Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:09 

0001 // module for producing a TTree with jet information for doing jet validation studies
0002 //  for questions/bugs please contact Virginia Bailey vbailey13@gsu.edu
0003 #include <JetValidation.h>
0004 // -- c++
0005 #include <cmath>
0006 #include <iomanip>
0007 #include <iostream>
0008 #include <sstream>
0009 #include <vector>
0010 // -- event / trigger
0011 #include <ffaobjects/EventHeader.h>
0012 #include <ffarawobjects/Gl1Packet.h>
0013 // -- fun4all
0014 #include <fun4all/Fun4AllBase.h>
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 // -- vertex
0017 #include <globalvertex/GlobalVertex.h>
0018 #include <globalvertex/GlobalVertexMap.h>
0019 // -- jet
0020 #include <jetbase/JetContainer.h>
0021 // -- DST node
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/getClass.h>
0024 // -- root
0025 #include <TFile.h>
0026 #include <TTree.h>
0027 // -- Tower stuff
0028 #include <calobase/TowerInfo.h>
0029 #include <calobase/RawTowerDefs.h>
0030 #include <calobase/TowerInfoDefs.h>
0031 #include <calobase/TowerInfoContainer.h>
0032 
0033 using std::cout;
0034 using std::endl;
0035 using std::string;
0036 using std::to_string;
0037 using std::stringstream;
0038 
0039 //____________________________________________________________________________..
0040 JetValidation::JetValidation()
0041   : SubsysReco("JetValidation")
0042   , m_recoJetName_r02("AntiKt_Tower_r02_Sub1")
0043   , m_recoJetName_r04("AntiKt_Tower_r04_Sub1")
0044   , m_recoJetName_r06("AntiKt_Tower_r06_Sub1")
0045   , m_outputTreeFile(nullptr)
0046   , m_outputQAFile(nullptr)
0047   , m_outputTreeFileName("test.root")
0048   , m_outputQAFileName("qa.root")
0049   , m_emcTowerNodeBase("TOWERINFO_CALIB_CEMC")
0050   , m_emcTowerNode("TOWERINFO_CALIB_CEMC_RETOWER")
0051   , m_ihcalTowerNode("TOWERINFO_CALIB_HCALIN")
0052   , m_ohcalTowerNode("TOWERINFO_CALIB_HCALOUT")
0053   , m_emcTowerNodeSub("TOWERINFO_CALIB_CEMC_RETOWER_SUB1")
0054   , m_ihcalTowerNodeSub("TOWERINFO_CALIB_HCALIN_SUB1")
0055   , m_ohcalTowerNodeSub("TOWERINFO_CALIB_HCALOUT_SUB1")
0056   , m_saveTree(false)
0057   , m_lowPtThreshold(10) /*GeV*/
0058   , m_highPtThreshold(40) /*GeV*/
0059   , m_subLeadPtThreshold(5) /*GeV*/
0060   , m_highPtJetCtr(0)
0061   , m_bins_phi(64)
0062   , m_bins_eta(24)
0063   , m_eta_low(-1)
0064   , m_eta_high(1)
0065   , m_bins_pt(200)
0066   , m_pt_low(0) /*GeV*/
0067   , m_pt_high(200) /*GeV*/
0068   , m_bins_events(18)
0069   , m_bins_zvtx(400)
0070   , m_zvtx_low(-100)
0071   , m_zvtx_high(100)
0072   , m_bkg_tower_energy(0.6) /*GeV*/
0073   , m_bkg_tower_neighbor_energy(0.06) /*GeV*/
0074   , m_bkg_towers(6)
0075   , m_T(nullptr)
0076   , m_run(0)
0077   , m_globalEvent(0)
0078   , m_event(0)
0079   , m_nJets_r02(0)
0080   , m_nJets_r04(0)
0081   , m_nJets_r06(0)
0082   , m_nJet_r02(0)
0083   , m_id_r02()
0084   , m_nComponent_r02()
0085   , m_eta_r02()
0086   , m_phi_r02()
0087   , m_e_r02()
0088   , m_pt_r02()
0089   , m_nJet_r04(0)
0090   , m_id_r04()
0091   , m_nComponent_r04()
0092   , m_eta_r04()
0093   , m_phi_r04()
0094   , m_e_r04()
0095   , m_pt_r04()
0096   , m_nJet_r06(0)
0097   , m_id_r06()
0098   , m_nComponent_r06()
0099   , m_eta_r06()
0100   , m_phi_r06()
0101   , m_e_r06()
0102   , m_pt_r06()
0103 {
0104   cout << "JetValidation::JetValidation(const std::string &name) Calling ctor" << endl;
0105 }
0106 
0107 //____________________________________________________________________________..
0108 JetValidation::~JetValidation()
0109 {
0110   cout << "JetValidation::~JetValidation() Calling dtor" << endl;
0111 }
0112 
0113 //____________________________________________________________________________..
0114 Int_t JetValidation::Init(PHCompositeNode *topNode)
0115 {
0116   cout << "JetValidation::Init(PHCompositeNode *topNode) Initializing" << endl;
0117 
0118   if(m_saveTree) {
0119     m_outputTreeFile = new TFile(m_outputTreeFileName.c_str(),"RECREATE");
0120     cout << "JetValidation::Init - Output to " << m_outputTreeFileName << endl;
0121 
0122     // configure Tree
0123     m_T = new TTree("T", "T");
0124     m_T->Branch("event", &m_globalEvent, "event/I");
0125     m_T->Branch("run", &m_run, "run/I");
0126     m_T->Branch("zvtx", &m_zvtx);
0127     m_T->Branch("scaledVector", &m_scaledVector);
0128 
0129     m_T->Branch("towersCEMCBase_isGood", &m_towersCEMCBase_isGood);
0130     m_T->Branch("towersCEMCBase_energy", &m_towersCEMCBase_energy);
0131     m_T->Branch("towersCEMCBase_time",   &m_towersCEMCBase_time);
0132 
0133     m_T->Branch("towersCEMC_isGood",  &m_towersCEMC_isGood);
0134     m_T->Branch("towersIHCal_isGood", &m_towersIHCal_isGood);
0135     m_T->Branch("towersOHCal_isGood", &m_towersOHCal_isGood);
0136 
0137     m_T->Branch("towersCEMC_energy",  &m_towersCEMC_energy);
0138     m_T->Branch("towersIHCal_energy", &m_towersIHCal_energy);
0139     m_T->Branch("towersOHCal_energy", &m_towersOHCal_energy);
0140 
0141     m_T->Branch("towersCEMCSub_energy",  &m_towersCEMCSub_energy);
0142     m_T->Branch("towersIHCalSub_energy", &m_towersIHCalSub_energy);
0143     m_T->Branch("towersOHCalSub_energy", &m_towersOHCalSub_energy);
0144 
0145     m_T->Branch("nJet_r02", &m_nJet_r02, "nJet_r02/I");
0146     m_T->Branch("nComponent_r02", &m_nComponent_r02);
0147     m_T->Branch("id_r02", &m_id_r02);
0148     m_T->Branch("eta_r02", &m_eta_r02);
0149     m_T->Branch("phi_r02", &m_phi_r02);
0150     m_T->Branch("e_r02", &m_e_r02);
0151     m_T->Branch("pt_r02", &m_pt_r02);
0152 
0153     m_T->Branch("nJet_r04", &m_nJet_r04, "nJet_r04/I");
0154     m_T->Branch("nComponent_r04", &m_nComponent_r04);
0155     m_T->Branch("id_r04", &m_id_r04);
0156     m_T->Branch("eta_r04", &m_eta_r04);
0157     m_T->Branch("phi_r04", &m_phi_r04);
0158     m_T->Branch("e_r04", &m_e_r04);
0159     m_T->Branch("pt_r04", &m_pt_r04);
0160 
0161     m_T->Branch("nJet_r06", &m_nJet_r06, "nJet_r06/I");
0162     m_T->Branch("nComponent_r04", &m_nComponent_r04);
0163     m_T->Branch("id_r06", &m_id_r06);
0164     m_T->Branch("eta_r06", &m_eta_r06);
0165     m_T->Branch("phi_r06", &m_phi_r06);
0166     m_T->Branch("e_r06", &m_e_r06);
0167     m_T->Branch("pt_r06", &m_pt_r06);
0168   }
0169 
0170   m_outputQAFile = new TFile(m_outputQAFileName.c_str(),"RECREATE");
0171   cout << "JetValidation::Init - Output to " << m_outputQAFileName << endl;
0172 
0173   for(auto suffix : JetEvent_Status_vec) {
0174     // pt
0175     hJetPt_r02.push_back(new TH1F(("hJetPt_r02_"+suffix).c_str(), "Jet: R = 0.2; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0176     hJetPt_r04.push_back(new TH1F(("hJetPt_r04_"+suffix).c_str(), "Jet: R = 0.4; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0177     hJetPt_r06.push_back(new TH1F(("hJetPt_r06_"+suffix).c_str(), "Jet: R = 0.6; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0178 
0179     hJetPt_r02_bkg.push_back(new TH1F(("hJetPt_r02_"+suffix+"_bkg").c_str(), "Jet: R = 0.2; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0180     hJetPt_r04_bkg.push_back(new TH1F(("hJetPt_r04_"+suffix+"_bkg").c_str(), "Jet: R = 0.4; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0181     hJetPt_r06_bkg.push_back(new TH1F(("hJetPt_r06_"+suffix+"_bkg").c_str(), "Jet: R = 0.6; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0182 
0183     hJetPt_r02_nobkg.push_back(new TH1F(("hJetPt_r02_"+suffix+"_nobkg").c_str(), "Jet: R = 0.2; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0184     hJetPt_r04_nobkg.push_back(new TH1F(("hJetPt_r04_"+suffix+"_nobkg").c_str(), "Jet: R = 0.4; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0185     hJetPt_r06_nobkg.push_back(new TH1F(("hJetPt_r06_"+suffix+"_nobkg").c_str(), "Jet: R = 0.6; p_{T} [GeV]; Counts", m_bins_pt, m_pt_low, m_pt_high));
0186 
0187     // deltaPhi
0188     hJetDeltaPhi_r02.push_back(new TH1F(("hJetDeltaPhi_r02_"+suffix).c_str(), "Jet: R = 0.2; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0189     hJetDeltaPhi_r04.push_back(new TH1F(("hJetDeltaPhi_r04_"+suffix).c_str(), "Jet: R = 0.4; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0190     hJetDeltaPhi_r06.push_back(new TH1F(("hJetDeltaPhi_r06_"+suffix).c_str(), "Jet: R = 0.6; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0191 
0192     hJetDeltaPhi_r02_bkg.push_back(new TH1F(("hJetDeltaPhi_r02_"+suffix+"_bkg").c_str(), "Jet: R = 0.2; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0193     hJetDeltaPhi_r04_bkg.push_back(new TH1F(("hJetDeltaPhi_r04_"+suffix+"_bkg").c_str(), "Jet: R = 0.4; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0194     hJetDeltaPhi_r06_bkg.push_back(new TH1F(("hJetDeltaPhi_r06_"+suffix+"_bkg").c_str(), "Jet: R = 0.6; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0195 
0196     hJetDeltaPhi_r02_nobkg.push_back(new TH1F(("hJetDeltaPhi_r02_"+suffix+"_nobkg").c_str(), "Jet: R = 0.2; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0197     hJetDeltaPhi_r04_nobkg.push_back(new TH1F(("hJetDeltaPhi_r04_"+suffix+"_nobkg").c_str(), "Jet: R = 0.4; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0198     hJetDeltaPhi_r06_nobkg.push_back(new TH1F(("hJetDeltaPhi_r06_"+suffix+"_nobkg").c_str(), "Jet: R = 0.6; |#Delta #phi|; Counts", m_bins_phi, 0, M_PI));
0199 
0200     // eta_phi
0201     h2JetEtaPhi_r02.push_back(new TH2F(("h2JetEtaPhi_r02_"+suffix).c_str(), "Jet: R = 0.2; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0202     h2JetEtaPhi_r04.push_back(new TH2F(("h2JetEtaPhi_r04_"+suffix).c_str(), "Jet: R = 0.4; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0203     h2JetEtaPhi_r06.push_back(new TH2F(("h2JetEtaPhi_r06_"+suffix).c_str(), "Jet: R = 0.6; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0204 
0205     h2JetEtaPhi_r02_bkg.push_back(new TH2F(("h2JetEtaPhi_r02_"+suffix+"_bkg").c_str(), "Jet: R = 0.2; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0206     h2JetEtaPhi_r04_bkg.push_back(new TH2F(("h2JetEtaPhi_r04_"+suffix+"_bkg").c_str(), "Jet: R = 0.4; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0207     h2JetEtaPhi_r06_bkg.push_back(new TH2F(("h2JetEtaPhi_r06_"+suffix+"_bkg").c_str(), "Jet: R = 0.6; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0208 
0209     h2JetEtaPhi_r02_nobkg.push_back(new TH2F(("h2JetEtaPhi_r02_"+suffix+"_nobkg").c_str(), "Jet: R = 0.2; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0210     h2JetEtaPhi_r04_nobkg.push_back(new TH2F(("h2JetEtaPhi_r04_"+suffix+"_nobkg").c_str(), "Jet: R = 0.4; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0211     h2JetEtaPhi_r06_nobkg.push_back(new TH2F(("h2JetEtaPhi_r06_"+suffix+"_nobkg").c_str(), "Jet: R = 0.6; #phi; #eta", m_bins_phi, -M_PI, M_PI, m_bins_eta, m_eta_low, m_eta_high));
0212   }
0213 
0214   hEvents    = new TH1F("hEvents","Events; Status; Counts", m_bins_events, 0, m_bins_events);
0215   hEventsBkg = new TH1F("hEventsBkg","Background Events; Status; Counts", m_bins_events, 0, m_bins_events);
0216   hZVtx      = new TH1F("hZVtx","Z Vertex; z [cm]; Counts", m_bins_zvtx, m_zvtx_low, m_zvtx_high);
0217 
0218   return Fun4AllReturnCodes::EVENT_OK;
0219 }
0220 
0221 //____________________________________________________________________________..
0222 Int_t JetValidation::process_event(PHCompositeNode *topNode)
0223 {
0224   EventHeader* eventInfo = findNode::getClass<EventHeader>(topNode,"EventHeader");
0225   if(!eventInfo)
0226   {
0227     cout << PHWHERE << "JetValidation::process_event - Fatal Error - EventHeader node is missing. " << endl;
0228     return Fun4AllReturnCodes::ABORTEVENT;
0229   }
0230 
0231   m_globalEvent = eventInfo->get_EvtSequence();
0232   m_run         = eventInfo->get_RunNumber();
0233 
0234   if (m_event % 20 == 0) cout << "Progress: " << m_event << ", Global: " << m_globalEvent << endl;
0235   ++m_event;
0236 
0237   // zvertex
0238   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0239   if (!vertexmap) {
0240     cout << "JetValidation::process_event - Error can not find global vertex node " << endl;
0241     return Fun4AllReturnCodes::ABORTRUN;
0242   }
0243   if(!vertexmap->empty()) {
0244     GlobalVertex* vtx = vertexmap->begin()->second;
0245     m_zvtx = vtx->get_z();
0246   }
0247   else m_zvtx = -9999;
0248 
0249   hZVtx->Fill(m_zvtx);
0250 
0251   // grab the gl1 data
0252   Gl1Packet* gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0253   if (!gl1PacketInfo)
0254   {
0255     cout << PHWHERE << "JetValidation::process_event: GL1Packet node is missing" << endl;
0256     return Fun4AllReturnCodes::ABORTRUN;
0257   }
0258 
0259   uint64_t scaledvec = gl1PacketInfo->getScaledVector();
0260   // uint64_t scaledvec = gl1PacketInfo->getTriggerVector();
0261   for (Int_t i = 0; i < 64; ++i)
0262   {
0263     Bool_t trig_decision = ((scaledvec & 0x1U) == 0x1U);
0264     m_scaledVector.push_back(trig_decision);
0265     scaledvec = (scaledvec >> 1U) & 0xffffffffU;
0266   }
0267 
0268   // Get TowerInfoContainer
0269   // Base EMCal Towers
0270   TowerInfoContainer* towersCEMCBase = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNodeBase.c_str());
0271 
0272   // unsubtracted
0273   TowerInfoContainer* towersCEMC  = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0274   TowerInfoContainer* towersIHCal = findNode::getClass<TowerInfoContainer>(topNode, m_ihcalTowerNode.c_str());
0275   TowerInfoContainer* towersOHCal = findNode::getClass<TowerInfoContainer>(topNode, m_ohcalTowerNode.c_str());
0276 
0277   // subtracted
0278   TowerInfoContainer* towersCEMCSub  = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNodeSub.c_str());
0279   TowerInfoContainer* towersIHCalSub = findNode::getClass<TowerInfoContainer>(topNode, m_ihcalTowerNodeSub.c_str());
0280   TowerInfoContainer* towersOHCalSub = findNode::getClass<TowerInfoContainer>(topNode, m_ohcalTowerNodeSub.c_str());
0281 
0282   if (!towersCEMCBase || !towersCEMC    || !towersIHCal    || !towersOHCal    ||
0283       !towersCEMCSub  || !towersIHCalSub || !towersOHCalSub) {
0284     cout << PHWHERE << "JetValidation::process_event Could not find one of "
0285          << m_emcTowerNodeBase  << ", "
0286          << m_emcTowerNode      << ", "
0287          << m_ihcalTowerNode    << ", "
0288          << m_ohcalTowerNode    << ", "
0289          << m_emcTowerNodeSub   << ", "
0290          << m_ihcalTowerNodeSub << ", "
0291          << m_ohcalTowerNodeSub
0292          << endl;
0293     return Fun4AllReturnCodes::ABORTEVENT;
0294   }
0295 
0296   std::vector<Float_t> towerEnergy(towersCEMC->size());
0297   std::vector<Float_t> towerEnergySub(towersCEMC->size());
0298 
0299   // loop over towers
0300   for(UInt_t towerIndex = 0; towerIndex < towersCEMC->size(); ++towerIndex) {
0301     // unsubtracted
0302     TowerInfo* tower = towersCEMC->get_tower_at_channel(towerIndex);
0303 
0304     Float_t energy = (tower->get_isGood()) ? tower->get_energy() : 0;
0305 
0306     m_towersCEMC_energy.push_back(tower->get_energy());
0307     m_towersCEMC_isGood.push_back(tower->get_isGood());
0308 
0309     tower = towersIHCal->get_tower_at_channel(towerIndex);
0310     m_towersIHCal_energy.push_back(tower->get_energy());
0311     m_towersIHCal_isGood.push_back(tower->get_isGood());
0312 
0313     energy += (tower->get_isGood()) ? tower->get_energy() : 0;
0314 
0315     tower = towersOHCal->get_tower_at_channel(towerIndex);
0316     m_towersOHCal_energy.push_back(tower->get_energy());
0317     m_towersOHCal_isGood.push_back(tower->get_isGood());
0318 
0319     energy += (tower->get_isGood()) ? tower->get_energy() : 0;
0320 
0321     towerEnergy[towerIndex] = energy;
0322 
0323     // subtracted
0324     tower = towersCEMCSub->get_tower_at_channel(towerIndex);
0325     m_towersCEMCSub_energy.push_back(tower->get_energy());
0326 
0327     energy = (tower->get_isGood()) ? tower->get_energy() : 0;
0328 
0329     tower = towersIHCalSub->get_tower_at_channel(towerIndex);
0330     m_towersIHCalSub_energy.push_back(tower->get_energy());
0331 
0332     energy += (tower->get_isGood()) ? tower->get_energy() : 0;
0333 
0334     tower = towersOHCalSub->get_tower_at_channel(towerIndex);
0335     m_towersOHCalSub_energy.push_back(tower->get_energy());
0336 
0337     energy += (tower->get_isGood()) ? tower->get_energy() : 0;
0338 
0339     towerEnergySub[towerIndex] = energy;
0340   }
0341 
0342   Bool_t m_hasBkg = isBackgroundEvent(towerEnergy);
0343 
0344   hEvents->Fill((Int_t)EventStatus::ALL);
0345   if(m_hasBkg) {
0346     hEventsBkg->Fill((Int_t)EventStatus::ALL);
0347   }
0348 
0349   if(m_zvtx != -9999) {
0350     hEvents->Fill((Int_t)EventStatus::ZVTX);
0351     if(m_hasBkg) {
0352       hEventsBkg->Fill((Int_t)EventStatus::ZVTX);
0353     }
0354   }
0355   if(fabs(m_zvtx) < m_zvtx_max[0]) {
0356     hEvents->Fill((Int_t)EventStatus::ZVTX60);
0357     if(m_hasBkg) {
0358       hEventsBkg->Fill((Int_t)EventStatus::ZVTX60);
0359     }
0360   }
0361   if(fabs(m_zvtx) < m_zvtx_max[1]) {
0362     hEvents->Fill((Int_t)EventStatus::ZVTX50);
0363     if(m_hasBkg) {
0364       hEventsBkg->Fill((Int_t)EventStatus::ZVTX50);
0365     }
0366   }
0367   if(fabs(m_zvtx) < m_zvtx_max[2]) {
0368     hEvents->Fill((Int_t)EventStatus::ZVTX30);
0369     if(m_hasBkg) {
0370       hEventsBkg->Fill((Int_t)EventStatus::ZVTX30);
0371     }
0372   }
0373   if(fabs(m_zvtx) < m_zvtx_max[3]) {
0374     hEvents->Fill((Int_t)EventStatus::ZVTX20);
0375     if(m_hasBkg) {
0376       hEventsBkg->Fill((Int_t)EventStatus::ZVTX20);
0377     }
0378   }
0379   if(fabs(m_zvtx) < m_zvtx_max[4]) {
0380     hEvents->Fill((Int_t)EventStatus::ZVTX10);
0381     if(m_hasBkg) {
0382       hEventsBkg->Fill((Int_t)EventStatus::ZVTX10);
0383     }
0384   }
0385 
0386   std::vector<Bool_t> eventTrig_status(JetEvent_Status_vec.size());
0387 
0388   // Fill Event Status
0389   if(m_scaledVector[(Int_t)Trigger::MBD_NS_1]) {
0390     // ALL
0391     hEvents->Fill((Int_t)EventStatus::ALL_MBDNS1);
0392     eventTrig_status[(Int_t)JetEvent_Status::ALL_MBDNS1] = true;
0393     if(m_hasBkg) {
0394       hEventsBkg->Fill((Int_t)EventStatus::ALL_MBDNS1);
0395     }
0396 
0397     // ZVTX60
0398     if(fabs(m_zvtx) < m_zvtx_max[0]) {
0399       hEvents->Fill((Int_t)EventStatus::ZVTX60_MBDNS1);
0400       eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_MBDNS1] = true;
0401       if(m_hasBkg) {
0402         hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_MBDNS1);
0403       }
0404     }
0405 
0406     // Jet 8 GeV
0407     if(m_scaledVector[(Int_t)Trigger::JET_8]) {
0408       hEvents->Fill((Int_t)EventStatus::ALL_MBDNS1_JET8);
0409       eventTrig_status[(Int_t)JetEvent_Status::ALL_MBDNS1_JET8] = true;
0410 
0411       if(m_hasBkg) {
0412         hEventsBkg->Fill((Int_t)EventStatus::ALL_MBDNS1_JET8);
0413       }
0414 
0415       if(fabs(m_zvtx) < m_zvtx_max[0]) {
0416         hEvents->Fill((Int_t)EventStatus::ZVTX60_MBDNS1_JET8);
0417         eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_MBDNS1_JET8] = true;
0418 
0419         if(m_hasBkg) {
0420           hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_MBDNS1_JET8);
0421         }
0422       }
0423     }
0424 
0425     // Jet 10 GeV
0426     if(m_scaledVector[(Int_t)Trigger::JET_10]) {
0427       hEvents->Fill((Int_t)EventStatus::ALL_MBDNS1_JET10);
0428       eventTrig_status[(Int_t)JetEvent_Status::ALL_MBDNS1_JET10] = true;
0429 
0430       if(m_hasBkg) {
0431         hEventsBkg->Fill((Int_t)EventStatus::ALL_MBDNS1_JET10);
0432       }
0433 
0434       if(fabs(m_zvtx) < m_zvtx_max[0]) {
0435         hEvents->Fill((Int_t)EventStatus::ZVTX60_MBDNS1_JET10);
0436         eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_MBDNS1_JET10] = true;
0437 
0438         if(m_hasBkg) {
0439           hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_MBDNS1_JET10);
0440         }
0441       }
0442     }
0443 
0444     // Jet 12 GeV
0445     if(m_scaledVector[(Int_t)Trigger::JET_12]) {
0446       hEvents->Fill((Int_t)EventStatus::ALL_MBDNS1_JET12);
0447       eventTrig_status[(Int_t)JetEvent_Status::ALL_MBDNS1_JET12] = true;
0448 
0449       if(m_hasBkg) {
0450         hEventsBkg->Fill((Int_t)EventStatus::ALL_MBDNS1_JET12);
0451       }
0452 
0453       if(fabs(m_zvtx) < m_zvtx_max[0]) {
0454         hEvents->Fill((Int_t)EventStatus::ZVTX60_MBDNS1_JET12);
0455         eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_MBDNS1_JET12] = true;
0456 
0457         if(m_hasBkg) {
0458           hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_MBDNS1_JET12);
0459         }
0460       }
0461     }
0462   }
0463 
0464   if(fabs(m_zvtx) < m_zvtx_max[0]) {
0465     eventTrig_status[(Int_t)JetEvent_Status::ZVTX60] = true;
0466 
0467     // Jet 8 GeV
0468     if(m_scaledVector[(Int_t)Trigger::JET_8]) {
0469       eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_JET8] = true;
0470       hEvents->Fill((Int_t)EventStatus::ZVTX60_JET8);
0471 
0472       if(m_hasBkg) {
0473         hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_JET8);
0474       }
0475     }
0476 
0477     // Jet 10 GeV
0478     if(m_scaledVector[(Int_t)Trigger::JET_10]) {
0479       eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_JET10] = true;
0480       hEvents->Fill((Int_t)EventStatus::ZVTX60_JET10);
0481 
0482       if(m_hasBkg) {
0483         hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_JET10);
0484       }
0485     }
0486 
0487     // Jet 12 GeV
0488     if(m_scaledVector[(Int_t)Trigger::JET_12]) {
0489       eventTrig_status[(Int_t)JetEvent_Status::ZVTX60_JET12] = true;
0490       hEvents->Fill((Int_t)EventStatus::ZVTX60_JET12);
0491 
0492       if(m_hasBkg) {
0493         hEventsBkg->Fill((Int_t)EventStatus::ZVTX60_JET12);
0494       }
0495     }
0496   }
0497 
0498   // interface to reco jets
0499   JetContainer* jets_r02 = findNode::getClass<JetContainer>(topNode, m_recoJetName_r02);
0500   JetContainer* jets_r04 = findNode::getClass<JetContainer>(topNode, m_recoJetName_r04);
0501   JetContainer* jets_r06 = findNode::getClass<JetContainer>(topNode, m_recoJetName_r06);
0502   if (!jets_r02 || !jets_r04 || !jets_r06) {
0503     cout << "JetValidation::process_event - Error can not find one of DST Reco JetContainer node "
0504          << m_recoJetName_r02 << ", "
0505          << m_recoJetName_r04 << ", "
0506          << m_recoJetName_r06 << ", " << endl;
0507     return Fun4AllReturnCodes::ABORTRUN;
0508   }
0509 
0510   // get reco jets
0511   m_nJet_r02 = 0;
0512   m_nJet_r04 = 0;
0513   m_nJet_r06 = 0;
0514 
0515   // R = 0.2
0516   for (auto jet : *jets_r02) {
0517     Float_t jet_pt = jet->get_pt();
0518     if (jet_pt < m_lowPtThreshold || fabs(jet->get_eta()) > m_eta_high-0.2) continue;  // to remove noise jets
0519 
0520     m_id_r02.push_back(jet->get_id());
0521     m_nComponent_r02.push_back(jet->size_comp());
0522     m_eta_r02.push_back(jet->get_eta());
0523     m_phi_r02.push_back(jet->get_phi());
0524     m_e_r02.push_back(jet->get_e());
0525     m_pt_r02.push_back(jet_pt);
0526 
0527     for (UInt_t i = 0; i < eventTrig_status.size(); ++i) {
0528       if (eventTrig_status[i]) {
0529         hJetPt_r02[i]->Fill(jet_pt);
0530         if(m_hasBkg) hJetPt_r02_bkg[i]->Fill(jet_pt);
0531         else         hJetPt_r02_nobkg[i]->Fill(jet_pt);
0532       }
0533     }
0534 
0535     ++m_nJet_r02;
0536   }
0537 
0538   Bool_t hasHighPtJet = false;
0539   Float_t leadPt = 0;
0540   Float_t subLeadPt = 0;
0541   Float_t leadPhi = 0;
0542   Float_t subLeadPhi = 0;
0543   // R = 0.4
0544   for (auto jet : *jets_r04) {
0545 
0546     if (fabs(jet->get_eta()) > m_eta_high-0.4) continue;  // to remove noise jets
0547 
0548     Float_t jet_pt = jet->get_pt();
0549 
0550     if(jet_pt > subLeadPt) {
0551         if(jet_pt > leadPt) {
0552             subLeadPt = leadPt;
0553             subLeadPhi = leadPhi;
0554             leadPt = jet_pt;
0555             leadPhi = jet->get_phi();
0556         }
0557         else {
0558           subLeadPt = jet_pt;
0559           subLeadPhi = jet->get_phi();
0560         }
0561     }
0562 
0563     if (jet_pt < m_lowPtThreshold) continue;  // to remove noise jets
0564 
0565     m_id_r04.push_back(jet->get_id());
0566     m_nComponent_r04.push_back(jet->size_comp());
0567     m_eta_r04.push_back(jet->get_eta());
0568     m_phi_r04.push_back(jet->get_phi());
0569     m_e_r04.push_back(jet->get_e());
0570     m_pt_r04.push_back(jet_pt);
0571 
0572     for (UInt_t i = 0; i < eventTrig_status.size(); ++i) {
0573       if (eventTrig_status[i]) {
0574         hJetPt_r04[i]->Fill(jet_pt);
0575         h2JetEtaPhi_r04[i]->Fill(jet->get_phi(), jet->get_eta());
0576         if(m_hasBkg) {
0577           hJetPt_r04_bkg[i]->Fill(jet_pt);
0578           h2JetEtaPhi_r04_bkg[i]->Fill(jet->get_phi(), jet->get_eta());
0579         }
0580         else{
0581           hJetPt_r04_nobkg[i]->Fill(jet_pt);
0582           h2JetEtaPhi_r04_nobkg[i]->Fill(jet->get_phi(), jet->get_eta());
0583         }
0584       }
0585     }
0586 
0587     if(jet_pt >= m_highPtThreshold) {
0588       hasHighPtJet = true;
0589       ++m_highPtJetCtr;
0590     }
0591 
0592     ++m_nJet_r04;
0593   }
0594 
0595   if(subLeadPt >= m_subLeadPtThreshold && leadPt >= m_lowPtThreshold) {
0596     // ensure that |leadPhi - subLeadPhi|  <= pi
0597     Float_t deltaPhi = (fabs(leadPhi-subLeadPhi) > M_PI) ? 2*M_PI-fabs(leadPhi-subLeadPhi) : fabs(leadPhi-subLeadPhi);
0598     for (UInt_t i = 0; i < eventTrig_status.size(); ++i) {
0599       if (eventTrig_status[i]) {
0600         hJetDeltaPhi_r04[i]->Fill(deltaPhi);
0601         if(m_hasBkg) hJetDeltaPhi_r04_bkg[i]->Fill(deltaPhi);
0602         else         hJetDeltaPhi_r04_nobkg[i]->Fill(deltaPhi);
0603       }
0604     }
0605   }
0606 
0607   if(hasHighPtJet && !m_hasBkg) {
0608     cout << "High Pt Jet Found, Event: " << m_globalEvent << ", Jet Pt: " << (Int_t)leadPt << " GeV" << endl;
0609 
0610     stringstream name;
0611     stringstream title;
0612 
0613     name << "h2TowerEnergy_" << m_run << "_" << m_globalEvent << "_" << (Int_t)leadPt;
0614     title << "Tower Energy, Run: " << m_run << ", Event: " << m_globalEvent << ", Z: " << (Int_t)(m_zvtx*10)/10. << " cm, Jet p_{T} = " <<  (Int_t)leadPt << " GeV; Tower Index #phi; Tower Index #eta";
0615 
0616     auto h2 = new TH2F(name.str().c_str(), title.str().c_str(), m_bins_phi, -0.5, m_bins_phi-0.5, m_bins_eta, -0.5, m_bins_eta-0.5);
0617 
0618     name.str("");
0619     title.str("");
0620 
0621     name << "h2TowerEnergySub_" << m_run << "_" << m_globalEvent << "_" << (Int_t)leadPt;
0622     title << "Tower Energy, Run: " << m_run << ", Event: " << m_globalEvent << ", Z: " << (Int_t)(m_zvtx*10)/10. << " cm, Jet p_{T} = " <<  (Int_t)leadPt << " GeV; Tower Index #phi; Tower Index #eta";
0623 
0624     auto h2Sub = new TH2F(name.str().c_str(), title.str().c_str(), m_bins_phi, -0.5, m_bins_phi-0.5, m_bins_eta, -0.5, m_bins_eta-0.5);
0625 
0626     for(UInt_t towerIndex = 0; towerIndex < towersCEMC->size(); ++towerIndex) {
0627       UInt_t key       = TowerInfoDefs::encode_hcal(towerIndex);
0628       UInt_t iphi      = TowerInfoDefs::getCaloTowerPhiBin(key);
0629       UInt_t ieta      = TowerInfoDefs::getCaloTowerEtaBin(key);
0630 
0631       h2->SetBinContent(iphi+1,ieta+1, towerEnergy[towerIndex]);
0632       h2Sub->SetBinContent(iphi+1,ieta+1, towerEnergySub[towerIndex]);
0633     }
0634 
0635     h2TowerEnergy.push_back(h2);
0636     h2TowerEnergySub.push_back(h2Sub);
0637   }
0638 
0639   if(m_saveTree && hasHighPtJet) {
0640     // loop over base towers
0641     for(UInt_t towerIndex = 0; towerIndex < towersCEMCBase->size(); ++towerIndex) {
0642       TowerInfo* tower = towersCEMCBase->get_tower_at_channel(towerIndex);
0643       m_towersCEMCBase_isGood.push_back(tower->get_isGood());
0644       m_towersCEMCBase_time.push_back(tower->get_time_float());
0645       m_towersCEMCBase_energy.push_back(tower->get_energy());
0646     }
0647   }
0648 
0649   // R = 0.6
0650   for (auto jet : *jets_r06) {
0651     Float_t jet_pt = jet->get_pt();
0652     if (jet_pt < m_lowPtThreshold || fabs(jet->get_eta()) > m_eta_high-0.6) continue;  // to remove noise jets
0653 
0654     m_id_r06.push_back(jet->get_id());
0655     m_nComponent_r06.push_back(jet->size_comp());
0656     m_eta_r06.push_back(jet->get_eta());
0657     m_phi_r06.push_back(jet->get_phi());
0658     m_e_r06.push_back(jet->get_e());
0659     m_pt_r06.push_back(jet_pt);
0660 
0661     for (UInt_t i = 0; i < eventTrig_status.size(); ++i) {
0662       if (eventTrig_status[i]) {
0663         hJetPt_r06[i]->Fill(jet_pt);
0664         if(m_hasBkg) hJetPt_r06_bkg[i]->Fill(jet_pt);
0665         else         hJetPt_r06_nobkg[i]->Fill(jet_pt);
0666       }
0667     }
0668 
0669     ++m_nJet_r06;
0670   }
0671 
0672   if(m_saveTree) {
0673     // fill the tree
0674     m_T->Fill();
0675   }
0676 
0677   m_nJets_r02 += m_nJet_r02;
0678   m_nJets_r04 += m_nJet_r04;
0679   m_nJets_r06 += m_nJet_r06;
0680 
0681   return Fun4AllReturnCodes::EVENT_OK;
0682 }
0683 
0684 //____________________________________________________________________________..
0685 Int_t JetValidation::ResetEvent(PHCompositeNode *topNode)
0686 {
0687   m_id_r02.clear();
0688   m_nComponent_r02.clear();
0689   m_eta_r02.clear();
0690   m_phi_r02.clear();
0691   m_e_r02.clear();
0692   m_pt_r02.clear();
0693 
0694   m_id_r04.clear();
0695   m_nComponent_r04.clear();
0696   m_eta_r04.clear();
0697   m_phi_r04.clear();
0698   m_e_r04.clear();
0699   m_pt_r04.clear();
0700 
0701   m_id_r06.clear();
0702   m_nComponent_r06.clear();
0703   m_eta_r06.clear();
0704   m_phi_r06.clear();
0705   m_e_r06.clear();
0706   m_pt_r06.clear();
0707 
0708   m_scaledVector.clear();
0709 
0710   m_towersCEMCBase_isGood.clear();
0711   m_towersCEMCBase_energy.clear();
0712   m_towersCEMCBase_time.clear();
0713 
0714   m_towersCEMC_isGood.clear();
0715   m_towersIHCal_isGood.clear();
0716   m_towersOHCal_isGood.clear();
0717 
0718   m_towersCEMC_energy.clear();
0719   m_towersIHCal_energy.clear();
0720   m_towersOHCal_energy.clear();
0721 
0722   m_towersCEMCSub_energy.clear();
0723   m_towersIHCalSub_energy.clear();
0724   m_towersOHCalSub_energy.clear();
0725 
0726   return Fun4AllReturnCodes::EVENT_OK;
0727 }
0728 
0729 //____________________________________________________________________________..
0730 Int_t JetValidation::End(PHCompositeNode *topNode)
0731 {
0732   cout << "Events" << endl;
0733   cout << "All: " << hEvents->GetBinContent((Int_t)EventStatus::ALL+1) << endl;
0734   cout << endl;
0735 
0736   cout << "ZVTX: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ALL+1) << " %" << endl;
0737   cout << "ZVTX60: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << " %" << endl;
0738   cout << "ZVTX50: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX50+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX50+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << " %" << endl;
0739   cout << "ZVTX30: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX30+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX30+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << " %" << endl;
0740   cout << "ZVTX20: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX20+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX20+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << " %" << endl;
0741   cout << "ZVTX10: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX10+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX10+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << " %" << endl;
0742   cout << endl;
0743 
0744   cout << "ALL_MBDNS1: " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ALL+1) << " %" << endl;
0745   cout << "ALL_MBDNS1_JET8: " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1_JET8+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1_JET8+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1+1) << " %" << endl;
0746   cout << "ALL_MBDNS1_JET10: " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1_JET10+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1_JET10+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1+1) << " %" << endl;
0747   cout << "ALL_MBDNS1_JET12: " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1_JET12+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1_JET12+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ALL_MBDNS1+1) << " %" << endl;
0748   cout << endl;
0749 
0750   cout << "ZVTX60_MBDNS1: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0751   cout << "ZVTX60_MBDNS1_JET8: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1_JET8+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1_JET8+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0752   cout << "ZVTX60_MBDNS1_JET10: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1_JET10+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1_JET10+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0753   cout << "ZVTX60_MBDNS1_JET12: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1_JET12+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_MBDNS1_JET12+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0754   cout << endl;
0755 
0756   cout << "ZVTX60: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX+1) << " %" << endl;
0757   cout << "ZVTX60_JET8: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_JET8+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_JET8+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0758   cout << "ZVTX60_JET10: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_JET10+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_JET10+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0759   cout << "ZVTX60_JET12: " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_JET12+1) << ", " << hEvents->GetBinContent((Int_t)EventStatus::ZVTX60_JET12+1)*100./hEvents->GetBinContent((Int_t)EventStatus::ZVTX60+1) << " %" << endl;
0760   cout << endl;
0761 
0762   cout << "Jets" << endl;
0763   cout << "R=0.2: " << m_nJets_r02 << endl;
0764   cout << "R=0.4: " << m_nJets_r04 << endl;
0765   cout << "R=0.4 with pT >= " << m_highPtThreshold << ": " << m_highPtJetCtr << endl;
0766   cout << "R=0.6: " << m_nJets_r06 << endl;
0767 
0768   if(m_saveTree) {
0769     cout << "JetValidation::End - Output to " << m_outputTreeFileName << endl;
0770     m_outputTreeFile->cd();
0771 
0772     m_T->Write();
0773 
0774     m_outputTreeFile->Close();
0775     delete m_outputTreeFile;
0776   }
0777   cout << "JetValidation::End - Output to " << m_outputQAFileName << endl;
0778 
0779   m_outputQAFile->cd();
0780 
0781   hEvents->Write();
0782   hEventsBkg->Write();
0783   hZVtx->Write();
0784 
0785   m_outputQAFile->mkdir("hJet_r02/pt");
0786   m_outputQAFile->mkdir("hJet_r04/pt");
0787   m_outputQAFile->mkdir("hJet_r06/pt");
0788 
0789   m_outputQAFile->mkdir("hJet_r02/deltaPhi");
0790   m_outputQAFile->mkdir("hJet_r04/deltaPhi");
0791   m_outputQAFile->mkdir("hJet_r06/deltaPhi");
0792 
0793   m_outputQAFile->mkdir("hJet_r02/eta_phi");
0794   m_outputQAFile->mkdir("hJet_r04/eta_phi");
0795   m_outputQAFile->mkdir("hJet_r06/eta_phi");
0796 
0797   m_outputQAFile->mkdir("TowerEnergy/base");
0798   m_outputQAFile->mkdir("TowerEnergy/sub");
0799 
0800   for(UInt_t i = 0; i < JetEvent_Status_vec.size(); ++i) {
0801     m_outputQAFile->cd("hJet_r02/pt");
0802     hJetPt_r02[i]->Write();
0803     hJetPt_r02_bkg[i]->Write();
0804     hJetPt_r02_nobkg[i]->Write();
0805 
0806     m_outputQAFile->cd("hJet_r04/pt");
0807     hJetPt_r04[i]->Write();
0808     hJetPt_r04_bkg[i]->Write();
0809     hJetPt_r04_nobkg[i]->Write();
0810 
0811     m_outputQAFile->cd("hJet_r06/pt");
0812     hJetPt_r06[i]->Write();
0813     hJetPt_r06_bkg[i]->Write();
0814     hJetPt_r06_nobkg[i]->Write();
0815 
0816     m_outputQAFile->cd("hJet_r04/deltaPhi");
0817     hJetDeltaPhi_r04[i]->Write();
0818     hJetDeltaPhi_r04_bkg[i]->Write();
0819     hJetDeltaPhi_r04_nobkg[i]->Write();
0820 
0821     m_outputQAFile->cd("hJet_r04/eta_phi");
0822     h2JetEtaPhi_r04[i]->Write();
0823     h2JetEtaPhi_r04_bkg[i]->Write();
0824     h2JetEtaPhi_r04_nobkg[i]->Write();
0825   }
0826 
0827   for(UInt_t i = 0; i < h2TowerEnergy.size(); ++i) {
0828     m_outputQAFile->cd("TowerEnergy/base");
0829     h2TowerEnergy[i]->Write();
0830 
0831     m_outputQAFile->cd("TowerEnergy/sub");
0832     h2TowerEnergySub[i]->Write();
0833   }
0834 
0835   m_outputQAFile->Close();
0836   delete m_outputQAFile;
0837 
0838   cout << "JetValidation::End(PHCompositeNode *topNode) This is the End..." << endl;
0839   return Fun4AllReturnCodes::EVENT_OK;
0840 }
0841 
0842 
0843 //____________________________________________________________________________..
0844 Bool_t JetValidation::isBackgroundEvent(std::vector<Float_t> &towerEnergy) {
0845   UInt_t ctr = 0;
0846 
0847   for(UInt_t i = 0; i < m_bins_phi; ++i) {
0848     for(UInt_t j = 0; j < m_bins_eta; ++j) {
0849       UInt_t key        = TowerInfoDefs::encode_hcal(j,i);
0850       UInt_t towerIndex = TowerInfoDefs::decode_hcal(key);
0851       Float_t energy    = towerEnergy[towerIndex];
0852 
0853       if(energy < m_bkg_tower_energy) {
0854         ctr = 0;
0855         continue;
0856       }
0857 
0858       // left side
0859       key        = (i == 0) ? TowerInfoDefs::encode_hcal(j,m_bins_phi-1) : TowerInfoDefs::encode_hcal(j,i-1);
0860       towerIndex = TowerInfoDefs::decode_hcal(key);
0861       energy     = towerEnergy[towerIndex];
0862 
0863       if(energy >= m_bkg_tower_neighbor_energy) {
0864         ctr = 0;
0865         continue;
0866       }
0867 
0868       // right side
0869       key        = (i == m_bins_phi-1) ? TowerInfoDefs::encode_hcal(j,0) : TowerInfoDefs::encode_hcal(j,i+1);
0870       towerIndex = TowerInfoDefs::decode_hcal(key);
0871       energy     = towerEnergy[towerIndex];
0872 
0873       if(energy < m_bkg_tower_neighbor_energy){
0874         ++ctr;
0875         if(ctr >= m_bkg_towers) {
0876           cout << "Found Background: (" << i << "," << j-m_bkg_towers+1 << ")" << " to (" << i << "," << j << ")" << endl;
0877           return true;
0878         }
0879       }
0880       else ctr = 0;
0881     }
0882     ctr = 0;
0883   }
0884 
0885   return false;
0886 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!