JetUnderlyingEvent/src/JetUnderlyingEvent.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 JetUnderlyingEvent.h.
0007 //
0008 // JetUnderlyingEvent(const std::string &name = "JetUnderlyingEvent")
0009 // everything is keyed to JetUnderlyingEvent, 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 // JetUnderlyingEvent::~JetUnderlyingEvent()
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 JetUnderlyingEvent::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 JetUnderlyingEvent::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 JetUnderlyingEvent::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 JetUnderlyingEvent::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 JetUnderlyingEvent::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 JetUnderlyingEvent::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 JetUnderlyingEvent::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void JetUnderlyingEvent::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 "JetUnderlyingEvent.h"
0065
0066 #include <TH3.h>
0067 #include <fun4all/Fun4AllReturnCodes.h>
0068
0069 #include <phool/PHCompositeNode.h>
0070
0071 #include <fun4all/Fun4AllReturnCodes.h>
0072 #include <fun4all/PHTFileServer.h>
0073
0074 #include <phool/PHCompositeNode.h>
0075 #include <phool/getClass.h>
0076
0077 #include <g4jets/JetMap.h>
0078 #include <g4jets/Jetv1.h>
0079
0080 #include <centrality/CentralityInfo.h>
0081 #include <TFile.h>
0082 #include <TH3F.h>
0083 #include <TH2.h>
0084 #include <TH2F.h>
0085 #include <TString.h>
0086
0087 #include <calobase/RawTower.h>
0088 #include <calobase/RawTowerContainer.h>
0089 #include <calobase/RawTowerGeom.h>
0090 #include <calobase/RawTowerGeomContainer.h>
0091
0092 #include <jetbackground/TowerBackground.h>
0093
0094
0095
0096
0097 //____________________________________________________________________________..
0098 JetUnderlyingEvent::JetUnderlyingEvent(const std::string& name, const std::string& outputfilename):
0099  SubsysReco(name)
0100  , m_outputFileName(outputfilename)
0101 {
0102   std::cout << "JetUnderlyingEvent::JetUnderlyingEvent(const std::string &name) Calling ctor" << std::endl;
0103 }
0104
0105 //____________________________________________________________________________..
0106 JetUnderlyingEvent::~JetUnderlyingEvent()
0107 {
0108   std::cout << "JetUnderlyingEvent::~JetUnderlyingEvent() Calling dtor" << std::endl;
0109 }
0110
0111 //____________________________________________________________________________..
0112 int JetUnderlyingEvent::Init(PHCompositeNode *topNode)
0113 {
0114   std::cout << "JetUnderlyingEvent::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0115   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0116   std::cout << "JetUnderlyingEvent::Init - Output to " << m_outputFileName << std::endl;
0117
0118   // Histograms
0119
0120   hsubtractedE = new TH3F("HsubtractedE", "", 100, 0.0, 100, 100, 0, 100, 10, 0, 100);
0121    hsubtractedE->GetXaxis()->SetTitle("E_{T} [GeV]");
0122    hsubtractedE->GetYaxis()->SetTitle("Subtracted E_{T} [GeV]");
0123    hsubtractedE->GetZaxis()->SetTitle("Centrality");
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126
0127 //____________________________________________________________________________..
0128 int JetUnderlyingEvent::InitRun(PHCompositeNode *topNode)
0129 {
0130   std::cout << "JetUnderlyingEvent::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0131   return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133
0134 //____________________________________________________________________________..
0135 int JetUnderlyingEvent::process_event(PHCompositeNode *topNode)
0136 {
0137   // interface to reco jets
0138   JetMap* jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_r04_Sub1");
0139   if (!jets){
0140       std::cout
0141     << "MyJetAnalysis::process_event - Error can not find DST Reco JetMap node "
0142         << std::endl;
0143       exit(-1);
0144     }
0145
0146
0147   //calorimeter towers
0148   RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0149   RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0150   RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0151   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0152   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0153   if(!towersEM3 || !towersIH3 || !towersOH3){
0154     std::cout
0155       <<"MyJetAnalysis::process_event - Error can not find raw tower node "
0156       << std::endl;
0157     exit(-1);
0158   }
0159
0160   if(!tower_geom || !tower_geomOH){
0161     std::cout
0162       <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
0163       << std::endl;
0164     exit(-1);
0165   }
0166
0167   //centrality
0168
0169   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0170   if (!cent_node)
0171     {
0172       std::cout
0173         << "MyJetAnalysis::process_event - Error can not find centrality node "
0174     << std::endl;
0175       exit(-1);
0176     }
0177
0178   int m_centrality =  cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0179
0180   //underlying event
0181   TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub2");
0182   if(!background){
0183     std::cout<<"Can't get background. Exiting"<<std::endl;
0184     return Fun4AllReturnCodes::EVENT_OK;
0185   }
0186
0187
0188   //get reco jets
0189   float background_v2 = 0;
0190   float background_Psi2 = 0;
0191
0192     {
0193       background_v2 = background->get_v2();
0194       background_Psi2 = background->get_Psi2();
0195     }
0196   for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0197     {
0198
0199       Jet* jet = iter->second;
0200
0201
0202
0203       Jet* unsubjet = new Jetv1();
0204       float totalPx = 0;
0205       float totalPy = 0;
0206       float totalPz = 0;
0207       float totalE = 0;
0208       int nconst = 0;
0209
0210       for (Jet::ConstIter comp = jet->begin_comp(); comp != jet->end_comp(); ++comp)
0211     {
0212       RawTower *tower;
0213       nconst++;
0214
0215       if ((*comp).first == 15)
0216         {
0217           tower = towersIH3->getTower((*comp).second);
0218           if(!tower || !tower_geom){
0219         continue;
0220           }
0221           float UE = background->get_UE(1).at(tower->get_bineta());
0222           float tower_phi = tower_geom->get_tower_geometry(tower->get_key())->get_phi();
0223           float tower_eta = tower_geom->get_tower_geometry(tower->get_key())->get_eta();
0224           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0225           totalE += tower->get_energy() + UE;
0226           double pt = tower->get_energy() / cosh(tower_eta);
0227           totalPx += pt * cos(tower_phi);
0228           totalPy += pt * sin(tower_phi);
0229           totalPz += pt * sinh(tower_eta);
0230         }
0231       else if ((*comp).first == 16)
0232         {
0233           tower = towersOH3->getTower((*comp).second);
0234           if(!tower || !tower_geomOH)
0235         {
0236           continue;
0237         }
0238
0239           float UE = background->get_UE(2).at(tower->get_bineta());
0240           float tower_phi = tower_geomOH->get_tower_geometry(tower->get_key())->get_phi();
0241           float tower_eta = tower_geomOH->get_tower_geometry(tower->get_key())->get_eta();
0242           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0243           totalE +=tower->get_energy() + UE;
0244           double pt = tower->get_energy() / cosh(tower_eta);
0245           totalPx += pt * cos(tower_phi);
0246           totalPy += pt * sin(tower_phi);
0247           totalPz += pt * sinh(tower_eta);
0248         }
0249       else if ((*comp).first == 14)
0250         {
0251           tower = towersEM3->getTower((*comp).second);
0252           if(!tower || !tower_geom)
0253         {
0254           continue;
0255         }
0256           float UE = background->get_UE(0).at(tower->get_bineta());
0257           float tower_phi = tower_geom->get_tower_geometry(tower->get_key())->get_phi();
0258           float tower_eta = tower_geom->get_tower_geometry(tower->get_key())->get_eta();
0259           UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
0260           totalE +=tower->get_energy() + UE;
0261           double pt = tower->get_energy() / cosh(tower_eta);
0262           totalPx += pt * cos(tower_phi);
0263           totalPy += pt * sin(tower_phi);
0264           totalPz += pt * sinh(tower_eta);
0265         }
0266     }
0267
0268       //get unsubtracted jet
0269       unsubjet->set_px(totalPx);
0270       unsubjet->set_py(totalPy);
0271       unsubjet->set_pz(totalPz);
0272       unsubjet->set_e(totalE);
0273
0274       // fill histograms
0275       assert(hsubtractedE);
0276       hsubtractedE->Fill(jet->get_pt(),unsubjet->get_pt() - jet->get_pt(),m_centrality);
0277
0278
0279     }
0280
0281
0282
0283
0284    return Fun4AllReturnCodes::EVENT_OK;
0285 }
0286
0287 //____________________________________________________________________________..
0288  int JetUnderlyingEvent::ResetEvent(PHCompositeNode *topNode)
0289  {
0290    return Fun4AllReturnCodes::EVENT_OK;
0291  }
0292
0293 //____________________________________________________________________________..
0294  int JetUnderlyingEvent::EndRun(const int runnumber)
0295  {
0296   std::cout << "JetUnderlyingEvent::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0297    return Fun4AllReturnCodes::EVENT_OK;
0298  }
0299
0300 //____________________________________________________________________________..
0301 int JetUnderlyingEvent::End(PHCompositeNode *topNode)
0302 {
0303   std::cout << "JetValidation::End - Output to " << m_outputFileName << std::endl;
0304   PHTFileServer::get().cd(m_outputFileName);
0305
0306   TH2 *Centrality;
0307   for(int i = 0; i < hsubtractedE ->GetNbinsZ(); i++)
0308     {
0309       hsubtractedE -> GetZaxis()->SetRange(i+1,i+1);
0310       Centrality = (TH2*)  hsubtractedE -> Project3D("yx");
0311       Centrality -> Write(Form("hsubtractedE%1.0f",hsubtractedE->GetZaxis()->GetBinLowEdge(i+1)));
0312     }
0313
0314   return Fun4AllReturnCodes::EVENT_OK;
0315 }
0316
0317 //____________________________________________________________________________..
0318  int JetUnderlyingEvent::Reset(PHCompositeNode *topNode)
0319  {
0320   std::cout << "JetUnderlyingEvent::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0321   return Fun4AllReturnCodes::EVENT_OK;
0322  }
0323
0324 //____________________________________________________________________________..
0325  void JetUnderlyingEvent::Print(const std::string &what) const
0326  {
0327   std::cout << "JetUnderlyingEvent::Print(const std::string &what) const Printing info for " << what << std::endl;
0328  }