14 using namespace fastjet;
18 namespace SColdQcdCorrelatorAnalysis {
22 void SEnergyCorrelator::InitializeMembers() {
25 if (m_inDebugMode) PrintDebug(0);
28 m_inFileNames.clear();
29 m_eecLongSide.clear();
30 m_subEvtsToUse.clear();
31 m_jetCstVector.clear();
32 m_outHistVarDrAxis.clear();
33 m_outHistErrDrAxis.clear();
34 m_outHistVarLnDrAxis.clear();
35 m_outHistErrLnDrAxis.clear();
42 void SEnergyCorrelator::InitializeTree() {
45 if (m_inDebugMode) PrintDebug(4);
53 m_inChain -> SetMakeClass(1);
56 if (m_isInputTreeTruth) {
57 m_inChain ->
SetBranchAddress(
"Parton3_ID", &m_partonID[0], &m_brPartonID[0]);
58 m_inChain ->
SetBranchAddress(
"Parton4_ID", &m_partonID[1], &m_brPartonID[1]);
59 m_inChain ->
SetBranchAddress(
"Parton3_MomX", &m_partonMomX[0], &m_brPartonMomX[0]);
60 m_inChain ->
SetBranchAddress(
"Parton3_MomY", &m_partonMomY[0], &m_brPartonMomY[0]);
61 m_inChain ->
SetBranchAddress(
"Parton3_MomZ", &m_partonMomZ[0], &m_brPartonMomZ[0]);
62 m_inChain ->
SetBranchAddress(
"Parton4_MomX", &m_partonMomX[1], &m_brPartonMomX[1]);
63 m_inChain ->
SetBranchAddress(
"Parton4_MomY", &m_partonMomY[1], &m_brPartonMomY[1]);
64 m_inChain ->
SetBranchAddress(
"Parton4_MomZ", &m_partonMomZ[1], &m_brPartonMomZ[1]);
70 m_inChain ->
SetBranchAddress(
"EvtSumECalEne", &m_evtSumECal, &m_brEvtSumECal);
71 m_inChain ->
SetBranchAddress(
"EvtSumHCalEne", &m_evtSumHCal, &m_brEvtSumHCal);
95 if (m_inStandaloneMode) PrintMessage(2);
102 void SEnergyCorrelator::InitializeHists() {
105 if (m_inDebugMode) PrintDebug(5);
107 for (
size_t iPtBin = 0; iPtBin < m_nBinsJetPt; iPtBin++) {
108 TH1D* hInitialVarDrAxis = NULL;
109 TH1D* hInitialErrDrAxis = NULL;
110 TH1D* hInitialVarLnDrAxis = NULL;
111 TH1D* hInitialErrLnDrAxis = NULL;
112 m_outHistVarDrAxis.push_back(hInitialVarDrAxis);
113 m_outHistVarLnDrAxis.push_back(hInitialVarLnDrAxis);
114 m_outHistErrDrAxis.push_back(hInitialErrDrAxis);
115 m_outHistErrLnDrAxis.push_back(hInitialErrLnDrAxis);
119 if (m_doSecondCstLoop) {
120 vector<double> drBinEdges = m_eecLongSide[0] -> bin_edges();
121 size_t nDrBinEdges = drBinEdges.size();
123 double drBinEdgeArray[nDrBinEdges];
124 for (
size_t iDrEdge = 0; iDrEdge < nDrBinEdges; iDrEdge++) {
125 drBinEdgeArray[iDrEdge] = drBinEdges.at(iDrEdge);
127 hCstPtOneVsDr =
new TH2D(
"hCstPtOneVsDr",
"", m_nBinsDr, drBinEdgeArray, 200, 0., 100.);
128 hCstPtTwoVsDr =
new TH2D(
"hCstPtTwoVsDr",
"", m_nBinsDr, drBinEdgeArray, 200, 0., 100.);
129 hCstPtFracVsDr =
new TH2D(
"hCstPtFracVsDr",
"", m_nBinsDr, drBinEdgeArray, 500, 0., 5.);
130 hCstPhiOneVsDr =
new TH2D(
"hCstPhiOneVsDr",
"", m_nBinsDr, drBinEdgeArray, 360, -3.15, 3.15);;
131 hCstPhiTwoVsDr =
new TH2D(
"hCstPhiTwoVsDr",
"", m_nBinsDr, drBinEdgeArray, 360, -3.15, 3.15);
132 hCstEtaOneVsDr =
new TH2D(
"hCstEtaOneVsDr",
"", m_nBinsDr, drBinEdgeArray, 400, -2., 2.);
133 hCstEtaTwoVsDr =
new TH2D(
"hCstEtaTwoVsDr",
"", m_nBinsDr, drBinEdgeArray, 400, -2., 2.);
134 hDeltaPhiOneVsDr =
new TH2D(
"hDeltaPhiOneVsDr",
"", m_nBinsDr, drBinEdgeArray, 720, -6.30, 6.30);
135 hDeltaPhiTwoVsDr =
new TH2D(
"hDeltaPhiTwoVsDr",
"", m_nBinsDr, drBinEdgeArray, 720, -6.30, 6.30);
136 hDeltaEtaOneVsDr =
new TH2D(
"hDeltaEtaOneVsDr",
"", m_nBinsDr, drBinEdgeArray, 800, -4., 4.);
137 hDeltaEtaTwoVsDr =
new TH2D(
"hDeltaEtaTwoVsDr",
"", m_nBinsDr, drBinEdgeArray, 800, -4., 4.);
138 hJetPtFracOneVsDr =
new TH2D(
"hJetPtFracOneVsDr",
"", m_nBinsDr, drBinEdgeArray, 500, 0., 5.);
139 hJetPtFracTwoVsDr =
new TH2D(
"hJetPtFracTwoVsDr",
"", m_nBinsDr, drBinEdgeArray, 500, 0., 5.);
140 hCstPairWeightVsDr =
new TH2D(
"hCstPairWeightVsDr",
"", m_nBinsDr, drBinEdgeArray, 100, 0., 1.);
144 if (m_inStandaloneMode) PrintMessage(3);
151 void SEnergyCorrelator::InitializeCorrs() {
154 if (m_inDebugMode) PrintDebug(6);
157 for (
size_t iPtBin = 0; iPtBin < m_nBinsJetPt; iPtBin++) {
158 m_eecLongSide.push_back(
new contrib::eec::EECLongestSide<contrib::eec::hist::axis::log>(m_nPointCorr, m_nBinsDr, {m_drBinRange.first, m_drBinRange.second}));
162 if (m_inStandaloneMode) PrintMessage(4);
169 void SEnergyCorrelator::PrintMessage(
const uint32_t code,
const uint64_t nEvts,
const uint64_t
event) {
172 if (m_inDebugMode && (m_verbosity > 5)) PrintDebug(22);
176 cout <<
"\n Running standalone correlator calculation...\n"
177 <<
" Set name & modes:\n"
178 <<
" module name = " << m_moduleName.data() <<
"\n"
179 <<
" complex mode? = " << m_inComplexMode <<
"\n"
180 <<
" standalone mode? = " << m_inStandaloneMode <<
"\n"
181 <<
" debug mode? = " << m_inDebugMode <<
"\n"
182 <<
" batch mode? = " << m_inBatchMode
186 cout <<
" Opened files:\n"
187 <<
" output = " << m_outFileName.data() <<
"\n"
190 for (
const string&
inFileName : m_inFileNames) {
193 cout <<
" }" << endl;
196 cout <<
" Initialized input chain:\n"
197 <<
" tree name = " << m_inTreeName.data()
201 cout <<
" Initialized output histograms." << endl;
204 cout <<
" Initialized correlators." << endl;
207 cout <<
" Set correlator parameters:\n"
208 <<
" n-point = " << m_nPointCorr <<
", number of dR bins = " << m_nBinsDr <<
"\n"
209 <<
" dR bin range = (" << m_drBinRange.first <<
", " << m_drBinRange.second <<
")"
213 cout <<
" Set jet parameters:\n"
214 <<
" eta range = (" << m_etaJetRange.first <<
", " << m_etaJetRange.second <<
")\n"
215 <<
" pt range = (" << m_ptJetRange.first <<
", " << m_ptJetRange.second <<
")\n"
216 <<
" Set pTjet bins:"
218 for (uint32_t iPtBin = 0; iPtBin < m_nBinsJetPt; iPtBin++) {
219 cout <<
" bin[" << iPtBin <<
"] = (" << m_ptJetBins.at(iPtBin).first <<
", " << m_ptJetBins.at(iPtBin).second <<
")" << endl;
223 cout <<
" Beginning event loop: " << nEvts <<
" events to process..." << endl;
227 cout <<
" processing event " << (
event + 1) <<
"/" << nEvts <<
"..." << endl;
229 cout <<
" processing event " << (
event + 1) <<
"/" << nEvts <<
"...\r" <<
flush;
230 if ((event + 1) == nEvts) cout << endl;
234 cout <<
" Analysis finished!" << endl;
237 cout <<
" Saved output histograms." << endl;
240 cout <<
" Finished correlator calculation!\n" << endl;
243 cout <<
" Set constituent parameters:\n"
244 <<
" apply constituent cuts? = " << m_applyCstCuts <<
"\n"
245 <<
" momentum range = (" << m_momCstRange.first <<
", " << m_momCstRange.second <<
")\n"
246 <<
" dr range = (" << m_drCstRange.first <<
", " << m_drCstRange.second <<
")"
250 cout <<
" Finished event loop!" << endl;
253 cout <<
" Extracted output histograms from correlators." << endl;
256 cout <<
" Set which sub-events to use:" << endl;
257 switch (m_subEvtOpt) {
259 cout <<
" Option " << m_subEvtOpt <<
": use only signal event" << endl;
262 cout <<
" Option " << m_subEvtOpt <<
": use only background events" << endl;
265 cout <<
" Option " << m_subEvtOpt <<
": use only primary background event" << endl;
268 cout <<
" Option " << m_subEvtOpt <<
": use only pileup events" << endl;
271 cout <<
" Option " << m_subEvtOpt <<
": use events only with these embedding IDs: ";
272 for (
size_t iEvtToUse = 0; iEvtToUse < m_subEvtsToUse.size(); iEvtToUse++) {
273 cout << m_subEvtsToUse[iEvtToUse];
274 if ((iEvtToUse + 1) < m_subEvtsToUse.size()) {
282 cout <<
" Option " << m_subEvtOpt <<
": use everything (check what you entered)" << endl;
292 void SEnergyCorrelator::PrintDebug(
const uint32_t code) {
295 if (m_inDebugMode && (m_verbosity > 7)) {
296 cout <<
"SEnergyCorrelator::PrintDebug(uint32_t) printing a debugging statement..." << endl;
301 cout <<
"SEnergyCorrelator::InitializeMembers() initializing internal variables..." << endl;
304 cout <<
"SEnergyCorrelator::SEnergyCorrelator(string, bool, bool) calling ctor..." << endl;
307 cout <<
"SEnergyCorrelator::Init(PHCompositeNode*) initializing..." << endl;
310 cout <<
"SEnergyCorrelator::GrabInputNode() grabbing input node..." << endl;
313 cout <<
"SEnergyCorrelator::InitializeTree() initializing input tree..." << endl;
316 cout <<
"SEnergyCorrelator::InitializeHists() initializing histograms..." << endl;
319 cout <<
"SEnergyCorrelator::InitializeCorrs() initializing correlators" << endl;
322 cout <<
"SEnergyCorrelator::process_event(PHCompositeNode*) processing event..." << endl;
325 cout <<
"SEnergyCorrelator::End(PHCompositeNode*) this is the end..." << endl;
328 cout <<
"SEnergyCorrelator::SaveOutput() saving output..." << endl;
331 cout <<
"SEnergyCorrelator::Init() initializing..." << endl;
334 cout <<
"SenergyCorrelator::OpenInputFile() opening input file..." << endl;
337 cout <<
"SEnergyCorrelator::Analyze() analyzing input..." << endl;
340 cout <<
"SEnergyCorrelator::End() this is the end..." << endl;
343 cout <<
"SEnergyCorrelator::~SEnergyCorrelator() calling dtor..." << endl;
346 cout <<
"SEnergyCorrelator::OpenOutputFile() opening output file..." << endl;
349 cout <<
"SEnergyCorrelator::GetEntry(uint64_t) getting tree entry..." << endl;
352 cout <<
"SEnergyCorrelator::LoadTree(uint64_t) loading tree..." << endl;
355 cout <<
"SEnergyCorrelator::SetInputTree(string, bool) setting input tree name..." << endl;
358 cout <<
"SEnergyCorrelator::SetCorrelatorParameters(uint32_t, uint64_t, pair<double, double>) setting correlator parameters..." << endl;
361 cout <<
"SEnergyCorrelator::SetJetParameters(vector<pair<double, double>>, pair<double, double>) setting jet parameters..." << endl;
364 cout <<
"SEnergyCorrelators:CheckCriticalParameters() checking critical parameters..." << endl;
367 cout <<
"SEnergyCorrelator::PrintMessage(uint32_t, uint64_t, uint64_t) printing a message..." << endl;
370 cout <<
"SEnergyCorrelator::PrintError(uint32_t) printing an error..." << endl;
373 cout <<
"SEnergyCorrelator::SetConstituentParameters(pair<double, double>, pair<double, double>) setting constituent parameters..." << endl;
376 cout <<
"SEnergyCorrelator::ExtractHistsFromCorr() extracting output histograms..." << endl;
379 cout <<
"SEnergyCorrelator::ApplyJetCuts(double, double) applying jet cuts..." << endl;
382 cout <<
"SEnergyCorrelator::ApplyCstCuts(double, double) applying constituent cuts..." << endl;
385 cout <<
"SEnergyCorrelator::GetJetPtBin(double) getting jet pT bin..." << endl;
388 cout <<
"SEnergyCorrelator::CloseInputFile() closing input file..." << endl;
391 cout <<
"SEnergyCorrelator::CloseOutputFile() closing output file..." << endl;
394 cout <<
"SEnergyCorrelator::DoCorrelatorCalculation() looping over events and calculating correlators..." << endl;
397 cout <<
"SEnergyCorrelator::SetSubEventsToUse(uint16_t, vector<int>) setting sub-events to use..." << endl;
400 cout <<
"SEnergyCorrelator::CheckIfSubEvtGood(int) checking if sub-event is good..." << endl;
409 void SEnergyCorrelator::PrintError(
const uint32_t code,
const size_t nDrBinEdges,
const size_t iDrBin,
const string sInFileName) {
412 if (m_inDebugMode && (m_verbosity > 5)) PrintDebug(23);
416 if (m_inComplexMode) {
417 cerr <<
"SEnergyCorrelator::Init(PHCompositeNode*) PANIC: calling complex method in standalone mode! Aborting!" << endl;
419 cerr <<
"PANIC: calling complex method in standalone mode! Aborting!" << endl;
423 if (m_inComplexMode) {
424 cerr <<
"SEnergyCorrelator::GrabInputNode() PANIC: couldn't grab node \"" << m_inNodeName <<
"\"! Aborting!" << endl;
426 cerr <<
"PANIC: couldn't grab node \"" << m_inNodeName <<
"\"! Aborting!" << endl;
430 if (m_inComplexMode) {
431 cerr <<
"SEnergyCorrelator::GrabInputNode() PANIC: couldn't grab tree \"" << m_inTreeName <<
"\" from node \"" << m_inNodeName <<
"\"! Aborting!" << endl;
433 cerr <<
"PANIC: couldn't grab tree \"" << m_inTreeName <<
"\" from node \"" << m_inNodeName <<
"\"! Aborting!" << endl;
437 if (m_inComplexMode) {
438 cerr <<
"SEnergyCorrelator::process_event(PHCompositeNode*) PANIC: calling complex method in standalone mode! Aborting!" << endl;
440 cerr <<
"PANIC: calling complex method in standalone mode! Aborting!" << endl;
444 if (m_inComplexMode) {
445 cerr <<
"SEnergyCorrelator::End(PHCompositeNode*) PANIC: calling complex method in standalone mode! Aborting!" << endl;
447 cerr <<
"PANIC: calling complex method in standalone mode! Aborting!" << endl;
451 if (m_inComplexMode) {
452 cerr <<
"SEnergyCorrelator::Init() PANIC: calling standalone method in complex mode! Aborting!" << endl;
454 cerr <<
"PANIC: calling standalone method in complex mode! Aborting!" << endl;
458 if (m_inComplexMode) {
459 cerr <<
"SEnergyCorrelator::OpenInputFiles() PANIC: couldn't create input TChain! Aborting" << endl;
461 cerr <<
"PANIC: couldn't create input TChain! Aborting!" << endl;
465 if (m_inComplexMode) {
466 cerr <<
"SEnergyCorrelator::OpenInputFiles() PANIC: couldn't grab tree \"" << m_inTreeName <<
"\" from file \"" << sInFileName <<
"\"! Aborting!" << endl;
468 cerr <<
"PANIC: couldn't grab tree \"" << m_inTreeName <<
"\" from file \"" << sInFileName <<
"\"! Aborting!" << endl;
472 if (m_inComplexMode) {
473 cerr <<
"SEnergyCorrelator::Analyze() PANIC: calling standalone method in complex mode! Aborting!" << endl;
475 cerr <<
"PANIC: calling standalone method in complex mode! Aborting!" << endl;
479 if (m_inComplexMode) {
480 cerr <<
"SEnergyCorrelator::End() PANIC: calling standalone method in complex mode! Aborting!" << endl;
482 cerr <<
"PANIC: calling standalone method in complex mode! Aborting!" << endl;
486 if (m_inComplexMode) {
487 cerr <<
"SEnergyCorrelator::InitializeTree() PANIC: no TTree! Aborting!" << endl;
489 cerr <<
"PANIC: no TTree! Aborting!" << endl;
493 if (m_inComplexMode) {
494 cerr <<
"SEnergyCorrelator::OpenOutputFile() PANIC: couldn't open output file! Aborting!" << endl;
496 cerr <<
"PANIC: couldn't open output file! Aborting!" << endl;
500 if (m_inComplexMode) {
501 cerr <<
"SEnergyCorrelator::ExtraHistsFromCorr() PANIC: number of dR bin edges is no good! Aborting!" << endl;
503 cerr <<
"PANIC: number of dR bin edges is no good! Aborting!\n"
504 <<
" nDrBinEdges = " << nDrBinEdges <<
", nDrBins = " << m_nBinsDr
509 if (m_inStandaloneMode) {
510 cerr <<
"WARNING: dR bin #" << iDrBin <<
" with variance has a NAN as content or error..." << endl;
514 if (m_inStandaloneMode) {
515 cerr <<
"WARNING: dR bin #" << iDrBin <<
" with statistical error has a NAN as content or error..." << endl;
519 if (m_inComplexMode) {
520 cerr <<
"SEnergyCorrelatorFile::End() PANIC: calling standalone method in complex mode! Aborting!" << endl;
522 cerr <<
"PANIC: calling standalone method in complex mode! Aborting!" << endl;
532 bool SEnergyCorrelator::CheckCriticalParameters() {
535 if (m_inDebugMode) PrintDebug(21);
547 if (m_inDebugMode && (m_verbosity > 5)) PrintDebug(16);
549 int64_t entryStatus(-1);
553 entryStatus = m_inChain ->
GetEntry(entry);
564 if (m_inDebugMode && (m_verbosity > 5)) PrintDebug(17);
568 int64_t treeStatus(-1);
572 treeNumber = m_inChain -> GetTreeNumber();
573 treeStatus = m_inChain ->
LoadTree(entry);
577 const bool isTreeStatusGood = (treeStatus >= 0);
578 const bool isNotCurrentTree = (treeNumber != m_fCurrent);
579 if (isTreeStatusGood && isNotCurrentTree) {
580 m_fCurrent = m_inChain -> GetTreeNumber();