calotriggeremulator/src/FindOutlier.cc

0001 #include "FindOutlier.h"
0002 #include <calobase/RawTowerGeomContainer.h>
0003 #include <calobase/RawTowerGeom.h>
0004 // Calo
0005 #include <calobase/TowerInfo.h>
0006 #include <fun4all/Fun4AllHistoManager.h>
0007 #include <jetbackground/TowerBackground.h>
0008
0009 #include <calobase/TowerInfoContainer.h>
0010 #include <ffarawobjects/CaloPacketContainer.h>
0011 #include <ffarawobjects/CaloPacket.h>
0012 // jet includes
0013 #include <jetbackground/TowerBackground.h>
0014 #include <jetbase/Jet.h>
0015 #include <jetbase/JetContainer.h>
0016 #include <jetbase/JetContainerv1.h>
0017 #include <jetbase/Jetv2.h>
0018
0019 #include <qautils/QAHistManagerDef.h>
0020
0021 #include <fun4all/Fun4AllHistoManager.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>  // for PHWHERE
0026
0027 #include <ffarawobjects/Gl1Packet.h>
0028 #include <TNtuple.h>
0029 #include <TFile.h>
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <TProfile2D.h>
0033 #include <TSystem.h>
0034
0035 #include <boost/format.hpp>
0036 #include <phool/PHCompositeNode.h>
0037 #include <phool/PHIODataNode.h>
0038 #include <phool/PHNode.h>
0039 #include <phool/PHNodeIterator.h>
0040 #include <phool/PHObject.h>
0041 #include <phool/getClass.h>
0042 #include <phool/phool.h>
0043 #include"TrashInfov1.h"
0044
0045 #include <cassert>
0046 #include <cmath>  // for log10, pow, sqrt, abs, M_PI
0047 #include <cstdint>
0048 #include <iostream>  // for operator<<, endl, basic_...
0049 #include <limits>
0050 #include <map>  // for operator!=, _Rb_tree_con...
0051 #include <string>
0052 #include <utility>  // for pair
0053
0054 FindOutlier::FindOutlier(const std::string& name, const std::string& outfile)
0055 : SubsysReco(name)
0056 , m_outfilename(outfile)
0057
0058 {
0059
0060 }
0061
0062 FindOutlier::~FindOutlier()
0063 {
0064
0065 }
0066
0067 int FindOutlier::Init(PHCompositeNode* /*unused*/)
0068 {
0069   std::cout << "Leaving FindOutlier::Init" << std::endl;
0070
0071   tn = new TNtuple("tn","tn","event:trash:r1:energy:emcal_energy:spread");
0072   hm = new Fun4AllHistoManager("FindOutlierHistos");
0073   h_emcal = new TH2F("EMCAL","", 12, -0.5, 95.5, 32, -0.5, 255.5);
0074   h_hcalout = new TH2F("HCALOUT","", 24, -0.5, 23.5, 64, -0.5, 63.5);
0075   h_hcalout_time = new TH2F("HCALOUT_TIME","", 24, -0.5, 23.5, 64, -0.5, 63.5);
0076   return Fun4AllReturnCodes::EVENT_OK;
0077 }
0078
0079 int FindOutlier::process_event(PHCompositeNode* topNode)
0080 {
0081
0082   _eventcounter++;
0083   if (Verbosity())  std::cout << "In process_event" << std::endl;
0084
0085   Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0086   if (!gl1PacketInfo)
0087   {
0088     std::cout << PHWHERE << "FindOutlier::process_event: GL1Packet node is missing" << std::endl;
0089     return Fun4AllReturnCodes::EVENT_OK;
0090   }
0091
0092   trashinfo = findNode::getClass<TrashInfo>(topNode, "TrashInfo");
0093   if (!trashinfo)
0094   {
0095     std::cout << PHWHERE << "FindOutlier::process_event: Trashinfo node is missing" << std::endl;
0096     return Fun4AllReturnCodes::EVENT_OK;
0097   }
0098   //  if (!trashinfo->isTrash()) return Fun4AllReturnCodes::EVENT_OK;
0099
0100   TowerInfoContainer* emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0101   h_emcal->Reset();
0102   int size = emcal_towers->size();
0103   for (int i = 0; i < size; i++)
0104     {
0105       TowerInfo* tower = emcal_towers->get_tower_at_channel(i);
0106       float offlineenergy = tower->get_energy();
0107
0108       unsigned int towerkey = emcal_towers->encode_key(i);
0109       int ieta = emcal_towers->getTowerEtaBin(towerkey);
0110       int iphi = emcal_towers->getTowerPhiBin(towerkey);
0111       if (offlineenergy < 0.03) continue;
0112       h_emcal->Fill(ieta, iphi, offlineenergy);
0113     }
0114   if (Verbosity())  std::cout << "emcal" << std::endl;
0115   char hname[100];
0116   TH2F *h = (TH2F*) h_emcal->Clone();
0117   sprintf(hname, "h_emcal_%d_%d", (trashinfo->isTrash()? 1 : 0), _eventcounter);
0118   h->SetName(hname);
0119   hm->registerHisto(h);
0120
0121   if (useWave)
0122     {
0123       CaloPacketContainer *emcalcont = findNode::getClass<CaloPacketContainer>(topNode, "CEMCPackets");
0124       for (int packetid = 6001; packetid <= 6128; packetid++)
0125     {
0126       CaloPacket *packet = emcalcont->getPacketbyId(packetid);
0127       if (!packet)
0128         return Fun4AllReturnCodes::ABORTRUN;
0129
0130       for (int i = 0; i < 192; i++)
0131         {
0132           int cha = 192*(packetid-6001) + i;
0133           if (emcal_towers->get_tower_at_channel(cha)->get_energy() < 0.03 ) continue;
0134           char name[100];
0135           sprintf(name, "wave_%d_%d", _eventcounter, cha);
0136           auto *h_waveform = new TH1F(name, ";sample; adc", 12, -0.5, 11.5);
0137
0138           int channel = i;
0139           for (int is = 0; is < 12; is++)
0140         {
0141           if (packet->iValue(channel, "SUPPRESSED"))
0142             {
0143               if (is < 6)
0144             {
0145               h_waveform->Fill(is,packet->iValue(channel, "PRE"));
0146             }
0147               else
0148             {
0149               h_waveform->Fill(is,packet->iValue(channel, "POST"));
0150             }
0151             }
0152           else
0153             {
0154               h_waveform->Fill(is,packet->iValue(is, channel));
0155             }
0156
0157         }
0158           h_waveform->SetDirectory(nullptr);
0159           hm->registerHisto(h_waveform);
0160         }
0161     }
0162     }
0163
0164   TowerInfoContainer* hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0165   h_hcalout->Reset();
0166   h_hcalout_time->Reset();
0167   size = hcalout_towers->size();
0168   for (int i = 0; i < size; i++)
0169     {
0170       TowerInfo* tower = hcalout_towers->get_tower_at_channel(i);
0171       float offlineenergy = tower->get_energy();
0172       float offlinetime = tower->get_time();
0173       unsigned int towerkey = hcalout_towers->encode_key(i);
0174       int ieta = hcalout_towers->getTowerEtaBin(towerkey);
0175       int iphi = hcalout_towers->getTowerPhiBin(towerkey);
0176       if (offlineenergy < 0.03) continue;
0177       h_hcalout->Fill(ieta, iphi, offlineenergy);
0178       h_hcalout_time->Fill(ieta, iphi, offlinetime);
0179     }
0180   if (Verbosity())  std::cout << "hcal" << std::endl;
0181   h = (TH2F*) h_hcalout->Clone();
0182   sprintf(hname, "h_hcalout_%d_%d", (trashinfo->isTrash()? 1 : 0), _eventcounter);
0183   h->SetName(hname);
0184   hm->registerHisto(h);
0185   TH2F *ht = (TH2F*) h_hcalout_time->Clone();
0186   sprintf(hname, "h_hcalout_time_%d_%d", (trashinfo->isTrash()? 1 : 0), _eventcounter);
0187   ht->SetName(hname);
0188   hm->registerHisto(ht);
0189
0190   tn->Fill(_eventcounter, (trashinfo->isTrash()? 1:0), trashinfo->getR1(), trashinfo->getEnergy(), trashinfo->getEMCALEnergy(), trashinfo->getSpread());
0191
0192   if (useWave)
0193     {
0194
0195       CaloPacketContainer *emcalcont = findNode::getClass<CaloPacketContainer>(topNode, "HCALPackets");
0196       for (int packetid = 8001; packetid <= 8008; packetid++)
0197     {
0198       CaloPacket *packet = emcalcont->getPacketbyId(packetid);
0199       if (!packet)
0200         return Fun4AllReturnCodes::ABORTRUN;
0201
0202       std::cout << "FindOutlier::Found one at " << _eventcounter <<std::endl;
0203       for (int i = 0; i < 192; i++)
0204         {
0205           int cha = 192*(packetid-8001) + i;
0206           if (hcalout_towers->get_tower_at_channel(cha)->get_energy() < 0.03 ) continue;
0207           char name[100];
0208           sprintf(name, "wave_%d_%d", _eventcounter, cha);
0209           auto *h_waveform = new TH1F(name, ";sample; adc", 12, -0.5, 11.5);
0210
0211           int channel = i;
0212           for (int is = 0; is < 12; is++)
0213         {
0214           if (packet->iValue(channel, "SUPPRESSED"))
0215             {
0216               if (is < 6)
0217             {
0218               h_waveform->Fill(is,packet->iValue(channel, "PRE"));
0219             }
0220               else
0221             {
0222               h_waveform->Fill(is,packet->iValue(channel, "POST"));
0223             }
0224             }
0225           else
0226             {
0227               h_waveform->Fill(is,packet->iValue(is, channel));
0228             }
0229
0230         }
0231           h_waveform->SetDirectory(nullptr);
0232           hm->registerHisto(h_waveform);
0233         }
0234     }
0235     }
0236   // uint64_t triggervec = gl1PacketInfo->getScaledVector();
0237       // bool mbd_event = (((triggervec >> 10U) & 0x1U) == 0x1U);
0238       // if (!mbd_event)   return Fun4AllReturnCodes::EVENT_OK;
0239
0240
0241   // uint64_t rawvec = gl1PacketInfo->lValue(0, "TriggerInput");
0242   // bool rare_event = ((rawvec >> 18U) > 0);
0243   // if (rare_event)   return Fun4AllReturnCodes::EVENT_OK;
0244
0245   // // calculate max
0246   // TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0247   // int size = towers->size();  // online towers should be the same!
0248   // h_emcal->Reset();
0249
0250   // for (int i = 0; i < size; i++)
0251   //   {
0252   //       TowerInfo* tower = towers->get_tower_at_channel(i);
0253   //       float offlineenergy = tower->get_energy();
0254   //       unsigned int towerkey = towers->encode_key(i);
0255   //       int ieta = towers->getTowerEtaBin(towerkey);
0256   //       int iphi = towers->getTowerPhiBin(towerkey);
0257   //       bool isGood = tower->get_isGood();
0258   //    if (!isGood) continue;
0259   //    h_emcal->Fill(ieta, iphi, offlineenergy);
0260   //    if (offlineenergy > max_cut) {
0261
0262   //      std::cout << i << " " <<ieta<< " " << iphi <<std::endl;
0263   //    }
0264   //   }
0265   // int binnumber = h_emcal->GetMaximumBin();
0266   // if (h_emcal->GetBinContent(binnumber) < max_cut)  return Fun4AllReturnCodes::EVENT_OK;
0267
0268   // int binx = -1;
0269   // int binz = -1;
0270   // int biny = -1;
0271   // h_emcal->GetBinXYZ(binnumber, binx, biny, binz);
0272   // if (binx < 0)
0273   //   return Fun4AllReturnCodes::ABORTRUN;
0274
0275   // binx--;
0276   // biny--;
0277   // int packetid = 4 * biny;
0278   // int channeloffset = 0;
0279   // if (binx < 3)
0280   //   {
0281   //     channeloffset = 2 - binx;
0282   //     packetid += 1;
0283   //   }
0284   // else if (binx < 6)
0285   //   {
0286   //     channeloffset = 2 - (binx%3);
0287   //     packetid += 0;
0288   //   }
0289   // else if (binx < 9)
0290   //   {
0291   //     channeloffset = (binx%3);
0292   //     packetid += 2;
0293   //   }
0294   // else
0295   //   {
0296   //     channeloffset = (binx%3);
0297   //     packetid += 3;
0298   //   }
0299   // channeloffset *= 64;
0300   // packetid += 6001;
0301
0302   // CaloPacketContainer *emcalcont = findNode::getClass<CaloPacketContainer>(topNode, "CEMCPackets");
0303   // CaloPacket *packet = emcalcont->getPacketbyId(packetid);
0304   // if (!packet)
0305   //   return Fun4AllReturnCodes::ABORTRUN;
0306
0307   // std::cout << "FindOutlier::Found one at " << _eventcounterer <<" ch " << binx << " " << biny << " " << std::endl;
0308   // for (int i = 0; i < 192; i++)
0309   //   {
0310   //     char name[100];
0311   //     sprintf(name, "wave_%d_%d", i, _eventcounteryer);
0312   //     auto *h_waveform = new TH1F(name, ";sample; adc", 12, -0.5, 11.5);
0313
0314   //     int channel = i;
0315   //     for (int is = 0; is < 12; is++)
0316   //    {
0317   //      if (packet->iValue(channel, "SUPPRESSED"))
0318   //        {
0319   //          if (is < 6)
0320   //        {
0321   //          h_waveform->Fill(is,packet->iValue(channel, "PRE"));
0322   //        }
0323   //          else
0324   //        {
0325   //          h_waveform->Fill(is,packet->iValue(channel, "POST"));
0326   //        }
0327   //        }
0328   //      else
0329   //        {
0330   //          h_waveform->Fill(is,packet->iValue(is, channel));
0331   //        }
0332
0333   //    }
0334   //     h_waveform->SetDirectory(nullptr);
0335   //     hm->registerHisto(h_waveform);
0336   //   }
0337
0338
0339   return Fun4AllReturnCodes::EVENT_OK;
0340
0341
0342 }
0343
0344 int FindOutlier::End(PHCompositeNode* /*unused*/)
0345 {
0346   hm->dumpHistos(m_outfilename.c_str());
0347   std::string fail=m_outfilename;
0348   auto it = fail.find("HIST");
0349   fail.replace(it, 4, "NTUPLE");
0350   TFile *fout = new TFile(fail.c_str(),"recreate");
0351   tn->Write();
0352   fout->Close();
0353   return Fun4AllReturnCodes::EVENT_OK;
0354 }