JetSeedCount/src/JetSeedCount.cc

0001 #include <iostream>
0002 #include "JetSeedCount.h"
0003
0004 #include <fun4all/Fun4AllBase.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <fun4all/PHTFileServer.h>
0007
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/getClass.h>
0010
0011 #include <jetbase/JetMap.h>
0012 #include <jetbase/JetContainer.h>
0013 #include <jetbase/Jetv2.h>
0014 #include <jetbase/Jetv1.h>
0015
0016 #include <centrality/CentralityInfov2.h>
0017 #include <globalvertex/GlobalVertex.h>
0018 #include <globalvertex/GlobalVertexMap.h>
0019
0020 #include <ffaobjects/EventHeader.h>
0021 #include <ffaobjects/EventHeaderv1.h>
0022
0023 #include <calobase/RawTower.h>
0024 #include <calobase/RawTowerContainer.h>
0025 #include <calobase/RawTowerGeom.h>
0026 #include <calobase/RawTowerGeomContainer.h>
0027 #include <calobase/TowerInfoContainer.h>
0028 #include <calobase/TowerInfo.h>
0029
0030 #include <jetbackground/TowerBackground.h>
0031 #include <TCanvas.h>
0032 #include <TH2.h>
0033
0034 JetSeedCount::JetSeedCount(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename):
0035  SubsysReco("JetSeedCount_" + recojetname + "_" + truthjetname)
0036  , m_recoJetName(recojetname)
0037  , m_truthJetName(truthjetname)
0038  , m_outputFileName(outputfilename)
0039  , m_etaRange(-1,1)
0040  , m_ptRange(5, 100)
0041 {
0042     //std::cout << "JetSeedCount::JetSeedCount(const std::string &name) Calling ctor" << std::endl;
0043   }
0044
0045 JetSeedCount::~JetSeedCount()
0046 {
0047   std::cout << "JetSeedCount::~JetSeedCount() Calling dtor" << std::endl;
0048 }
0049
0050 int JetSeedCount::Init(PHCompositeNode *topNode)
0051 {
0052   std::cout << "JetSeedCount::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0053   std::cout << "Opening output file named " << m_outputFileName << std::endl;
0054   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0055   return Fun4AllReturnCodes::EVENT_OK;
0056 }
0057
0058 //____________________________________________________________________________..
0059 int JetSeedCount::InitRun(PHCompositeNode *topNode)
0060 {
0061   std::cout << "JetSeedCount::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0062   return Fun4AllReturnCodes::EVENT_OK;
0063 }
0064
0065 //____________________________________________________________________________..
0066 int JetSeedCount::process_event(PHCompositeNode *topNode)
0067 {
0068   //std::cout << "JetSeedCount::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0069   ++m_event;
0070   //Calling Raw Jet Seeds
0071   JetContainer* seedjetsraw = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0072   if (!seedjetsraw){
0073      std::cout << "JetSeedCount::process_event - Error can not find DST raw seed jets" << std::endl;
0074      exit(-1);
0075   }
0076
0077   //Calling Sub jet seeds
0078   JetContainer* seedjetssub = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
0079   if (!seedjetssub){
0080      std::cout << "JetSeedCount::process_event - Error can not find DST sub seed jets" << std::endl;
0081      exit(-1);
0082   }
0083
0084   //Calling Centrality Info
0085   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0086   if (!cent_node){
0087      std::cout << "JetSeedCount::process_event - Error can not find CentralityInfo" << std::endl;
0088      exit(-1);
0089   }
0090
0091   //Z vertex info
0092   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0093   if (!vertexmap){
0094      std::cout
0095      << "MyJetAnalysis::process_event - Error can not find global vertex  node "
0096      << std::endl;
0097      exit(-1);
0098    }
0099    if (vertexmap->empty()){
0100       std::cout
0101       << "MyJetAnalysis::process_event - global vertex node is empty "
0102       << std::endl;
0103       return Fun4AllReturnCodes::ABORTEVENT;
0104     }
0105
0106    GlobalVertex *vtx = vertexmap->begin()->second;
0107    z_vtx = vtx->get_z();
0108
0109   //Saving Centrailty info
0110   float cent_mbd = cent_node->get_centile(CentralityInfo::PROP::bimp);
0111   m_centrality.push_back(cent_mbd);
0112
0113   //Raw Seed Count
0114   m_seed_raw = 0;
0115   float Counter = 0;
0116   //for (JetMap::Iter iter = seedjetsraw->begin(); iter != seedjetsraw->end(); ++iter){
0117   for(auto jet : *seedjetsraw){
0118      //Jet* jet = iter->second;
0119      int passesCut = jet->get_property(seedjetsraw->property_index(Jet::PROPERTY::prop_SeedItr));
0120      Counter += 1;
0121      m_rawpt_all.push_back(jet->get_pt());
0122      if (passesCut == 1){
0123         m_rawpt.push_back(jet->get_pt());
0124         m_rawenergy.push_back(jet->get_e());
0125         m_RawEta.push_back(jet->get_eta());
0126         m_RawPhi.push_back(jet->get_phi());
0127         m_rawcent.push_back(cent_mbd);
0128         m_seed_raw++;
0129      }
0130   }
0131   m_raw_counts.push_back(m_seed_raw);
0132
0133   //Sub Seed Count
0134   m_seed_sub = 0;
0135   Counter = 0;
0136   //for (unsigned int iter = 0; iter < seedjetssub->size(); ++iter){
0137   //Jet* jet = seedjetsub->get_jet(iter);
0138   for(auto jet : *seedjetssub){
0139      //Jet* jet = iter->second;
0140      int passesCut = jet->get_property(seedjetssub->property_index(Jet::PROPERTY::prop_SeedItr));
0141      Counter += 1;
0142      m_subpt_all.push_back(jet->get_pt());
0143      if (passesCut == 2){
0144         m_subpt.push_back(jet->get_pt());
0145         m_subenergy.push_back(jet->get_e());
0146         m_SubEta.push_back(jet->get_eta());
0147         m_SubPhi.push_back(jet->get_phi());
0148         m_subcent.push_back(cent_mbd);
0149         m_seed_sub++;
0150      }
0151   }
0152   m_sub_counts.push_back(m_seed_sub);
0153   return Fun4AllReturnCodes::EVENT_OK;
0154 }
0155
0156 //____________________________________________________________________________..
0157 int JetSeedCount::ResetEvent(PHCompositeNode *topNode)
0158 {
0159   //std::cout << "JetSeedCount::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0160   return Fun4AllReturnCodes::EVENT_OK;
0161 }
0162
0163 //____________________________________________________________________________..
0164 int JetSeedCount::EndRun(const int runnumber)
0165 {
0166   std::cout << "JetSeedCount::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0167   return Fun4AllReturnCodes::EVENT_OK;
0168 }
0169
0170 //____________________________________________________________________________..
0171 int JetSeedCount::End(PHCompositeNode *topNode)
0172 {
0173   std::cout << "JetSeedCount::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0174   PHTFileServer::get().cd(m_outputFileName);
0175
0176   TH1F *hRawSeedCount = new TH1F("hRawSeedCount", "Raw Seed Count per Event", 100, 0.00, 50.00);
0177   hRawSeedCount->GetXaxis()->SetTitle("Raw Seed Count per Event");
0178   hRawSeedCount->GetYaxis()->SetTitle("Number of Entries");
0179   for (int j = 0; j < (int)m_raw_counts.size(); j++){
0180      hRawSeedCount->Fill(m_raw_counts.at(j));
0181   }
0182
0183   TH1F *hRawPt = new TH1F("hRawPt", "Raw p_{T}", 1000, 0.00, 50.00);
0184   hRawPt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0185   hRawPt->GetYaxis()->SetTitle("Number of Entries");
0186   for (int j = 0; j < (int)m_rawpt.size(); j++){
0187      hRawPt->Fill(m_rawpt.at(j));
0188   }
0189
0190   TH1F *hRawPt_All = new TH1F("hRawPt_All", "Raw p_{T} (all jet seeds)", 1000, 0.00, 50.00);
0191   hRawPt_All->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0192   hRawPt_All->GetYaxis()->SetTitle("Number of Entries");
0193   for (int j = 0; j < (int)m_rawpt_all.size(); j++){
0194      hRawPt_All->Fill(m_rawpt_all.at(j));
0195   }
0196
0197   TH2F *hRawEtaVsPhi = new TH2F("hRawEtaVsPhi", "Raw Seed Eta Vs Phi", 220, -1.1, 1.1, 628, -3.14, 3.14);
0198   hRawEtaVsPhi->GetXaxis()->SetTitle("Jet #eta [Rads.]");
0199   hRawEtaVsPhi->GetYaxis()->SetTitle("Jet #phi [Rads.]");
0200   for (int j = 0; j < (int)m_RawEta.size(); j++){
0201      hRawEtaVsPhi->Fill(m_RawEta.at(j), m_RawPhi.at(j));
0202   }
0203
0204   TH1F *hSubSeedCount = new TH1F("hSubSeedCount", "Sub Seed Count per Event", 100, 0.00, 50.00);
0205   hSubSeedCount->GetXaxis()->SetTitle("Sub Seed Count per Event");
0206   hSubSeedCount->GetYaxis()->SetTitle("Number of Entries");
0207   for (int j = 0; j < (int)m_sub_counts.size(); j++){
0208      hSubSeedCount->Fill(m_sub_counts.at(j));
0209   }
0210
0211   TH1F *hSubPt = new TH1F("hSubPt", "Sub. p_{T}", 1000, 0.00, 50.00);
0212   hSubPt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0213   hSubPt->GetYaxis()->SetTitle("Number of Entries");
0214   for (int j = 0; j < (int)m_subpt.size(); j++){
0215      hSubPt->Fill(m_subpt.at(j));
0216   }
0217
0218   TH1F *hSubPt_All = new TH1F("hSubPt_All", "Sub. p_{T} (all jet seeds)", 1000, 0.00, 50.00);
0219   hSubPt_All->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0220   hSubPt_All->GetYaxis()->SetTitle("Number of Entries");
0221   for (int j = 0; j < (int)m_subpt_all.size(); j++){
0222      hSubPt->Fill(m_subpt_all.at(j));
0223   }
0224
0225   TH2F *hSubEtaVsPhi = new TH2F("hSubEtaVsPhi", "Sub. Seed Eta Vs Phi", 220, -1.1, 1.1, 628, -3.14, 3.14);
0226   hSubEtaVsPhi->GetXaxis()->SetTitle("Jet #eta [Rads.]");
0227   hSubEtaVsPhi->GetYaxis()->SetTitle("Jet #phi [Rads.]");
0228   for (int j = 0; j < (int)m_SubEta.size(); j++){
0229      hSubEtaVsPhi->Fill(m_SubEta.at(j), m_SubPhi.at(j));
0230   }
0231
0232   TH2F *hRawSeedEnergyVsCent = new TH2F("hRawSeedEnergyVsCent", "Raw Seed Energy Vs Centrality", 10.00, 0.00, 100.00, 100, 0.00, 50.00);
0233   hRawSeedEnergyVsCent->GetXaxis()->SetTitle("Centrality");
0234   hRawSeedEnergyVsCent->GetYaxis()->SetTitle("RawSeedEnergy");
0235   for (int j = 0; j < (int)m_rawenergy.size(); j++){
0236      hRawSeedEnergyVsCent->Fill(m_rawcent.at(j), m_rawenergy.at(j));
0237   }
0238
0239   TH2F *hSubSeedEnergyVsCent = new TH2F("hSubSeedEnergyVsCent", "Sub Seed Energy Vs Centrality", 10.00, 0.00, 100.00, 100, 0.00, 50.00);
0240   hSubSeedEnergyVsCent->GetXaxis()->SetTitle("Centrality");
0241   hSubSeedEnergyVsCent->GetYaxis()->SetTitle("SubSeedEnergy");
0242   for (int j = 0; j < (int)m_subenergy.size(); j++){
0243      hSubSeedEnergyVsCent->Fill(m_subcent.at(j), m_subenergy.at(j));
0244   }
0245
0246    TH1F *hCentMbd = new TH1F("hCentMbd", "hCentMbd", 10, 0.00, 100.00);
0247   hCentMbd->GetXaxis()->SetTitle("Centrality (Mbd)");
0248   hCentMbd->GetYaxis()->SetTitle("Number of Entries");
0249   for (int j = 0; j < (int)m_centrality.size(); j++){
0250      hCentMbd->Fill(m_centrality.at(j));
0251   }
0252
0253   TH2F *hRawSeedVsCent = new TH2F("hRawSeedVsCent", "Raw Seed Vs Centrality", 10, 0.00, 100.00, 101, -0.5, 100.5);
0254   hRawSeedVsCent->GetXaxis()->SetTitle("Centrality");
0255   hRawSeedVsCent->GetYaxis()->SetTitle("Raw Seed Count");
0256   for (int j = 0; j < (int)m_raw_counts.size(); j++){
0257      hRawSeedVsCent->Fill(m_centrality.at(j), m_raw_counts.at(j));
0258   }
0259
0260    TH2F *hSubSeedVsCent = new TH2F("hSubSeedVsCent", "Sub Seed Vs Centrality", 10, 0.00, 100.00, 101, -0.5, 100.5);
0261   hSubSeedVsCent->GetXaxis()->SetTitle("Centrality");
0262   hSubSeedVsCent->GetYaxis()->SetTitle("Sub Seed Count");
0263   for (int j = 0; j < (int)m_sub_counts.size(); j++){
0264      hSubSeedVsCent->Fill(m_centrality.at(j), m_sub_counts.at(j));
0265   }
0266
0267   hRawPt->Write();
0268   hSubPt->Write();
0269   hRawPt_All->Write();
0270   hSubPt_All->Write();
0271   hRawEtaVsPhi->Write();
0272   hSubEtaVsPhi->Write();
0273   hRawSeedCount->Write();
0274   hSubSeedCount->Write();
0275   hRawSeedEnergyVsCent->Write();
0276   hSubSeedEnergyVsCent->Write();
0277   hCentMbd->Write();
0278   hRawSeedVsCent->Write();
0279   hSubSeedVsCent->Write();
0280
0281   return Fun4AllReturnCodes::EVENT_OK;
0282 }
0283
0284 //____________________________________________________________________________..
0285 int JetSeedCount::Reset(PHCompositeNode *topNode)
0286 {
0287  std::cout << "JetSeedCount::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0288  return Fun4AllReturnCodes::EVENT_OK;
0289 }
0290
0291 //____________________________________________________________________________..
0292 void JetSeedCount::Print(const std::string &what) const
0293 {
0294   std::cout << "JetSeedCount::Print(const std::string &what) const Printing info for " << what << std::endl;
0295 }