calotriggeremulator/src/FindRejects.cc

0001 #include "FindRejects.h"
0002 #include <calobase/RawTowerGeomContainer.h>
0003 #include <calobase/RawTowerGeom.h>
0004 #include <calotrigger/TriggerAnalyzer.h>
0005 // Calo
0006 #include <calobase/TowerInfo.h>
0007 #include <fun4all/Fun4AllHistoManager.h>
0008 #include <jetbackground/TowerBackground.h>
0009
0010 #include <calobase/TowerInfoContainer.h>
0011 #include <ffarawobjects/CaloPacketContainer.h>
0012 #include <ffarawobjects/CaloPacket.h>
0013 // jet includes
0014 #include <jetbackground/TowerBackground.h>
0015 #include <jetbase/Jet.h>
0016 #include <jetbase/JetContainer.h>
0017 #include <jetbase/JetContainerv1.h>
0018 #include <jetbase/Jetv2.h>
0019
0020 #include <qautils/QAHistManagerDef.h>
0021
0022 #include <fun4all/Fun4AllHistoManager.h>
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>  // for PHWHERE
0027
0028 #include <ffarawobjects/Gl1Packet.h>
0029 #include <TNtuple.h>
0030 #include <TFile.h>
0031 #include <TH1.h>
0032 #include <TH2.h>
0033 #include <TProfile2D.h>
0034 #include <TSystem.h>
0035
0036 #include <boost/format.hpp>
0037 #include <phool/PHCompositeNode.h>
0038 #include <phool/PHIODataNode.h>
0039 #include <phool/PHNode.h>
0040 #include <phool/PHNodeIterator.h>
0041 #include <phool/PHObject.h>
0042 #include <phool/getClass.h>
0043 #include <phool/phool.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 FindRejects::FindRejects(const std::string& name, const std::string& outfile)
0055 : SubsysReco(name)
0056 , m_outfilename(outfile)
0057
0058 {
0059   m_triggeranalyzer = new TriggerAnalyzer();
0060 }
0061
0062 FindRejects::~FindRejects()
0063 {
0064
0065 }
0066
0067 int FindRejects::Init(PHCompositeNode* /*unused*/)
0068 {
0069   std::cout << "Leaving FindRejects::Init" << std::endl;
0070
0071   tn = new TNtuple("tn","tn","event:trash:r1:energy:emcal_energy:spread");
0072   hm = new Fun4AllHistoManager("FindRejectsHistos");
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 FindRejects::process_event(PHCompositeNode* topNode)
0080 {
0081
0082   _eventcounter++;
0083   std::cout << "In process_event" << std::endl;
0084
0085   std::string recoJetName = "AntiKt_Tower_r04_Sub1";
0086   JetContainer *jets_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0087   TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0088   TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0089   TowerInfoContainer *emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0090   TowerInfoContainer *emcalre_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0091   RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0092   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0093
0094   int b_njet_4 = 0;
0095   std::vector<float> b_jet_et_4{};
0096
0097   if (jets_4)
0098      {
0099        std::cout << "In jets" << std::endl;
0100        for (auto jet : *jets_4)
0101      {
0102        if (jet->get_pt() < 3) continue;
0103        std::cout << "In b_njet " << b_njet_4 << std::endl;
0104
0105        float jet_total_eT = 0;
0106
0107        b_njet_4++;
0108
0109
0110        for (auto comp : jet->get_comp_vec())
0111          {
0112            unsigned int channel = comp.second;
0113            TowerInfo *tower;
0114            float tower_eT = 0;
0115            if (comp.first == 26 || comp.first == 30)
0116          {  // IHcal
0117            tower = hcalin_towers->get_tower_at_channel(channel);
0118            if (!tower || !tower_geomIH)
0119              {
0120                continue;
0121              }
0122            if(!tower->get_isGood()) continue;
0123
0124            unsigned int calokey = hcalin_towers->encode_key(channel);
0125
0126            int ieta = hcalin_towers->getTowerEtaBin(calokey);
0127            int iphi = hcalin_towers->getTowerPhiBin(calokey);
0128            const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0129            float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0130            tower_eT = tower->get_energy() / std::cosh(tower_eta);
0131            jet_total_eT += tower_eT;
0132          }
0133            else if (comp.first == 27 || comp.first == 31)
0134          {  // IHcal
0135            tower = hcalout_towers->get_tower_at_channel(channel);
0136
0137            if (!tower || !tower_geomOH)
0138              {
0139                continue;
0140              }
0141            unsigned int calokey = hcalout_towers->encode_key(channel);
0142            int ieta = hcalout_towers->getTowerEtaBin(calokey);
0143            int iphi = hcalout_towers->getTowerPhiBin(calokey);
0144            const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0145            float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0146            tower_eT = tower->get_energy() / std::cosh(tower_eta);
0147            jet_total_eT += tower_eT;
0148          }
0149            else if (comp.first == 28 || comp.first == 29)
0150          {  // IHcal
0151            tower = emcalre_towers->get_tower_at_channel(channel);
0152
0153            if (!tower || !tower_geomIH)
0154              {
0155                continue;
0156              }
0157
0158            unsigned int calokey = emcalre_towers->encode_key(channel);
0159            int ieta = emcalre_towers->getTowerEtaBin(calokey);
0160            int iphi = emcalre_towers->getTowerPhiBin(calokey);
0161            const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0162            float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0163            tower_eT = tower->get_energy() / std::cosh(tower_eta);
0164            jet_total_eT += tower_eT;
0165          }
0166          }
0167        b_jet_et_4.push_back(jet_total_eT);
0168        std::cout << "   total et = " << jet_total_eT  << std::endl;
0169      }
0170      }
0171    // find if there is a jet over 20 gev
0172
0173    std::cout << b_njet_4 << " jets"<<std::endl;
0174    if (b_njet_4 < 1) return Fun4AllReturnCodes::EVENT_OK;
0175
0176    auto is_greater_than_cut = [] (float et) { return et >= 15; };
0177
0178    auto it = std::find_if(b_jet_et_4.begin(), b_jet_et_4.end(), is_greater_than_cut);
0179
0180    if (it == b_jet_et_4.end())
0181      {
0182        return Fun4AllReturnCodes::EVENT_OK;
0183      }
0184
0185
0186    std::cout << " ET > 15 FOUND " << (*it) << " GeV" << std::endl;
0187
0188    m_triggeranalyzer->decodeTriggers(topNode);
0189
0190    if (m_triggeranalyzer->didTriggerFire(19))
0191      {
0192        return Fun4AllReturnCodes::EVENT_OK;
0193      }
0194
0195    std::cout << " REJECT FOUND " << (*it) << " GeV" << std::endl;
0196
0197    {
0198      h_emcal->Reset();
0199      int size = emcal_towers->size();
0200      for (int i = 0; i < size; i++)
0201        {
0202      TowerInfo* tower = emcal_towers->get_tower_at_channel(i);
0203      float offlineenergy = tower->get_energy();
0204
0205      unsigned int towerkey = emcal_towers->encode_key(i);
0206      int ieta = emcal_towers->getTowerEtaBin(towerkey);
0207      int iphi = emcal_towers->getTowerPhiBin(towerkey);
0208      if (offlineenergy < 0.03) continue;
0209      h_emcal->Fill(ieta, iphi, offlineenergy);
0210        }
0211
0212      char hname[100];
0213      TH2F *h = (TH2F*) h_emcal->Clone();
0214      sprintf(hname, "h_emcal_%d", _eventcount);
0215      h->SetName(hname);
0216      hm->registerHisto(h);
0217
0218      CaloPacketContainer *emcalcont = findNode::getClass<CaloPacketContainer>(topNode, "CEMCPackets");
0219      for (int packetid = 6001; packetid <= 6128; packetid++)
0220        {
0221      CaloPacket *packet = emcalcont->getPacketbyId(packetid);
0222      if (!packet)
0223        return Fun4AllReturnCodes::ABORTRUN;
0224
0225      for (int i = 0; i < 192; i++)
0226        {
0227          int cha = 192*(packetid-6001) + i;
0228          if (emcal_towers->get_tower_at_channel(cha)->get_energy() < 0.03 ) continue;
0229          char name[100];
0230          sprintf(name, "ewave_%d_%d", _eventcount, cha);
0231          auto *h_waveform = new TH1F(name, ";sample; adc", 12, -0.5, 11.5);
0232
0233          int channel = i;
0234          for (int is = 0; is < 12; is++)
0235            {
0236          if (packet->iValue(channel, "SUPPRESSED"))
0237            {
0238              if (is < 6)
0239                {
0240              h_waveform->Fill(is,packet->iValue(channel, "PRE"));
0241                }
0242              else
0243                {
0244              h_waveform->Fill(is,packet->iValue(channel, "POST"));
0245                }
0246            }
0247          else
0248            {
0249              h_waveform->Fill(is,packet->iValue(is, channel));
0250            }
0251
0252            }
0253          h_waveform->SetDirectory(nullptr);
0254          hm->registerHisto(h_waveform);
0255        }
0256        }
0257
0258    }
0259    {
0260      h_hcalout->Reset();
0261
0262      int size = hcalout_towers->size();
0263      for (int i = 0; i < size; i++)
0264        {
0265      TowerInfo* tower = hcalout_towers->get_tower_at_channel(i);
0266      float offlineenergy = tower->get_energy();
0267      float offlinetime = tower->get_time();
0268      unsigned int towerkey = hcalout_towers->encode_key(i);
0269      int ieta = hcalout_towers->getTowerEtaBin(towerkey);
0270      int iphi = hcalout_towers->getTowerPhiBin(towerkey);
0271      if (offlineenergy < 0.03) continue;
0272      h_hcalout->Fill(ieta, iphi, offlineenergy);
0273      h_hcalout_time->Fill(ieta, iphi, offlinetime);
0274        }
0275
0276      char hname[100];
0277      TH2F *h = (TH2F*) h_hcalout->Clone();
0278      sprintf(hname, "h_hcalout_%d", _eventcount);
0279      h->SetName(hname);
0280      hm->registerHisto(h);
0281      TH2F *ht = (TH2F*) h_hcalout_time->Clone();
0282      sprintf(hname, "h_hcalout_time_%d", _eventcount);
0283      ht->SetName(hname);
0284      hm->registerHisto(ht);
0285
0286      CaloPacketContainer *emcalcont = findNode::getClass<CaloPacketContainer>(topNode, "HCALPackets");
0287      for (int packetid = 8001; packetid <= 8008; packetid++)
0288        {
0289      CaloPacket *packet = emcalcont->getPacketbyId(packetid);
0290      if (!packet)
0291        return Fun4AllReturnCodes::ABORTRUN;
0292
0293      std::cout << "FindRejects::Found one at " << _eventcounter <<std::endl;
0294      for (int i = 0; i < 192; i++)
0295        {
0296          int cha = 192*(packetid-8001) + i;
0297          if (hcalout_towers->get_tower_at_channel(cha)->get_energy() < 0.03 ) continue;
0298          char name[100];
0299          sprintf(name, "hwave_%d_%d", _eventcount, cha);
0300          auto *h_waveform = new TH1F(name, ";sample; adc", 12, -0.5, 11.5);
0301
0302          int channel = i;
0303          for (int is = 0; is < 12; is++)
0304            {
0305          if (packet->iValue(channel, "SUPPRESSED"))
0306            {
0307              if (is < 6)
0308                {
0309              h_waveform->Fill(is,packet->iValue(channel, "PRE"));
0310                }
0311              else
0312                {
0313              h_waveform->Fill(is,packet->iValue(channel, "POST"));
0314                }
0315            }
0316          else
0317            {
0318              h_waveform->Fill(is,packet->iValue(is, channel));
0319            }
0320
0321            }
0322          h_waveform->SetDirectory(nullptr);
0323          hm->registerHisto(h_waveform);
0324        }
0325        }
0326    }
0327
0328    _eventcount++;
0329      if (_eventcount > 1000)
0330        {
0331      return Fun4AllReturnCodes::ABORTRUN;
0332        }
0333      return Fun4AllReturnCodes::EVENT_OK;
0334
0335
0336 }
0337
0338 int FindRejects::End(PHCompositeNode* /*unused*/)
0339 {
0340   hm->dumpHistos(m_outfilename.c_str());
0341   TFile *fout = new TFile("failtest.root","recreate");
0342   tn->Write();
0343   fout->Close();
0344   return Fun4AllReturnCodes::EVENT_OK;
0345 }