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