Jet-Study/src/EventValidation.cc

0001 #include <EventValidation.h>
0002 // -- c++
0003 #include <cmath>
0004 #include <iomanip>
0005 #include <iostream>
0006 #include <sstream>
0007 #include <fstream>
0008 // -- event
0009 #include <ffaobjects/EventHeader.h>
0010 #include <fun4all/Fun4AllServer.h>
0011 // -- fun4all
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 // -- vertex
0014 #include <globalvertex/GlobalVertex.h>
0015 #include <globalvertex/GlobalVertexMap.h>
0016 // -- jet
0017 #include <jetbase/JetContainer.h>
0018 // -- DST node
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/getClass.h>
0021 // -- root
0022 #include <TFile.h>
0023 // -- Tower stuff
0024 #include <calobase/RawTowerDefs.h>
0025 #include <calobase/TowerInfoContainer.h>
0026 #include <calobase/TowerInfo.h>
0027 #include <calobase/RawTowerGeomContainer.h>
0028 #include <calobase/RawTowerGeom.h>
0029
0030 // Jet Utils
0031 #include "JetUtils.h"
0032
0033 using std::cout;
0034 using std::endl;
0035 using std::to_string;
0036 using std::make_pair;
0037 using std::stringstream;
0038 using std::ofstream;
0039
0040 //____________________________________________________________________________..
0041 EventValidation::EventValidation()
0042   : SubsysReco("EventValidation")
0043   , m_emcTowerNode("TOWERINFO_CALIB_CEMC_RETOWER")
0044   , m_ihcTowerNode("TOWERINFO_CALIB_HCALIN")
0045   , m_ohcTowerNode("TOWERINFO_CALIB_HCALOUT")
0046   , m_emcGeomNode("TOWERGEOM_CEMC")
0047   , m_ihcGeomNode("TOWERGEOM_HCALIN")
0048   , m_ohcGeomNode("TOWERGEOM_HCALOUT")
0049   , m_recoJetName_r04("AntiKt_unsubtracted_r04")
0050   , m_bins_phi(64)
0051   , m_phi_low(-M_PI)
0052   , m_phi_high(M_PI)
0053   , m_bins_eta(24)
0054   , m_eta_low(-1.1)
0055   , m_eta_high(1.1)
0056   , m_R(0.4)
0057 {
0058   cout << "EventValidation::EventValidation(const std::string &name) Calling ctor" << endl;
0059 }
0060
0061 //____________________________________________________________________________..
0062 EventValidation::~EventValidation()
0063 {
0064   cout << "EventValidation::~EventValidation() Calling dtor" << endl;
0065 }
0066
0067 //____________________________________________________________________________..
0068 Int_t EventValidation::Init(PHCompositeNode *topNode)
0069 {
0070   cout << "EventValidation::Init(PHCompositeNode *topNode) Initializing" << endl;
0071
0072   // so that the histos actually get written out
0073   Fun4AllServer *se = Fun4AllServer::instance();
0074   se->Print("NODETREE");
0075
0076   return Fun4AllReturnCodes::EVENT_OK;
0077 }
0078
0079 //____________________________________________________________________________..
0080 Int_t EventValidation::process_event(PHCompositeNode *topNode)
0081 {
0082   EventHeader* eventInfo = findNode::getClass<EventHeader>(topNode,"EventHeader");
0083   if(!eventInfo)
0084   {
0085     cout << PHWHERE << "EventValidation::process_event - Fatal Error - EventHeader node is missing. " << endl;
0086     return Fun4AllReturnCodes::ABORTEVENT;
0087   }
0088
0089   int m_globalEvent = eventInfo->get_EvtSequence();
0090   int m_run         = eventInfo->get_RunNumber();
0091
0092   // zvertex
0093   Double_t m_zvtx = -9999;
0094   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0095
0096   if (!vertexmap) {
0097     cout << "JetValidationv2::process_event - Error can not find global vertex node " << endl;
0098     return Fun4AllReturnCodes::ABORTRUN;
0099   }
0100
0101   if(!vertexmap->empty()) {
0102     GlobalVertex* vtx = vertexmap->begin()->second;
0103     m_zvtx = vtx->get_z();
0104   }
0105
0106   // interface to reco jets
0107   JetContainer* jets_r04 = findNode::getClass<JetContainer>(topNode, m_recoJetName_r04);
0108   if (!jets_r04) {
0109     cout << "EventValidation::process_event - Error can not find DST Reco JetContainer node " << m_recoJetName_r04 << endl;
0110     return Fun4AllReturnCodes::ABORTRUN;
0111   }
0112
0113   Int_t jetPtLead    = 0;
0114   Int_t jetPtSubLead = 0;
0115
0116   Double_t jetPhiLead = 0;
0117   Double_t jetPhiSubLead = 0;
0118
0119   Double_t jetEtaLead = 0;
0120   Double_t jetEtaSubLead = 0;
0121
0122   Int_t jetIdxLead = 0;
0123
0124   for (UInt_t i = 0; i < jets_r04->size(); ++i) {
0125     Jet* jet = jets_r04->get_jet(i);
0126
0127     Double_t phi = jet->get_phi();
0128     Double_t eta = jet->get_eta();
0129     Double_t pt = jet->get_pt();
0130
0131     // exclude jets near the edge of the detector
0132     if(JetUtils::check_bad_jet_eta(eta, m_zvtx, m_R)) continue;
0133
0134     if((Int_t)pt > jetPtSubLead) {
0135       if((Int_t)pt > jetPtLead) {
0136         jetPtSubLead = jetPtLead;
0137         jetPtLead = pt;
0138
0139         jetPhiSubLead = jetPhiLead;
0140         jetPhiLead = phi;
0141
0142         jetEtaSubLead = jetEtaLead;
0143         jetEtaLead = eta;
0144
0145         jetIdxLead = i;
0146       }
0147       else {
0148         jetPtSubLead  = pt;
0149         jetPhiSubLead = phi;
0150         jetEtaSubLead = eta;
0151       }
0152     }
0153   }
0154
0155   cout << "==============" << endl;
0156   cout << "Event Jet Info" << endl;
0157   cout << "Z: " << m_zvtx << endl;
0158   cout << "Lead Jet pT: " << jetPtLead << ", phi: " << jetPhiLead << ", eta: " << jetEtaLead << ", idx: " << jetIdxLead << endl;
0159   cout << "SubLead Jet pT: " << jetPtSubLead << ", phi: " << jetPhiSubLead << ", eta: " << jetEtaSubLead << endl;
0160   cout << "==============" << endl;
0161
0162   // Get TowerInfoContainer
0163   TowerInfoContainer* towersCEMC  = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0164   TowerInfoContainer* towersIHCal = findNode::getClass<TowerInfoContainer>(topNode, m_ihcTowerNode.c_str());
0165   TowerInfoContainer* towersOHCal = findNode::getClass<TowerInfoContainer>(topNode, m_ohcTowerNode.c_str());
0166
0167   RawTowerGeomContainer* geomCEMC  = findNode::getClass<RawTowerGeomContainer>(topNode, m_emcGeomNode.c_str());
0168   RawTowerGeomContainer* geomIHCAL = findNode::getClass<RawTowerGeomContainer>(topNode, m_ihcGeomNode.c_str());
0169   RawTowerGeomContainer* geomOHCAL = findNode::getClass<RawTowerGeomContainer>(topNode, m_ohcGeomNode.c_str());
0170
0171   if (!towersCEMC || !towersIHCal || !towersOHCal) {
0172     cout << PHWHERE << "EventValidation::process_event Could not find one of "
0173          << m_emcTowerNode << ", "
0174          << m_ihcTowerNode << ", "
0175          << m_ohcTowerNode
0176          << endl;
0177     return Fun4AllReturnCodes::ABORTEVENT;
0178   }
0179
0180   if (!geomCEMC || !geomIHCAL || !geomOHCAL) {
0181     cout << PHWHERE << "EventValidation::process_event Could not find one of "
0182          << m_emcGeomNode << ", "
0183          << m_ihcGeomNode << ", "
0184          << m_ohcGeomNode
0185          << endl;
0186     return Fun4AllReturnCodes::ABORTEVENT;
0187   }
0188
0189   Jet::TYPE_comp_vec comp_vec = jets_r04->get_jet(jetIdxLead)->get_comp_vec();
0190
0191   Double_t totalCEMCEnergy = 0;
0192   Double_t totalIHCalEnergy = 0;
0193   Double_t totalOHCalEnergy = 0;
0194
0195   Double_t totalCEMCPx = 0;
0196   Double_t totalIHCalPx = 0;
0197   Double_t totalOHCalPx = 0;
0198
0199   Double_t totalCEMCPy = 0;
0200   Double_t totalIHCalPy = 0;
0201   Double_t totalOHCalPy = 0;
0202
0203   Double_t totalCEMCPz = 0;
0204   Double_t totalIHCalPz = 0;
0205   Double_t totalOHCalPz = 0;
0206
0207   UInt_t nTowers = 0;
0208
0209   stringstream name;
0210   stringstream nameSuffix;
0211   stringstream title;
0212   stringstream titleSuffix;
0213
0214   titleSuffix << "Run: " << m_run << ", Event: " << m_globalEvent
0215               << "; #phi; #eta";
0216
0217   nameSuffix << m_run << "_" << m_globalEvent;
0218
0219   name << "CEMC_h2TowerPx_" << nameSuffix.str();
0220   title << "CEMC Tower Px: " << titleSuffix.str();
0221
0222   auto h2TowerPx_ = new TH2F(name.str().c_str(), title.str().c_str(), m_bins_phi, m_phi_low, m_phi_high, m_bins_eta, m_eta_low, m_eta_high);
0223
0224   name.str("");
0225   title.str("");
0226   name << "CEMC_h2TowerPy_" << nameSuffix.str();
0227   title << "CEMC Tower Py: " << titleSuffix.str();
0228
0229   auto h2TowerPy_ = new TH2F(name.str().c_str(), title.str().c_str(), m_bins_phi, m_phi_low, m_phi_high, m_bins_eta, m_eta_low, m_eta_high);
0230
0231   name.str("");
0232   title.str("");
0233   name << "CEMC_h2TowerEnergy_" << nameSuffix.str();
0234   title << "CEMC Tower Energy: " << titleSuffix.str();
0235
0236   auto h2TowerEnergy_ = new TH2F(name.str().c_str(), title.str().c_str(), m_bins_phi, m_phi_low, m_phi_high, m_bins_eta, m_eta_low, m_eta_high);
0237
0238
0239   name.str("");
0240   name << "event-log-" << m_run << "-" << m_globalEvent << ".csv";
0241
0242   ofstream event_log(name.str());
0243   event_log << "detector,iphi,ieta,phi,eta,energy,pt,px,py,pz" << endl;
0244
0245   for (auto comp : comp_vec) {
0246     TowerInfoContainer* towers = nullptr;
0247     UInt_t towerIndex = comp.second;
0248     string det = "";
0249     Double_t* totalEnergy = nullptr;
0250     Double_t* totalPx = nullptr;
0251     Double_t* totalPy = nullptr;
0252     Double_t* totalPz = nullptr;
0253     RawTowerGeomContainer* geom = nullptr;
0254     RawTowerDefs::CalorimeterId geocaloid{RawTowerDefs::CalorimeterId::NONE};
0255
0256     if (comp.first == Jet::CEMC_TOWERINFO_RETOWER) {
0257       towers = towersCEMC;
0258       det = "CEMC";
0259       totalEnergy = &totalCEMCEnergy;
0260       totalPx = &totalCEMCPx;
0261       totalPy = &totalCEMCPy;
0262       totalPz = &totalCEMCPz;
0263       geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0264       geom = geomIHCAL;
0265     }
0266     if (comp.first == Jet::HCALIN_TOWERINFO) {
0267       towers = towersIHCal;
0268       det = "IHCal";
0269       totalEnergy = &totalIHCalEnergy;
0270       totalPx = &totalIHCalPx;
0271       totalPy = &totalIHCalPy;
0272       totalPz = &totalIHCalPz;
0273       geocaloid = RawTowerDefs::CalorimeterId::HCALIN;
0274       geom = geomIHCAL;
0275     }
0276     if (comp.first == Jet::HCALOUT_TOWERINFO) {
0277       towers = towersIHCal;
0278       det = "OHCal";
0279       totalEnergy = &totalOHCalEnergy;
0280       totalPx = &totalOHCalPx;
0281       totalPy = &totalOHCalPy;
0282       totalPz = &totalOHCalPz;
0283       geocaloid = RawTowerDefs::CalorimeterId::HCALOUT;
0284       geom = geomOHCAL;
0285     }
0286     if(towers == nullptr) {
0287       cout << "Unknown Detector: " << comp.first << ", Index: " << towerIndex << endl;
0288       continue;
0289     }
0290
0291     UInt_t calokey = towers->encode_key(towerIndex);
0292     Int_t ieta     = towers->getTowerEtaBin(calokey);
0293     Int_t iphi     = towers->getTowerPhiBin(calokey);
0294
0295     // get the radius of the detector
0296     const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(geocaloid, ieta, iphi);
0297     RawTowerGeom* tower_geom = geom->get_tower_geometry(key);
0298     Double_t r = tower_geom->get_center_radius();
0299
0300     if(det == "CEMC") {
0301         const RawTowerDefs::keytype EMCal_key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, 0, 0);
0302         RawTowerGeom *EMCal_tower_geom = geomCEMC->get_tower_geometry(EMCal_key);
0303         r = EMCal_tower_geom->get_center_radius();
0304     }
0305
0306     TowerInfo* tower = towers->get_tower_at_channel(towerIndex);
0307
0308     Double_t energy = tower->get_energy();
0309     *totalEnergy = *totalEnergy + energy;
0310
0311     Double_t phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0312     Double_t towereta = tower_geom->get_eta();
0313     Double_t z0 = sinh(towereta) * r;
0314     Double_t z = z0 - m_zvtx;
0315     Double_t eta = asinh(z / r);  // eta after shift from vertex
0316     Double_t pt = energy / cosh(eta);
0317     Double_t px = pt * cos(phi);
0318     Double_t py = pt * sin(phi);
0319     Double_t pz = pt * sinh(eta);
0320
0321     *totalPx = *totalPx + px;
0322     *totalPy = *totalPy + py;
0323     *totalPz = *totalPz + pz;
0324
0325     if(det == "CEMC") {
0326       h2TowerEnergy_->Fill(phi, eta, energy);
0327       h2TowerPx_->Fill(phi, eta, px);
0328       h2TowerPy_->Fill(phi, eta, py);
0329     }
0330
0331     // ensure tower is flagged as good
0332     if(!tower->get_isGood()) {
0333       cout << "Bad Tower in Jet" << endl;
0334     }
0335
0336     if(det == "CEMC") {
0337       cout << "Detector: " << det << ", iphi: " << iphi << ", ieta: " << ieta << ", energy: " << energy
0338            << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
0339       ++nTowers;
0340     }
0341     event_log << det << "," << iphi << "," << ieta << "," << phi << "," << eta << "," << energy << "," << pt << "," << px << "," << py << "," << pz << endl;
0342   }
0343
0344   event_log.close();
0345
0346   h2TowerEnergy.push_back(h2TowerEnergy_);
0347   h2TowerPx.push_back(h2TowerPx_);
0348   h2TowerPy.push_back(h2TowerPy_);
0349
0350   cout << "============" << endl;
0351   cout << "Total Energy" << endl;
0352   cout << "CEMC: " << totalCEMCEnergy << endl;
0353   cout << "IHCal: " << totalIHCalEnergy << endl;
0354   cout << "OHCal: " << totalOHCalEnergy << endl;
0355   cout << "============" << endl;
0356
0357   Double_t totalPx = totalCEMCPx + totalIHCalPx + totalOHCalPx;
0358   Double_t totalPy = totalCEMCPy + totalIHCalPy + totalOHCalPy;
0359   Double_t totalPz = totalCEMCPz + totalIHCalPz + totalOHCalPz;
0360
0361   Double_t totalPt = sqrt(totalPx*totalPx + totalPy*totalPy);
0362   Double_t phi     = atan2(totalPy, totalPx);
0363   Double_t eta     = asinh(totalPz / totalPt);
0364
0365   cout << "============" << endl;
0366   cout << "Lead Jet: pt: " << (Int_t)totalPt << ", px: " << (Int_t)totalPx << ", py: " << (Int_t)totalPy << ", pz: " << (Int_t)totalPz
0367                            << ", phi: " << phi << ", eta: " << eta << ", Towers: " << nTowers << endl;
0368   cout << "============" << endl;
0369
0370
0371   return Fun4AllReturnCodes::EVENT_OK;
0372 }
0373
0374 //____________________________________________________________________________..
0375 Int_t EventValidation::ResetEvent(PHCompositeNode *topNode)
0376 {
0377   return Fun4AllReturnCodes::EVENT_OK;
0378 }
0379
0380 //____________________________________________________________________________..
0381 Int_t EventValidation::End(PHCompositeNode *topNode)
0382 {
0383   cout << "EventValidation::End(PHCompositeNode *topNode) This is the End..." << endl;
0384   TFile output(m_outputFile.c_str(),"recreate");
0385
0386   output.cd();
0387
0388   output.mkdir("energy");
0389   output.mkdir("px");
0390   output.mkdir("py");
0391
0392   for(UInt_t i = 0; i < h2TowerEnergy.size(); ++i) {
0393     output.cd("energy");
0394     h2TowerEnergy[i]->Write();
0395
0396     output.cd("px");
0397     h2TowerPx[i]->Write();
0398
0399     output.cd("py");
0400     h2TowerPy[i]->Write();
0401   }
0402
0403   output.Close();
0404
0405   cout << "EventValidation::End(PHCompositeNode *topNode) This is the End..." << endl;
0406   return Fun4AllReturnCodes::EVENT_OK;
0407 }
0408