QA/Jet/ConstituentsinJets.cc

0001 #include "ConstituentsinJets.h"
0002
0003 // calobase includes
0004 #include <calobase/RawTower.h>
0005 #include <calobase/RawTowerContainer.h>
0006 #include <calobase/RawTowerGeom.h>
0007 #include <calobase/RawTowerGeomContainer.h>
0008 #include <calobase/TowerInfo.h>
0009 #include <calobase/TowerInfoContainer.h>
0010
0011 // calotrigger includes
0012 #include <calotrigger/TriggerAnalyzer.h>
0013
0014 // jet includes
0015 #include <jetbase/Jet.h>
0016 #include <jetbase/JetContainer.h>
0017 #include <jetbase/JetContainerv1.h>
0018 #include <jetbase/Jetv2.h>
0019
0020 // fun4all includes
0021 #include <fun4all/Fun4AllHistoManager.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023
0024 // jetbackground includes
0025 #include <jetbackground/TowerBackground.h>
0026
0027 // phool includes
0028 #include <phool/PHCompositeNode.h>
0029 #include <phool/getClass.h>
0030
0031 // qautils includes
0032 #include <qautils/QAHistManagerDef.h>
0033
0034 #include <TH1.h>
0035 #include <TH2.h>
0036 #include <TStyle.h>
0037
0038 #include <algorithm>
0039 #include <cassert>
0040 #include <cmath>
0041 #include <cstdlib>
0042 #include <iostream>
0043 #include <map>
0044 #include <string>
0045 #include <utility>
0046 #include <vector>
0047
0048 ConstituentsinJets::ConstituentsinJets(const std::string &moduleName, const std::string &recojetname, const std::string &towBkgdName, const std::string &histTag, const std::string &towCEMCName, const std::string &towIHCALName, const std::string &towOHCALName)
0049   : SubsysReco(moduleName)
0050   , m_moduleName(moduleName)
0051   , m_recoJetName(recojetname)
0052   , m_towBkgdName(towBkgdName)
0053   , m_histTag(histTag)
0054   , m_towCEMCName(towCEMCName)
0055   , m_towIHCALName(towIHCALName)
0056   , m_towOHCALName(towOHCALName)
0057 {
0058 }
0059
0060 int ConstituentsinJets::Init(PHCompositeNode * /*topNode*/)
0061 {
0062   // create output file
0063   delete m_analyzer; // make cppcheck happy
0064   m_analyzer = new TriggerAnalyzer();
0065
0066   gStyle->SetOptTitle(0);
0067   m_manager = QAHistManagerDef::getHistoManager();
0068   if (!m_manager)
0069   {
0070     std::cout << PHWHERE << ": PANIC: couldn't grab histogram manager!" << std::endl;
0071     assert(m_manager);
0072   }
0073   // Initialize histograms
0074   const int N_const = 200;
0075   Double_t N_const_bins[N_const + 1];
0076   for (int i = 0; i <= N_const; i++)
0077   {
0078     N_const_bins[i] = 1.0 * i;
0079   }
0080
0081   const int N_frac = 100;
0082   double frac_max = 1.0;
0083   Double_t frac_bins[N_frac + 1];
0084   for (int i = 0; i <= N_frac; i++)
0085   {
0086     frac_bins[i] = (frac_max / N_frac) * i;
0087   }
0088
0089   const int N_calo_layers = 3;
0090   Double_t calo_layers_bins[N_calo_layers + 1] = {0.5, 1.5, 2.5, 3.5};  // emcal, ihcal, ohcal
0091
0092   // make sure module name is lower case
0093   std::string smallModuleName = m_moduleName;
0094   std::transform(
0095       smallModuleName.begin(),
0096       smallModuleName.end(),
0097       smallModuleName.begin(),
0098       ::tolower);
0099
0100   // construct histogram names
0101   std::vector<std::string> vecHistNames = {
0102       "ncsts_total",
0103       "ncsts_ihcal",
0104       "ncsts_ohcal",
0105       "ncsts_cemc",
0106       "ncstsvscalolayer",
0107       "efracjet_ihcal",
0108       "efracjet_ohcal",
0109       "efracjet_cemc",
0110       "efracjetvscalolayer"};
0111   for (auto &vecHistName : vecHistNames)
0112   {
0113     vecHistName.insert(0, "h_" + smallModuleName + "_");
0114     if (!m_histTag.empty())
0115     {
0116       vecHistName.append("_" + m_histTag);
0117     }
0118   }
0119
0120   // declare histograms and include x and y axis labels in the constructor
0121   h1_ConstituentsinJets_total = new TH1D(vecHistNames[0].data(), "", N_const, N_const_bins);
0122   h1_ConstituentsinJets_total->GetXaxis()->SetTitle("N constituents");
0123   h1_ConstituentsinJets_total->GetYaxis()->SetTitle("Counts");
0124
0125   h1_ConstituentsinJets_IHCAL = new TH1D(vecHistNames[1].data(), "", N_const, N_const_bins);
0126   h1_ConstituentsinJets_IHCAL->GetXaxis()->SetTitle("N constituents");
0127   h1_ConstituentsinJets_IHCAL->GetYaxis()->SetTitle("Counts");
0128
0129   h1_ConstituentsinJets_OHCAL = new TH1D(vecHistNames[2].data(), "", N_const, N_const_bins);
0130   h1_ConstituentsinJets_OHCAL->GetXaxis()->SetTitle("N constituents");
0131   h1_ConstituentsinJets_OHCAL->GetYaxis()->SetTitle("Counts");
0132
0133   h1_ConstituentsinJets_CEMC = new TH1D(vecHistNames[3].data(), "", N_const, N_const_bins);
0134   h1_ConstituentsinJets_CEMC->GetXaxis()->SetTitle("N constituents");
0135   h1_ConstituentsinJets_CEMC->GetYaxis()->SetTitle("Counts");
0136
0137   h2_ConstituentsinJets_vs_caloLayer = new TH2D(vecHistNames[4].data(), "", N_calo_layers, calo_layers_bins, N_const, N_const_bins);
0138   h2_ConstituentsinJets_vs_caloLayer->GetXaxis()->SetTitle("Calo Layer");
0139   h2_ConstituentsinJets_vs_caloLayer->GetYaxis()->SetTitle("N constituents");
0140
0141   h1_jetFracE_IHCAL = new TH1D(vecHistNames[5].data(), "", N_frac, frac_bins);
0142   h1_jetFracE_IHCAL->GetXaxis()->SetTitle("E fraction");
0143   h1_jetFracE_IHCAL->GetYaxis()->SetTitle("Counts");
0144
0145   h1_jetFracE_OHCAL = new TH1D(vecHistNames[6].data(), "", N_frac, frac_bins);
0146   h1_jetFracE_OHCAL->GetXaxis()->SetTitle("E fraction");
0147   h1_jetFracE_OHCAL->GetYaxis()->SetTitle("Counts");
0148
0149   h1_jetFracE_CEMC = new TH1D(vecHistNames[7].data(), "", N_frac, frac_bins);
0150   h1_jetFracE_CEMC->GetXaxis()->SetTitle("E fraction");
0151   h1_jetFracE_CEMC->GetYaxis()->SetTitle("Counts");
0152
0153   h2_jetFracE_vs_caloLayer = new TH2D(vecHistNames[8].data(), "", N_calo_layers, calo_layers_bins, N_frac, frac_bins);
0154   h2_jetFracE_vs_caloLayer->GetXaxis()->SetTitle("Calo Layer");
0155   h2_jetFracE_vs_caloLayer->GetYaxis()->SetTitle("E fraction");
0156
0157   // register histograms
0158   // m_manager->registerHisto(h1_ConstituentsinJets_total);
0159   // m_manager->registerHisto(h1_ConstituentsinJets_IHCAL);
0160   // m_manager->registerHisto(h1_ConstituentsinJets_OHCAL);
0161   // m_manager->registerHisto(h1_ConstituentsinJets_CEMC);
0162   // m_manager->registerHisto(h2_ConstituentsinJets_vs_caloLayer);
0163   // m_manager->registerHisto(h1_jetFracE_IHCAL);
0164   // m_manager->registerHisto(h1_jetFracE_OHCAL);
0165   // m_manager->registerHisto(h1_jetFracE_CEMC);
0166   // m_manager->registerHisto(h2_jetFracE_vs_caloLayer);
0167
0168   if (Verbosity() > 0)
0169   {
0170     std::cout << "ConstituentsinJets::Init - Histograms created" << std::endl;
0171   }
0172
0173   return Fun4AllReturnCodes::EVENT_OK;
0174 }
0175
0176 int ConstituentsinJets::process_event(PHCompositeNode *topNode)
0177 {
0178   if (Verbosity() > 1)
0179   {
0180     std::cout << "ConstituentsinJets::process_event - Process event..." << std::endl;
0181   }
0182
0183   // if needed, check if selected trigger fired
0184   if (m_doTrgSelect)
0185   {
0186     m_analyzer->decodeTriggers(topNode);
0187     bool hasTrigger = JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer);
0188     if (!hasTrigger)
0189     {
0190       return Fun4AllReturnCodes::EVENT_OK;
0191     }
0192   }
0193
0194   // get the jets
0195   JetContainer *jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0196   if (!jets)
0197   {
0198     std::cout << "ConstituentsinJets::process_event - Error can not find jet node " << m_recoJetName << std::endl;
0199     return Fun4AllReturnCodes::EVENT_OK;
0200   }
0201
0202   // get unsub towers
0203   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, m_towCEMCName);
0204   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, m_towIHCALName);
0205   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, m_towOHCALName);
0206   if (!towersEM3 || !towersIH3 || !towersOH3)
0207   {
0208     std::cout << "ConstituentsinJets::process_event - Error can not find tower node " << std::endl;
0209     return Fun4AllReturnCodes::EVENT_OK;
0210   }
0211
0212   // get tower geometry
0213   RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0214   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0215   if (!tower_geomIH || !tower_geomOH)
0216   {
0217     std::cout << "ConstituentsinJets::process_event - Error can not find tower geometry node " << std::endl;
0218     return Fun4AllReturnCodes::EVENT_OK;
0219   }
0220
0221   // for tower background if needed
0222   TowerBackground *towBack = nullptr;
0223
0224   // get underlying event
0225   float background_v2 = 0;
0226   float background_Psi2 = 0;
0227   bool has_tower_background = false;
0228   if (!m_inPPMode)
0229   {
0230     towBack = findNode::getClass<TowerBackground>(topNode, m_towBkgdName);
0231     if (!towBack)
0232     {
0233       std::cout << "ConstituentsinJets::process_event - Error can not find tower background node " << std::endl;
0234     }
0235     else
0236     {
0237       has_tower_background = true;
0238       background_v2 = towBack->get_v2();
0239       background_Psi2 = towBack->get_Psi2();
0240     }
0241   }
0242
0243   // loop over jets
0244   for (auto *jet : *jets)
0245   {
0246     // remove noise
0247     if (jet->get_pt() < m_ptRange.first)//change from <1 to <m_ptRange.first
0248     {
0249       continue;
0250     }
0251
0252     // apply eta and pt cuts
0253     bool eta_cut = (jet->get_eta() >= m_etaRange.first) and (jet->get_eta() <= m_etaRange.second);
0254     bool pt_cut = (jet->get_pt() >= m_ptRange.first) and (jet->get_pt() <= m_ptRange.second);
0255     if ((not eta_cut) or (not pt_cut))
0256     {
0257       continue;
0258     }
0259
0260     // zero out counters
0261     int n_comp_total = 0;
0262     int n_comp_ihcal = 0;
0263     int n_comp_ohcal = 0;
0264     int n_comp_emcal = 0;
0265
0266     float jet_total_eT = 0;
0267     float eFrac_ihcal = 0;
0268     float eFrac_ohcal = 0;
0269     float eFrac_emcal = 0;
0270
0271     // loop over jet constituents
0272     for (auto comp : jet->get_comp_vec())
0273     {
0274       // get tower
0275       unsigned int channel = comp.second;
0276       TowerInfo *tower;
0277
0278       float tower_eT = 0;
0279
0280       // get tower info
0281       if (comp.first == 26 || comp.first == 30)
0282       {  // IHcal
0283
0284         tower = towersIH3->get_tower_at_channel(channel);
0285
0286         if (!tower || !tower_geomIH)
0287         {
0288           continue;
0289         }
0290
0291         unsigned int calokey = towersIH3->encode_key(channel);
0292         int ieta = towersIH3->getTowerEtaBin(calokey);
0293         int iphi = towersIH3->getTowerPhiBin(calokey);
0294         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0295         float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0296         float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0297         tower_eT = tower->get_energy() / std::cosh(tower_eta);
0298
0299         if (comp.first == 30)
0300         {  // is sub1
0301           if (has_tower_background)
0302           {
0303             float UE = towBack->get_UE(1).at(ieta);
0304             float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0305             tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0306           }
0307         }
0308
0309         eFrac_ihcal += tower_eT;
0310         jet_total_eT += tower_eT;
0311         n_comp_ihcal++;
0312         n_comp_total++;
0313       }
0314       else if (comp.first == 27 || comp.first == 31)
0315       {  // OHcal
0316
0317         tower = towersOH3->get_tower_at_channel(channel);
0318
0319         if (!tower || !tower_geomOH)
0320         {
0321           continue;
0322         }
0323
0324         unsigned int calokey = towersOH3->encode_key(channel);
0325         int ieta = towersOH3->getTowerEtaBin(calokey);
0326         int iphi = towersOH3->getTowerPhiBin(calokey);
0327         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0328         float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0329         float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0330         tower_eT = tower->get_energy() / std::cosh(tower_eta);
0331
0332         if (comp.first == 31)
0333         {  // is sub1
0334           if (has_tower_background)
0335           {
0336             float UE = towBack->get_UE(2).at(ieta);
0337             float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0338             tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0339           }
0340         }
0341
0342         eFrac_ohcal += tower_eT;
0343         jet_total_eT += tower_eT;
0344         n_comp_ohcal++;
0345         n_comp_total++;
0346       }
0347       else if (comp.first == 28 || comp.first == 29)
0348       {  // EMCal (retowered)
0349
0350         tower = towersEM3->get_tower_at_channel(channel);
0351
0352         if (!tower || !tower_geomIH)
0353         {
0354           continue;
0355         }
0356
0357         unsigned int calokey = towersEM3->encode_key(channel);
0358         int ieta = towersEM3->getTowerEtaBin(calokey);
0359         int iphi = towersEM3->getTowerPhiBin(calokey);
0360         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0361         float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0362         float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0363         tower_eT = tower->get_energy() / std::cosh(tower_eta);
0364
0365         if (comp.first == 29)
0366         {  // is sub1
0367           if (has_tower_background)
0368           {
0369             float UE = towBack->get_UE(0).at(ieta);
0370             float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0371             tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0372           }
0373         }
0374
0375         eFrac_emcal += tower_eT;
0376         jet_total_eT += tower_eT;
0377         n_comp_emcal++;
0378         n_comp_total++;
0379       }
0380     }
0381
0382     // normalize energy fractions
0383     eFrac_ihcal /= jet_total_eT;
0384     eFrac_ohcal /= jet_total_eT;
0385     eFrac_emcal /= jet_total_eT;
0386
0387     // fill histograms
0388     assert(h1_ConstituentsinJets_total);
0389     assert(h1_ConstituentsinJets_IHCAL);
0390     assert(h1_ConstituentsinJets_OHCAL);
0391     assert(h1_ConstituentsinJets_CEMC);
0392     assert(h2_ConstituentsinJets_vs_caloLayer);
0393     assert(h1_jetFracE_IHCAL);
0394     assert(h1_jetFracE_OHCAL);
0395     assert(h1_jetFracE_CEMC);
0396     assert(h2_jetFracE_vs_caloLayer);
0397
0398     h1_ConstituentsinJets_total->Fill(1.0 * n_comp_total);
0399     h1_ConstituentsinJets_IHCAL->Fill(1.0 * n_comp_ihcal);
0400     h1_ConstituentsinJets_OHCAL->Fill(1.0 * n_comp_ohcal);
0401     h1_ConstituentsinJets_CEMC->Fill(1.0 * n_comp_emcal);
0402     h2_ConstituentsinJets_vs_caloLayer->Fill(1.0, 1.0 * n_comp_emcal);
0403     h2_ConstituentsinJets_vs_caloLayer->Fill(2.0, 1.0 * n_comp_ihcal);
0404     h2_ConstituentsinJets_vs_caloLayer->Fill(3.0, 1.0 * n_comp_ohcal);
0405
0406     h1_jetFracE_IHCAL->Fill(eFrac_ihcal);
0407     h1_jetFracE_OHCAL->Fill(eFrac_ohcal);
0408     h1_jetFracE_CEMC->Fill(eFrac_emcal);
0409
0410     h2_jetFracE_vs_caloLayer->Fill(1.0, eFrac_emcal);
0411     h2_jetFracE_vs_caloLayer->Fill(2.0, eFrac_ihcal);
0412     h2_jetFracE_vs_caloLayer->Fill(3.0, eFrac_ohcal);
0413   }
0414
0415   if (Verbosity() > 1)
0416   {
0417     std::cout << "ConstituentsinJets::process_event - Event processed" << std::endl;
0418   }
0419
0420   return Fun4AllReturnCodes::EVENT_OK;
0421 }
0422
0423 int ConstituentsinJets::End(PHCompositeNode * /*topNode*/)
0424 {
0425   if (Verbosity() > 0)
0426   {
0427     std::cout << "ConstituentsinJets::EndRun - End run " << std::endl;
0428     // std::cout << "ConstituentsinJets::EndRun - Writing to " << m_outputFileName << std::endl;
0429   }
0430
0431   m_manager->registerHisto(h1_ConstituentsinJets_total);
0432   m_manager->registerHisto(h1_ConstituentsinJets_IHCAL);
0433   m_manager->registerHisto(h1_ConstituentsinJets_OHCAL);
0434   m_manager->registerHisto(h1_ConstituentsinJets_CEMC);
0435   m_manager->registerHisto(h2_ConstituentsinJets_vs_caloLayer);
0436   m_manager->registerHisto(h1_jetFracE_IHCAL);
0437   m_manager->registerHisto(h1_jetFracE_OHCAL);
0438   m_manager->registerHisto(h1_jetFracE_CEMC);
0439   m_manager->registerHisto(h2_jetFracE_vs_caloLayer);
0440
0441   // write histograms to root file
0442   // h1_ConstituentsinJets_total->Write();
0443   // h1_ConstituentsinJets_IHCAL->Write();
0444   // h1_ConstituentsinJets_OHCAL->Write();
0445   // h1_ConstituentsinJets_CEMC->Write();
0446   // h2_ConstituentsinJets_vs_caloLayer->Write();
0447
0448   // h1_jetFracE_IHCAL->Write();
0449   // h1_jetFracE_OHCAL->Write();
0450   // h1_jetFracE_CEMC->Write();
0451   // h2_jetFracE_vs_caloLayer->Write();
0452
0453   if (Verbosity() > 0)
0454   {
0455     std::cout << "ConstituentsinJets::EndRun - Done" << std::endl;
0456   }
0457
0458   return Fun4AllReturnCodes::EVENT_OK;
0459 }