JetNconstituents/src/JetNconstituents.cc

0001 #include "JetNconstituents.h"
0002
0003 // fun4all includes
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/PHTFileServer.h>
0006
0007 // phool includes
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/getClass.h>
0010
0011 // jet includes
0012 #include <jetbase/Jet.h>
0013 #include <jetbase/Jetv2.h>
0014 #include <jetbase/JetContainer.h>
0015 #include <jetbase/JetContainerv1.h>
0016
0017 // calobase includes
0018 #include <calobase/RawTower.h>
0019 #include <calobase/RawTowerContainer.h>
0020 #include <calobase/RawTowerGeom.h>
0021 #include <calobase/RawTowerGeomContainer.h>
0022 #include <calobase/TowerInfoContainer.h>
0023 #include <calobase/TowerInfo.h>
0024
0025 // jetbackground includes
0026 #include <jetbackground/TowerBackground.h>
0027
0028 #include <TH2D.h>
0029 #include <TH1D.h>
0030
0031 #include <iostream>
0032 #include <string>
0033 #include <vector>
0034
0035 JetNconstituents::JetNconstituents(const std::string &outputfilename)
0036   : SubsysReco("JetNconstituents")
0037   , m_outputFileName(outputfilename)
0038   // these are all initialized but included here for clarity
0039   , m_etaRange(-1.1, 1.1)
0040   , m_ptRange(1.0, 1000)
0041   , m_recoJetName("AntiKt_Tower_r04")
0042   , h1_jetNconstituents_total(nullptr)
0043   , h1_jetNconstituents_IHCAL(nullptr)
0044   , h1_jetNconstituents_OHCAL(nullptr)
0045   , h1_jetNconstituents_CEMC(nullptr)
0046   , h2_jetNconstituents_vs_caloLayer(nullptr)
0047   , h1_jetFracE_IHCAL(nullptr)
0048   , h1_jetFracE_OHCAL(nullptr)
0049   , h1_jetFracE_CEMC(nullptr)
0050   , h2_jetFracE_vs_caloLayer(nullptr)
0051 {
0052 }
0053
0054
0055 int JetNconstituents::Init(PHCompositeNode *topNode)
0056 {
0057   if (Verbosity() > 0)
0058   {
0059     std::cout << "JetNconstituents::Init - Output to " << m_outputFileName << std::endl;
0060   }
0061
0062   // create output file
0063   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0064
0065
0066   // Initialize histograms
0067   const int N_const = 200;
0068   Double_t N_const_bins[N_const+1];
0069   for (int i = 0; i <= N_const; i++)
0070   {
0071     N_const_bins[i] = 1.0*i;
0072   }
0073
0074   const int N_frac = 100;
0075   double frac_max = 1.0;
0076   Double_t frac_bins[N_frac+1];
0077   for (int i = 0; i <= N_frac; i++)
0078   {
0079     frac_bins[i] = (frac_max/N_frac)*i;
0080   }
0081
0082   const int N_calo_layers = 3;
0083   Double_t calo_layers_bins[N_calo_layers+1] = {0.5, 1.5, 2.5, 3.5}; // emcal, ihcal, ohcal
0084
0085   h1_jetNconstituents_total = new TH1D("h1_jetNconstituents_total", "Jet N constituents", N_const, N_const_bins);
0086   h1_jetNconstituents_IHCAL = new TH1D("h1_jetNconstituents_IHCAL", "Jet N constituents in IHCal", N_const, N_const_bins);
0087   h1_jetNconstituents_OHCAL = new TH1D("h1_jetNconstituents_OHCAL", "Jet N constituents in OHCal", N_const, N_const_bins);
0088   h1_jetNconstituents_CEMC = new TH1D("h1_jetNconstituents_CEMC", "Jet N constituents in CEMC", N_const, N_const_bins);
0089   h2_jetNconstituents_vs_caloLayer = new TH2D("h2_jetNconstituents_vs_caloLayer", "Jet N constituents vs Calo Layer", N_calo_layers, calo_layers_bins, N_const, N_const_bins);
0090
0091   h1_jetFracE_IHCAL = new TH1D("h1_jetFracE_IHCAL", "Jet E fraction in IHCal", N_frac, frac_bins);
0092   h1_jetFracE_OHCAL = new TH1D("h1_jetFracE_OHCAL", "Jet E fraction in OHCal", N_frac, frac_bins);
0093   h1_jetFracE_CEMC = new TH1D("h1_jetFracE_CEMC", "Jet E fraction in CEMC", N_frac, frac_bins);
0094   h2_jetFracE_vs_caloLayer = new TH2D("h2_jetFracE_vs_caloLayer", "Jet E fraction vs Calo Layer", N_calo_layers, calo_layers_bins, N_frac, frac_bins);
0095
0096   if(Verbosity() > 0)
0097   {
0098     std::cout << "JetNconstituents::Init - Histograms created" << std::endl;
0099   }
0100
0101   return Fun4AllReturnCodes::EVENT_OK;
0102 }
0103
0104 int JetNconstituents::process_event(PHCompositeNode *topNode)
0105 {
0106
0107   if(Verbosity() > 1)
0108   {
0109     std::cout << "JetNconstituents::process_event - Process event..." << std::endl;
0110   }
0111
0112   // get the jets
0113   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0114   if(!jets)
0115   {
0116     std::cout <<"JetNconstituents::process_event - Error can not find jet node " << m_recoJetName << std::endl;
0117     exit(-1); // fatal
0118   }
0119
0120   // get unsub towers
0121   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0122   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0123   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0124   if(!towersEM3 || !towersIH3 || !towersOH3)
0125   {
0126     std::cout <<"JetNconstituents::process_event - Error can not find tower node " << std::endl;
0127     exit(-1); // fatal
0128   }
0129
0130   // get tower geometry
0131   RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0132   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0133   if(!tower_geomIH || !tower_geomOH)
0134   {
0135     std::cout <<"JetNconstituents::process_event - Error can not find tower geometry node " << std::endl;
0136     exit(-1); // fatal
0137   }
0138
0139   // get underlying event
0140   float background_v2 = 0;
0141   float background_Psi2 = 0;
0142   bool has_tower_background = false;
0143   TowerBackground *towBack = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0144   if(!towBack)
0145   {
0146     std::cout <<"JetNconstituents::process_event - Error can not find tower background node " << std::endl;
0147   }
0148   else
0149   {
0150     has_tower_background = true;
0151     background_v2 = towBack->get_v2();
0152     background_Psi2 = towBack->get_Psi2();
0153   }
0154
0155
0156   // loop over jets
0157   for (auto jet : *jets)
0158   {
0159
0160     // remove noise
0161     if(jet->get_pt() < 1) continue;
0162
0163     // apply eta and pt cuts
0164     bool eta_cut = (jet->get_eta() >= m_etaRange.first) and (jet->get_eta() <= m_etaRange.second);
0165     bool pt_cut = (jet->get_pt() >= m_ptRange.first) and (jet->get_pt() <= m_ptRange.second);
0166     if ((not eta_cut) or (not pt_cut)) continue;
0167
0168
0169     // zero out counters
0170     int n_comp_total = 0;
0171     int n_comp_ihcal = 0;
0172     int n_comp_ohcal = 0;
0173     int n_comp_emcal = 0;
0174
0175     float jet_total_eT = 0;
0176     float eFrac_ihcal = 0;
0177     float eFrac_ohcal = 0;
0178     float eFrac_emcal = 0;
0179
0180     // loop over jet constituents
0181     for (auto comp: jet->get_comp_vec())
0182       {
0183
0184       // get tower
0185       unsigned int channel = comp.second;
0186       TowerInfo *tower;
0187
0188       float tower_eT = 0;
0189
0190       // get tower info
0191       if(comp.first == 26 || comp.first == 30)
0192       { // IHcal
0193
0194         tower = towersIH3->get_tower_at_channel(channel);
0195
0196         if(!tower || !tower_geomIH){ continue; }
0197
0198         unsigned int calokey = towersIH3->encode_key(channel);
0199         int ieta = towersIH3->getTowerEtaBin(calokey);
0200           int iphi = towersIH3->getTowerPhiBin(calokey);
0201           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0202           float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0203           float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0204         tower_eT = tower->get_energy()/cosh(tower_eta);
0205
0206         if(comp.first == 30)
0207         { // is sub1
0208           if(has_tower_background)
0209           {
0210             float UE = towBack->get_UE(1).at(ieta);
0211             float tower_UE = UE *( 1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0212             tower_eT = (tower->get_energy() - tower_UE)/cosh(tower_eta);
0213           }
0214         }
0215
0216         eFrac_ihcal+= tower_eT;
0217         jet_total_eT += tower_eT;
0218         n_comp_ihcal++;
0219         n_comp_total++;
0220
0221       }
0222       else if(comp.first == 27 || comp.first == 31)
0223       { // OHcal
0224
0225         tower = towersOH3->get_tower_at_channel(channel);
0226
0227         if(!tower || !tower_geomOH){ continue; }
0228
0229         unsigned int calokey = towersOH3->encode_key(channel);
0230         int ieta = towersOH3->getTowerEtaBin(calokey);
0231         int iphi = towersOH3->getTowerPhiBin(calokey);
0232         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0233         float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0234         float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0235         tower_eT = tower->get_energy()/cosh(tower_eta);
0236
0237         if(comp.first == 31)
0238         { // is sub1
0239           if(has_tower_background)
0240           {
0241             float UE = towBack->get_UE(2).at(ieta);
0242             float tower_UE = UE *( 1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0243             tower_eT = (tower->get_energy() - tower_UE)/cosh(tower_eta);
0244           }
0245         }
0246
0247         eFrac_ohcal+= tower_eT;
0248         jet_total_eT += tower_eT;
0249         n_comp_ohcal++;
0250         n_comp_total++;
0251
0252       }
0253       else if(comp.first == 28 || comp.first == 29)
0254       { // EMCal (retowered)
0255
0256         tower = towersEM3->get_tower_at_channel(channel);
0257
0258         if(!tower || !tower_geomIH){ continue; }
0259
0260
0261         unsigned int calokey = towersEM3->encode_key(channel);
0262         int ieta = towersEM3->getTowerEtaBin(calokey);
0263         int iphi = towersEM3->getTowerPhiBin(calokey);
0264         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0265         float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0266         float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0267         tower_eT = tower->get_energy()/cosh(tower_eta);
0268
0269         if(comp.first == 29)
0270         { // is sub1
0271           if(has_tower_background)
0272           {
0273             float UE = towBack->get_UE(0).at(ieta);
0274             float tower_UE = UE *( 1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0275             tower_eT = (tower->get_energy() - tower_UE)/cosh(tower_eta);
0276           }
0277         }
0278
0279         eFrac_emcal+= tower_eT;
0280         jet_total_eT += tower_eT;
0281         n_comp_emcal++;
0282         n_comp_total++;
0283       }
0284
0285
0286     }
0287
0288     // normalize energy fractions
0289     eFrac_ihcal/= jet_total_eT;
0290     eFrac_ohcal/= jet_total_eT;
0291     eFrac_emcal/= jet_total_eT;
0292
0293     // fill histograms
0294     h1_jetNconstituents_total->Fill(1.0*n_comp_total);
0295     h1_jetNconstituents_IHCAL->Fill(1.0*n_comp_ihcal);
0296     h1_jetNconstituents_OHCAL->Fill(1.0*n_comp_ohcal);
0297     h1_jetNconstituents_CEMC->Fill(1.0*n_comp_emcal);
0298     h2_jetNconstituents_vs_caloLayer->Fill(1.0, 1.0*n_comp_emcal);
0299     h2_jetNconstituents_vs_caloLayer->Fill(2.0, 1.0*n_comp_ihcal);
0300     h2_jetNconstituents_vs_caloLayer->Fill(3.0, 1.0*n_comp_ohcal);
0301
0302     h1_jetFracE_IHCAL->Fill(eFrac_ihcal);
0303     h1_jetFracE_OHCAL->Fill(eFrac_ohcal);
0304     h1_jetFracE_CEMC->Fill(eFrac_emcal);
0305
0306     h2_jetFracE_vs_caloLayer->Fill(1.0, eFrac_emcal);
0307     h2_jetFracE_vs_caloLayer->Fill(2.0, eFrac_ihcal);
0308     h2_jetFracE_vs_caloLayer->Fill(3.0, eFrac_ohcal);
0309
0310
0311
0312
0313   }
0314
0315   if(Verbosity() > 1)
0316   {
0317     std::cout << "JetNconstituents::process_event - Event processed" << std::endl;
0318   }
0319
0320   return Fun4AllReturnCodes::EVENT_OK;
0321 }
0322
0323
0324 //____________________________________________________________________________..
0325 int JetNconstituents::End(PHCompositeNode *topNode)
0326 {
0327
0328   if(Verbosity() > 0)
0329   {
0330     std::cout << "JetNconstituents::EndRun - End run " << std::endl;
0331     std::cout << "JetNconstituents::EndRun - Writing to " << m_outputFileName << std::endl;
0332   }
0333
0334   PHTFileServer::get().cd(m_outputFileName);
0335
0336   // write histograms to root file
0337   h1_jetNconstituents_total->Write();
0338   h1_jetNconstituents_IHCAL->Write();
0339   h1_jetNconstituents_OHCAL->Write();
0340   h1_jetNconstituents_CEMC->Write();
0341   h2_jetNconstituents_vs_caloLayer->Write();
0342
0343   h1_jetFracE_IHCAL->Write();
0344   h1_jetFracE_OHCAL->Write();
0345   h1_jetFracE_CEMC->Write();
0346   h2_jetFracE_vs_caloLayer->Write();
0347
0348   if(Verbosity() > 0)
0349   {
0350     std::cout << "JetNconstituents::EndRun - Done" << std::endl;
0351   }
0352
0353   return Fun4AllReturnCodes::EVENT_OK;
0354 }
0355