QA/Jet/PhotonJetsKinematics.cc

0001 //____________________________________________________________________________..
0002
0003 #include "PhotonJetsKinematics.h"
0004
0005 // calobase includes
0006 #include <calobase/RawCluster.h>
0007 #include <calobase/RawClusterContainer.h>
0008 #include <calobase/RawClusterUtility.h>
0009 #include <calobase/RawTower.h>
0010 #include <calobase/RawTowerContainer.h>
0011 #include <calobase/RawTowerGeom.h>
0012 #include <calobase/RawTowerGeomContainer.h>
0013 #include <calobase/TowerInfo.h>
0014 #include <calobase/TowerInfoContainer.h>
0015 #include <calobase/TowerInfoDefs.h>
0016
0017 // calotrigger includes
0018 #include <calotrigger/TriggerAnalyzer.h>
0019
0020 // fun4all includes
0021 #include <fun4all/Fun4AllHistoManager.h>
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023
0024 // phool includes
0025 #include <phool/getClass.h>
0026 #include <phool/PHCompositeNode.h>
0027
0028 // qautils includes
0029 #include <qautils/QAHistManagerDef.h>
0030
0031 // root includes
0032 #include <TH1.h>
0033 #include <TH2.h>
0034 #include <TStyle.h>
0035
0036 // c++ includes
0037 #include <algorithm>
0038 #include <cassert>
0039 #include <iostream>
0040 #include <vector>
0041 //____________________________________________________________________________..
0042
0043 PhotonJetsKinematics::PhotonJetsKinematics(const std::string &modulename, const std::string &inputnode, const std::string &histtag) :
0044   SubsysReco(modulename)
0045   , m_modulename(modulename)
0046   , m_inputnode(inputnode)
0047   , m_histtag(histtag)
0048   , m_trgToSelect(JetQADefs::GL1::MBDNSJet1)
0049   , m_doTrgSelect(false)
0050   , m_doOptHist(false)
0051 {
0052   if (Verbosity() > 1)
0053     {
0054       std::cout << "PhotonJetsKinematics::PhotonJetsKinematics(const std::string &name, const std::string&outputfilename) Calling ctor" << std::endl;
0055     }
0056 }
0057 //____________________________________________________________________________..
0058 PhotonJetsKinematics::~PhotonJetsKinematics()
0059 {
0060   if (Verbosity() > 1)
0061     {
0062       std::cout << "PhotonJetsKinematics::~PhotonJetsKinematics() Calling dtor" << std::endl;
0063     }
0064   delete m_analyzer;
0065 }
0066
0067 //____________________________________________________________________________..
0068 int PhotonJetsKinematics::Init(PHCompositeNode* /*topNode*/)
0069 {
0070   if (Verbosity() > 1)
0071     {
0072       std::cout << "PhotonJetsKinematics::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0073     }
0074
0075   // initialize trigger analyzer and hist manager
0076   delete m_analyzer;
0077   m_analyzer = new TriggerAnalyzer();
0078
0079   gStyle->SetOptTitle(0);
0080   m_manager = QAHistManagerDef::getHistoManager();
0081   if (!m_manager)
0082     {
0083       std::cerr << PHWHERE << "PANIC: couldn't grab histogram manager!" << std::endl;
0084       assert(m_manager);
0085     }
0086
0087   // make sure module name is lower case
0088   std::string smallModuleName = m_modulename;
0089   std::transform(
0090          smallModuleName.begin(),
0091          smallModuleName.end(),
0092          smallModuleName.begin(),
0093          ::tolower);
0094   // construct histogram names
0095   std::vector<std::string> vecHistNames = {
0096     "emcal_cluster_chi2",
0097     "emcal_cluster_energy",
0098     "emcal_cluster_eta_phi", // optional
0099     "emcal_cluster_eta", // optional
0100     "emcal_cluster_phi", // optional
0101     "emcal_cluster_energy_eta", //see what eta distribution looks like with different energy
0102     "emcal_cluster_chi2_eta", //see what eta distribution looks like with different chi2
0103     "emcal_cluster_eta_with_cuts",//chi2<3,Et>500 MeV
0104     "emcal_cluster_energy_phi",
0105     "emcal_cluster_chi2_phi",
0106     "emcal_cluster_phi_with_cuts",//chi2 < 3, Et > 500 MeV
0107     "emcal_cluster_eta_phi_with_cuts",//chi2 < 3, Et > 500 MeV
0108     "emcal_cluster_eta_with_energy_cut",//optional: Et > 500 MeV only
0109     "emcal_cluster_chi2_energy"
0110   };
0111
0112   for (auto& histName : vecHistNames)
0113     {
0114       histName.insert(0, "h_" + smallModuleName + "_");
0115       if (!m_histtag.empty())
0116     {
0117       histName.append("_" + m_histtag);
0118     }
0119     }
0120
0121   //initializing histograms
0122   if (m_doOptHist)
0123   {
0124     h_emcal_cluster_eta_phi = new TH2F(vecHistNames[2].data(), "", 48, -1.2, 1.2, 64, -M_PI, M_PI);
0125     h_emcal_cluster_eta_phi->GetXaxis()->SetTitle("#eta");
0126     h_emcal_cluster_eta_phi->GetYaxis()->SetTitle("#Phi");
0127     h_emcal_cluster_eta_phi->SetOption("COLZ");
0128
0129     h_emcal_cluster_eta = new TH1D(vecHistNames[3].data(), "", 48, -1.2, 1.2);
0130     h_emcal_cluster_eta->GetXaxis()->SetTitle("#eta");
0131
0132     h_emcal_cluster_phi = new TH1D(vecHistNames[4].data(), "", 64, -M_PI, M_PI);
0133     h_emcal_cluster_phi->GetXaxis()->SetTitle("#phi");
0134
0135     h_emcal_cluster_eta_with_energy_cut = new TH1D(vecHistNames[12].data(), "", 48, -1.2, 1.2);
0136     h_emcal_cluster_eta_with_energy_cut->GetXaxis()->SetTitle("#eta");
0137   }
0138   h_emcal_cluster_chi2 = new TH1D(vecHistNames[0].data(), "", 30, 0, 150);
0139   h_emcal_cluster_chi2->GetXaxis()->SetTitle("#chi^{2}");
0140
0141   h_emcal_cluster_energy = new TH1D(vecHistNames[1].data(), "", 100, 0, 10);
0142   h_emcal_cluster_energy->GetXaxis()->SetTitle("E_{T} [GeV]");
0143
0144   h_emcal_cluster_energy_eta = new TH2D(vecHistNames[5].data(), "",48,-1.2,1.2,100,0,10);
0145   h_emcal_cluster_energy_eta->GetXaxis()->SetTitle("#eta");
0146   h_emcal_cluster_energy_eta->GetYaxis()->SetTitle("E_{T} [GeV]");
0147
0148   h_emcal_cluster_chi2_eta = new TH2D(vecHistNames[6].data(), "",48,-1.2,1.2,150,0,150);
0149   h_emcal_cluster_chi2_eta->GetXaxis()->SetTitle("#eta");
0150   h_emcal_cluster_chi2_eta->GetYaxis()->SetTitle("#chi^{2}");
0151
0152   h_emcal_cluster_eta_with_cuts = new TH1D(vecHistNames[7].data(), "", 48, -1.2, 1.2);
0153   h_emcal_cluster_eta_with_cuts->GetXaxis()->SetTitle("#eta");
0154
0155   h_emcal_cluster_energy_phi = new TH2D(vecHistNames[8].data(), "",64,-M_PI,M_PI,100,0,10);
0156   h_emcal_cluster_energy_phi->GetXaxis()->SetTitle("#phi");
0157   h_emcal_cluster_energy_phi->GetYaxis()->SetTitle("E_{T} [GeV]");
0158
0159   h_emcal_cluster_chi2_phi = new TH2D(vecHistNames[9].data(), "",64,-M_PI,M_PI,150,0,150);
0160   h_emcal_cluster_chi2_phi->GetXaxis()->SetTitle("#phi");
0161   h_emcal_cluster_chi2_phi->GetYaxis()->SetTitle("#chi2^{2}");
0162
0163   h_emcal_cluster_phi_with_cuts = new TH1D(vecHistNames[10].data(), "", 64,-M_PI, M_PI);
0164   h_emcal_cluster_phi_with_cuts->GetXaxis()->SetTitle("#phi");
0165
0166   h_emcal_cluster_eta_phi_with_cuts = new TH2F(vecHistNames[11].data(), "", 48, -1.2, 1.2, 64, -M_PI, M_PI);
0167   h_emcal_cluster_eta_phi_with_cuts->GetXaxis()->SetTitle("#eta");
0168   h_emcal_cluster_eta_phi_with_cuts->GetYaxis()->SetTitle("#phi");
0169   h_emcal_cluster_eta_phi_with_cuts->SetOption("COLZ");
0170
0171   h_emcal_cluster_chi2_energy = new TH2F(vecHistNames[13].data(),"",100,0,10,150,0,150);
0172   h_emcal_cluster_chi2_energy->GetXaxis()->SetTitle("E_{T} [GeV]");
0173   h_emcal_cluster_chi2_energy->GetYaxis()->SetTitle("#chi2^{2}");
0174
0175   return Fun4AllReturnCodes::EVENT_OK;
0176
0177 }
0178
0179 //____________________________________________________________________________..
0180 int PhotonJetsKinematics::InitRun(PHCompositeNode* /*topNode*/)
0181 {
0182
0183   if (Verbosity() > 1)
0184     {
0185       std::cout << "PhotonJetsKinematics::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0186     }
0187   return Fun4AllReturnCodes::EVENT_OK;
0188 }
0189
0190 //____________________________________________________________________________..
0191 int PhotonJetsKinematics::process_event(PHCompositeNode *topNode)
0192 {
0193
0194   RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_inputnode);
0195   if (!clusterContainer)
0196     {
0197       std::cout << PHWHERE << "PhotonJetsKinematics::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0198       return 0;
0199     }
0200
0201   // if needed, check if trigger fired
0202   if (m_doTrgSelect)
0203     {
0204       m_analyzer->decodeTriggers(topNode);
0205       bool hasTrigger = JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer);
0206       if (!hasTrigger)
0207     {
0208       return Fun4AllReturnCodes::EVENT_OK;
0209     }
0210     }
0211
0212   RawClusterContainer::ConstIterator clusterIter;
0213   RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0214
0215   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0216     {
0217       RawCluster* recoCluster = clusterIter->second;
0218       CLHEP::Hep3Vector vertex(0, 0, 0);
0219       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0220
0221       float clusE = E_vec_cluster.mag();
0222       float clus_eta = E_vec_cluster.pseudoRapidity();
0223       float clus_phi = E_vec_cluster.phi();
0224       float clus_chisq = recoCluster->get_chi2();
0225       float clusEt = clusE/std::cosh(clus_eta);//this is what really wanted to be plotted on histogram related to energy
0226
0227       //Filling Histograms, only taking into account E vec cluster, not reco cluster
0228
0229        h_emcal_cluster_chi2->Fill(clus_chisq);
0230        h_emcal_cluster_energy->Fill(clusEt);
0231        if (m_doOptHist)
0232        {
0233          h_emcal_cluster_eta_phi->Fill(clus_eta, clus_phi);
0234          h_emcal_cluster_eta->Fill(clus_eta);  // 1D eta plot
0235          h_emcal_cluster_phi->Fill(clus_phi);  // 1D phi plot
0236        }
0237        h_emcal_cluster_energy_eta->Fill(clus_eta, clusEt);//2D energy eta plot
0238        h_emcal_cluster_chi2_eta->Fill(clus_eta,clus_chisq);//2D chi2 eta plot
0239        h_emcal_cluster_energy_phi->Fill(clus_phi, clusEt);//2D energy phi plot
0240        h_emcal_cluster_chi2_phi->Fill(clus_phi,clus_chisq);//2D chi2 phi plot
0241        h_emcal_cluster_chi2_energy->Fill(clusEt,clus_chisq);
0242        if(clus_chisq<3 && clusEt> 0.5){
0243        h_emcal_cluster_eta_with_cuts->Fill(clus_eta);
0244        h_emcal_cluster_phi_with_cuts->Fill(clus_phi);
0245        h_emcal_cluster_eta_phi_with_cuts->Fill(clus_eta, clus_phi);
0246        }
0247        if (m_doOptHist)
0248        {
0249          if (clusEt>0.5)
0250          {
0251            h_emcal_cluster_eta_with_energy_cut->Fill(clus_eta);
0252          }
0253        }
0254      }
0255
0256   if (Verbosity() > 1)
0257     {
0258       std::cout << "PhotonJetsKinematics::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0259     }
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263 //____________________________________________________________________________..
0264 int PhotonJetsKinematics::ResetEvent(PHCompositeNode* /*topNode*/)
0265 {
0266   if (Verbosity() > 1)
0267     {
0268       std::cout << "PhotonJetsKinematics::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0269     }
0270   return Fun4AllReturnCodes::EVENT_OK;
0271 }
0272
0273 //____________________________________________________________________________..
0274 int PhotonJetsKinematics::EndRun(const int runnumber)
0275 {
0276   if (Verbosity() > 1)
0277     {
0278       std::cout << "PhotonJetsKinematics::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0279     }
0280   return Fun4AllReturnCodes::EVENT_OK;
0281 }
0282
0283 //____________________________________________________________________________..
0284 int PhotonJetsKinematics::End(PHCompositeNode* /*topNode*/)
0285 {
0286   // Commenting out the following line because Jenkins keeps failing the build test :(
0287   // if (Verbosity() > 1) std::cout << "PhotonJetsKinematics::End - Output to " << outfilename << std::endl;
0288
0289   //Outputting the histograms
0290   if (m_doOptHist)
0291   {
0292     m_manager->registerHisto(h_emcal_cluster_eta_phi);
0293     m_manager->registerHisto(h_emcal_cluster_eta);
0294     m_manager->registerHisto(h_emcal_cluster_phi);
0295     m_manager->registerHisto(h_emcal_cluster_eta_with_energy_cut);
0296   }
0297   m_manager->registerHisto(h_emcal_cluster_energy);
0298   m_manager->registerHisto(h_emcal_cluster_energy_eta);
0299   m_manager->registerHisto(h_emcal_cluster_chi2);
0300   m_manager->registerHisto(h_emcal_cluster_chi2_eta);
0301   m_manager->registerHisto(h_emcal_cluster_eta_with_cuts);
0302   m_manager->registerHisto(h_emcal_cluster_phi_with_cuts);
0303   m_manager->registerHisto(h_emcal_cluster_energy_phi);
0304   m_manager->registerHisto(h_emcal_cluster_chi2_phi);
0305   m_manager->registerHisto(h_emcal_cluster_eta_phi_with_cuts);
0306   m_manager->registerHisto(h_emcal_cluster_chi2_energy);
0307
0308   if (Verbosity() > 1)
0309     {
0310       std::cout << "PhotonJetsKinematics::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0311     }
0312   return Fun4AllReturnCodes::EVENT_OK;
0313
0314 }
0315
0316 //____________________________________________________________________________..
0317 int PhotonJetsKinematics::Reset(PHCompositeNode* /*topNode*/)
0318 {
0319   if (Verbosity() > 1)
0320     {
0321       std::cout << "PhotonJetsKinematics::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0322     }
0323   return Fun4AllReturnCodes::EVENT_OK;
0324 }
0325
0326 //____________________________________________________________________________..
0327 void PhotonJetsKinematics::Print(const std::string &what) const
0328 {
0329   if (Verbosity() > 1)
0330     {
0331       std::cout << "PhotonJetsKinematics::Print(const std::string &what) const Printing info for " << what << std::endl;
0332     }
0333 }