Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SCorrelatorJetTree.sys.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SCorrelatorJetTree.sys.h
1 // ----------------------------------------------------------------------------
2 // 'SCorrelatorJetTree.system.h'
3 // Derek Anderson
4 // 01.18.2023
5 //
6 // A module to produce a tree of jets for the sPHENIX
7 // Cold QCD Energy-Energy Correlator analysis.
8 //
9 // Derived from code by Antonio Silva (thanks!!)
10 // ----------------------------------------------------------------------------
11 
12 #pragma once
13 
14 using namespace std;
15 using namespace findNode;
16 
17 
18 
19 namespace SColdQcdCorrelatorAnalysis {
20 
21  // system methods -----------------------------------------------------------
22 
23  void SCorrelatorJetTree::InitVariables() {
24 
25  // print debug statement
26  if (m_doDebug) {
27  cout << "SCorrelatorJetTree::InitVariables() Initializing class members..." << endl;
28  }
29 
30  // initialize parton and other variables
31  m_partonID[0] = -9999;
32  m_partonID[1] = -9999;
33  m_partonMom[0] = CLHEP::Hep3Vector(-9999., -9999., -9999.);
34  m_partonMom[1] = CLHEP::Hep3Vector(-9999., -9999., -9999.);
35  m_vecEvtsToGrab.clear();
36  m_mapCstToEmbedID.clear();
37 
38  // initialize truth (inclusive) tree address members
39  m_trueVtx = CLHEP::Hep3Vector(-9999., -9999., -9999.);
40  m_trueNumJets = 0;
41  m_trueNumChrgPars = -9999;
42  m_trueSumPar = -9999.;
43  m_truePartonID[0] = -9999;
44  m_truePartonID[1] = -9999;
45  m_truePartonMomX[0] = -9999.;
46  m_truePartonMomX[1] = -9999.;
47  m_truePartonMomY[0] = -9999.;
48  m_truePartonMomY[1] = -9999.;
49  m_truePartonMomZ[0] = -9999.;
50  m_truePartonMomZ[1] = -9999.;
51  m_trueJets.clear();
52  m_trueJetNCst.clear();
53  m_trueJetID.clear();
54  m_trueJetE.clear();
55  m_trueJetPt.clear();
56  m_trueJetEta.clear();
57  m_trueJetPhi.clear();
58  m_trueJetArea.clear();
59  m_trueCstID.clear();
60  m_trueCstEmbedID.clear();
61  m_trueCstZ.clear();
62  m_trueCstDr.clear();
63  m_trueCstE.clear();
64  m_trueCstPt.clear();
65  m_trueCstEta.clear();
66  m_trueCstPhi.clear();
67 
68  // initialize reco tree address members
69  m_recoVtx = CLHEP::Hep3Vector(-9999., -9999., -9999.);
70  m_trueVtxX = -9999.;
71  m_trueVtxY = -9999.;
72  m_trueVtxZ = -9999.;
73  m_recoNumJets = 0;
74  m_recoVtxX = -9999.;
75  m_recoVtxY = -9999.;
76  m_recoVtxZ = -9999.;
77  m_recoNumTrks = -9999;
78  m_recoSumECal = -9999.;
79  m_recoSumHCal = -9999.;
80  m_recoJets.clear();
81  m_recoJetNCst.clear();
82  m_recoJetID.clear();
83  m_recoJetE.clear();
84  m_recoJetPt.clear();
85  m_recoJetEta.clear();
86  m_recoJetPhi.clear();
87  m_recoJetArea.clear();
88  m_recoCstMatchID.clear();
89  m_recoCstZ.clear();
90  m_recoCstDr.clear();
91  m_recoCstE.clear();
92  m_recoCstPt.clear();
93  m_recoCstEta.clear();
94  m_recoCstPhi.clear();
95  return;
96 
97  } // end 'InitVariables()'
98 
99 
100 
101  void SCorrelatorJetTree::InitHists() {
102 
103  // print debug statement
104  if (m_doDebug) {
105  cout << "SCorrelatorJetTree::InitHists() Initializing QA histograms..." << endl;
106  }
107 
108  // binning
109  const uint64_t nNumBins(500);
110  const uint64_t nEneBins(200);
111  const uint64_t nPhiBins(63);
112  const uint64_t nEtaBins(40);
113  const uint64_t nPtBins(200);
114  const uint64_t nDcaBins(10000);
115  const uint64_t nQualBins(200);
116  const uint64_t nDeltaBins(1000);
117  const uint64_t nAreaBins(1000);
118  const double rNumBins[CONST::NRange] = {0., 500.};
119  const double rPhiBins[CONST::NRange] = {-3.15, 3.15};
120  const double rEtaBins[CONST::NRange] = {-2., 2.};
121  const double rEneBins[CONST::NRange] = {0., 100.};
122  const double rPtBins[CONST::NRange] = {0., 100.};
123  const double rDcaBins[CONST::NRange] = {-5., 5.};
124  const double rQualBins[CONST::NRange] = {0., 20.};
125  const double rDeltaBins[CONST::NRange] = {0., 5.};
126  const double rAreaBins[CONST::NRange] = {0., 10.};
127 
128  m_outFile -> cd();
129  // no. of objects in acceptance
130  m_hNumObject[OBJECT::TRACK] = new TH1D("hNumTrks", "N_{accept} (tracks)", nNumBins, rNumBins[0], rNumBins[1]);
131  m_hNumObject[OBJECT::ECLUST] = new TH1D("hNumEClust", "N_{accept} (EMCal clust.)", nNumBins, rNumBins[0], rNumBins[1]);
132  m_hNumObject[OBJECT::HCLUST] = new TH1D("hNumHClust", "N_{accept} (HCal clust.)", nNumBins, rNumBins[0], rNumBins[1]);
133  m_hNumObject[OBJECT::FLOW] = new TH1D("hNumFlow", "N_{accept} (flow)", nNumBins, rNumBins[0], rNumBins[1]);
134  m_hNumObject[OBJECT::PART] = new TH1D("hNumPart", "N_{accept} (par.s)", nNumBins, rNumBins[0], rNumBins[1]);
135  m_hNumObject[OBJECT::TJET] = new TH1D("hNumTruthJet", "N_{jet} (truth)", nNumBins, rNumBins[0], rNumBins[1]);
136  m_hNumObject[OBJECT::RJET] = new TH1D("hNumRecoJets", "N_{jet} (reco.)", nNumBins, rNumBins[0], rNumBins[1]);
137  m_hNumObject[OBJECT::TCST] = new TH1D("hNumTruthCst", "N_{cst} (truth)", nNumBins, rNumBins[0], rNumBins[1]);
138  m_hNumObject[OBJECT::RCST] = new TH1D("hNumRecoCst", "N_{cst} (reco.)", nNumBins, rNumBins[0], rNumBins[1]);
139  // sum of cst. energies
140  m_hSumCstEne[CST_TYPE::PART_CST] = new TH1D("hSumPartEne", "#SigmaE (cst. par.s)", nEneBins, rEneBins[0], rEneBins[1]);
141  m_hSumCstEne[CST_TYPE::FLOW_CST] = new TH1D("hSumFlowEne", "#SigmaE (cst. flow)", nEneBins, rEneBins[0], rEneBins[1]);
142  m_hSumCstEne[CST_TYPE::TRACK_CST] = new TH1D("hSumTrackEne", "#SigmaE (cst. track)", nEneBins, rEneBins[0], rEneBins[1]);
143  m_hSumCstEne[CST_TYPE::ECAL_CST] = new TH1D("hSumECalEne", "#SigmaE (cst. ecal)", nEneBins, rEneBins[0], rEneBins[1]);
144  m_hSumCstEne[CST_TYPE::HCAL_CST] = new TH1D("hSumHCalEne", "#SigmaE (cst. hcal)", nEneBins, rEneBins[0], rEneBins[1]);
145  // particle QA
146  m_hObjectQA[OBJECT::PART][INFO::PT] = new TH1D("hPartPt", "p_{T} (par.s)", nPtBins, rPtBins[0], rPtBins[1]);
147  m_hObjectQA[OBJECT::PART][INFO::ETA] = new TH1D("hPartEta", "#eta (par.s)", nEtaBins, rEtaBins[0], rEtaBins[1]);
148  m_hObjectQA[OBJECT::PART][INFO::PHI] = new TH1D("hPartPhi", "#phi (par.s)", nPhiBins, rPhiBins[0], rPhiBins[1]);
149  m_hObjectQA[OBJECT::PART][INFO::ENE] = new TH1D("hPartEne", "E (par.s)", nEneBins, rEneBins[0], rEneBins[1]);
150  m_hObjectQA[OBJECT::PART][INFO::QUAL] = new TH1D("hPartQuality", "Quality (par.s, N/A)", nQualBins, rQualBins[0], rQualBins[1]);
151  m_hObjectQA[OBJECT::PART][INFO::DCAXY] = new TH1D("hPartDcaXY", "DCA_{xy} (par.s, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
152  m_hObjectQA[OBJECT::PART][INFO::DCAZ] = new TH1D("hPartDcaZ", "DCA_{z} (par.s, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
153  m_hObjectQA[OBJECT::PART][INFO::DELTAPT] = new TH1D("hPartDeltaPt", "#deltap_{T}/p_{T} (par.s, N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
154  m_hObjectQA[OBJECT::PART][INFO::NTPC] = new TH1D("hPartNumTpc", "N_{clust}^{tpc} (par.s, N/A)", nNumBins, rNumBins[0], rNumBins[1]);
155  // track QA
156  m_hObjectQA[OBJECT::TRACK][INFO::PT] = new TH1D("hTrackPt", "p_{T} (tracks)", nPtBins, rPtBins[0], rPtBins[1]);
157  m_hObjectQA[OBJECT::TRACK][INFO::ETA] = new TH1D("hTrackEta", "#eta (tracks)", nEtaBins, rEtaBins[0], rEtaBins[1]);
158  m_hObjectQA[OBJECT::TRACK][INFO::PHI] = new TH1D("hTrackPhi", "#phi (tracks)", nPhiBins, rPhiBins[0], rPhiBins[1]);
159  m_hObjectQA[OBJECT::TRACK][INFO::ENE] = new TH1D("hTrackEne", "E (tracks)", nEneBins, rEneBins[0], rEneBins[1]);
160  m_hObjectQA[OBJECT::TRACK][INFO::QUAL] = new TH1D("hTrackQuality", "Quality (tracks)", nQualBins, rQualBins[0], rQualBins[1]);
161  m_hObjectQA[OBJECT::TRACK][INFO::DCAXY] = new TH1D("hTrackDcaXY", "DCA_{xy} (tracks)", nDcaBins, rDcaBins[0], rDcaBins[1]);
162  m_hObjectQA[OBJECT::TRACK][INFO::DCAZ] = new TH1D("hTrackDcaZ", "DCA_{z} (tracks)", nDcaBins, rDcaBins[0], rDcaBins[1]);
163  m_hObjectQA[OBJECT::TRACK][INFO::DELTAPT] = new TH1D("hTrackDeltaPt", "#deltap_{T}/p_{T} (tracks)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
164  m_hObjectQA[OBJECT::TRACK][INFO::NTPC] = new TH1D("hTrackNumTpc", "N_{clust}^{tpc} (tracks)", nNumBins, rNumBins[0], rNumBins[1]);
165  // particle flow QA
166  m_hObjectQA[OBJECT::FLOW][INFO::PT] = new TH1D("hFlowPt", "p_{T} (flow)", nPtBins, rPtBins[0], rPtBins[1]);
167  m_hObjectQA[OBJECT::FLOW][INFO::ETA] = new TH1D("hFlowEta", "#eta (flow)", nEtaBins, rEtaBins[0], rEtaBins[1]);
168  m_hObjectQA[OBJECT::FLOW][INFO::PHI] = new TH1D("hFlowPhi", "#phi (flow)", nPhiBins, rPhiBins[0], rPhiBins[1]);
169  m_hObjectQA[OBJECT::FLOW][INFO::ENE] = new TH1D("hFlowEne", "E (flow)", nEneBins, rEneBins[0], rEneBins[1]);
170  m_hObjectQA[OBJECT::FLOW][INFO::QUAL] = new TH1D("hFlowQuality", "Quality (flow, N/A)", nQualBins, rQualBins[0], rQualBins[1]);
171  m_hObjectQA[OBJECT::FLOW][INFO::DCAXY] = new TH1D("hFlowDcaXY", "DCA_{xy} (flow, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
172  m_hObjectQA[OBJECT::FLOW][INFO::DCAZ] = new TH1D("hFlowDcaZ", "DCA_{z} (flow, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
173  m_hObjectQA[OBJECT::FLOW][INFO::DELTAPT] = new TH1D("hFlowDeltaPt", "#deltap_{T}/p_{T} (flow, N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
174  m_hObjectQA[OBJECT::FLOW][INFO::NTPC] = new TH1D("hFlowNumTpc", "N_{clust}^{tpc} (flow, N/A)", nNumBins, rNumBins[0], rNumBins[1]);
175  // emcal cluster QA
176  m_hObjectQA[OBJECT::ECLUST][INFO::PT] = new TH1D("hEClustPt", "p_{T} (EMCal clust.)", nPtBins, rPtBins[0], rPtBins[1]);
177  m_hObjectQA[OBJECT::ECLUST][INFO::ETA] = new TH1D("hEClustEta", "#eta (EMCal clust.)", nEtaBins, rEtaBins[0], rEtaBins[1]);
178  m_hObjectQA[OBJECT::ECLUST][INFO::PHI] = new TH1D("hEClustPhi", "#phi (EMCal clust.)", nPhiBins, rPhiBins[0], rPhiBins[1]);
179  m_hObjectQA[OBJECT::ECLUST][INFO::ENE] = new TH1D("hEClustEne", "E (EMCal clust.)", nEneBins, rEneBins[0], rEneBins[1]);
180  m_hObjectQA[OBJECT::ECLUST][INFO::QUAL] = new TH1D("hEClustQuality", "Quality (EMCal, N/A)", nQualBins, rQualBins[0], rQualBins[1]);
181  m_hObjectQA[OBJECT::ECLUST][INFO::DCAXY] = new TH1D("hEClustDcaXY", "DCA_{xy} (EMCal, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
182  m_hObjectQA[OBJECT::ECLUST][INFO::DCAZ] = new TH1D("hEClustDcaZ", "DCA_{z} (EMCal, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
183  m_hObjectQA[OBJECT::ECLUST][INFO::DELTAPT] = new TH1D("hEClustDeltaPt", "#deltap_{T}/p_{T} (EMCal, N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
184  m_hObjectQA[OBJECT::ECLUST][INFO::NTPC] = new TH1D("hEClustNumTpc", "N_{clust}^{tpc} (EMCal, N/A)", nNumBins, rNumBins[0], rNumBins[1]);
185  // hcal cluster QA
186  m_hObjectQA[OBJECT::HCLUST][INFO::PT] = new TH1D("hHClustPt", "p_{T} (HCal clust.)", nPtBins, rPtBins[0], rPtBins[1]);
187  m_hObjectQA[OBJECT::HCLUST][INFO::ETA] = new TH1D("hHClustEta", "#eta (HCal clust.)", nEtaBins, rEtaBins[0], rEtaBins[1]);
188  m_hObjectQA[OBJECT::HCLUST][INFO::PHI] = new TH1D("hHClustPhi", "#phi (HCal clust.)", nPhiBins, rPhiBins[0], rPhiBins[1]);
189  m_hObjectQA[OBJECT::HCLUST][INFO::ENE] = new TH1D("hHClustEne", "E (HCal clust.)", nEneBins, rEneBins[0], rEneBins[1]);
190  m_hObjectQA[OBJECT::HCLUST][INFO::QUAL] = new TH1D("hHClustQuality", "Quality (HCal, N/A)", nQualBins, rQualBins[0], rQualBins[1]);
191  m_hObjectQA[OBJECT::HCLUST][INFO::DCAXY] = new TH1D("hHClustDcaXY", "DCA_{xy} (HCal, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
192  m_hObjectQA[OBJECT::HCLUST][INFO::DCAZ] = new TH1D("hHClustDcaZ", "DCA_{z} (HCal, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
193  m_hObjectQA[OBJECT::HCLUST][INFO::DELTAPT] = new TH1D("hHClustDeltaPt", "#deltap_{T}/p_{T} (HCal, N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
194  m_hObjectQA[OBJECT::HCLUST][INFO::NTPC] = new TH1D("hHClustNumTpc", "N_{clust}^{tpc} (HCal, N/A)", nNumBins, rNumBins[0], rNumBins[1]);
195  // truth jet QA
196  m_hObjectQA[OBJECT::TJET][INFO::PT] = new TH1D("hTruthJetPt", "p_{T} (truth jet)", nPtBins, rPtBins[0], rPtBins[1]);
197  m_hObjectQA[OBJECT::TJET][INFO::ETA] = new TH1D("hTruthJetEta", "#eta (truth jet)", nEtaBins, rEtaBins[0], rEtaBins[1]);
198  m_hObjectQA[OBJECT::TJET][INFO::PHI] = new TH1D("hTruthJetPhi", "#phi (truth jet)", nPhiBins, rPhiBins[0], rPhiBins[1]);
199  m_hObjectQA[OBJECT::TJET][INFO::ENE] = new TH1D("hTruthJetEne", "E (truth jet)", nEneBins, rEneBins[0], rEneBins[1]);
200  m_hObjectQA[OBJECT::TJET][INFO::QUAL] = new TH1D("hTruthJetQuality", "Quality (truth jet, N/A)", nQualBins, rQualBins[0], rQualBins[1]);
201  m_hObjectQA[OBJECT::TJET][INFO::DCAXY] = new TH1D("hTruthJetDcaXY", "DCA_{xy} (truth jet, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
202  m_hObjectQA[OBJECT::TJET][INFO::DCAZ] = new TH1D("hTruthJetDcaZ", "DCA_{z} (truth jet, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
203  m_hObjectQA[OBJECT::TJET][INFO::DELTAPT] = new TH1D("hTruthJetDeltaPt", "#deltap_{T}/p_{T} (truth jet, N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
204  m_hObjectQA[OBJECT::TJET][INFO::NTPC] = new TH1D("hTruthJetNumTpc", "N_{clust}^{tpc} (truth jet, N/A)", nNumBins, rNumBins[0], rNumBins[1]);
205  m_hJetArea[0] = new TH1D("hTruthJetArea", "Area (truth jet)", nAreaBins, rAreaBins[0], rAreaBins[1]);
206  m_hJetNumCst[0] = new TH1D("hTruthJetNCst", "N_{cst} (truth jet)", nNumBins, rNumBins[0], rNumBins[1]);
207  // reco jet QA
208  m_hObjectQA[OBJECT::RJET][INFO::PT] = new TH1D("hRecoJetPt", "p_{T} (reco. jet)", nPtBins, rPtBins[0], rPtBins[1]);
209  m_hObjectQA[OBJECT::RJET][INFO::ETA] = new TH1D("hRecoJetEta", "#eta (reco. jet)", nEtaBins, rEtaBins[0], rEtaBins[1]);
210  m_hObjectQA[OBJECT::RJET][INFO::PHI] = new TH1D("hRecoJetPhi", "#phi (reco. jet)", nPhiBins, rPhiBins[0], rPhiBins[1]);
211  m_hObjectQA[OBJECT::RJET][INFO::ENE] = new TH1D("hRecoJetEne", "E (reco. jet)", nEneBins, rEneBins[0], rEneBins[1]);
212  m_hObjectQA[OBJECT::RJET][INFO::QUAL] = new TH1D("hRecoJetQuality", "Quality (reco. jet, N/A)", nQualBins, rQualBins[0], rQualBins[1]);
213  m_hObjectQA[OBJECT::RJET][INFO::DCAXY] = new TH1D("hRecoJetDcaXY", "DCA_{xy} (reco. jet, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
214  m_hObjectQA[OBJECT::RJET][INFO::DCAZ] = new TH1D("hRecoJetDcaZ", "DCA_{z} (reco. jet, N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
215  m_hObjectQA[OBJECT::RJET][INFO::DELTAPT] = new TH1D("hRecoJetDeltaPt", "#deltap_{T}/p_{T} (reco. jet, N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
216  m_hObjectQA[OBJECT::RJET][INFO::NTPC] = new TH1D("hRecoJetNumTpc", "N_{clust}^{tpc} (reco. jet, N/A)", nNumBins, rNumBins[0], rNumBins[1]);
217  m_hJetArea[1] = new TH1D("hRecoJetArea", "Area (reco. jet)", nAreaBins, rAreaBins[0], rAreaBins[1]);
218  m_hJetNumCst[1] = new TH1D("hRecoJetNCst", "N_{cst} (reco. jet)", nNumBins, rNumBins[0], rNumBins[1]);
219  // truth cst. QA
220  m_hObjectQA[OBJECT::TCST][INFO::PT] = new TH1D("hTruthCstPt", "p_{T} (truth cst.)", nPtBins, rPtBins[0], rPtBins[1]);
221  m_hObjectQA[OBJECT::TCST][INFO::ETA] = new TH1D("hTruthCstEta", "#eta (truth cst.)", nEtaBins, rEtaBins[0], rEtaBins[1]);
222  m_hObjectQA[OBJECT::TCST][INFO::PHI] = new TH1D("hTruthCstPhi", "#phi (truth cst.)", nPhiBins, rPhiBins[0], rPhiBins[1]);
223  m_hObjectQA[OBJECT::TCST][INFO::ENE] = new TH1D("hTruthCstEne", "E (truth cst.)", nEneBins, rEneBins[0], rEneBins[1]);
224  m_hObjectQA[OBJECT::TCST][INFO::QUAL] = new TH1D("hTruthCstQuality", "Quality (truth cst., N/A)", nQualBins, rQualBins[0], rQualBins[1]);
225  m_hObjectQA[OBJECT::TCST][INFO::DCAXY] = new TH1D("hTruthCstDcaXY", "DCA_{xy} (truth cst., N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
226  m_hObjectQA[OBJECT::TCST][INFO::DCAZ] = new TH1D("hTruthCstDcaZ", "DCA_{z} (truth cst., N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
227  m_hObjectQA[OBJECT::TCST][INFO::DELTAPT] = new TH1D("hTruthCstDeltaPt", "#deltap_{T}/p_{T} (truth cst., N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
228  m_hObjectQA[OBJECT::TCST][INFO::NTPC] = new TH1D("hTruthCstNumTpc", "N_{clust}^{tpc} (truth cst., N/A)", nNumBins, rNumBins[0], rNumBins[1]);
229  // reco cst. QA
230  m_hObjectQA[OBJECT::RCST][INFO::PT] = new TH1D("hRecoCstPt", "p_{T} (reco. cst.)", nPtBins, rPtBins[0], rPtBins[1]);
231  m_hObjectQA[OBJECT::RCST][INFO::ETA] = new TH1D("hRecoCstEta", "#eta (reco. cst.)", nEtaBins, rEtaBins[0], rEtaBins[1]);
232  m_hObjectQA[OBJECT::RCST][INFO::PHI] = new TH1D("hRecoCstPhi", "#phi (reco. cst.)", nPhiBins, rPhiBins[0], rPhiBins[1]);
233  m_hObjectQA[OBJECT::RCST][INFO::ENE] = new TH1D("hRecoCstEne", "E (reco. cst.)", nEneBins, rEneBins[0], rEneBins[1]);
234  m_hObjectQA[OBJECT::RCST][INFO::QUAL] = new TH1D("hRecoCstQuality", "Quality (reco. cst., N/A)", nQualBins, rQualBins[0], rQualBins[1]);
235  m_hObjectQA[OBJECT::RCST][INFO::DCAXY] = new TH1D("hRecoCstDcaXY", "DCA_{xy} (reco. cst., N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
236  m_hObjectQA[OBJECT::RCST][INFO::DCAZ] = new TH1D("hRecoCstDcaZ", "DCA_{z} (reco. cst., N/A)", nDcaBins, rDcaBins[0], rDcaBins[1]);
237  m_hObjectQA[OBJECT::RCST][INFO::DELTAPT] = new TH1D("hRecoCstDeltaPt", "#deltap_{T}/p_{T} (reco. cst., N/A)", nDeltaBins, rDeltaBins[0], rDeltaBins[1]);
238  m_hObjectQA[OBJECT::RCST][INFO::NTPC] = new TH1D("hRecoCstNumTpc", "N_{clust}^{tpc} (reco. cst., N/A)", nNumBins, rNumBins[0], rNumBins[1]);
239  // no. of cst.s
240  m_hNumCstAccept[CST_TYPE::PART_CST][0] = new TH1D("hNumPartCstTot", "N_{cst}^{par} total", nNumBins, rNumBins[0], rNumBins[1]);
241  m_hNumCstAccept[CST_TYPE::PART_CST][1] = new TH1D("hNumPartCstAcc", "N_{cst}^{par} accepted", nNumBins, rNumBins[0], rNumBins[1]);
242  m_hNumCstAccept[CST_TYPE::TRACK_CST][0] = new TH1D("hNumTrkCstTot", "N_{cst}^{trk} total", nNumBins, rNumBins[0], rNumBins[1]);
243  m_hNumCstAccept[CST_TYPE::TRACK_CST][1] = new TH1D("hNumTrkCstAcc", "N_{cst}^{trk} accepted", nNumBins, rNumBins[0], rNumBins[1]);
244  m_hNumCstAccept[CST_TYPE::FLOW_CST][0] = new TH1D("hNumFlowCstTot", "N_{cst}^{flow} total", nNumBins, rNumBins[0], rNumBins[1]);
245  m_hNumCstAccept[CST_TYPE::FLOW_CST][1] = new TH1D("hNumFlowCstAcc", "N_{cst}^{flow} accepted", nNumBins, rNumBins[0], rNumBins[1]);
246  m_hNumCstAccept[CST_TYPE::ECAL_CST][0] = new TH1D("hNumECalCstTot", "N_{cst}^{clust} total", nNumBins, rNumBins[0], rNumBins[1]);
247  m_hNumCstAccept[CST_TYPE::ECAL_CST][1] = new TH1D("hNumECalCstAcc", "N_{cst}^{clust} accepted", nNumBins, rNumBins[0], rNumBins[1]);
248  m_hNumCstAccept[CST_TYPE::HCAL_CST][0] = new TH1D("hNumHCalCstTot", "N_{cst}^{clust} total", nNumBins, rNumBins[0], rNumBins[1]);
249  m_hNumCstAccept[CST_TYPE::HCAL_CST][1] = new TH1D("hNumHCalCstAcc", "N_{cst}^{clust} accepted", nNumBins, rNumBins[0], rNumBins[1]);
250 
251  // set errors
252  for (size_t iObj = OBJECT::TRACK; iObj < CONST::NObjType; iObj++) {
253  m_hNumObject[iObj] -> Sumw2();
254  for (size_t iInfo = INFO::PT; iInfo < CONST::NInfoQA; iInfo++) {
255  m_hObjectQA[iObj][iInfo] -> Sumw2();
256  } // end info loop
257  } // end object loop
258  for (size_t iCst = CST_TYPE::TRACK_CST; iCst < CONST::NCstType; iCst++) {
259  m_hNumCstAccept[iCst][0] -> Sumw2();
260  m_hNumCstAccept[iCst][1] -> Sumw2();
261  m_hSumCstEne[iCst] -> Sumw2();
262  } // end cst loop
263  m_hJetArea[0] -> Sumw2();
264  m_hJetNumCst[0] -> Sumw2();
265  m_hJetArea[1] -> Sumw2();
266  m_hJetNumCst[1] -> Sumw2();
267  return;
268 
269  } // end 'InitHists()'
270 
271 
272 
273  void SCorrelatorJetTree::InitFuncs() {
274 
275  // print debug statement
276  if (m_doDebug) {
277  cout << "SCorrelatorJetTree::InitFuncs() Initializing functions for internal calculations..." << endl;
278  }
279 
280  // pt range of functions
281  const pair<float, float> ptRange = {0., 100.};
282 
283  // initialize functions
284  m_fSigDcaXY = new TF1("fSigmaDcaXY", "[0]+[1]/x+[2]/(x*x)", ptRange.first, ptRange.second);
285  m_fSigDcaZ = new TF1("fSigmaDcaZ", "[0]+[1]/x+[2]/(x*x)", ptRange.first, ptRange.second);
286  for (uint8_t iParam = 0; iParam < CONST::NParam; iParam++) {
287  m_fSigDcaXY -> SetParameter(iParam, m_parSigDcaXY[iParam]);
288  m_fSigDcaZ -> SetParameter(iParam, m_parSigDcaZ[iParam]);
289  }
290  return;
291 
292  } // end 'InitFuncs()'
293 
294 
295 
296  void SCorrelatorJetTree::InitTuples() {
297 
298  // print debug statement
299  if (m_doDebug) {
300  cout << "SCorrelatorJetTree::InitTuples() Initializing output tuples..." << endl;
301  }
302 
303  // track QA leaves
304  const vector<string> vecTrkQALeaves = {
305  "pt",
306  "eta",
307  "phi",
308  "energy",
309  "dcaxy",
310  "dcaz",
311  "deltapt",
312  "quality",
313  "nmvtxlayer",
314  "ninttlayer",
315  "ntpclayer",
316  "vtxx",
317  "vtxy",
318  "vtxz"
319  };
320 
321  // flatten leaf list
322  string argTrkQALeaves("");
323  for (size_t iLeaf = 0; iLeaf < vecTrkQALeaves.size(); iLeaf++) {
324  argTrkQALeaves.append(vecTrkQALeaves[iLeaf]);
325  if ((iLeaf + 1) != vecTrkQALeaves.size()) {
326  argTrkQALeaves.append(":");
327  }
328  }
329  m_ntTrkQA = new TNtuple("ntTrkQA", "Track QA", argTrkQALeaves.data());
330 
331  // leaves for weird track check
332  const vector<string> vecWeirdTrkLeaves = {
333  "trkid_a",
334  "trkid_b",
335  "pt_a",
336  "pt_b",
337  "eta_a",
338  "eta_b",
339  "phi_a",
340  "phi_b",
341  "ene_a",
342  "ene_b",
343  "dcaxy_a",
344  "dcaxy_b",
345  "dcaz_a",
346  "dcaz_b",
347  "vtxx_a",
348  "vtxx_b",
349  "vtxy_a",
350  "vtxy_b",
351  "vtxz_a",
352  "vtxz_b",
353  "quality_a",
354  "quality_b",
355  "deltapt_a",
356  "deltapt_b",
357  "nmvtxlayers_a",
358  "nmvtxlayers_b",
359  "ninttlayers_a",
360  "ninttlayers_b",
361  "ntpclayers_a",
362  "ntpclayers_b",
363  "nmvtxclusts_a",
364  "nmvtxclusts_b",
365  "ninttclusts_a",
366  "ninttclusts_b",
367  "ntpcclusts_a",
368  "ntpcclusts_b",
369  "nclustkey_a",
370  "nclustkey_b",
371  "nsameclustkey",
372  "deltartrack"
373  };
374 
375  string argWeirdTrkLeaves("");
376  for (size_t iLeaf = 0; iLeaf < vecWeirdTrkLeaves.size(); iLeaf++) {
377  argWeirdTrkLeaves.append(vecWeirdTrkLeaves[iLeaf]);
378  if ((iLeaf + 1) != vecWeirdTrkLeaves.size()) {
379  argWeirdTrkLeaves.append(":");
380  }
381  }
382  m_ntWeirdTracks = new TNtuple("ntWeirdTracks", "Weird Tracks", argWeirdTrkLeaves.data());
383  return;
384 
385  } // end 'InitTuples()'
386 
387 
388 
389  void SCorrelatorJetTree::InitTrees() {
390 
391  // print debug statement
392  if (m_doDebug) {
393  cout << "SCorrelatorJetTree::InitTrees() Initializing output trees..." << endl;
394  }
395 
396  // initialize true (inclusive) jet tree
397  m_trueTree = new TTree("TruthJetTree", "A tree of truth jets");
398  m_trueTree -> Branch("EvtNumJets", &m_trueNumJets, "EvtNumJets/I");
399  m_trueTree -> Branch("EvtNumChrgPars", &m_trueNumChrgPars, "EvtNumChrgPars/I");
400  m_trueTree -> Branch("Parton3_ID", &m_truePartonID[0], "Parton3_ID/I");
401  m_trueTree -> Branch("Parton4_ID", &m_truePartonID[1], "Parton4_ID/I");
402  m_trueTree -> Branch("Parton3_MomX", &m_truePartonMomX[0], "Parton3_MomX/D");
403  m_trueTree -> Branch("Parton3_MomY", &m_truePartonMomY[0], "Parton3_MomY/D");
404  m_trueTree -> Branch("Parton3_MomZ", &m_truePartonMomZ[0], "Parton3_MomZ/D");
405  m_trueTree -> Branch("Parton4_MomX", &m_truePartonMomX[1], "Parton4_MomX/D");
406  m_trueTree -> Branch("Parton4_MomY", &m_truePartonMomY[1], "Parton4_MomY/D");
407  m_trueTree -> Branch("Parton4_MomZ", &m_truePartonMomZ[1], "Parton4_MomZ/D");
408  m_trueTree -> Branch("EvtVtxX", &m_trueVtxX, "EvtVtxX/D");
409  m_trueTree -> Branch("EvtVtxY", &m_trueVtxY, "EvtVtxY/D");
410  m_trueTree -> Branch("EvtVtxZ", &m_trueVtxZ, "EvtVtxZ/D");
411  m_trueTree -> Branch("EvtSumParEne", &m_trueSumPar, "EvtSumParEne/D");
412  m_trueTree -> Branch("JetNumCst", &m_trueJetNCst);
413  m_trueTree -> Branch("JetID", &m_trueJetID);
414  m_trueTree -> Branch("JetEnergy", &m_trueJetE);
415  m_trueTree -> Branch("JetPt", &m_trueJetPt);
416  m_trueTree -> Branch("JetEta", &m_trueJetEta);
417  m_trueTree -> Branch("JetPhi", &m_trueJetPhi);
418  m_trueTree -> Branch("JetArea", &m_trueJetArea);
419  m_trueTree -> Branch("CstID", &m_trueCstID);
420  m_trueTree -> Branch("CstEmbedID", &m_trueCstEmbedID);
421  m_trueTree -> Branch("CstZ", &m_trueCstZ);
422  m_trueTree -> Branch("CstDr", &m_trueCstDr);
423  m_trueTree -> Branch("CstEnergy", &m_trueCstE);
424  m_trueTree -> Branch("CstPt", &m_trueCstPt);
425  m_trueTree -> Branch("CstEta", &m_trueCstEta);
426  m_trueTree -> Branch("CstPhi", &m_trueCstPhi);
427 
428  // initialize reco jet tree
429  m_recoTree = new TTree("RecoJetTree", "A tree of reconstructed jets");
430  m_recoTree -> Branch("EvtNumJets", &m_recoNumJets, "EvtNumJets/I");
431  m_recoTree -> Branch("EvtNumTrks", &m_recoNumTrks, "EvtNumTrks/I");
432  m_recoTree -> Branch("EvtVtxX", &m_recoVtxX, "EvtVtxX/D");
433  m_recoTree -> Branch("EvtVtxY", &m_recoVtxY, "EvtVtxY/D");
434  m_recoTree -> Branch("EvtVtxZ", &m_recoVtxZ, "EvtVtxZ/D");
435  m_recoTree -> Branch("EvtSumECalEne", &m_recoSumECal, "EvtSumECalEne/D");
436  m_recoTree -> Branch("EvtSumHCalEne", &m_recoSumHCal, "EvtSumHCalEne/D");
437  m_recoTree -> Branch("JetNumCst", &m_recoJetNCst);
438  m_recoTree -> Branch("JetID", &m_recoJetID);
439  m_recoTree -> Branch("JetEnergy", &m_recoJetE);
440  m_recoTree -> Branch("JetPt", &m_recoJetPt);
441  m_recoTree -> Branch("JetEta", &m_recoJetEta);
442  m_recoTree -> Branch("JetPhi", &m_recoJetPhi);
443  m_recoTree -> Branch("JetArea", &m_recoJetArea);
444  m_recoTree -> Branch("CstMatchID", &m_recoCstMatchID);
445  m_recoTree -> Branch("CstZ", &m_recoCstZ);
446  m_recoTree -> Branch("CstDr", &m_recoCstDr);
447  m_recoTree -> Branch("CstEnergy", &m_recoCstE);
448  m_recoTree -> Branch("CstPt", &m_recoCstPt);
449  m_recoTree -> Branch("CstEta", &m_recoCstEta);
450  m_recoTree -> Branch("CstPhi", &m_recoCstPhi);
451  return;
452 
453  } // end 'InitTrees()'
454 
455 
456 
457  void SCorrelatorJetTree::InitEvals(PHCompositeNode* topNode) {
458 
459  // print debug statement
460  if (m_doDebug) {
461  cout << "SCorrelatorJetTree::InitEvals(PHCompositeNode*) Initializing evaluators..." << endl;
462  }
463 
464  m_evalStack = new SvtxEvalStack(topNode);
465  if (!m_evalStack) {
466  cerr << "SCorrelatorJetTree::InitEvals(PHCompositeNode*) PANIC: couldn't grab SvtxEvalStack! Aborting!" << endl;
467  assert(m_evalStack);
468  } else {
469  m_evalStack -> next_event(topNode);
470  }
471 
472  m_trackEval = m_evalStack -> get_track_eval();
473  if (!m_trackEval) {
474  cerr << "SCorrelatorJetTree::InitEvals(PHCompositeNode*) PANIC: couldn't grab track evaluator! Aborting!" << endl;
475  assert(m_trackEval);
476  }
477  return;
478 
479  } // end 'InitEvals(PHCompositeNode*)'
480 
481 
482 
483  void SCorrelatorJetTree::FillTrueTree() {
484 
485  // print debug statement
486  if (m_doDebug) {
487  cout << "SCorrelatorJetTree::FillTrueTree() Filling truth jet tree..." << endl;
488  }
489 
490  // prepare vectors for filling
491  m_trueJetNCst.clear();
492  m_trueJetID.clear();
493  m_trueJetE.clear();
494  m_trueJetPt.clear();
495  m_trueJetEta.clear();
496  m_trueJetPhi.clear();
497  m_trueJetArea.clear();
498  m_trueCstID.clear();
499  m_trueCstEmbedID.clear();
500  m_trueCstZ.clear();
501  m_trueCstDr.clear();
502  m_trueCstE.clear();
503  m_trueCstPt.clear();
504  m_trueCstEta.clear();
505  m_trueCstPhi.clear();
506 
507  // declare vectors to storing constituents
508  vector<int> vecTruCstID;
509  vector<int> vecTruCstEmbedID;
510  vector<double> vecTruCstZ;
511  vector<double> vecTruCstDr;
512  vector<double> vecTruCstE;
513  vector<double> vecTruCstPt;
514  vector<double> vecTruCstEta;
515  vector<double> vecTruCstPhi;
516  vecTruCstID.clear();
517  vecTruCstZ.clear();
518  vecTruCstDr.clear();
519  vecTruCstE.clear();
520  vecTruCstPt.clear();
521  vecTruCstEta.clear();
522  vecTruCstPhi.clear();
523 
524  // fill jets/constituent variables
525  unsigned int nTruJet(0);
526  unsigned int nTruCst(0);
527  for (unsigned int iTruJet = 0; iTruJet < m_trueJets.size(); ++iTruJet) {
528 
529  // get jet info
530  const unsigned int jetNCst = m_trueJets[iTruJet].constituents().size();
531  const unsigned int jetTruID = iTruJet;
532  const double jetPhi = m_trueJets[iTruJet].phi_std();
533  const double jetEta = m_trueJets[iTruJet].pseudorapidity();
534  const double jetArea = 0.; // FIXME: jet area needs to be defined
535  const double jetE = m_trueJets[iTruJet].E();
536  const double jetPt = m_trueJets[iTruJet].perp();
537  const double jetPx = m_trueJets[iTruJet].px();
538  const double jetPy = m_trueJets[iTruJet].py();
539  const double jetPz = m_trueJets[iTruJet].pz();
540  const double jetP = sqrt((jetPx * jetPx) + (jetPy * jetPy) + (jetPz * jetPz));
541 
542  // clear constituent vectors
543  vecTruCstID.clear();
544  vecTruCstEmbedID.clear();
545  vecTruCstZ.clear();
546  vecTruCstDr.clear();
547  vecTruCstE.clear();
548  vecTruCstPt.clear();
549  vecTruCstEta.clear();
550  vecTruCstPhi.clear();
551 
552  // loop over constituents
553  vector<fastjet::PseudoJet> trueCsts = m_trueJets[iTruJet].constituents();
554  for (unsigned int iTruCst = 0; iTruCst < trueCsts.size(); ++iTruCst) {
555 
556  // get constituent info
557  const double cstPhi = trueCsts[iTruCst].phi_std();
558  const double cstEta = trueCsts[iTruCst].pseudorapidity();
559  const double cstE = trueCsts[iTruCst].E();
560  const double cstPt = trueCsts[iTruCst].perp();
561  const double cstPx = trueCsts[iTruCst].px();
562  const double cstPy = trueCsts[iTruCst].py();
563  const double cstPz = trueCsts[iTruCst].pz();
564  const double cstP = ((cstPx * cstPx) + (cstPy * cstPy) + (cstPz * cstPz));
565  const double cstZ = cstP / jetP;
566  const double cstDf = cstPhi - jetPhi;
567  const double cstDh = cstEta - jetEta;
568  const double cstDr = sqrt((cstDf * cstDf) + (cstDh * cstDh));
569 
570  // get barcode and embedding ID
571  const int cstID = trueCsts[iTruCst].user_index();
572  const int embedID = m_mapCstToEmbedID[cstID];
573 
574  // add csts to vectors
575  vecTruCstID.push_back(abs(cstID));
576  vecTruCstEmbedID.push_back(embedID);
577  vecTruCstZ.push_back(cstZ);
578  vecTruCstDr.push_back(cstDr);
579  vecTruCstE.push_back(cstE);
580  vecTruCstPt.push_back(cstPt);
581  vecTruCstEta.push_back(cstEta);
582  vecTruCstPhi.push_back(cstPhi);
583 
584  // fill QA histograms and increment counters
585  m_hObjectQA[OBJECT::TCST][INFO::PT] -> Fill(cstPt);
586  m_hObjectQA[OBJECT::TCST][INFO::ETA] -> Fill(cstEta);
587  m_hObjectQA[OBJECT::TCST][INFO::PHI] -> Fill(cstPhi);
588  m_hObjectQA[OBJECT::TCST][INFO::ENE] -> Fill(cstE);
589  ++nTruCst;
590  } // end constituent loop
591 
592  // store jet/cst output
593  m_trueJetNCst.push_back(jetNCst);
594  m_trueJetID.push_back(jetTruID);
595  m_trueJetE.push_back(jetE);
596  m_trueJetPt.push_back(jetPt);
597  m_trueJetEta.push_back(jetEta);
598  m_trueJetPhi.push_back(jetPhi);
599  m_trueJetArea.push_back(jetArea);
600  m_trueCstID.push_back(vecTruCstID);
601  m_trueCstEmbedID.push_back(vecTruCstEmbedID);
602  m_trueCstZ.push_back(vecTruCstZ);
603  m_trueCstDr.push_back(vecTruCstDr);
604  m_trueCstE.push_back(vecTruCstE);
605  m_trueCstPt.push_back(vecTruCstPt);
606  m_trueCstEta.push_back(vecTruCstEta);
607  m_trueCstPhi.push_back(vecTruCstPhi);
608 
609  // fill QA histograms and increment counters
610  m_hJetArea[0] -> Fill(jetArea);
611  m_hJetNumCst[0] -> Fill(jetNCst);
612  m_hObjectQA[OBJECT::TJET][INFO::PT] -> Fill(jetPt);
613  m_hObjectQA[OBJECT::TJET][INFO::ETA] -> Fill(jetEta);
614  m_hObjectQA[OBJECT::TJET][INFO::PHI] -> Fill(jetPhi);
615  m_hObjectQA[OBJECT::TJET][INFO::ENE] -> Fill(jetE);
616  ++nTruJet;
617  } // end jet loop
618 
619  // fill QA histograms
620  m_hNumObject[OBJECT::TJET] -> Fill(nTruJet);
621  m_hNumObject[OBJECT::TCST] -> Fill(nTruCst);
622 
623  // store evt info
624  m_trueNumJets = nTruJet;
625  m_truePartonID[0] = m_partonID[0];
626  m_truePartonID[1] = m_partonID[1];
627  m_truePartonMomX[0] = m_partonMom[0].x();
628  m_truePartonMomX[1] = m_partonMom[1].x();
629  m_truePartonMomY[0] = m_partonMom[0].y();
630  m_truePartonMomY[1] = m_partonMom[1].y();
631  m_truePartonMomZ[0] = m_partonMom[0].z();
632  m_truePartonMomZ[1] = m_partonMom[1].z();
633  m_trueVtxX = m_trueVtx.x();
634  m_trueVtxY = m_trueVtx.y();
635  m_trueVtxZ = m_trueVtx.z();
636 
637  // fill output tree
638  m_trueTree -> Fill();
639  return;
640 
641  } // end 'FillTrueTree()'
642 
643 
644 
645  void SCorrelatorJetTree::FillRecoTree() {
646 
647  // print debug statement
648  if (m_doDebug) {
649  cout << "SCorrelatorJetTree::FillRecoTree() Filling reco jet tree..." << endl;
650  }
651 
652  // prepare vectors for filling
653  m_recoJetNCst.clear();
654  m_recoJetID.clear();
655  m_recoJetE.clear();
656  m_recoJetPt.clear();
657  m_recoJetEta.clear();
658  m_recoJetPhi.clear();
659  m_recoJetArea.clear();
660  m_recoCstMatchID.clear();
661  m_recoCstZ.clear();
662  m_recoCstDr.clear();
663  m_recoCstE.clear();
664  m_recoCstPt.clear();
665  m_recoCstEta.clear();
666  m_recoCstPhi.clear();
667 
668  // declare vectors for storing constituents
669  vector<int> vecRecCstMatchID;
670  vector<double> vecRecCstZ;
671  vector<double> vecRecCstDr;
672  vector<double> vecRecCstE;
673  vector<double> vecRecCstPt;
674  vector<double> vecRecCstEta;
675  vector<double> vecRecCstPhi;
676  vecRecCstMatchID.clear();
677  vecRecCstZ.clear();
678  vecRecCstDr.clear();
679  vecRecCstE.clear();
680  vecRecCstPt.clear();
681  vecRecCstEta.clear();
682  vecRecCstPhi.clear();
683 
684  // fill jet/constituent variables
685  unsigned long nRecJet(0);
686  unsigned long nRecCst(0);
687  for (unsigned int iJet = 0; iJet < m_recoJets.size(); ++iJet) {
688 
689  // get jet info
690  const unsigned int jetNCst = m_recoJets[iJet].constituents().size();
691  const unsigned int jetRecID = iJet;
692  const double jetPhi = m_recoJets[iJet].phi_std();
693  const double jetEta = m_recoJets[iJet].pseudorapidity();
694  const double jetArea = 0.; // FIXME: jet area needs to be defined
695  const double jetE = m_recoJets[iJet].E();
696  const double jetPt = m_recoJets[iJet].perp();
697  const double jetPx = m_recoJets[iJet].px();
698  const double jetPy = m_recoJets[iJet].py();
699  const double jetPz = m_recoJets[iJet].pz();
700  const double jetP = sqrt((jetPx * jetPx) + (jetPy * jetPy) + (jetPz * jetPz));
701 
702  // clear constituent vectors
703  vecRecCstMatchID.clear();
704  vecRecCstZ.clear();
705  vecRecCstDr.clear();
706  vecRecCstE.clear();
707  vecRecCstPt.clear();
708  vecRecCstEta.clear();
709  vecRecCstPhi.clear();
710 
711  // loop over constituents
712  vector<fastjet::PseudoJet> recoCsts = m_recoJets[iJet].constituents();
713  for (unsigned int iCst = 0; iCst < recoCsts.size(); ++iCst) {
714 
715  // get constituent info
716  const double cstMatchID = recoCsts[iCst].user_index();
717  const double cstPhi = recoCsts[iCst].phi_std();
718  const double cstEta = recoCsts[iCst].pseudorapidity();
719  const double cstE = recoCsts[iCst].E();
720  const double cstPt = recoCsts[iCst].perp();
721  const double cstPx = recoCsts[iCst].px();
722  const double cstPy = recoCsts[iCst].py();
723  const double cstPz = recoCsts[iCst].pz();
724  const double cstP = ((cstPx * cstPx) + (cstPy * cstPy) + (cstPz * cstPz));
725  const double cstZ = cstP / jetP;
726  const double cstDf = cstPhi - jetPhi;
727  const double cstDh = cstEta - jetEta;
728  const double cstDr = sqrt((cstDf * cstDf) + (cstDh * cstDh));
729 
730  // add csts to vectors
731  vecRecCstMatchID.push_back(cstMatchID);
732  vecRecCstZ.push_back(cstZ);
733  vecRecCstDr.push_back(cstDr);
734  vecRecCstE.push_back(cstE);
735  vecRecCstPt.push_back(cstPt);
736  vecRecCstEta.push_back(cstEta);
737  vecRecCstPhi.push_back(cstPhi);
738 
739  // fill QA histograms and increment counters
740  m_hObjectQA[OBJECT::RCST][INFO::PT] -> Fill(cstPt);
741  m_hObjectQA[OBJECT::RCST][INFO::ETA] -> Fill(cstEta);
742  m_hObjectQA[OBJECT::RCST][INFO::PHI] -> Fill(cstPhi);
743  m_hObjectQA[OBJECT::RCST][INFO::ENE] -> Fill(cstE);
744  ++nRecCst;
745  } // end constituent loop
746 
747  // store jet/cst output
748  m_recoJetNCst.push_back(jetNCst);
749  m_recoJetID.push_back(jetRecID);
750  m_recoJetE.push_back(jetE);
751  m_recoJetPt.push_back(jetPt);
752  m_recoJetEta.push_back(jetEta);
753  m_recoJetPhi.push_back(jetPhi);
754  m_recoJetArea.push_back(jetArea);
755  m_recoCstMatchID.push_back(vecRecCstMatchID);
756  m_recoCstZ.push_back(vecRecCstZ);
757  m_recoCstDr.push_back(vecRecCstDr);
758  m_recoCstE.push_back(vecRecCstE);
759  m_recoCstPt.push_back(vecRecCstPt);
760  m_recoCstEta.push_back(vecRecCstEta);
761  m_recoCstPhi.push_back(vecRecCstPhi);
762 
763  // fill QA histograms and increment counters
764  m_hJetArea[1] -> Fill(jetArea);
765  m_hJetNumCst[1] -> Fill(jetNCst);
766  m_hObjectQA[OBJECT::RJET][INFO::PT] -> Fill(jetPt);
767  m_hObjectQA[OBJECT::RJET][INFO::ETA] -> Fill(jetEta);
768  m_hObjectQA[OBJECT::RJET][INFO::PHI] -> Fill(jetPhi);
769  m_hObjectQA[OBJECT::RJET][INFO::ENE] -> Fill(jetE);
770  ++nRecJet;
771  } // end jet loop
772 
773  // fill QA histograms
774  m_hNumObject[OBJECT::RJET] -> Fill(nRecJet);
775  m_hNumObject[OBJECT::RCST] -> Fill(nRecCst);
776 
777  // store event info
778  m_recoNumJets = nRecJet;
779  m_recoVtxX = m_recoVtx.x();
780  m_recoVtxY = m_recoVtx.y();
781  m_recoVtxZ = m_recoVtx.z();
782 
783  // fill object tree
784  m_recoTree -> Fill();
785  return;
786 
787  } // end 'FillRecoTree()'
788 
789 
790 
791  void SCorrelatorJetTree::SaveOutput() {
792 
793  // print debug statement
794  if (m_doDebug) {
795  cout << "SCorrelatorJetTree::SaveOutput() Saving output trees and histograms..." << endl;
796  }
797 
798  // save QA histograms if need be
799  const string sQuality[CONST::NDirectory + 1] = {"Tracks", "CaloClusters", "ParticleFlow", "Particles", "TruthJets", "RecoJets", "QA"};
800  TDirectory* dQuality[CONST::NDirectory + 1];
801  if (m_doQualityPlots) {
802 
803  // create QA directories
804  dQuality[CONST::NDirectory] = (TDirectory*) m_outFile -> mkdir(sQuality[CONST::NDirectory].data());
805  for (size_t iDirect = 0; iDirect < CONST::NDirectory; iDirect++) {
806  dQuality[iDirect] = (TDirectory*) dQuality[CONST::NDirectory] -> mkdir(sQuality[iDirect].data());
807  }
808 
809  // save object-specific QA hists
810  for (size_t iObj = OBJECT::TRACK; iObj < CONST::NObjType; iObj++) {
811  switch (iObj) {
812  case OBJECT::TRACK:
813  dQuality[0] -> cd();
814  break;
815  case OBJECT::HCLUST:
816  dQuality[1] -> cd();
817  break;
818  case OBJECT::ECLUST:
819  dQuality[1] -> cd();
820  break;
821  case OBJECT::FLOW:
822  dQuality[2] -> cd();
823  break;
824  case OBJECT::PART:
825  dQuality[3] -> cd();
826  break;
827  case OBJECT::TJET:
828  dQuality[4] -> cd();
829  break;
830  case OBJECT::RJET:
831  dQuality[5] -> cd();
832  break;
833  case OBJECT::TCST:
834  dQuality[4] -> cd();
835  break;
836  case OBJECT::RCST:
837  dQuality[5] -> cd();
838  break;
839  }
840  m_hNumObject[iObj] -> Write();
841  for (size_t iInfo = INFO::PT; iInfo < CONST::NInfoQA; iInfo++) {
842  m_hObjectQA[iObj][iInfo] -> Write();
843  }
844  } // end object loop
845 
846  // save cst-specific histograms
847  for (size_t iCst = CST_TYPE::TRACK_CST; iCst < CONST::NCstType; iCst++) {
848  switch (iCst) {
849  case CST_TYPE::TRACK_CST:
850  dQuality[0] -> cd();
851  break;
852  case CST_TYPE::ECAL_CST:
853  dQuality[1] -> cd();
854  break;
855  case CST_TYPE::HCAL_CST:
856  dQuality[1] -> cd();
857  break;
858  case CST_TYPE::FLOW_CST:
859  dQuality[2] -> cd();
860  break;
861  case CST_TYPE::PART_CST:
862  dQuality[3] -> cd();
863  break;
864  }
865  m_hNumCstAccept[iCst][0] -> Write();
866  m_hNumCstAccept[iCst][1] -> Write();
867  m_hSumCstEne[iCst] -> Write();
868  } // end cst loop
869 
870  // save jet-specific histograms
871  dQuality[4] -> cd();
872  m_hJetArea[0] -> Write();
873  m_hJetNumCst[0] -> Write();
874  dQuality[5] -> cd();
875  m_hJetArea[1] -> Write();
876  m_hJetNumCst[1] -> Write();
877  }
878 
879  // save QA tuples
880  dQuality[0] -> cd();
881  m_ntTrkQA -> Write();
882  if (m_checkWeirdTrks) {
883  m_ntWeirdTracks -> Write();
884  }
885 
886  // save output trees
887  m_outFile -> cd();
888  m_recoTree -> Write();
889  if (m_isMC) {
890  m_trueTree -> Write();
891  }
892  return;
893 
894  } // end 'SaveOutput()'
895 
896 
897 
898  void SCorrelatorJetTree::ResetVariables() {
899 
900  // print debug statement
901  if (m_doDebug) {
902  cout << "SCorrelatorJetTree::ResetTreeVariables() Resetting tree variables..." << endl;
903  }
904 
905  // reset fastjet members
906  m_trueJetDef = NULL;
907  m_recoJetDef = NULL;
908  m_trueClust = NULL;
909  m_recoClust = NULL;
910 
911  // reset parton and other variables
912  m_partonID[0] = -9999;
913  m_partonID[1] = -9999;
914  m_partonMom[0] = CLHEP::Hep3Vector(-9999., -9999., -9999.);
915  m_partonMom[1] = CLHEP::Hep3Vector(-9999., -9999., -9999.);
916  m_vecEvtsToGrab.clear();
917  m_mapCstToEmbedID.clear();
918 
919  // reset truth (inclusive) tree variables
920  m_trueVtx = CLHEP::Hep3Vector(-9999., -9999., -9999.);
921  m_trueNumJets = 0;
922  m_trueNumChrgPars = -9999;
923  m_trueSumPar = -9999.;
924  m_truePartonID[0] = -9999;
925  m_truePartonID[1] = -9999;
926  m_truePartonMomX[0] = -9999.;
927  m_truePartonMomX[1] = -9999.;
928  m_truePartonMomY[0] = -9999.;
929  m_truePartonMomY[1] = -9999.;
930  m_truePartonMomZ[0] = -9999.;
931  m_truePartonMomZ[1] = -9999.;
932  m_trueJets.clear();
933  m_trueJetNCst.clear();
934  m_trueJetID.clear();
935  m_trueJetE.clear();
936  m_trueJetPt.clear();
937  m_trueJetEta.clear();
938  m_trueJetPhi.clear();
939  m_trueJetArea.clear();
940  m_trueCstID.clear();
941  m_trueCstEmbedID.clear();
942  m_trueCstZ.clear();
943  m_trueCstDr.clear();
944  m_trueCstE.clear();
945  m_trueCstPt.clear();
946  m_trueCstEta.clear();
947  m_trueCstPhi.clear();
948 
949  // reset reco tree variables
950  m_recoVtx = CLHEP::Hep3Vector(-9999., -9999., -9999.);
951  m_recoNumJets = 0;
952  m_recoNumTrks = -9999;
953  m_recoSumECal = -9999.;
954  m_recoSumHCal = -9999.;
955  m_recoJets.clear();
956  m_recoJetNCst.clear();
957  m_recoJetID.clear();
958  m_recoJetE.clear();
959  m_recoJetPt.clear();
960  m_recoJetEta.clear();
961  m_recoJetPhi.clear();
962  m_recoJetArea.clear();
963  m_recoCstZ.clear();
964  m_recoCstDr.clear();
965  m_recoCstE.clear();
966  m_recoCstPt.clear();
967  m_recoCstEta.clear();
968  m_recoCstPhi.clear();
969  return;
970 
971  } // end 'ResetTreeVariables()
972 
973 
974 
975  void SCorrelatorJetTree::DetermineEvtsToGrab(PHCompositeNode* topNode) {
976 
977  // print debug statement
978  if (m_doDebug) {
979  cout << "SCorrelatorJetTree::DetermineEvtsToGrab() Determining which subevents to grab..." << endl;
980  }
981 
982  // make sure vector is clear
983  m_vecEvtsToGrab.clear();
984 
985  // if not embedding, grab signal event
986  // otherwise add all subevents
987  if (!m_isEmbed) {
988  m_vecEvtsToGrab.push_back(1);
989  } else {
990  PHHepMCGenEventMap* mapMcEvts = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
991  for (PHHepMCGenEventMap::ConstIter itEvt = mapMcEvts -> begin(); itEvt != mapMcEvts -> end(); ++itEvt) {
992  m_vecEvtsToGrab.push_back(itEvt -> second -> get_embedding_id());
993  }
994  }
995  return;
996 
997  } // end 'DetermineEvtsToGrab()'
998 
999 
1000 
1001  int SCorrelatorJetTree::CreateJetNode(PHCompositeNode* topNode) {
1002 
1003  // print debug statement
1004  if (m_doDebug) {
1005  cout << "SCorrelatorJetTree::CreateJetNode(PHCompositeNode*) Creating jet node..." << endl;
1006  }
1007 
1008  // create iterator & DST node
1009  PHNodeIterator iter(topNode);
1010  PHCompositeNode* lowerNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
1011  if (!lowerNode) {
1012  lowerNode = new PHCompositeNode("DST");
1013  topNode -> addNode(lowerNode);
1014  cout << "DST node added" << endl;
1015  }
1016 
1017  // construct jet tree name
1018  string baseName;
1019  string recoNodeName;
1020  string trueNodeName;
1021  if (m_jetTreeName.empty()) {
1022  baseName = "JetTree";
1023  } else {
1024  baseName = m_jetTreeName;
1025  }
1026 
1027  // can't have forward slashes in DST or else you make a subdirectory on save!!!
1028  string undrscr = "_";
1029  string nothing = "";
1030 
1031  // define good strings to replace bad ones
1032  map<string, string> forbiddenStrings;
1033  forbiddenStrings["/"] = undrscr;
1034  forbiddenStrings["("] = undrscr;
1035  forbiddenStrings[")"] = nothing;
1036  forbiddenStrings["+"] = "plus";
1037  forbiddenStrings["-"] = "minus";
1038  forbiddenStrings["*"] = "star";
1039  for (auto const& [badString, goodString] : forbiddenStrings) {
1040  size_t pos;
1041  while ((pos = baseName.find(badString)) != string::npos) {
1042  baseName.replace(pos, 1, goodString);
1043  }
1044  }
1045 
1046  // construct jet node name
1047  recoNodeName = baseName + "_RecoJets";
1048  trueNodeName = baseName + "_TruthJets";
1049 
1050  // construct jet maps
1051  m_recoJetMap = new JetMapv1();
1052  if (m_isMC && m_saveDST) {
1053  m_trueJetMap = new JetMapv1();
1054  }
1055 
1056  // add jet node
1057  if (m_saveDST) {
1058  PHIODataNode<PHObject>* recoJetNode = new PHIODataNode<PHObject>(m_recoJetMap, recoNodeName.c_str(), "PHObject");
1059  lowerNode -> addNode(recoJetNode);
1060  cout << recoNodeName << " node added" << endl;
1061  }
1062 
1063  // save truth DSTs if needed
1064  if(m_isMC && m_saveDST) {
1065  PHIODataNode<PHObject> *trueJetNode = new PHIODataNode<PHObject>(m_trueJetMap, trueNodeName.c_str(), "PHObject");
1066  lowerNode -> addNode(trueJetNode);
1067  cout << trueNodeName << " node added" << endl;
1068  }
1069  return Fun4AllReturnCodes::EVENT_OK;
1070 
1071  } // end 'CreateJetNode(PHCompositeNode*)'
1072 
1073 
1074 
1075  int SCorrelatorJetTree::GetEmbedID(PHCompositeNode* topNode, const int iEvtToGrab) {
1076 
1077  // print debug statement
1078  if (m_doDebug) {
1079  cout << "SCorrelatorJetTree::GetEmbedID(PHCompositeNode*) Grabbing embedding ID from MC subevent #" << iEvtToGrab << "..." << endl;
1080  }
1081 
1082  // grab mc event map
1083  PHHepMCGenEventMap* mapMcEvts = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
1084  if (!mapMcEvts) {
1085  cerr << PHWHERE
1086  << "PANIC: HEPMC event map node is missing!"
1087  << endl;
1088  assert(mapMcEvts);
1089  }
1090 
1091  // grab mc event & return embedding id
1092  PHHepMCGenEvent* mcEvtStart = mapMcEvts -> get(iEvtToGrab);
1093  if (!mcEvtStart) {
1094  cerr << PHWHERE
1095  << "PANIC: Couldn't grab start of mc events!"
1096  << endl;
1097  assert(mcEvtStart);
1098  }
1099  return mcEvtStart -> get_embedding_id();
1100 
1101  } // end 'GetEmbedID(PHCompositeNode*, int)'
1102 
1103 
1104 
1105  SvtxTrackMap* SCorrelatorJetTree::GetTrackMap(PHCompositeNode* topNode) {
1106 
1107  // print debug statement
1108  if (m_doDebug) {
1109  cout << "SCorrelatorJetTree::GetTrackMap(PHCompositeNode*) Grabbing track map..." << endl;
1110  }
1111 
1112  // grab track map
1113  SvtxTrackMap* mapTrks = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
1114  if (!mapTrks) {
1115  cerr << PHWHERE
1116  << "PANIC: SvtxTrackMap node is missing!"
1117  << endl;
1118  assert(mapTrks);
1119  }
1120  return mapTrks;
1121 
1122  } // end 'GetTrackMap(PHCompositeNode*)'
1123 
1124 
1125 
1126  GlobalVertex* SCorrelatorJetTree::GetGlobalVertex(PHCompositeNode* topNode, const int iVtxToGrab) {
1127 
1128  // print debug statement
1129  if (m_doDebug) {
1130  cout << "SCorrelatorJetTree::GetGlobalVertex(PHCompositeNode*) Getting global vertex..." << endl;
1131  }
1132 
1133  // get vertex map
1134  GlobalVertexMap* mapVtx = GetVertexMap(topNode);
1135 
1136  // get specified vertex
1137  GlobalVertex* vtx = NULL;
1138  if (iVtxToGrab < 0) {
1139  vtx = mapVtx -> begin() -> second;
1140  } else {
1141  vtx = mapVtx -> get(iVtxToGrab);
1142  }
1143 
1144  // check if good
1145  if (!vtx) {
1146  cerr << PHWHERE
1147  << "PANIC: no vertex!"
1148  << endl;
1149  assert(vtx);
1150  }
1151  return vtx;
1152 
1153  } // end 'GetGlobalVertex(PHCompositeNode*, int)'
1154 
1155 
1156 
1157  GlobalVertexMap* SCorrelatorJetTree::GetVertexMap(PHCompositeNode* topNode) {
1158 
1159  // print debug statement
1160  if (m_doDebug) {
1161  cout << "SCorrelatorJetTree::GetVertexMap(PHCompositeNode*) Getting global vertex map..." << endl;
1162  }
1163 
1164  // get vertex map
1165  GlobalVertexMap* mapVtx = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1166 
1167  // check if good
1168  const bool isVtxMapGood = (mapVtx && !(mapVtx -> empty()));
1169  if (!isVtxMapGood) {
1170  cerr << PHWHERE
1171  << "PANIC: GlobalVertexMap node is missing or empty!\n"
1172  << " Please turn on the do_global flag in the main macro in order to reconstruct the global vertex!"
1173  << endl;
1174  assert(isVtxMapGood);
1175  }
1176  return mapVtx;
1177 
1178  } // end 'GetVertexMap(PHCompositeNode*)'
1179 
1180 
1181 
1182  HepMC::GenEvent* SCorrelatorJetTree::GetMcEvent(PHCompositeNode* topNode, const int iEvtToGrab) {
1183 
1184  // print debug statement
1185  if (m_doDebug) {
1186  cout << "SCorrelatorJetTree::GetMcEvent(PHCompositeNode*, int) Grabbing mc subevent #" << iEvtToGrab << "..." << endl;
1187  }
1188 
1189  // grab mc event map
1190  PHHepMCGenEventMap* mapMcEvts = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
1191  if (!mapMcEvts) {
1192  cerr << PHWHERE
1193  << "PANIC: HEPMC event map node is missing!"
1194  << endl;
1195  assert(mapMcEvts);
1196  }
1197 
1198  // grab mc event & check if good
1199  PHHepMCGenEvent* mcEvtStart = mapMcEvts -> get(iEvtToGrab);
1200  if (!mcEvtStart) {
1201  cerr << PHWHERE
1202  << "PANIC: Couldn't grab start of mc events!"
1203  << endl;
1204  assert(mcEvtStart);
1205  }
1206 
1207  HepMC::GenEvent* mcEvt = mcEvtStart -> getEvent();
1208  if (!mcEvt) {
1209  cerr << PHWHERE
1210  << "PANIC: Couldn't grab HepMC event!"
1211  << endl;
1212  assert(mcEvt);
1213  }
1214  return mcEvt;
1215 
1216  } // end 'GetMcEvent(PHCompositeNode*, int)'
1217 
1218 
1219 
1220  RawClusterContainer* SCorrelatorJetTree::GetClusterStore(PHCompositeNode* topNode, const TString sNodeName) {
1221 
1222  // print debug statement
1223  if (m_doDebug) {
1224  cout << "SCorrelatorJetTree::GetClusterStore(PHCompositeNode*, TString) Grabbing calorimeter cluster container..." << endl;
1225  }
1226 
1227  // grab clusters
1228  RawClusterContainer *clustStore = findNode::getClass<RawClusterContainer>(topNode, sNodeName.Data());
1229  if (!clustStore) {
1230  cout << PHWHERE
1231  << "PANIC: " << sNodeName.Data() << " node is missing!"
1232  << endl;
1233  assert(clustStore);
1234  }
1235  return clustStore;
1236 
1237  } // end 'GetClusterStore(PHCompositeNode*, TString)'
1238 
1239 
1240 
1241  ParticleFlowElementContainer* SCorrelatorJetTree::GetFlowStore(PHCompositeNode* topNode) {
1242 
1243  // print debug statement
1244  if (m_doDebug) {
1245  cout << "SCorrelatorJetTree::GetFlowStore(PHCompositeNode*) Grabbing particle flow container..." << endl;
1246  }
1247 
1248  // declare pf objects
1249  ParticleFlowElementContainer* flowStore = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
1250  if (!flowStore) {
1251  cerr << PHWHERE
1252  << "PANIC: Couldn't grab particle flow container!"
1253  << endl;
1254  assert(flowStore);
1255  }
1256  return flowStore;
1257 
1258  } // end 'GetFlowStore(PHCompositeNode*)'
1259 
1260 } // end SColdQcdCorrelatorAnalysis namespace
1261 
1262 // end ------------------------------------------------------------------------