Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DHistogramActions_Reaction.cc
Go to the documentation of this file.
2 
3 void DHistogramAction_PID::Initialize(JEventLoop* locEventLoop)
4 {
5  string locHistName, locHistTitle;
6  string locParticleName, locParticleName2, locParticleROOTName, locParticleROOTName2;
7  Particle_t locPID, locPID2;
8 
9  vector<const DParticleID*> locParticleIDs;
10  locEventLoop->Get(locParticleIDs);
11  auto locDesiredPIDs = Get_Reaction()->Get_FinalPIDs(-1, false, false, d_AllCharges, false);
12 
13  vector<const DMCThrown*> locMCThrowns;
14  locEventLoop->Get(locMCThrowns);
15 
16  vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;
17  locEventLoop->Get(locAnalysisUtilitiesVector);
18 
19  //CREATE THE HISTOGRAMS
20  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
21  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
22  {
24  dParticleID = locParticleIDs[0];
25  dAnalysisUtilities = locAnalysisUtilitiesVector[0];
26  for(size_t loc_i = 0; loc_i < locDesiredPIDs.size(); ++loc_i)
27  {
28  locPID = locDesiredPIDs[loc_i];
29  locParticleName = ParticleType(locPID);
30  locParticleROOTName = ParticleName_ROOT(locPID);
31  CreateAndChangeTo_Directory(locParticleName, locParticleName);
32 
33  if(ParticleCharge(locPID) == 0)
34  {
35  //BCAL
36  CreateAndChangeTo_Directory("BCAL", "BCAL");
37  locHistName = "Beta";
38  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;#beta");
39  dHistMap_Beta[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBetaBins, dMinBeta, dMaxBeta);
40  gDirectory->cd("..");
41 
42  //FCAL
43  CreateAndChangeTo_Directory("FCAL", "FCAL");
44  locHistName = "Beta";
45  locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;#beta");
46  dHistMap_Beta[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBetaBins, dMinBeta, dMaxBeta);
47  gDirectory->cd("..");
48  }
49 
50  //q = 0
51  if(locPID == Gamma)
52  {
53  //BCAL
54  CreateAndChangeTo_Directory("BCAL", "BCAL");
55 
56  locHistName = "BetaVsP";
57  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#beta");
58  dHistMap_BetaVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
59 
60  locHistName = "DeltaTVsShowerE_Photon";
61  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#Deltat_{BCAL - RF}");
62  dHistMap_DeltaTVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
63 
64  locHistName = "TimePullVsShowerE_Photon";
65  locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
66  dHistMap_TimePullVsP[Gamma][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
67 
68  locHistName = "TimeFOMVsShowerE_Photon";
69  locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
70  dHistMap_TimeFOMVsP[Gamma][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
71 
72  gDirectory->cd("..");
73 
74  //FCAL
75  CreateAndChangeTo_Directory("FCAL", "FCAL");
76 
77  locHistName = "BetaVsP";
78  locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#beta");
79  dHistMap_BetaVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
80 
81  locHistName = "DeltaTVsShowerE_Photon";
82  locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat_{FCAL - RF}");
83  dHistMap_DeltaTVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
84 
85  locHistName = "TimePullVsShowerE_Photon";
86  locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
87  dHistMap_TimePullVsP[Gamma][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
88 
89  locHistName = "TimeFOMVsShowerE_Photon";
90  locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
91  dHistMap_TimeFOMVsP[Gamma][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
92 
93  gDirectory->cd("..");
94  }
95  else if(ParticleCharge(locPID) != 0)
96  {
97  //SC
98  CreateAndChangeTo_Directory("SC", "SC");
99 
100  locHistName = "dEdXVsP";
101  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC dE/dX (MeV/cm)");
102  dHistMap_dEdXVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
103 
104  locHistName = "DeltadEdXVsP";
105  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC #Delta(dE/dX) (MeV/cm)");
106  dHistMap_DeltadEdXVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
107 
108  locHistName = "DeltaBetaVsP";
109  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC #Delta#beta");
110  dHistMap_DeltaBetaVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
111 
112  locHistName = "BetaVsP";
113  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC #beta");
114  dHistMap_BetaVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
115 
116  locHistName = "DeltaTVsP";
117  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{SC - RF}");
118  dHistMap_DeltaTVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
119 
120  locHistName = string("TimePullVsP_") + locParticleName;
121  locHistTitle = string("SC ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
122  dHistMap_TimePullVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
123 
124  locHistName = string("TimeFOMVsP_") + locParticleName;
125  locHistTitle = string("SC") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");
126  dHistMap_TimeFOMVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
127 
128  gDirectory->cd("..");
129 
130  //TOF
131  CreateAndChangeTo_Directory("TOF", "TOF");
132 
133  locHistName = "dEdXVsP";
134  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF dE/dX (MeV/cm)");
135  dHistMap_dEdXVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
136 
137  locHistName = "DeltadEdXVsP";
138  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF #Delta(dE/dX) (MeV/cm)");
139  dHistMap_DeltadEdXVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
140 
141  locHistName = "DeltaBetaVsP";
142  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF #Delta#beta");
143  dHistMap_DeltaBetaVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
144 
145  locHistName = "BetaVsP";
146  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF #beta");
147  dHistMap_BetaVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
148 
149  locHistName = "DeltaTVsP";
150  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{TOF - RF}");
151  dHistMap_DeltaTVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
152 
153  locHistName = string("TimePullVsP_") + locParticleName;
154  locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
155  dHistMap_TimePullVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
156 
157  locHistName = string("TimeFOMVsP_") + locParticleName;
158  locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");
159  dHistMap_TimeFOMVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
160 
161  gDirectory->cd("..");
162 
163  //BCAL
164  CreateAndChangeTo_Directory("BCAL", "BCAL");
165 
166  locHistName = "BetaVsP";
167  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);BCAL #beta");
168  dHistMap_BetaVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
169 
170  locHistName = "DeltaBetaVsP";
171  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);BCAL #Delta#beta");
172  dHistMap_DeltaBetaVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
173 
174  locHistName = "DeltaTVsP";
175  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{BCAL - RF}");
176  dHistMap_DeltaTVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
177 
178  locHistName = string("TimePullVsP_") + locParticleName;
179  locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
180  dHistMap_TimePullVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
181 
182  locHistName = string("TimeFOMVsP_") + locParticleName;
183  locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");
184  dHistMap_TimeFOMVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
185 
186  locHistName = "EOverPVsP";
187  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);BCAL E_{Shower}/p_{Track} (c);");
188  dHistMap_EOverPVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
189 
190  locHistName = "EOverPVsTheta";
191  locHistTitle = locParticleROOTName + string(" Candidates;#theta#circ;BCAL E_{Shower}/p_{Track} (c);");
192  dHistMap_EOverPVsTheta[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALThetaBins, dMinBCALTheta, dMaxBCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
193 
194  gDirectory->cd("..");
195 
196  //FCAL
197  CreateAndChangeTo_Directory("FCAL", "FCAL");
198 
199  locHistName = "BetaVsP";
200  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FCAL #beta");
201  dHistMap_BetaVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
202 
203  locHistName = "DeltaBetaVsP";
204  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FCAL #Delta#beta");
205  dHistMap_DeltaBetaVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
206 
207  locHistName = "DeltaTVsP";
208  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{FCAL - RF}");
209  dHistMap_DeltaTVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
210 
211  locHistName = string("TimePullVsP_") + locParticleName;
212  locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");
213  dHistMap_TimePullVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
214 
215  locHistName = string("TimeFOMVsP_") + locParticleName;
216  locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");
217  dHistMap_TimeFOMVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
218 
219  locHistName = "EOverPVsP";
220  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FCAL E_{Shower}/p_{Track} (c);");
221  dHistMap_EOverPVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
222 
223  locHistName = "EOverPVsTheta";
224  locHistTitle = locParticleROOTName + string(" Candidates;#theta#circ;FCAL E_{Shower}/p_{Track} (c);");
225  dHistMap_EOverPVsTheta[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DFCALThetaBins, dMinFCALTheta, dMaxFCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
226 
227  gDirectory->cd("..");
228 
229  //CDC
230  CreateAndChangeTo_Directory("CDC", "CDC");
231 
232  locHistName = "dEdXVsP_Int";
233  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx [Integral-based] (keV/cm)");
234  dHistMap_dEdXVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
235 
236  locHistName = "dEdXVsP_Amp";
237  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx [Amplitude-based] (keV/cm)");
238  dHistMap_dEdXVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
239 
240  locHistName = "DeltadEdXVsP_Int";
241  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) [Integral-based] (keV/cm)");
242  dHistMap_DeltadEdXVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
243 
244  locHistName = "DeltadEdXVsP_Amp";
245  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) [Amplitude-based] (keV/cm)");
246  dHistMap_DeltadEdXVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
247 
248  locHistName = string("dEdXPullVsP_Int_") + locParticleName;
249  locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)} [Integral-based]");
250  dHistMap_dEdXPullVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
251 
252  locHistName = string("dEdXPullVsP_Amp_") + locParticleName;
253  locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)} [Amplitude-based]");
254  dHistMap_dEdXPullVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
255 
256  locHistName = string("dEdXFOMVsP_Int_") + locParticleName;
257  locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx [Integral-based] PID Confidence Level");
258  dHistMap_dEdXFOMVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
259 
260  locHistName = string("dEdXFOMVsP_Amp_") + locParticleName;
261  locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx [Amplitude-based] PID Confidence Level");
262  dHistMap_dEdXFOMVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
263 
264  gDirectory->cd("..");
265 
266  //FDC
267  CreateAndChangeTo_Directory("FDC", "FDC");
268 
269  locHistName = "dEdXVsP";
270  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC dE/dx (keV/cm)");
271  dHistMap_dEdXVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
272 
273  locHistName = "DeltadEdXVsP";
274  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX) (keV/cm)");
275  dHistMap_DeltadEdXVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
276 
277  locHistName = string("dEdXPullVsP_") + locParticleName;
278  locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
279  dHistMap_dEdXPullVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
280 
281  locHistName = string("dEdXFOMVsP_") + locParticleName;
282  locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC dE/dx PID Confidence Level");
283  dHistMap_dEdXFOMVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
284 
285  gDirectory->cd("..");
286 
287  // DIRC
288  CreateAndChangeTo_Directory("DIRC", "DIRC");
289 
290  locHistName = string("NumPhotons_") + locParticleName;
291  locHistTitle = locParticleROOTName + string("; DIRC NumPhotons");
292  dHistMap_NumPhotons_DIRC[locPID] = new TH1I(locHistName.c_str(), locHistTitle.c_str(), dDIRCNumPhotonsBins, dDIRCMinNumPhotons, dDIRCMaxNumPhotons);
293 
294  locHistName = string("ThetaCVsP_") + locParticleName;
295  locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC #theta_{C}");
296  dHistMap_ThetaCVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCThetaCBins, dDIRCMinThetaC, dDIRCMaxThetaC);
297 
298  locHistName = string("Ldiff_kpiVsP_") + locParticleName;
299  locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC L_{K}-L_{#pi}");
300  dHistMap_Ldiff_kpiVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCLikelihoodBins, -1*dDIRCMaxLikelihood, dDIRCMaxLikelihood);
301 
302  locHistName = string("Ldiff_pkVsP_") + locParticleName;
303  locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC L_{p}-L_{K}");
304  dHistMap_Ldiff_pkVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCLikelihoodBins, -1*dDIRCMaxLikelihood, dDIRCMaxLikelihood);
305 
306  gDirectory->cd("..");
307  } //end of charged
308 
309  //one per thrown pid:
310  if(!locMCThrowns.empty())
311  {
312  CreateAndChangeTo_Directory("Throwns", "Throwns");
313  for(size_t loc_j = 0; loc_j < dThrownPIDs.size(); ++loc_j)
314  {
315  locPID2 = dThrownPIDs[loc_j];
316  if((ParticleCharge(locPID2) != ParticleCharge(locPID)) && (locPID2 != Unknown))
317  continue;
318  locParticleName2 = ParticleType(locPID2);
319  locParticleROOTName2 = ParticleName_ROOT(locPID2);
320 
321  //Confidence Level for Thrown PID
322  pair<Particle_t, Particle_t> locPIDPair(locPID, locPID2);
323  locHistName = string("PIDConfidenceLevel_ForThrownPID_") + locParticleName2;
324  locHistTitle = locParticleROOTName + string(" PID, Thrown PID = ") + locParticleROOTName2 + string(";PID Confidence Level");
325  dHistMap_PIDFOMForTruePID[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFOMBins, 0.0, 1.0);
326  }
327  gDirectory->cd("..");
328  }
329 
330  //no FOM for massive neutrals
331  if((ParticleCharge(locPID) != 0) || (locPID == Gamma))
332  {
333  // Overall Confidence Level
334  locHistName = "PIDConfidenceLevel";
335  locHistTitle = locParticleROOTName + string(" PID;PID Confidence Level");
336  dHistMap_PIDFOM[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFOMBins, 0.0, 1.0);
337 
338  // P Vs Theta, PID FOM = NaN
339  locHistName = "PVsTheta_NaNPIDFOM";
340  locHistTitle = locParticleROOTName + string(", PID FOM = NaN;#theta#circ;p (GeV/c)");
341  dHistMap_PVsTheta_NaNPIDFOM[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
342  }
343 
344  gDirectory->cd("..");
345  } //end of PID loop
346 
347  //Return to the base directory
349  }
350  japp->RootUnLock(); //RELEASE ROOT LOCK!!
351 }
352 
353 bool DHistogramAction_PID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
354 {
355  vector<const DMCThrownMatching*> locMCThrownMatchingVector;
356  locEventLoop->Get(locMCThrownMatchingVector);
357  const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector.empty() ? NULL : locMCThrownMatchingVector[0];
358 
359  const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();
360 
361  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
362  {
363  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);
364  auto locParticles = Get_UseKinFitResultsFlag() ? locParticleComboStep->Get_FinalParticles(Get_Reaction()->Get_ReactionStep(loc_i), false, false, d_AllCharges) : locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_AllCharges);
365  for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
366  {
367  //check if will be duplicate
368  pair<Particle_t, const JObject*> locParticleInfo(locParticles[loc_j]->PID(), Get_FinalParticle_SourceObject(locParticles[loc_j]));
369  pair<const DEventRFBunch*, pair<Particle_t, const JObject*> > locHistInfo(locEventRFBunch, locParticleInfo);
371  continue; //previously histogrammed
372 
373  dPreviouslyHistogrammedParticles.insert(locHistInfo);
374  if(ParticleCharge(locParticles[loc_j]->PID()) != 0) //charged
375  Fill_ChargedHists(static_cast<const DChargedTrackHypothesis*>(locParticles[loc_j]), locMCThrownMatching, locEventRFBunch);
376  else //neutral
377  Fill_NeutralHists(static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_j]), locMCThrownMatching, locEventRFBunch);
378  }
379  }
380 
381  return true;
382 }
383 
384 void DHistogramAction_PID::Fill_ChargedHists(const DChargedTrackHypothesis* locChargedTrackHypothesis, const DMCThrownMatching* locMCThrownMatching, const DEventRFBunch* locEventRFBunch)
385 {
386  Particle_t locPID = locChargedTrackHypothesis->PID();
387  double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
388  double locDeltaBeta = locBeta_Timing - locChargedTrackHypothesis->lorentzMomentum().Beta();
389  double locDeltaT = (locChargedTrackHypothesis->Get_TimeAtPOCAToVertex() - locChargedTrackHypothesis->t0());
390  double locFOM_Timing = (locChargedTrackHypothesis->Get_NDF_Timing() > 0) ? TMath::Prob(locChargedTrackHypothesis->Get_ChiSq_Timing(), locChargedTrackHypothesis->Get_NDF_Timing()) : numeric_limits<double>::quiet_NaN();
391 
392  double locP = locChargedTrackHypothesis->momentum().Mag();
393  double locTheta = locChargedTrackHypothesis->momentum().Theta()*180.0/TMath::Pi();
394  double locMatchFOM = 0.0;
395  const DMCThrown* locMCThrown = (locMCThrownMatching != NULL) ? locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM) : NULL;
396 
397  auto locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
398  auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
399  auto locTOFHitMatchParams = locChargedTrackHypothesis->Get_TOFHitMatchParams();
400  auto locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
401  auto locDIRCMatchParams = locChargedTrackHypothesis->Get_DIRCMatchParams();
402 
403  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
404 
405  double locTimePull = 0.0;
406  unsigned int locTimeNDF = 0;
407  dParticleID->Calc_TimingChiSq(locChargedTrackHypothesis, locTimeNDF, locTimePull);
408  DetectorSystem_t locTimeDetector = locChargedTrackHypothesis->t1_detector();
409 
410  //FILL HISTOGRAMS
411  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
412  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
413  Lock_Action();
414  {
415  dHistMap_PIDFOM[locPID]->Fill(locChargedTrackHypothesis->Get_FOM());
416 
417  //SC dE/dx
418  if(locSCHitMatchParams != NULL)
419  {
420  dHistMap_dEdXVsP[locPID][SYS_START]->Fill(locP, locSCHitMatchParams->dEdx*1.0E3);
421  double locdx = locSCHitMatchParams->dHitEnergy/locSCHitMatchParams->dEdx;
422  double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
423  dParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
424  dHistMap_DeltadEdXVsP[locPID][SYS_START]->Fill(locP, (locSCHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
425  }
426 
427  //TOF dE/dx
428  if(locTOFHitMatchParams != NULL)
429  {
430  dHistMap_dEdXVsP[locPID][SYS_TOF]->Fill(locP, locTOFHitMatchParams->dEdx*1.0E3);
431  double locdx = locTOFHitMatchParams->dHitEnergy/locTOFHitMatchParams->dEdx;
432  double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
433  dParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
434  dHistMap_DeltadEdXVsP[locPID][SYS_TOF]->Fill(locP, (locTOFHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
435  }
436 
437  //BCAL E/p
438  if(locBCALShowerMatchParams != NULL)
439  {
440  const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
441  double locEOverP = locBCALShower->E/locP;
442  dHistMap_EOverPVsP[locPID][SYS_BCAL]->Fill(locP, locEOverP);
443  dHistMap_EOverPVsTheta[locPID][SYS_BCAL]->Fill(locTheta, locEOverP);
444  }
445 
446  //FCAL E/p
447  if(locFCALShowerMatchParams != NULL)
448  {
449  const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
450  double locEOverP = locFCALShower->getEnergy()/locP;
451  dHistMap_EOverPVsP[locPID][SYS_FCAL]->Fill(locP, locEOverP);
452  dHistMap_EOverPVsTheta[locPID][SYS_FCAL]->Fill(locTheta, locEOverP);
453  }
454 
455  //DIRC
456  if(locDIRCMatchParams != NULL) {
457 
458  int locNumPhotons_DIRC = locDIRCMatchParams->dNPhotons;
459  double locThetaC_DIRC = locDIRCMatchParams->dThetaC;
460  dHistMap_NumPhotons_DIRC[locPID]->Fill(locNumPhotons_DIRC);
461  dHistMap_ThetaCVsP_DIRC[locPID]->Fill(locP, locThetaC_DIRC);
462  double locLpi_DIRC = locDIRCMatchParams->dLikelihoodPion;
463  double locLk_DIRC = locDIRCMatchParams->dLikelihoodKaon;
464  double locLp_DIRC = locDIRCMatchParams->dLikelihoodProton;
465  dHistMap_Ldiff_kpiVsP_DIRC[locPID]->Fill(locP, locLk_DIRC-locLpi_DIRC);
466  dHistMap_Ldiff_pkVsP_DIRC[locPID]->Fill(locP, locLp_DIRC-locLk_DIRC);
467  }
468 
469  //Timing
470  if(dHistMap_BetaVsP[locPID].find(locTimeDetector) != dHistMap_BetaVsP[locPID].end())
471  {
472  dHistMap_BetaVsP[locPID][locTimeDetector]->Fill(locP, locBeta_Timing);
473  dHistMap_DeltaBetaVsP[locPID][locTimeDetector]->Fill(locP, locDeltaBeta);
474  dHistMap_DeltaTVsP[locPID][locTimeDetector]->Fill(locP, locDeltaT);
475  dHistMap_TimePullVsP[locPID][locTimeDetector]->Fill(locP, locTimePull);
476  dHistMap_TimeFOMVsP[locPID][locTimeDetector]->Fill(locP, locFOM_Timing);
477  }
478 
479  //CDC dE/dx
480  if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC > 0)
481  {
482  dHistMap_dEdXVsP[locPID][SYS_CDC]->Fill(locP, locTrackTimeBased->ddEdx_CDC*1.0E6);
483  dHistMap_dEdXVsP[locPID][SYS_CDC_AMP]->Fill(locP, locTrackTimeBased->ddEdx_CDC_amp*1.0E6);
484 
485  double locProbabledEdx = dParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_CDC, true);
486  double locDeltadEdx = locTrackTimeBased->ddEdx_CDC - locProbabledEdx;
487  dHistMap_DeltadEdXVsP[locPID][SYS_CDC]->Fill(locP, 1.0E6*locDeltadEdx);
488  double locDeltadEdx_amp = locTrackTimeBased->ddEdx_CDC_amp - locProbabledEdx;
489  dHistMap_DeltadEdXVsP[locPID][SYS_CDC_AMP]->Fill(locP, 1.0E6*locDeltadEdx_amp);
490 
491  double locMeandx = locTrackTimeBased->ddx_CDC/locTrackTimeBased->dNumHitsUsedFordEdx_CDC;
492  double locSigmadEdx = dParticleID->GetdEdxSigma_DC(locTrackTimeBased->dNumHitsUsedFordEdx_CDC, locP, locChargedTrackHypothesis->mass(), locMeandx, true);
493  double locdEdXPull = locDeltadEdx/locSigmadEdx;
494  dHistMap_dEdXPullVsP[locPID][SYS_CDC]->Fill(locP, locDeltadEdx/locSigmadEdx);
495  double locdEdXPull_amp = locDeltadEdx_amp/locSigmadEdx;
496  dHistMap_dEdXPullVsP[locPID][SYS_CDC_AMP]->Fill(locP, locDeltadEdx_amp/locSigmadEdx);
497 
498  double locdEdXChiSq = locdEdXPull*locdEdXPull;
499  double locdEdXFOM = TMath::Prob(locdEdXChiSq, locTrackTimeBased->dNumHitsUsedFordEdx_CDC);
500  dHistMap_dEdXFOMVsP[locPID][SYS_CDC]->Fill(locP, locdEdXFOM);
501  double locdEdXChiSq_amp = locdEdXPull_amp*locdEdXPull_amp;
502  double locdEdXFOM_amp = TMath::Prob(locdEdXChiSq_amp, locTrackTimeBased->dNumHitsUsedFordEdx_CDC);
503  dHistMap_dEdXFOMVsP[locPID][SYS_CDC_AMP]->Fill(locP, locdEdXFOM_amp);
504  }
505 
506  //FDC dE/dx
507  if(locTrackTimeBased->dNumHitsUsedFordEdx_FDC > 0)
508  {
509  dHistMap_dEdXVsP[locPID][SYS_FDC]->Fill(locP, locTrackTimeBased->ddEdx_FDC*1.0E6);
510  double locProbabledEdx = dParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_FDC, false);
511  double locDeltadEdx = locTrackTimeBased->ddEdx_FDC - locProbabledEdx;
512  dHistMap_DeltadEdXVsP[locPID][SYS_FDC]->Fill(locP, 1.0E6*locDeltadEdx);
513  double locMeandx = locTrackTimeBased->ddx_FDC/locTrackTimeBased->dNumHitsUsedFordEdx_FDC;
514  double locSigmadEdx = dParticleID->GetdEdxSigma_DC(locTrackTimeBased->dNumHitsUsedFordEdx_FDC, locP, locChargedTrackHypothesis->mass(), locMeandx, false);
515  double locdEdXPull = locDeltadEdx/locSigmadEdx;
516  dHistMap_dEdXPullVsP[locPID][SYS_FDC]->Fill(locP, locDeltadEdx/locSigmadEdx);
517  double locdEdXChiSq = locdEdXPull*locdEdXPull;
518  double locdEdXFOM = TMath::Prob(locdEdXChiSq, locTrackTimeBased->dNumHitsUsedFordEdx_FDC);
519  dHistMap_dEdXFOMVsP[locPID][SYS_FDC]->Fill(locP, locdEdXFOM);
520  }
521 
522  pair<Particle_t, Particle_t> locPIDPair(locPID, Unknown); //default unless matched
523  if(locMCThrown != NULL) //else bogus track (not matched to any thrown tracks)
524  locPIDPair.second = (Particle_t)(locMCThrown->type); //matched
525  if(dHistMap_PIDFOMForTruePID.find(locPIDPair) != dHistMap_PIDFOMForTruePID.end()) //else hist not created or PID is weird
526  dHistMap_PIDFOMForTruePID[locPIDPair]->Fill(locChargedTrackHypothesis->Get_FOM());
527 
528  if(locChargedTrackHypothesis->Get_NDF() == 0) //NaN
529  dHistMap_PVsTheta_NaNPIDFOM[locPID]->Fill(locTheta, locP);
530  }
531  Unlock_Action();
532 }
533 
534 void DHistogramAction_PID::Fill_NeutralHists(const DNeutralParticleHypothesis* locNeutralParticleHypothesis, const DMCThrownMatching* locMCThrownMatching, const DEventRFBunch* locEventRFBunch)
535 {
536  Particle_t locPID = locNeutralParticleHypothesis->PID();
537 
538  double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
539  double locDeltaT = locNeutralParticleHypothesis->time() - locNeutralParticleHypothesis->t0();
540 
541  double locP = locNeutralParticleHypothesis->momentum().Mag();
542  double locMatchFOM = 0.0;
543  const DMCThrown* locMCThrown = (locMCThrownMatching != NULL) ? locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM) : NULL;
544 
545  double locTimePull = 0.0;
546  unsigned int locTimeNDF = 0;
547  dParticleID->Calc_TimingChiSq(locNeutralParticleHypothesis, locTimeNDF, locTimePull);
548 
549  DetectorSystem_t locSystem = locNeutralParticleHypothesis->t1_detector();
550 
551  //FILL HISTOGRAMS
552  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
553  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
554  Lock_Action();
555  {
556  //Beta (good for all PIDs)
557  dHistMap_Beta[locPID][locSystem]->Fill(locBeta_Timing);
558  if(locPID != Gamma)
559  {
560  Unlock_Action();
561  return;
562  }
563 
564  dHistMap_PIDFOM[locPID]->Fill(locNeutralParticleHypothesis->Get_FOM());
565  dHistMap_BetaVsP[locPID][locSystem]->Fill(locP, locBeta_Timing);
566  dHistMap_DeltaTVsP[locPID][locSystem]->Fill(locP, locDeltaT);
567  dHistMap_TimePullVsP[locPID][locSystem]->Fill(locP, locTimePull);
568  dHistMap_TimeFOMVsP[locPID][locSystem]->Fill(locP, locNeutralParticleHypothesis->Get_FOM());
569 
570  pair<Particle_t, Particle_t> locPIDPair(locPID, Unknown); //default unless matched
571  if(locMCThrown != NULL) //else bogus track (not matched to any thrown tracks)
572  locPIDPair.second = (Particle_t)(locMCThrown->type); //matched
573  if(dHistMap_PIDFOMForTruePID.find(locPIDPair) != dHistMap_PIDFOMForTruePID.end()) //else hist not created or PID is weird
574  dHistMap_PIDFOMForTruePID[locPIDPair]->Fill(locNeutralParticleHypothesis->Get_FOM());
575  }
576  Unlock_Action();
577 }
578 
580 {
581  string locHistName, locHistTitle, locStepROOTName, locParticleName, locParticleROOTName;
582  Particle_t locPID, locHigherMassPID, locLowerMassPID;
583  string locHigherMassParticleName, locLowerMassParticleName, locHigherMassParticleROOTName, locLowerMassParticleROOTName;
584 
585  vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;
586  locEventLoop->Get(locAnalysisUtilitiesVector);
587 
588  size_t locNumSteps = Get_Reaction()->Get_NumReactionSteps();
589  dHistDeque_TrackZToCommon.resize(locNumSteps);
590  dHistDeque_TrackTToCommon.resize(locNumSteps);
591  dHistDeque_TrackDOCAToCommon.resize(locNumSteps);
592  dHistDeque_MaxTrackDeltaZ.resize(locNumSteps);
593  dHistDeque_MaxTrackDeltaT.resize(locNumSteps);
594  dHistDeque_MaxTrackDOCA.resize(locNumSteps);
595  dHistDeque_TrackDeltaTVsP.resize(locNumSteps);
596 
597  //CREATE THE HISTOGRAMS
598  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
599  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
600  {
601  dAnalysisUtilities = locAnalysisUtilitiesVector[0];
603  for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)
604  {
605  auto locDetectedChargedPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_Charged, false);
606  if(locDetectedChargedPIDs.empty())
607  continue;
608 
609  const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
610  ostringstream locStepName;
611  locStepName << "Step" << loc_i << "__" << DAnalysis::Get_StepName(locReactionStep, true, false);
612  locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);
613  CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());
614 
615  // Max Track DeltaZ
616  locHistName = "MaxTrackDeltaZ";
617  locHistTitle = locStepROOTName + string(";Largest Track #DeltaVertex-Z (cm)");
618  dHistDeque_MaxTrackDeltaZ[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexZBins, 0.0, dMaxDeltaVertexZ);
619 
620  // Max Track DeltaT
621  locHistName = "MaxTrackDeltaT";
622  locHistTitle = locStepROOTName + string(";Largest Track #DeltaVertex-T (ns)");
623  dHistDeque_MaxTrackDeltaT[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexTBins, 0.0, dMaxDeltaVertexT);
624 
625  // Max Track DOCA
626  locHistName = "MaxTrackDOCA";
627  locHistTitle = locStepROOTName + string(";Largest Track DOCA (cm)");
628  dHistDeque_MaxTrackDOCA[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDOCABins, 0.0, dMaxDOCA);
629 
630  for(size_t loc_j = 0; loc_j < locDetectedChargedPIDs.size(); ++loc_j)
631  {
632  locPID = locDetectedChargedPIDs[loc_j];
633  if(dHistDeque_TrackZToCommon[loc_i].find(locPID) != dHistDeque_TrackZToCommon[loc_i].end())
634  continue; //already created for this pid
635  locParticleName = ParticleType(locPID);
636  locParticleROOTName = ParticleName_ROOT(locPID);
637 
638  // TrackZ To Common
639  locHistName = string("TrackZToCommon_") + locParticleName;
640  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#DeltaVertex-Z (Track, Common) (cm)");
641  dHistDeque_TrackZToCommon[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexZBins, dMinDeltaVertexZ, dMaxDeltaVertexZ);
642 
643  // TrackT To Common
644  locHistName = string("TrackTToCommon_") + locParticleName;
645  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#DeltaVertex-T (Track, Common) (ns)");
646  dHistDeque_TrackTToCommon[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexTBins, dMinDeltaVertexT, dMaxDeltaVertexT);
647 
648  // TrackDOCA To Common
649  locHistName = string("TrackDOCAToCommon_") + locParticleName;
650  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";DOCA (Track, Common) (cm)");
651  dHistDeque_TrackDOCAToCommon[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDOCABins, dMinDOCA, dMaxDOCA);
652 
653  // DeltaT Vs P against beam photon
655  {
656  locHistName = string("TrackDeltaTVsP_") + ParticleType(locPID) + string("_Beam") + ParticleType(locPID);
657  locHistTitle = locStepROOTName + string(";") + ParticleName_ROOT(locPID) + string(" Momentum (GeV/c);t_{") + ParticleName_ROOT(locPID) + string("} - t_{Beam ") + ParticleName_ROOT(locPID) + string("} (ns)");
658  dHistMap_BeamTrackDeltaTVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaVertexTBins, dMinDeltaVertexT, dMaxDeltaVertexT);
659  }
660  }
661 
662  //delta-t vs p
663  auto locDetectedChargedPIDs_HasDupes = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_Charged, true);
664  for(int loc_j = 0; loc_j < int(locDetectedChargedPIDs_HasDupes.size()) - 1; ++loc_j)
665  {
666  locPID = locDetectedChargedPIDs_HasDupes[loc_j];
667  for(size_t loc_k = loc_j + 1; loc_k < locDetectedChargedPIDs_HasDupes.size(); ++loc_k)
668  {
669  if(ParticleMass(locDetectedChargedPIDs_HasDupes[loc_k]) > ParticleMass(locPID))
670  {
671  locHigherMassPID = locDetectedChargedPIDs_HasDupes[loc_k];
672  locLowerMassPID = locPID;
673  }
674  else
675  {
676  locHigherMassPID = locPID;
677  locLowerMassPID = locDetectedChargedPIDs_HasDupes[loc_k];
678  }
679  pair<Particle_t, Particle_t> locParticlePair(locHigherMassPID, locLowerMassPID);
680  if(dHistDeque_TrackDeltaTVsP[loc_i].find(locParticlePair) != dHistDeque_TrackDeltaTVsP[loc_i].end())
681  continue; //already created for this pair
682 
683  locHigherMassParticleName = ParticleType(locHigherMassPID);
684  locHigherMassParticleROOTName = ParticleName_ROOT(locHigherMassPID);
685  locLowerMassParticleName = ParticleType(locLowerMassPID);
686  locLowerMassParticleROOTName = ParticleName_ROOT(locLowerMassPID);
687 
688  // DeltaT Vs P
689  locHistName = string("TrackDeltaTVsP_") + locHigherMassParticleName + string("_") + locLowerMassParticleName;
690  locHistTitle = locStepROOTName + string(";") + locHigherMassParticleROOTName + string(" Momentum (GeV/c);t_{") + locHigherMassParticleROOTName + string("} - t_{") + locLowerMassParticleROOTName + string("} (ns)");
691  dHistDeque_TrackDeltaTVsP[loc_i][locParticlePair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaVertexTBins, dMinDeltaVertexT, dMaxDeltaVertexT);
692  }
693  }
694  gDirectory->cd("..");
695  }
696  //Return to the base directory
698  }
699  japp->RootUnLock(); //RELEASE ROOT LOCK!!
700 }
701 
702 bool DHistogramAction_TrackVertexComparison::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
703 {
704  Particle_t locPID;
705  double locDOCA, locDeltaZ, locDeltaT, locMaxDOCA, locMaxDeltaZ, locMaxDeltaT;
706 
707  //note: very difficult to tell when results will be duplicate: just histogram all combos
708  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
709  {
710  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);
711  auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_Charged);
712  if(locParticles.empty())
713  continue;
714 
715  //Grab/Find common vertex & time
716  DVector3 locVertex = locParticleComboStep->Get_Position();
717  double locVertexTime = locParticleComboStep->Get_Time();
718 
719  for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
720  {
721  locPID = locParticles[loc_j]->PID();
722 
723  //find max's
724  locMaxDOCA = -1.0;
725  locMaxDeltaZ = -1.0;
726  locMaxDeltaT = -1.0;
727  for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
728  {
729  locDeltaZ = fabs(locParticles[loc_j]->position().Z() - locParticles[loc_k]->position().Z());
730  if(locDeltaZ > locMaxDeltaZ)
731  locMaxDeltaZ = locDeltaZ;
732 
733  locDeltaT = fabs(locParticles[loc_j]->time() - locParticles[loc_k]->time());
734  if(locDeltaT > locMaxDeltaT)
735  locMaxDeltaT = locDeltaT;
736 
737  locDOCA = dAnalysisUtilities->Calc_DOCA(locParticles[loc_j], locParticles[loc_k]);
738  if(locDOCA > locMaxDOCA)
739  locMaxDOCA = locDOCA;
740  }
741 
742  //delta-t vs p
743  deque<pair<const DKinematicData*, size_t> > locParticlePairs;
744  size_t locHigherMassParticleIndex, locLowerMassParticleIndex;
745  //minimize number of locks by keeping track of the results and saving them at the end
746  deque<pair<Particle_t, Particle_t> > locPIDPairs;
747  deque<double> locPs;
748  deque<double> locDeltaTs;
749  for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)
750  {
751  locParticlePairs.clear();
752  locParticlePairs.push_back(pair<const DKinematicData*, size_t>(locParticles[loc_j], loc_i));
753  locParticlePairs.push_back(pair<const DKinematicData*, size_t>(locParticles[loc_k], loc_i));
754 
755  if(locParticles[loc_k]->mass() > locParticles[loc_j]->mass())
756  {
757  locHigherMassParticleIndex = loc_k;
758  locLowerMassParticleIndex = loc_j;
759  }
760  else
761  {
762  locHigherMassParticleIndex = loc_j;
763  locLowerMassParticleIndex = loc_k;
764  }
765 
766  locDeltaTs.push_back(locParticles[locHigherMassParticleIndex]->time() - locParticles[locLowerMassParticleIndex]->time());
767  locPs.push_back(locParticles[locHigherMassParticleIndex]->momentum().Mag());
768  locPIDPairs.push_back(pair<Particle_t, Particle_t>(locParticles[locHigherMassParticleIndex]->PID(), locParticles[locLowerMassParticleIndex]->PID()));
769  }
770 
771  const DKinematicData* locBeamParticle = ((loc_i == 0) && DAnalysis::Get_IsFirstStepBeam(Get_Reaction())) ? locParticleComboStep->Get_InitialParticle() : NULL;
772  double locBeamDeltaT = (locBeamParticle != NULL) ? locParticles[loc_j]->time() - locBeamParticle->time() : numeric_limits<double>::quiet_NaN();
773  locDOCA = dAnalysisUtilities->Calc_DOCAToVertex(locParticles[loc_j], locVertex);
774 
775  //do all at once to reduce #locks & amount of time within the lock
776  //FILL HISTOGRAMS
777  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
778  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
779  Lock_Action();
780  {
781  //comparison to common vertex/time
782  dHistDeque_TrackZToCommon[loc_i][locPID]->Fill(locParticles[loc_j]->position().Z() - locVertex.Z());
783  dHistDeque_TrackTToCommon[loc_i][locPID]->Fill(locParticles[loc_j]->time() - locVertexTime);
784  dHistDeque_TrackDOCAToCommon[loc_i][locPID]->Fill(locDOCA);
785  //hist max's
786  if(locMaxDeltaZ > 0.0) //else none found (e.g. only 1 detected charged track)
787  {
788  dHistDeque_MaxTrackDeltaZ[loc_i]->Fill(locMaxDeltaZ);
789  dHistDeque_MaxTrackDeltaT[loc_i]->Fill(locMaxDeltaT);
790  dHistDeque_MaxTrackDOCA[loc_i]->Fill(locMaxDOCA);
791  }
792  //delta-t's
793  if(locBeamParticle != NULL)
794  dHistMap_BeamTrackDeltaTVsP[locPID]->Fill(locParticles[loc_j]->momentum().Mag(), locBeamDeltaT);
795  for(size_t loc_k = 0; loc_k < locPIDPairs.size(); ++loc_k)
796  {
797  if(dHistDeque_TrackDeltaTVsP[loc_i].find(locPIDPairs[loc_k]) == dHistDeque_TrackDeltaTVsP[loc_i].end())
798  {
799  //pair not found: equal masses and order switched somehow //e.g. mass set differently between REST and reconstruction
800  pair<Particle_t, Particle_t> locTempPIDPair(locPIDPairs[loc_k]);
801  locPIDPairs[loc_k].first = locTempPIDPair.second;
802  locPIDPairs[loc_k].second = locTempPIDPair.first;
803  }
804  dHistDeque_TrackDeltaTVsP[loc_i][locPIDPairs[loc_k]]->Fill(locPs[loc_k], locDeltaTs[loc_k]);
805  }
806  }
807  Unlock_Action();
808  } //end of particle loop
809  } //end of step loop
810  return true;
811 }
812 
814 {
815  if(Get_UseKinFitResultsFlag() && (Get_Reaction()->Get_KinFitType() == d_NoFit))
816  {
817  cout << "WARNING: REQUESTED HISTOGRAM OF KINEMATIC FIT RESULTS WHEN NO KINEMATIC FIT!!!" << endl;
818  return; //no fit performed, but kinfit data requested!!
819  }
820 
821  vector<const DParticleID*> locParticleIDs;
822  locEventLoop->Get(locParticleIDs);
823 
824  string locHistName, locHistTitle, locStepROOTName, locParticleName, locParticleROOTName;
825  Particle_t locPID;
826 
827  size_t locNumSteps = Get_Reaction()->Get_NumReactionSteps();
828  dHistDeque_PVsTheta.resize(locNumSteps);
829  dHistDeque_PhiVsTheta.resize(locNumSteps);
830  dHistDeque_BetaVsP.resize(locNumSteps);
831  dHistDeque_DeltaBetaVsP.resize(locNumSteps);
832  dHistDeque_P.resize(locNumSteps);
833  dHistDeque_Theta.resize(locNumSteps);
834  dHistDeque_Phi.resize(locNumSteps);
835  dHistDeque_VertexZ.resize(locNumSteps);
836  dHistDeque_VertexYVsX.resize(locNumSteps);
837 
838  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
839  DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
840  vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;
841  locEventLoop->Get(locAnalysisUtilitiesVector);
842 
843  auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());
844 
845  //CREATE THE HISTOGRAMS
846  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
847  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
848  {
850 
851  dParticleID = locParticleIDs[0];
852  dAnalysisUtilities = locAnalysisUtilitiesVector[0];
853  locGeometry->GetTargetZ(dTargetZCenter);
854 
855  //beam particle
856  locPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID();
857  if(locIsFirstStepBeam)
858  {
859  locParticleName = string("Beam_") + ParticleType(locPID);
860  CreateAndChangeTo_Directory(locParticleName, locParticleName);
861  locParticleROOTName = ParticleName_ROOT(locPID);
862 
863  // Momentum
864  locHistName = "Momentum";
865  locHistTitle = string("Beam ") + locParticleROOTName + string(";p (GeV/c)");
866  dBeamParticleHist_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
867 
868  // Theta
869  locHistName = "Theta";
870  locHistTitle = string("Beam ") + locParticleROOTName + string(";#theta#circ");
871  dBeamParticleHist_Theta = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
872 
873  // Phi
874  locHistName = "Phi";
875  locHistTitle = string("Beam ") + locParticleROOTName + string(";#phi#circ");
876  dBeamParticleHist_Phi = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
877 
878  // P Vs Theta
879  locHistName = "PVsTheta";
880  locHistTitle = string("Beam ") + locParticleROOTName + string(";#theta#circ;p (GeV/c)");
881  dBeamParticleHist_PVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
882 
883  // Phi Vs Theta
884  locHistName = "PhiVsTheta";
885  locHistTitle = string("Beam ") + locParticleROOTName + string(";#theta#circ;#phi#circ");
886  dBeamParticleHist_PhiVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
887 
888  // Vertex-Z
889  locHistName = "VertexZ";
890  locHistTitle = string("Beam ") + locParticleROOTName + string(";Vertex-Z (cm)");
891  dBeamParticleHist_VertexZ = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
892 
893  // Vertex-Y Vs Vertex-X
894  locHistName = "VertexYVsX";
895  locHistTitle = string("Beam ") + locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
896  dBeamParticleHist_VertexYVsX = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
897 
898  // Delta-T (Beam, RF)
899  locHistName = "DeltaTRF";
900  locHistTitle = string("Beam ") + locParticleROOTName + string(";#Deltat_{Beam - RF} (ns)");
901  dBeamParticleHist_DeltaTRF = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTRFBins, dMinDeltaTRF, dMaxDeltaTRF);
902 
903  // Delta-T (Beam, RF) Vs Beam E
904  locHistName = "DeltaTRFVsBeamE";
905  locHistTitle = string("Beam ") + locParticleROOTName + string(";E (GeV);#Deltat_{Beam - RF} (ns)");
906  dBeamParticleHist_DeltaTRFVsBeamE = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaTRFBins, dMinDeltaTRF, dMaxDeltaTRF);
907 
908  gDirectory->cd("..");
909  }
910 
911  //Steps
912  for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)
913  {
914  const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
915  ostringstream locStepName;
916  locStepName << "Step" << loc_i << "__" << DAnalysis::Get_StepName(locReactionStep, true, false);
917  locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);
918 
919  Particle_t locInitialPID = locReactionStep->Get_InitialPID();
920 
921 
922  //get PIDs
923  vector<Particle_t> locPIDs;
924  bool locLastPIDMissingFlag = false;
925  if(!Get_UseKinFitResultsFlag()) //measured, ignore missing & decaying particles (ignore target anyway)
926  locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_AllCharges, false);
927  else //kinematic fit: decaying & missing particles are reconstructed
928  {
929  locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, true, d_AllCharges, false);
930  if((!locIsFirstStepBeam) || (loc_i != 0)) //add decaying parent particle //skip if on beam particle!
931  locPIDs.insert(locPIDs.begin(), locInitialPID);
932 
933  Particle_t locMissingPID = locReactionStep->Get_MissingPID();
934  if(locMissingPID != Unknown)
935  {
936  locPIDs.push_back(locMissingPID);
937  locLastPIDMissingFlag = true;
938  }
939  }
940 
941  bool locDirectoryCreatedFlag = false;
942  for(size_t loc_j = 0; loc_j < locPIDs.size(); ++loc_j)
943  {
944  locPID = locPIDs[loc_j];
945  bool locIsMissingFlag = false;
946  if((loc_j == locPIDs.size() - 1) && locLastPIDMissingFlag)
947  {
948  locParticleName = string("(") + ParticleType(locPID) + string(")");
949  locParticleROOTName = string("(") + ParticleName_ROOT(locPID) + string(")");
950  locIsMissingFlag = true;
951  }
952  else
953  {
954  locParticleName = ParticleType(locPID);
955  locParticleROOTName = ParticleName_ROOT(locPID);
956  if(dHistDeque_P[loc_i].find(locPID) != dHistDeque_P[loc_i].end())
957  continue; //pid already done
958  }
959 
960  if(!locDirectoryCreatedFlag)
961  {
962  CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());
963  locDirectoryCreatedFlag = true;
964  }
965 
966  CreateAndChangeTo_Directory(locParticleName, locParticleName);
967 
968  // Momentum
969  locHistName = "Momentum";
970  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";p (GeV/c)");
971  dHistDeque_P[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
972 
973  // Theta
974  locHistName = "Theta";
975  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#theta#circ");
976  dHistDeque_Theta[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
977 
978  // Phi
979  locHistName = "Phi";
980  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#phi#circ");
981  dHistDeque_Phi[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
982 
983  // P Vs Theta
984  locHistName = "PVsTheta";
985  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#theta#circ;p (GeV/c)");
986  dHistDeque_PVsTheta[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
987 
988  // Phi Vs Theta
989  locHistName = "PhiVsTheta";
990  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#theta#circ;#phi#circ");
991  dHistDeque_PhiVsTheta[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
992 
993  //beta vs p
994  locHistName = "BetaVsP";
995  locHistTitle = locParticleROOTName + string(";p (GeV/c);#beta");
996  dHistDeque_BetaVsP[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
997 
998  //delta-beta vs p
999  locHistName = "DeltaBetaVsP";
1000  locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#beta");
1001  dHistDeque_DeltaBetaVsP[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1002 
1003  // Vertex-Z
1004  locHistName = "VertexZ";
1005  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";Vertex-Z (cm)");
1006  dHistDeque_VertexZ[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
1007 
1008  // Vertex-Y Vs Vertex-X
1009  locHistName = "VertexYVsX";
1010  locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
1011  dHistDeque_VertexYVsX[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
1012 
1013  gDirectory->cd("..");
1014  } //end of particle loop
1015 
1016  // Vertex
1017  string locInitParticleROOTName = ParticleName_ROOT(locInitialPID);
1018  string locInitParticleName = ParticleType(locInitialPID);
1019  if((loc_i == 0) || IsDetachedVertex(locInitialPID))
1020  {
1021  if(!locDirectoryCreatedFlag)
1022  {
1023  CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());
1024  locDirectoryCreatedFlag = true;
1025  }
1026 
1027  locHistName = "StepVertexZ";
1028  locHistTitle = (loc_i == 0) ? ";Production Vertex-Z (cm)" : string(";") + locInitParticleROOTName + string(" Decay Vertex-Z (cm)");
1029  dHistMap_StepVertexZ[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
1030 
1031  locHistName = "StepVertexYVsX";
1032  locHistTitle = (loc_i == 0) ? "Production Vertex" : locInitParticleROOTName + string(" Decay Vertex");
1033  locHistTitle += string(";Vertex-X (cm);Vertex-Y (cm)");
1034  dHistMap_StepVertexYVsX[loc_i] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
1035  }
1036 
1037  if((loc_i != 0) && IsDetachedVertex(locInitialPID))
1038  {
1039  if(!locDirectoryCreatedFlag)
1040  {
1041  CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());
1042  locDirectoryCreatedFlag = true;
1043  }
1044 
1045  locHistName = locInitParticleName + string("PathLength");
1046  locHistTitle = string(";") + locInitParticleROOTName + string(" Path Length (cm)");
1047  dHistMap_DetachedPathLength[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPathLengthBins, 0.0, dMaxPathLength);
1048 
1049  locHistName = locInitParticleName + string("Lifetime");
1050  locHistTitle = string(";") + locInitParticleROOTName + string(" Lifetime (ns)");
1051  dHistMap_DetachedLifetime[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumLifetimeBins, 0.0, dMaxLifetime);
1052 
1054  {
1055  locHistName = locInitParticleName + string("Lifetime_RestFrame");
1056  locHistTitle = string(";") + locInitParticleROOTName + string(" Rest Frame Lifetime (ns)");
1057  dHistMap_DetachedLifetime_RestFrame[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumLifetimeBins, 0.0, dMaxLifetime);
1058  }
1059  }
1060 
1061  if(locDirectoryCreatedFlag)
1062  gDirectory->cd("..");
1063  } //end of step loop
1064 
1065  //Return to the base directory
1067  }
1068  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1069 }
1070 
1071 bool DHistogramAction_ParticleComboKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1072 {
1073  if(Get_UseKinFitResultsFlag() && (Get_Reaction()->Get_KinFitType() == d_NoFit))
1074  {
1075  cout << "WARNING: REQUESTED HISTOGRAM OF KINEMATIC FIT RESULTS WHEN NO KINEMATIC FIT!!! Skipping histogram." << endl;
1076  return true; //no fit performed, but kinfit data requested!!
1077  }
1078 
1079  const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();
1080  auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());
1081 
1082  const DKinematicData* locKinematicData;
1083  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
1084  {
1085  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);
1086  auto locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
1087 
1088  //initial particle
1089  Particle_t locInitialPID = locReactionStep->Get_InitialPID();
1091  locKinematicData = locParticleComboStep->Get_InitialParticle();
1092  else
1093  locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();
1094  if(locKinematicData != NULL)
1095  {
1096  if(locIsFirstStepBeam && (loc_i == 0))
1097  {
1098  //check if will be duplicate
1099  const JObject* locSourceObject = locParticleComboStep->Get_InitialParticle_Measured();
1101  {
1102  dPreviouslyHistogrammedBeamParticles.insert(locSourceObject);
1103  Fill_BeamHists(locKinematicData, locEventRFBunch);
1104  }
1105  }
1106  else if(Get_UseKinFitResultsFlag()) //decaying particle, but kinfit so can hist
1107  Fill_Hists(locEventLoop, locKinematicData, false, loc_i); //has many source object, and is unique to this combo: no dupes to check against: let it ride
1108  }
1109 
1110  //VERTEX INFORMATION
1111  //other than first, skipped if not detached vertex
1112  DLorentzVector locStepSpacetimeVertex = locParticleComboStep->Get_SpacetimeVertex();
1113  if((loc_i == 0) || IsDetachedVertex(locInitialPID))
1114  {
1115  dHistMap_StepVertexZ[loc_i]->Fill(locStepSpacetimeVertex.Z());
1116  dHistMap_StepVertexYVsX[loc_i]->Fill(locStepSpacetimeVertex.X(), locStepSpacetimeVertex.Y());
1117  }
1118 
1119  //DETACHED VERTEX INFORMATION
1120  if((loc_i != 0) && IsDetachedVertex(locInitialPID))
1121  {
1122  int locFromStepIndex = DAnalysis::Get_InitialParticleDecayFromIndices(Get_Reaction(), loc_i).first;
1123  DLorentzVector locFromSpacetimeVertex = locParticleCombo->Get_ParticleComboStep(locFromStepIndex)->Get_SpacetimeVertex();
1124  DLorentzVector locDeltaSpacetime = locStepSpacetimeVertex - locFromSpacetimeVertex;
1125 
1126  double locPathLength = locDeltaSpacetime.Vect().Mag();
1127  dHistMap_DetachedPathLength[loc_i]->Fill(locPathLength);
1128  dHistMap_DetachedLifetime[loc_i]->Fill(locDeltaSpacetime.T());
1129 
1130  if(locKinematicData != NULL)
1131  {
1132  DLorentzVector locInitialP4 = locKinematicData->lorentzMomentum();
1133  //below, t, x, and tau are really delta-t, delta-x, and delta-tau
1134  //tau = gamma*(t - beta*x/c) //plug in: t = x/(beta*c), gamma = 1/sqrt(1 - beta^2)
1135  //tau = (x/(beta*c) - beta*x/c)/sqrt(1 - beta^2) //plug in beta = p*c/E
1136  //tau = (x*E/(p*c^2) - p*x/E)/sqrt(1 - p^2*c^2/E^2) //multiply num & denom by E*P, factor x/c^2 out of num
1137  //tau = (x/c^2)*(E^2 - p^2*c^2)/(p*sqrt(E^2 - p^2*c^2)) //plug in m^2*c^4 = E^2 - p^2*c^2
1138  //tau = (x/c^2)*m^2*c^4/(p*m*c^2) //cancel c's & m's
1139  //tau = x*m/p
1140  //however, in data, p & m are in units with c = 1, so need an extra 1/c: tau = x*m/(c*p)
1141  double locRestFrameLifetime = locPathLength*ParticleMass(locInitialPID)/(29.9792458*locInitialP4.P()); //tau
1142  dHistMap_DetachedLifetime_RestFrame[loc_i]->Fill(locRestFrameLifetime);
1143  //note that tau = hbar / Gamma, hbar = 6.582119E-22 MeV*s, Gamma = Resonance FWHM
1144  }
1145  }
1146 
1147  //final particles
1148  for(size_t loc_j = 0; loc_j < locParticleComboStep->Get_NumFinalParticles(); ++loc_j)
1149  {
1151  locKinematicData = locParticleComboStep->Get_FinalParticle(loc_j);
1152  else
1153  locKinematicData = locParticleComboStep->Get_FinalParticle_Measured(loc_j);
1154  if(locKinematicData == NULL)
1155  continue; //e.g. a decaying or missing particle whose params aren't set yet
1156 
1157  //check if duplicate
1158  const JObject* locSourceObject = Get_FinalParticle_SourceObject(locKinematicData);
1159  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
1160  {
1161  pair<Particle_t, const JObject*> locParticleInfo(locKinematicData->PID(), locSourceObject);
1162  pair<size_t, pair<Particle_t, const JObject*> > locHistInfo(loc_i, locParticleInfo);
1164  continue; //previously histogrammed
1165  dPreviouslyHistogrammedParticles.insert(locHistInfo);
1166  }
1167 
1168  bool locIsMissingFlag = (Get_Reaction()->Get_ReactionStep(loc_i)->Get_MissingParticleIndex() == int(loc_j));
1169  Fill_Hists(locEventLoop, locKinematicData, locIsMissingFlag, loc_i);
1170  } //end of particle loop
1171  } //end of step loop
1172  return true;
1173 }
1174 
1175 void DHistogramAction_ParticleComboKinematics::Fill_Hists(JEventLoop* locEventLoop, const DKinematicData* locKinematicData, bool locIsMissingFlag, size_t locStepIndex)
1176 {
1177  Particle_t locPID = locKinematicData->PID();
1178  DVector3 locMomentum = locKinematicData->momentum();
1179  DVector3 locPosition = locKinematicData->position();
1180 
1181  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
1182  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
1183  double locP = locMomentum.Mag();
1184 
1185  double locBeta_Timing = 0.0;
1186  auto locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locKinematicData);
1187  if(locNeutralHypo != nullptr)
1188  locBeta_Timing = locNeutralHypo->measuredBeta();
1189  else
1190  {
1191  auto locChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locKinematicData);
1192  if(locChargedHypo != nullptr)
1193  locBeta_Timing = locChargedHypo->measuredBeta();
1194  }
1195  double locDeltaBeta = locBeta_Timing - locKinematicData->lorentzMomentum().Beta();
1196 
1197  //FILL HISTOGRAMS
1198  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1199  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1200  Lock_Action();
1201  {
1202  dHistDeque_P[locStepIndex][locPID][locIsMissingFlag]->Fill(locP);
1203  dHistDeque_Phi[locStepIndex][locPID][locIsMissingFlag]->Fill(locPhi);
1204  dHistDeque_Theta[locStepIndex][locPID][locIsMissingFlag]->Fill(locTheta);
1205  dHistDeque_PVsTheta[locStepIndex][locPID][locIsMissingFlag]->Fill(locTheta, locP);
1206  dHistDeque_PhiVsTheta[locStepIndex][locPID][locIsMissingFlag]->Fill(locTheta, locPhi);
1207  dHistDeque_BetaVsP[locStepIndex][locPID][locIsMissingFlag]->Fill(locP, locBeta_Timing);
1208  dHistDeque_DeltaBetaVsP[locStepIndex][locPID][locIsMissingFlag]->Fill(locP, locDeltaBeta);
1209  dHistDeque_VertexZ[locStepIndex][locPID][locIsMissingFlag]->Fill(locKinematicData->position().Z());
1210  dHistDeque_VertexYVsX[locStepIndex][locPID][locIsMissingFlag]->Fill(locKinematicData->position().X(), locKinematicData->position().Y());
1211  }
1212  Unlock_Action();
1213 }
1214 
1216 {
1217  DVector3 locMomentum = locKinematicData->momentum();
1218  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
1219  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
1220  double locP = locMomentum.Mag();
1221  double locDeltaTRF = locKinematicData->time() - (locEventRFBunch->dTime + (locKinematicData->z() - dTargetZCenter)/29.9792458);
1222 
1223  //FILL HISTOGRAMS
1224  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1225  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1226  Lock_Action();
1227  {
1228  dBeamParticleHist_P->Fill(locP);
1229  dBeamParticleHist_Phi->Fill(locPhi);
1230  dBeamParticleHist_Theta->Fill(locTheta);
1231  dBeamParticleHist_PVsTheta->Fill(locTheta, locP);
1232  dBeamParticleHist_PhiVsTheta->Fill(locTheta, locPhi);
1233  dBeamParticleHist_VertexZ->Fill(locKinematicData->position().Z());
1234  dBeamParticleHist_VertexYVsX->Fill(locKinematicData->position().X(), locKinematicData->position().Y());
1235  dBeamParticleHist_DeltaTRF->Fill(locDeltaTRF);
1236  dBeamParticleHist_DeltaTRFVsBeamE->Fill(locKinematicData->energy(), locDeltaTRF);
1237  }
1238  Unlock_Action();
1239 }
1240 
1241 void DHistogramAction_InvariantMass::Initialize(JEventLoop* locEventLoop)
1242 {
1243  string locHistName, locHistTitle;
1244  double locMassPerBin = 1000.0*(dMaxMass - dMinMass)/((double)dNumMassBins);
1245  locEventLoop->GetSingle(dAnalysisUtilities);
1246 
1247  string locParticleNamesForHist = "";
1248  if(dInitialPID != Unknown)
1249  {
1251  locParticleNamesForHist = DAnalysis::Convert_PIDsToROOTName(locChainPIDs);
1252  }
1253  else
1254  {
1255  for(size_t loc_i = 0; loc_i < dToIncludePIDs.size(); ++loc_i)
1256  locParticleNamesForHist += ParticleName_ROOT(dToIncludePIDs[loc_i]);
1257  }
1258 
1259  bool locBeamPresent = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());
1260 
1261  //CREATE THE HISTOGRAMS
1262  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1263  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1264  {
1266 
1267  locHistName = "InvariantMass";
1268  ostringstream locStream;
1269  locStream << locMassPerBin;
1270  locHistTitle = string(";") + locParticleNamesForHist + string(" Invariant Mass (GeV/c^{2});# Combos / ") + locStream.str() + string(" MeV/c^{2}");
1271  dHist_InvariantMass = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMassBins, dMinMass, dMaxMass);
1272 
1273  if(locBeamPresent)
1274  {
1275  locHistName = "InvariantMassVsBeamE";
1276  locHistTitle = string(";Beam Energy (GeV);") + locParticleNamesForHist + string(" Invariant Mass (GeV/c^{2})");
1277  dHist_InvariantMassVsBeamE = GetOrCreate_Histogram<TH2D>(locHistName, locHistTitle, dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DMassBins, dMinMass, dMaxMass);
1278  }
1279 
1280  //Return to the base directory
1282  }
1283  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1284 }
1285 
1286 bool DHistogramAction_InvariantMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1287 {
1288  auto locFirstStep = locParticleCombo->Get_ParticleComboStep(0);
1289  bool locBeamPresent = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());
1290  auto locBeam = Get_UseKinFitResultsFlag() ? locFirstStep->Get_InitialParticle() : locFirstStep->Get_InitialParticle_Measured();
1291 
1292  vector<double> locMassesToFill;
1293  vector<double> loc2DMassesToFill;
1294  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
1295  {
1296  const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
1297  if((dInitialPID != Unknown) && (locReactionStep->Get_InitialPID() != dInitialPID))
1298  continue;
1299  if((dStepIndex != -1) && (int(loc_i) != dStepIndex))
1300  continue;
1301 
1302  //build all possible combinations of the included pids
1303  set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dToIncludePIDs);
1304 
1305  //loop over them
1306  set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin();
1307  for(; locComboIterator != locIndexCombos.end(); ++locComboIterator)
1308  {
1309  set<pair<const JObject*, unsigned int> > locSourceObjects;
1310  DLorentzVector locFinalStateP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locComboIterator, locSourceObjects, Get_UseKinFitResultsFlag());
1311  if(dPreviousSourceObjects.find(locSourceObjects) == dPreviousSourceObjects.end())
1312  {
1313  dPreviousSourceObjects.insert(locSourceObjects);
1314  locMassesToFill.push_back(locFinalStateP4.M());
1315  }
1316 
1317  auto locBeamSourceObjects = std::make_pair(locSourceObjects, locBeam);
1318  if(locBeamPresent && (dPreviousSourceObjects_Beam.find(locBeamSourceObjects) == dPreviousSourceObjects_Beam.end()))
1319  {
1320  dPreviousSourceObjects_Beam.insert(locBeamSourceObjects);
1321  loc2DMassesToFill.push_back(locFinalStateP4.M());
1322  }
1323  }
1324  //don't break: e.g. if multiple pi0's, histogram invariant mass of each one
1325  }
1326 
1327  //FILL HISTOGRAMS
1328  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1329  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1330  Lock_Action();
1331  {
1332  for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)
1333  dHist_InvariantMass->Fill(locMassesToFill[loc_i]);
1334  for(size_t loc_i = 0; loc_i < loc2DMassesToFill.size(); ++loc_i)
1335  dHist_InvariantMassVsBeamE->Fill(locBeam->energy(), loc2DMassesToFill[loc_i]);
1336  }
1337  Unlock_Action();
1338 
1339  return true;
1340 }
1341 
1342 void DHistogramAction_MissingMass::Initialize(JEventLoop* locEventLoop)
1343 {
1344  double locMassPerBin = 1000.0*(dMaxMass - dMinMass)/((double)dNumMassBins);
1345  string locInitialParticlesROOTName = DAnalysis::Get_InitialParticlesName(Get_Reaction()->Get_ReactionStep(0), true);
1346  vector<Particle_t> locMissingMassOffOfPIDs(dMissingMassOffOfPIDs.begin(), dMissingMassOffOfPIDs.end());
1347  auto locChainPIDs = DAnalysis::Get_ChainPIDs(Get_Reaction(), Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID(), dMissingMassOffOfStepIndex, locMissingMassOffOfPIDs, !Get_UseKinFitResultsFlag(), true);
1348  string locFinalParticlesROOTName = DAnalysis::Convert_PIDsToROOTName(locChainPIDs);
1349 
1350  locEventLoop->GetSingle(dAnalysisUtilities);
1351 
1352  //CREATE THE HISTOGRAMS
1353  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1354  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1355  {
1357 
1358  string locHistName = "MissingMass";
1359  ostringstream locStream;
1360  locStream << locMassPerBin;
1361  string locHistTitle = string(";") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass (GeV/c^{2});# Combos / ") + locStream.str() + string(" MeV/c^{2}");
1362  dHist_MissingMass = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMassBins, dMinMass, dMaxMass);
1363 
1364  locHistName = "MissingMassVsBeamE";
1365  locMassPerBin *= ((double)dNumMassBins)/((double)dNum2DMassBins);
1366  locStream.str("");
1367  locStream << locMassPerBin;
1368  locHistTitle = string(";Beam Energy (GeV);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass (GeV/c^{2})");
1369  dHist_MissingMassVsBeamE = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DMassBins, dMinMass, dMaxMass);
1370 
1371  locHistName = "MissingMassVsMissingP";
1372  locStream.str("");
1373  locStream << locMassPerBin;
1374  locHistTitle = string(";Missing P (GeV/c);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass (GeV/c^{2})");
1375  dHist_MissingMassVsMissingP = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DMissPBins, dMinMissP, dMaxMissP, dNum2DMassBins, dMinMass, dMaxMass);
1376 
1377  //Return to the base directory
1379  }
1380  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1381 }
1382 
1383 bool DHistogramAction_MissingMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1384 {
1385  double locBeamEnergy = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->energy();
1386 
1387  //build all possible combinations of the included pids
1388  set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(Get_Reaction()->Get_ReactionStep(dMissingMassOffOfStepIndex), dMissingMassOffOfPIDs);
1389 
1390  //loop over them
1391  vector<pair<double, double> > locMassesToFill; //first is missing mass, second is missing p
1392  set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin();
1393  for(; locComboIterator != locIndexCombos.end(); ++locComboIterator)
1394  {
1395  set<pair<const JObject*, unsigned int> > locSourceObjects;
1396  DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo, 0, dMissingMassOffOfStepIndex, *locComboIterator, locSourceObjects, Get_UseKinFitResultsFlag());
1397 
1398  if(dPreviousSourceObjects.find(locSourceObjects) != dPreviousSourceObjects.end())
1399  continue; //dupe: already histed!
1400  dPreviousSourceObjects.insert(locSourceObjects);
1401 
1402  locMassesToFill.push_back(pair<double, double>(locMissingP4.M(), locMissingP4.P()));
1403  }
1404 
1405  //FILL HISTOGRAMS
1406  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1407  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1408  Lock_Action();
1409  {
1410  for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)
1411  {
1412  dHist_MissingMass->Fill(locMassesToFill[loc_i].first);
1413  dHist_MissingMassVsBeamE->Fill(locBeamEnergy, locMassesToFill[loc_i].first);
1414  dHist_MissingMassVsMissingP->Fill(locMassesToFill[loc_i].second, locMassesToFill[loc_i].first);
1415  }
1416  }
1417  Unlock_Action();
1418 
1419  return true;
1420 }
1421 
1423 {
1424  double locMassSqPerBin = 1000.0*1000.0*(dMaxMassSq - dMinMassSq)/((double)dNumMassBins);
1425  string locInitialParticlesROOTName = DAnalysis::Get_InitialParticlesName(Get_Reaction()->Get_ReactionStep(0), true);
1426  vector<Particle_t> locMissingMassOffOfPIDs(dMissingMassOffOfPIDs.begin(), dMissingMassOffOfPIDs.end());
1427  auto locChainPIDs = DAnalysis::Get_ChainPIDs(Get_Reaction(), Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID(), dMissingMassOffOfStepIndex, locMissingMassOffOfPIDs, !Get_UseKinFitResultsFlag(), true);
1428  string locFinalParticlesROOTName = DAnalysis::Convert_PIDsToROOTName(locChainPIDs);
1429 
1430  locEventLoop->GetSingle(dAnalysisUtilities);
1431 
1432  //CREATE THE HISTOGRAMS
1433  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1434  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1435  {
1437 
1438  string locHistName = "MissingMassSquared";
1439  ostringstream locStream;
1440  locStream << locMassSqPerBin;
1441  string locHistTitle = string(";") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass Squared (GeV/c^{2})^{2};# Combos / ") + locStream.str() + string(" (MeV/c^{2})^{2}");
1442  dHist_MissingMassSquared = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMassBins, dMinMassSq, dMaxMassSq);
1443 
1444  locHistName = "MissingMassSquaredVsBeamE";
1445  locMassSqPerBin *= ((double)dNumMassBins)/((double)dNum2DMassBins);
1446  locStream.str();
1447  locStream << locMassSqPerBin;
1448  locHistTitle = string(";Beam Energy (GeV);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass Squared (GeV/c^{2})^{2};");
1449  dHist_MissingMassSquaredVsBeamE = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DMassBins, dMinMassSq, dMaxMassSq);
1450 
1451  locHistName = "MissingMassSquaredVsMissingP";
1452  locStream.str("");
1453  locStream << locMassSqPerBin;
1454  locHistTitle = string(";Missing P (GeV/c);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass Squared (GeV/c^{2})^{2}");
1455  dHist_MissingMassSquaredVsMissingP = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DMissPBins, dMinMissP, dMaxMissP, dNum2DMassBins, dMinMassSq, dMaxMassSq);
1456 
1457  //Return to the base directory
1459  }
1460  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1461 }
1462 
1463 bool DHistogramAction_MissingMassSquared::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1464 {
1465  double locBeamEnergy = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->energy();
1466 
1467  //build all possible combinations of the included pids
1468  set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(Get_Reaction()->Get_ReactionStep(dMissingMassOffOfStepIndex), dMissingMassOffOfPIDs);
1469 
1470  //loop over them
1471  vector<pair<double, double> > locMassesToFill; //first is missing mass, second is missing p
1472  set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin();
1473  for(; locComboIterator != locIndexCombos.end(); ++locComboIterator)
1474  {
1475  set<pair<const JObject*, unsigned int> > locSourceObjects;
1476  DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo, 0, dMissingMassOffOfStepIndex, *locComboIterator, locSourceObjects, Get_UseKinFitResultsFlag());
1477 
1478  if(dPreviousSourceObjects.find(locSourceObjects) != dPreviousSourceObjects.end())
1479  continue; //dupe: already histed!
1480  dPreviousSourceObjects.insert(locSourceObjects);
1481 
1482  locMassesToFill.push_back(pair<double, double>(locMissingP4.M2(), locMissingP4.P()));
1483  }
1484 
1485  //FILL HISTOGRAMS
1486  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1487  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1488  Lock_Action();
1489  {
1490  for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)
1491  {
1492  dHist_MissingMassSquared->Fill(locMassesToFill[loc_i].first);
1493  dHist_MissingMassSquaredVsBeamE->Fill(locBeamEnergy, locMassesToFill[loc_i].first);
1494  dHist_MissingMassSquaredVsMissingP->Fill(locMassesToFill[loc_i].second, locMassesToFill[loc_i].first);
1495  }
1496  }
1497  Unlock_Action();
1498 
1499  return true;
1500 }
1501 
1502 void DHistogramAction_2DInvariantMass::Initialize(JEventLoop* locEventLoop)
1503 {
1504  string locHistName, locHistTitle;
1505  locEventLoop->GetSingle(dAnalysisUtilities);
1506 
1507  string locParticleNamesForHistX = "";
1508  for(size_t loc_i = 0; loc_i < dXPIDs.size(); ++loc_i)
1509  locParticleNamesForHistX += ParticleName_ROOT(dXPIDs[loc_i]);
1510 
1511  string locParticleNamesForHistY = "";
1512  for(size_t loc_i = 0; loc_i < dYPIDs.size(); ++loc_i)
1513  locParticleNamesForHistY += ParticleName_ROOT(dYPIDs[loc_i]);
1514 
1515  string locMassString = " Invariant Mass (GeV/c^{2});";
1516 
1517  //CREATE THE HISTOGRAMS
1518  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1519  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1520  {
1522 
1523  locHistName = "2DInvariantMass";
1524  locHistTitle = string(";") + locParticleNamesForHistX + locMassString + locParticleNamesForHistY + locMassString;
1525  dHist_2DInvaraintMass = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumXBins, dMinX, dMaxX, dNumYBins, dMinY, dMaxY);
1526 
1527  //Return to the base directory
1529  }
1530  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1531 }
1532 
1533 bool DHistogramAction_2DInvariantMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1534 {
1535  vector<pair<double, double> > locMassesToFill;
1536  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
1537  {
1538  const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
1539  if((dStepIndex != -1) && (int(loc_i) != dStepIndex))
1540  continue;
1541 
1542  //build all possible combinations of the included pids
1543  set<set<size_t> > locXIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dXPIDs);
1544  set<set<size_t> > locYIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dYPIDs);
1545 
1546  //(double) loop over them
1547  set<set<size_t> >::iterator locXComboIterator = locXIndexCombos.begin();
1548  for(; locXComboIterator != locXIndexCombos.end(); ++locXComboIterator)
1549  {
1550  set<pair<const JObject*, unsigned int> > locXSourceObjects;
1551  DLorentzVector locXP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locXComboIterator, locXSourceObjects, Get_UseKinFitResultsFlag());
1552 
1553  set<set<size_t> >::iterator locYComboIterator = locYIndexCombos.begin();
1554  for(; locYComboIterator != locYIndexCombos.end(); ++locYComboIterator)
1555  {
1556  set<pair<const JObject*, unsigned int> > locYSourceObjects;
1557  DLorentzVector locYP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locYComboIterator, locYSourceObjects, Get_UseKinFitResultsFlag());
1558 
1559  if(locXSourceObjects == locYSourceObjects)
1560  continue; //the same!
1561 
1562  set<set<pair<const JObject*, unsigned int> > > locAllSourceObjects;
1563  locAllSourceObjects.insert(locXSourceObjects);
1564  locAllSourceObjects.insert(locYSourceObjects);
1565  if(dPreviousSourceObjects.find(locAllSourceObjects) != dPreviousSourceObjects.end())
1566  continue; //dupe: already histed! (also, could be that the X/Y swapped combo has already been histed: don't double-count!
1567  dPreviousSourceObjects.insert(locAllSourceObjects);
1568 
1569  locMassesToFill.push_back(pair<double, double>(locXP4.M(), locYP4.M()));
1570  }
1571  }
1572  }
1573 
1574  //FILL HISTOGRAMS
1575  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1576  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1577  Lock_Action();
1578  {
1579  for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)
1580  dHist_2DInvaraintMass->Fill(locMassesToFill[loc_i].first, locMassesToFill[loc_i].second);
1581  }
1582  Unlock_Action();
1583 
1584  return true;
1585 }
1586 
1587 void DHistogramAction_Dalitz::Initialize(JEventLoop* locEventLoop)
1588 {
1589  string locHistName, locHistTitle;
1590  locEventLoop->GetSingle(dAnalysisUtilities);
1591 
1592  string locParticleNamesForHistX = "";
1593  for(size_t loc_i = 0; loc_i < dXPIDs.size(); ++loc_i)
1594  locParticleNamesForHistX += ParticleName_ROOT(dXPIDs[loc_i]);
1595 
1596  string locParticleNamesForHistY = "";
1597  for(size_t loc_i = 0; loc_i < dYPIDs.size(); ++loc_i)
1598  locParticleNamesForHistY += ParticleName_ROOT(dYPIDs[loc_i]);
1599 
1600  string locMassString = " Invariant Mass Squared (GeV/c^{2})^{2};";
1601 
1602  //CREATE THE HISTOGRAMS
1603  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1604  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1605  {
1607 
1608  locHistName = "DalitzPlot";
1609  locHistTitle = string(";") + locParticleNamesForHistX + locMassString + locParticleNamesForHistY + locMassString;
1610  dHist_DalitzPlot = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumXBins, dMinX, dMaxX, dNumYBins, dMinY, dMaxY);
1611 
1612  //Return to the base directory
1614  }
1615  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1616 }
1617 
1618 bool DHistogramAction_Dalitz::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1619 {
1620  vector<pair<double, double> > locMassesToFill;
1621  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
1622  {
1623  const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
1624  if((dStepIndex != -1) && (int(loc_i) != dStepIndex))
1625  continue;
1626 
1627  //build all possible combinations of the included pids
1628  set<set<size_t> > locXIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dXPIDs);
1629  set<set<size_t> > locYIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dYPIDs);
1630 
1631  //(double) loop over them
1632  set<set<size_t> >::iterator locXComboIterator = locXIndexCombos.begin();
1633  for(; locXComboIterator != locXIndexCombos.end(); ++locXComboIterator)
1634  {
1635  set<pair<const JObject*, unsigned int> > locXSourceObjects;
1636  DLorentzVector locXP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locXComboIterator, locXSourceObjects, Get_UseKinFitResultsFlag());
1637 
1638  set<set<size_t> >::iterator locYComboIterator = locYIndexCombos.begin();
1639  for(; locYComboIterator != locYIndexCombos.end(); ++locYComboIterator)
1640  {
1641  set<pair<const JObject*, unsigned int> > locYSourceObjects;
1642  DLorentzVector locYP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locYComboIterator, locYSourceObjects, Get_UseKinFitResultsFlag());
1643 
1644  if(locXSourceObjects == locYSourceObjects)
1645  continue; //the same!
1646 
1647  set<set<pair<const JObject*, unsigned int> > > locAllSourceObjects;
1648  locAllSourceObjects.insert(locXSourceObjects);
1649  locAllSourceObjects.insert(locYSourceObjects);
1650  if(dPreviousSourceObjects.find(locAllSourceObjects) != dPreviousSourceObjects.end())
1651  continue; //dupe: already histed! (also, could be that the X/Y swapped combo has already been histed: don't double-count!
1652  dPreviousSourceObjects.insert(locAllSourceObjects);
1653 
1654  locMassesToFill.push_back(pair<double, double>(locXP4.M2(), locYP4.M2()));
1655  }
1656  }
1657  }
1658 
1659  //FILL HISTOGRAMS
1660  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1661  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1662  Lock_Action();
1663  {
1664  for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)
1665  dHist_DalitzPlot->Fill(locMassesToFill[loc_i].first, locMassesToFill[loc_i].second);
1666  }
1667  Unlock_Action();
1668 
1669  return true;
1670 }
1671 
1672 void DHistogramAction_KinFitResults::Initialize(JEventLoop* locEventLoop)
1673 {
1674  auto locReaction = Get_Reaction();
1675  DKinFitType locKinFitType = locReaction->Get_KinFitType();
1676  if(locKinFitType == d_NoFit)
1677  return;
1678 
1679  locEventLoop->GetSingle(dAnalysisUtilities);
1680 
1681  //Get the DReactionVertexInfo for this reaction
1682  vector<const DReactionVertexInfo*> locReactionVertexInfos;
1683  locEventLoop->Get(locReactionVertexInfos);
1684  const DReactionVertexInfo* locReactionVertexInfo = nullptr;
1685  for(auto& locVertexInfo : locReactionVertexInfos)
1686  {
1687  auto locReactions = locVertexInfo->Get_Reactions();
1688  std::sort(locReactions.begin(), locReactions.end());
1689  if(!std::binary_search(locReactions.begin(), locReactions.end(), Get_Reaction()))
1690  continue;
1691  locReactionVertexInfo = locVertexInfo;
1692  break;
1693  }
1694 
1695  dKinFitUtils = new DKinFitUtils_GlueX(locEventLoop);
1696 
1697  size_t locNumConstraints = 0, locNumUnknowns = 0;
1698  string locConstraintString = dKinFitUtils->Get_ConstraintInfo(locReactionVertexInfo, Get_Reaction(), locNumConstraints, locNumUnknowns);
1699 
1700  size_t locNDF = locNumConstraints - locNumUnknowns;
1701  bool locIncludeBeamlineInVertexFitFlag = dKinFitUtils->Get_IncludeBeamlineInVertexFitFlag();
1702 
1703  bool locP4IsFit = ((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit));
1704  bool locSpactimeIsFitFlag = (locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit);
1705  bool locVertexIsFitFlag = (locSpactimeIsFitFlag || (locKinFitType == d_VertexFit) || (locKinFitType == d_P4AndVertexFit));
1706 
1707  //bool locIsInclusiveChannelFlag = Get_Reaction()->Get_IsInclusiveChannelFlag();
1708  //Below, should in theory check on whether to create pxyz pull histograms in the inclusive channel case
1709  //But, this is tricky: can have inclusive (no p4) but still have mass constraints (p4)
1710 
1711  //CREATE THE HISTOGRAMS
1712  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1713  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1714  {
1716 
1717  // Confidence Level
1718  string locHistName = "ConfidenceLevel";
1719  ostringstream locHistTitle;
1720  locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";Confidence Level (" << locNumConstraints;
1721  locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit);# Combos";
1722  dHist_ConfidenceLevel = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle.str(), dNumConfidenceLevelBins, 0.0, 1.0);
1723 
1724  //beam (pulls & conlev)
1725  Particle_t locInitialPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID();
1726  bool locBeamFlag = Get_IsFirstStepBeam(Get_Reaction());
1727  if(locBeamFlag)
1728  {
1729  auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(0);
1730  auto locFullConstrainParticles = locStepVertexInfo->Get_FullConstrainParticles(true, d_InitialState, d_AllCharges, false);
1731  bool locIsInVertexFitFlag = locVertexIsFitFlag && locIncludeBeamlineInVertexFitFlag && !locFullConstrainParticles.empty();
1732  bool locIsChargedFlag = (ParticleCharge(locInitialPID) != 0);
1733 
1734  if(locP4IsFit || locIsInVertexFitFlag)
1735  {
1736  string locFullROOTName = string("Beam ") + ParticleName_ROOT(locInitialPID);
1737 
1739  {
1740  //Conlev
1741  locHistName = "ConfidenceLevel_VsBeamEnergy";
1742  locHistTitle.str("");
1743  locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";Beam Energy (GeV);Confidence Level (" << locNumConstraints;
1744  locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";
1745  dHistMap_ConfidenceLevel_VsP[pair<int, Particle_t>(-1, locInitialPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DConfidenceLevelBins, 0.0, 1.0);
1746  }
1747 
1748  //Pulls
1749  CreateAndChangeTo_Directory("Beam", "Beam");
1750  map<DKinFitPullType, TH1I*> locParticlePulls;
1751  map<DKinFitPullType, TH2I*> locParticlePullsVsP, locParticlePullsVsTheta, locParticlePullsVsPhi;
1752  Create_ParticlePulls(locFullROOTName, locIsChargedFlag, locIsInVertexFitFlag, false, -1, locInitialPID);
1753 
1754  gDirectory->cd("..");
1755  }
1756  }
1757 
1758  //final particles
1759  for(size_t loc_i = 0; loc_i < Get_Reaction()->Get_NumReactionSteps(); ++loc_i)
1760  {
1761  const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);
1762  auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(loc_i);
1763  auto locFullConstrainParticles = locStepVertexInfo->Get_FullConstrainParticles(true, d_FinalState, d_AllCharges, false);
1764  ostringstream locStepName;
1765  locStepName << "Step" << loc_i << "__" << DAnalysis::Get_StepName(locReactionStep, true, false);
1766  string locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);
1767 
1768  //Conlev
1769  auto locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_AllCharges, false);
1770  for(auto locPID : locPIDs)
1771  {
1772  bool locIsInVertexFitFlag = locVertexIsFitFlag && locStepVertexInfo->Get_FittableVertexFlag();
1773 
1774  bool locIsNeutralShowerFlag = (locIsInVertexFitFlag && (ParticleCharge(locPID) == 0));
1775  if((ParticleMass(locPID) > 0.0) && !locSpactimeIsFitFlag)
1776  locIsNeutralShowerFlag = false; //massive shower momentum is defined by t, which isn't fit: use particle
1777  if(!locP4IsFit && !locIsInVertexFitFlag)
1778  continue; //p4 is not fit, and this is not in a vertex fit: no pulls
1779  if(locIsNeutralShowerFlag && !locP4IsFit)
1780  continue; //vertex-only fit: neutral shower does not constrain: no pulls
1781 
1782  string locParticleName = ParticleType(locPID);
1783  string locParticleROOTName = ParticleName_ROOT(locPID);
1784 
1786  {
1787  //vs p
1788  locHistName = string("ConfidenceLevel_Vs") + locParticleName + string("P");
1789  locHistTitle.str("");
1790  locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";" << locParticleROOTName << " p (GeV/c);Confidence Level (" << locNumConstraints;
1791  locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";
1792  dHistMap_ConfidenceLevel_VsP[pair<int, Particle_t>(loc_i, locPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DPBins, dMinP, dMaxP, dNum2DConfidenceLevelBins, 0.0, 1.0);
1793 
1794  //vs theta
1795  locHistName = string("ConfidenceLevel_Vs") + locParticleName + string("Theta");
1796  locHistTitle.str("");
1797  locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";" << locParticleROOTName << " #theta#circ;Confidence Level (" << locNumConstraints;
1798  locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";
1799  dHistMap_ConfidenceLevel_VsTheta[pair<int, Particle_t>(loc_i, locPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DConfidenceLevelBins, 0.0, 1.0);
1800 
1801  //vs phi
1802  locHistName = string("ConfidenceLevel_Vs") + locParticleName + string("Phi");
1803  locHistTitle.str("");
1804  locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";" << locParticleROOTName << " #phi#circ;Confidence Level (" << locNumConstraints;
1805  locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";
1806  dHistMap_ConfidenceLevel_VsPhi[pair<int, Particle_t>(loc_i, locPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DConfidenceLevelBins, 0.0, 1.0);
1807  }
1808  }
1809 
1810  //pulls
1811  bool locFilledFlag = false;
1812  for(auto locPID : locPIDs)
1813  {
1814  bool locIsInVertexFitFlag = locVertexIsFitFlag && locStepVertexInfo->Get_FittableVertexFlag();
1815  bool locIsChargedFlag = (ParticleCharge(locPID) != 0);
1816 
1817  bool locIsNeutralShowerFlag = (locIsInVertexFitFlag && (ParticleCharge(locPID) == 0));
1818  if((ParticleMass(locPID) > 0.0) && !locSpactimeIsFitFlag)
1819  locIsNeutralShowerFlag = false; //massive shower momentum is defined by t, which isn't fit: use particle
1820  if(!locP4IsFit && !locIsInVertexFitFlag)
1821  continue; //p4 is not fit, and this is not in a vertex fit: no pulls
1822  if(locIsNeutralShowerFlag && !locP4IsFit)
1823  continue; //vertex-only fit: neutral shower does not constrain: no pulls
1824 
1825  if(!locFilledFlag) //first call
1826  {
1827  locFilledFlag = true;
1828  CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());
1829  }
1830 
1831  string locParticleName = ParticleType(locPID);
1832  CreateAndChangeTo_Directory(locParticleName, locParticleName);
1833 
1834  string locParticleROOTName = ParticleName_ROOT(locPID);
1835  string locFullROOTName = locParticleROOTName + string(", ") + locStepROOTName;
1836  Create_ParticlePulls(locFullROOTName, locIsChargedFlag, locIsInVertexFitFlag, locIsNeutralShowerFlag, loc_i, locPID);
1837 
1838  gDirectory->cd("..");
1839  } //end of particle loop
1840 
1841  if(locFilledFlag)
1842  gDirectory->cd("..");
1843  } //end of step loop
1844 
1845  //Return to the base directory
1847  }
1848  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1849 }
1850 
1851 void DHistogramAction_KinFitResults::Create_ParticlePulls(string locFullROOTName, bool locIsChargedFlag, bool locIsInVertexFitFlag, bool locIsNeutralShowerFlag, int locStepIndex, Particle_t locPID)
1852 {
1853  auto locParticlePair = std::make_pair(locStepIndex, locPID);
1854 
1855  DKinFitType locKinFitType = Get_Reaction()->Get_KinFitType();
1856  bool locP4IsFit = ((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit));
1857 
1858  //Determine pull types
1859  set<pair<DKinFitPullType, pair<string, string> > > locPullTypes;
1860 
1861  //p4 pulls:
1862  string locHistName, locHistTitle;
1863  if(locIsNeutralShowerFlag && locP4IsFit) //neutral shower not in a p4-only fit
1864  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_EPull, pair<string, string>("Pull_E", "E Pull")));
1865  else if(!locIsNeutralShowerFlag && (locP4IsFit || locIsInVertexFitFlag))
1866  {
1867  //all detected particles (except neutral showers) have p3 used in p4 fits and vertex fits
1868  //however, don't include if the particles aren't actually in the fit (vertex-only, and too few particles at that vertex to constrain)
1869  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_PxPull, pair<string, string>("Pull_Px", "p_{x} Pull")));
1870  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_PyPull, pair<string, string>("Pull_Py", "p_{y} Pull")));
1871  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_PzPull, pair<string, string>("Pull_Pz", "p_{z} Pull")));
1872  }
1873 
1874  //vertex pulls:
1875  if((locIsNeutralShowerFlag && locP4IsFit) || (locIsChargedFlag && locIsInVertexFitFlag))
1876  {
1877  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_XxPull, pair<string, string>("Pull_Xx", "x_{x} Pull")));
1878  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_XyPull, pair<string, string>("Pull_Xy", "x_{y} Pull")));
1879  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_XzPull, pair<string, string>("Pull_Xz", "x_{z} Pull")));
1880  }
1881 
1882  //time pulls:
1883  if(locIsInVertexFitFlag && ((locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit)))
1884  locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_TPull, pair<string, string>("Pull_T", "t Pull")));
1885 
1886  bool locIsBeamFlag = (locFullROOTName.substr(0, 4) == "Beam");
1887  int locNum2DPBins = locIsBeamFlag ? dNum2DBeamEBins : dNum2DPBins;
1888  double locMinP = locIsBeamFlag ? dMinBeamE : dMinP;
1889  double locMaxP = locIsBeamFlag ? dMaxBeamE : dMaxP;
1890 
1891  //Create histograms
1892  set<pair<DKinFitPullType, pair<string, string> > >::iterator locIterator = locPullTypes.begin();
1893  for(; locIterator != locPullTypes.end(); ++locIterator)
1894  {
1895  pair<string, string> locStrings = (*locIterator).second;
1896  auto locPullType = (*locIterator).first;
1897 
1898  //1D
1899  locHistName = locStrings.first;
1900  locHistTitle = locFullROOTName + string(";") + locStrings.second + string(";# Combos");
1901  dHistMap_Pulls[locParticlePair][locPullType] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
1902 
1903  int locNumDeltaBins;
1904  double locMaxDeltaRange;
1905  Get_DeltaBinningParams(locPullType, false, locNumDeltaBins, locMaxDeltaRange);
1906  Get_HistNameTitle(locPullType, locFullROOTName, 0, locHistName, locHistTitle);
1907  dHistMap_Deltas[locParticlePair][locPullType] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);
1908 
1909  if(!dHistDependenceFlag)
1910  continue;
1911  Get_DeltaBinningParams(locPullType, true, locNumDeltaBins, locMaxDeltaRange);
1912 
1913  //vs p
1914  locHistName = locStrings.first + string("_VsP");
1915  locHistTitle = locFullROOTName + string(";p (GeV/c);") + locStrings.second;
1916  dHistMap_PullsVsP[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, locNum2DPBins, locMinP, locMaxP, dNumPullBins, dMinPull, dMaxPull);
1917  Get_HistNameTitle(locPullType, locFullROOTName, 1, locHistName, locHistTitle);
1918  dHistMap_DeltasVsP[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, locNum2DPBins, locMinP, locMaxP, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);
1919  if(locIsBeamFlag)
1920  continue;
1921 
1922  //vs theta
1923  locHistName = locStrings.first + string("_VsTheta");
1924  locHistTitle = locFullROOTName + string(";#theta#circ;") + locStrings.second;
1925  dHistMap_PullsVsTheta[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumPullBins, dMinPull, dMaxPull);
1926  Get_HistNameTitle(locPullType, locFullROOTName, 2, locHistName, locHistTitle);
1927  dHistMap_DeltasVsTheta[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);
1928 
1929  //vs phi
1930  locHistName = locStrings.first + string("_VsPhi");
1931  locHistTitle = locFullROOTName + string(";#phi#circ;") + locStrings.second;
1932  dHistMap_PullsVsPhi[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNumPullBins, dMinPull, dMaxPull);
1933  Get_HistNameTitle(locPullType, locFullROOTName, 3, locHistName, locHistTitle);
1934  dHistMap_DeltasVsPhi[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);
1935  }
1936 }
1937 
1938 void DHistogramAction_KinFitResults::Get_HistNameTitle(DKinFitPullType locPullType, string locFullROOTName, int locVsKey, string& locHistName, string& locHistTitle)
1939 {
1940  string locNameParam, locTitleParam;
1941  if(locPullType == d_EPull)
1942  {
1943  locNameParam = "ShowerE";
1944  locTitleParam = "E";
1945  }
1946  else if(locPullType == d_PxPull) //delta-theta!!
1947  {
1948  locNameParam = "Theta";
1949  locTitleParam = "#theta#circ";
1950  }
1951  else if(locPullType == d_PyPull) //delta-phi!!
1952  {
1953  locNameParam = "Phi";
1954  locTitleParam = "#phi#circ";
1955  }
1956  else if(locPullType == d_PzPull) //delta-p/p!!
1957  {
1958  locNameParam = "POverP";
1959  locTitleParam = "p/p";
1960  }
1961  else if(locPullType == d_XxPull)
1962  {
1963  locNameParam = "X";
1964  locTitleParam = "Vertex-X (cm)";
1965  }
1966  else if(locPullType == d_XyPull)
1967  {
1968  locNameParam = "Y";
1969  locTitleParam = "Vertex-Y (cm)";
1970  }
1971  else if(locPullType == d_XzPull)
1972  {
1973  locNameParam = "Z";
1974  locTitleParam = "Vertex-Z (cm)";
1975  }
1976  else if(locPullType == d_TPull)
1977  {
1978  locNameParam = "T";
1979  locTitleParam = "T (ns)";
1980  }
1981 
1982  if(locVsKey == 0) //1D
1983  {
1984  locHistName = string("Delta") + locNameParam;
1985  locHistTitle = locFullROOTName + string(";#Delta") + locTitleParam;
1986  }
1987  else if(locVsKey == 1) //p
1988  {
1989  locHistName = string("Delta") + locNameParam + string("_VsP");
1990  locHistTitle = locFullROOTName + string(";p (GeV/c);#Delta") + locTitleParam;
1991  }
1992  else if(locVsKey == 2) //theta
1993  {
1994  locHistName = string("Delta") + locNameParam + string("_VsTheta");
1995  locHistTitle = locFullROOTName + string(";#theta#circ;#Delta") + locTitleParam;
1996  }
1997  else if(locVsKey == 3) //phi
1998  {
1999  locHistName = string("Delta") + locNameParam + string("_VsPhi");
2000  locHistTitle = locFullROOTName + string(";#phi#circ;#Delta") + locTitleParam;
2001  }
2002 }
2003 
2004 void DHistogramAction_KinFitResults::Get_DeltaBinningParams(DKinFitPullType locPullType, bool loc2DFlag, int& locNumBins, double& locMax)
2005 {
2006  if(locPullType == d_EPull)
2007  {
2008  locNumBins = loc2DFlag ? dNum2DDeltaShowerEBins : dNumDeltaShowerEBins;
2009  locMax = dMaxDeltaShowerE;
2010  }
2011  else if(locPullType == d_PxPull) //delta-theta!!
2012  {
2013  locNumBins = loc2DFlag ? dNum2DDeltaThetaBins : dNumDeltaThetaBins;
2014  locMax = dMaxDeltaTheta;
2015  }
2016  else if(locPullType == d_PyPull) //delta-phi!!
2017  {
2018  locNumBins = loc2DFlag ? dNum2DDeltaPhiBins : dNumDeltaPhiBins;
2019  locMax = dMaxDeltaPhi;
2020  }
2021  else if(locPullType == d_PzPull) //delta-p/p!!
2022  {
2023  locNumBins = loc2DFlag ? dNum2DDeltaPOverPBins : dNumDeltaPOverPBins;
2024  locMax = dMaxDeltaPOverP;
2025  }
2026  else if((locPullType == d_XxPull) || (locPullType == d_XyPull))
2027  {
2028  locNumBins = loc2DFlag ? dNum2DDeltaVertexXYBins : dNumDeltaVertexXYBins;
2029  locMax = dMaxDeltaVertexXY;
2030  }
2031  else if(locPullType == d_XzPull)
2032  {
2033  locNumBins = loc2DFlag ? dNum2DDeltaVertexZBins : dNumDeltaVertexZBins;
2034  locMax = dMaxDeltaVertexZ;
2035  }
2036  else if(locPullType == d_TPull)
2037  {
2038  locNumBins = loc2DFlag ? dNum2DDeltaTBins : dNumDeltaTBins;
2039  locMax = dMaxDeltaT;
2040  }
2041 }
2042 
2043 bool DHistogramAction_KinFitResults::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2044 {
2045  //kinfit results are unique for each DParticleCombo: no need to check for duplicates
2046  const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults();
2047  if(locKinFitResults == NULL)
2048  return true;
2049 
2050  // Get Confidence Level
2051  double locConfidenceLevel = locKinFitResults->Get_ConfidenceLevel();
2052 
2053  // Get Pulls
2054  map<const JObject*, map<DKinFitPullType, double> > locPulls; //DKinematicData is the MEASURED particle
2055  locKinFitResults->Get_Pulls(locPulls);
2056 
2058  bool locBeamFlag = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());
2059 
2060  //FILL HISTOGRAMS
2061  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2062  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2063  Lock_Action();
2064  {
2065  dHist_ConfidenceLevel->Fill(locConfidenceLevel);
2066 
2067  // beam
2068  if(locBeamFlag)
2069  {
2070  const DKinematicData* locKinematicData = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured();
2071  map<DKinFitPullType, double> locParticlePulls = locPulls[locKinematicData];
2072  if(!locParticlePulls.empty())
2073  {
2074  pair<int, Particle_t> locParticlePair(-1, locInitPID);
2075  DVector3 locMomentum = locKinematicData->momentum();
2076  double locP = locMomentum.Mag();
2077 
2078  //con lev
2080  dHistMap_ConfidenceLevel_VsP[locParticlePair]->Fill(locP, locConfidenceLevel);
2081 
2082  //pulls
2083  if(locConfidenceLevel >= dPullHistConfidenceLevelCut)
2084  {
2085  map<DKinFitPullType, double>::iterator locIterator = locParticlePulls.begin();
2086  for(; locIterator != locParticlePulls.end(); ++locIterator)
2087  {
2088  dHistMap_Pulls[locParticlePair][locIterator->first]->Fill(locIterator->second);
2090  dHistMap_PullsVsP[locParticlePair][locIterator->first]->Fill(locP, locIterator->second);
2091  }
2092  }
2093 
2094  //deltas
2095  double locDeltaPOverP = (locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->pmag() - locP)/locP;
2096  dHistMap_Deltas[locParticlePair][d_PzPull]->Fill(locDeltaPOverP);
2098  dHistMap_DeltasVsP[locParticlePair][d_PzPull]->Fill(locP, locDeltaPOverP);
2099  }
2100  }
2101 
2102  // final particle pulls
2103  for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)
2104  {
2105  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);
2106  auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i));
2107  auto locKinFitParticles = locParticleComboStep->Get_FinalParticles(Get_Reaction()->Get_ReactionStep(loc_i), false, false);
2108  for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)
2109  {
2110  auto& locParticle = locParticles[loc_j];
2111  auto& locKinFitParticle = locKinFitParticles[loc_j];
2112 
2113  //get pulls for this particle
2114  map<DKinFitPullType, double> locParticlePulls;
2115  map<const JObject*, map<DKinFitPullType, double> >::iterator locParticleIterator = locPulls.find(locParticle);
2116  if(locParticleIterator != locPulls.end())
2117  locParticlePulls = locParticleIterator->second;
2118  else //is neutral shower
2119  {
2120  auto locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locParticle);
2121  locParticlePulls = locPulls[locNeutralHypo->Get_NeutralShower()];
2122  }
2123  if(locParticlePulls.empty())
2124  continue;
2125 
2126  pair<int, Particle_t> locParticlePair(loc_i, locParticle->PID());
2127  DVector3 locMomentum = locParticle->momentum();
2128  double locP = locMomentum.Mag();
2129  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2130  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2131 
2132  //deltas //loop used to easily figure out what params changed
2133  for(auto locIterator = locParticlePulls.begin(); locIterator != locParticlePulls.end(); ++locIterator)
2134  {
2135  auto locPullType = locIterator->first;
2136  auto locDelta = Get_Delta(locPullType, locParticle, locKinFitParticle);
2137  dHistMap_Deltas[locParticlePair][locPullType]->Fill(locDelta);
2138  if(!dHistDependenceFlag)
2139  continue;
2140 
2141  dHistMap_DeltasVsP[locParticlePair][locPullType]->Fill(locP, locDelta);
2142  dHistMap_DeltasVsTheta[locParticlePair][locPullType]->Fill(locTheta, locDelta);
2143  dHistMap_DeltasVsPhi[locParticlePair][locPullType]->Fill(locPhi, locDelta);
2144  }
2145 
2146  //con lev
2148  {
2149  dHistMap_ConfidenceLevel_VsP[locParticlePair]->Fill(locP, locConfidenceLevel);
2150  dHistMap_ConfidenceLevel_VsTheta[locParticlePair]->Fill(locTheta, locConfidenceLevel);
2151  dHistMap_ConfidenceLevel_VsPhi[locParticlePair]->Fill(locPhi, locConfidenceLevel);
2152  }
2153  if(locConfidenceLevel < dPullHistConfidenceLevelCut)
2154  continue;
2155 
2156  //fill histograms
2157  map<DKinFitPullType, double>::iterator locIterator = locParticlePulls.begin();
2158  for(; locIterator != locParticlePulls.end(); ++locIterator)
2159  {
2160  dHistMap_Pulls[locParticlePair][locIterator->first]->Fill(locIterator->second);
2161  if(!dHistDependenceFlag)
2162  continue;
2163  dHistMap_PullsVsP[locParticlePair][locIterator->first]->Fill(locP, locIterator->second);
2164  dHistMap_PullsVsTheta[locParticlePair][locIterator->first]->Fill(locTheta, locIterator->second);
2165  dHistMap_PullsVsPhi[locParticlePair][locIterator->first]->Fill(locPhi, locIterator->second);
2166  }
2167  }
2168  }
2169  }
2170  Unlock_Action();
2171 
2172  return true;
2173 }
2174 
2175 double DHistogramAction_KinFitResults::Get_Delta(DKinFitPullType locPullType, const DKinematicData* locMeasured, const DKinematicData* locKinFit)
2176 {
2177  //kinfit - measured
2178  if(locPullType == d_EPull)
2179  return locKinFit->energy() - locMeasured->energy();
2180  else if(locPullType == d_PxPull) //delta-theta!!
2181  return (locKinFit->momentum().Theta() - locMeasured->momentum().Theta())*180.0/TMath::Pi();
2182  else if(locPullType == d_PyPull) //delta-phi!!
2183  {
2184  auto locDeltaPhi = (locKinFit->momentum().Phi() - locMeasured->momentum().Phi())*180.0/TMath::Pi();
2185  while(locDeltaPhi > 180.0)
2186  locDeltaPhi -= 360.0;
2187  while(locDeltaPhi < -180.0)
2188  locDeltaPhi += 360.0;
2189  return locDeltaPhi;
2190  }
2191  else if(locPullType == d_PzPull) //delta-p-over-p!!
2192  {
2193  auto locP = locMeasured->momentum().Mag();
2194  return (locKinFit->momentum().Mag() - locP)/locP;
2195  }
2196  else if(locPullType == d_XxPull)
2197  return locKinFit->position().X() - locMeasured->position().X();
2198  else if(locPullType == d_XyPull)
2199  return locKinFit->position().Y() - locMeasured->position().Y();
2200  else if(locPullType == d_XzPull)
2201  return locKinFit->position().Z() - locMeasured->position().Z();
2202  else if(locPullType == d_TPull)
2203  return locKinFit->time() - locMeasured->time();
2204  else
2205  return 0.0;
2206 }
2207 
2209 {
2210  string locHistName, locHistTitle;
2211  double locPtPerBin = 1000.0*(dMaxPt - dMinPt)/((double)dNumPtBins);
2212 
2213  vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;
2214  locEventLoop->Get(locAnalysisUtilitiesVector);
2215 
2216  //CREATE THE HISTOGRAMS
2217  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2218  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2219  {
2220  dAnalysisUtilities = locAnalysisUtilitiesVector[0];
2222 
2223  locHistName = "MissingTransverseMomentum";
2224  ostringstream locStream;
2225  locStream << locPtPerBin;
2226  locHistTitle = string(";") + string(" Missing Transverse Momentum (GeV/c);# Combos / ") + locStream.str() + string(" MeV/c");
2227  dHist_MissingTransverseMomentum = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPtBins, dMinPt, dMaxPt);
2228 
2229  //Return to the base directory
2231  }
2232  japp->RootUnLock(); //RELEASE ROOT LOCK!!
2233 }
2234 
2235 bool DHistogramAction_MissingTransverseMomentum::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2236 {
2237  set<pair<const JObject*, unsigned int> > locSourceObjects;
2238  DLorentzVector locFinalStateP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, 0, locSourceObjects, Get_UseKinFitResultsFlag()); // Use step '0'
2239 
2240  if(dPreviousSourceObjects.find(locSourceObjects) != dPreviousSourceObjects.end())
2241  return true; //dupe: already histed!
2242  dPreviousSourceObjects.insert(locSourceObjects);
2243 
2244  double locMissingTransverseMomentum = locFinalStateP4.Pt();
2245 
2246  //FILL HISTOGRAMS
2247  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2248  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2249  Lock_Action();
2250  {
2251  dHist_MissingTransverseMomentum->Fill(locMissingTransverseMomentum);
2252  }
2253  Unlock_Action();
2254 
2255  return true;
2256 }
2257 
vector< const DKinematicData * > Get_FinalParticles_Measured(void) const
DetectorSystem_t t1_detector(void) const
vector< pair< int, int > > Get_FullConstrainParticles(bool locFitFlag, DReactionState_t locState=d_EitherState, Charge_t locCharge=d_AllCharges, bool locIncludeDecayingFlag=true) const
TDirectoryFile * ChangeTo_BaseDirectory(void)
set< pair< set< pair< const JObject *, unsigned int > >, const JObject * > > dPreviousSourceObjects_Beam
bool Get_IsFirstStepBeam(const DReaction *locReaction)
Definition: DReaction.h:178
void Initialize(JEventLoop *locEventLoop)
double getEnergy() const
Definition: DFCALShower.h:156
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
DKinFitPullType
const DKinematicData * Get_FinalParticle_Measured(size_t locFinalParticleIndex) const
deque< map< Particle_t, map< bool, TH2I * > > > dHistDeque_DeltaBetaVsP
double Get_ChiSq_Timing(void) const
const DAnalysisUtilities * dAnalysisUtilities
static char * ParticleName_ROOT(Particle_t p)
Definition: particleType.h:851
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_dEdXFOMVsP
bool Get_IncludeBeamlineInVertexFitFlag(void) const
pair< int, int > Get_InitialParticleDecayFromIndices(const DReaction *locReaction, int locStepIndex)
Definition: DReaction.cc:93
void Fill_ChargedHists(const DChargedTrackHypothesis *locChargedTrackHypothesis, const DMCThrownMatching *locMCThrownMatching, const DEventRFBunch *locEventRFBunch)
size_t Get_NumParticleComboSteps(void) const
void Create_ParticlePulls(string locFullROOTName, bool locIsChargedFlag, bool locIsInVertexFitFlag, bool locIsNeutralShowerFlag, int locStepIndex, Particle_t locPID)
DKinFitType Get_KinFitType(void) const
Definition: DReaction.h:79
deque< map< Particle_t, map< bool, TH2I * > > > dHistDeque_VertexYVsX
const DAnalysisUtilities * dAnalysisUtilities
double pmag(void) const
double z(void) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
TVector3 DVector3
Definition: DVector3.h:14
string Convert_PIDsToROOTName(const vector< Particle_t > &locPIDs)
Definition: DReaction.h:230
double energy(void) const
const DReaction * Get_Reaction(void) const
deque< map< Particle_t, map< bool, TH2I * > > > dHistDeque_PhiVsTheta
const DAnalysisUtilities * dAnalysisUtilities
const DReactionStepVertexInfo * Get_StepVertexInfo(size_t locStepIndex) const
char string[256]
axes Fill(100, 100)
TDirectoryFile * CreateAndChangeTo_Directory(TDirectoryFile *locBaseDirectory, string locDirName, string locDirTitle)
double Get_TimeAtPOCAToVertex(void) const
bool Get(string xpath, string &sval) const
Definition: DGeometry.h:50
const DVector3 & position(void) const
map< pair< int, Particle_t >, TH2I * > dHistMap_ConfidenceLevel_VsPhi
set< set< pair< const JObject *, unsigned int > > > dPreviousSourceObjects
const DTrackTimeBased * Get_TrackTimeBased(void) const
void Fill_BeamHists(const DKinematicData *locKinematicData, const DEventRFBunch *locEventRFBunch)
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
TDirectoryFile * CreateAndChangeTo_ActionDirectory(void)
const JObject * Get_FinalParticle_SourceObject(const DKinematicData *locParticle)
static unsigned short int IsDetachedVertex(Particle_t p)
Definition: particleType.h:817
static char * ParticleType(Particle_t p)
Definition: particleType.h:142
DetectorSystem_t
Definition: GlueX.h:15
set< set< size_t > > Build_IndexCombos(const DReactionStep *locReactionStep, deque< Particle_t > locToIncludePIDs) const
double measuredBeta(void) const
double Calc_DOCAToVertex(const DVector3 &locUnitDir, const DVector3 &locPosition, const DVector3 &locVertex) const
const DKinematicData * Get_InitialParticle_Measured(void) const
TLorentzVector DLorentzVector
string Get_StepName(const DReactionStep *locStep, bool locIncludeMissingFlag, bool locTLatexFlag)
unsigned int Get_NDF_Timing(void) const
const DEventRFBunch * Get_EventRFBunch(void) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
static int ParticleCharge(Particle_t p)
double Get_Delta(DKinFitPullType locPullType, const DKinematicData *locMeasured, const DKinematicData *locKinFit)
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
double Calc_DOCA(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2) const
Definition: GlueX.h:17
bool Get_UseKinFitResultsFlag(void) const
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_EOverPVsTheta
Definition: GlueX.h:19
shared_ptr< const DTOFHitMatchParams > Get_TOFHitMatchParams(void) const
map< pair< int, Particle_t >, map< DKinFitPullType, TH1I * > > dHistMap_Pulls
JApplication * japp
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
vector< Particle_t > Get_ChainPIDs(const DReaction *locReaction, size_t locStepIndex, int locUpToStepIndex, vector< Particle_t > locUpThroughPIDs, bool locExpandDecayingFlag, bool locExcludeMissingFlag)
Definition: DReaction.cc:162
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
set< pair< size_t, pair< Particle_t, const JObject * > > > dPreviouslyHistogrammedParticles
void Initialize(JEventLoop *locEventLoop)
void Initialize(JEventLoop *locEventLoop)
void Get_HistNameTitle(DKinFitPullType locPullType, string locFullROOTName, int locVsKey, string &locHistName, string &locHistTitle)
double time(void) const
string Get_ConstraintInfo(const DReactionVertexInfo *locReactionVertexInfo, const DReaction *locReaction, size_t &locNumConstraints, size_t &locNumUnknowns) const
DLorentzVector Calc_FinalStateP4(const DReaction *locReaction, const DParticleCombo *locParticleCombo, size_t locStepIndex, bool locUseKinFitDataFlag) const
const DAnalysisUtilities * dAnalysisUtilities
deque< map< Particle_t, map< bool, TH1I * > > > dHistDeque_P
deque< map< Particle_t, TH1I * > > dHistDeque_TrackTToCommon
map< pair< int, Particle_t >, map< DKinFitPullType, TH2I * > > dHistMap_PullsVsTheta
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_EOverPVsP
Particle_t Get_InitialPID(void) const
Definition: DReactionStep.h:81
void Initialize(JEventLoop *locEventLoop)
deque< map< Particle_t, TH1I * > > dHistDeque_TrackDOCAToCommon
void Fill_Hists(JEventLoop *locEventLoop, const DKinematicData *locKinematicData, bool locIsMissingFlag, size_t locStepIndex)
Definition: GlueX.h:20
map< pair< int, Particle_t >, map< DKinFitPullType, TH2I * > > dHistMap_PullsVsP
Definition: GlueX.h:18
map< Particle_t, TH2I * > dHistMap_Ldiff_kpiVsP_DIRC
DLorentzVector Calc_MissingP4(const DReaction *locReaction, const DParticleCombo *locParticleCombo, bool locUseKinFitDataFlag) const
map< pair< int, Particle_t >, map< DKinFitPullType, TH1I * > > dHistMap_Deltas
double GetMostProbabledEdx_DC(double p, double mass, double dx, bool locIsCDCFlag) const
Definition: DParticleID.cc:657
set< set< pair< const JObject *, unsigned int > > > dPreviousSourceObjects
DGeometry * GetDGeometry(unsigned int run_number)
Definition: GlueX.h:22
map< Particle_t, TH2I * > dHistMap_ThetaCVsP_DIRC
deque< map< pair< Particle_t, Particle_t >, TH2I * > > dHistDeque_TrackDeltaTVsP
vector< const DKinematicData * > Get_FinalParticles(void) const
map< Particle_t, TH2I * > dHistMap_Ldiff_pkVsP_DIRC
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
size_t Get_NumFinalParticles(void) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
map< pair< int, Particle_t >, map< DKinFitPullType, TH2I * > > dHistMap_DeltasVsP
shared_ptr< const DDIRCMatchParams > Get_DIRCMatchParams(void) const
double Get_ConfidenceLevel(void) const
DLorentzVector Get_SpacetimeVertex(void) const
double GetdEdxSigma_DC(double num_hits, double p, double mass, double mean_path_length, bool locIsCDCFlag) const
Definition: DParticleID.cc:679
void Initialize(JEventLoop *locEventLoop)
deque< map< Particle_t, map< bool, TH1I * > > > dHistDeque_Theta
int Get_MissingParticleIndex(void) const
Definition: DReactionStep.h:93
shared_ptr< const DSCHitMatchParams > Get_SCHitMatchParams(void) const
void Get_DeltaBinningParams(DKinFitPullType locPullType, bool loc2DFlag, int &locNumBins, double &locMax)
set< pair< const DEventRFBunch *, pair< Particle_t, const JObject * > > > dPreviouslyHistogrammedParticles
map< Particle_t, TH1I * > dHistMap_PIDFOM
set< set< pair< const JObject *, unsigned int > > > dPreviousSourceObjects
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_dEdXPullVsP
DLorentzVector lorentzMomentum(void) const
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_dEdXVsP
static double ParticleMass(Particle_t p)
map< Particle_t, TH1I * > dHistMap_NumPhotons_DIRC
double Calc_TimingChiSq(const DChargedTrackHypothesis *locChargedHypo, unsigned int &locNDF, double &locTimingPull) const
shared_ptr< const DFCALShowerMatchParams > Get_FCALShowerMatchParams(void) const
const DKinematicData * Get_InitialParticle(void) const
deque< map< Particle_t, map< bool, TH1I * > > > dHistDeque_Phi
map< pair< int, Particle_t >, map< DKinFitPullType, TH2I * > > dHistMap_DeltasVsTheta
const DAnalysisUtilities * dAnalysisUtilities
unsigned int Get_NDF(void) const
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_DeltaTVsP
map< pair< int, Particle_t >, map< DKinFitPullType, TH2I * > > dHistMap_DeltasVsPhi
const DVector3 & momentum(void) const
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_DeltadEdXVsP
set< set< set< pair< const JObject *, unsigned int > > > > dPreviousSourceObjects
void Fill_NeutralHists(const DNeutralParticleHypothesis *locNeutralParticleHypothesis, const DMCThrownMatching *locMCThrownMatching, const DEventRFBunch *locEventRFBunch)
void Initialize(JEventLoop *locEventLoop)
const DAnalysisUtilities * dAnalysisUtilities
const DReactionStep * Get_ReactionStep(size_t locStepIndex) const
Definition: DReaction.h:84
map< Particle_t, map< DetectorSystem_t, TH1I * > > dHistMap_Beta
set< set< pair< const JObject *, unsigned int > > > dPreviousSourceObjects
const DKinFitResults * Get_KinFitResults(void) const
deque< map< Particle_t, map< bool, TH1I * > > > dHistDeque_VertexZ
set< set< set< pair< const JObject *, unsigned int > > > > dPreviousSourceObjects
DVector3 Get_Position(void) const
shared_ptr< const DBCALShowerMatchParams > Get_BCALShowerMatchParams(void) const
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_TimePullVsP
map< Particle_t, TH2I * > dHistMap_PVsTheta_NaNPIDFOM
const DParticleComboStep * Get_ParticleComboStep(size_t locStepIndex) const
string Get_InitialParticlesName(const DReactionStep *locStep, bool locTLatexFlag)
void GetScintMPdEandSigma(double p, double M, double x, double &most_prob_dE, double &sigma_dE) const
Definition: DParticleID.cc:622
deque< map< Particle_t, map< bool, TH2I * > > > dHistDeque_BetaVsP
const DKinematicData * Get_FinalParticle(size_t locFinalParticleIndex) const
size_t Get_NumReactionSteps(void) const
Definition: DReaction.h:83
map< pair< int, Particle_t >, TH2I * > dHistMap_ConfidenceLevel_VsTheta
int type
GEANT particle ID.
Definition: DMCThrown.h:20
deque< map< Particle_t, TH1I * > > dHistDeque_TrackZToCommon
DetectorSystem_t t1_detector(void) const
deque< map< Particle_t, map< bool, TH2I * > > > dHistDeque_PVsTheta
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_BetaVsP
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
map< pair< int, Particle_t >, TH2I * > dHistMap_ConfidenceLevel_VsP
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
Particle_t Get_MissingPID(void) const
Particle_t PID(void) const
const DMCThrown * Get_MatchingMCThrown(const DChargedTrackHypothesis *locChargedTrackHypothesis, double &locMatchFOM) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_DeltaBetaVsP
void Get_Pulls(map< const JObject *, map< DKinFitPullType, double > > &locPulls) const
double mass(void) const
map< pair< Particle_t, Particle_t >, TH1I * > dHistMap_PIDFOMForTruePID
map< Particle_t, map< DetectorSystem_t, TH2I * > > dHistMap_TimeFOMVsP
Particle_t
Definition: particleType.h:12
vector< const DReaction * > Get_Reactions(void) const
map< pair< int, Particle_t >, map< DKinFitPullType, TH2I * > > dHistMap_PullsVsPhi