Bug Summary

File:libraries/ANALYSIS/DHistogramActions_Thrown.cc
Location:line 909, column 3
Description:Value stored to 'locEventRFBunch' is never read

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