JetUESize-v2Psi2Centrality/src/UEvsEtaCentrality.cc

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in UEvsEtaCentrality.h.
0007 //
0008 // UEvsEtaCentrality(const std::string &name = "UEvsEtaCentrality")
0009 // everything is keyed to UEvsEtaCentrality, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // UEvsEtaCentrality::~UEvsEtaCentrality()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int UEvsEtaCentrality::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int UEvsEtaCentrality::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int UEvsEtaCentrality::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT;
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT;
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int UEvsEtaCentrality::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int UEvsEtaCentrality::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int UEvsEtaCentrality::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int UEvsEtaCentrality::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void UEvsEtaCentrality::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 //____________________________________________________________________________..
0063
0064 #include "UEvsEtaCentrality.h"
0065
0066 #include <fun4all/Fun4AllReturnCodes.h>
0067 #include <fun4all/PHTFileServer.h>
0068
0069 #include <phool/PHCompositeNode.h>
0070
0071 #include <phool/PHCompositeNode.h>
0072 #include <phool/getClass.h>
0073
0074 #include <g4jets/JetMap.h>
0075 #include <g4jets/Jetv1.h>
0076
0077 #include <centrality/CentralityInfo.h>
0078
0079 #include <calobase/RawTower.h>
0080 #include <calobase/RawTowerContainer.h>
0081 #include <calobase/RawTowerGeom.h>
0082 #include <calobase/RawTowerGeomContainer.h>
0083 #include <calobase/TowerInfoContainer.h>
0084 #include <calobase/TowerInfo.h>
0085
0086 #include <jetbackground/TowerBackground.h>
0087
0088 #include <TTree.h>
0089 #include <TH1F.h>
0090 #include <TH2F.h>
0091 //____________________________________________________________________________..
0092 UEvsEtaCentrality::UEvsEtaCentrality(const std::string& name, const std::string& outputfilename)
0093 : SubsysReco(name)
0094   , m_outputFileName(outputfilename)
0095 {
0096   std::cout << "UEvsEtaCentrality::UEvsEtaCentrality(const std::string &name) Calling ctor" << std::endl;
0097 }
0098
0099 //____________________________________________________________________________..
0100 UEvsEtaCentrality::~UEvsEtaCentrality()
0101 {
0102   std::cout << "UEvsEtaCentrality::~UEvsEtaCentrality() Calling dtor" << std::endl;
0103 }
0104
0105 //____________________________________________________________________________..
0106 int UEvsEtaCentrality::Init(PHCompositeNode *topNode)
0107 {
0108   std::cout << "UEvsEtaCentrality::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0109
0110   if (Verbosity() >=  UEvsEtaCentrality::VERBOSITY_SOME)
0111     std::cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << std::endl;
0112
0113   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0114
0115   hv2_cent = new TH2F("hv2_cent","",10,0,100,50,0,0.5);
0116   hPsi2_cent = new TH2F("hPsi2_cent","",10,0,100,50,-1.57,1.57);
0117   hUEiHcalEta_Cent0_20 = new TH2F("hUEiHcalEta_Cent0_20", "", 24, -1.1,1.1, 48, 0, 0.25);
0118   hUEoHcalEta_Cent0_20 = new TH2F("hUEoHcalEta_Cent0_20", "", 24, -1.1,1.1, 48, 0, 0.5);
0119   hUEemcalEta_Cent0_20 = new TH2F("hUEemcalEta_Cent0_20", "", 24, -1.1,1.1, 48, 0, 1.5);
0120
0121   hUEiHcalEta_Cent20_50 = new TH2F("hUEiHcalEta_Cent20_50", "", 24, -1.1,1.1, 48, 0, 0.25);
0122   hUEoHcalEta_Cent20_50 = new TH2F("hUEoHcalEta_Cent20_50", "", 24, -1.1,1.1, 48, 0, 0.5);
0123   hUEemcalEta_Cent20_50 = new TH2F("hUEemcalEta_Cent20_50", "", 24, -1.1,1.1, 48, 0, 1.5);
0124
0125   hUEiHcalEta_Cent50_100 = new TH2F("hUEiHcalEta_Cent50_100", "", 24, -1.1,1.1, 48, 0, 0.25);
0126   hUEoHcalEta_Cent50_100 = new TH2F("hUEoHcalEta_Cent50_100", "", 24, -1.1,1.1, 48, 0, 0.5);
0127   hUEemcalEta_Cent50_100 = new TH2F("hUEemcalEta_Cent50_100", "", 24, -1.1,1.1, 48, 0, 1.5);
0128
0129   return Fun4AllReturnCodes::EVENT_OK;
0130 }
0131
0132 //____________________________________________________________________________..
0133 int UEvsEtaCentrality::InitRun(PHCompositeNode *topNode)
0134 {
0135   std::cout << "UEvsEtaCentrality::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0136   return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138
0139 //____________________________________________________________________________..
0140 int UEvsEtaCentrality::process_event(PHCompositeNode *topNode)
0141 {
0142
0143   //centrality
0144   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0145   if (!cent_node)
0146     {
0147       std::cout
0148         << "MyJetAnalysis::process_event - Error can not find centrality node "
0149         << std::endl;
0150       exit(-1);
0151     }
0152
0153
0154
0155   //calorimeter towers
0156   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0157   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0158   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0159   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0160   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0161   if(!towersEM3 || !towersIH3 || !towersOH3){
0162     std::cout
0163       <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0164       << std::endl;
0165     exit(-1);
0166   }
0167
0168   if(!tower_geom || !tower_geomOH){
0169     std::cout
0170       <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0171       << std::endl;
0172     exit(-1);
0173   }
0174
0175   //underlying event
0176   TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0177   if(!background){
0178     std::cout<<"Can't get background. Exiting"<<std::endl;
0179     return Fun4AllReturnCodes::EVENT_OK;
0180   }
0181
0182     float background_v2 = 0;
0183     float background_Psi2 = 0;
0184     float m_centrality = 0;
0185
0186     m_centrality =  cent_node->get_centile(CentralityInfo::PROP::epd_NS);
0187     background_v2 = background->get_v2();
0188     background_Psi2 = background->get_Psi2();
0189
0190     hv2_cent->Fill(m_centrality, background_v2);
0191     hPsi2_cent->Fill(m_centrality, background_Psi2);
0192
0193     int i=0;
0194     for (i=0;i<=23;i++){
0195
0196
0197
0198     float UEi = background->get_UE(1).at(i); //24 eta bins for each detector Inner Hcal
0199     float UEo = background->get_UE(2).at(i); //24 eta bins for each detector, Outer Hcal
0200     float UEe = background->get_UE(0).at(i); //24 eta bins for each detector, EMCal
0201
0202     double eta = tower_geom->get_etacenter(i);
0203
0204     if (m_centrality > 0 && m_centrality <= 20){
0205     hUEiHcalEta_Cent0_20->Fill(eta, UEi);
0206     hUEoHcalEta_Cent0_20->Fill(eta, UEo);
0207     hUEemcalEta_Cent0_20->Fill(eta, UEe);
0208       }
0209     if (m_centrality > 20 && m_centrality <= 50){
0210     hUEiHcalEta_Cent20_50->Fill(eta, UEi);
0211     hUEoHcalEta_Cent20_50->Fill(eta, UEo);
0212     hUEemcalEta_Cent20_50->Fill(eta, UEe);
0213     }
0214
0215     if (m_centrality > 50 && m_centrality <= 100){
0216     hUEiHcalEta_Cent50_100->Fill(eta, UEi);
0217     hUEoHcalEta_Cent50_100->Fill(eta, UEo);
0218     hUEemcalEta_Cent50_100->Fill(eta, UEe);
0219     }
0220
0221
0222     }
0223
0224
0225
0226   // std::cout << "UEvsEtaCentrality::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229
0230 //____________________________________________________________________________..
0231 int UEvsEtaCentrality::ResetEvent(PHCompositeNode *topNode)
0232 {
0233   // std::cout << "UEvsEtaCentrality::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0234
0235   // IMPORTANT!! CLEAR YOUR VECTORS AND RESET YOUR TREE VARIABLES HERE!!
0236
0237   return Fun4AllReturnCodes::EVENT_OK;
0238 }
0239
0240 //____________________________________________________________________________..
0241 int UEvsEtaCentrality::EndRun(const int runnumber)
0242 {
0243   std::cout << "UEvsEtaCentrality::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0244   return Fun4AllReturnCodes::EVENT_OK;
0245 }
0246
0247 //____________________________________________________________________________..
0248 int UEvsEtaCentrality::End(PHCompositeNode *topNode)
0249 {
0250   std::cout << "MyJetAnalysis::End - Output to " << m_outputFileName << std::endl;
0251   PHTFileServer::get().cd(m_outputFileName);
0252
0253   hv2_cent->Write();
0254   hPsi2_cent->Write();
0255   hUEiHcalEta_Cent0_20->Write();
0256   hUEoHcalEta_Cent0_20->Write();
0257   hUEemcalEta_Cent0_20->Write();
0258   hUEiHcalEta_Cent20_50->Write();
0259   hUEoHcalEta_Cent20_50->Write();
0260   hUEemcalEta_Cent20_50->Write();
0261   hUEiHcalEta_Cent50_100->Write();
0262   hUEoHcalEta_Cent50_100->Write();
0263   hUEemcalEta_Cent50_100->Write();
0264   std::cout << "UEvsEtaCentrality::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0265   return Fun4AllReturnCodes::EVENT_OK;
0266 }
0267
0268 //____________________________________________________________________________..
0269 int UEvsEtaCentrality::Reset(PHCompositeNode *topNode)
0270 {
0271  std::cout << "UEvsEtaCentrality::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0272   return Fun4AllReturnCodes::EVENT_OK;
0273 }
0274
0275 //____________________________________________________________________________..
0276 void UEvsEtaCentrality::Print(const std::string &what) const
0277 {
0278   std::cout << "UEvsEtaCentrality::Print(const std::string &what) const Printing info for " << what << std::endl;
0279 }