calotriggeremulator/src/JetPerformance.cc

0001 #include "JetPerformance.h"
0002 //for emc clusters
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TH2D.h>
0006 #include <TH1D.h>
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008 #include <g4main/PHG4Particle.h>
0009 #include <TMath.h>
0010 #include <calotrigger/TriggerRunInfov1.h>
0011
0012 #include <calobase/RawTowerGeom.h>
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 #include <calobase/TowerInfoContainer.h>
0019 #include <calobase/TowerInfoContainerv1.h>
0020 #include <calobase/TowerInfoContainerv2.h>
0021 #include <calobase/TowerInfoContainerv3.h>
0022 #include <calobase/TowerInfo.h>
0023 #include <calobase/TowerInfov1.h>
0024 #include <calobase/TowerInfov2.h>
0025 #include <calobase/TowerInfov3.h>
0026
0027 #include <calobase/RawCluster.h>
0028 #include <calobase/RawClusterContainer.h>
0029 #include <calobase/RawClusterUtility.h>
0030 #include <calobase/RawTowerGeomContainer.h>
0031 #include <calobase/RawTower.h>
0032 #include <calobase/RawTowerContainer.h>
0033 #include <fun4all/Fun4AllHistoManager.h>
0034 #include <HepMC/SimpleVector.h>
0035 //for vetex information
0036 #include <globalvertex/GlobalVertex.h>
0037 #include <globalvertex/GlobalVertexMap.h>
0038 #include "TrashInfov1.h"
0039 #include <vector>
0040
0041 #include <fun4all/Fun4AllReturnCodes.h>
0042 #include <phool/PHCompositeNode.h>
0043 #include <phool/PHIODataNode.h>
0044 #include <phool/PHNode.h>
0045 #include <phool/PHNodeIterator.h>
0046 #include <phool/PHObject.h>
0047 #include <phool/getClass.h>
0048 // G4Cells includes
0049
0050 #include <iostream>
0051
0052 #include <map>
0053
0054 //____________________________________________________________________________..
0055 JetPerformance::JetPerformance(const std::string &name, const std::string &outfilename):
0056   SubsysReco(name)
0057
0058 {
0059   isSim = false;
0060   _foutname = outfilename;
0061   m_calo_nodename = "TOWERINFO_CALIB";
0062   m_jet_triggernames = {"Jet 6 GeV + MBD NS >= 1",
0063             "Jet 8 GeV + MBD NS >= 1",
0064             "Jet 10 GeV + MBD NS >= 1",
0065             "Jet 12 GeV + MBD NS >= 1",
0066             "Jet 6 GeV",
0067             "Jet 8 GeV",
0068             "Jet 10 GeV",
0069             "Jet 12 GeV"};
0070 }
0071
0072 //____________________________________________________________________________..
0073 JetPerformance::~JetPerformance()
0074 {
0075
0076 }
0077
0078 //____________________________________________________________________________..
0079 int JetPerformance::Init(PHCompositeNode *topNode)
0080 {
0081
0082   triggeranalyzer = new TriggerAnalyzer();
0083
0084   hm = new Fun4AllHistoManager("DIJET_HISTOGRAMS");
0085
0086   TH1D *h;
0087
0088   h = new TH1D("h_jet_et_mbd", ";et;N", 50, 0, 50);
0089   hm->registerHisto(h);
0090   h = new TH1D("h_jet_pt_mbd", ";pt;N", 50, 0, 50);
0091   hm->registerHisto(h);
0092
0093   h = new TH1D("h_jet_et_scaled_mbd", ";et;N", 50, 0, 50);
0094   hm->registerHisto(h);
0095   h = new TH1D("h_jet_pt_scaled_mbd", ";pt;N", 50, 0, 50);
0096   hm->registerHisto(h);
0097
0098   for (int i = 0 ; i < 8; i++)
0099     {
0100       h = new TH1D(Form("h_jet_et_scaled_gl1_%d", i), ";et;N", 50, 0, 50);
0101       hm->registerHisto(h);
0102
0103       h = new TH1D(Form("h_jet_pt_scaled_gl1_%d", i), ";pt;N", 50, 0, 50);
0104       hm->registerHisto(h);
0105
0106       h = new TH1D(Form("h_jet_et_raw_gl1_%d", i), ";et;N", 50, 0, 50);
0107       hm->registerHisto(h);
0108
0109       h = new TH1D(Form("h_jet_pt_raw_gl1_%d", i), ";pt;N", 50, 0, 50);
0110       hm->registerHisto(h);
0111     }
0112
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 //____________________________________________________________________________..
0118 int JetPerformance::InitRun(PHCompositeNode *topNode)
0119 {
0120   return Fun4AllReturnCodes::EVENT_OK;
0121 }
0122
0123 //____________________________________________________________________________..
0124
0125 int JetPerformance::process_event(PHCompositeNode *topNode)
0126 {
0127
0128
0129   _i_event++;
0130
0131   triggeranalyzer->decodeTriggers(topNode);
0132
0133
0134   if (mbd_prescale == -2)
0135     {
0136       mbd_prescale = triggeranalyzer->getPrescaleByName("MBD N&S >= 1");
0137     }
0138   for (int i = 0; i < 8 ; i++)
0139     {
0140       if (jet_prescale[i] == -2)
0141     {
0142       jet_prescale[i] = triggeranalyzer->getPrescaleByName(m_jet_triggernames[i]);
0143     }
0144     }
0145
0146   bool mbd_fired = triggeranalyzer->didTriggerFire("MBD N&S >= 1");
0147   bool jet_fired[8];
0148   bool jet_raw[8];
0149   for (int i = 0; i < 8; i++)
0150     {
0151       jet_fired[i] = triggeranalyzer->didTriggerFire(m_jet_triggernames);
0152       jet_raw[i] = triggeranalyzer->checkRawTrigger(m_jet_triggernames);
0153     }
0154   RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0155   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0156
0157
0158   std::string recoJetName = "AntiKt_Tower_r04_Sub1";
0159
0160   JetContainer *jets_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0161
0162   TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALIN");
0163   TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALOUT");
0164   TowerInfoContainer *emcalre_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_CEMC_RETOWER");
0165
0166   float background_v2 = 0;
0167   float background_Psi2 = 0;
0168   bool has_tower_background = false;
0169   TowerBackground *towBack = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0170   if (!towBack)
0171     {
0172
0173     }
0174   else
0175     {
0176       has_tower_background = true;
0177       background_v2 = towBack->get_v2();
0178       background_Psi2 = towBack->get_Psi2();
0179     }
0180
0181   if (!jets_4)
0182     {
0183       return Fun4AllReturnCodes::ABORTRUN;
0184     }
0185
0186   auto hem = dynamic_cast<TH1*>(hm->getHisto("h_jet_et_mbd"));// ";et;N", 50, 0, 50);
0187   auto hpm = dynamic_cast<TH1*>(hm->getHisto("h_jet_pt_mbd"));// ";pt;N", 50, 0, 50);
0188   auto hesm = dynamic_cast<TH1*>(hm->getHisto("h_jet_et_scaled_mbd"));// ";et;N", 50, 0, 50);
0189   auto hpsm = dynamic_cast<TH1*>(hm->getHisto("h_jet_pt_scaled_mbd"));// ";pt;N", 50, 0, 50);
0190   TH1D *hesg[8];
0191   TH1D *hpsg[8];
0192   TH1D *herg[8];
0193   TH1D *hprg[8];
0194   for (int i = 0 ; i < 8; i++)
0195     {
0196       hesg[i] = dynamic_cast<TH1*>(hm->getHisto(Form("h_jet_et_scaled_gl1_%d", i)));// ";et;N", 50, 0, 50);
0197       hpsg[i] = dynamic_cast<TH1*>(hm->getHisto(Form("h_jet_pt_scaled_gl1_%d", i)));// ";pt;N", 50, 0, 50);
0198       herg[i] = dynamic_cast<TH1*>(hm->getHisto(Form("h_jet_et_raw_gl1_%d", i)));// ";et;N", 50, 0, 50);
0199       hprg[i] = dynamic_cast<TH1*>(hm->getHisto(Form("h_jet_pt_raw_gl1_%d", i)));// ";pt;N", 50, 0, 50);
0200     }
0201
0202
0203
0204
0205   for (auto jet : *jets_4)
0206     {
0207       float jet_total_eT = 0;
0208       float eFrac_ihcal = 0;
0209       float eFrac_ohcal = 0;
0210       float eFrac_emcal = 0;
0211
0212       for (auto comp : jet->get_comp_vec())
0213     {
0214
0215       unsigned int channel = comp.second;
0216       TowerInfo *tower;
0217       float tower_eT = 0;
0218       if (comp.first == 26 || comp.first == 30)
0219         {  // IHcal
0220           tower = hcalin_towers->get_tower_at_channel(channel);
0221           if (!tower || !tower_geomIH)
0222         {
0223           continue;
0224         }
0225           if(!tower->get_isGood()) continue;
0226
0227           unsigned int calokey = hcalin_towers->encode_key(channel);
0228
0229           int ieta = hcalin_towers->getTowerEtaBin(calokey);
0230           int iphi = hcalin_towers->getTowerPhiBin(calokey);
0231           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0232           float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0233           float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0234           tower_eT = tower->get_energy() / std::cosh(tower_eta);
0235
0236           if (comp.first == 30)
0237         {  // is sub1
0238           if (has_tower_background)
0239             {
0240               float UE = towBack->get_UE(1).at(ieta);
0241               float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0242               tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0243             }
0244         }
0245
0246           eFrac_ihcal += tower_eT;
0247           jet_total_eT += tower_eT;
0248         }
0249       else if (comp.first == 27 || comp.first == 31)
0250         {  // IHcal
0251           tower = hcalout_towers->get_tower_at_channel(channel);
0252
0253           if (!tower || !tower_geomOH)
0254         {
0255           continue;
0256         }
0257
0258           unsigned int calokey = hcalout_towers->encode_key(channel);
0259           int ieta = hcalout_towers->getTowerEtaBin(calokey);
0260           int iphi = hcalout_towers->getTowerPhiBin(calokey);
0261           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0262           float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0263           float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0264           tower_eT = tower->get_energy() / std::cosh(tower_eta);
0265
0266           if (comp.first == 31)
0267         {  // is sub1
0268           if (has_tower_background)
0269             {
0270               float UE = towBack->get_UE(2).at(ieta);
0271               float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0272               tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0273             }
0274         }
0275
0276           eFrac_ohcal += tower_eT;
0277           jet_total_eT += tower_eT;
0278         }
0279       else if (comp.first == 28 || comp.first == 29)
0280         {  // IHcal
0281           tower = emcalre_towers->get_tower_at_channel(channel);
0282
0283           if (!tower || !tower_geomIH)
0284         {
0285           continue;
0286         }
0287
0288           unsigned int calokey = emcalre_towers->encode_key(channel);
0289           int ieta = emcalre_towers->getTowerEtaBin(calokey);
0290           int iphi = emcalre_towers->getTowerPhiBin(calokey);
0291           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0292           float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0293           float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0294           tower_eT = tower->get_energy() / std::cosh(tower_eta);
0295
0296           if (comp.first == 29)
0297         {  // is sub1
0298           if (has_tower_background)
0299             {
0300               float UE = towBack->get_UE(0).at(ieta);
0301               float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0302               tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0303             }
0304         }
0305
0306
0307           eFrac_emcal += tower_eT;
0308           jet_total_eT += tower_eT;
0309         }
0310     }
0311
0312       if (jet->get_pt() > pt_cut)
0313     }
0314   return Fun4AllReturnCodes::EVENT_OK;
0315 }
0316
0317
0318
0319 void JetPerformance::GetNodes(PHCompositeNode* topNode)
0320 {
0321
0322
0323 }
0324
0325 int JetPerformance::ResetEvent(PHCompositeNode *topNode)
0326 {
0327
0328   return Fun4AllReturnCodes::EVENT_OK;
0329 }
0330
0331 //____________________________________________________________________________..
0332 int JetPerformance::EndRun(const int runnumber)
0333 {
0334   return Fun4AllReturnCodes::EVENT_OK;
0335 }
0336
0337 //____________________________________________________________________________..
0338 int JetPerformance::End(PHCompositeNode *topNode)
0339 {
0340   if (Verbosity() > 0)
0341     {
0342       std::cout << "JetPerformance::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0343     }
0344   std::cout<<"Total events: "<<_i_event<<std::endl;
0345   hm->dumpHistos(_foutname.c_str());
0346
0347   return Fun4AllReturnCodes::EVENT_OK;
0348 }
0349
0350 //____________________________________________________________________________..
0351 int JetPerformance::Reset(PHCompositeNode *topNode)
0352 {
0353   if (Verbosity() > 0)
0354     {
0355       std::cout << "JetPerformance::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0356     }
0357   return Fun4AllReturnCodes::EVENT_OK;
0358 }