Bug Summary

File:/volatile/halld/gluex/nightly/2024-04-23/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/ANALYSIS/DHistogramActions_Thrown.cc
Location:line 483, column 79
Description:Called C++ object pointer is uninitialized

Annotated Source Code

1#include "ANALYSIS/DHistogramActions.h"
2
3
4void DHistogramAction_ParticleComboGenReconComparison::Initialize(JEventLoop* locEventLoop)
5{
6 if(Get_UseKinFitResultsFlag() && (Get_Reaction()->Get_KinFitType() == d_NoFit))
7 {
8 cout << "WARNING: REQUESTED HISTOGRAM OF KINEMAITIC FIT RESULTS WHEN NO KINEMATIC FIT!!!" << endl;
9 return; //no fit performed, but kinfit data requested!!
10 }
11
12 Run_Update(locEventLoop);
13
14 vector<const DParticleID*> locParticleIDs;
15 locEventLoop->Get(locParticleIDs);
16
17 string locHistName, locHistTitle, locStepName, locStepROOTName, locParticleName, locParticleROOTName;
18 Particle_t locPID;
19
20 size_t locNumSteps = Get_Reaction()->Get_NumReactionSteps();
21
22 dHistDeque_DeltaPOverP.resize(locNumSteps);
23 dHistDeque_DeltaTheta.resize(locNumSteps);
24 dHistDeque_DeltaPhi.resize(locNumSteps);
25 dHistDeque_DeltaT.resize(locNumSteps);
26 dHistDeque_DeltaT_TOF.resize(locNumSteps);
27 dHistDeque_DeltaT_BCAL.resize(locNumSteps);
28 dHistDeque_DeltaT_FCAL.resize(locNumSteps);
29 dHistDeque_DeltaVertexZ.resize(locNumSteps);
30 dHistDeque_DeltaPOverPVsP.resize(locNumSteps);
31 dHistDeque_DeltaPOverPVsTheta.resize(locNumSteps);
32 dHistDeque_DeltaThetaVsP.resize(locNumSteps);
33 dHistDeque_DeltaThetaVsTheta.resize(locNumSteps);
34 dHistDeque_DeltaPhiVsP.resize(locNumSteps);
35 dHistDeque_DeltaPhiVsTheta.resize(locNumSteps);
36 dHistDeque_DeltaTVsTheta.resize(locNumSteps);
37 dHistDeque_DeltaTVsP.resize(locNumSteps);
38 dHistDeque_DeltaVertexZVsTheta.resize(locNumSteps);
39
40 dHistDeque_Pulls.resize(locNumSteps);
41 dHistDeque_PullsVsP.resize(locNumSteps);
42 dHistDeque_PullsVsTheta.resize(locNumSteps);
43
44 dHistDeque_TimePull_CDC.resize(locNumSteps);
45 dHistDeque_TimePull_ST.resize(locNumSteps);
46 dHistDeque_TimePull_BCAL.resize(locNumSteps);
47 dHistDeque_TimePull_TOF.resize(locNumSteps);
48 dHistDeque_TimePull_FCAL.resize(locNumSteps);
49
50 dHistDeque_TimePullVsTheta_CDC.resize(locNumSteps);
51 dHistDeque_TimePullVsTheta_BCAL.resize(locNumSteps);
52 dHistDeque_TimePullVsTheta_ST.resize(locNumSteps);
53
54 dHistDeque_TimePullVsP_CDC.resize(locNumSteps);
55 dHistDeque_TimePullVsP_BCAL.resize(locNumSteps);
56 dHistDeque_TimePullVsP_ST.resize(locNumSteps);
57 dHistDeque_TimePullVsP_TOF.resize(locNumSteps);
58 dHistDeque_TimePullVsP_FCAL.resize(locNumSteps);
59
60 //CREATE THE HISTOGRAMS
61 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
62 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
63 {
64 CreateAndChangeTo_ActionDirectory();
65
66 //RF
67 locHistName = "DeltaT_RFBeamBunch";
68 dRFBeamBunchDeltaT_Hist = GetOrCreate_Histogram<TH1I>(locHistName, ";RF #Deltat (Reconstructed - Thrown)", dNumRFDeltaTBins, dMinRFDeltaT, dMaxRFDeltaT);
69
70 //beam particle
71 locPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID();
72 bool locBeamParticleUsed = (locPID == Gamma);
73 if(locBeamParticleUsed)
74 {
75 locParticleName = string("Beam_") + ParticleType(locPID);
76 CreateAndChangeTo_Directory(locParticleName, locParticleName);
77 locParticleROOTName = ParticleName_ROOT(locPID);
78
79 // DeltaP/P
80 locHistName = string("DeltaPOverP");
81 locHistTitle = locParticleROOTName + string(";#Deltap/p (Reconstructed - Thrown)");
82 dBeamParticleHist_DeltaPOverP = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
83
84 // DeltaP/P Vs P
85 locHistName = string("DeltaPOverPVsP");
86 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltap/p (Reconstructed - Thrown)");
87 dBeamParticleHist_DeltaPOverPVsP = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
88
89 // DeltaT
90 locHistName = string("DeltaT");
91 locHistTitle = locParticleROOTName + string(";#Deltat (ns) (Reconstructed - Thrown)");
92 dBeamParticleHist_DeltaT = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
93
94 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
95 }
96
97 deque<string> locPullNames(8, "");
98 locPullNames[0] = "E"; locPullNames[1] = "Px"; locPullNames[2] = "Py"; locPullNames[3] = "Pz";
99 locPullNames[4] = "Xx"; locPullNames[5] = "Xy"; locPullNames[6] = "Xz"; locPullNames[7] = "T";
100
101 deque<string> locPullTitles(8, "");
102 locPullTitles[0] = "E"; locPullTitles[1] = "p_{x}"; locPullTitles[2] = "p_{y}"; locPullTitles[3] = "p_{z}";
103 locPullTitles[4] = "x_{x}"; locPullTitles[5] = "x_{y}"; locPullTitles[6] = "x_{z}"; locPullTitles[7] = "t";
104
105 //CREATE THE HISTOGRAMS
106 deque<Particle_t> locPIDs;
107 for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)
108 {
109 const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
110 locStepName = DAnalysis::Get_StepName(locReactionStep, true, false);
111 locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);
112
113 Particle_t locInitialPID = locReactionStep->Get_InitialPID();
114
115 //get PIDs
116 auto locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, Get_UseKinFitResultsFlag(), Get_UseKinFitResultsFlag(), d_AllCharges, false);
117 if(Get_UseKinFitResultsFlag())
118 {
119 if((!locBeamParticleUsed) || (loc_i != 0)) //add decaying parent particle //skip if on beam particle!
120 locPIDs.insert(locPIDs.begin(), locInitialPID);
121 }
122 if(locPIDs.empty())
123 continue;
124
125 CreateAndChangeTo_Directory(locStepName, locStepName);
126 for(size_t loc_j = 0; loc_j < locPIDs.size(); ++loc_j)
127 {
128 locPID = locPIDs[loc_j];
129 locParticleName = ParticleType(locPID);
130 locParticleROOTName = ParticleName_ROOT(locPID);
131 if(dHistDeque_DeltaPOverP[loc_i].find(locPID) != dHistDeque_DeltaPOverP[loc_i].end())
132 continue; //pid already done
133
134 CreateAndChangeTo_Directory(locParticleName, locParticleName);
135
136 // DeltaP/P
137 locHistName = string("DeltaPOverP");
138 locHistTitle = locParticleROOTName + string(";#Deltap/p (Reconstructed - Thrown)");
139 dHistDeque_DeltaPOverP[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
140
141 // DeltaTheta
142 locHistName = string("DeltaTheta");
143 locHistTitle = locParticleROOTName + string(";#Delta#theta#circ (Reconstructed - Thrown)");
144 dHistDeque_DeltaTheta[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);
145
146 // DeltaPhi
147 locHistName = string("DeltaPhi");
148 locHistTitle = locParticleROOTName + string(";#Delta#phi#circ (Reconstructed - Thrown)");
149 dHistDeque_DeltaPhi[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
150
151 // DeltaT
152 locHistName = string("DeltaT");
153 locHistTitle = locParticleROOTName + string(";#Deltat (ns) (Reconstructed - Thrown)");
154 dHistDeque_DeltaT[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
155
156 // DeltaT - BCAL
157 locHistName = string("DeltaT_BCAL");
158 locHistTitle = locParticleROOTName + string(" in BCAL;#Deltat (ns) (Reconstructed - Thrown)");
159 dHistDeque_DeltaT_BCAL[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
160
161 // DeltaT - TOF (charged only)
162 if(ParticleCharge(locPID) != 0)
163 {
164 locHistName = string("DeltaT_TOF");
165 locHistTitle = locParticleROOTName + string(" in TOF;#Deltat (ns) (Reconstructed - Thrown)");
166 dHistDeque_DeltaT_TOF[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
167 }
168
169 // DeltaT - FCAL (neutral only)
170 if(ParticleCharge(locPID) == 0)
171 {
172 locHistName = string("DeltaT_FCAL");
173 locHistTitle = locParticleROOTName + string(" in FCAL;#Deltat (ns) (Reconstructed - Thrown)");
174 dHistDeque_DeltaT_FCAL[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
175 }
176
177 // DeltaVertexZ
178 locHistName = string("DeltaVertexZ");
179 locHistTitle = locParticleROOTName + string(";#DeltaVertex-Z (cm) (Reconstructed - Thrown)");
180 dHistDeque_DeltaVertexZ[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexZBins, dMinDeltaVertexZ, dMaxDeltaVertexZ);
181
182 // DeltaP/P Vs P
183 locHistName = string("DeltaPOverPVsP");
184 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltap/p (Reconstructed - Thrown)");
185 dHistDeque_DeltaPOverPVsP[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
186
187 // DeltaP/P Vs Theta
188 locHistName = string("DeltaPOverPVsTheta");
189 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltap/p (Reconstructed - Thrown)");
190 dHistDeque_DeltaPOverPVsTheta[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
191
192 // DeltaTheta Vs P
193 locHistName = string("DeltaThetaVsP");
194 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#theta#circ (Reconstructed - Thrown)");
195 dHistDeque_DeltaThetaVsP[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);
196
197 // DeltaTheta Vs Theta
198 locHistName = string("DeltaThetaVsTheta");
199 locHistTitle = locParticleROOTName + string(";#theta#circ;#Delta#theta#circ (Reconstructed - Thrown)");
200 dHistDeque_DeltaThetaVsTheta[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);
201
202 // DeltaPhi Vs P
203 locHistName = string("DeltaPhiVsP");
204 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#phi#circ (Reconstructed - Thrown)");
205 dHistDeque_DeltaPhiVsP[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
206
207 // DeltaPhi Vs Theta
208 locHistName = string("DeltaPhiVsTheta");
209 locHistTitle = locParticleROOTName + string(";#theta#circ;#Delta#phi#circ (Reconstructed - Thrown)");
210 dHistDeque_DeltaPhiVsTheta[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
211
212 // DeltaT Vs Theta
213 locHistName = string("DeltaTVsTheta");
214 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat (ns) (Reconstructed - Thrown)");
215 dHistDeque_DeltaTVsTheta[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
216
217 // DeltaT Vs P
218 locHistName = string("DeltaTVsP");
219 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat (ns) (Reconstructed - Thrown)");
220 dHistDeque_DeltaTVsP[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
221
222 // DeltaVertexZ Vs Theta
223 locHistName = string("DeltaVertexZVsTheta");
224 locHistTitle = locParticleROOTName + string(";#theta#circ;#DeltaVertex-Z (cm) (Reconstructed - Thrown)");
225 dHistDeque_DeltaVertexZVsTheta[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaVertexZBins, dMinDeltaVertexZ, dMaxDeltaVertexZ);
226
227 /************************************************************************ Pulls ************************************************************************/
228
229 CreateAndChangeTo_Directory("Pulls", "Pulls");
230 for(size_t loc_j = 0; loc_j < dPullTypes.size(); ++loc_j)
231 {
232 if((ParticleCharge(locPID) != 0) && (dPullTypes[loc_j] == d_EPull))
233 continue;
234 if((ParticleCharge(locPID) == 0) && ((dPullTypes[loc_j] >= d_PxPull) && (dPullTypes[loc_j] <= d_PzPull)))
235 continue;
236
237 //Pull 1D
238 locHistName = locPullNames[loc_j] + string("Pull");
239 locHistTitle = locParticleROOTName + string(";#Delta") + locPullTitles[loc_j] + string("/#sigma_{") + locPullTitles[loc_j] + string("} (Reconstructed - Thrown)");
240 dHistDeque_Pulls[loc_i][locPID][dPullTypes[loc_j]] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
241
242 //Pull vs P
243 locHistName = locPullNames[loc_j] + string("PullVsP");
244 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta") + locPullTitles[loc_j] + string("/#sigma_{") + locPullTitles[loc_j] + string("} (Reconstructed - Thrown)");
245 dHistDeque_PullsVsP[loc_i][locPID][dPullTypes[loc_j]] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
246
247 //Pull vs Theta
248 locHistName = locPullNames[loc_j] + string("PullVsTheta");
249 locHistTitle = locParticleROOTName + string(";#theta#circ;#Delta") + locPullTitles[loc_j] + string("/#sigma_{") + locPullTitles[loc_j] + string("} (Reconstructed - Thrown)");
250 dHistDeque_PullsVsTheta[loc_i][locPID][dPullTypes[loc_j]] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
251 }
252
253 //Delta-t Pulls - CDC & ST
254 if(ParticleCharge(locPID) != 0)
255 {
256 //CDC
257 locHistName = "TimePull_CDC";
258 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
259 dHistDeque_TimePull_CDC[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
260
261 locHistName = "TimePullVsTheta_CDC";
262 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat/#sigma_{#Deltat}");
263 dHistDeque_TimePullVsTheta_CDC[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
264
265 locHistName = "TimePullVsP_CDC";
266 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
267 dHistDeque_TimePullVsP_CDC[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
268
269 //ST
270 locHistName = "TimePull_ST";
271 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
272 dHistDeque_TimePull_ST[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
273
274 locHistName = "TimePullVsTheta_ST";
275 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat/#sigma_{#Deltat}");
276 dHistDeque_TimePullVsTheta_ST[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
277
278 locHistName = "TimePullVsP_ST";
279 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
280 dHistDeque_TimePullVsP_ST[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
281 }
282
283 //Delta-t Pulls - BCAL
284 locHistName = "TimePull_BCAL";
285 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
286 dHistDeque_TimePull_BCAL[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
287
288 locHistName = "TimePullVsTheta_BCAL";
289 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat/#sigma_{#Deltat}");
290 dHistDeque_TimePullVsTheta_BCAL[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
291
292 locHistName = "TimePullVsP_BCAL";
293 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
294 dHistDeque_TimePullVsP_BCAL[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
295
296 //Delta-t Pulls - TOF
297 if(ParticleCharge(locPID) != 0) //TOF
298 {
299 locHistName = "TimePull_TOF";
300 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
301 dHistDeque_TimePull_TOF[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
302
303 locHistName = "TimePullVsP_TOF";
304 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
305 dHistDeque_TimePullVsP_TOF[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
306 }
307
308 //Delta-t Pulls - FCAL
309 locHistName = "TimePull_FCAL";
310 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
311 dHistDeque_TimePull_FCAL[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
312
313 locHistName = "TimePullVsP_FCAL";
314 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
315 dHistDeque_TimePullVsP_FCAL[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
316
317 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
318
319 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
320 } //end of particle loop
321 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
322 } //end of step loop
323
324 //Return to the base directory
325 ChangeTo_BaseDirectory();
326 }
327 japp->RootUnLock(); //RELEASE ROOT LOCK!!
328}
329
330void DHistogramAction_ParticleComboGenReconComparison::Run_Update(JEventLoop* locEventLoop)
331{
332 DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
333 DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
334 locGeometry->GetTargetZ(dTargetZCenter);
335}
336
337bool DHistogramAction_ParticleComboGenReconComparison::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
338{
339 if(Get_UseKinFitResultsFlag() && (Get_Reaction()->Get_KinFitType() == d_NoFit))
340 {
341 cout << "WARNING: REQUESTED HISTOGRAM OF KINEMAITIC FIT RESULTS WHEN NO KINEMATIC FIT!!! Skipping histogram." << endl;
342 return true; //no fit performed, but kinfit data requested!!
343 }
344
345 vector<const DMCThrown*> locMCThrowns;
346 locEventLoop->Get(locMCThrowns);
347 if(locMCThrowns.empty())
348 return true; //e.g. non-simulated event
349
350 vector<const DBeamPhoton*> locBeamPhotons;
351 locEventLoop->Get(locBeamPhotons, "MCGEN");
352 const DBeamPhoton* locThrownBeamPhoton = locBeamPhotons.empty() ? NULL__null : locBeamPhotons[0];
353
354 vector<const DMCThrownMatching*> locMCThrownMatchingVector;
355 locEventLoop->Get(locMCThrownMatchingVector);
356 if(locMCThrownMatchingVector.empty())
357 return true;
358 const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0];
359
360 const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();
361 const DEventRFBunch* locThrownEventRFBunch = NULL__null;
362 locEventLoop->GetSingle(locThrownEventRFBunch, "Thrown");
363
364 //RF time difference
365 double locRFTime = locEventRFBunch->dTime;
366 double locRFDeltaT = locRFTime - locThrownEventRFBunch->dTime;
367 //FILL HISTOGRAMS
368 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
369 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
370 Lock_Action();
371 {
372 dRFBeamBunchDeltaT_Hist->Fill(locRFDeltaT);
373 }
374 Unlock_Action();
375
376 const DKinematicData* locKinematicData;
377 for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
378 {
379 const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);
380
381 //initial particle
382 if(Get_UseKinFitResultsFlag())
383 locKinematicData = locParticleComboStep->Get_InitialParticle();
384 else
385 locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();
386 if(locKinematicData != NULL__null)
387 {
388 if(locKinematicData->PID() == Gamma)
389 {
390 //check if will be duplicate
391 const JObject* locSourceObject = locParticleComboStep->Get_InitialParticle_Measured();
392 if(dPreviouslyHistogrammedBeamParticles.find(locSourceObject) == dPreviouslyHistogrammedBeamParticles.end())
393 {
394 dPreviouslyHistogrammedBeamParticles.insert(locSourceObject);
395 Fill_BeamHists(locKinematicData, locThrownBeamPhoton);
396 }
397 }
398 }
399
400 //final particles
401 for(size_t loc_j = 0; loc_j < locParticleComboStep->Get_NumFinalParticles(); ++loc_j)
402 {
403 if(Get_UseKinFitResultsFlag())
404 locKinematicData = locParticleComboStep->Get_FinalParticle(loc_j);
405 else
406 locKinematicData = locParticleComboStep->Get_FinalParticle_Measured(loc_j);
407 if(locKinematicData == NULL__null)
408 continue; //e.g. a decaying or missing particle whose params aren't set yet
409 if(DAnalysis::Get_DecayStepIndex(Get_Reaction(), loc_i, loc_j) > 0)
410 continue; //decaying: not supported
411 if(Get_Reaction()->Get_ReactionStep(loc_i)->Get_MissingParticleIndex() == int(loc_j))
412 continue; //missing: not supported
413
414 //check if duplicate
415 const JObject* locSourceObject = locParticleComboStep->Get_FinalParticle_SourceObject(loc_j);
416 if(locSourceObject != NULL__null) //else is reconstructed missing/decaying particle: has many source object, and is unique to this combo: no dupes to check against: let it ride
417 {
418 pair<Particle_t, const JObject*> locParticleInfo(locKinematicData->PID(), locSourceObject);
419 pair<size_t, pair<Particle_t, const JObject*> > locHistInfo(loc_i, locParticleInfo);
420 if(dPreviouslyHistogrammedParticles.find(locHistInfo) != dPreviouslyHistogrammedParticles.end())
421 continue; //previously histogrammed
422 dPreviouslyHistogrammedParticles.insert(locHistInfo);
423 }
424
425 if(ParticleCharge(locKinematicData->PID()) != 0)
426 {
427 const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locKinematicData);
428 double locMatchFOM = 0.0;
429 const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);
430 if(locMCThrown != NULL__null)
431 Fill_ChargedHists(locChargedTrackHypothesis, locMCThrown, locThrownEventRFBunch, loc_i);
432 }
433 else
434 {
435 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locKinematicData);
436 double locMatchFOM = 0.0;
437 const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM);
438 if(locMCThrown != NULL__null)
439 Fill_NeutralHists(locNeutralParticleHypothesis, locMCThrown, locThrownEventRFBunch, loc_i);
440 }
441 } //end of particle loop
442 } //end of step loop
443 return true;
444}
445
446void DHistogramAction_ParticleComboGenReconComparison::Fill_BeamHists(const DKinematicData* locKinematicData, const DKinematicData* locThrownKinematicData)
447{
448 if(locThrownKinematicData == NULL__null)
449 return;
450
451 DVector3 locMomentum = locKinematicData->momentum();
452 DVector3 locThrownMomentum = locThrownKinematicData->momentum();
453
454 double locThrownP = locThrownMomentum.Mag();
455 double locDeltaPOverP = (locMomentum.Mag() - locThrownP)/locThrownP;
456 double locDeltaT = locKinematicData->time() - locThrownKinematicData->time();
457
458 //FILL HISTOGRAMS
459 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
460 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
461 Lock_Action();
462 {
463 dBeamParticleHist_DeltaPOverP->Fill(locDeltaPOverP);
464 dBeamParticleHist_DeltaPOverPVsP->Fill(locThrownP, locDeltaPOverP);
465 dBeamParticleHist_DeltaT->Fill(locDeltaT);
466 }
467 Unlock_Action();
468}
469
470void DHistogramAction_ParticleComboGenReconComparison::Fill_ChargedHists(const DChargedTrackHypothesis* locChargedTrackHypothesis, const DMCThrown* locMCThrown, const DEventRFBunch* locThrownEventRFBunch, size_t locStepIndex)
471{
472 Particle_t locPID = locChargedTrackHypothesis->PID();
473 double locThrownP = locMCThrown->momentum().Mag();
474 double locThrownTheta = locMCThrown->momentum().Theta()*180.0/TMath::Pi();
475 double locDeltaPOverP = (locChargedTrackHypothesis->momentum().Mag() - locThrownP)/locThrownP;
476 double locDeltaTheta = locChargedTrackHypothesis->momentum().Theta()*180.0/TMath::Pi() - locThrownTheta;
477 double locDeltaPhi = locChargedTrackHypothesis->momentum().Phi()*180.0/TMath::Pi() - locMCThrown->momentum().Phi()*180.0/TMath::Pi();
478 double locDeltaT = locChargedTrackHypothesis->time() - locMCThrown->time(); //time comparison isn't fair if track comes from a detached vertex!!!
479 double locDeltaVertexZ = locChargedTrackHypothesis->position().Z() - locMCThrown->position().Z();
480 const TMatrixDSym& locCovarianceMatrix = *(locChargedTrackHypothesis->errorMatrix().get());
1
'locCovarianceMatrix' initialized to a garbage value
→
481
482 double locStartTime = locThrownEventRFBunch->dTime + (locMCThrown->z() - dTargetZCenter)/29.9792458;
483 double locTimePull = (locStartTime - locChargedTrackHypothesis->time())/sqrt(locCovarianceMatrix(6, 6));
2
←
Called C++ object pointer is uninitialized
484 double locT0Pull = (locStartTime - locChargedTrackHypothesis->t0())/locChargedTrackHypothesis->t0_err();
485
486 //FILL HISTOGRAMS
487 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
488 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
489 Lock_Action();
490 {
491 dHistDeque_DeltaPOverP[locStepIndex][locPID]->Fill(locDeltaPOverP);
492 dHistDeque_DeltaTheta[locStepIndex][locPID]->Fill(locDeltaTheta);
493 dHistDeque_DeltaPhi[locStepIndex][locPID]->Fill(locDeltaPhi);
494 dHistDeque_DeltaT[locStepIndex][locPID]->Fill(locDeltaT);
495 if(locChargedTrackHypothesis->t0_detector() == SYS_START)
496 {
497 dHistDeque_TimePull_ST[locStepIndex][locPID]->Fill(locT0Pull);
498 dHistDeque_TimePullVsTheta_ST[locStepIndex][locPID]->Fill(locThrownTheta, locT0Pull);
499 dHistDeque_TimePullVsP_ST[locStepIndex][locPID]->Fill(locThrownP, locT0Pull);
500 }
501 if(locChargedTrackHypothesis->t0_detector() == SYS_CDC)
502 {
503 dHistDeque_TimePull_CDC[locStepIndex][locPID]->Fill(locT0Pull);
504 dHistDeque_TimePullVsTheta_CDC[locStepIndex][locPID]->Fill(locThrownTheta, locT0Pull);
505 dHistDeque_TimePullVsP_CDC[locStepIndex][locPID]->Fill(locThrownP, locT0Pull);
506 }
507 else if(locChargedTrackHypothesis->t1_detector() == SYS_CDC)
508 {
509 dHistDeque_TimePull_CDC[locStepIndex][locPID]->Fill(locTimePull);
510 dHistDeque_TimePullVsTheta_CDC[locStepIndex][locPID]->Fill(locThrownTheta, locTimePull);
511 dHistDeque_TimePullVsP_CDC[locStepIndex][locPID]->Fill(locThrownP, locTimePull);
512 }
513 if(locChargedTrackHypothesis->t1_detector() == SYS_BCAL)
514 {
515 dHistDeque_DeltaT_BCAL[locStepIndex][locPID]->Fill(locDeltaT);
516 dHistDeque_TimePull_BCAL[locStepIndex][locPID]->Fill(locTimePull);
517 dHistDeque_TimePullVsTheta_BCAL[locStepIndex][locPID]->Fill(locThrownTheta, locTimePull);
518 dHistDeque_TimePullVsP_BCAL[locStepIndex][locPID]->Fill(locThrownP, locTimePull);
519 }
520 else if(locChargedTrackHypothesis->t1_detector() == SYS_TOF)
521 {
522 dHistDeque_DeltaT_TOF[locStepIndex][locPID]->Fill(locDeltaT);
523 dHistDeque_TimePull_TOF[locStepIndex][locPID]->Fill(locTimePull);
524 dHistDeque_TimePullVsP_TOF[locStepIndex][locPID]->Fill(locThrownP, locTimePull);
525 }
526 else if(locChargedTrackHypothesis->t1_detector() == SYS_FCAL)
527 {
528 dHistDeque_TimePull_FCAL[locStepIndex][locPID]->Fill(locTimePull);
529 dHistDeque_TimePullVsP_FCAL[locStepIndex][locPID]->Fill(locThrownP, locTimePull);
530 }
531 dHistDeque_DeltaVertexZ[locStepIndex][locPID]->Fill(locDeltaVertexZ);
532 dHistDeque_DeltaPOverPVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaPOverP);
533 dHistDeque_DeltaPOverPVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaPOverP);
534 dHistDeque_DeltaThetaVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaTheta);
535 dHistDeque_DeltaThetaVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaTheta);
536 dHistDeque_DeltaPhiVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaPhi);
537 dHistDeque_DeltaPhiVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaPhi);
538 dHistDeque_DeltaTVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaT);
539 dHistDeque_DeltaTVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaT);
540 dHistDeque_DeltaVertexZVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaVertexZ);
541
542 for(size_t loc_j = 0; loc_j < dPullTypes.size(); ++loc_j)
543 {
544 if(dPullTypes[loc_j] == d_EPull)
545 continue;
546 double locPull = 0.0;
547 if((dPullTypes[loc_j] >= d_PxPull) && (dPullTypes[loc_j] <= d_PzPull))
548 {
549 int locIndex = int(dPullTypes[loc_j] - d_PxPull);
550 locPull = (locChargedTrackHypothesis->momentum()(locIndex) - locMCThrown->momentum()(locIndex))/sqrt(locCovarianceMatrix(locIndex, locIndex));
551 }
552 else if((dPullTypes[loc_j] >= d_XxPull) && (dPullTypes[loc_j] <= d_XzPull))
553 {
554 int locIndex = int(dPullTypes[loc_j] - d_XxPull);
555 locPull = (locChargedTrackHypothesis->position()(locIndex) - locMCThrown->position()(locIndex))/sqrt(locCovarianceMatrix(locIndex + 3, locIndex + 3));
556 }
557 else if(dPullTypes[loc_j] == d_TPull)
558 locPull = (locChargedTrackHypothesis->time() - locMCThrown->time())/sqrt(locCovarianceMatrix(6, 6));
559 (dHistDeque_Pulls[locStepIndex][locPID])[dPullTypes[loc_j]]->Fill(locPull);
560 (dHistDeque_PullsVsP[locStepIndex][locPID])[dPullTypes[loc_j]]->Fill(locThrownP, locPull);
561 (dHistDeque_PullsVsTheta[locStepIndex][locPID])[dPullTypes[loc_j]]->Fill(locThrownTheta, locPull);
562 }
563 }
564 Unlock_Action();
565}
566
567void DHistogramAction_ParticleComboGenReconComparison::Fill_NeutralHists(const DNeutralParticleHypothesis* locNeutralParticleHypothesis, const DMCThrown* locMCThrown, const DEventRFBunch* locThrownEventRFBunch, size_t locStepIndex)
568{
569 Particle_t locPID = locNeutralParticleHypothesis->PID();
570 const DNeutralShower* locNeutralShower = locNeutralParticleHypothesis->Get_NeutralShower();
571
572 double locThrownP = locMCThrown->momentum().Mag();
573 double locThrownTheta = locMCThrown->momentum().Theta()*180.0/TMath::Pi();
574 DVector3 locNeutralP3 = locNeutralParticleHypothesis->momentum();
575 double locPMag = locNeutralP3.Mag();
576 double locDeltaPOverP = (locPMag - locThrownP)/locThrownP;
577 double locDeltaTheta = locNeutralP3.Theta()*180.0/TMath::Pi() - locThrownTheta;
578 double locDeltaPhi = locNeutralP3.Phi()*180.0/TMath::Pi() - locMCThrown->momentum().Phi()*180.0/TMath::Pi();
579 double locDeltaT = locNeutralParticleHypothesis->time() - locMCThrown->time(); //time comparison isn't fair if track comes from a detached vertex!!!
580 double locDeltaVertexZ = locNeutralParticleHypothesis->position().Z() - locMCThrown->position().Z();
581 auto locCovarianceMatrix = locNeutralParticleHypothesis->errorMatrix();
582
583 double locStartTime = locThrownEventRFBunch->dTime + (locMCThrown->z() - dTargetZCenter)/29.9792458;
584 double locTimePull = (locCovarianceMatrix != nullptr) ? (locStartTime - locNeutralParticleHypothesis->time())/sqrt((*locCovarianceMatrix)(6, 6)) : std::numeric_limits<double>::quiet_NaN();
585
586 double locEPull = 0.0;
587 if(Get_UseKinFitResultsFlag())
588 {
589 DMatrix locJacobian(1, 7);
590 locJacobian(0, 0) = locNeutralP3.Px()/locPMag;
591 locJacobian(0, 1) = locNeutralP3.Py()/locPMag;
592 locJacobian(0, 2) = locNeutralP3.Pz()/locPMag;
593 for(unsigned int loc_i = 0; loc_i < 4; ++loc_i)
594 locJacobian(0, 3 + loc_i) = 0.0;
595
596 if(locCovarianceMatrix != nullptr)
597 {
598 TMatrixDSym locCovCopy = *locCovarianceMatrix;
599 locCovCopy.Similarity(locJacobian);
600 double locEUncertainty = sqrt(locCovCopy(0, 0));
601 locEPull = (locNeutralParticleHypothesis->energy() - locMCThrown->energy())/locEUncertainty;
602 }
603 else
604 locEPull = std::numeric_limits<double>::quiet_NaN();
605 }
606 else
607 locEPull = (locNeutralShower->dEnergy - locMCThrown->energy())/sqrt((*(locNeutralShower->dCovarianceMatrix))(0, 0));
608
609 //FILL HISTOGRAMS
610 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
611 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
612 Lock_Action();
613 {
614 dHistDeque_DeltaPOverP[locStepIndex][locPID]->Fill(locDeltaPOverP);
615 dHistDeque_DeltaTheta[locStepIndex][locPID]->Fill(locDeltaTheta);
616 dHistDeque_DeltaPhi[locStepIndex][locPID]->Fill(locDeltaPhi);
617 dHistDeque_DeltaT[locStepIndex][locPID]->Fill(locDeltaT);
618 if(locNeutralParticleHypothesis->t1_detector() == SYS_BCAL)
619 {
620 dHistDeque_DeltaT_BCAL[locStepIndex][locPID]->Fill(locDeltaT);
621 dHistDeque_TimePull_BCAL[locStepIndex][locPID]->Fill(locTimePull);
622 dHistDeque_TimePullVsTheta_BCAL[locStepIndex][locPID]->Fill(locThrownTheta, locTimePull);
623 dHistDeque_TimePullVsP_BCAL[locStepIndex][locPID]->Fill(locThrownP, locTimePull);
624 }
625 else if(locNeutralParticleHypothesis->t1_detector() == SYS_FCAL)
626 {
627 dHistDeque_DeltaT_FCAL[locStepIndex][locPID]->Fill(locDeltaT);
628 dHistDeque_TimePull_FCAL[locStepIndex][locPID]->Fill(locTimePull);
629 dHistDeque_TimePullVsP_FCAL[locStepIndex][locPID]->Fill(locThrownP, locTimePull);
630 }
631
632 dHistDeque_DeltaVertexZ[locStepIndex][locPID]->Fill(locDeltaVertexZ);
633 dHistDeque_DeltaPOverPVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaPOverP);
634 dHistDeque_DeltaPOverPVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaPOverP);
635 dHistDeque_DeltaThetaVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaTheta);
636 dHistDeque_DeltaThetaVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaTheta);
637 dHistDeque_DeltaPhiVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaPhi);
638 dHistDeque_DeltaPhiVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaPhi);
639 dHistDeque_DeltaTVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaT);
640 dHistDeque_DeltaTVsP[locStepIndex][locPID]->Fill(locThrownP, locDeltaT);
641 dHistDeque_DeltaVertexZVsTheta[locStepIndex][locPID]->Fill(locThrownTheta, locDeltaVertexZ);
642
643 if(locCovarianceMatrix != nullptr)
644 {
645 for(size_t loc_j = 0; loc_j < dPullTypes.size(); ++loc_j)
646 {
647 if((dPullTypes[loc_j] >= d_PxPull) && (dPullTypes[loc_j] <= d_PzPull))
648 continue;
649 double locPull = 0.0;
650 if(dPullTypes[loc_j] == d_EPull)
651 locPull = locEPull;
652 else if((dPullTypes[loc_j] >= d_XxPull) && (dPullTypes[loc_j] <= d_XzPull))
653 {
654 int locIndex = int(dPullTypes[loc_j] - d_XxPull);
655 locPull = (locNeutralParticleHypothesis->position()(locIndex) - locMCThrown->position()(locIndex))/sqrt((*locCovarianceMatrix)(locIndex + 3, locIndex + 3));
656 }
657 else if(dPullTypes[loc_j] == d_TPull)
658 locPull = (locNeutralParticleHypothesis->time() - locMCThrown->time())/sqrt((*locCovarianceMatrix)(6, 6));
659 (dHistDeque_Pulls[locStepIndex][locPID])[dPullTypes[loc_j]]->Fill(locPull);
660 (dHistDeque_PullsVsP[locStepIndex][locPID])[dPullTypes[loc_j]]->Fill(locThrownP, locPull);
661 (dHistDeque_PullsVsTheta[locStepIndex][locPID])[dPullTypes[loc_j]]->Fill(locThrownTheta, locPull);
662 }
663 }
664 }
665 Unlock_Action();
666}
667
668void DHistogramAction_ThrownParticleKinematics::Initialize(JEventLoop* locEventLoop)
669{
670 string locHistName, locHistTitle, locParticleName, locParticleROOTName;
671 Particle_t locPID;
672
673 //CREATE THE HISTOGRAMS
674 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
675 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
676 {
677 CreateAndChangeTo_ActionDirectory();
678
679 // Beam Particle: MCGEN
680 {
681 locPID = Gamma;
682 locParticleName = string("MCGENBeamParticle_") + ParticleType(locPID);
683 locParticleROOTName = ParticleName_ROOT(locPID);
684 CreateAndChangeTo_Directory(locParticleName, locParticleName);
685
686 locHistName = "Momentum";
687 locHistTitle = string("MCGEN Thrown Beam ") + locParticleROOTName + string(";p (GeV/c)");
688 dMCGENBeamParticle_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
689
690 locHistName = "Time";
691 locHistTitle = string("MCGEN Thrown Beam ") + locParticleROOTName + string(";t (ns)");
692 dMCGENBeamParticle_Time = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
693
694 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
695 }
696
697 // Beam Particle: All
698 {
699 locPID = Gamma;
700 locParticleName = string("TRUTHBeamParticles_") + ParticleType(locPID);
701 locParticleROOTName = ParticleName_ROOT(locPID);
702 CreateAndChangeTo_Directory(locParticleName, locParticleName);
703
704 locHistName = "Momentum";
705 locHistTitle = string("TRUTH Thrown Beam ") + locParticleROOTName + string(";p (GeV/c)");
706 dAllBeamParticle_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
707
708 locHistName = "Time";
709 locHistTitle = string("TRUTH Thrown Beam ") + locParticleROOTName + string(";t (ns)");
710 dAllBeamParticle_Time = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
711
712 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
713 }
714
715 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
716 {
717 locPID = dFinalStatePIDs[loc_i];
718 locParticleName = ParticleType(locPID);
719 locParticleROOTName = ParticleName_ROOT(locPID);
720 CreateAndChangeTo_Directory(locParticleName, locParticleName);
721
722 // Momentum
723 locHistName = "Momentum";
724 locHistTitle = string("Thrown ") + locParticleROOTName + string(";p (GeV/c)");
725 dHistMap_P[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
726
727 // Theta
728 locHistName = "Theta";
729 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#theta#circ");
730 dHistMap_Theta[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
731
732 // Phi
733 locHistName = "Phi";
734 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#phi#circ");
735 dHistMap_Phi[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
736
737 // P Vs Theta
738 locHistName = "PVsTheta";
739 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#theta#circ;p (GeV/c)");
740 dHistMap_PVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
741
742 // Phi Vs Theta
743 locHistName = "PhiVsTheta";
744 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#theta#circ;#phi#circ");
745 dHistMap_PhiVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
746
747 // Vertex-Z
748 locHistName = "VertexZ";
749 locHistTitle = string("Thrown ") + locParticleROOTName + string(";Vertex-Z (cm)");
750 dHistMap_VertexZ[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
751
752 // Vertex-Y Vs Vertex-X
753 locHistName = "VertexYVsX";
754 locHistTitle = string("Thrown ") + locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
755 dHistMap_VertexYVsX[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
756
757 // Vertex-T
758 locHistName = "VertexT";
759 locHistTitle = string("Thrown ") + locParticleROOTName + string(";Vertex-T (ns)");
760 dHistMap_VertexT[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
761
762 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
763 }
764
765 //Return to the base directory
766 ChangeTo_BaseDirectory();
767 }
768 japp->RootUnLock(); //RELEASE ROOT LOCK!!
769}
770
771bool DHistogramAction_ThrownParticleKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
772{
773 vector<const DMCThrown*> locMCThrowns, locMCThrowns_Decaying;
774 locEventLoop->Get(locMCThrowns, "FinalState");
775 locEventLoop->Get(locMCThrowns_Decaying, "Decaying");
776 locMCThrowns.insert(locMCThrowns.begin(), locMCThrowns_Decaying.begin(), locMCThrowns_Decaying.end());
777 if(locMCThrowns.empty())
778 return true; //e.g. non-simulated event
779
780 if(Get_CalledPriorWithComboFlag())
781 return true; //else double-counting!
782
783 Particle_t locPID;
784 const DMCThrown* locMCThrown;
785
786 vector<const DBeamPhoton*> locMCGENBeamPhotons;
787 locEventLoop->Get(locMCGENBeamPhotons, "MCGEN");
788
789 vector<const DBeamPhoton*> locBeamPhotons;
790 locEventLoop->Get(locBeamPhotons, "TRUTH");
791
792 //FILL HISTOGRAMS
793 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
794 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
795 Lock_Action();
796 {
797 for(size_t loc_i = 0; loc_i < locMCGENBeamPhotons.size(); ++loc_i)
798 {
799 dMCGENBeamParticle_P->Fill(locMCGENBeamPhotons[loc_i]->energy());
800 dMCGENBeamParticle_Time->Fill(locMCGENBeamPhotons[loc_i]->time());
801 }
802 for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
803 {
804 dAllBeamParticle_P->Fill(locBeamPhotons[loc_i]->energy());
805 dAllBeamParticle_Time->Fill(locBeamPhotons[loc_i]->time());
806 }
807 }
808 Unlock_Action();
809
810 for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
811 {
812 locMCThrown = locMCThrowns[loc_i];
813 locPID = (Particle_t)locMCThrown->type;
814 if(dHistMap_P.find(locPID) == dHistMap_P.end())
815 continue; //not interested in histogramming
816
817 DVector3 locMomentum = locMCThrown->momentum();
818 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
819 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
820 double locP = locMomentum.Mag();
821
822 //FILL HISTOGRAMS
823 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
824 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
825 Lock_Action();
826 {
827 dHistMap_P[locPID]->Fill(locP);
828 dHistMap_Phi[locPID]->Fill(locPhi);
829 dHistMap_Theta[locPID]->Fill(locTheta);
830 dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
831 dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
832 dHistMap_VertexZ[locPID]->Fill(locMCThrown->position().Z());
833 dHistMap_VertexYVsX[locPID]->Fill(locMCThrown->position().X(), locMCThrown->position().Y());
834 dHistMap_VertexT[locPID]->Fill(locMCThrown->time());
835 }
836 Unlock_Action();
837 }
838 return true;
839}
840
841void DHistogramAction_ReconnedThrownKinematics::Initialize(JEventLoop* locEventLoop)
842{
843 string locHistName, locHistTitle, locParticleName, locParticleROOTName;
844 Particle_t locPID;
845
846 Run_Update(locEventLoop);
847
848 //CREATE THE HISTOGRAMS
849 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
850 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
851 {
852
853 CreateAndChangeTo_ActionDirectory();
854
855 // Beam Particle
856 locPID = Gamma;
857 locParticleName = string("Beam_") + ParticleType(locPID);
858 locParticleROOTName = ParticleName_ROOT(locPID);
859 CreateAndChangeTo_Directory(locParticleName, locParticleName);
860
861 locHistName = "Momentum";
862 locHistTitle = string("Thrown Beam ") + locParticleROOTName + string(";p (GeV/c)");
863 dBeamParticle_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
864 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
865
866 //PID
867 CreateAndChangeTo_Directory("PID", "PID");
868 {
869 //beta vs p
870 locHistName = "BetaVsP_Q+";
871 locHistTitle = "q^{+};p (GeV/c);#beta";
872 dHistMap_QBetaVsP[1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
873
874 locHistName = "BetaVsP_Q-";
875 locHistTitle = "q^{-};p (GeV/c);#beta";
876 dHistMap_QBetaVsP[-1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
877 }
878 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
879
880 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
881 {
882 locPID = dFinalStatePIDs[loc_i];
883 locParticleName = ParticleType(locPID);
884 locParticleROOTName = ParticleName_ROOT(locPID);
885 CreateAndChangeTo_Directory(locParticleName, locParticleName);
886
887 // Momentum
888 locHistName = "Momentum";
889 locHistTitle = string("Thrown ") + locParticleROOTName + string(";p (GeV/c)");
890 dHistMap_P[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
891
892 // Theta
893 locHistName = "Theta";
894 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#theta#circ");
895 dHistMap_Theta[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
896
897 // Phi
898 locHistName = "Phi";
899 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#phi#circ");
900 dHistMap_Phi[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
901
902 // P Vs Theta
903 locHistName = "PVsTheta";
904 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#theta#circ;p (GeV/c)");
905 dHistMap_PVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
906
907 // Phi Vs Theta
908 locHistName = "PhiVsTheta";
909 locHistTitle = string("Thrown ") + locParticleROOTName + string(";#theta#circ;#phi#circ");
910 dHistMap_PhiVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
911
912 // Vertex-Z
913 locHistName = "VertexZ";
914 locHistTitle = string("Thrown ") + locParticleROOTName + string(";Vertex-Z (cm)");
915 dHistMap_VertexZ[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
916
917 // Vertex-Y Vs Vertex-X
918 locHistName = "VertexYVsX";
919 locHistTitle = string("Thrown ") + locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
920 dHistMap_VertexYVsX[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
921
922 // Vertex-T
923 locHistName = "VertexT";
924 locHistTitle = string("Thrown ") + locParticleROOTName + string(";Vertex-T (ns)");
925 dHistMap_VertexT[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
926
927 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
928 }
929
930 //Return to the base directory
931 ChangeTo_BaseDirectory();
932 }
933 japp->RootUnLock(); //RELEASE ROOT LOCK!!
934}
935
936void DHistogramAction_ReconnedThrownKinematics::Run_Update(JEventLoop* locEventLoop)
937{
938 const DAnalysisUtilities* locAnalysisUtilities = NULL__null;
939 locEventLoop->GetSingle(locAnalysisUtilities);
940 dAnalysisUtilities = locAnalysisUtilities;
941}
942
943bool DHistogramAction_ReconnedThrownKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
944{
945 vector<const DMCThrown*> locMCThrowns, locMCThrowns_Decaying;
946 locEventLoop->Get(locMCThrowns, "FinalState");
947 locEventLoop->Get(locMCThrowns_Decaying, "Decaying");
948 locMCThrowns.insert(locMCThrowns.begin(), locMCThrowns_Decaying.begin(), locMCThrowns_Decaying.end());
949 if(locMCThrowns.empty())
950 return true; //e.g. non-simulated event
951
952 if(Get_CalledPriorWithComboFlag())
953 return true; //else double-counting!
954
955 const DMCThrownMatching* locMCThrownMatching = NULL__null;
956 locEventLoop->GetSingle(locMCThrownMatching);
957
958 vector<const DBeamPhoton*> locBeamPhotons;
959 locEventLoop->Get(locBeamPhotons);
960
961 //FILL HISTOGRAMS
962 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
963 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
964 Lock_Action();
965 {
966 for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
967 dBeamParticle_P->Fill(locBeamPhotons[loc_i]->energy());
968 }
969 Unlock_Action();
970
971 for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
972 {
973 const DMCThrown* locMCThrown = locMCThrowns[loc_i];
974 Particle_t locPID = (Particle_t)locMCThrown->type;
975
976 double locMatchFOM = 0.0;
977 double locBeta_Timing = 0.0;
978 if(ParticleCharge(locPID) != 0)
979 {
980 const DChargedTrackHypothesis* locChargedTrackHypothesis = locMCThrownMatching->Get_MatchingChargedHypothesis(locMCThrown, locMatchFOM);
981 if(locChargedTrackHypothesis == NULL__null)
982 continue; //not reconstructed
983 locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
984 }
985 else
986 {
987 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locMCThrownMatching->Get_MatchingNeutralHypothesis(locMCThrown, locMatchFOM);
988 if(locNeutralParticleHypothesis == NULL__null)
989 continue; //not reconstructed
990 locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
991 }
992
993 DVector3 locMomentum = locMCThrown->momentum();
994 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
995 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
996 double locP = locMomentum.Mag();
997 int locCharge = ParticleCharge(locPID);
998
999 //FILL HISTOGRAMS
1000 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1001 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1002 Lock_Action();
1003 {
1004 if(dHistMap_QBetaVsP.find(locCharge) != dHistMap_QBetaVsP.end())
1005 dHistMap_QBetaVsP[locCharge]->Fill(locP, locBeta_Timing);
1006 if(dHistMap_P.find(locPID) != dHistMap_P.end())
1007 {
1008 dHistMap_P[locPID]->Fill(locP);
1009 dHistMap_Phi[locPID]->Fill(locPhi);
1010 dHistMap_Theta[locPID]->Fill(locTheta);
1011 dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
1012 dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
1013 dHistMap_VertexZ[locPID]->Fill(locMCThrown->position().Z());
1014 dHistMap_VertexYVsX[locPID]->Fill(locMCThrown->position().X(), locMCThrown->position().Y());
1015 dHistMap_VertexT[locPID]->Fill(locMCThrown->time());
1016 }
1017 }
1018 Unlock_Action();
1019 }
1020 return true;
1021}
1022
1023void DHistogramAction_GenReconTrackComparison::Initialize(JEventLoop* locEventLoop)
1024{
1025 string locHistName, locHistTitle, locParticleName, locParticleROOTName;
1026 Particle_t locPID;
1027
1028 Run_Update(locEventLoop);
1029
1030 //CREATE THE HISTOGRAMS
1031 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1032 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1033 {
1034 CreateAndChangeTo_ActionDirectory();
1035
1036 locHistName = "DeltaT_RFBeamBunch";
1037 dRFBeamBunchDeltaT_Hist = GetOrCreate_Histogram<TH1I>(locHistName, ";RF #Deltat (Reconstructed - Thrown)", dNumRFDeltaTBins, dMinRFDeltaT, dMaxRFDeltaT);
1038
1039 deque<string> locPullNames(8, "");
1040 locPullNames[0] = "E"; locPullNames[1] = "Px"; locPullNames[2] = "Py"; locPullNames[3] = "Pz";
1041 locPullNames[4] = "Xx"; locPullNames[5] = "Xy"; locPullNames[6] = "Xz"; locPullNames[7] = "T";
1042
1043 deque<string> locPullTitles(8, "");
1044 locPullTitles[0] = "E"; locPullTitles[1] = "p_{x}"; locPullTitles[2] = "p_{y}"; locPullTitles[3] = "p_{z}";
1045 locPullTitles[4] = "x_{x}"; locPullTitles[5] = "x_{y}"; locPullTitles[6] = "x_{z}"; locPullTitles[7] = "t";
1046
1047 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
1048 {
1049 locPID = dFinalStatePIDs[loc_i];
1050 locParticleName = ParticleType(locPID);
1051 locParticleROOTName = ParticleName_ROOT(locPID);
1052 CreateAndChangeTo_Directory(locParticleName, locParticleName);
1053
1054 // MatchChiSqPerDF
1055 locHistName = string("MatchFOM");
1056 locHistTitle = locParticleROOTName + string(";Thrown/Reconstructed Matching FOM");
1057 dHistMap_MatchFOM[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMCMatchingFOMBins, 0.0, 1.0);
1058
1059 // DeltaP/P
1060 locHistName = string("DeltaPOverP");
1061 locHistTitle = locParticleROOTName + string(";#Deltap/p (Reconstructed - Thrown)");
1062 dHistMap_DeltaPOverP[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
1063
1064 // DeltaTheta
1065 locHistName = string("DeltaTheta");
1066 locHistTitle = locParticleROOTName + string(";#Delta#theta#circ (Reconstructed - Thrown)");
1067 dHistMap_DeltaTheta[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);
1068
1069 // DeltaPhi
1070 locHistName = string("DeltaPhi");
1071 locHistTitle = locParticleROOTName + string(";#Delta#phi#circ (Reconstructed - Thrown)");
1072 dHistMap_DeltaPhi[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
1073
1074 // DeltaT
1075 locHistName = string("DeltaT");
1076 locHistTitle = locParticleROOTName + string(";#Deltat (ns) (Reconstructed - Thrown)");
1077 dHistMap_DeltaT[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1078
1079 // DeltaT - BCAL
1080 locHistName = string("DeltaT_BCAL");
1081 locHistTitle = locParticleROOTName + string(" in BCAL;#Deltat (ns) (Reconstructed - Thrown)");
1082 dHistMap_DeltaT_BCAL[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1083
1084 // DeltaT - TOF (charged only)
1085 if(ParticleCharge(locPID) != 0)
1086 {
1087 locHistName = string("DeltaT_TOF");
1088 locHistTitle = locParticleROOTName + string(" in TOF;#Deltat (ns) (Reconstructed - Thrown)");
1089 dHistMap_DeltaT_TOF[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1090 }
1091
1092 // DeltaT - FCAL (neutral only)
1093 if(ParticleCharge(locPID) == 0)
1094 {
1095 locHistName = string("DeltaT_FCAL");
1096 locHistTitle = locParticleROOTName + string(" in FCAL;#Deltat (ns) (Reconstructed - Thrown)");
1097 dHistMap_DeltaT_FCAL[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1098 }
1099
1100 // DeltaVertexZ
1101 locHistName = string("DeltaVertexZ");
1102 locHistTitle = locParticleROOTName + string(";#DeltaVertex-Z (cm) (Reconstructed - Thrown)");
1103 dHistMap_DeltaVertexZ[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexZBins, dMinDeltaVertexZ, dMaxDeltaVertexZ);
1104
1105 // DeltaP/P Vs P
1106 locHistName = string("DeltaPOverPVsP");
1107 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltap/p (Reconstructed - Thrown)");
1108 dHistMap_DeltaPOverPVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
1109
1110 // DeltaP/P Vs Theta
1111 locHistName = string("DeltaPOverPVsTheta");
1112 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltap/p (Reconstructed - Thrown)");
1113 dHistMap_DeltaPOverPVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaPOverPBins, dMinDeltaPOverP, dMaxDeltaPOverP);
1114
1115 // DeltaTheta Vs P
1116 locHistName = string("DeltaThetaVsP");
1117 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#theta#circ (Reconstructed - Thrown)");
1118 dHistMap_DeltaThetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);
1119
1120 // DeltaTheta Vs Theta
1121 locHistName = string("DeltaThetaVsTheta");
1122 locHistTitle = locParticleROOTName + string(";#theta#circ;#Delta#theta#circ (Reconstructed - Thrown)");
1123 dHistMap_DeltaThetaVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaThetaBins, dMinDeltaTheta, dMaxDeltaTheta);
1124
1125 // DeltaPhi Vs P
1126 locHistName = string("DeltaPhiVsP");
1127 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#phi#circ (Reconstructed - Thrown)");
1128 dHistMap_DeltaPhiVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
1129
1130 // DeltaPhi Vs Theta
1131 locHistName = string("DeltaPhiVsTheta");
1132 locHistTitle = locParticleROOTName + string(";#theta#circ;#Delta#phi#circ (Reconstructed - Thrown)");
1133 dHistMap_DeltaPhiVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
1134
1135 // DeltaT Vs Theta
1136 locHistName = string("DeltaTVsTheta");
1137 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat (ns) (Reconstructed - Thrown)");
1138 dHistMap_DeltaTVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1139
1140 // DeltaT Vs P
1141 locHistName = string("DeltaTVsP");
1142 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat (ns) (Reconstructed - Thrown)");
1143 dHistMap_DeltaTVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1144
1145 // DeltaVertexZ Vs Theta
1146 locHistName = string("DeltaVertexZVsTheta");
1147 locHistTitle = locParticleROOTName + string(";#theta#circ;#DeltaVertex-Z (cm) (Reconstructed - Thrown)");
1148 dHistMap_DeltaVertexZVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumDeltaVertexZBins, dMinDeltaVertexZ, dMaxDeltaVertexZ);
1149
1150 // P Vs Theta
1151 locHistName = "PVsTheta_LargeDeltaT";
1152 locHistTitle = locParticleROOTName + string(";#theta#circ;p (GeV/c)");
1153 dHistMap_PVsTheta_LargeDeltaT[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
1154
1155 /************************************************************************ Pulls ************************************************************************/
1156
1157 CreateAndChangeTo_Directory("Pulls", "Pulls");
1158 for(size_t loc_j = 0; loc_j < dPullTypes.size(); ++loc_j)
1159 {
1160 if((ParticleCharge(locPID) != 0) && (dPullTypes[loc_j] == d_EPull))
1161 continue;
1162 if((ParticleCharge(locPID) == 0) && ((dPullTypes[loc_j] >= d_PxPull) && (dPullTypes[loc_j] <= d_PzPull)))
1163 continue;
1164
1165 //Pull 1D
1166 locHistName = locPullNames[loc_j] + string("Pull");
1167 locHistTitle = locParticleROOTName + string(";#Delta") + locPullTitles[loc_j] + string("/#sigma_{") + locPullTitles[loc_j] + string("} (Reconstructed - Thrown)");
1168 dHistMap_Pulls[locPID][dPullTypes[loc_j]] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
1169
1170 //Pull vs P
1171 locHistName = locPullNames[loc_j] + string("PullVsP");
1172 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta") + locPullTitles[loc_j] + string("/#sigma_{") + locPullTitles[loc_j] + string("} (Reconstructed - Thrown)");
1173 dHistMap_PullsVsP[locPID][dPullTypes[loc_j]] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
1174
1175 //Pull vs Theta
1176 locHistName = locPullNames[loc_j] + string("PullVsTheta");
1177 locHistTitle = locParticleROOTName + string(";#theta#circ;#Delta") + locPullTitles[loc_j] + string("/#sigma_{") + locPullTitles[loc_j] + string("} (Reconstructed - Thrown)");
1178 dHistMap_PullsVsTheta[locPID][dPullTypes[loc_j]] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
1179 }
1180
1181 //Delta-t Pulls - CDC & ST
1182 if(ParticleCharge(locPID) != 0)
1183 {
1184 //CDC
1185 locHistName = "TimePull_CDC";
1186 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
1187 dHistMap_TimePull_CDC[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
1188
1189 locHistName = "TimePullVsTheta_CDC";
1190 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat/#sigma_{#Deltat}");
1191 dHistMap_TimePullVsTheta_CDC[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
1192
1193 locHistName = "TimePullVsP_CDC";
1194 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
1195 dHistMap_TimePullVsP_CDC[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
1196
1197 //ST
1198 locHistName = "TimePull_ST";
1199 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
1200 dHistMap_TimePull_ST[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
1201
1202 locHistName = "TimePullVsTheta_ST";
1203 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat/#sigma_{#Deltat}");
1204 dHistMap_TimePullVsTheta_ST[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
1205
1206 locHistName = "TimePullVsP_ST";
1207 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
1208 dHistMap_TimePullVsP_ST[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
1209 }
1210
1211 //Delta-t Pulls - BCAL
1212 locHistName = "TimePull_BCAL";
1213 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
1214 dHistMap_TimePull_BCAL[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
1215
1216 locHistName = "TimePullVsTheta_BCAL";
1217 locHistTitle = locParticleROOTName + string(";#theta#circ;#Deltat/#sigma_{#Deltat}");
1218 dHistMap_TimePullVsTheta_BCAL[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPullBins, -10.0, 10.0);
1219
1220 locHistName = "TimePullVsP_BCAL";
1221 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
1222 dHistMap_TimePullVsP_BCAL[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
1223
1224 //Delta-t Pulls - TOF
1225 if(ParticleCharge(locPID) != 0) //TOF
1226 {
1227 locHistName = "TimePull_TOF";
1228 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
1229 dHistMap_TimePull_TOF[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
1230
1231 locHistName = "TimePullVsP_TOF";
1232 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
1233 dHistMap_TimePullVsP_TOF[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
1234 }
1235
1236 //Delta-t Pulls - FCAL
1237 locHistName = "TimePull_FCAL";
1238 locHistTitle = locParticleROOTName + string(";#Deltat/#sigma_{#Deltat}");
1239 dHistMap_TimePull_FCAL[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, -10.0, 10.0);
1240
1241 locHistName = "TimePullVsP_FCAL";
1242 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
1243 dHistMap_TimePullVsP_FCAL[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, -10.0, 10.0);
1244
1245 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
1246
1247 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
1248 }
1249
1250 //Return to the base directory
1251 ChangeTo_BaseDirectory();
1252 }
1253 japp->RootUnLock(); //RELEASE ROOT LOCK!!
1254}
1255
1256void DHistogramAction_GenReconTrackComparison::Run_Update(JEventLoop* locEventLoop)
1257{
1258 DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
1259 DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
1260 locGeometry->GetTargetZ(dTargetZCenter);
1261}
1262
1263bool DHistogramAction_GenReconTrackComparison::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1264{
1265 vector<const DMCThrown*> locMCThrowns;
1266 locEventLoop->Get(locMCThrowns);
1267 if(locMCThrowns.empty())
1268 return true; //e.g. non-simulated event
1269
1270 if(Get_CalledPriorWithComboFlag())
1271 return true; //else double-counting!
1272
1273 Particle_t locPID;
1274 double locDeltaPOverP, locDeltaTheta, locDeltaPhi, locDeltaVertexZ;
1275 double locThrownP, locThrownTheta, locDeltaT;
1276
1277 vector<const DMCThrownMatching*> locMCThrownMatchingVector;
1278 locEventLoop->Get(locMCThrownMatchingVector);
1279 if(locMCThrownMatchingVector.empty())
1280 return true;
1281 const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0];
1282
1283 const DEventRFBunch* locThrownEventRFBunch = NULL__null;
1284 locEventLoop->GetSingle(locThrownEventRFBunch, "Thrown");
1285
1286 //RF time difference
1287 vector<const DEventRFBunch*> locEventRFBunches;
1288 locEventLoop->Get(locEventRFBunches);
1289 const DEventRFBunch* locEventRFBunch = locEventRFBunches[0];
1290 double locRFTime = locEventRFBunch->dTime;
1291 double locRFDeltaT = locRFTime - locThrownEventRFBunch->dTime;
1292
1293 //FILL HISTOGRAMS
1294 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1295 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1296 Lock_Action();
1297 {
1298 dRFBeamBunchDeltaT_Hist->Fill(locRFDeltaT);
1299 }
1300 Unlock_Action();
1301
1302 //charged particles
1303 map<const DMCThrown*, pair<const DChargedTrack*, double> > locThrownToChargedMap;
1304 locMCThrownMatching->Get_ThrownToChargedMap(locThrownToChargedMap);
1305 map<const DMCThrown*, pair<const DChargedTrack*, double> >::iterator locChargedIterator = locThrownToChargedMap.begin();
1306 for(; locChargedIterator != locThrownToChargedMap.end(); ++locChargedIterator)
1307 {
1308 const DMCThrown* locMCThrown = locChargedIterator->first;
1309 locPID = (Particle_t)locMCThrown->type;
1310 if(dHistMap_DeltaPOverP.find(locPID) == dHistMap_DeltaPOverP.end())
1311 continue; //e.g. not interested in histogramming
1312
1313 double locMatchFOM = locChargedIterator->second.second;
1314 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedIterator->second.first->Get_Hypothesis(locPID);
1315 if(locChargedTrackHypothesis == NULL__null)
1316 locChargedTrackHypothesis = locChargedIterator->second.first->Get_BestFOM();
1317
1318 locThrownP = locMCThrown->momentum().Mag();
1319 locThrownTheta = locMCThrown->momentum().Theta()*180.0/TMath::Pi();
1320 locDeltaPOverP = (locChargedTrackHypothesis->momentum().Mag() - locThrownP)/locThrownP;
1321 locDeltaTheta = locChargedTrackHypothesis->momentum().Theta()*180.0/TMath::Pi() - locThrownTheta;
1322 locDeltaPhi = locChargedTrackHypothesis->momentum().Phi()*180.0/TMath::Pi() - locMCThrown->momentum().Phi()*180.0/TMath::Pi();
1323 locDeltaT = locChargedTrackHypothesis->time() - locMCThrown->time(); //time comparison isn't fair if track comes from a detached vertex!!!
1324 locDeltaVertexZ = locChargedTrackHypothesis->position().Z() - locMCThrown->position().Z();
1325 const TMatrixDSym& locCovarianceMatrix = *(locChargedTrackHypothesis->errorMatrix().get());
1326
1327 const DTrackTimeBased* locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
1328
1329 double locStartTime = locThrownEventRFBunch->dTime + (locMCThrown->z() - dTargetZCenter)/29.9792458;
1330 double locTimePull = (locStartTime - locChargedTrackHypothesis->time())/sqrt(locCovarianceMatrix(6, 6));
1331 double locT0Pull = (locStartTime - locChargedTrackHypothesis->t0())/locChargedTrackHypothesis->t0_err();
1332
1333 //FILL HISTOGRAMS
1334 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1335 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1336 Lock_Action();
1337 {
1338 dHistMap_MatchFOM[locPID]->Fill(locMatchFOM);
1339 dHistMap_DeltaPOverP[locPID]->Fill(locDeltaPOverP);
1340 dHistMap_DeltaTheta[locPID]->Fill(locDeltaTheta);
1341 dHistMap_DeltaPhi[locPID]->Fill(locDeltaPhi);
1342 dHistMap_DeltaT[locPID]->Fill(locDeltaT);
1343 if(locChargedTrackHypothesis->t0_detector() == SYS_START)
1344 {
1345 dHistMap_TimePull_ST[locPID]->Fill(locT0Pull);
1346 dHistMap_TimePullVsTheta_ST[locPID]->Fill(locThrownTheta, locT0Pull);
1347 dHistMap_TimePullVsP_ST[locPID]->Fill(locThrownP, locT0Pull);
1348 }
1349 if(locChargedTrackHypothesis->t0_detector() == SYS_CDC)
1350 {
1351 dHistMap_TimePull_CDC[locPID]->Fill(locT0Pull);
1352 dHistMap_TimePullVsTheta_CDC[locPID]->Fill(locThrownTheta, locT0Pull);
1353 dHistMap_TimePullVsP_CDC[locPID]->Fill(locThrownP, locT0Pull);
1354 }
1355 else if(locChargedTrackHypothesis->t1_detector() == SYS_CDC)
1356 {
1357 dHistMap_TimePull_CDC[locPID]->Fill(locTimePull);
1358 dHistMap_TimePullVsTheta_CDC[locPID]->Fill(locThrownTheta, locTimePull);
1359 dHistMap_TimePullVsP_CDC[locPID]->Fill(locThrownP, locTimePull);
1360 }
1361 if(locChargedTrackHypothesis->t1_detector() == SYS_BCAL)
1362 {
1363 dHistMap_DeltaT_BCAL[locPID]->Fill(locDeltaT);
1364 dHistMap_TimePull_BCAL[locPID]->Fill(locTimePull);
1365 dHistMap_TimePullVsTheta_BCAL[locPID]->Fill(locThrownTheta, locTimePull);
1366 dHistMap_TimePullVsP_BCAL[locPID]->Fill(locThrownP, locTimePull);
1367 }
1368 else if(locChargedTrackHypothesis->t1_detector() == SYS_TOF)
1369 {
1370 dHistMap_DeltaT_TOF[locPID]->Fill(locDeltaT);
1371 dHistMap_TimePull_TOF[locPID]->Fill(locTimePull);
1372 dHistMap_TimePullVsP_TOF[locPID]->Fill(locThrownP, locTimePull);
1373 }
1374 else if(locChargedTrackHypothesis->t1_detector() == SYS_FCAL)
1375 {
1376 dHistMap_TimePull_FCAL[locPID]->Fill(locTimePull);
1377 dHistMap_TimePullVsP_FCAL[locPID]->Fill(locThrownP, locTimePull);
1378 }
1379 dHistMap_DeltaVertexZ[locPID]->Fill(locDeltaVertexZ);
1380 dHistMap_DeltaPOverPVsP[locPID]->Fill(locThrownP, locDeltaPOverP);
1381 dHistMap_DeltaPOverPVsTheta[locPID]->Fill(locThrownTheta, locDeltaPOverP);
1382 dHistMap_DeltaThetaVsP[locPID]->Fill(locThrownP, locDeltaTheta);
1383 dHistMap_DeltaThetaVsTheta[locPID]->Fill(locThrownTheta, locDeltaTheta);
1384 dHistMap_DeltaPhiVsP[locPID]->Fill(locThrownP, locDeltaPhi);
1385 dHistMap_DeltaPhiVsTheta[locPID]->Fill(locThrownTheta, locDeltaPhi);
1386 dHistMap_DeltaTVsTheta[locPID]->Fill(locThrownTheta, locDeltaT);
1387 dHistMap_DeltaTVsP[locPID]->Fill(locThrownP, locDeltaT);
1388 dHistMap_DeltaVertexZVsTheta[locPID]->Fill(locThrownTheta, locDeltaVertexZ);
1389 if((locTrackTimeBased->FOM > 0.01) && (locDeltaT >= 1.002))
1390 dHistMap_PVsTheta_LargeDeltaT[locPID]->Fill(locThrownTheta, locThrownP);
1391
1392 for(size_t loc_j = 0; loc_j < dPullTypes.size(); ++loc_j)
1393 {
1394 if(dPullTypes[loc_j] == d_EPull)
1395 continue;
1396 double locPull = 0.0;
1397 if((dPullTypes[loc_j] >= d_PxPull) && (dPullTypes[loc_j] <= d_PzPull))
1398 {
1399 int locIndex = int(dPullTypes[loc_j] - d_PxPull);
1400 locPull = (locChargedTrackHypothesis->momentum()(locIndex) - locMCThrown->momentum()(locIndex))/sqrt(locCovarianceMatrix(locIndex, locIndex));
1401 }
1402 else if((dPullTypes[loc_j] >= d_XxPull) && (dPullTypes[loc_j] <= d_XzPull))
1403 {
1404 int locIndex = int(dPullTypes[loc_j] - d_XxPull);
1405 locPull = (locChargedTrackHypothesis->position()(locIndex) - locMCThrown->position()(locIndex))/sqrt(locCovarianceMatrix(locIndex + 3, locIndex + 3));
1406 }
1407 else if(dPullTypes[loc_j] == d_TPull)
1408 locPull = (locChargedTrackHypothesis->time() - locMCThrown->time())/sqrt(locCovarianceMatrix(6, 6));
1409 (dHistMap_Pulls[locPID])[dPullTypes[loc_j]]->Fill(locPull);
1410 (dHistMap_PullsVsP[locPID])[dPullTypes[loc_j]]->Fill(locThrownP, locPull);
1411 (dHistMap_PullsVsTheta[locPID])[dPullTypes[loc_j]]->Fill(locThrownTheta, locPull);
1412 }
1413 }
1414 Unlock_Action();
1415 }
1416
1417 //neutral particles
1418 map<const DMCThrown*, pair<const DNeutralParticle*, double> > locThrownToNeutralMap;
1419 locMCThrownMatching->Get_ThrownToNeutralMap(locThrownToNeutralMap);
1420 map<const DMCThrown*, pair<const DNeutralParticle*, double> >::iterator locNeutralIterator = locThrownToNeutralMap.begin();
1421 for(; locNeutralIterator != locThrownToNeutralMap.end(); ++locNeutralIterator)
1422 {
1423 const DMCThrown* locMCThrown = locNeutralIterator->first;
1424 locPID = (Particle_t)locMCThrown->type;
1425 if(dHistMap_DeltaPOverP.find(locPID) == dHistMap_DeltaPOverP.end())
1426 continue; //e.g. not interested in histogramming
1427
1428 double locMatchFOM = locNeutralIterator->second.second;
1429 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralIterator->second.first->Get_Hypothesis(locPID);
1430 if(locNeutralParticleHypothesis == NULL__null)
1431 locNeutralParticleHypothesis = locNeutralIterator->second.first->Get_BestFOM();
1432
1433 const DNeutralShower* locNeutralShower = locNeutralParticleHypothesis->Get_NeutralShower();
1434
1435 locThrownP = locMCThrown->momentum().Mag();
1436 locThrownTheta = locMCThrown->momentum().Theta()*180.0/TMath::Pi();
1437 locDeltaPOverP = (locNeutralParticleHypothesis->momentum().Mag() - locThrownP)/locThrownP;
1438 locDeltaTheta = locNeutralParticleHypothesis->momentum().Theta()*180.0/TMath::Pi() - locThrownTheta;
1439 locDeltaPhi = locNeutralParticleHypothesis->momentum().Phi()*180.0/TMath::Pi() - locMCThrown->momentum().Phi()*180.0/TMath::Pi();
1440 locDeltaT = locNeutralParticleHypothesis->time() - locMCThrown->time(); //time comparison isn't fair if track comes from a detached vertex!!!
1441 locDeltaVertexZ = locNeutralParticleHypothesis->position().Z() - locMCThrown->position().Z();
1442 const TMatrixDSym& locCovarianceMatrix = *(locNeutralParticleHypothesis->errorMatrix().get());
1443
1444 double locStartTime = locThrownEventRFBunch->dTime + (locMCThrown->z() - dTargetZCenter)/29.9792458;
1445 double locTimePull = (locStartTime - locNeutralParticleHypothesis->time())/sqrt(locCovarianceMatrix(6, 6));
1446
1447 //FILL HISTOGRAMS
1448 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1449 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1450 Lock_Action();
1451 {
1452 dHistMap_MatchFOM[locPID]->Fill(locMatchFOM);
1453 dHistMap_DeltaPOverP[locPID]->Fill(locDeltaPOverP);
1454 dHistMap_DeltaTheta[locPID]->Fill(locDeltaTheta);
1455 dHistMap_DeltaPhi[locPID]->Fill(locDeltaPhi);
1456 dHistMap_DeltaT[locPID]->Fill(locDeltaT);
1457 if(locNeutralParticleHypothesis->t1_detector() == SYS_BCAL)
1458 {
1459 dHistMap_DeltaT_BCAL[locPID]->Fill(locDeltaT);
1460 dHistMap_TimePull_BCAL[locPID]->Fill(locTimePull);
1461 dHistMap_TimePullVsTheta_BCAL[locPID]->Fill(locThrownTheta, locTimePull);
1462 dHistMap_TimePullVsP_BCAL[locPID]->Fill(locThrownP, locTimePull);
1463 }
1464 else if(locNeutralParticleHypothesis->t1_detector() == SYS_FCAL)
1465 {
1466 dHistMap_DeltaT_FCAL[locPID]->Fill(locDeltaT);
1467 dHistMap_TimePull_FCAL[locPID]->Fill(locTimePull);
1468 dHistMap_TimePullVsP_FCAL[locPID]->Fill(locThrownP, locTimePull);
1469 }
1470
1471 dHistMap_DeltaVertexZ[locPID]->Fill(locDeltaVertexZ);
1472 dHistMap_DeltaPOverPVsP[locPID]->Fill(locThrownP, locDeltaPOverP);
1473 dHistMap_DeltaPOverPVsTheta[locPID]->Fill(locThrownTheta, locDeltaPOverP);
1474 dHistMap_DeltaThetaVsP[locPID]->Fill(locThrownP, locDeltaTheta);
1475 dHistMap_DeltaThetaVsTheta[locPID]->Fill(locThrownTheta, locDeltaTheta);
1476 dHistMap_DeltaPhiVsP[locPID]->Fill(locThrownP, locDeltaPhi);
1477 dHistMap_DeltaPhiVsTheta[locPID]->Fill(locThrownTheta, locDeltaPhi);
1478 dHistMap_DeltaTVsTheta[locPID]->Fill(locThrownTheta, locDeltaT);
1479 dHistMap_DeltaTVsP[locPID]->Fill(locThrownP, locDeltaT);
1480 dHistMap_DeltaVertexZVsTheta[locPID]->Fill(locThrownTheta, locDeltaVertexZ);
1481 if(locDeltaT >= 1.002)
1482 dHistMap_PVsTheta_LargeDeltaT[locPID]->Fill(locThrownTheta, locThrownP);
1483
1484 for(size_t loc_j = 0; loc_j < dPullTypes.size(); ++loc_j)
1485 {
1486 if((dPullTypes[loc_j] >= d_PxPull) && (dPullTypes[loc_j] <= d_PzPull))
1487 continue;
1488 double locPull = 0.0;
1489 if(dPullTypes[loc_j] == d_EPull)
1490 locPull = (locNeutralShower->dEnergy - locMCThrown->energy())/sqrt((*(locNeutralShower->dCovarianceMatrix))(0, 0));
1491 else if((dPullTypes[loc_j] >= d_XxPull) && (dPullTypes[loc_j] <= d_XzPull))
1492 {
1493 int locIndex = int(dPullTypes[loc_j] - d_XxPull);
1494 locPull = (locNeutralParticleHypothesis->position()(locIndex) - locMCThrown->position()(locIndex))/sqrt(locCovarianceMatrix(locIndex + 3, locIndex + 3));
1495 }
1496 else if(dPullTypes[loc_j] == d_TPull)
1497 locPull = (locNeutralParticleHypothesis->time() - locMCThrown->time())/sqrt(locCovarianceMatrix(6, 6));
1498 (dHistMap_Pulls[locPID])[dPullTypes[loc_j]]->Fill(locPull);
1499 (dHistMap_PullsVsP[locPID])[dPullTypes[loc_j]]->Fill(locThrownP, locPull);
1500 (dHistMap_PullsVsTheta[locPID])[dPullTypes[loc_j]]->Fill(locThrownTheta, locPull);
1501 }
1502
1503 }
1504 Unlock_Action();
1505 }
1506 return true;
1507}
1508
1509void DHistogramAction_TruePID::Initialize(JEventLoop* locEventLoop)
1510{
1511 string locStepName, locStepROOTName, locHistTitle, locHistName, locParticleName, locParticleROOTName;
1512
1513 size_t locNumSteps = Get_Reaction()->Get_NumReactionSteps();
1514 dHistDeque_P_CorrectID.resize(locNumSteps);
1515 dHistDeque_P_IncorrectID.resize(locNumSteps);
1516 dHistDeque_PVsTheta_CorrectID.resize(locNumSteps);
1517 dHistDeque_PVsTheta_IncorrectID.resize(locNumSteps);
1518
1519 Run_Update(locEventLoop);
1520
1521 //CREATE THE HISTOGRAMS
1522 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1523 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1524 {
1525 CreateAndChangeTo_ActionDirectory();
1526
1527 for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)
1528 {
1529 auto locDetectedPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_AllCharges, false);
1530 if(locDetectedPIDs.empty())
1531 continue;
1532
1533 const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
1534 locStepName = DAnalysis::Get_StepName(locReactionStep, true, false);
1535 locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);
1536 CreateAndChangeTo_Directory(locStepName, locStepName);
1537
1538 for(auto locPID : locDetectedPIDs)
1539 {
1540 locParticleName = ParticleType(locPID);
1541 locParticleROOTName = ParticleName_ROOT(locPID);
1542
1543 if(dHistDeque_P_CorrectID[loc_i].find(locPID) != dHistDeque_P_CorrectID[loc_i].end())
1544 continue; //hists already created for this pid
1545
1546 //P of Correct ID
1547 locHistName = string("Momentum_CorrectID_") + locParticleName;
1548 locHistTitle = string("Correct ") + locParticleROOTName + string(" ID, ") + locStepROOTName + string(";p (GeV/c)");
1549 dHistDeque_P_CorrectID[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
1550
1551 //P of Incorrect ID
1552 locHistName = string("Momentum_IncorrectID_") + locParticleName;
1553 locHistTitle = string("Incorrect ") + locParticleROOTName + string(" ID, ") + locStepROOTName + string(";p (GeV/c)");
1554 dHistDeque_P_IncorrectID[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
1555
1556 //P Vs Theta of Correct ID
1557 locHistName = string("PVsTheta_CorrectID_") + locParticleName;
1558 locHistTitle = string("Correct ") + locParticleROOTName + string(" ID, ") + locStepROOTName + string(";#theta#circ;p (GeV/c)");
1559 dHistDeque_PVsTheta_CorrectID[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
1560
1561 //P Vs Theta of Incorrect ID
1562 locHistName = string("PVsTheta_IncorrectID_") + locParticleName;
1563 locHistTitle = string("Incorrect ") + locParticleROOTName + string(" ID, ") + locStepROOTName + string(";#theta#circ;p (GeV/c)");
1564 dHistDeque_PVsTheta_IncorrectID[loc_i][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
1565 }
1566 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");
1567 } //end of step loop
1568
1569 //# Combos Pass/Fail All True PID
1570 locHistName = "Combo_TruePIDStatus";
1571 locHistTitle = Get_Reaction()->Get_ReactionName() + string(";# Combos;All Combo Particles True PID Status");
1572 dHist_TruePIDStatus = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 2, -0.5, 1.5);
1573
1574 //Return to the base directory
1575 ChangeTo_BaseDirectory();
1576 }
1577 japp->RootUnLock(); //RELEASE ROOT LOCK!!
1578}
1579
1580void DHistogramAction_TruePID::Run_Update(JEventLoop* locEventLoop)
1581{
1582 const DAnalysisUtilities* locAnalysisUtilities = NULL__null;
1583 locEventLoop->GetSingle(locAnalysisUtilities);
1584 dAnalysisUtilities = locAnalysisUtilities;
1585}
1586
1587bool DHistogramAction_TruePID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1588{
1589 vector<const DMCThrownMatching*> locMCThrownMatchingVector;
1590 locEventLoop->Get(locMCThrownMatchingVector);
1591 if(locMCThrownMatchingVector.empty())
1592 return true;
1593 const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0];
1594 double locP, locTheta;
1595 const DMCThrown* locMCThrown;
1596 Particle_t locPID;
1597
1598 int locComboTruePIDStatus = 1;
1599 for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
1600 {
1601 const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);
1602 auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i));
1603 for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
1604 {
1605 locPID = locParticles[loc_j]->PID();
1606
1607 double locMatchFOM = 0.0;
1608 if(ParticleCharge(locPID) == 0)
1609 locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_j]), locMatchFOM);
1610 else
1611 locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(static_cast<const DChargedTrackHypothesis*>(locParticles[loc_j]), locMatchFOM);
1612
1613 bool locCutResult = (locMCThrown == NULL__null) ? false : (((Particle_t)locMCThrown->type) == locPID);
1614 if(!locCutResult)
1615 locComboTruePIDStatus = 0;
1616
1617 //check if duplicate
1618 const JObject* locSourceObject = Get_FinalParticle_SourceObject(locParticles[loc_j]);
1619 pair<Particle_t, const JObject*> locParticleInfo(locParticles[loc_j]->PID(), locSourceObject);
1620 pair<size_t, pair<Particle_t, const JObject*> > locHistInfo(loc_i, locParticleInfo);
1621 if(dPreviouslyHistogrammedParticles.find(locHistInfo) != dPreviouslyHistogrammedParticles.end())
1622 continue; //previously histogrammed
1623 dPreviouslyHistogrammedParticles.insert(locHistInfo);
1624
1625 locP = locParticles[loc_j]->momentum().Mag();
1626 locTheta = locParticles[loc_j]->momentum().Theta()*180.0/TMath::Pi();
1627
1628 //FILL HISTOGRAMS
1629 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1630 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1631 Lock_Action();
1632 {
1633 if(locCutResult)
1634 {
1635 dHistDeque_P_CorrectID[loc_i][locPID]->Fill(locP);
1636 dHistDeque_PVsTheta_CorrectID[loc_i][locPID]->Fill(locTheta, locP);
1637 }
1638 else
1639 {
1640 dHistDeque_P_IncorrectID[loc_i][locPID]->Fill(locP);
1641 dHistDeque_PVsTheta_IncorrectID[loc_i][locPID]->Fill(locTheta, locP);
1642 }
1643 }
1644 Unlock_Action();
1645 }
1646 }
1647
1648 //FILL HISTOGRAMS
1649 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1650 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1651 Lock_Action();
1652 {
1653 dHist_TruePIDStatus->Fill(locComboTruePIDStatus);
1654 }
1655 Unlock_Action();
1656
1657 return true;
1658}
1659