Jet-Study/src/JetValidationv2.cc

0001 #include <JetValidationv2.h>
0002 // -- c++
0003 #include <cmath>
0004 #include <iomanip>
0005 #include <iostream>
0006 #include <sstream>
0007 #include <vector>
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 // -- flags
0022 #include <phool/recoConsts.h>
0023 // -- root
0024 #include <TFile.h>
0025 // -- Tower stuff
0026 #include <calobase/TowerInfoDefs.h>
0027 #include <calobase/TowerInfoContainer.h>
0028 #include <calobase/TowerInfo.h>
0029 // -- Trigger
0030 #include <calotrigger/TriggerRunInfo.h>
0031
0032 // Jet Utils
0033 #include "JetUtils.h"
0034
0035 using std::cout;
0036 using std::endl;
0037 using std::string;
0038 using std::to_string;
0039 using std::stringstream;
0040 using std::min;
0041 using std::max;
0042
0043 //____________________________________________________________________________..
0044 JetValidationv2::JetValidationv2()
0045   : SubsysReco("JetValidationv2")
0046   , m_recoJetName_r04("AntiKt_unsubtracted_r04")
0047   , m_emcTowerBaseNode("TOWERINFO_CALIB_CEMC")
0048   , m_emcTowerNode("TOWERINFO_CALIB_CEMC_RETOWER")
0049   , m_ihcalTowerNode("TOWERINFO_CALIB_HCALIN")
0050   , m_ohcalTowerNode("TOWERINFO_CALIB_HCALOUT")
0051   , m_triggerBit(17) /*Jet 8 GeV + MBD NS >= 1*/
0052   , m_triggerBits(42)
0053   , m_zvtx_max(30) /*cm*/
0054   , m_bins_phi(64)
0055   , m_bins_phi_cemc(256)
0056   , m_phi_low(-M_PI)
0057   , m_phi_high(M_PI)
0058   , m_bins_eta(24)
0059   , m_bins_eta_cemc(96)
0060   , m_eta_low(-1.1)
0061   , m_eta_high(1.1)
0062   , m_bins_pt(400)
0063   , m_pt_low(0) /*GeV*/
0064   , m_pt_high(200) /*GeV*/
0065   , m_pt_background(60) /*GeV*/
0066   , m_bins_zvtx(400)
0067   , m_zvtx_low(-100)
0068   , m_zvtx_high(100)
0069   , m_bins_frac(140)
0070   , m_frac_low(-0.2)
0071   , m_frac_high(1.2)
0072   , m_bins_ET(190)
0073   , m_ET_low(10)
0074   , m_ET_high(200)
0075   , m_bins_constituents(300)
0076   , m_constituents_low(0)
0077   , m_constituents_high(300)
0078   , m_bins_nJets(100)
0079   , m_nJets_low(0)
0080   , m_nJets_high(100)
0081   , m_event(0)
0082   , m_R(0.4)
0083   , m_nJets_min(9999)
0084   , m_nJets_max(0)
0085   , m_constituents_min(9999)
0086   , m_constituents_max(0)
0087   , m_pt_min(9999)
0088   , m_pt_max(0)
0089 {
0090   cout << "JetValidationv2::JetValidationv2(const std::string &name) Calling ctor" << endl;
0091 }
0092
0093 //____________________________________________________________________________..
0094 JetValidationv2::~JetValidationv2()
0095 {
0096   cout << "JetValidationv2::~JetValidationv2() Calling dtor" << endl;
0097 }
0098
0099 //____________________________________________________________________________..
0100 Int_t JetValidationv2::Init(PHCompositeNode *topNode)
0101 {
0102   cout << "JetValidationv2::Init(PHCompositeNode *topNode) Initializing" << endl;
0103
0104   // so that the histos actually get written out
0105   Fun4AllServer *se = Fun4AllServer::instance();
0106   se->Print("NODETREE");
0107
0108   hEvents = new TH1F("hEvents","Events; Status; Counts", m_eventStatus.size(), 0, m_eventStatus.size());
0109
0110   for(UInt_t i = 0; i < m_eventStatus.size(); ++i) {
0111     hEvents->GetXaxis()->SetBinLabel(i+1, m_eventStatus[i].c_str());
0112   }
0113
0114   hzvtxAll = new TH1F("zvtxAll","Z Vertex; z [cm]; Counts", m_bins_zvtx, m_zvtx_low, m_zvtx_high);
0115
0116   hjetPhiEtaPt = new TH3F("hjetPhiEtaPt", "Jet; #phi; #eta; p_{T} [GeV]", m_bins_phi, m_phi_low, m_phi_high
0117                                                                         , m_bins_eta, m_eta_low, m_eta_high
0118                                                                         , m_bins_pt, m_pt_low, m_pt_high);
0119
0120   hjetConstituentsVsPt = new TH2F("hjetConstituentsVsPt", "Jet; p_{T} [GeV]; Constituents", m_bins_pt, m_pt_low, m_pt_high
0121                                                             , m_bins_constituents, m_constituents_low, m_constituents_high);
0122
0123   hNJetsVsLeadPt = new TH2F("hNJetsVsLeadPt", "Event; Lead p_{T} [GeV]; # of Jets", m_bins_pt, m_pt_low, m_pt_high
0124                                                                                   , m_bins_nJets, m_nJets_low, m_nJets_high);
0125
0126   h2ETVsFracCEMC = new TH2F("h2ETVsFracCEMC","Jet; Fraction of E_{T,Lead Jet} EMCal; E_{T,Lead Jet} [GeV]"
0127                             , m_bins_frac, m_frac_low, m_frac_high, m_bins_ET, m_ET_low, m_ET_high);
0128
0129   h2ETVsFracCEMC_miss = new TH2F("h2ETVsFracCEMC_miss","Jet; Fraction of E_{T,Lead Jet} EMCal; E_{T,Lead Jet} [GeV]"
0130                                  , m_bins_frac, m_frac_low, m_frac_high, m_bins_ET, m_ET_low, m_ET_high);
0131
0132   h2FracOHCalVsFracCEMC = new TH2F("h2FracOHCalVsFracCEMC","Jet; Fraction of E_{T,Lead Jet} EMCal; Fraction of E_{T,Lead Jet} OHCal"
0133                                    , m_bins_frac, m_frac_low, m_frac_high, m_bins_frac, m_frac_low, m_frac_high);
0134
0135   h2FracOHCalVsFracCEMC_miss = new TH2F("h2FracOHCalVsFracCEMC_miss","Jet; Fraction of E_{T,Lead Jet} EMCal; Fraction of E_{T,Lead Jet} OHCal"
0136                                    , m_bins_frac, m_frac_low, m_frac_high, m_bins_frac, m_frac_low, m_frac_high);
0137
0138   m_triggeranalyzer = new TriggerAnalyzer();
0139
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143 //____________________________________________________________________________..
0144 Int_t JetValidationv2::process_event(PHCompositeNode *topNode)
0145 {
0146   EventHeader* eventInfo = findNode::getClass<EventHeader>(topNode,"EventHeader");
0147   if(!eventInfo)
0148   {
0149     cout << PHWHERE << "JetValidationv2::process_event - Fatal Error - EventHeader node is missing. " << endl;
0150     return Fun4AllReturnCodes::ABORTEVENT;
0151   }
0152
0153   Int_t m_globalEvent = eventInfo->get_EvtSequence();
0154   Int_t m_run         = eventInfo->get_RunNumber();
0155
0156   if (m_event % 20 == 0) cout << "Progress: " << m_event << ", Global: " << m_globalEvent << endl;
0157   ++m_event;
0158
0159   TriggerRunInfo* triggerruninfo = findNode::getClass<TriggerRunInfo>(topNode, "TriggerRunInfo");
0160   if (!triggerruninfo) {
0161     cout << "JetValidationv2::process_event - Error can not find TriggerRunInfo node " << endl;
0162     return Fun4AllReturnCodes::ABORTRUN;
0163   }
0164   if(m_event == 1) {
0165     triggerruninfo->identify();
0166   }
0167
0168   // zvertex
0169   Float_t m_zvtx = -9999;
0170   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0171
0172   if (!vertexmap) {
0173     cout << "JetValidationv2::process_event - Error can not find global vertex node " << endl;
0174     return Fun4AllReturnCodes::ABORTRUN;
0175   }
0176
0177   if(!vertexmap->empty()) {
0178     GlobalVertex* vtx = vertexmap->begin()->second;
0179     m_zvtx = vtx->get_z();
0180   }
0181
0182   hzvtxAll->Fill(m_zvtx);
0183   hEvents->Fill(m_status::ALL);
0184
0185   // skip event if zvtx is too large
0186   if(fabs(m_zvtx) >= m_zvtx_max) {
0187     return Fun4AllReturnCodes::ABORTEVENT;
0188   }
0189
0190   hEvents->Fill(m_status::ZVTX30);
0191
0192   m_triggeranalyzer->decodeTriggers(topNode);
0193
0194   // skip event if it did not fire the desired trigger
0195   if(!m_triggeranalyzer->didTriggerFire(m_triggerBit)) {
0196     return Fun4AllReturnCodes::ABORTEVENT;
0197   }
0198
0199   vector<Int_t> triggerIdx;
0200   for(Int_t i = 0; i < m_triggerBits; ++i) {
0201     if(m_triggeranalyzer->didTriggerFire(i)) {
0202       triggerIdx.push_back(i);
0203     }
0204   }
0205
0206   // interface to reco jets
0207   JetContainer* jets_r04 = findNode::getClass<JetContainer>(topNode, m_recoJetName_r04);
0208   if (!jets_r04) {
0209     cout << "JetValidationv2::process_event - Error can not find DST Reco JetContainer node " << m_recoJetName_r04 << endl;
0210     return Fun4AllReturnCodes::ABORTRUN;
0211   }
0212
0213   Int_t jetPtLead    = 0;
0214   Int_t jetPtSubLead = 0;
0215
0216   Float_t jetPhiLead = 0;
0217   Float_t jetPhiSubLead = 0;
0218
0219   Float_t jetEtaLead = 0;
0220   Float_t jetEtaSubLead = 0;
0221
0222   Int_t nJets = 0;
0223
0224   // trigger
0225   Bool_t hasBkg = false;
0226   for (auto jet : *jets_r04) {
0227     Float_t phi = jet->get_phi();
0228     Float_t eta = jet->get_eta();
0229     Float_t pt = jet->get_pt();
0230     Int_t constituents = jet->get_comp_vec().size();
0231
0232     // exclude jets near the edge of the detector
0233     if(JetUtils::check_bad_jet_eta(eta, m_zvtx, m_R)) continue;
0234     ++nJets;
0235
0236     hjetPhiEtaPt->Fill(phi, eta, pt);
0237
0238     hjetConstituentsVsPt->Fill(pt, constituents);
0239
0240     m_constituents_min = min(m_constituents_min, constituents);
0241     m_constituents_max = max(m_constituents_max, constituents);
0242
0243     m_pt_min = min(m_pt_min, (Int_t)pt);
0244     m_pt_max = max(m_pt_max, (Int_t)pt);
0245
0246     if(pt >= m_pt_background) {
0247       hasBkg = true;
0248     }
0249
0250     if((Int_t)pt > jetPtSubLead) {
0251       if((Int_t)pt > jetPtLead) {
0252         jetPtSubLead = jetPtLead;
0253         jetPtLead = pt;
0254
0255         jetPhiSubLead = jetPhiLead;
0256         jetPhiLead = phi;
0257
0258         jetEtaSubLead = jetEtaLead;
0259         jetEtaLead = eta;
0260       }
0261       else {
0262         jetPtSubLead  = pt;
0263         jetPhiSubLead = phi;
0264         jetEtaSubLead = eta;
0265       }
0266     }
0267   }
0268
0269   hNJetsVsLeadPt->Fill(jetPtLead, nJets);
0270
0271   m_nJets_min = min(m_nJets_min, nJets);
0272   m_nJets_max = max(m_nJets_max, nJets);
0273
0274   recoConsts* rc = recoConsts::instance();
0275
0276   Bool_t failsLoEmJetCut = rc->get_IntFlag("failsLoEmJetCut");
0277   Bool_t failsHiEmJetCut = rc->get_IntFlag("failsHiEmJetCut");
0278   Bool_t failsIhJetCut   = rc->get_IntFlag("failsIhJetCut");
0279   Bool_t failsAnyJetCut  = rc->get_IntFlag("failsAnyJetCut");
0280
0281   Bool_t isDijet   = rc->get_IntFlag("isDijet");
0282   Float_t frcem    = rc->get_FloatFlag("frcem");
0283   Float_t frcoh    = rc->get_FloatFlag("frcoh");
0284   Float_t maxJetET = rc->get_FloatFlag("maxJetET");
0285   Float_t dPhi     = rc->get_FloatFlag("dPhi");
0286
0287   h2ETVsFracCEMC->Fill(frcem, maxJetET);
0288
0289   // if event doesn't have a high pt jet then skip to the next event
0290   if(!hasBkg) return Fun4AllReturnCodes::ABORTEVENT;
0291
0292   h2FracOHCalVsFracCEMC->Fill(frcem, frcoh);
0293
0294   cout << "Background Jet: " << jetPtLead << " GeV, Event: " << m_globalEvent << endl;
0295
0296   cout << "isDijet: " << isDijet
0297        << ", frcem: " << frcem
0298        << ", frcoh: " << frcoh
0299        << ", maxJetET: " << maxJetET
0300        << ", dPhi: " << dPhi << endl;
0301
0302   hEvents->Fill(m_status::ZVTX30_BKG);
0303
0304   if(failsLoEmJetCut) {
0305     hEvents->Fill(m_status::ZVTX30_BKG_failsLoEmJetCut);
0306   }
0307   if(failsHiEmJetCut) {
0308     hEvents->Fill(m_status::ZVTX30_BKG_failsHiEmJetCut);
0309   }
0310   if(failsIhJetCut) {
0311     hEvents->Fill(m_status::ZVTX30_BKG_failsIhJetCut);
0312   }
0313   if(!failsAnyJetCut) {
0314     hEvents->Fill(m_status::ZVTX30_BKG_failsNoneJetCut);
0315   }
0316   else {
0317     // okay to skip the event since the background jet is correctly flagged
0318     return Fun4AllReturnCodes::ABORTEVENT;
0319   }
0320
0321   h2ETVsFracCEMC_miss->Fill(frcem, maxJetET);
0322   h2FracOHCalVsFracCEMC_miss->Fill(frcem, frcoh);
0323
0324   cout << "Background Jet Not Flagged!" << endl;
0325
0326   // round nearest 0.1
0327   jetPhiLead    = (Int_t)(jetPhiLead*10)/10.;
0328   jetPhiSubLead = (Int_t)(jetPhiSubLead*10)/10.;
0329
0330   jetEtaLead    = (Int_t)(jetEtaLead*10)/10.;
0331   jetEtaSubLead = (Int_t)(jetEtaSubLead*10)/10.;
0332
0333   // round the bkg check variables
0334   frcem = (Int_t)(frcem*100)/100.;
0335   frcoh = (Int_t)(frcoh*100)/100.;
0336   dPhi  = (Int_t)(dPhi*1000)/1000.;
0337
0338   // Get TowerInfoContainer
0339   TowerInfoContainer* towersCEMCBase = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerBaseNode.c_str());
0340   TowerInfoContainer* towersCEMC     = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0341   TowerInfoContainer* towersIHCal    = findNode::getClass<TowerInfoContainer>(topNode, m_ihcalTowerNode.c_str());
0342   TowerInfoContainer* towersOHCal    = findNode::getClass<TowerInfoContainer>(topNode, m_ohcalTowerNode.c_str());
0343
0344   if (!towersCEMCBase || !towersCEMC || !towersIHCal || !towersOHCal) {
0345     cout << PHWHERE << "JetValidation::process_event Could not find one of "
0346          << m_emcTowerBaseNode  << ", "
0347          << m_emcTowerNode      << ", "
0348          << m_ihcalTowerNode    << ", "
0349          << m_ohcalTowerNode    << ", "
0350          << endl;
0351     return Fun4AllReturnCodes::ABORTEVENT;
0352   }
0353
0354   stringstream name;
0355   stringstream nameSuffix;
0356   stringstream title;
0357   stringstream titleSuffix;
0358
0359   titleSuffix << "Run: " << m_run << ", Event: " << m_globalEvent
0360               << ", Jet 1 p_{T}: " << jetPtLead << " GeV, (" << jetPhiLead << "," << jetEtaLead << ")"
0361               << ", Jet 2 p_{T}: " << jetPtSubLead << " GeV, (" << jetPhiSubLead << "," << jetEtaSubLead << ")"
0362               << "; #phi; #eta";
0363
0364   string s_triggerIdx = "";
0365
0366   for(auto idx : triggerIdx) {
0367     if(!s_triggerIdx.empty()) s_triggerIdx += "-";
0368     s_triggerIdx += to_string(idx);
0369   }
0370
0371   nameSuffix << m_run << "_" << m_globalEvent << "_" << s_triggerIdx << "_" << (Int_t)m_zvtx << "_" << jetPtLead << "_" << jetPtSubLead
0372             << "_" << frcem << "_" << frcoh << "_" << (Int_t)maxJetET << "_" << dPhi << "_" << isDijet;
0373
0374   name << "hCEMCBase_" << nameSuffix.str();
0375   title << "CEMC Base: " << titleSuffix.str();
0376
0377   auto hCEMCBase_ = new TH2F(name.str().c_str(), title.str().c_str(), m_bins_phi_cemc, m_phi_low, m_phi_high, m_bins_eta_cemc, m_eta_low, m_eta_high);
0378
0379   name.str("");
0380   title.str("");
0381   name << "hCEMC_" << nameSuffix.str();
0382   title << "CEMC: " << titleSuffix.str();
0383
0384   auto hCEMC_ = 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);
0385
0386   name.str("");
0387   title.str("");
0388   name << "hIHCal_" << nameSuffix.str();
0389   title << "IHCal: " << titleSuffix.str();
0390
0391   auto hIHCal_ = 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);
0392
0393   name.str("");
0394   title.str("");
0395   name << "hOHCal_" << nameSuffix.str();
0396   title << "OHCal: " << titleSuffix.str();
0397
0398   auto hOHCal_ = 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);
0399
0400   // loop over CEMC towers
0401   for(UInt_t towerIndex = 0; towerIndex < towersCEMCBase->size(); ++towerIndex) {
0402     UInt_t key  = TowerInfoDefs::encode_emcal(towerIndex);
0403     UInt_t iphi = TowerInfoDefs::getCaloTowerPhiBin(key);
0404     UInt_t ieta = TowerInfoDefs::getCaloTowerEtaBin(key);
0405
0406     // map: [0, 127]  -> [128,255], [0,pi]
0407     // map: [128, 255] -> [0,127], [-pi,0]
0408     iphi = (iphi + m_bins_phi_cemc/2) % m_bins_phi_cemc;
0409
0410     // EMCal
0411     TowerInfo* tower = towersCEMCBase->get_tower_at_channel(towerIndex);
0412     Float_t energy = (tower->get_isGood()) ? tower->get_energy() : 0;
0413
0414     hCEMCBase_->SetBinContent(iphi+1, ieta+1, energy);
0415   }
0416
0417   // loop over towers
0418   for(UInt_t towerIndex = 0; towerIndex < towersCEMC->size(); ++towerIndex) {
0419     UInt_t key  = TowerInfoDefs::encode_hcal(towerIndex);
0420     UInt_t iphi = TowerInfoDefs::getCaloTowerPhiBin(key);
0421     UInt_t ieta = TowerInfoDefs::getCaloTowerEtaBin(key);
0422
0423     // map: [0, 31]  -> [32,63], [0,pi]
0424     // map: [32, 63] -> [0,31], [-pi,0]
0425     iphi = (iphi + m_bins_phi/2) % m_bins_phi;
0426
0427     // EMCal
0428     TowerInfo* tower = towersCEMC->get_tower_at_channel(towerIndex);
0429     Float_t energy = (tower->get_isGood()) ? tower->get_energy() : 0;
0430
0431     hCEMC_->SetBinContent(iphi+1, ieta+1, energy);
0432
0433     // HCal
0434     tower = towersIHCal->get_tower_at_channel(towerIndex);
0435     energy = (tower->get_isGood()) ? tower->get_energy() : 0;
0436
0437     hIHCal_->SetBinContent(iphi+1, ieta+1, energy);
0438
0439     tower = towersOHCal->get_tower_at_channel(towerIndex);
0440     energy = (tower->get_isGood()) ? tower->get_energy() : 0;
0441
0442     hOHCal_->SetBinContent(iphi+1, ieta+1, energy);
0443   }
0444
0445   hCEMCBase.push_back(hCEMCBase_);
0446   hCEMC.push_back(hCEMC_);
0447   hIHCal.push_back(hIHCal_);
0448   hOHCal.push_back(hOHCal_);
0449
0450   return Fun4AllReturnCodes::EVENT_OK;
0451 }
0452
0453 //____________________________________________________________________________..
0454 Int_t JetValidationv2::ResetEvent(PHCompositeNode *topNode)
0455 {
0456   return Fun4AllReturnCodes::EVENT_OK;
0457 }
0458
0459 //____________________________________________________________________________..
0460 Int_t JetValidationv2::End(PHCompositeNode *topNode)
0461 {
0462   cout << "JetValidationv2::End(PHCompositeNode *topNode) This is the End..." << endl;
0463   cout << "Events Summary" << endl;
0464   for(UInt_t i = 0; i < m_eventStatus.size(); ++i) {
0465     cout << m_eventStatus[i] << ": " << hEvents->GetBinContent(i+1) << " Events" << endl;
0466   }
0467   cout << "=======================" << endl;
0468   cout << "Jets Summary" << endl;
0469   cout << "Constituents min: " << m_constituents_min << ", max: " << m_constituents_max << endl;
0470   cout << "nJets min: " << m_nJets_min << ", max: " << m_nJets_max << endl;
0471   cout << "Jet pT min: " << m_pt_min << " GeV, max: " << m_pt_max << " GeV" << endl;
0472
0473   TFile output(m_outputFile.c_str(),"recreate");
0474
0475   output.cd();
0476
0477   output.mkdir("event");
0478   output.mkdir("zvtx");
0479   output.mkdir("jets");
0480   output.mkdir("bkg_checks");
0481
0482   output.mkdir("CEMCBase");
0483   output.mkdir("CEMC");
0484   output.mkdir("IHCal");
0485   output.mkdir("OHCal");
0486
0487   output.cd("event");
0488   hEvents->Write();
0489
0490   output.cd("zvtx");
0491   hzvtxAll->Write();
0492
0493   output.cd("jets");
0494   hjetPhiEtaPt->Write();
0495   hjetConstituentsVsPt->Write();
0496   hNJetsVsLeadPt->Write();
0497
0498   output.cd("bkg_checks");
0499   h2ETVsFracCEMC->Write();
0500   h2FracOHCalVsFracCEMC->Write();
0501   h2ETVsFracCEMC_miss->Write();
0502   h2FracOHCalVsFracCEMC_miss->Write();
0503
0504   for(UInt_t i = 0; i < hCEMC.size(); ++i) {
0505     output.cd("CEMCBase");
0506     hCEMCBase[i]->Write();
0507
0508     output.cd("CEMC");
0509     hCEMC[i]->Write();
0510
0511     output.cd("IHCal");
0512     hIHCal[i]->Write();
0513
0514     output.cd("OHCal");
0515     hOHCal[i]->Write();
0516   }
0517
0518   output.Close();
0519
0520   cout << "JetValidationv2::End(PHCompositeNode *topNode) This is the End..." << endl;
0521   return Fun4AllReturnCodes::EVENT_OK;
0522 }
0523