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