Bug Summary

File:DHistogramActions_Thrown.cc
Location:line 585, column 82
Description:Called C++ object pointer is uninitialized

Annotated Source Code

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