QA/Jet/JetSeedCount.cc

0001 #include "JetSeedCount.h"
0002
0003 #include <calotrigger/TriggerAnalyzer.h>
0004
0005 #include <centrality/CentralityInfo.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <fun4all/PHTFileServer.h>
0009
0010 #include <globalvertex/GlobalVertex.h>
0011 #include <globalvertex/GlobalVertexMap.h>
0012
0013 #include <jetbase/JetContainer.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/getClass.h>
0017
0018 #include <TH1.h>
0019 #include <TH2.h>
0020 #include <TStyle.h>
0021
0022 #include <algorithm>
0023 #include <cstdlib>
0024 #include <iostream>
0025
0026 JetSeedCount::JetSeedCount(const std::string &moduleName, const std::string &recojetname, const std::string &rawSeedName, const std::string &subSeedName, const std::string &truthjetname, const std::string &outputfilename)
0027   : SubsysReco(moduleName)
0028   , m_moduleName(moduleName)
0029   , m_recoJetName(recojetname)
0030   , m_rawSeedName(rawSeedName)
0031   , m_subSeedName(subSeedName)
0032   , m_truthJetName(truthjetname)
0033   , m_outputFileName(outputfilename)
0034 {
0035   // std::cout << "JetSeedCount::JetSeedCount(const std::string &name) Calling ctor" << std::endl;
0036
0037 }
0038
0039 int JetSeedCount::Init(PHCompositeNode * /*topNode*/)
0040 {
0041   if (Verbosity() > 1)
0042   {
0043     std::cout << "JetSeedCount::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0044   }
0045   if (m_writeToOutputFile)
0046   {
0047     std::cout << "Opening output file named " << m_outputFileName << std::endl;
0048     PHTFileServer::open(m_outputFileName, "RECREATE");
0049   }
0050   delete m_analyzer;
0051   m_analyzer = new TriggerAnalyzer();
0052
0053   gStyle->SetOptTitle(0);
0054   m_manager = QAHistManagerDef::getHistoManager();
0055   if (!m_manager)
0056   {
0057     std::cout << "No m_manager" << std::endl;
0058     assert(m_manager);
0059   }
0060
0061   // make sure module name is lower case
0062   std::string smallModuleName = m_moduleName;
0063   std::transform(
0064       smallModuleName.begin(),
0065       smallModuleName.end(),
0066       smallModuleName.begin(),
0067       ::tolower);
0068
0069   // construct histogram names
0070   std::vector<std::string> vecHistNames = {
0071       "rawseedcount",
0072       "rawpt",
0073       "rawptall",
0074       "rawetavsphi",
0075       "subseedcount",
0076       "subpt",
0077       "subptall",
0078       "subetavsphi",
0079       "rawseedenergyvscent",
0080       "subseedenergyvscent",
0081       "centmbd",
0082       "rawseedvscent",
0083       "subseedvscent"};
0084   for (auto &vecHistName : vecHistNames)
0085   {
0086     vecHistName.insert(0, "h_" + smallModuleName + "_");
0087     if (!m_histTag.empty())
0088     {
0089       vecHistName.append("_" + m_histTag);
0090     }
0091   }
0092
0093   // make histograms
0094   m_hRawSeedCount = new TH1F(vecHistNames[0].data(), "", 100, 0.00, 50.00);
0095   m_hRawSeedCount->GetXaxis()->SetTitle("Raw Seed Count per Event");
0096   m_hRawSeedCount->GetYaxis()->SetTitle("Number of Entries");
0097   m_manager->registerHisto(m_hRawSeedCount);
0098
0099   m_hRawPt = new TH1F(vecHistNames[1].data(), "", 1000, 0.00, 50.00);
0100   m_hRawPt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0101   m_hRawPt->GetYaxis()->SetTitle("Number of Entries");
0102   m_manager->registerHisto(m_hRawPt);
0103
0104   m_hRawPt_All = new TH1F(vecHistNames[2].data(), "", 1000, 0.00, 50.00);
0105   m_hRawPt_All->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0106   m_hRawPt_All->GetYaxis()->SetTitle("Number of Entries");
0107   m_manager->registerHisto(m_hRawPt_All);
0108
0109   m_hRawEtaVsPhi = new TH2F(vecHistNames[3].data(), "", 220, -1.1, 1.1, 628, -3.14, 3.14);
0110   m_hRawEtaVsPhi->GetXaxis()->SetTitle("Jet #eta [Rads.]");
0111   m_hRawEtaVsPhi->GetYaxis()->SetTitle("Jet #phi [Rads.]");
0112   m_manager->registerHisto(m_hRawEtaVsPhi);
0113
0114   m_hSubSeedCount = new TH1F(vecHistNames[4].data(), "", 100, 0.00, 50.00);
0115   m_hSubSeedCount->GetXaxis()->SetTitle("Sub Seed Count per Event");
0116   m_hSubSeedCount->GetYaxis()->SetTitle("Number of Entries");
0117   m_manager->registerHisto(m_hSubSeedCount);
0118
0119   m_hSubPt = new TH1F(vecHistNames[5].data(), "", 1000, 0.00, 50.00);
0120   m_hSubPt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0121   m_hSubPt->GetYaxis()->SetTitle("Number of Entries");
0122   m_manager->registerHisto(m_hSubPt);
0123
0124   m_hSubPt_All = new TH1F(vecHistNames[6].data(), "", 1000, 0.00, 50.00);
0125   m_hSubPt_All->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0126   m_hSubPt_All->GetYaxis()->SetTitle("Number of Entries");
0127   m_manager->registerHisto(m_hSubPt_All);
0128
0129   m_hSubEtaVsPhi = new TH2F(vecHistNames[7].data(), "", 220, -1.1, 1.1, 628, -3.14, 3.14);
0130   m_hSubEtaVsPhi->GetXaxis()->SetTitle("Jet #eta [Rads.]");
0131   m_hSubEtaVsPhi->GetYaxis()->SetTitle("Jet #phi [Rads.]");
0132   m_manager->registerHisto(m_hSubEtaVsPhi);
0133
0134   // If not in pp mode, plot quantities vs. centrality
0135   if (!m_inPPMode)
0136   {
0137     m_hRawSeedEnergyVsCent = new TH2F(vecHistNames[8].data(), "", 10.00, 0.00, 100.00, 100, 0.00, 50.00);
0138     m_hRawSeedEnergyVsCent->GetXaxis()->SetTitle("Centrality");
0139     m_hRawSeedEnergyVsCent->GetYaxis()->SetTitle("RawSeedEnergy");
0140     m_manager->registerHisto(m_hRawSeedEnergyVsCent);
0141
0142     m_hSubSeedEnergyVsCent = new TH2F(vecHistNames[9].data(), "", 10.00, 0.00, 100.00, 100, 0.00, 50.00);
0143     m_hSubSeedEnergyVsCent->GetXaxis()->SetTitle("Centrality");
0144     m_hSubSeedEnergyVsCent->GetYaxis()->SetTitle("SubSeedEnergy");
0145     m_manager->registerHisto(m_hSubSeedEnergyVsCent);
0146
0147     m_hCentMbd = new TH1F(vecHistNames[10].data(), "", 10, 0.00, 100.00);
0148     m_hCentMbd->GetXaxis()->SetTitle("Centrality (Mbd)");
0149     m_hCentMbd->GetYaxis()->SetTitle("Number of Entries");
0150     m_manager->registerHisto(m_hCentMbd);
0151
0152     m_hRawSeedVsCent = new TH2F(vecHistNames[11].data(), "", 10, 0.00, 100.00, 101, -0.5, 100.5);
0153     m_hRawSeedVsCent->GetXaxis()->SetTitle("Centrality");
0154     m_hRawSeedVsCent->GetYaxis()->SetTitle("Raw Seed Count");
0155     m_manager->registerHisto(m_hRawSeedVsCent);
0156
0157     m_hSubSeedVsCent = new TH2F(vecHistNames[12].data(), "", 10, 0.00, 100.00, 101, -0.5, 100.5);
0158     m_hSubSeedVsCent->GetXaxis()->SetTitle("Centrality");
0159     m_hSubSeedVsCent->GetYaxis()->SetTitle("Sub Seed Count");
0160     m_manager->registerHisto(m_hSubSeedVsCent);
0161   }  // end if not in pp mode
0162   return Fun4AllReturnCodes::EVENT_OK;
0163 }
0164
0165 //____________________________________________________________________________..
0166 int JetSeedCount::InitRun(PHCompositeNode * /*topNode*/)
0167 {
0168   if (Verbosity() > 1)
0169   {
0170     std::cout << "JetSeedCount::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0171   }
0172   return Fun4AllReturnCodes::EVENT_OK;
0173 }
0174
0175 //____________________________________________________________________________..
0176 int JetSeedCount::process_event(PHCompositeNode *topNode)
0177 {
0178   // if needed, check if selected trigger fired
0179   if (m_doTrgSelect)
0180   {
0181     m_analyzer->decodeTriggers(topNode);
0182     bool hasTrigger = JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer);
0183     if (!hasTrigger)
0184     {
0185       return Fun4AllReturnCodes::EVENT_OK;
0186     }
0187   }
0188
0189   // std::cout << "JetSeedCount::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0190   ++m_event;
0191   // Calling Raw Jet Seeds
0192   JetContainer *seedjetsraw = findNode::getClass<JetContainer>(topNode, m_rawSeedName);
0193   if (!seedjetsraw)
0194   {
0195     std::cout << "JetSeedCount::process_event - Error can not find DST raw seed jets" << std::endl;
0196     return Fun4AllReturnCodes::EVENT_OK;
0197   }
0198
0199   // Calling Sub jet seeds
0200   JetContainer *seedjetssub = findNode::getClass<JetContainer>(topNode, m_subSeedName);
0201   if (!seedjetssub)
0202   {
0203     std::cout << "JetSeedCount::process_event - Error can not find DST sub seed jets" << std::endl;
0204     return Fun4AllReturnCodes::EVENT_OK;
0205   }
0206
0207   // If not in pp mode, call Centrality Info
0208   CentralityInfo *cent_node = nullptr;
0209   if (!m_inPPMode)
0210   {
0211     cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0212     if (!cent_node)
0213     {
0214       std::cout << "JetSeedCount::process_event - Error can not find CentralityInfo" << std::endl;
0215       return Fun4AllReturnCodes::EVENT_OK;
0216     }
0217   }
0218
0219   // Z vertex info
0220   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0221   if (!vertexmap)
0222   {
0223     std::cout
0224         << "JetSeedCount::process_event - Error can not find global vertex  node "
0225         << std::endl;
0226     return Fun4AllReturnCodes::EVENT_OK;
0227   }
0228   if (vertexmap->empty())
0229   {
0230     if (Verbosity() > 1)
0231     {
0232       std::cout
0233           << "JetSeedCount::process_event - global vertex node is empty "
0234           << std::endl;
0235     }
0236     return Fun4AllReturnCodes::EVENT_OK;
0237   }
0238
0239   // If not in pp mode, save centrailty info
0240   float cent_mbd = std::numeric_limits<float>::max();
0241   if (!m_inPPMode)
0242   {
0243     cent_mbd = cent_node->get_centile(CentralityInfo::PROP::bimp);
0244     m_hCentMbd->Fill(cent_mbd);
0245   }
0246
0247   // Raw Seed Count
0248   uint64_t n_seed_raw = 0;
0249   //  float Counter = 0;
0250   // for (JetMap::Iter iter = seedjetsraw->begin(); iter != seedjetsraw->end(); ++iter){
0251   for (auto *jet : *seedjetsraw)
0252   {
0253     // Jet* jet = iter->second;
0254     int passesCut = jet->get_property(seedjetsraw->property_index(Jet::PROPERTY::prop_SeedItr));
0255     // Counter += 1;
0256     m_hRawPt_All->Fill(jet->get_pt());
0257     if (passesCut == 1)
0258     {
0259       m_hRawPt->Fill(jet->get_pt());
0260       m_hRawEtaVsPhi->Fill(jet->get_eta(), jet->get_phi());
0261       if (!m_inPPMode)
0262       {
0263         m_hRawSeedEnergyVsCent->Fill(cent_mbd, jet->get_e());
0264       }
0265       n_seed_raw++;
0266     }
0267   }
0268   m_hRawSeedCount->Fill(n_seed_raw);
0269   if (!m_inPPMode)
0270   {
0271     m_hRawSeedVsCent->Fill(cent_mbd, n_seed_raw);
0272   }
0273
0274   // Sub Seed Count
0275   uint64_t n_seed_sub = 0;
0276   //  Counter = 0;
0277   // for (unsigned int iter = 0; iter < seedjetssub->size(); ++iter){
0278   // Jet* jet = seedjetsub->get_jet(iter);
0279   for (auto *jet : *seedjetssub)
0280   {
0281     // Jet* jet = iter->second;
0282     int passesCut = jet->get_property(seedjetssub->property_index(Jet::PROPERTY::prop_SeedItr));
0283     //  Counter += 1;
0284     m_hSubPt_All->Fill(jet->get_pt());
0285     if (passesCut == 2)
0286     {
0287       m_hSubPt->Fill(jet->get_pt());
0288       m_hSubEtaVsPhi->Fill(jet->get_eta(), jet->get_phi());
0289       if (!m_inPPMode)
0290       {
0291         m_hSubSeedEnergyVsCent->Fill(cent_mbd, jet->get_e());
0292       }
0293       n_seed_sub++;
0294     }
0295   }
0296   m_hSubSeedCount->Fill(n_seed_sub);
0297   if (!m_inPPMode)
0298   {
0299     m_hSubSeedVsCent->Fill(cent_mbd, n_seed_sub);
0300   }
0301   return Fun4AllReturnCodes::EVENT_OK;
0302 }
0303
0304 //____________________________________________________________________________..
0305 int JetSeedCount::End(PHCompositeNode * /*topNode*/)
0306 {
0307   if (Verbosity() > 1)
0308   {
0309     std::cout << "JetSeedCount::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0310   }
0311   if (m_writeToOutputFile)
0312   {
0313     PHTFileServer::cd(m_outputFileName);
0314     m_hRawSeedCount->Write();
0315     m_hRawPt->Write();
0316     m_hRawPt_All->Write();
0317     m_hRawEtaVsPhi->Write();
0318     m_hSubSeedCount->Write();
0319     m_hSubPt->Write();
0320     m_hSubPt_All->Write();
0321     m_hSubEtaVsPhi->Write();
0322     if (!m_inPPMode)
0323     {
0324       m_hRawSeedEnergyVsCent->Write();
0325       m_hSubSeedEnergyVsCent->Write();
0326       m_hCentMbd->Write();
0327       m_hRawSeedVsCent->Write();
0328       m_hSubSeedVsCent->Write();
0329     }
0330   }
0331   return Fun4AllReturnCodes::EVENT_OK;
0332 }
0333
0334 //____________________________________________________________________________..
0335 void JetSeedCount::Print(const std::string &what) const
0336 {
0337   std::cout << "JetSeedCount::Print(const std::string &what) const Printing info for " << what << std::endl;
0338 }