Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DHistogramActions_Independent.cc
Go to the documentation of this file.
1 
2 #include <unistd.h>
3 
5 
6 void DHistogramAction_ObjectMemory::Initialize(JEventLoop* locEventLoop)
7 {
8  //setup binning
9  vector<string> locBinLabels = {"TMatrixFSym", "DKinematicInfo", "Charged DTimingInfo", "DTrackingInfo", "Neutral DTimingInfo", "KinematicDatas", "Charged Hypos", "Neutral Hypos", "Beam Photons",
10  "Combo RF Bunches", "Source Combos", "Source Combo Vectors", "Particle Combos", "Particle Combo Steps",
11  "DKinFitParticle", "DKinFitChainStep", "DKinFitChain", "DKinFitResults", "DKinFitConstraint_Mass", "DKinFitConstraint_P4", "DKinFitConstraint_Vertex", "DKinFitConstraint_Spacetime"};
12  for(size_t loc_i = 0; loc_i < locBinLabels.size(); ++loc_i)
13  dBinMap[locBinLabels[loc_i]] = loc_i + 1;
14 
15  //CREATE THE HISTOGRAMS
16  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
17  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
18  {
20 
21  dVirtualMemoryVsEventNumber = new TH1F("VirtualMemoryVsEventNumber", ";Event Counter;Virtual Memory (MB)", dMaxNumEvents, 0.5, (double)dMaxNumEvents + 0.5);
22  dResidentMemoryVsEventNumber = new TH1F("ResidentMemoryVsEventNumber", ";Event Counter;Resident Memory (MB)", dMaxNumEvents, 0.5, (double)dMaxNumEvents + 0.5);
23 
24  // Total Memory
25  string locHistName = "TotalMemory";
26  string locHistTitle = ";Event Counter;Total Memory (MB)";
27  dHist_TotalMemory = GetOrCreate_Histogram<TH1F>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5);
28 
29  // # Objects
30  locHistName = "NumObjects2D";
31  locHistTitle = "# Objects;Event Counter";
32  dHist_NumObjects = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5, locBinLabels.size(), 0.5, float(locBinLabels.size()) + 0.5);
33  for(size_t loc_i = 0; loc_i < locBinLabels.size(); ++loc_i)
34  dHist_NumObjects->GetYaxis()->SetBinLabel(1 + loc_i, locBinLabels[loc_i].c_str());
35 
36  // Object Memory
37  locHistName = "Memory2D";
38  locHistTitle = "Memory (MB);Event Counter";
39  dHist_Memory = GetOrCreate_Histogram<TH2F>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5, locBinLabels.size(), 0.5, float(locBinLabels.size()) + 0.5);
40  for(size_t loc_i = 0; loc_i < locBinLabels.size(); ++loc_i)
41  dHist_Memory->GetYaxis()->SetBinLabel(1 + loc_i, locBinLabels[loc_i].c_str());
42 
43  //Return to the base directory
45  }
46  japp->RootUnLock(); //RELEASE ROOT LOCK!!
47 }
48 
49 bool DHistogramAction_ObjectMemory::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
50 {
52  return true; //else double-counting!
53 
55  return true;
56 
57  if(locParticleCombo != nullptr)
58  return true; // Protect against infinite recursion (see below)
59 
60  //THIS IS EXTREMELY DANGEROUS, AND SHOULD BE AVOIDED UNLESS YOU KNOW !!!!EXACTLY!!!! WHAT YOU ARE DOING
61  //And if you're doing this, you probably don't
62  //This will result in a infinite-recursion crash if it is called by DAnalysisResults_factory.
63  //This action should only be used directly in an event processsor
64  //This casuses the analysis to run, generating the objects needed for histogramming below.
65  vector<const DAnalysisResults*> locAnalysisResults;
66  locEventLoop->Get(locAnalysisResults);
67 
68  map<int, size_t> locNumObjectsMap; //int is bin
69  map<int, double> locMemoryMap; //int is bin
70  double locTotalMemory = 0.0;
71 
72  /******************************************************* COMPONENT OBJECTS *******************************************************/
73 
74  //TMatrixFSym
75  auto locBin = dBinMap["TMatrixFSym"];
76  locNumObjectsMap[locBin] = dResourcePool_TMatrixFSym.Get_NumObjectsAllThreads();
77  auto locMemory = (sizeof(TMatrixDSym) + 7*7*4)*locNumObjectsMap[locBin]; //assume 7x7 matrix of floats (4)
78  locMemoryMap[locBin] = locMemory;
79  locTotalMemory += locMemory;
80 
81  //DKinematicData::DKinematicInfo
82  locBin = dBinMap["DKinematicInfo"];
83  locNumObjectsMap[locBin] = dResourcePool_KinematicInfo.Get_NumObjectsAllThreads();
84  locMemory = sizeof(DKinematicData::DKinematicInfo)*locNumObjectsMap[locBin];
85  locMemoryMap[locBin] = locMemory;
86  locTotalMemory += locMemory;
87 
88  //DChargedTrackHypothesis::DTimingInfo
89  locBin = dBinMap["Charged DTimingInfo"];
91  locMemory = sizeof(DChargedTrackHypothesis::DTimingInfo)*locNumObjectsMap[locBin];
92  locMemoryMap[locBin] = locMemory;
93  locTotalMemory += locMemory;
94 
95  //DChargedTrackHypothesis::DTrackingInfo
96  locBin = dBinMap["DTrackingInfo"];
98  locMemory = sizeof(DChargedTrackHypothesis::DTrackingInfo)*locNumObjectsMap[locBin];
99  locMemoryMap[locBin] = locMemory;
100  locTotalMemory += locMemory;
101 
102  //DNeutralParticleHypothesis::DTimingInfo
103  locBin = dBinMap["Neutral DTimingInfo"];
105  locMemory = sizeof(DNeutralParticleHypothesis::DTimingInfo)*locNumObjectsMap[locBin];
106  locMemoryMap[locBin] = locMemory;
107  locTotalMemory += locMemory;
108 
109  /******************************************************* KINEMATIC DATA OBJECTS *******************************************************/
110 
111  //DKinematicData
112  locBin = dBinMap["KinematicDatas"];
113  locNumObjectsMap[locBin] = dResourcePool_KinematicData.Get_NumObjectsAllThreads();
114  locMemory = sizeof(DKinematicData)*locNumObjectsMap[locBin];
115  locMemoryMap[locBin] = locMemory;
116  locTotalMemory += locMemory;
117 
118  //DChargedTrackHypothesis
119  locBin = dBinMap["Charged Hypos"];
121  locMemory = sizeof(DChargedTrackHypothesis)*locNumObjectsMap[locBin];
122  locMemoryMap[locBin] = locMemory;
123  locTotalMemory += locMemory;
124 
125  //DNeutralParticleHypothesis
126  locBin = dBinMap["Neutral Hypos"];
128  locMemory = sizeof(DNeutralParticleHypothesis)*locNumObjectsMap[locBin];
129  locMemoryMap[locBin] = locMemory;
130  locTotalMemory += locMemory;
131 
132  //DBeamPhoton
133  locBin = dBinMap["Beam Photons"];
134  locNumObjectsMap[locBin] = dResourcePool_BeamPhotons.Get_NumObjectsAllThreads();
135  locMemory = sizeof(DBeamPhoton)*locNumObjectsMap[locBin];
136  locMemoryMap[locBin] = locMemory;
137  locTotalMemory += locMemory;
138 
139  /******************************************************* COMBO OBJECTS *******************************************************/
140 
141  //DEventRFBunch
142  locBin = dBinMap["Combo RF Bunches"];
143  locNumObjectsMap[locBin] = dResourcePool_EventRFBunch.Get_NumObjectsAllThreads();
144  locMemory = sizeof(DEventRFBunch)*locNumObjectsMap[locBin];
145  locMemoryMap[locBin] = locMemory;
146  locTotalMemory += locMemory;
147 
148  //DSourceCombo
149  locBin = dBinMap["Source Combos"];
150  locNumObjectsMap[locBin] = dResourcePool_SourceCombo.Get_NumObjectsAllThreads();
151  locMemory = sizeof(DSourceCombo)*locNumObjectsMap[locBin];
152  locMemoryMap[locBin] = locMemory;
153  locTotalMemory += locMemory;
154 
155  //DSourceCombo Vectors
156  locBin = dBinMap["Source Combo Vectors"];
157  locNumObjectsMap[locBin] = dResourcePool_SourceComboVector.Get_NumObjectsAllThreads();
158  locMemory = sizeof(vector<const DSourceCombo*>)*locNumObjectsMap[locBin];
159  locMemoryMap[locBin] = locMemory;
160  locTotalMemory += locMemory;
161 
162  //DParticleCombo
163  locBin = dBinMap["Particle Combos"];
164  locNumObjectsMap[locBin] = dResourcePool_ParticleCombo.Get_NumObjectsAllThreads();
165  locMemory = sizeof(DParticleCombo)*locNumObjectsMap[locBin];
166  locMemoryMap[locBin] = locMemory;
167  locTotalMemory += locMemory;
168 
169  //DParticleComboStep
170  locBin = dBinMap["Particle Combo Steps"];
171  locNumObjectsMap[locBin] = dResourcePool_ParticleComboStep.Get_NumObjectsAllThreads();
172  locMemory = sizeof(DParticleComboStep)*locNumObjectsMap[locBin];
173  locMemoryMap[locBin] = locMemory;
174  locTotalMemory += locMemory;
175 
176  /******************************************************* KINFIT OBJECTS *******************************************************/
177 
178  //DKinFitParticle
179  locBin = dBinMap["DKinFitParticle"];
180  locNumObjectsMap[locBin] = dResourcePool_KinFitParticle.Get_NumObjectsAllThreads();
181  locMemory = sizeof(DKinFitParticle)*locNumObjectsMap[locBin];
182  locMemoryMap[locBin] = locMemory;
183  locTotalMemory += locMemory;
184 
185  //DKinFitChainStep
186  locBin = dBinMap["DKinFitChainStep"];
187  locNumObjectsMap[locBin] = dResourcePool_KinFitChainStep.Get_NumObjectsAllThreads();
188  locMemory = sizeof(DKinFitChainStep)*locNumObjectsMap[locBin];
189  locMemoryMap[locBin] = locMemory;
190  locTotalMemory += locMemory;
191 
192  //DKinFitChain
193  locBin = dBinMap["DKinFitChain"];
194  locNumObjectsMap[locBin] = dResourcePool_KinFitChain.Get_NumObjectsAllThreads();
195  locMemory = sizeof(DKinFitChain)*locNumObjectsMap[locBin];
196  locMemoryMap[locBin] = locMemory;
197  locTotalMemory += locMemory;
198 
199  //DKinFitResults
200  locBin = dBinMap["DKinFitResults"];
201  locNumObjectsMap[locBin] = dResourcePool_KinFitResults.Get_NumObjectsAllThreads();
202  locMemory = sizeof(DKinFitResults)*locNumObjectsMap[locBin];
203  locMemoryMap[locBin] = locMemory;
204  locTotalMemory += locMemory;
205 
206  /******************************************************* KINFIT CONSTRAINTS *******************************************************/
207 
208  //DKinFitConstraint_Mass
209  locBin = dBinMap["DKinFitConstraint_Mass"];
210  locNumObjectsMap[locBin] = dResourcePool_MassConstraint.Get_NumObjectsAllThreads();
211  locMemory = sizeof(DKinFitConstraint_Mass)*locNumObjectsMap[locBin];
212  locMemoryMap[locBin] = locMemory;
213  locTotalMemory += locMemory;
214 
215  //DKinFitConstraint_P4
216  locBin = dBinMap["DKinFitConstraint_P4"];
217  locNumObjectsMap[locBin] = dResourcePool_P4Constraint.Get_NumObjectsAllThreads();
218  locMemory = sizeof(DKinFitConstraint_P4)*locNumObjectsMap[locBin];
219  locMemoryMap[locBin] = locMemory;
220  locTotalMemory += locMemory;
221 
222  //DKinFitConstraint_Vertex
223  locBin = dBinMap["DKinFitConstraint_Vertex"];
224  locNumObjectsMap[locBin] = dResourcePool_VertexConstraint.Get_NumObjectsAllThreads();
225  locMemory = sizeof(DKinFitConstraint_Vertex)*locNumObjectsMap[locBin];
226  locMemoryMap[locBin] = locMemory;
227  locTotalMemory += locMemory;
228 
229  //DKinFitConstraint_Spacetime
230  locBin = dBinMap["DKinFitConstraint_Spacetime"];
231  locNumObjectsMap[locBin] = dResourcePool_SpacetimeConstraint.Get_NumObjectsAllThreads();
232  locMemory = sizeof(DKinFitConstraint_Spacetime)*locNumObjectsMap[locBin];
233  locMemoryMap[locBin] = locMemory;
234  locTotalMemory += locMemory;
235 
236  //Convert to MB
237  for(auto& locMemoryPair : locMemoryMap)
238  locMemoryPair.second /= (1024.0*1024.0);
239  locTotalMemory /= (1024.0*1024.0); //convert to MB
240 
241  //FILL HISTOGRAMS
242  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
243  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
244  Lock_Action(); //ACQUIRE ROOT LOCK!!
245  {
246  ++dEventCounter;
248  {
249  for(auto& locNumObjectsPair : locNumObjectsMap)
250  {
251  int locObjectBin = locNumObjectsPair.first;
252  dHist_NumObjects->SetBinContent(dEventCounter, locObjectBin, locNumObjectsMap[locObjectBin]);
253  dHist_Memory->SetBinContent(dEventCounter, locObjectBin, locMemoryMap[locObjectBin]);
254  }
255  dHist_TotalMemory->SetBinContent(dEventCounter, locTotalMemory);
256 
257  double vm, rss;
258  Read_MemoryUsage(vm, rss);
259  dVirtualMemoryVsEventNumber->SetBinContent(dEventCounter, vm / 1024.0);
260  dResidentMemoryVsEventNumber->SetBinContent(dEventCounter, rss / 1024.0);
261  }
262  }
263  Unlock_Action(); //RELEASE ROOT LOCK!!
264 
265  return true;
266 }
267 
268 void DHistogramAction_ObjectMemory::Read_MemoryUsage(double& vm_usage, double& resident_set)
269 {
270  vm_usage = 0.0;
271  resident_set = 0.0;
272 
273  // 'file' stat seems to give the most reliable results
274  ifstream stat_stream("/proc/self/stat",ios_base::in);
275 
276  // dummy vars for leading entries in stat that we don't care about
277  string pid, comm, state, ppid, pgrp, session, tty_nr;
278  string tpgid, flags, minflt, cminflt, majflt, cmajflt;
279  string utime, stime, cutime, cstime, priority, nice;
280  string O, itrealvalue, starttime;
281 
282  // the two fields we want
283  unsigned long vsize;
284  long rss;
285 
286  stat_stream >> pid >> comm >> state >> ppid >> pgrp >> session >> tty_nr
287  >> tpgid >> flags >> minflt >> cminflt >> majflt >> cmajflt
288  >> utime >> stime >> cutime >> cstime >> priority >> nice
289  >> O >> itrealvalue >> starttime >> vsize >> rss; // don't care about the rest
290 
291  stat_stream.close();
292 
293  long page_size_kb = sysconf(_SC_PAGE_SIZE) / 1024; // in case x86-64 is configured to use 2MB pages
294  vm_usage = vsize / 1024.0;
295  resident_set = rss * page_size_kb;
296 }
297 
298 void DHistogramAction_Reconstruction::Initialize(JEventLoop* locEventLoop)
299 {
300  //Create any histograms/trees/etc. within a ROOT lock.
301  //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
302 
303  //When creating a reaction-independent action, only modify member variables within a ROOT lock.
304  //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
305 
306  string locHistName, locHistTitle;
307 
308  //Check if is REST event (high-level objects only)
309  bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
310 
311  vector<const DMCThrown*> locMCThrowns;
312  locEventLoop->Get(locMCThrowns);
313 
314  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
315  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
316  double locTargetZCenter = 0.0;
317  locGeometry->GetTargetZ(locTargetZCenter);
318 
319  //CREATE THE HISTOGRAMS
320  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
321  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
322  {
323  if(dTargetCenter.Z() < -9.8E9)
324  dTargetCenter.SetXYZ(0.0, 0.0, locTargetZCenter);
325 
326  //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
327  //If another thread has already created the folder, it just changes to it.
329 
330  //FCAL
331  locHistName = "FCALShowerYVsX";
332  dHist_FCALShowerYVsX = GetOrCreate_Histogram<TH2I>(locHistName, ";FCAL Shower X (cm);FCAL Shower Y (cm)", dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
333  locHistName = "FCALShowerEnergy";
334  dHist_FCALShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
335 
336  //BCAL
337  locHistName = "BCALShowerEnergy";
338  dHist_BCALShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
339  locHistName = "BCALShowerPhi";
340  dHist_BCALShowerPhi = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #phi#circ", dNumPhiBins, dMinPhi, dMaxPhi);
341  locHistName = "BCALShowerPhiVsZ";
342  dHist_BCALShowerPhiVsZ = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Shower Z (cm);BCAL Shower #phi#circ", dNum2DBCALZBins, 0.0, 450.0, dNum2DPhiBins, dMinPhi, dMaxPhi);
343 
344  //TOF
345  locHistName = "TOFPointEnergy";
346  dHist_TOFPointEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";TOF Point Energy (MeV)", dNumHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
347  locHistName = "TOFPointYVsX";
348  dHist_TOFPointYVsX = GetOrCreate_Histogram<TH2I>(locHistName, ";TOF Point X (cm);TOF Point Y (cm)", dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
349 
350  //SC
351  locHistName = "SCHitSector";
352  dHist_SCHitSector = GetOrCreate_Histogram<TH1I>(locHistName, ";SC Hit Sector", 30, 0.5, 30.5);
353  locHistName = "SCHitEnergy";
354  dHist_SCHitEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";SC Hit Energy (MeV)", dNumHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
355  locHistName = "SCHitEnergyVsSector";
356  dHist_SCHitEnergyVsSector = GetOrCreate_Histogram<TH2I>(locHistName, ";SC Hit Sector;SC Hit Energy (MeV)", 30, 0.5, 30.5, dNum2DHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
357  locHistName = "SCRFDeltaTVsSector";
358  dHist_SCRFDeltaTVsSector = GetOrCreate_Histogram<TH2I>(locHistName, ";SC Hit Sector;SC - RF #Deltat (ns)", 30, 0.5, 30.5, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
359 
360  //TAGH, TAGM
361  locHistName = "TAGMRFDeltaTVsColumn";
362  dHist_TAGMRFDeltaTVsColumn = GetOrCreate_Histogram<TH2I>(locHistName, ";TAGM Column;TAGM - RF #Deltat (ns)", 102, 0.5, 102.5, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
363  locHistName = "TAGHRFDeltaTVsCounter";
364  dHist_TAGHRFDeltaTVsCounter = GetOrCreate_Histogram<TH2I>(locHistName, ";TAGH Counter;TAGH - RF #Deltat (ns)", 274, 0.5, 274.5, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
365 
366  //TRACKING
367  CreateAndChangeTo_Directory("Tracking", "Tracking");
368  locHistName = "NumDCHitsPerTrack";
369  dHist_NumDCHitsPerTrack = GetOrCreate_Histogram<TH1I>(locHistName, ";# Track Hits", 50, 0.5, 50.5);
370  locHistName = "NumPossDCHitsPerTrack";
371  dHist_NumPossDCHitsPerTrack = GetOrCreate_Histogram<TH1I>(locHistName, ";# Possible Track Hits", 50, 0.5, 50.5);
372  locHistName = "TrackHitFraction";
373  dHist_TrackHitFraction = GetOrCreate_Histogram<TH1I>(locHistName, ";# Track Hits/# Possible Track Hits", 50, 0.0, 1.0);
374  locHistName = "NumDCHitsPerTrackVsTheta";
375  dHist_NumDCHitsPerTrackVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, ";#theta#circ;# Track Hits", dNum2DThetaBins, dMinTheta, dMaxTheta, 46, 4.5, 50.5);
376  locHistName = "NumPossDCHitsPerTrackVsTheta";
377  dHist_NumPossDCHitsPerTrackVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, ";#theta#circ;# Possible Track Hits", dNum2DThetaBins, dMinTheta, dMaxTheta, 46, 4.5, 50.5);
378  locHistName = "TrackHitFractionVsTheta";
379  dHist_TrackHitFractionVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, ";#theta#circ;# Track Hits/# Possible Track Hits", dNum2DThetaBins, dMinTheta, dMaxTheta, 50, 0.0, 1.0);
380  locHistName = "TrackingFOM_WireBased";
381  dHist_TrackingFOM_WireBased = GetOrCreate_Histogram<TH1I>(locHistName, ";Confidence Level", dNumFOMBins, 0.0, 1.0);
382  locHistName = "TrackingFOM";
383  dHist_TrackingFOM = GetOrCreate_Histogram<TH1I>(locHistName, ";Confidence Level", dNumFOMBins, 0.0, 1.0);
384  locHistName = "TrackingFOMVsP";
385  dHist_TrackingFOMVsP = GetOrCreate_Histogram<TH2I>(locHistName, ";p (GeV/c);Confidence Level", dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
386  locHistName = "TrackingFOMVsTheta";
387  dHist_TrackingFOMVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, ";#theta#circ;Confidence Level", dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DFOMBins, 0.0, 1.0);
388  locHistName = "TrackingFOMVsNumHits";
389  dHist_TrackingFOMVsNumHits = GetOrCreate_Histogram<TH2I>(locHistName, ";# Track Hits;Confidence Level", 46, 4.5, 50.5, dNum2DFOMBins, 0.0, 1.0);
390 
391  if(!locIsRESTEvent)
392  {
393  locHistName = "CDCRingVsTheta_Candidates";
394  dHist_CDCRingVsTheta_Candidates = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Track Candidates;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
395  locHistName = "CDCRingVsTheta_WireBased";
396  dHist_CDCRingVsTheta_WireBased = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Wire-Based Tracks;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
397  }
398 
399  locHistName = "CDCRingVsTheta_TimeBased";
400  dHist_CDCRingVsTheta_TimeBased = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Time-Based Tracks;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
401  locHistName = "CDCRingVsTheta_TimeBased_GoodTrackFOM";
402  dHist_CDCRingVsTheta_TimeBased_GoodTrackFOM = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Good FOM Time-Based Tracks;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
403 
404  if(!locIsRESTEvent)
405  {
406  locHistName = "FDCPlaneVsTheta_Candidates";
407  dHist_FDCPlaneVsTheta_Candidates = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Track Candidates;p (GeV/c);FDC Plane", dNum2DThetaBins, dMinTheta, dMaxTheta, 24, 0.5, 24.5);
408  locHistName = "FDCPlaneVsTheta_WireBased";
409  dHist_FDCPlaneVsTheta_WireBased = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Wire-Based Tracks;p (GeV/c);FDC Plane", dNum2DThetaBins, dMinTheta, dMaxTheta, 24, 0.5, 24.5);
410  }
411 
412  locHistName = "FDCPlaneVsTheta_TimeBased";
413  dHist_FDCPlaneVsTheta_TimeBased = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Time-Based Tracks;p (GeV/c);FDC Plane", dNum2DThetaBins, dMinTheta, dMaxTheta, 24, 0.5, 24.5);
414  locHistName = "FDCPlaneVsTheta_TimeBased_GoodTrackFOM";
415  dHist_FDCPlaneVsTheta_TimeBased_GoodTrackFOM = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Good FOM Time-Based Tracks;p (GeV/c);FDC Plane", dNum2DThetaBins, dMinTheta, dMaxTheta, 24, 0.5, 24.5);
416 
417  if(!locMCThrowns.empty())
418  {
419  locHistName = "MCMatchedHitsVsTheta";
420  dHist_MCMatchedHitsVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, "Fraction of Track Hits Matched to MC;#theta#circ;", dNum2DThetaBins, dMinTheta, dMaxTheta, 100, 0.0, 1.0);
421  locHistName = "MCMatchedHitsVsP";
422  dHist_MCMatchedHitsVsP = GetOrCreate_Histogram<TH2I>(locHistName, "Fraction of Track Hits Matched to MC;p (GeV/c);", dNum2DPBins, dMinP, dMaxP, 100, 0.0, 1.0);
423  }
424 
425  for(int locCharge = -1; locCharge <= 1; locCharge += 2)
426  {
427  string locParticleROOTName = (locCharge == -1) ? "#it{q}^{-}" : "#it{q}^{+}";
428  string locParticleName = (locCharge == -1) ? "q-" : "q+";
429  CreateAndChangeTo_Directory(locParticleName, locParticleName);
430 
431  if(!locIsRESTEvent)
432  {
433  // PVsTheta Track Candidates
434  locHistName = string("PVsTheta_Candidates_") + locParticleName;
435  locHistTitle = locParticleROOTName + string(" Track Candidates;#theta#circ;p (GeV/c)");
436  dHistMap_PVsTheta_Candidates[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
437 
438  // PVsTheta Wire-Based Tracks
439  locHistName = string("PVsTheta_WireBased_") + locParticleName;
440  locHistTitle = locParticleROOTName + string(" Wire-Based Tracks;#theta#circ;p (GeV/c)");
441  dHistMap_PVsTheta_WireBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
442  }
443 
444  // PVsTheta Time-Based Tracks
445  locHistName = string("PVsTheta_TimeBased_") + locParticleName;
446  locHistTitle = locParticleROOTName + string(" Time-Based Tracks;#theta#circ;p (GeV/c)");
447  dHistMap_PVsTheta_TimeBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
448 
449  // PVsTheta Time-Based Tracks Good Track FOM
450  locHistName = string("PVsTheta_TimeBased_GoodTrackFOM_") + locParticleName;
451  locHistTitle = locParticleROOTName + string(" Time-Based Tracks, Good Tracking FOM;#theta#circ;p (GeV/c)");
452  dHistMap_PVsTheta_TimeBased_GoodTrackFOM[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
453 
454  // PVsTheta Time-Based Tracks Low Track FOM
455  locHistName = string("PVsTheta_TimeBased_LowTrackFOM_") + locParticleName;
456  locHistTitle = locParticleROOTName + string(" Time-Based Tracks, Low Tracking FOM;#theta#circ;p (GeV/c)");
457  dHistMap_PVsTheta_TimeBased_LowTrackFOM[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
458 
459  // PVsTheta Time-Based Tracks High Track FOM
460  locHistName = string("PVsTheta_TimeBased_HighTrackFOM_") + locParticleName;
461  locHistTitle = locParticleROOTName + string(" Time-Based Tracks, High Tracking FOM;#theta#circ;p (GeV/c)");
462  dHistMap_PVsTheta_TimeBased_HighTrackFOM[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
463 
464  if(!locIsRESTEvent)
465  {
466  // PVsTheta: Good Wire-Based, Good Time-Based
467  locHistName = string("PVsTheta_GoodWireBased_GoodTimeBased_") + locParticleName;
468  locHistTitle = locParticleROOTName + string(" Good Wire-Based, Good Time-Based;#theta#circ;p (GeV/c)");
469  dHistMap_PVsTheta_GoodWireBased_GoodTimeBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
470 
471  // PVsTheta: Good Wire-Based, Bad Time-Based
472  locHistName = string("PVsTheta_GoodWireBased_BadTimeBased_") + locParticleName;
473  locHistTitle = locParticleROOTName + string(" Good Wire-Based, Bad Time-Based;#theta#circ;p (GeV/c)");
474  dHistMap_PVsTheta_GoodWireBased_BadTimeBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
475  }
476 
477  gDirectory->cd(".."); //end of charge
478  }
479  gDirectory->cd(".."); //End of "Tracking"
480 
481  //Return to the base directory
483  }
484  japp->RootUnLock(); //RELEASE ROOT LOCK!!
485 }
486 
487 bool DHistogramAction_Reconstruction::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
488 {
489  //Expect locParticleCombo to be NULL since this is a reaction-independent action.
490 
492  return true; //else double-counting!
493 
494  bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
495 
496  vector<const DBCALShower*> locBCALShowers;
497  locEventLoop->Get(locBCALShowers);
498 
499  vector<const DFCALShower*> locFCALShowers;
500  locEventLoop->Get(locFCALShowers);
501 
502  vector<const DTOFPoint*> locTOFPoints;
503  locEventLoop->Get(locTOFPoints);
504 
505  vector<const DSCHit*> locSCHits;
506  locEventLoop->Get(locSCHits);
507 
508  vector<const DBeamPhoton*> locBeamPhotons;
509  locEventLoop->Get(locBeamPhotons);
510 
511  const DDetectorMatches* locDetectorMatches = NULL;
512  locEventLoop->GetSingle(locDetectorMatches);
513 
514  vector<const DTrackTimeBased*> locTrackTimeBasedVector;
515  locEventLoop->Get(locTrackTimeBasedVector);
516 
517  vector<const DMCThrownMatching*> locMCThrownMatchingVector;
518  locEventLoop->Get(locMCThrownMatchingVector);
519 
520  vector<const DMCThrown*> locMCThrowns;
521  locEventLoop->Get(locMCThrowns, "FinalState");
522 
523  const DParticleID* locParticleID = NULL;
524  locEventLoop->GetSingle(locParticleID);
525 
526  const DEventRFBunch* locEventRFBunch = NULL;
527  locEventLoop->GetSingle(locEventRFBunch);
528 
529  const DDetectorMatches* locDetectorMatches_WireBased = NULL;
530  vector<const DTrackCandidate*> locTrackCandidates;
531  vector<const DTrackWireBased*> locTrackWireBasedVector;
532  if(!locIsRESTEvent)
533  {
534  vector<const DDetectorMatches*> locDetectorMatchesVector_WireBased;
535  locEventLoop->Get(locDetectorMatchesVector_WireBased, "WireBased");
536  if(!locDetectorMatchesVector_WireBased.empty())
537  locDetectorMatches_WireBased = locDetectorMatchesVector_WireBased[0];
538  locEventLoop->Get(locTrackCandidates);
539  locEventLoop->Get(locTrackWireBasedVector);
540  }
541 
542  //select the best DTrackWireBased for each track: use best tracking FOM
543  map<JObject::oid_t, const DTrackWireBased*> locBestTrackWireBasedMap; //lowest tracking FOM for each candidate id
544  for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
545  {
546  JObject::oid_t locCandidateID = locTrackWireBasedVector[loc_i]->candidateid;
547  if(locBestTrackWireBasedMap.find(locCandidateID) == locBestTrackWireBasedMap.end())
548  locBestTrackWireBasedMap[locCandidateID] = locTrackWireBasedVector[loc_i];
549  else if(locTrackWireBasedVector[loc_i]->FOM > locBestTrackWireBasedMap[locCandidateID]->FOM)
550  locBestTrackWireBasedMap[locCandidateID] = locTrackWireBasedVector[loc_i];
551  }
552 
553  //select the best DTrackTimeBased for each track: use best tracking FOM
554  //also, make map from WBT -> TBT (if not rest)
555  //also, select best sc matches for each track
556  map<JObject::oid_t, const DTrackTimeBased*> locBestTrackTimeBasedMap; //lowest tracking FOM for each candidate id
557  map<const DTrackWireBased*, const DTrackTimeBased*> locWireToTimeBasedTrackMap;
558  map<const DTrackTimeBased*, shared_ptr<const DSCHitMatchParams>> locTimeBasedToBestSCMatchMap;
559  for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
560  {
561  //Best SC Match Params
562  shared_ptr<const DSCHitMatchParams> locSCHitMatchParams;
563  if(locParticleID->Get_BestSCMatchParams(locTrackTimeBasedVector[loc_i], locDetectorMatches, locSCHitMatchParams))
564  locTimeBasedToBestSCMatchMap[locTrackTimeBasedVector[loc_i]] = locSCHitMatchParams;
565 
566  JObject::oid_t locCandidateID = locTrackTimeBasedVector[loc_i]->candidateid;
567  if(locBestTrackTimeBasedMap.find(locCandidateID) == locBestTrackTimeBasedMap.end())
568  locBestTrackTimeBasedMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
569  else if(locTrackTimeBasedVector[loc_i]->FOM > locBestTrackTimeBasedMap[locCandidateID]->FOM)
570  locBestTrackTimeBasedMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
571  if(locIsRESTEvent)
572  continue;
573  const DTrackWireBased* locTrackWireBased = NULL;
574  locTrackTimeBasedVector[loc_i]->GetSingle(locTrackWireBased);
575  locWireToTimeBasedTrackMap[locTrackWireBased] = locTrackTimeBasedVector[loc_i];
576  }
577 
578  //FILL HISTOGRAMS
579  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
580  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
581  Lock_Action(); //ACQUIRE ROOT LOCK!!
582  {
583  //FCAL
584  for(size_t loc_i = 0; loc_i < locFCALShowers.size(); ++loc_i)
585  {
586  dHist_FCALShowerEnergy->Fill(locFCALShowers[loc_i]->getEnergy());
587  dHist_FCALShowerYVsX->Fill(locFCALShowers[loc_i]->getPosition().X(), locFCALShowers[loc_i]->getPosition().Y());
588  }
589 
590  //BCAL
591  for(size_t loc_i = 0; loc_i < locBCALShowers.size(); ++loc_i)
592  {
593  dHist_BCALShowerEnergy->Fill(locBCALShowers[loc_i]->E);
594 
595  DVector3 locBCALPosition(locBCALShowers[loc_i]->x, locBCALShowers[loc_i]->y, locBCALShowers[loc_i]->z);
596  double locBCALPhi = locBCALPosition.Phi()*180.0/TMath::Pi();
597  dHist_BCALShowerPhi->Fill(locBCALPhi);
598  dHist_BCALShowerPhiVsZ->Fill(locBCALPosition.Z(), locBCALPhi);
599  }
600 
601  //TOF
602  for(size_t loc_i = 0; loc_i < locTOFPoints.size(); ++loc_i)
603  {
604  dHist_TOFPointEnergy->Fill(locTOFPoints[loc_i]->dE*1.0E3);
605  dHist_TOFPointYVsX->Fill(locTOFPoints[loc_i]->pos.X(), locTOFPoints[loc_i]->pos.Y());
606  }
607 
608  //SC
609  for(size_t loc_i = 0; loc_i < locSCHits.size(); ++loc_i)
610  {
611  dHist_SCHitSector->Fill(locSCHits[loc_i]->sector);
612  dHist_SCHitEnergy->Fill(locSCHits[loc_i]->dE*1.0E3);
613  dHist_SCHitEnergyVsSector->Fill(locSCHits[loc_i]->sector, locSCHits[loc_i]->dE*1.0E3);
614  }
615 
616  //TAGM, TAGH
617  for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
618  {
619  double locDeltaT = locBeamPhotons[loc_i]->time() - locEventRFBunch->dTime;
620  if(locBeamPhotons[loc_i]->dSystem == SYS_TAGM)
621  dHist_TAGMRFDeltaTVsColumn->Fill(locBeamPhotons[loc_i]->dCounter, locDeltaT);
622  else
623  dHist_TAGHRFDeltaTVsCounter->Fill(locBeamPhotons[loc_i]->dCounter, locDeltaT);
624  }
625 
626  //TRACK CANDIDATES
627  for(size_t loc_i = 0; loc_i < locTrackCandidates.size(); ++loc_i)
628  {
629  int locCharge = (locTrackCandidates[loc_i]->charge() > 0.0) ? 1 : -1;
630  double locTheta = locTrackCandidates[loc_i]->momentum().Theta()*180.0/TMath::Pi();
631  double locP = locTrackCandidates[loc_i]->momentum().Mag();
632  dHistMap_PVsTheta_Candidates[locCharge]->Fill(locTheta, locP);
633 
634  set<int> locCDCRings;
635  locParticleID->Get_CDCRings(locTrackCandidates[loc_i]->dCDCRings, locCDCRings);
636  for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator)
637  dHist_CDCRingVsTheta_Candidates->Fill(locTheta, *locIterator);
638 
639  set<int> locFDCPlanes;
640  locParticleID->Get_FDCPlanes(locTrackCandidates[loc_i]->dFDCPlanes, locFDCPlanes);
641  for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator)
642  dHist_FDCPlaneVsTheta_Candidates->Fill(locTheta, *locIterator);
643  }
644 
645  //WIRE-BASED TRACKS
646  map<JObject::oid_t, const DTrackWireBased*>::iterator locWireBasedIterator = locBestTrackWireBasedMap.begin();
647  for(; locWireBasedIterator != locBestTrackWireBasedMap.end(); ++locWireBasedIterator)
648  {
649  const DTrackWireBased* locTrackWireBased = locWireBasedIterator->second;
650  int locCharge = (locTrackWireBased->charge() > 0.0) ? 1 : -1;
651  double locTheta = locTrackWireBased->momentum().Theta()*180.0/TMath::Pi();
652  double locP = locTrackWireBased->momentum().Mag();
653  dHistMap_PVsTheta_WireBased[locCharge]->Fill(locTheta, locP);
654 
655  set<int> locCDCRings;
656  locParticleID->Get_CDCRings(locTrackWireBased->dCDCRings, locCDCRings);
657  for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator)
658  dHist_CDCRingVsTheta_WireBased->Fill(locTheta, *locIterator);
659 
660  set<int> locFDCPlanes;
661  locParticleID->Get_FDCPlanes(locTrackWireBased->dFDCPlanes, locFDCPlanes);
662  for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator)
663  dHist_FDCPlaneVsTheta_WireBased->Fill(locTheta, *locIterator);
664 
665  dHist_TrackingFOM_WireBased->Fill(locTrackWireBased->FOM);
666  }
667 
668  //TIME-BASED TRACKS
669  map<JObject::oid_t, const DTrackTimeBased*>::iterator locTimeBasedIterator = locBestTrackTimeBasedMap.begin();
670  for(; locTimeBasedIterator != locBestTrackTimeBasedMap.end(); ++locTimeBasedIterator)
671  {
672  const DTrackTimeBased* locTrackTimeBased = locTimeBasedIterator->second;
673  int locCharge = (locTrackTimeBased->charge() > 0.0) ? 1 : -1;
674  double locTheta = locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi();
675  double locP = locTrackTimeBased->momentum().Mag();
676  int locPossibleHits = locTrackTimeBased->potential_cdc_hits_on_track + locTrackTimeBased->potential_fdc_hits_on_track;
677  double locTrackHitFraction = (locTrackTimeBased->Ndof + 5)/(double)locPossibleHits;
678 
679  dHistMap_PVsTheta_TimeBased[locCharge]->Fill(locTheta, locP);
680  dHist_NumDCHitsPerTrack->Fill(locTrackTimeBased->Ndof + 5);
681  dHist_NumPossDCHitsPerTrack->Fill(locPossibleHits);
682  dHist_TrackHitFraction->Fill(locTrackHitFraction);
683  dHist_NumDCHitsPerTrackVsTheta->Fill(locTheta, locTrackTimeBased->Ndof + 5);
684  dHist_NumPossDCHitsPerTrackVsTheta->Fill(locTheta, locPossibleHits);
685  dHist_TrackHitFractionVsTheta->Fill(locTheta, locTrackHitFraction);
686 
687  dHist_TrackingFOM->Fill(locTrackTimeBased->FOM);
688  dHist_TrackingFOMVsTheta->Fill(locTheta, locTrackTimeBased->FOM);
689  dHist_TrackingFOMVsP->Fill(locP, locTrackTimeBased->FOM);
690  dHist_TrackingFOMVsNumHits->Fill(locTrackTimeBased->Ndof + 5, locTrackTimeBased->FOM);
691 
692  //CDC
693  set<int> locCDCRings;
694  locParticleID->Get_CDCRings(locTrackTimeBased->dCDCRings, locCDCRings);
695  for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator)
696  {
697  dHist_CDCRingVsTheta_TimeBased->Fill(locTheta, *locIterator);
698  if(locTrackTimeBased->FOM > dGoodTrackFOM)
699  dHist_CDCRingVsTheta_TimeBased_GoodTrackFOM->Fill(locTheta, *locIterator);
700  }
701 
702  //FDC
703  set<int> locFDCPlanes;
704  locParticleID->Get_FDCPlanes(locTrackTimeBased->dFDCPlanes, locFDCPlanes);
705  for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator)
706  {
707  dHist_FDCPlaneVsTheta_TimeBased->Fill(locTheta, *locIterator);
708  if(locTrackTimeBased->FOM > dGoodTrackFOM)
709  dHist_FDCPlaneVsTheta_TimeBased_GoodTrackFOM->Fill(locTheta, *locIterator);
710  }
711 
712  //FOM
713  if(locTrackTimeBased->FOM > dGoodTrackFOM)
714  dHistMap_PVsTheta_TimeBased_GoodTrackFOM[locCharge]->Fill(locTheta, locP);
715  else
716  dHistMap_PVsTheta_TimeBased_LowTrackFOM[locCharge]->Fill(locTheta, locP);
717  if(locTrackTimeBased->FOM > dHighTrackFOM)
718  dHistMap_PVsTheta_TimeBased_HighTrackFOM[locCharge]->Fill(locTheta, locP);
719 
720  //SC/RF DELTA-T
721  auto locSCIterator = locTimeBasedToBestSCMatchMap.find(locTrackTimeBased);
722  if(locSCIterator != locTimeBasedToBestSCMatchMap.end())
723  {
724  auto& locSCHitMatchParams = locSCIterator->second;
725  double locPropagatedSCTime = locSCHitMatchParams->dHitTime - locSCHitMatchParams->dFlightTime + (dTargetCenter.Z() - locTrackTimeBased->z())/29.9792458;
726  double locDeltaT = locPropagatedSCTime - locEventRFBunch->dTime;
727  dHist_SCRFDeltaTVsSector->Fill(locSCHitMatchParams->dSCHit->sector, locDeltaT);
728  }
729  }
730 
731  // If "Good" WBT, see if TBT is good
732  locWireBasedIterator = locBestTrackWireBasedMap.begin();
733  for(; locWireBasedIterator != locBestTrackWireBasedMap.end(); ++locWireBasedIterator)
734  {
735  if(locDetectorMatches_WireBased == NULL)
736  continue;
737  const DTrackWireBased* locTrackWireBased = locWireBasedIterator->second;
738  if(locTrackWireBased->FOM < dGoodTrackFOM)
739  continue; //no good
740  if(!locDetectorMatches_WireBased->Get_IsMatchedToHit(locTrackWireBased))
741  continue; //no good
742 
743  int locCharge = (locTrackWireBased->charge() > 0.0) ? 1 : -1;
744  double locTheta = locTrackWireBased->momentum().Theta()*180.0/TMath::Pi();
745  double locP = locTrackWireBased->momentum().Mag();
746 
747  map<const DTrackWireBased*, const DTrackTimeBased*>::iterator locReconIterator = locWireToTimeBasedTrackMap.find(locTrackWireBased);
748  if(locReconIterator == locWireToTimeBasedTrackMap.end())
749  {
750  dHistMap_PVsTheta_GoodWireBased_BadTimeBased[locCharge]->Fill(locTheta, locP);
751  continue; //no time-based
752  }
753 
754  const DTrackTimeBased* locTrackTimeBased = locReconIterator->second;
755  if((locTrackTimeBased->FOM < dGoodTrackFOM) || (!locDetectorMatches->Get_IsMatchedToHit(locTrackTimeBased)))
756  dHistMap_PVsTheta_GoodWireBased_BadTimeBased[locCharge]->Fill(locTheta, locP);
757  else
758  dHistMap_PVsTheta_GoodWireBased_GoodTimeBased[locCharge]->Fill(locTheta, locP);
759  }
760 
761  //THROWN
762  for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
763  {
764  if(fabs(locMCThrowns[loc_i]->charge()) < 0.9)
765  continue;
766 
767  double locMatchFOM;
768  const DChargedTrackHypothesis* locChargedTrackHypothesis = locMCThrownMatchingVector[0]->Get_MatchingChargedHypothesis(locMCThrowns[loc_i], locMatchFOM);
769  if(locChargedTrackHypothesis == NULL)
770  continue;
771 
772  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
773  double locHitFraction = 1.0*locTrackTimeBased->dNumHitsMatchedToThrown/(locTrackTimeBased->Ndof + 5);
774  dHist_MCMatchedHitsVsTheta->Fill(locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi(), locHitFraction);
775  dHist_MCMatchedHitsVsP->Fill(locTrackTimeBased->momentum().Mag(), locHitFraction);
776  }
777  }
778  Unlock_Action(); //RELEASE ROOT LOCK!!
779 
780  return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
781 }
782 
783 void DHistogramAction_DetectorMatching::Initialize(JEventLoop* locEventLoop)
784 {
785  //Create any histograms/trees/etc. within a ROOT lock.
786  //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
787 
788  //When creating a reaction-independent action, only modify member variables within a ROOT lock.
789  //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
790  string locHistName, locHistTitle;
791 
792  bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
793 
794  map<string, double> tofparms;
795  locEventLoop->GetCalib("TOF/tof_parms", tofparms);
796 
797  //CREATE THE HISTOGRAMS
798  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
799  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
800  {
801  TOF_E_THRESHOLD = tofparms["TOF_E_THRESHOLD"];
802 
803  //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
804  //If another thread has already created the folder, it just changes to it.
806 
807  //Loop over filling for time-based and wire-based tracks
808  for(unsigned int locDummy = 0; locDummy < 2; ++locDummy)
809  {
810  bool locIsTimeBased = (locDummy == 0);
811  if(locIsRESTEvent && (!locIsTimeBased))
812  continue;
813 
814  string locDirectoryName = locIsTimeBased ? "TimeBased" : "WireBased";
815  CreateAndChangeTo_Directory(locDirectoryName.c_str(), locDirectoryName.c_str());
816  string locTrackString = locIsTimeBased ? "Time-Based Tracks" : "Wire-Based Tracks";
817 
818  //Kinematics of has (no) hit
819  vector<DetectorSystem_t> locDetectorSystems;
820  locDetectorSystems.push_back(SYS_START); locDetectorSystems.push_back(SYS_BCAL);
821  locDetectorSystems.push_back(SYS_TOF); locDetectorSystems.push_back(SYS_FCAL);
822  for(size_t loc_i = 0; loc_i < locDetectorSystems.size(); ++loc_i)
823  {
824  DetectorSystem_t locSystem = locDetectorSystems[loc_i];
825 
826  double locMaxTheta = ((locSystem == SYS_FCAL) || (locSystem == SYS_TOF)) ? 12.0 : dMaxTheta;
827  double locMaxP = (locSystem == SYS_BCAL) ? 3.0 : dMaxP;
828 
829  string locSystemName = SystemName(locSystem);
830  if(locSystemName == "ST")
831  locSystemName = "SC";
832  string locDirName = locSystemName;
833  if(locSystemName == "TOF")
834  locDirName = "TOFPoint";
835 
836  CreateAndChangeTo_Directory(locDirName, locDirName);
837 
838  // PVsTheta Has Hit
839  locHistName = "PVsTheta_HasHit";
840  locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" Has Hit;#theta#circ;p (GeV/c)");
841  dHistMap_PVsTheta_HasHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPBins, dMinP, locMaxP);
842 
843  // PVsTheta Has No Hit
844  locHistName = "PVsTheta_NoHit";
845  locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" No Hit;#theta#circ;p (GeV/c)");
846  dHistMap_PVsTheta_NoHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPBins, dMinP, locMaxP);
847 
848  // PhiVsTheta Has Hit
849  locHistName = "PhiVsTheta_HasHit";
850  locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" Has Hit;#theta#circ;#phi#circ");
851  dHistMap_PhiVsTheta_HasHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
852 
853  // PhiVsTheta Has No Hit
854  locHistName = "PhiVsTheta_NoHit";
855  locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" No Hit;#theta#circ;#phi#circ");
856  dHistMap_PhiVsTheta_NoHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
857 
858  gDirectory->cd("..");
859  }
860 
861  //SC
862  CreateAndChangeTo_Directory("SC", "SC");
863  locHistName = "SCPaddleVsTheta_HasHit";
864  locHistTitle = locTrackString + string(", Has Other Match, SC Has Hit;#theta#circ;Projected SC Paddle");
865  dHistMap_SCPaddleVsTheta_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, 30, 0.5, 30.5);
866 
867  locHistName = "SCPaddleVsTheta_NoHit";
868  locHistTitle = locTrackString + string(", Has Other Match, SC No Hit;#theta#circ;Projected SC Paddle");
869  dHistMap_SCPaddleVsTheta_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, 30, 0.5, 30.5);
870 
871  locHistName = "SCPaddleVsZ_HasHit";
872  locHistTitle = locTrackString + string(", Has Other Match, SC Has Hit;Projected SC Hit-Z (cm);Projected SC Paddle");
873  dHistMap_SCPaddleVsZ_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DSCZBins, 0.0, 120.0, 30, 0.5, 30.5);
874 
875  locHistName = "SCPaddleVsZ_NoHit";
876  locHistTitle = locTrackString + string(", Has Other Match, SC No Hit;Projected SC Hit-Z (cm);Projected SC Paddle");
877  dHistMap_SCPaddleVsZ_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DSCZBins, 0.0, 120.0, 30, 0.5, 30.5);
878 
879  locHistName = "SCPaddle_BarrelRegion_HasHit";
880  locHistTitle = locTrackString + string(", Has Other Match, SC Barrel Region Has Hit;Projected SC Paddle");
881  dHistMap_SCPaddle_BarrelRegion_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
882 
883  locHistName = "SCPaddle_BarrelRegion_NoHit";
884  locHistTitle = locTrackString + string(", Has Other Match, SC Barrel Region No Hit;Projected SC Paddle");
885  dHistMap_SCPaddle_BarrelRegion_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
886 
887  locHistName = "SCPaddle_NoseRegion_HasHit";
888  locHistTitle = locTrackString + string(", Has Other Match, SC Front Region Has Hit;Projected SC Paddle");
889  dHistMap_SCPaddle_NoseRegion_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
890 
891  locHistName = "SCPaddle_NoseRegion_NoHit";
892  locHistTitle = locTrackString + string(", Has Other Match, SC Front Region No Hit;Projected SC Paddle");
893  dHistMap_SCPaddle_NoseRegion_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
894 
895  locHistName = "SCTrackDeltaPhiVsP";
896  locHistTitle = locTrackString + string(";p (GeV/c);SC / Track #Delta#phi#circ");
897  dHistMap_SCTrackDeltaPhiVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaPhiBins, dSCMatchMinDeltaPhi, dSCMatchMaxDeltaPhi);
898 
899  locHistName = "SCTrackDeltaPhiVsTheta";
900  locHistTitle = locTrackString + string(";#theta#circ;SC / Track #Delta#phi#circ");
901  dHistMap_SCTrackDeltaPhiVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaPhiBins, dSCMatchMinDeltaPhi, dSCMatchMaxDeltaPhi);
902 
903  locHistName = "SCTrackDeltaPhiVsZ";
904  locHistTitle = locTrackString + string(";Projected SC Hit-Z (cm);SC / Track #Delta#phi#circ");
905  dHistMap_SCTrackDeltaPhiVsZ[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DSCZBins, 0.0, 120.0, dNum2DDeltaPhiBins, dSCMatchMinDeltaPhi, dSCMatchMaxDeltaPhi);
906  gDirectory->cd("..");
907 
908  //TOFPaddle
909  CreateAndChangeTo_Directory("TOFPaddle", "TOFPaddle");
910  locHistName = "VerticalPaddleTrackDeltaX";
911  locHistTitle = locTrackString + string(";Vertical TOF Paddle / Track |#DeltaX| (cm)");
912  dHistMap_TOFPaddleTrackDeltaX[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins/2, 0.0, dMaxTrackMatchDOCA);
913 
914  locHistName = "HorizontalPaddleTrackDeltaY";
915  locHistTitle = locTrackString + string(";Horizontal TOF Paddle / Track |#DeltaY| (cm)");
916  dHistMap_TOFPaddleTrackDeltaY[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins/2, 0.0, dMaxTrackMatchDOCA);
917 
918  locHistName = "TrackYVsVerticalPaddle_HasHit";
919  locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle Has Hit;Projected Vertical Paddle;Projected TOF Hit Y (cm)");
920  dHistMap_TOFPaddleTrackYVsVerticalPaddle_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNumFCALTOFXYBins, -130.0, 130.0);
921 
922  locHistName = "TrackYVsVerticalPaddle_NoHit";
923  locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle No Hit;Projected Vertical Paddle;Projected TOF Hit Y (cm)");
924  dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNumFCALTOFXYBins, -130.0, 130.0);
925 
926  locHistName = "HorizontalPaddleVsTrackX_HasHit";
927  locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle Has Hit;Projected TOF Hit X (cm);Projected Horizontal Paddle");
928  dHistMap_TOFPaddleHorizontalPaddleVsTrackX_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, 44, 0.5, 44.5);
929 
930  locHistName = "HorizontalPaddleVsTrackX_NoHit";
931  locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle No Hit;Projected TOF Hit X (cm);Projected Horizontal Paddle");
932  dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, 44, 0.5, 44.5);
933  gDirectory->cd("..");
934 
935  //TOFPoint
936  CreateAndChangeTo_Directory("TOFPoint", "TOFPoint");
937  locHistName = "TrackTOFYVsX_HasHit";
938  locHistTitle = locTrackString + string(", Has Other Match, TOF Has Hit;Projected TOF Hit X (cm);Projected TOF Hit Y (cm)");
939  dHistMap_TrackTOFYVsX_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
940 
941  locHistName = "TrackTOFYVsX_NoHit";
942  locHistTitle = locTrackString + string(", Has Other Match, TOF No Hit;Projected TOF Hit X (cm);Projected TOF Hit Y (cm)");
943  dHistMap_TrackTOFYVsX_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
944 
945  locHistName = "TrackTOF2DPaddles_HasHit";
946  locHistTitle = locTrackString + string(", Has Other Match, TOF Has Hit;Projected Vertical TOF Paddle;Projected Horizontal TOF Paddle");
947  dHistMap_TrackTOF2DPaddles_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, 44, 0.5, 44.5);
948 
949  locHistName = "TrackTOF2DPaddles_NoHit";
950  locHistTitle = locTrackString + string(", Has Other Match, TOF No Hit;Projected Vertical TOF Paddle;Projected Horizontal TOF Paddle");
951  dHistMap_TrackTOF2DPaddles_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, 44, 0.5, 44.5);
952 
953  locHistName = "TrackTOFP_HasHit";
954  locHistTitle = locTrackString + string(", Has Other Match, TOF Has Hit;p (GeV/c)");
955  dHistMap_TrackTOFP_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
956 
957  locHistName = "TrackTOFP_NoHit";
958  locHistTitle = locTrackString + string(", Has Other Match, TOF No Hit;p (GeV/c)");
959  dHistMap_TrackTOFP_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
960 
961  locHistName = "TrackTOFR_HasHit";
962  locHistTitle = locTrackString + string(", Has Other Match, TOF Has Hit;Projected TOF Hit R (cm)");
963  dHistMap_TrackTOFR_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTOFRBins, 0.0, 180);
964 
965  locHistName = "TrackTOFR_NoHit";
966  locHistTitle = locTrackString + string(", Has Other Match, TOF No Hit;Projected TOF Hit R (cm)");
967  dHistMap_TrackTOFR_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTOFRBins, 0.0, 180);
968 
969  locHistName = "TOFTrackDistanceVsP";
970  locHistTitle = locTrackString + string(";p (GeV/c);TOF / Track Distance (cm)");
971  dHistMap_TOFPointTrackDistanceVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
972 
973  locHistName = "TOFTrackDistanceVsTheta";
974  locHistTitle = locTrackString + string(";#theta#circ;TOF / Track Distance (cm)");
975  dHistMap_TOFPointTrackDistanceVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, 20.0, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
976 
977  locHistName = "TOFTrackDeltaXVsHorizontalPaddle";
978  locHistTitle = locTrackString + string(";TOF Horizontal Paddle;TOF / Track #DeltaX (cm)");
979  dHistMap_TOFPointTrackDeltaXVsHorizontalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
980 
981  locHistName = "TOFTrackDeltaXVsVerticalPaddle";
982  locHistTitle = locTrackString + string(";TOF Vertical Paddle;TOF / Track #DeltaX (cm)");
983  dHistMap_TOFPointTrackDeltaXVsVerticalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
984 
985  locHistName = "TOFTrackDeltaYVsHorizontalPaddle";
986  locHistTitle = locTrackString + string(";TOF Horizontal Paddle;TOF / Track #DeltaY (cm)");
987  dHistMap_TOFPointTrackDeltaYVsHorizontalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
988 
989  locHistName = "TOFTrackDeltaYVsVerticalPaddle";
990  locHistTitle = locTrackString + string(";TOF Vertical Paddle;TOF / Track #DeltaY (cm)");
991  dHistMap_TOFPointTrackDeltaYVsVerticalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
992 
993  locHistName = "TOFTrackDistance_BothPlanes";
994  locHistTitle = locTrackString + string("TOF Hit in Both Planes;TOF / Track Distance (cm)");
995  dHistMap_TOFPointTrackDistance_BothPlanes[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
996 
997  locHistName = "TOFTrackDistance_OnePlane";
998  locHistTitle = locTrackString + string("TOF Hit in One Plane;TOF / Track Distance (cm)");
999  dHistMap_TOFPointTrackDistance_OnePlane[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
1000  gDirectory->cd("..");
1001 
1002  //FCAL
1003  CreateAndChangeTo_Directory("FCAL", "FCAL");
1004  locHistName = "TrackFCALYVsX_HasHit";
1005  locHistTitle = locTrackString + string(", Has Other Match, FCAL Has Hit;Projected FCAL Hit X (cm);Projected FCAL Hit Y (cm)");
1006  dHistMap_TrackFCALYVsX_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
1007 
1008  locHistName = "TrackFCALYVsX_NoHit";
1009  locHistTitle = locTrackString + string(", Has Other Match, FCAL No Hit;Projected FCAL Hit X (cm);Projected FCAL Hit Y (cm)");
1010  dHistMap_TrackFCALYVsX_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
1011 
1012  locHistName = "TrackFCALP_HasHit";
1013  locHistTitle = locTrackString + string(", Has Other Match, FCAL Has Hit;p (GeV/c)");
1014  dHistMap_TrackFCALP_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
1015 
1016  locHistName = "TrackFCALP_NoHit";
1017  locHistTitle = locTrackString + string(", Has Other Match, FCAL No Hit;p (GeV/c)");
1018  dHistMap_TrackFCALP_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
1019 
1020  locHistName = "TrackFCALR_HasHit";
1021  locHistTitle = locTrackString + string(", Has Other Match, FCAL Has Hit;Projected FCAL Hit R (cm)");
1022  dHistMap_TrackFCALR_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFCALTOFXYBins, 0.0, 130);
1023 
1024  locHistName = "TrackFCALR_NoHit";
1025  locHistTitle = locTrackString + string(", Has Other Match, FCAL No Hit;Projected FCAL Hit R (cm)");
1026  dHistMap_TrackFCALR_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFCALTOFXYBins, 0.0, 130);
1027 
1028  locHistName = "TrackFCALRowVsColumn_HasHit";
1029  locHistTitle = locTrackString + string(", Has Other Match, FCAL Has Hit;Projected FCAL Hit Column;Projected FCAL Hit Row");
1030  dHistMap_TrackFCALRowVsColumn_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 59, -0.5, 58.5, 59, -0.5, 58.5);
1031 
1032  locHistName = "TrackFCALRowVsColumn_NoHit";
1033  locHistTitle = locTrackString + string(", Has Other Match, FCAL No Hit;Projected FCAL Hit Column;Projected FCAL Hit Row");
1034  dHistMap_TrackFCALRowVsColumn_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 59, -0.5, 58.5, 59, -0.5, 58.5);
1035 
1036  locHistName = "FCALTrackDistanceVsP";
1037  locHistTitle = locTrackString + string(";p (GeV/c);FCAL / Track Distance (cm)");
1038  dHistMap_FCALTrackDistanceVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
1039 
1040  locHistName = "FCALTrackDistanceVsTheta";
1041  locHistTitle = locTrackString + string(";#theta#circ;FCAL / Track Distance (cm)");
1042  dHistMap_FCALTrackDistanceVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, 20.0, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
1043  gDirectory->cd("..");
1044 
1045  //BCAL
1046  CreateAndChangeTo_Directory("BCAL", "BCAL");
1047  locHistName = "TrackBCALModuleVsZ_HasHit";
1048  locHistTitle = locTrackString + string(", Has Other Match, BCAL Has Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit Module");
1049  dHistMap_TrackBCALModuleVsZ_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, 48, 0.5, 48.5);
1050 
1051  locHistName = "TrackBCALModuleVsZ_NoHit";
1052  locHistTitle = locTrackString + string(", Has Other Match, BCAL No Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit Module");
1053  dHistMap_TrackBCALModuleVsZ_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, 48, 0.5, 48.5);
1054 
1055  locHistName = "TrackBCALPhiVsZ_HasHit";
1056  locHistTitle = locTrackString + string(", Has Other Match, BCAL Has Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit #phi#circ");
1057  dHistMap_TrackBCALPhiVsZ_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, dNum2DPhiBins, dMinPhi, dMaxPhi);
1058 
1059  locHistName = "TrackBCALPhiVsZ_NoHit";
1060  locHistTitle = locTrackString + string(", Has Other Match, BCAL No Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit #phi#circ");
1061  dHistMap_TrackBCALPhiVsZ_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, dNum2DPhiBins, dMinPhi, dMaxPhi);
1062 
1063  locHistName = "BCALDeltaPhiVsP";
1064  locHistTitle = locTrackString + string(";p (GeV/c);BCAL / Track #Delta#phi#circ");
1065  dHistMap_BCALDeltaPhiVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, 4.0, dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
1066 
1067  locHistName = "BCALDeltaPhiVsZ";
1068  locHistTitle = locTrackString + string(";Projected BCAL Hit-Z (cm);BCAL / Track #Delta#phi#circ");
1069  dHistMap_BCALDeltaPhiVsZ[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
1070 
1071  locHistName = "BCALDeltaZVsTheta";
1072  locHistTitle = locTrackString + string(";#theta#circ;BCAL / Track #Deltaz (cm)");
1073  dHistMap_BCALDeltaZVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaZBins, dMinDeltaZ, dMaxDeltaZ);
1074 
1075  locHistName = "BCALDeltaZVsZ";
1076  locHistTitle = locTrackString + string(";Projected BCAL Hit-Z (cm);BCAL / Track #Deltaz (cm)");
1077  dHistMap_BCALDeltaZVsZ[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, dNum2DDeltaZBins, dMinDeltaZ, dMaxDeltaZ);
1078  gDirectory->cd("..");
1079 
1080  //TRACKING
1081  locHistName = "PVsTheta_NoHitMatch";
1082  locHistTitle = locTrackString + string(", No Hit Match;#theta#circ;p (GeV/c)");
1083  dHistMap_TrackPVsTheta_NoHitMatch[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
1084 
1085  locHistName = "PVsTheta_HitMatch";
1086  locHistTitle = locTrackString + string(", Hit Match;#theta#circ;p (GeV/c)");
1087  dHistMap_TrackPVsTheta_HitMatch[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
1088 
1089  gDirectory->cd("..");
1090  }
1091 
1092  //Return to the base directory
1094  }
1095  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1096 }
1097 
1098 bool DHistogramAction_DetectorMatching::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1099 {
1100  //Expect locParticleCombo to be NULL since this is a reaction-independent action.
1101 
1103  return true; //else double-counting!
1104 
1105  bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
1106 
1107  Fill_MatchingHists(locEventLoop, true); //Time-based tracks
1108  if(!locIsRESTEvent)
1109  Fill_MatchingHists(locEventLoop, false); //Wire-based tracks
1110 
1111  return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
1112 }
1113 
1114 void DHistogramAction_DetectorMatching::Fill_MatchingHists(JEventLoop* locEventLoop, bool locIsTimeBased)
1115 {
1116  const DParticleID* locParticleID = NULL;
1117  locEventLoop->GetSingle(locParticleID);
1118 
1119  //can't make this a class member: may cause race condition
1121  locCutAction_TrackHitPattern.Initialize(locEventLoop);
1122 
1123  //get the best tracks for each candidate id, based on good hit pattern & tracking FOM
1124  map<JObject::oid_t, const DTrackingData*> locBestTrackMap; //lowest tracking FOM for each candidate id
1125  if(locIsTimeBased)
1126  {
1127  vector<const DTrackTimeBased*> locTrackTimeBasedVector;
1128  locEventLoop->Get(locTrackTimeBasedVector);
1129 
1130  //select the best DTrackTimeBased for each track: of tracks with good hit pattern, use best tracking FOM
1131  for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
1132  {
1133  if(locTrackTimeBasedVector[loc_i]->FOM < dMinTrackingFOM)
1134  continue;
1135  if(!locCutAction_TrackHitPattern.Cut_TrackHitPattern(locParticleID, locTrackTimeBasedVector[loc_i]))
1136  continue;
1137  JObject::oid_t locCandidateID = locTrackTimeBasedVector[loc_i]->candidateid;
1138  if(locBestTrackMap.find(locCandidateID) == locBestTrackMap.end())
1139  locBestTrackMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
1140  else if(locTrackTimeBasedVector[loc_i]->FOM > (dynamic_cast<const DTrackTimeBased*>(locBestTrackMap[locCandidateID]))->FOM)
1141  locBestTrackMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
1142  }
1143  }
1144  else
1145  {
1146  vector<const DTrackWireBased*> locTrackWireBasedVector;
1147  locEventLoop->Get(locTrackWireBasedVector);
1148 
1149  //select the best DTrackWireBased for each track: of tracks with good hit pattern, use best tracking FOM
1150  for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
1151  {
1152  if(locTrackWireBasedVector[loc_i]->FOM < dMinTrackingFOM)
1153  continue;
1154  if(!locCutAction_TrackHitPattern.Cut_TrackHitPattern(locParticleID, locTrackWireBasedVector[loc_i]))
1155  continue;
1156  JObject::oid_t locCandidateID = locTrackWireBasedVector[loc_i]->candidateid;
1157  if(locBestTrackMap.find(locCandidateID) == locBestTrackMap.end())
1158  locBestTrackMap[locCandidateID] = locTrackWireBasedVector[loc_i];
1159  else if(locTrackWireBasedVector[loc_i]->FOM > (dynamic_cast<const DTrackWireBased*>(locBestTrackMap[locCandidateID]))->FOM)
1160  locBestTrackMap[locCandidateID] = locTrackWireBasedVector[loc_i];
1161  }
1162  }
1163 
1164  vector<const DBCALShower*> locBCALShowers;
1165  locEventLoop->Get(locBCALShowers);
1166 
1167  vector<const DFCALShower*> locFCALShowers;
1168  locEventLoop->Get(locFCALShowers);
1169 
1170  vector<const DTOFPoint*> locTOFPoints;
1171  locEventLoop->Get(locTOFPoints);
1172 
1173  vector<const DTOFPaddleHit*> locTOFPaddleHits;
1174  locEventLoop->Get(locTOFPaddleHits);
1175 
1176  vector<const DSCHit*> locSCHits;
1177  locEventLoop->Get(locSCHits);
1178 
1179  const DEventRFBunch* locEventRFBunch = nullptr;
1180  locEventLoop->GetSingle(locEventRFBunch);
1181 
1182  string locDetectorMatchesTag = locIsTimeBased ? "" : "WireBased";
1183  const DDetectorMatches* locDetectorMatches = NULL;
1184  locEventLoop->GetSingle(locDetectorMatches, locDetectorMatchesTag.c_str());
1185 
1186  //TRACK / BCAL CLOSEST MATCHES
1187  map<const DTrackingData*, pair<shared_ptr<const DBCALShowerMatchParams>, double> > locBCALTrackDistanceMap; //double = z
1188  for(auto locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1189  {
1190  map<DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >extrapolations;
1191  Get_Extrapolations(locTrackIterator->second,extrapolations);
1192 
1193  if(extrapolations.size()==0)
1194  break; //e.g. REST data: no trajectory
1195 
1196  double locStartTime = locTrackIterator->second->t0();
1197  double locStartTimeVariance = 0.0;
1198  DVector3 locProjPos, locProjMom;
1199 
1200  shared_ptr<const DBCALShowerMatchParams> locBestMatchParams;
1201  if(locParticleID->Get_ClosestToTrack(extrapolations.at(SYS_BCAL), locBCALShowers, false, locStartTime, locBestMatchParams, &locStartTimeVariance, &locProjPos, &locProjMom))
1202  locBCALTrackDistanceMap[locTrackIterator->second] = std::make_pair(locBestMatchParams, locProjPos.Z());
1203  }
1204 
1205  //TRACK / FCAL CLOSEST MATCHES
1206  map<const DTrackingData*, shared_ptr<const DFCALShowerMatchParams>> locFCALTrackDistanceMap;
1207  for(auto locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1208  {
1209  map<DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >extrapolations;
1210  Get_Extrapolations(locTrackIterator->second,extrapolations);
1211 
1212  if(extrapolations.size()==0)
1213  break; //e.g. REST data: no trajectory
1214 
1215  double locStartTime = locTrackIterator->second->t0();
1216  shared_ptr<const DFCALShowerMatchParams> locBestMatchParams;
1217  if(locParticleID->Get_ClosestToTrack(extrapolations.at(SYS_FCAL), locFCALShowers, false, locStartTime, locBestMatchParams))
1218  locFCALTrackDistanceMap.emplace(locTrackIterator->second, locBestMatchParams);
1219  }
1220 
1221  //TRACK / SC CLOSEST MATCHES
1222  map<const DTrackingData*, pair<shared_ptr<const DSCHitMatchParams>, double> > locSCTrackDistanceMap; //double = z
1223  for(auto locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1224  {
1225  map<DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >extrapolations;
1226  Get_Extrapolations(locTrackIterator->second,extrapolations);
1227 
1228  if(extrapolations.size()==0)
1229  break; //e.g. REST data: no trajectory
1230 
1231  double locStartTime = locTrackIterator->second->t0();
1232  double locStartTimeVariance = 0.0;
1233  DVector3 locProjPos, locProjMom;
1234 
1235  shared_ptr<const DSCHitMatchParams> locBestMatchParams;
1236  if(locParticleID->Get_ClosestToTrack(extrapolations.at(SYS_START), locSCHits, locIsTimeBased, false, locStartTime, locBestMatchParams, &locStartTimeVariance, &locProjPos, &locProjMom))
1237  locSCTrackDistanceMap[locTrackIterator->second] = std::make_pair(locBestMatchParams, locProjPos.Z());
1238  }
1239 
1240  //TRACK / TOF POINT CLOSEST MATCHES
1241  map<const DTrackingData*, shared_ptr<const DTOFHitMatchParams>> locTOFPointTrackDistanceMap;
1242  for(auto locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1243  {
1244  map<DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >extrapolations;
1245  Get_Extrapolations(locTrackIterator->second,extrapolations);
1246 
1247  if(extrapolations.size()==0)
1248  break; //e.g. REST data: no trajectory
1249 
1250  double locStartTime = locTrackIterator->second->t0();
1251  shared_ptr<const DTOFHitMatchParams> locBestMatchParams;
1252  if(locParticleID->Get_ClosestToTrack(extrapolations.at(SYS_TOF), locTOFPoints, false, locStartTime, locBestMatchParams))
1253  locTOFPointTrackDistanceMap.emplace(locTrackIterator->second, locBestMatchParams);
1254  }
1255 
1256  //TRACK / TOF PADDLE CLOSEST MATCHES
1257  map<const DTrackingData*, pair<const DTOFPaddleHit*, pair<double, double> > > locHorizontalTOFPaddleTrackDistanceMap; //doubles: delta-y, distance
1258  map<const DTrackingData*, pair<const DTOFPaddleHit*, pair<double, double> > > locVerticalTOFPaddleTrackDistanceMap; //doubles: delta-x, distance
1259  for(auto locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1260  {
1261  const DKinematicData* locKinematicData = locTrackIterator->second;
1262  map<DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >extrapolations;
1263  Get_Extrapolations(locTrackIterator->second,extrapolations);
1264 
1265  if(extrapolations.size()==0)
1266  break; //e.g. REST data: no trajectory
1267 
1268  double locBestDeltaX = 999.9, locBestDeltaY = 999.9, locBestDistance_Vertical = 999.9, locBestDistance_Horizontal = 999.9;
1269  double locStartTime = locParticleID->Calc_PropagatedRFTime(locKinematicData, locEventRFBunch);
1270 
1271  const DTOFPaddleHit* locClosestTOFPaddleHit_Vertical = locParticleID->Get_ClosestTOFPaddleHit_Vertical(extrapolations.at(SYS_TOF), locTOFPaddleHits, locStartTime, locBestDeltaX, locBestDistance_Vertical);
1272  pair<double, double> locDistancePair_Vertical(locBestDeltaX, locBestDistance_Vertical);
1273  if(locClosestTOFPaddleHit_Vertical != NULL)
1274  locVerticalTOFPaddleTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPaddleHit*, pair<double, double> >(locClosestTOFPaddleHit_Vertical, locDistancePair_Vertical);
1275 
1276  const DTOFPaddleHit* locClosestTOFPaddleHit_Horizontal = locParticleID->Get_ClosestTOFPaddleHit_Horizontal(extrapolations.at(SYS_TOF), locTOFPaddleHits, locStartTime, locBestDeltaY, locBestDistance_Horizontal);
1277  pair<double, double> locDistancePair_Horizontal(locBestDeltaY, locBestDistance_Horizontal);
1278  if(locClosestTOFPaddleHit_Horizontal != NULL)
1279  locHorizontalTOFPaddleTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPaddleHit*, pair<double, double> >(locClosestTOFPaddleHit_Horizontal, locDistancePair_Horizontal);
1280  }
1281 
1282  //PROJECTED HIT POSITIONS
1283  map<const DTrackingData*, pair<int, bool> > locProjectedSCPaddleMap; //pair: paddle, hit-barrel-flag (false if bend/nose)
1284  map<const DTrackingData*, pair<int, int> > locProjectedTOF2DPaddlesMap; //pair: vertical, horizontal
1285  map<const DTrackingData*, pair<float, float> > locProjectedTOFXYMap; //pair: x, y
1286  map<const DTrackingData*, pair<int, int> > locProjectedFCALRowColumnMap; //pair: column, row
1287  map<const DTrackingData*, pair<float, float> > locProjectedFCALXYMap; //pair: x, y
1288  map<const DTrackingData*, pair<float, int> > locProjectedBCALModuleSectorMap; //pair: z, module
1289  map<const DTrackingData*, pair<float, float> > locProjectedBCALPhiZMap; //pair: z, phi
1290  for(auto locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1291  {
1292  const DTrackingData* locTrack = locTrackIterator->second;
1293  map<DetectorSystem_t,vector<DTrackFitter::Extrapolation_t> >extrapolations;
1294  Get_Extrapolations(locTrackIterator->second,extrapolations);
1295 
1296  if(extrapolations.size()==0)
1297  break; //e.g. REST data: no trajectory
1298 
1299  //SC
1300  DVector3 locSCIntersection;
1301  bool locProjBarrelFlag = false;
1302  unsigned int locProjectedSCPaddle = locParticleID->PredictSCSector(extrapolations.at(SYS_START), &locSCIntersection, &locProjBarrelFlag);
1303  if(locProjectedSCPaddle != 0)
1304  locProjectedSCPaddleMap[locTrack] = pair<int, bool>(locProjectedSCPaddle, locProjBarrelFlag);
1305 
1306  //TOF
1307  DVector3 locTOFIntersection;
1308  unsigned int locHorizontalBar = 0, locVerticalBar = 0;
1309  if(locParticleID->PredictTOFPaddles(extrapolations.at(SYS_TOF), locHorizontalBar, locVerticalBar, &locTOFIntersection))
1310  {
1311  locProjectedTOF2DPaddlesMap[locTrack] = pair<int, int>(locVerticalBar, locHorizontalBar);
1312  locProjectedTOFXYMap[locTrack] = pair<float, float>(locTOFIntersection.X(), locTOFIntersection.Y());
1313  }
1314 
1315  //FCAL
1316  DVector3 locFCALIntersection;
1317  unsigned int locRow = 0, locColumn = 0;
1318  if(locParticleID->PredictFCALHit(extrapolations.at(SYS_FCAL), locRow, locColumn, &locFCALIntersection))
1319  {
1320  locProjectedFCALRowColumnMap[locTrack] = pair<int, int>(locColumn, locRow);
1321  locProjectedFCALXYMap[locTrack] = pair<float, float>(locFCALIntersection.X(), locFCALIntersection.Y());
1322  }
1323 
1324  //BCAL
1325  DVector3 locBCALIntersection;
1326  unsigned int locModule = 0, locSector = 0;
1327  if(locParticleID->PredictBCALWedge(extrapolations.at(SYS_BCAL), locModule, locSector, &locBCALIntersection))
1328  {
1329  locProjectedBCALModuleSectorMap[locTrack] = pair<float, int>(locBCALIntersection.Z(), locModule);
1330  locProjectedBCALPhiZMap[locTrack] = pair<float, float>(locBCALIntersection.Z(), locBCALIntersection.Phi()*180.0/TMath::Pi());
1331  }
1332  }
1333 
1334  //FILL HISTOGRAMS
1335  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1336  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1337  Lock_Action(); //ACQUIRE ROOT LOCK!!
1338  {
1339  /********************************************************** MATCHING DISTANCE **********************************************************/
1340 
1341  //BCAL
1342  for(auto locMapPair : locBCALTrackDistanceMap)
1343  {
1344  auto locTrack = locMapPair.first;
1345  auto locMatchParams = locMapPair.second.first;
1346  double locDeltaPhi = locMatchParams->dDeltaPhiToShower*180.0/TMath::Pi();
1347  double locProjectedZ = locMapPair.second.second;
1348  dHistMap_BCALDeltaPhiVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDeltaPhi);
1349  dHistMap_BCALDeltaZVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locMatchParams->dDeltaZToShower);
1350 
1351  dHistMap_BCALDeltaPhiVsZ[locIsTimeBased]->Fill(locProjectedZ, locDeltaPhi);
1352  dHistMap_BCALDeltaZVsZ[locIsTimeBased]->Fill(locProjectedZ, locMatchParams->dDeltaZToShower);
1353  }
1354 
1355  //FCAL
1356  for(auto locMapPair : locFCALTrackDistanceMap)
1357  {
1358  auto locTrack = locMapPair.first;
1359  dHistMap_FCALTrackDistanceVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locMapPair.second->dDOCAToShower);
1360  dHistMap_FCALTrackDistanceVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locMapPair.second->dDOCAToShower);
1361  }
1362 
1363  //TOF Paddle
1364  //Horizontal
1365  auto locTOFPaddleIterator = locHorizontalTOFPaddleTrackDistanceMap.begin();
1366  for(; locTOFPaddleIterator != locHorizontalTOFPaddleTrackDistanceMap.end(); ++locTOFPaddleIterator)
1367  {
1368  double locDeltaY = locTOFPaddleIterator->second.second.first;
1369  dHistMap_TOFPaddleTrackDeltaY[locIsTimeBased]->Fill(locDeltaY);
1370  }
1371  //Vertical
1372  locTOFPaddleIterator = locVerticalTOFPaddleTrackDistanceMap.begin();
1373  for(; locTOFPaddleIterator != locVerticalTOFPaddleTrackDistanceMap.end(); ++locTOFPaddleIterator)
1374  {
1375  double locDeltaX = locTOFPaddleIterator->second.second.first;
1376  dHistMap_TOFPaddleTrackDeltaX[locIsTimeBased]->Fill(locDeltaX);
1377  }
1378 
1379  //TOF Point
1380  for(auto locMapPair : locTOFPointTrackDistanceMap)
1381  {
1382  auto locTrack = locMapPair.first;
1383  const DTOFPoint* locTOFPoint = locMapPair.second->dTOFPoint;
1384  double locDeltaX = locMapPair.second->dDeltaXToHit;
1385  double locDeltaY = locMapPair.second->dDeltaYToHit;
1386 
1387  double locDistance = sqrt(locDeltaX*locDeltaX + locDeltaY*locDeltaY);
1388  if((fabs(locDeltaX) < 500.0) && (fabs(locDeltaY) < 500.0)) //else position not well-defined
1389  {
1390  dHistMap_TOFPointTrackDistanceVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDistance);
1391  dHistMap_TOFPointTrackDistanceVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locDistance);
1392  if((locTOFPoint->dHorizontalBar != 0) && (locTOFPoint->dVerticalBar != 0))
1393  dHistMap_TOFPointTrackDistance_BothPlanes[locIsTimeBased]->Fill(locDistance);
1394  else
1395  dHistMap_TOFPointTrackDistance_OnePlane[locIsTimeBased]->Fill(locDistance);
1396  }
1397 
1398  dHistMap_TOFPointTrackDeltaXVsHorizontalPaddle[locIsTimeBased]->Fill(locTOFPoint->dHorizontalBar, locDeltaX);
1399  dHistMap_TOFPointTrackDeltaXVsVerticalPaddle[locIsTimeBased]->Fill(locTOFPoint->dVerticalBar, locDeltaX);
1400 
1401  dHistMap_TOFPointTrackDeltaYVsHorizontalPaddle[locIsTimeBased]->Fill(locTOFPoint->dHorizontalBar, locDeltaY);
1402  dHistMap_TOFPointTrackDeltaYVsVerticalPaddle[locIsTimeBased]->Fill(locTOFPoint->dVerticalBar, locDeltaY);
1403  }
1404 
1405  //SC
1406  if(locSCHits.size() <= 4) //don't fill if every paddle fired!
1407  {
1408  for(auto locMapPair : locSCTrackDistanceMap)
1409  {
1410  auto locTrack = locMapPair.first;
1411  auto locMatchParams = locMapPair.second.first;
1412  double locDeltaPhi = locMatchParams->dDeltaPhiToHit*180.0/TMath::Pi();
1413  double locProjectedZ = locMapPair.second.second;
1414  dHistMap_SCTrackDeltaPhiVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDeltaPhi);
1415  dHistMap_SCTrackDeltaPhiVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locDeltaPhi);
1416  dHistMap_SCTrackDeltaPhiVsZ[locIsTimeBased]->Fill(locProjectedZ, locDeltaPhi);
1417  }
1418  }
1419 
1420  /********************************************************* MATCHING EFFICINECY *********************************************************/
1421 
1422  //Does-it-match, by detector
1423  for(auto locMapPair : locBestTrackMap)
1424  {
1425  auto locTrack = locMapPair.second;
1426  double locTheta = locTrack->momentum().Theta()*180.0/TMath::Pi();
1427  double locPhi = locTrack->momentum().Phi()*180.0/TMath::Pi();
1428  double locP = locTrack->momentum().Mag();
1429 
1430  //BCAL
1431  if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_START))
1432  {
1433  if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_BCAL))
1434  {
1435  dHistMap_PVsTheta_HasHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locP);
1436  dHistMap_PhiVsTheta_HasHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1437  if(locProjectedBCALModuleSectorMap.find(locTrack) != locProjectedBCALModuleSectorMap.end())
1438  {
1439  pair<float, float>& locPositionPair = locProjectedBCALPhiZMap[locTrack];
1440  dHistMap_TrackBCALPhiVsZ_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1441  pair<float, int>& locElementPair = locProjectedBCALModuleSectorMap[locTrack];
1442  dHistMap_TrackBCALModuleVsZ_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1443  }
1444  }
1445  else
1446  {
1447  dHistMap_PVsTheta_NoHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locP);
1448  dHistMap_PhiVsTheta_NoHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1449  if(locProjectedBCALModuleSectorMap.find(locTrack) != locProjectedBCALModuleSectorMap.end())
1450  {
1451  pair<float, float>& locPositionPair = locProjectedBCALPhiZMap[locTrack];
1452  dHistMap_TrackBCALPhiVsZ_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1453  pair<float, int>& locElementPair = locProjectedBCALModuleSectorMap[locTrack];
1454  dHistMap_TrackBCALModuleVsZ_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1455  }
1456  }
1457  }
1458 
1459  //FCAL
1460  if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1461  {
1462  dHistMap_PVsTheta_HasHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locP);
1463  if(locP > 1.0)
1464  dHistMap_PhiVsTheta_HasHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1465  if(locProjectedFCALRowColumnMap.find(locTrack) != locProjectedFCALRowColumnMap.end())
1466  {
1467  dHistMap_TrackFCALP_HasHit[locIsTimeBased]->Fill(locP);
1468  if(locP > 1.0)
1469  {
1470  pair<float, float>& locPositionPair = locProjectedFCALXYMap[locTrack];
1471  dHistMap_TrackFCALYVsX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1472  float locFCALR = sqrt(locPositionPair.first*locPositionPair.first + locPositionPair.second*locPositionPair.second);
1473  dHistMap_TrackFCALR_HasHit[locIsTimeBased]->Fill(locFCALR);
1474  pair<int, int>& locElementPair = locProjectedFCALRowColumnMap[locTrack];
1475  dHistMap_TrackFCALRowVsColumn_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1476  }
1477  }
1478  }
1479  else
1480  {
1481  dHistMap_PVsTheta_NoHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locP);
1482  if(locP > 1.0)
1483  dHistMap_PhiVsTheta_NoHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1484  if(locProjectedFCALRowColumnMap.find(locTrack) != locProjectedFCALRowColumnMap.end())
1485  {
1486  dHistMap_TrackFCALP_NoHit[locIsTimeBased]->Fill(locP);
1487  if(locP > 1.0)
1488  {
1489  pair<float, float>& locPositionPair = locProjectedFCALXYMap[locTrack];
1490  dHistMap_TrackFCALYVsX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1491  float locFCALR = sqrt(locPositionPair.first*locPositionPair.first + locPositionPair.second*locPositionPair.second);
1492  dHistMap_TrackFCALR_NoHit[locIsTimeBased]->Fill(locFCALR);
1493  pair<int, int>& locElementPair = locProjectedFCALRowColumnMap[locTrack];
1494  dHistMap_TrackFCALRowVsColumn_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1495  }
1496  }
1497  }
1498 
1499  //TOF Paddle
1500  if((locP > 1.0) && (locProjectedTOFXYMap.find(locTrack) != locProjectedTOFXYMap.end()))
1501  {
1502  pair<float, float>& locPositionPair = locProjectedTOFXYMap[locTrack];
1503  pair<int, int>& locPaddlePair = locProjectedTOF2DPaddlesMap[locTrack]; //vertical, horizontal
1504 
1505  //Horizontal
1506  if(locHorizontalTOFPaddleTrackDistanceMap.find(locTrack) != locHorizontalTOFPaddleTrackDistanceMap.end())
1507  {
1508  auto& locMatch = locHorizontalTOFPaddleTrackDistanceMap[locTrack];
1509  const DTOFPaddleHit* locTOFPaddleHit = locMatch.first;
1510  bool locDoubleEndedHitFlag = ((locTOFPaddleHit->E_north > TOF_E_THRESHOLD) && (locTOFPaddleHit->E_south > TOF_E_THRESHOLD));
1511  double locDistance = locDoubleEndedHitFlag ? locMatch.second.second : locMatch.second.first;
1512  if(locDistance <= dMinTOFPaddleMatchDistance) //match
1513  dHistMap_TOFPaddleHorizontalPaddleVsTrackX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1514  else //no match
1515  dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1516  }
1517  else // no match
1518  dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1519 
1520  //Vertical
1521  if(locVerticalTOFPaddleTrackDistanceMap.find(locTrack) != locVerticalTOFPaddleTrackDistanceMap.end())
1522  {
1523  auto& locMatch = locVerticalTOFPaddleTrackDistanceMap[locTrack];
1524  const DTOFPaddleHit* locTOFPaddleHit = locMatch.first;
1525  bool locDoubleEndedHitFlag = ((locTOFPaddleHit->E_north > TOF_E_THRESHOLD) && (locTOFPaddleHit->E_south > TOF_E_THRESHOLD));
1526  double locDistance = locDoubleEndedHitFlag ? locMatch.second.second : locMatch.second.first;
1527  if(locDistance <= dMinTOFPaddleMatchDistance) //match
1528  dHistMap_TOFPaddleTrackYVsVerticalPaddle_HasHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1529  else //no match
1530  dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1531  }
1532  else // no match
1533  dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1534  }
1535 
1536  //TOF Point
1537  if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1538  {
1539  dHistMap_PVsTheta_HasHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locP);
1540  if(locP > 1.0)
1541  dHistMap_PhiVsTheta_HasHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locPhi);
1542  if(locProjectedTOFXYMap.find(locTrack) != locProjectedTOFXYMap.end())
1543  {
1544  dHistMap_TrackTOFP_HasHit[locIsTimeBased]->Fill(locP);
1545  if(locP > 1.0)
1546  {
1547  pair<float, float>& locPositionPair = locProjectedTOFXYMap[locTrack];
1548  dHistMap_TrackTOFYVsX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1549  float locTOFR = sqrt(locPositionPair.first*locPositionPair.first + locPositionPair.second*locPositionPair.second);
1550  dHistMap_TrackTOFR_HasHit[locIsTimeBased]->Fill(locTOFR);
1551  pair<int, int>& locElementPair = locProjectedTOF2DPaddlesMap[locTrack];
1552  dHistMap_TrackTOF2DPaddles_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1553  }
1554  }
1555  }
1556  else
1557  {
1558  dHistMap_PVsTheta_NoHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locP);
1559  if(locP > 1.0)
1560  dHistMap_PhiVsTheta_NoHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locPhi);
1561  if(locProjectedTOFXYMap.find(locTrack) != locProjectedTOFXYMap.end())
1562  {
1563  dHistMap_TrackTOFP_NoHit[locIsTimeBased]->Fill(locP);
1564  if(locP > 1.0)
1565  {
1566  pair<float, float>& locPositionPair = locProjectedTOFXYMap[locTrack];
1567  dHistMap_TrackTOFYVsX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1568  float locTOFR = sqrt(locPositionPair.first*locPositionPair.first + locPositionPair.second*locPositionPair.second);
1569  dHistMap_TrackTOFR_NoHit[locIsTimeBased]->Fill(locTOFR);
1570  pair<int, int>& locElementPair = locProjectedTOF2DPaddlesMap[locTrack];
1571  dHistMap_TrackTOF2DPaddles_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1572  }
1573  }
1574  }
1575 
1576  //SC
1577  if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_BCAL) || locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL) || locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1578  {
1579  if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_START))
1580  {
1581  dHistMap_PVsTheta_HasHit[SYS_START][locIsTimeBased]->Fill(locTheta, locP);
1582  dHistMap_PhiVsTheta_HasHit[SYS_START][locIsTimeBased]->Fill(locTheta, locPhi);
1583  if(locProjectedSCPaddleMap.find(locTrack) != locProjectedSCPaddleMap.end())
1584  {
1585  dHistMap_SCPaddleVsTheta_HasHit[locIsTimeBased]->Fill(locTheta, locProjectedSCPaddleMap[locTrack].first);
1586  if(locProjectedSCPaddleMap[locTrack].second)
1587  dHistMap_SCPaddle_BarrelRegion_HasHit[locIsTimeBased]->Fill(locProjectedSCPaddleMap[locTrack].first);
1588  else
1589  dHistMap_SCPaddle_NoseRegion_HasHit[locIsTimeBased]->Fill(locProjectedSCPaddleMap[locTrack].first);
1590  if(locSCTrackDistanceMap.find(locTrack) != locSCTrackDistanceMap.end())
1591  {
1592  double locProjectedZ = locSCTrackDistanceMap[locTrack].second;
1593  dHistMap_SCPaddleVsZ_HasHit[locIsTimeBased]->Fill(locProjectedZ, locProjectedSCPaddleMap[locTrack].first);
1594  }
1595  }
1596  }
1597  else
1598  {
1599  dHistMap_PVsTheta_NoHit[SYS_START][locIsTimeBased]->Fill(locTheta, locP);
1600  dHistMap_PhiVsTheta_NoHit[SYS_START][locIsTimeBased]->Fill(locTheta, locPhi);
1601  if(locProjectedSCPaddleMap.find(locTrack) != locProjectedSCPaddleMap.end())
1602  {
1603  dHistMap_SCPaddleVsTheta_NoHit[locIsTimeBased]->Fill(locTheta, locProjectedSCPaddleMap[locTrack].first);
1604  if(locProjectedSCPaddleMap[locTrack].second)
1605  dHistMap_SCPaddle_BarrelRegion_NoHit[locIsTimeBased]->Fill(locProjectedSCPaddleMap[locTrack].first);
1606  else
1607  dHistMap_SCPaddle_NoseRegion_NoHit[locIsTimeBased]->Fill(locProjectedSCPaddleMap[locTrack].first);
1608  if(locSCTrackDistanceMap.find(locTrack) != locSCTrackDistanceMap.end())
1609  {
1610  double locProjectedZ = locSCTrackDistanceMap[locTrack].second;
1611  dHistMap_SCPaddleVsZ_NoHit[locIsTimeBased]->Fill(locProjectedZ, locProjectedSCPaddleMap[locTrack].first);
1612  }
1613  }
1614  }
1615  }
1616  }
1617 
1618  //Is-Matched to Something
1619  for(auto locMapPair : locBestTrackMap)
1620  {
1621  auto locTrack = locMapPair.second;
1622  double locTheta = locTrack->momentum().Theta()*180.0/TMath::Pi();
1623  double locP = locTrack->momentum().Mag();
1624  if(locDetectorMatches->Get_IsMatchedToHit(locTrack))
1625  dHistMap_TrackPVsTheta_HitMatch[locIsTimeBased]->Fill(locTheta, locP);
1626  else
1627  dHistMap_TrackPVsTheta_NoHitMatch[locIsTimeBased]->Fill(locTheta, locP);
1628  }
1629  }
1630  Unlock_Action(); //RELEASE ROOT LOCK!!
1631 }
1632 
1633 void DHistogramAction_DetectorPID::Initialize(JEventLoop* locEventLoop)
1634 {
1635  //Create any histograms/trees/etc. within a ROOT lock.
1636  //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
1637 
1638  //When creating a reaction-independent action, only modify member variables within a ROOT lock.
1639  //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
1640 
1641  string locHistName, locHistTitle, locParticleROOTName;
1642 
1643  string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
1644  if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
1645  gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
1646  if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
1647  gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
1648 
1649  //CREATE THE HISTOGRAMS
1650  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1651  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1652  {
1653  //So: Default tag is "", User can set it to something else
1654  //In here, if tag is "", get from gparms, if not, leave it alone
1655  //If gparms value does not exist, set it to (and use) "PreSelect"
1656  if(dTrackSelectionTag == "NotATag")
1657  dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
1658  if(dShowerSelectionTag == "NotATag")
1659  dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
1660 
1661  //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
1662  //If another thread has already created the folder, it just changes to it.
1664 
1665  //q = 0
1666  locParticleROOTName = ParticleName_ROOT(Gamma);
1667  //BCAL
1668  CreateAndChangeTo_Directory("BCAL", "BCAL");
1669 
1670  locHistName = "BetaVsP_q0";
1671  locHistTitle = "BCAL q^{0};Shower Energy (GeV);#beta";
1672  dHistMap_BetaVsP[SYS_BCAL][0] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1673 
1674  locHistName = "DeltaTVsShowerE_Photon";
1675  locHistTitle = string("BCAL q^{0}") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat_{BCAL - RF}");
1676  dHistMap_DeltaTVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1677 
1678 /*
1679  //Uncomment when ready!
1680  locHistName = "TimePullVsShowerE_Photon";
1681  locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
1682  dHistMap_TimePullVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1683 
1684  locHistName = "TimeFOMVsShowerE_Photon";
1685  locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
1686  dHistMap_TimeFOMVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1687 */
1688  gDirectory->cd("..");
1689 
1690  //FCAL
1691  CreateAndChangeTo_Directory("FCAL", "FCAL");
1692 
1693  locHistName = "BetaVsP_q0";
1694  locHistTitle = "FCAL q^{0};Shower Energy (GeV);#beta";
1695  dHistMap_BetaVsP[SYS_FCAL][0] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1696 
1697  locHistName = "DeltaTVsShowerE_Photon";
1698  locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat_{FCAL - RF}");
1699  dHistMap_DeltaTVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1700 
1701 /*
1702  //Uncomment when ready!
1703  locHistName = "TimePullVsShowerE_Photon";
1704  locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
1705  dHistMap_TimePullVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1706 
1707  locHistName = "TimeFOMVsShowerE_Photon";
1708  locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
1709  dHistMap_TimeFOMVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1710 */
1711  gDirectory->cd("..");
1712 
1713  //q +/-
1714  for(int locCharge = -1; locCharge <= 1; locCharge += 2)
1715  {
1716  string locParticleName = (locCharge == -1) ? "q-" : "q+";
1717  string locParticleROOTName = (locCharge == -1) ? "q^{-}" : "q^{+}";
1718 
1719  //SC
1720  CreateAndChangeTo_Directory("SC", "SC");
1721 
1722  locHistName = string("dEdXVsP_") + locParticleName;
1723  locHistTitle = locParticleROOTName + ";p (GeV/c);SC dE/dX (MeV/cm)";
1724  dHistMap_dEdXVsP[SYS_START][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1725 
1726  locHistName = string("BetaVsP_") + locParticleName;
1727  locHistTitle = string("SC ") + locParticleROOTName + string(";p (GeV/c);#beta");
1728  dHistMap_BetaVsP[SYS_START][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1729 
1730  gDirectory->cd("..");
1731 
1732  //TOF
1733  CreateAndChangeTo_Directory("TOF", "TOF");
1734 
1735  locHistName = string("dEdXVsP_") + locParticleName;
1736  locHistTitle = locParticleROOTName + ";p (GeV/c);TOF dE/dX (MeV/cm)";
1737  dHistMap_dEdXVsP[SYS_TOF][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1738 
1739  locHistName = string("BetaVsP_") + locParticleName;
1740  locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);#beta");
1741  dHistMap_BetaVsP[SYS_TOF][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1742 
1743  gDirectory->cd("..");
1744 
1745  //BCAL
1746  CreateAndChangeTo_Directory("BCAL", "BCAL");
1747 
1748  locHistName = string("BetaVsP_") + locParticleName;
1749  locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);#beta");
1750  dHistMap_BetaVsP[SYS_BCAL][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1751 
1752  locHistName = string("EOverPVsP_") + locParticleName;
1753  locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);E_{Shower}/p_{Track} (c);");
1754  dHistMap_BCALEOverPVsP[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1755 
1756  locHistName = string("EOverPVsTheta_") + locParticleName;
1757  locHistTitle = string("BCAL ") + locParticleROOTName + string(";#theta#circ;E_{Shower}/p_{Track} (c);");
1758  dHistMap_BCALEOverPVsTheta[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALThetaBins, dMinBCALTheta, dMaxBCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1759 
1760  gDirectory->cd("..");
1761 
1762  //FCAL
1763  CreateAndChangeTo_Directory("FCAL", "FCAL");
1764 
1765  locHistName = string("BetaVsP_") + locParticleName;
1766  locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#beta");
1767  dHistMap_BetaVsP[SYS_FCAL][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1768 
1769  locHistName = string("EOverPVsP_") + locParticleName;
1770  locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);E_{Shower}/p_{Track} (c);");
1771  dHistMap_FCALEOverPVsP[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1772 
1773  locHistName = string("EOverPVsTheta_") + locParticleName;
1774  locHistTitle = string("FCAL ") + locParticleROOTName + string(";#theta#circ;E_{Shower}/p_{Track} (c);");
1775  dHistMap_FCALEOverPVsTheta[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DFCALThetaBins, dMinFCALTheta, dMaxFCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1776 
1777  gDirectory->cd("..");
1778 
1779  //CDC
1780  CreateAndChangeTo_Directory("CDC", "CDC");
1781 
1782  locHistName = string("dEdXVsP_Int_") + locParticleName;
1783  locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx [Integral-based] (keV/cm)");
1784  dHistMap_dEdXVsP[SYS_CDC][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1785  locHistName = string("dEdXVsP_Amp_") + locParticleName;
1786  locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx [Amplitude-based] (keV/cm)");
1787  dHistMap_dEdXVsP[SYS_CDC_AMP][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1788 
1789  gDirectory->cd("..");
1790 
1791  //FDC
1792  CreateAndChangeTo_Directory("FDC", "FDC");
1793 
1794  locHistName = string("dEdXVsP_") + locParticleName;
1795  locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC dE/dx (keV/cm)");
1796  dHistMap_dEdXVsP[SYS_FDC][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1797 
1798  gDirectory->cd("..");
1799  }
1800 
1801  //delta's by PID
1802  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
1803  {
1804  Particle_t locPID = dFinalStatePIDs[loc_i];
1805  string locParticleName = ParticleType(locPID);
1806  string locParticleROOTName = ParticleName_ROOT(locPID);
1807 
1808  //SC
1809  CreateAndChangeTo_Directory("SC", "SC");
1810 
1811  locHistName = string("DeltadEdXVsP_") + locParticleName;
1812  locHistTitle = locParticleROOTName + " Candidates;p (GeV/c);SC #Delta(dE/dX) (MeV/cm)";
1813  dHistMap_DeltadEdXVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1814 
1815  locHistName = string("DeltaBetaVsP_") + locParticleName;
1816  locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1817  dHistMap_DeltaBetaVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1818 
1819  locHistName = string("DeltaTVsP_") + locParticleName;
1820  locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{SC - RF}");
1821  dHistMap_DeltaTVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1822 
1823 /*
1824  //Uncomment when ready!
1825  locHistName = string("TimePullVsP_") + locParticleName;
1826  locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1827  dHistMap_TimePullVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1828 
1829  locHistName = string("TimeFOMVsP_") + locParticleName;
1830  locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1831  dHistMap_TimeFOMVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1832 */
1833  gDirectory->cd("..");
1834 
1835  //TOF
1836  CreateAndChangeTo_Directory("TOF", "TOF");
1837 
1838  locHistName = string("DeltadEdXVsP_") + locParticleName;
1839  locHistTitle = locParticleROOTName + " Candidates;p (GeV/c);TOF #Delta(dE/dX) (MeV/cm)";
1840  dHistMap_DeltadEdXVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1841 
1842  locHistName = string("DeltaBetaVsP_") + locParticleName;
1843  locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1844  dHistMap_DeltaBetaVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1845 
1846  locHistName = string("DeltaTVsP_") + locParticleName;
1847  locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{TOF - RF}");
1848  dHistMap_DeltaTVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1849 
1850 /*
1851  //Uncomment when ready!
1852  locHistName = string("TimePullVsP_") + locParticleName;
1853  locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1854  dHistMap_TimePullVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1855 
1856  locHistName = string("TimeFOMVsP_") + locParticleName;
1857  locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1858  dHistMap_TimeFOMVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1859 */
1860  gDirectory->cd("..");
1861 
1862  //BCAL
1863  CreateAndChangeTo_Directory("BCAL", "BCAL");
1864 
1865  locHistName = string("DeltaBetaVsP_") + locParticleName;
1866  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1867  dHistMap_DeltaBetaVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1868 
1869  locHistName = string("DeltaTVsP_") + locParticleName;
1870  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{BCAL - RF}");
1871  dHistMap_DeltaTVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1872 
1873 /*
1874  //Uncomment when ready!
1875  locHistName = string("TimePullVsP_") + locParticleName;
1876  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1877  dHistMap_TimePullVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1878 
1879  locHistName = string("TimeFOMVsP_") + locParticleName;
1880  locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1881  dHistMap_TimeFOMVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1882 */
1883  gDirectory->cd("..");
1884 
1885  //FCAL
1886  CreateAndChangeTo_Directory("FCAL", "FCAL");
1887 
1888  locHistName = string("DeltaBetaVsP_") + locParticleName;
1889  locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1890  dHistMap_DeltaBetaVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1891 
1892  locHistName = string("DeltaTVsP_") + locParticleName;
1893  locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{FCAL - RF}");
1894  dHistMap_DeltaTVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1895 
1896 /*
1897  //Uncomment when ready!
1898  locHistName = string("TimePullVsP_") + locParticleName;
1899  locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1900  dHistMap_TimePullVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1901 
1902  locHistName = string("TimeFOMVsP_") + locParticleName;
1903  locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1904  dHistMap_TimeFOMVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1905 */
1906  gDirectory->cd("..");
1907 
1908  //CDC
1909  CreateAndChangeTo_Directory("CDC", "CDC");
1910 
1911  locHistName = string("DeltadEdXVsP_Int_") + locParticleName;
1912  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) [Integral-based] (keV/cm)");
1913  dHistMap_DeltadEdXVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1914 
1915  locHistName = string("DeltadEdXVsP_Amp_") + locParticleName;
1916  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) [Amplitude-based] (keV/cm)");
1917  dHistMap_DeltadEdXVsP[SYS_CDC_AMP][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1918 
1919 /*
1920  //Uncomment when ready!
1921  locHistName = string("dEdXPullVsP_") + locParticleName;
1922  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
1923  dHistMap_dEdXPullVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1924 
1925  locHistName = string("dEdXFOMVsP_") + locParticleName;
1926  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx PID Confidence Level");
1927  dHistMap_dEdXFOMVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1928 */
1929  gDirectory->cd("..");
1930 
1931  //FDC
1932  CreateAndChangeTo_Directory("FDC", "FDC");
1933 
1934  locHistName = string("DeltadEdXVsP_") + locParticleName;
1935  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX) (keV/cm)");
1936  dHistMap_DeltadEdXVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1937 
1938 /*
1939  //Uncomment when ready!
1940  locHistName = string("dEdXPullVsP_") + locParticleName;
1941  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
1942  dHistMap_dEdXPullVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1943 
1944  locHistName = string("dEdXFOMVsP_") + locParticleName;
1945  locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC dE/dx PID Confidence Level");
1946  dHistMap_dEdXFOMVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1947 */
1948  gDirectory->cd("..");
1949 
1950  // DIRC
1951  CreateAndChangeTo_Directory("DIRC", "DIRC");
1952 
1953  locHistName = string("NumPhotons_") + locParticleName;
1954  locHistTitle = locParticleROOTName + string("; DIRC NumPhotons");
1955  dHistMap_NumPhotons_DIRC[locPID] = new TH1I(locHistName.c_str(), locHistTitle.c_str(), dDIRCNumPhotonsBins, dDIRCMinNumPhotons, dDIRCMaxNumPhotons);
1956 
1957  locHistName = string("ThetaCVsP_") + locParticleName;
1958  locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC #theta_{C}");
1959  dHistMap_ThetaCVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCThetaCBins, dDIRCMinThetaC, dDIRCMaxThetaC);
1960 
1961  locHistName = string("Ldiff_kpiVsP_") + locParticleName;
1962  locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC L_{K}-L_{#pi}");
1963  dHistMap_Ldiff_kpiVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCLikelihoodBins, -1*dDIRCMaxLikelihood, dDIRCMaxLikelihood);
1964 
1965  locHistName = string("Ldiff_pkVsP_") + locParticleName;
1966  locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC L_{p}-L_{K}");
1967  dHistMap_Ldiff_pkVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCLikelihoodBins, -1*dDIRCMaxLikelihood, dDIRCMaxLikelihood);
1968 
1969  gDirectory->cd("..");
1970  }
1971 
1972  //Return to the base directory
1974  }
1975  japp->RootUnLock(); //RELEASE ROOT LOCK!!
1976 }
1977 
1978 bool DHistogramAction_DetectorPID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1979 {
1980  //Expect locParticleCombo to be NULL since this is a reaction-independent action.
1981 
1983  return true; //else double-counting!
1984 
1985  vector<const DChargedTrack*> locChargedTracks;
1986  locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
1987 
1988  vector<const DNeutralParticle*> locNeutralParticles;
1989  locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
1990 
1991  const DDetectorMatches* locDetectorMatches = NULL;
1992  locEventLoop->GetSingle(locDetectorMatches);
1993 
1994  const DEventRFBunch* locEventRFBunch = NULL;
1995  locEventLoop->GetSingle(locEventRFBunch);
1996 
1997  const DParticleID* locParticleID = NULL;
1998  locEventLoop->GetSingle(locParticleID);
1999 
2000  //FILL HISTOGRAMS
2001  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2002  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2003  Lock_Action(); //ACQUIRE ROOT LOCK!!
2004  {
2005  if(locEventRFBunch->dTimeSource != SYS_NULL) //only histogram beta for neutrals if the t0 is well known
2006  {
2007  for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
2008  {
2009  //doesn't matter which hypothesis you use for beta: t0 is from DEventVertex time
2010  const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
2011  double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
2012  const DNeutralShower* locNeutralShower = locNeutralParticles[loc_i]->dNeutralShower;
2013  double locShowerEnergy = locNeutralShower->dEnergy;
2014 
2015  double locDeltaT = locNeutralParticleHypothesis->time() - locEventRFBunch->dTime;
2016  if(locNeutralShower->dDetectorSystem == SYS_BCAL)
2017  {
2018  dHistMap_BetaVsP[SYS_BCAL][0]->Fill(locShowerEnergy, locBeta_Timing);
2019  dHistMap_DeltaTVsP[SYS_BCAL][Gamma]->Fill(locShowerEnergy, locDeltaT);
2020  }
2021  else
2022  {
2023  dHistMap_BetaVsP[SYS_FCAL][0]->Fill(locShowerEnergy, locBeta_Timing);
2024  dHistMap_DeltaTVsP[SYS_FCAL][Gamma]->Fill(locShowerEnergy, locDeltaT);
2025  }
2026  }
2027  }
2028 
2029  for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2030  {
2031  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestTrackingFOM();
2032  int locCharge = ParticleCharge(locChargedTrackHypothesis->PID());
2033  if(dHistMap_dEdXVsP[SYS_START].find(locCharge) == dHistMap_dEdXVsP[SYS_START].end())
2034  continue;
2035 
2036  double locStartTime = locParticleID->Calc_PropagatedRFTime(locChargedTrackHypothesis, locEventRFBunch);
2037  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
2038 
2039  Particle_t locPID = locChargedTrackHypothesis->PID();
2040  double locP = locTrackTimeBased->momentum().Mag();
2041  double locTheta = locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi();
2042 
2043  //if RF time is indeterminate, start time will be NaN
2044  auto locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
2045  auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
2046  auto locTOFHitMatchParams = locChargedTrackHypothesis->Get_TOFHitMatchParams();
2047  auto locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
2048  auto locDIRCMatchParams = locChargedTrackHypothesis->Get_DIRCMatchParams();
2049 
2050  if(locSCHitMatchParams != NULL)
2051  {
2052  dHistMap_dEdXVsP[SYS_START][locCharge]->Fill(locP, locSCHitMatchParams->dEdx*1.0E3);
2053  if((locEventRFBunch->dTimeSource != SYS_START) && (locEventRFBunch->dNumParticleVotes > 1))
2054  {
2055  //If SC was used for RF time, don't compute delta-beta
2056  double locBeta_Timing = locSCHitMatchParams->dPathLength/(29.9792458*(locSCHitMatchParams->dHitTime - locChargedTrackHypothesis->t0()));
2057  dHistMap_BetaVsP[SYS_START][locCharge]->Fill(locP, locBeta_Timing);
2059  {
2060  double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
2061  dHistMap_DeltaBetaVsP[SYS_START][locPID]->Fill(locP, locDeltaBeta);
2062  double locDeltaT = locSCHitMatchParams->dHitTime - locSCHitMatchParams->dFlightTime - locStartTime;
2063  dHistMap_DeltaTVsP[SYS_START][locPID]->Fill(locP, locDeltaT);
2064  }
2065  }
2067  {
2068  double locdx = locSCHitMatchParams->dHitEnergy/locSCHitMatchParams->dEdx;
2069  double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
2070  locParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
2071  dHistMap_DeltadEdXVsP[SYS_START][locPID]->Fill(locP, (locSCHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
2072  }
2073  }
2074  if(locTOFHitMatchParams != NULL)
2075  {
2076  dHistMap_dEdXVsP[SYS_TOF][locCharge]->Fill(locP, locTOFHitMatchParams->dEdx*1.0E3);
2077  if(locEventRFBunch->dNumParticleVotes > 1)
2078  {
2079  double locBeta_Timing = locTOFHitMatchParams->dPathLength/(29.9792458*(locTOFHitMatchParams->dHitTime - locChargedTrackHypothesis->t0()));
2080  dHistMap_BetaVsP[SYS_TOF][locCharge]->Fill(locP, locBeta_Timing);
2081  if(dHistMap_DeltaBetaVsP[SYS_TOF].find(locPID) != dHistMap_DeltaBetaVsP[SYS_TOF].end())
2082  {
2083  double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
2084  dHistMap_DeltaBetaVsP[SYS_TOF][locPID]->Fill(locP, locDeltaBeta);
2085  double locDeltaT = locTOFHitMatchParams->dHitTime - locTOFHitMatchParams->dFlightTime - locStartTime;
2086  dHistMap_DeltaTVsP[SYS_TOF][locPID]->Fill(locP, locDeltaT);
2087  }
2088  }
2089  if(dHistMap_DeltadEdXVsP[SYS_TOF].find(locPID) != dHistMap_DeltadEdXVsP[SYS_TOF].end())
2090  {
2091  double locdx = locTOFHitMatchParams->dHitEnergy/locTOFHitMatchParams->dEdx;
2092  double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
2093  locParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
2094  dHistMap_DeltadEdXVsP[SYS_TOF][locPID]->Fill(locP, (locTOFHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
2095  }
2096  }
2097  if(locBCALShowerMatchParams != NULL)
2098  {
2099  const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
2100  double locEOverP = locBCALShower->E/locP;
2101  dHistMap_BCALEOverPVsP[locCharge]->Fill(locP, locEOverP);
2102  dHistMap_BCALEOverPVsTheta[locCharge]->Fill(locTheta, locEOverP);
2103  if(locEventRFBunch->dNumParticleVotes > 1)
2104  {
2105  double locBeta_Timing = locBCALShowerMatchParams->dPathLength/(29.9792458*(locBCALShower->t - locChargedTrackHypothesis->t0()));
2106  dHistMap_BetaVsP[SYS_BCAL][locCharge]->Fill(locP, locBeta_Timing);
2107  if(dHistMap_DeltaBetaVsP[SYS_BCAL].find(locPID) != dHistMap_DeltaBetaVsP[SYS_BCAL].end())
2108  {
2109  double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
2110  dHistMap_DeltaBetaVsP[SYS_BCAL][locPID]->Fill(locP, locDeltaBeta);
2111  double locDeltaT = locBCALShower->t - locBCALShowerMatchParams->dFlightTime - locStartTime;
2112  dHistMap_DeltaTVsP[SYS_BCAL][locPID]->Fill(locP, locDeltaT);
2113  }
2114  }
2115  }
2116  if(locFCALShowerMatchParams != NULL)
2117  {
2118  const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
2119  double locEOverP = locFCALShower->getEnergy()/locP;
2120  dHistMap_FCALEOverPVsP[locCharge]->Fill(locP, locEOverP);
2121  dHistMap_FCALEOverPVsTheta[locCharge]->Fill(locTheta, locEOverP);
2122  if(locEventRFBunch->dNumParticleVotes > 1)
2123  {
2124  double locBeta_Timing = locFCALShowerMatchParams->dPathLength/(29.9792458*(locFCALShower->getTime() - locChargedTrackHypothesis->t0()));
2125  dHistMap_BetaVsP[SYS_FCAL][locCharge]->Fill(locP, locBeta_Timing);
2126  if(dHistMap_DeltaBetaVsP[SYS_FCAL].find(locPID) != dHistMap_DeltaBetaVsP[SYS_FCAL].end())
2127  {
2128  double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
2129  dHistMap_DeltaBetaVsP[SYS_FCAL][locPID]->Fill(locP, locDeltaBeta);
2130  double locDeltaT = locFCALShower->getTime() - locFCALShowerMatchParams->dFlightTime - locStartTime;
2131  dHistMap_DeltaTVsP[SYS_FCAL][locPID]->Fill(locP, locDeltaT);
2132  }
2133  }
2134  }
2135  if(locDIRCMatchParams != NULL && (dHistMap_NumPhotons_DIRC.find(locPID) != dHistMap_NumPhotons_DIRC.end())) {
2136  int locNumPhotons_DIRC = locDIRCMatchParams->dNPhotons;
2137  double locThetaC_DIRC = locDIRCMatchParams->dThetaC;
2138  dHistMap_NumPhotons_DIRC[locPID]->Fill(locNumPhotons_DIRC);
2139  dHistMap_ThetaCVsP_DIRC[locPID]->Fill(locP, locThetaC_DIRC);
2140  double locLpi_DIRC = locDIRCMatchParams->dLikelihoodPion;
2141  double locLk_DIRC = locDIRCMatchParams->dLikelihoodKaon;
2142  double locLp_DIRC = locDIRCMatchParams->dLikelihoodProton;
2143  dHistMap_Ldiff_kpiVsP_DIRC[locPID]->Fill(locP, locLk_DIRC-locLpi_DIRC);
2144  dHistMap_Ldiff_pkVsP_DIRC[locPID]->Fill(locP, locLp_DIRC-locLk_DIRC);
2145  }
2146 
2147  if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC > 0)
2148  {
2149  dHistMap_dEdXVsP[SYS_CDC][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_CDC*1.0E6);
2150  dHistMap_dEdXVsP[SYS_CDC_AMP][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_CDC_amp*1.0E6);
2151  if(dHistMap_DeltadEdXVsP[SYS_CDC].find(locPID) != dHistMap_DeltadEdXVsP[SYS_CDC].end())
2152  {
2153  double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_CDC, true);
2154  dHistMap_DeltadEdXVsP[SYS_CDC][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_CDC - locProbabledEdx)*1.0E6);
2155  }
2157  {
2158  double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_CDC_amp, true);
2159  dHistMap_DeltadEdXVsP[SYS_CDC_AMP][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_CDC_amp - locProbabledEdx)*1.0E6);
2160  }
2161  }
2162  if(locTrackTimeBased->dNumHitsUsedFordEdx_FDC > 0)
2163  {
2164  dHistMap_dEdXVsP[SYS_FDC][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_FDC*1.0E6);
2165  if(dHistMap_DeltadEdXVsP[SYS_FDC].find(locPID) != dHistMap_DeltadEdXVsP[SYS_FDC].end())
2166  {
2167  double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_FDC, false);
2168  dHistMap_DeltadEdXVsP[SYS_FDC][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_FDC - locProbabledEdx)*1.0E6);
2169  }
2170  }
2171  }
2172  }
2173  Unlock_Action(); //RELEASE ROOT LOCK!!
2174 
2175  return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2176 }
2177 
2178 void DHistogramAction_Neutrals::Initialize(JEventLoop* locEventLoop)
2179 {
2180  //Create any histograms/trees/etc. within a ROOT lock.
2181  //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
2182 
2183  //When creating a reaction-independent action, only modify member variables within a ROOT lock.
2184  //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
2185 
2186  string locHistName;
2187 
2188  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
2189  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
2190  double locTargetZCenter = 0.0;
2191  locGeometry->GetTargetZ(locTargetZCenter);
2192 
2193  //CREATE THE HISTOGRAMS
2194  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2195  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2196  {
2197  //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
2198  //If another thread has already created the folder, it just changes to it.
2200 
2201  if(dTargetCenter.Z() < -9.8E9)
2202  dTargetCenter.SetXYZ(0.0, 0.0, locTargetZCenter);
2203 
2204  //BCAL
2205  locHistName = "BCALTrackDOCA";
2206  dHist_BCALTrackDOCA = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower Distance to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
2207  locHistName = "BCALTrackDeltaPhi";
2208  dHist_BCALTrackDeltaPhi = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #Delta#phi#circ to Nearest Track", dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
2209  locHistName = "BCALTrackDeltaZ";
2210  dHist_BCALTrackDeltaZ = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #DeltaZ to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
2211  locHistName = "BCALNeutralShowerEnergy";
2212  dHist_BCALNeutralShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Neutral Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2213  locHistName = "BCALNeutralShowerDeltaT";
2214  dHist_BCALNeutralShowerDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Neutral Shower #Deltat (Propagated - RF) (ns)", dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
2215  locHistName = "BCALNeutralShowerDeltaTVsE";
2216  dHist_BCALNeutralShowerDeltaTVsE = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Neutral Shower Energy (GeV);BCAL Neutral Shower #Deltat (ns)", dNum2DShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
2217  locHistName = "BCALNeutralShowerDeltaTVsZ";
2218  dHist_BCALNeutralShowerDeltaTVsZ = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Neutral Shower Z (cm);BCAL Neutral Shower #Deltat (ns)", dNum2DBCALZBins, 0.0, 450.0, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
2219 
2220  //FCAL
2221  locHistName = "FCALTrackDOCA";
2222  dHist_FCALTrackDOCA = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Shower Distance to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
2223  locHistName = "FCALNeutralShowerEnergy";
2224  dHist_FCALNeutralShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Neutral Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2225  locHistName = "FCALNeutralShowerDeltaT";
2226  dHist_FCALNeutralShowerDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Neutral Shower #Deltat (Propagated - RF) (ns)", dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
2227  locHistName = "FCALNeutralShowerDeltaTVsE";
2228  dHist_FCALNeutralShowerDeltaTVsE = GetOrCreate_Histogram<TH2I>(locHistName, ";FCAL Neutral Shower Energy (GeV);FCAL Neutral Shower #Deltat (ns)", dNum2DShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
2229 
2230  //Return to the base directory
2232  }
2233  japp->RootUnLock(); //RELEASE ROOT LOCK!!
2234 }
2235 
2236 bool DHistogramAction_Neutrals::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2237 {
2238  //Expect locParticleCombo to be NULL since this is a reaction-independent action.
2239 
2241  return true; //else double-counting!
2242 
2243  vector<const DNeutralShower*> locNeutralShowers;
2244  locEventLoop->Get(locNeutralShowers);
2245 
2246  vector<const DTrackTimeBased*> locTrackTimeBasedVector;
2247  locEventLoop->Get(locTrackTimeBasedVector);
2248 
2249  const DDetectorMatches* locDetectorMatches = NULL;
2250  locEventLoop->GetSingle(locDetectorMatches);
2251 
2252  vector<const DEventRFBunch*> locEventRFBunches;
2253  locEventLoop->Get(locEventRFBunches);
2254  double locStartTime = locEventRFBunches.empty() ? 0.0 : locEventRFBunches[0]->dTime;
2255 
2256  //FILL HISTOGRAMS
2257  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2258  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2259  Lock_Action(); //ACQUIRE ROOT LOCK!!
2260  {
2261  for(size_t loc_i = 0; loc_i < locNeutralShowers.size(); ++loc_i)
2262  {
2263  //assume is photon
2264  double locPathLength = (locNeutralShowers[loc_i]->dSpacetimeVertex.Vect() - dTargetCenter).Mag();
2265  double locDeltaT = locNeutralShowers[loc_i]->dSpacetimeVertex.T() - locPathLength/29.9792458 - locStartTime;
2266 
2267  if(locNeutralShowers[loc_i]->dDetectorSystem == SYS_FCAL)
2268  {
2269  const DFCALShower* locFCALShower = NULL;
2270  locNeutralShowers[loc_i]->GetSingle(locFCALShower);
2271 
2272  double locDistance = 9.9E9;
2273  if(locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistance))
2274  dHist_FCALTrackDOCA->Fill(locDistance);
2275 
2276  dHist_FCALNeutralShowerEnergy->Fill(locNeutralShowers[loc_i]->dEnergy);
2277  dHist_FCALNeutralShowerDeltaT->Fill(locDeltaT);
2278  dHist_FCALNeutralShowerDeltaTVsE->Fill(locNeutralShowers[loc_i]->dEnergy, locDeltaT);
2279  }
2280  else
2281  {
2282  const DBCALShower* locBCALShower = NULL;
2283  locNeutralShowers[loc_i]->GetSingle(locBCALShower);
2284 
2285  double locDistance = 9.9E9, locDeltaPhi = 9.9E9, locDeltaZ = 9.9E9;
2286  if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDistance))
2287  dHist_BCALTrackDOCA->Fill(locDistance);
2288  if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDeltaPhi, locDeltaZ))
2289  {
2290  dHist_BCALTrackDeltaPhi->Fill(180.0*locDeltaPhi/TMath::Pi());
2291  dHist_BCALTrackDeltaZ->Fill(locDeltaZ);
2292  }
2293 
2294  dHist_BCALNeutralShowerEnergy->Fill(locNeutralShowers[loc_i]->dEnergy);
2295  dHist_BCALNeutralShowerDeltaT->Fill(locDeltaT);
2296  dHist_BCALNeutralShowerDeltaTVsE->Fill(locNeutralShowers[loc_i]->dEnergy, locDeltaT);
2297  dHist_BCALNeutralShowerDeltaTVsZ->Fill(locNeutralShowers[loc_i]->dSpacetimeVertex.Z(), locDeltaT);
2298  }
2299  }
2300  }
2301  Unlock_Action(); //RELEASE ROOT LOCK!!
2302 
2303  return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2304 }
2305 
2306 
2308 {
2309  //Create any histograms/trees/etc. within a ROOT lock.
2310  //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
2311 
2312  //When creating a reaction-independent action, only modify member variables within a ROOT lock.
2313  //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
2314 
2315  string locHistName, locHistTitle;
2316 
2317  vector<const DMCThrown*> locMCThrowns;
2318  locEventLoop->Get(locMCThrowns);
2319 
2320  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
2321  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
2322  double locTargetZCenter = 0.0;
2323  locGeometry->GetTargetZ(locTargetZCenter);
2324 
2325  string locTrackSelectionTag = "NotATag";
2326  if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2327  gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2328 
2329  //CREATE THE HISTOGRAMS
2330  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2331  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2332  {
2333  //So: Default tag is "", User can set it to something else
2334  //In here, if tag is "", get from gparms, if not, leave it alone
2335  //If gparms value does not exist, set it to (and use) "PreSelect"
2336  if(dTrackSelectionTag == "NotATag")
2337  dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2338 
2339  //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
2340  //If another thread has already created the folder, it just changes to it.
2342 
2343  if(dTargetCenterZ < -9.8E9)
2344  dTargetCenterZ = locTargetZCenter; //only set if not already set
2345 
2346  //Track Matched to Hit
2347  for(int locTruePIDFlag = 0; locTruePIDFlag < 2; ++locTruePIDFlag)
2348  {
2349  if(locMCThrowns.empty() && (locTruePIDFlag == 1))
2350  continue; //not a simulated event: don't histogram thrown info!
2351 
2352  string locDirName = (locTruePIDFlag == 1) ? "TruePID" : "ReconstructedPID";
2353  CreateAndChangeTo_Directory(locDirName.c_str(), locDirName.c_str());
2354 
2355  //By PID
2356  for(int loc_i = -2; loc_i < int(dTrackingPIDs.size()); ++loc_i) //-2 = q-, -1 = q+
2357  {
2358  string locParticleName, locParticleROOTName;
2359  int locPID = loc_i;
2360  if(loc_i == -2)
2361  {
2362  locParticleName = "q-";
2363  locParticleROOTName = "#it{q}^{-}";
2364  }
2365  else if(loc_i == -1)
2366  {
2367  locParticleName = "q+";
2368  locParticleROOTName = "#it{q}^{+}";
2369  }
2370  else
2371  {
2372  locParticleName = ParticleType(dTrackingPIDs[loc_i]);
2373  locParticleROOTName = ParticleName_ROOT(dTrackingPIDs[loc_i]);
2374  locPID = int(dTrackingPIDs[loc_i]);
2375  }
2376  CreateAndChangeTo_Directory(locParticleName, locParticleName);
2377  pair<int, bool> locPIDPair(locPID, bool(locTruePIDFlag));
2378 
2379  //BCAL
2380  locHistName = "BCALShowerEnergy";
2381  locHistTitle = locParticleROOTName + ";BCAL Shower Energy (GeV)";
2382  dHistMap_BCALShowerEnergy[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2383 
2384  locHistName = "BCALShowerTrackDepth";
2385  locHistTitle = locParticleROOTName + ";BCAL Shower Track Depth (cm)";
2386  dHistMap_BCALShowerTrackDepth[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2387 
2388  locHistName = "BCALShowerTrackDepthVsP";
2389  locHistTitle = locParticleROOTName + ";p (GeV/c);BCAL Shower Track Depth (cm)";
2390  dHistMap_BCALShowerTrackDepthVsP[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2391 
2392  //FCAL
2393  locHistName = "FCALShowerEnergy";
2394  locHistTitle = locParticleROOTName + ";FCAL Shower Energy (GeV)";
2395  dHistMap_FCALShowerEnergy[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2396 
2397  locHistName = "FCALShowerTrackDepth";
2398  locHistTitle = locParticleROOTName + ";FCAL Shower Track Depth (cm)";
2399  dHistMap_FCALShowerTrackDepth[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2400 
2401  locHistName = "FCALShowerTrackDepthVsP";
2402  locHistTitle = locParticleROOTName + ";p (GeV/c);FCAL Shower Track Depth (cm)";
2403  dHistMap_FCALShowerTrackDepthVsP[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2404 
2405  //SC
2406  locHistName = "SCEnergyVsTheta";
2407  locHistTitle = locParticleROOTName + ";#theta#circ;SC Point Energy (MeV)";
2408  dHistMap_SCEnergyVsTheta[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
2409 
2410  locHistName = "SCPhiVsTheta";
2411  locHistTitle = locParticleROOTName + ";#theta#circ;#phi#circ";
2412  dHistMap_SCPhiVsTheta[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
2413 
2414  gDirectory->cd(".."); //end of PID
2415  }
2416  gDirectory->cd(".."); //end of true/recon PID
2417  }
2418 
2419  //Return to the base directory
2421  }
2422  japp->RootUnLock(); //RELEASE ROOT LOCK!!
2423 }
2424 
2425 bool DHistogramAction_DetectorMatchParams::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2426 {
2427  //Expect locParticleCombo to be NULL since this is a reaction-independent action.
2428 
2430  return true; //else double-counting!
2431 
2432  vector<const DMCThrown*> locMCThrowns;
2433  locEventLoop->Get(locMCThrowns);
2434 
2435  Fill_Hists(locEventLoop, false);
2436  if(!locMCThrowns.empty())
2437  Fill_Hists(locEventLoop, true);
2438 
2439  return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2440 }
2441 
2442 void DHistogramAction_DetectorMatchParams::Fill_Hists(JEventLoop* locEventLoop, bool locUseTruePIDFlag)
2443 {
2444  vector<const DChargedTrack*> locChargedTracks;
2445  locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
2446 
2447  vector<const DMCThrownMatching*> locMCThrownMatchingVector;
2448  locEventLoop->Get(locMCThrownMatchingVector);
2449 
2450  const DEventRFBunch* locEventRFBunch = NULL;
2451  locEventLoop->GetSingle(locEventRFBunch);
2452 
2453  //FILL HISTOGRAMS
2454  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2455  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2456  Lock_Action(); //ACQUIRE ROOT LOCK!!
2457  {
2458  for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2459  {
2460  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2461 
2462  if(locUseTruePIDFlag && (!locMCThrownMatchingVector.empty()))
2463  {
2464  double locMatchFOM = 0.0;
2465  const DMCThrown* locMCThrown = locMCThrownMatchingVector[0]->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);
2466  if(locMCThrown == NULL)
2467  continue;
2468  //OK, have the thrown. Now, grab the best charged track hypothesis to get the best matching
2469  locChargedTrackHypothesis = locMCThrownMatchingVector[0]->Get_MatchingChargedHypothesis(locMCThrown, locMatchFOM);
2470  }
2471 
2472  pair<int, bool> locPIDPair(int(locChargedTrackHypothesis->PID()), locUseTruePIDFlag);
2473  bool locDisregardPIDFlag = (dHistMap_BCALShowerEnergy.find(locPIDPair) == dHistMap_BCALShowerEnergy.end());
2474  int locQIndex = (locChargedTrackHypothesis->charge() > 0.0) ? -1 : -2;
2475  pair<int, bool> locChargePair(locQIndex, locUseTruePIDFlag);
2476 
2477  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2478  auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
2479  auto locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
2480  auto locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
2481 
2482  //BCAL
2483  if(locBCALShowerMatchParams != NULL)
2484  {
2485  const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
2486  dHistMap_BCALShowerEnergy[locChargePair]->Fill(locBCALShower->E);
2487  dHistMap_BCALShowerTrackDepth[locChargePair]->Fill(locBCALShowerMatchParams->dx);
2488  dHistMap_BCALShowerTrackDepthVsP[locChargePair]->Fill(locMomentum.Mag(), locBCALShowerMatchParams->dx);
2489 
2490  if(!locDisregardPIDFlag)
2491  {
2492  dHistMap_BCALShowerEnergy[locPIDPair]->Fill(locBCALShower->E);
2493  dHistMap_BCALShowerTrackDepth[locPIDPair]->Fill(locBCALShowerMatchParams->dx);
2494  dHistMap_BCALShowerTrackDepthVsP[locPIDPair]->Fill(locMomentum.Mag(), locBCALShowerMatchParams->dx);
2495  }
2496  }
2497 
2498  //FCAL
2499  if(locFCALShowerMatchParams != NULL)
2500  {
2501  const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
2502  dHistMap_FCALShowerEnergy[locChargePair]->Fill(locFCALShower->getEnergy());
2503  dHistMap_FCALShowerTrackDepth[locChargePair]->Fill(locFCALShowerMatchParams->dx);
2504  dHistMap_FCALShowerTrackDepthVsP[locChargePair]->Fill(locMomentum.Mag(), locFCALShowerMatchParams->dx);
2505 
2506  if(!locDisregardPIDFlag)
2507  {
2508  dHistMap_FCALShowerEnergy[locPIDPair]->Fill(locFCALShower->getEnergy());
2509  dHistMap_FCALShowerTrackDepth[locPIDPair]->Fill(locFCALShowerMatchParams->dx);
2510  dHistMap_FCALShowerTrackDepthVsP[locPIDPair]->Fill(locMomentum.Mag(), locFCALShowerMatchParams->dx);
2511  }
2512  }
2513 
2514  //SC
2515  if(locSCHitMatchParams != NULL)
2516  {
2517  dHistMap_SCEnergyVsTheta[locChargePair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locSCHitMatchParams->dHitEnergy*1.0E3);
2518  dHistMap_SCPhiVsTheta[locChargePair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locMomentum.Phi()*180.0/TMath::Pi());
2519 
2520  if(!locDisregardPIDFlag)
2521  {
2522  dHistMap_SCEnergyVsTheta[locPIDPair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locSCHitMatchParams->dHitEnergy*1.0E3);
2523  dHistMap_SCPhiVsTheta[locPIDPair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locMomentum.Phi()*180.0/TMath::Pi());
2524  }
2525  }
2526  }
2527  }
2528  Unlock_Action(); //RELEASE ROOT LOCK!!
2529 }
2530 
2531 void DHistogramAction_EventVertex::Initialize(JEventLoop* locEventLoop)
2532 {
2533  string locHistName, locHistTitle;
2534 
2535  string locTrackSelectionTag = "NotATag";
2536  if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2537  gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2538 
2539  //CREATE THE HISTOGRAMS
2540  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2541  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2542  {
2543  //So: Default tag is "", User can set it to something else
2544  //In here, if tag is "", get from gparms, if not, leave it alone
2545  //If gparms value does not exist, set it to (and use) "PreSelect"
2546  if(dTrackSelectionTag == "NotATag")
2547  dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2548 
2550 
2551  // Event RF Bunch Time
2552  locHistName = "RFTrackDeltaT";
2553  locHistTitle = ";#Deltat_{RF - Track} (ns)";
2554  dRFTrackDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumRFTBins, -3.0, 3.0);
2555 
2556  // ALL EVENTS
2557  CreateAndChangeTo_Directory("AllEvents", "AllEvents");
2558 
2559  // Event Vertex-Z
2560  locHistName = "EventVertexZ";
2561  locHistTitle = ";Event Vertex-Z (cm)";
2562  dEventVertexZ_AllEvents = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2563 
2564  // Event Vertex-Y Vs Vertex-X
2565  locHistName = "EventVertexYVsX";
2566  locHistTitle = ";Event Vertex-X (cm);Event Vertex-Y (cm)";
2567  dEventVertexYVsX_AllEvents = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2568 
2569  // Event Vertex-T
2570  locHistName = "EventVertexT";
2571  locHistTitle = ";Event Vertex Time (ns)";
2572  dEventVertexT_AllEvents = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2573 
2574  gDirectory->cd("..");
2575 
2576 
2577  // 2+ Good Tracks
2578  CreateAndChangeTo_Directory("2+GoodTracks", "2+GoodTracks");
2579 
2580  // Event Vertex-Z
2581  locHistName = "EventVertexZ";
2582  locHistTitle = ";Event Vertex-Z (cm)";
2583  dEventVertexZ_2OrMoreGoodTracks = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2584 
2585  // Event Vertex-Y Vs Vertex-X
2586  locHistName = "EventVertexYVsX";
2587  locHistTitle = ";Event Vertex-X (cm);Event Vertex-Y (cm)";
2588  dEventVertexYVsX_2OrMoreGoodTracks = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2589 
2590  // Event Vertex-T
2591  locHistName = "EventVertexT";
2592  locHistTitle = ";Event Vertex Time (ns)";
2593  dEventVertexT_2OrMoreGoodTracks = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2594 
2595  // Confidence Level
2596  locHistName = "ConfidenceLevel";
2597  dHist_KinFitConfidenceLevel = GetOrCreate_Histogram<TH1I>(locHistName, "Event Vertex Kinematic Fit;Confidence Level;# Events", dNumConfidenceLevelBins, 0.0, 1.0);
2598 
2599  //final particle pulls
2600  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2601  {
2602  Particle_t locPID = dFinalStatePIDs[loc_i];
2603  string locParticleDirName = string("Pulls_") + string(ParticleType(locPID));
2604  string locParticleROOTName = ParticleName_ROOT(locPID);
2605  CreateAndChangeTo_Directory(locParticleDirName, locParticleDirName);
2606 
2607  //Px Pull
2608  locHistName = "Pull_Px";
2609  locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{x} Pull;# Events");
2610  dHistMap_KinFitPulls[locPID][d_PxPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2611 
2612  //Py Pull
2613  locHistName = "Pull_Py";
2614  locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{y} Pull;# Events");
2615  dHistMap_KinFitPulls[locPID][d_PyPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2616 
2617  //Pz Pull
2618  locHistName = "Pull_Pz";
2619  locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{z} Pull;# Events");
2620  dHistMap_KinFitPulls[locPID][d_PzPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2621 
2622  //Xx Pull
2623  locHistName = "Pull_Xx";
2624  locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{x} Pull;# Events");
2625  dHistMap_KinFitPulls[locPID][d_XxPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2626 
2627  //Xy Pull
2628  locHistName = "Pull_Xy";
2629  locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{y} Pull;# Events");
2630  dHistMap_KinFitPulls[locPID][d_XyPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2631 
2632  //Xz Pull
2633  locHistName = "Pull_Xz";
2634  locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{z} Pull;# Events");
2635  dHistMap_KinFitPulls[locPID][d_XzPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2636 
2637  gDirectory->cd("..");
2638  }
2639 
2640  gDirectory->cd("..");
2641 
2642  //Return to the base directory
2644  }
2645  japp->RootUnLock(); //RELEASE ROOT LOCK!!
2646 }
2647 
2648 bool DHistogramAction_EventVertex::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2649 {
2651  return true; //else double-counting!
2652 
2653  const DVertex* locVertex = NULL;
2654  locEventLoop->GetSingle(locVertex);
2655 
2656  vector<const DChargedTrack*> locChargedTracks;
2657  locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
2658 
2659  const DDetectorMatches* locDetectorMatches = NULL;
2660  locEventLoop->GetSingle(locDetectorMatches);
2661 
2662  const DParticleID* locParticleID = NULL;
2663  locEventLoop->GetSingle(locParticleID);
2664 
2665  const DEventRFBunch* locEventRFBunch = NULL;
2666  locEventLoop->GetSingle(locEventRFBunch);
2667 
2668  //Make sure that brun() is called (to get rf period) before using.
2669  //Cannot call JEventLoop->Get() because object may be in datastream (REST), bypassing factory brun() call.
2670  //Must do here rather than in Initialize() function because this object is shared by all threads (which each have their own factory)
2671  DRFTime_factory* locRFTimeFactory = static_cast<DRFTime_factory*>(locEventLoop->GetFactory("DRFTime"));
2672  if(!locRFTimeFactory->brun_was_called())
2673  {
2674  locRFTimeFactory->brun(locEventLoop, locEventLoop->GetJEvent().GetRunNumber());
2675  locRFTimeFactory->Set_brun_called();
2676  }
2677 
2678  //Get time-based tracks: use best PID FOM
2679  //Note that these may not be the PIDs that were used in the fit!!!
2680  //e.g. for a DTrackTimeBased the proton hypothesis has the highest tracking FOM, so it is used in the fit, but the pi+ PID has the highest PID FOM
2681  vector<const DTrackTimeBased*> locTrackTimeBasedVector;
2682  for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2683  {
2684  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2685  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
2686  if(locTrackTimeBased != NULL)
2687  locTrackTimeBasedVector.push_back(locTrackTimeBased);
2688  }
2689 
2690  //Event Vertex
2691  //FILL HISTOGRAMS
2692  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2693  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2694  Lock_Action();
2695  {
2696  for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2697  {
2698  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2699  double locPropagatedRFTime = locParticleID->Calc_PropagatedRFTime(locChargedTrackHypothesis, locEventRFBunch);
2700  double locShiftedRFTime = locRFTimeFactory->Step_TimeToNearInputTime(locPropagatedRFTime, locChargedTrackHypothesis->time());
2701  double locDeltaT = locShiftedRFTime - locChargedTrackHypothesis->time();
2702  dRFTrackDeltaT->Fill(locDeltaT);
2703  }
2704  dEventVertexZ_AllEvents->Fill(locVertex->dSpacetimeVertex.Z());
2705  dEventVertexYVsX_AllEvents->Fill(locVertex->dSpacetimeVertex.X(), locVertex->dSpacetimeVertex.Y());
2706  dEventVertexT_AllEvents->Fill(locVertex->dSpacetimeVertex.T());
2707 
2708  if(locChargedTracks.size() >= 2)
2709  {
2710  dEventVertexZ_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.Z());
2711  dEventVertexYVsX_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.X(), locVertex->dSpacetimeVertex.Y());
2712  dEventVertexT_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.T());
2713  }
2714  }
2715  Unlock_Action();
2716 
2717  if(locVertex->dKinFitNDF == 0)
2718  return true; //kin fit not performed or didn't converge: no results to histogram
2719 
2720  double locConfidenceLevel = TMath::Prob(locVertex->dKinFitChiSq, locVertex->dKinFitNDF);
2721 
2722  //FILL HISTOGRAMS
2723  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2724  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2725  Lock_Action();
2726  {
2727  dHist_KinFitConfidenceLevel->Fill(locConfidenceLevel);
2728 
2729  //pulls
2730  if(locConfidenceLevel > dPullHistConfidenceLevelCut)
2731  {
2732  for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
2733  {
2734  const DKinematicData* locKinematicData = static_cast<const DKinematicData*>(locTrackTimeBasedVector[loc_i]);
2735  Particle_t locPID = locKinematicData->PID();
2736  if(dHistMap_KinFitPulls.find(locPID) == dHistMap_KinFitPulls.end())
2737  continue; //PID not histogrammed
2738 
2739  map<const JObject*, map<DKinFitPullType, double> >::const_iterator locParticleIterator = locVertex->dKinFitPulls.find(locKinematicData);
2740  if(locParticleIterator == locVertex->dKinFitPulls.end())
2741  continue;
2742 
2743  const map<DKinFitPullType, double>& locPullMap = locParticleIterator->second;
2744  map<DKinFitPullType, double>::const_iterator locPullIterator = locPullMap.begin();
2745  for(; locPullIterator != locPullMap.end(); ++locPullIterator)
2746  dHistMap_KinFitPulls[locPID][locPullIterator->first]->Fill(locPullIterator->second);
2747  }
2748  }
2749  }
2750  Unlock_Action();
2751 
2752  return true;
2753 }
2754 
2756 {
2757  string locHistName, locHistTitle, locParticleName, locParticleROOTName;
2758  Particle_t locPID;
2759 
2760  string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
2761  if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2762  gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2763  if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
2764  gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
2765 
2766  //CREATE THE HISTOGRAMS
2767  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2768  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2769  {
2770  //So: Default tag is "", User can set it to something else
2771  //In here, if tag is "", get from gparms, if not, leave it alone
2772  //If gparms value does not exist, set it to (and use) "PreSelect"
2773  if(dTrackSelectionTag == "NotATag")
2774  dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2775  if(dShowerSelectionTag == "NotATag")
2776  dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
2777 
2779 
2780  // Beam Particle
2781  locPID = Gamma;
2782  locParticleName = string("Beam_") + ParticleType(locPID);
2783  locParticleROOTName = ParticleName_ROOT(locPID);
2784  CreateAndChangeTo_Directory(locParticleName, locParticleName);
2785  locHistName = "Momentum";
2786  locHistTitle = string("Beam ") + locParticleROOTName + string(";p (GeV/c)");
2787  dBeamParticle_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBeamEBins, dMinP, dMaxBeamE);
2788  gDirectory->cd("..");
2789 
2790  //PID
2791  CreateAndChangeTo_Directory("PID", "PID");
2792  {
2793  //beta vs p
2794  locHistName = "BetaVsP_Q+";
2795  locHistTitle = "q^{+};p (GeV/c);#beta";
2796  dHistMap_QBetaVsP[1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2797 
2798  locHistName = "BetaVsP_Q-";
2799  locHistTitle = "q^{-};p (GeV/c);#beta";
2800  dHistMap_QBetaVsP[-1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2801  }
2802  gDirectory->cd("..");
2803 
2804  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2805  {
2806  locPID = dFinalStatePIDs[loc_i];
2807  locParticleName = ParticleType(locPID);
2808  locParticleROOTName = ParticleName_ROOT(locPID);
2809  CreateAndChangeTo_Directory(locParticleName, locParticleName);
2810 
2811  // Momentum
2812  locHistName = "Momentum";
2813  locHistTitle = locParticleROOTName + string(";p (GeV/c)");
2814  dHistMap_P[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
2815 
2816  // Theta
2817  locHistName = "Theta";
2818  locHistTitle = locParticleROOTName + string(";#theta#circ");
2819  dHistMap_Theta[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
2820 
2821  // Phi
2822  locHistName = "Phi";
2823  locHistTitle = locParticleROOTName + string(";#phi#circ");
2824  dHistMap_Phi[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
2825 
2826  // P Vs Theta
2827  locHistName = "PVsTheta";
2828  locHistTitle = locParticleROOTName + string(";#theta#circ;p (GeV/c)");
2829  dHistMap_PVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
2830 
2831  // Phi Vs Theta
2832  locHistName = "PhiVsTheta";
2833  locHistTitle = locParticleROOTName + string(";#theta#circ;#phi#circ");
2834  dHistMap_PhiVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
2835 
2836  //beta vs p
2837  locHistName = "BetaVsP";
2838  locHistTitle = locParticleROOTName + string(";p (GeV/c);#beta");
2839  dHistMap_BetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2840 
2841  //delta-beta vs p
2842  locHistName = "DeltaBetaVsP";
2843  locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#beta");
2844  dHistMap_DeltaBetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
2845 
2846  // Vertex-Z
2847  locHistName = "VertexZ";
2848  locHistTitle = locParticleROOTName + string(";Vertex-Z (cm)");
2849  dHistMap_VertexZ[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2850 
2851  // Vertex-Y Vs Vertex-X
2852  locHistName = "VertexYVsX";
2853  locHistTitle = locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
2854  dHistMap_VertexYVsX[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2855 
2856  // Vertex-T
2857  locHistName = "VertexT";
2858  locHistTitle = locParticleROOTName + string(";Vertex-T (ns)");
2859  dHistMap_VertexT[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2860 
2861  gDirectory->cd("..");
2862  }
2863 
2864  //Return to the base directory
2866  }
2867  japp->RootUnLock(); //RELEASE ROOT LOCK!!
2868 }
2869 
2870 bool DHistogramAction_DetectedParticleKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2871 {
2873  return true; //else double-counting!
2874 
2875  vector<const DBeamPhoton*> locBeamPhotons;
2876  locEventLoop->Get(locBeamPhotons);
2877 
2878  //FILL HISTOGRAMS
2879  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2880  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2881  Lock_Action();
2882  {
2883  for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
2884  dBeamParticle_P->Fill(locBeamPhotons[loc_i]->energy());
2885  }
2886  Unlock_Action();
2887 
2888  vector<const DChargedTrack*> locPreSelectChargedTracks;
2889  locEventLoop->Get(locPreSelectChargedTracks, dTrackSelectionTag.c_str());
2890 
2891  for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
2892  {
2893  const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestTrackingFOM();
2894  int locCharge = ParticleCharge(locChargedTrackHypothesis->PID());
2895 
2896  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2897  double locP = locMomentum.Mag();
2898  double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
2899 
2900  if(dHistMap_QBetaVsP.find(locCharge) == dHistMap_QBetaVsP.end())
2901  continue;
2902 
2903  //FILL HISTOGRAMS
2904  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2905  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2906  Lock_Action();
2907  {
2908  //Extremely inefficient, I know ...
2909  dHistMap_QBetaVsP[locCharge]->Fill(locP, locBeta_Timing);
2910  }
2911  Unlock_Action();
2912  }
2913 
2914  for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
2915  {
2916  const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestFOM();
2917 
2918  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2919  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2920  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2921  double locP = locMomentum.Mag();
2922  double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
2923  double locDeltaBeta = locBeta_Timing - locChargedTrackHypothesis->lorentzMomentum().Beta();
2924 
2925  Particle_t locPID = (locChargedTrackHypothesis->Get_FOM() < dMinPIDFOM) ? Unknown : locChargedTrackHypothesis->PID();
2926  if(dHistMap_P.find(locPID) == dHistMap_P.end())
2927  continue; //not interested in histogramming
2928 
2929  //FILL HISTOGRAMS
2930  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2931  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2932  Lock_Action();
2933  {
2934  dHistMap_P[locPID]->Fill(locP);
2935  dHistMap_Phi[locPID]->Fill(locPhi);
2936  dHistMap_Theta[locPID]->Fill(locTheta);
2937  dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
2938  dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
2939  dHistMap_BetaVsP[locPID]->Fill(locP, locBeta_Timing);
2940  dHistMap_DeltaBetaVsP[locPID]->Fill(locP, locDeltaBeta);
2941  dHistMap_VertexZ[locPID]->Fill(locChargedTrackHypothesis->position().Z());
2942  dHistMap_VertexYVsX[locPID]->Fill(locChargedTrackHypothesis->position().X(), locChargedTrackHypothesis->position().Y());
2943  dHistMap_VertexT[locPID]->Fill(locChargedTrackHypothesis->time());
2944  }
2945  Unlock_Action();
2946  }
2947 
2948  vector<const DNeutralParticle*> locNeutralParticles;
2949  locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
2950 
2951  for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
2952  {
2953  const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(Gamma);
2954  if(locNeutralParticleHypothesis->Get_FOM() < dMinPIDFOM)
2955  continue;
2956 
2957  Particle_t locPID = locNeutralParticleHypothesis->PID();
2958  if(dHistMap_P.find(locPID) == dHistMap_P.end())
2959  continue; //e.g. a decaying particle, or not interested in histogramming
2960 
2961  DVector3 locMomentum = locNeutralParticleHypothesis->momentum();
2962  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2963  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2964  double locP = locMomentum.Mag();
2965 
2966  double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
2967  double locDeltaBeta = locBeta_Timing - locNeutralParticleHypothesis->lorentzMomentum().Beta();
2968 
2969  //FILL HISTOGRAMS
2970  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2971  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2972  Lock_Action();
2973  {
2974  dHistMap_P[locPID]->Fill(locP);
2975  dHistMap_Phi[locPID]->Fill(locPhi);
2976  dHistMap_Theta[locPID]->Fill(locTheta);
2977  dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
2978  dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
2979  dHistMap_BetaVsP[locPID]->Fill(locP, locBeta_Timing);
2980  dHistMap_DeltaBetaVsP[locPID]->Fill(locP, locDeltaBeta);
2981  dHistMap_VertexZ[locPID]->Fill(locNeutralParticleHypothesis->position().Z());
2982  dHistMap_VertexYVsX[locPID]->Fill(locNeutralParticleHypothesis->position().X(), locNeutralParticleHypothesis->position().Y());
2983  dHistMap_VertexT[locPID]->Fill(locNeutralParticleHypothesis->time());
2984  }
2985  Unlock_Action();
2986  }
2987  return true;
2988 }
2989 
2991 {
2992  string locHistName, locHistTitle, locParticleName, locParticleROOTName;
2993  Particle_t locPID;
2994 
2995  string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
2996  if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2997  gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2998  if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
2999  gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
3000 
3001  //CREATE THE HISTOGRAMS
3002  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
3003  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
3004  {
3005  //So: Default tag is "", User can set it to something else
3006  //In here, if tag is "", get from gparms, if not, leave it alone
3007  //If gparms value does not exist, set it to (and use) "PreSelect"
3008  if(dTrackSelectionTag == "NotATag")
3009  dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
3010  if(dShowerSelectionTag == "NotATag")
3011  dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
3012 
3014 
3015  //TRACKS
3016  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3017  {
3018  locPID = dFinalStatePIDs[loc_i];
3019  locParticleName = ParticleType(locPID);
3020  locParticleROOTName = ParticleName_ROOT(locPID);
3021  CreateAndChangeTo_Directory(locParticleName, locParticleName);
3022 
3023  // Px
3024  locHistName = "PxErrorVsP";
3025  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{p_{x}} (GeV/c)");
3026  dHistMap_TrackPxErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
3027 
3028  locHistName = "PxErrorVsTheta";
3029  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{p_{x}} (GeV/c)");
3030  dHistMap_TrackPxErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
3031 
3032  locHistName = "PxErrorVsPhi";
3033  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{p_{x}} (GeV/c)");
3034  dHistMap_TrackPxErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
3035 
3036  // Py
3037  locHistName = "PyErrorVsP";
3038  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{p_{y}} (GeV/c)");
3039  dHistMap_TrackPyErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
3040 
3041  locHistName = "PyErrorVsTheta";
3042  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{p_{y}} (GeV/c)");
3043  dHistMap_TrackPyErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
3044 
3045  locHistName = "PyErrorVsPhi";
3046  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{p_{y}} (GeV/c)");
3047  dHistMap_TrackPyErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
3048 
3049  // Pz
3050  locHistName = "PzErrorVsP";
3051  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{p_{z}} (GeV/c)");
3052  dHistMap_TrackPzErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPzErrorBins, 0.0, dMaxPzError);
3053 
3054  locHistName = "PzErrorVsTheta";
3055  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{p_{z}} (GeV/c)");
3056  dHistMap_TrackPzErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPzErrorBins, 0.0, dMaxPzError);
3057 
3058  locHistName = "PzErrorVsPhi";
3059  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{p_{z}} (GeV/c)");
3060  dHistMap_TrackPzErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DPzErrorBins, 0.0, dMaxPzError);
3061 
3062  // X
3063  locHistName = "XErrorVsP";
3064  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{x} (cm)");
3065  dHistMap_TrackXErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
3066 
3067  locHistName = "XErrorVsTheta";
3068  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{x} (cm)");
3069  dHistMap_TrackXErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
3070 
3071  locHistName = "XErrorVsPhi";
3072  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{x} (cm)");
3073  dHistMap_TrackXErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
3074 
3075  // Y
3076  locHistName = "YErrorVsP";
3077  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{y} (cm)");
3078  dHistMap_TrackYErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
3079 
3080  locHistName = "YErrorVsTheta";
3081  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{y} (cm)");
3082  dHistMap_TrackYErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
3083 
3084  locHistName = "YErrorVsPhi";
3085  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{y} (cm)");
3086  dHistMap_TrackYErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
3087 
3088  // Z
3089  locHistName = "ZErrorVsP";
3090  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{z} (cm)");
3091  dHistMap_TrackZErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DZErrorBins, 0.0, dMaxZError);
3092 
3093  locHistName = "ZErrorVsTheta";
3094  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{z} (cm)");
3095  dHistMap_TrackZErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DZErrorBins, 0.0, dMaxZError);
3096 
3097  locHistName = "ZErrorVsPhi";
3098  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{z} (cm)");
3099  dHistMap_TrackZErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DZErrorBins, 0.0, dMaxZError);
3100 
3101  gDirectory->cd("..");
3102  }
3103 
3104  //SHOWERS
3105  for(bool locIsBCALFlag : {false, true})
3106  {
3107  string locDirName = locIsBCALFlag ? "Photon_BCAL" : "Photon_FCAL";
3108  locParticleROOTName = ParticleName_ROOT(Gamma);
3109  CreateAndChangeTo_Directory(locDirName, locDirName);
3110 
3111  double locMaxP = locIsBCALFlag ? dMaxPBCAL : dMaxP;
3112  double locMinTheta = locIsBCALFlag ? dMinThetaBCAL : dMinTheta;
3113  double locMaxTheta = locIsBCALFlag ? dMaxTheta : dMaxThetaFCAL;
3114 
3115  // E
3116  locHistName = "EErrorVsP";
3117  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{E} (GeV)");
3118  dHistMap_ShowerEErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DEErrorBins, 0.0, dMaxEError);
3119 
3120  locHistName = "EErrorVsTheta";
3121  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{E} (GeV)");
3122  dHistMap_ShowerEErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DEErrorBins, 0.0, dMaxEError);
3123 
3124  locHistName = "EErrorVsPhi";
3125  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{E} (GeV)");
3126  dHistMap_ShowerEErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DEErrorBins, 0.0, dMaxEError);
3127 
3128  // X
3129  locHistName = "XErrorVsP";
3130  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{x} (cm)");
3131  dHistMap_ShowerXErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
3132 
3133  locHistName = "XErrorVsTheta";
3134  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{x} (cm)");
3135  dHistMap_ShowerXErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
3136 
3137  locHistName = "XErrorVsPhi";
3138  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{x} (cm)");
3139  dHistMap_ShowerXErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
3140 
3141  // Y
3142  locHistName = "YErrorVsP";
3143  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{y} (cm)");
3144  dHistMap_ShowerYErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
3145 
3146  locHistName = "YErrorVsTheta";
3147  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{y} (cm)");
3148  dHistMap_ShowerYErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
3149 
3150  locHistName = "YErrorVsPhi";
3151  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{y} (cm)");
3152  dHistMap_ShowerYErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
3153 
3154  // Z
3155  locHistName = "ZErrorVsP";
3156  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{z} (cm)");
3157  dHistMap_ShowerZErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DShowerZErrorBins, 0.0, dMaxShowerZError);
3158 
3159  locHistName = "ZErrorVsTheta";
3160  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{z} (cm)");
3161  dHistMap_ShowerZErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DShowerZErrorBins, 0.0, dMaxShowerZError);
3162 
3163  locHistName = "ZErrorVsPhi";
3164  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{z} (cm)");
3165  dHistMap_ShowerZErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DShowerZErrorBins, 0.0, dMaxShowerZError);
3166 
3167  // T
3168  locHistName = "TErrorVsP";
3169  locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{t} (ns)");
3170  dHistMap_ShowerTErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DTErrorBins, 0.0, dMaxTError);
3171 
3172  locHistName = "TErrorVsTheta";
3173  locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{t} (ns)");
3174  dHistMap_ShowerTErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DTErrorBins, 0.0, dMaxTError);
3175 
3176  locHistName = "TErrorVsPhi";
3177  locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{t} (ns)");
3178  dHistMap_ShowerTErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DTErrorBins, 0.0, dMaxTError);
3179 
3180  gDirectory->cd("..");
3181  }
3182 
3183  //Return to the base directory
3185  }
3186  japp->RootUnLock(); //RELEASE ROOT LOCK!!
3187 }
3188 
3189 bool DHistogramAction_TrackShowerErrors::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3190 {
3192  return true; //else double-counting!
3193 
3194  vector<const DChargedTrack*> locPreSelectChargedTracks;
3195  locEventLoop->Get(locPreSelectChargedTracks, dTrackSelectionTag.c_str());
3196 
3197  for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
3198  {
3199  const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestFOM();
3200 
3201  Particle_t locPID = (locChargedTrackHypothesis->Get_FOM() < dMinPIDFOM) ? Unknown : locChargedTrackHypothesis->PID();
3202  if(dHistMap_TrackPxErrorVsP.find(locPID) == dHistMap_TrackPxErrorVsP.end())
3203  continue; //not interested in histogramming
3204 
3205  DVector3 locMomentum = locChargedTrackHypothesis->momentum();
3206  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
3207  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
3208  double locP = locMomentum.Mag();
3209 
3210  const TMatrixFSym& locCovarianceMatrix = *(locChargedTrackHypothesis->errorMatrix().get());
3211  double locPxError = sqrt(locCovarianceMatrix(0, 0));
3212  double locPyError = sqrt(locCovarianceMatrix(1, 1));
3213  double locPzError = sqrt(locCovarianceMatrix(2, 2));
3214  double locXError = sqrt(locCovarianceMatrix(3, 3));
3215  double locYError = sqrt(locCovarianceMatrix(4, 4));
3216  double locZError = sqrt(locCovarianceMatrix(5, 5));
3217 
3218  //FILL HISTOGRAMS
3219  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3220  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3221  Lock_Action();
3222  {
3223  dHistMap_TrackPxErrorVsP[locPID]->Fill(locP, locPxError);
3224  dHistMap_TrackPxErrorVsTheta[locPID]->Fill(locTheta, locPxError);
3225  dHistMap_TrackPxErrorVsPhi[locPID]->Fill(locPhi, locPxError);
3226 
3227  dHistMap_TrackPyErrorVsP[locPID]->Fill(locP, locPyError);
3228  dHistMap_TrackPyErrorVsTheta[locPID]->Fill(locTheta, locPyError);
3229  dHistMap_TrackPyErrorVsPhi[locPID]->Fill(locPhi, locPyError);
3230 
3231  dHistMap_TrackPzErrorVsP[locPID]->Fill(locP, locPzError);
3232  dHistMap_TrackPzErrorVsTheta[locPID]->Fill(locTheta, locPzError);
3233  dHistMap_TrackPzErrorVsPhi[locPID]->Fill(locPhi, locPzError);
3234 
3235  dHistMap_TrackXErrorVsP[locPID]->Fill(locP, locXError);
3236  dHistMap_TrackXErrorVsTheta[locPID]->Fill(locTheta, locXError);
3237  dHistMap_TrackXErrorVsPhi[locPID]->Fill(locPhi, locXError);
3238 
3239  dHistMap_TrackYErrorVsP[locPID]->Fill(locP, locYError);
3240  dHistMap_TrackYErrorVsTheta[locPID]->Fill(locTheta, locYError);
3241  dHistMap_TrackYErrorVsPhi[locPID]->Fill(locPhi, locYError);
3242 
3243  dHistMap_TrackZErrorVsP[locPID]->Fill(locP, locZError);
3244  dHistMap_TrackZErrorVsTheta[locPID]->Fill(locTheta, locZError);
3245  dHistMap_TrackZErrorVsPhi[locPID]->Fill(locPhi, locZError);
3246  }
3247  Unlock_Action();
3248  }
3249 
3250  vector<const DNeutralParticle*> locNeutralParticles;
3251  locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
3252 
3253  for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
3254  {
3255  const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(Gamma);
3256  if(locNeutralParticleHypothesis->Get_FOM() < dMinPIDFOM)
3257  continue;
3258 
3259  DVector3 locMomentum = locNeutralParticleHypothesis->momentum();
3260  double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
3261  double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
3262  double locP = locMomentum.Mag();
3263 
3264  const DNeutralShower* locNeutralShower = locNeutralParticles[loc_i]->dNeutralShower;
3265  const TMatrixFSym& locCovarianceMatrix = *(locNeutralShower->dCovarianceMatrix);
3266  bool locIsBCALFlag = (locNeutralShower->dDetectorSystem == SYS_BCAL);
3267  double locEError = sqrt(locCovarianceMatrix(0, 0));
3268  double locXError = sqrt(locCovarianceMatrix(1, 1));
3269  double locYError = sqrt(locCovarianceMatrix(2, 2));
3270  double locZError = sqrt(locCovarianceMatrix(3, 3));
3271  double locTError = sqrt(locCovarianceMatrix(4, 4));
3272 
3273  //FILL HISTOGRAMS
3274  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3275  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3276  Lock_Action();
3277  {
3278  dHistMap_ShowerEErrorVsP[locIsBCALFlag]->Fill(locP, locEError);
3279  dHistMap_ShowerEErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locEError);
3280  dHistMap_ShowerEErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locEError);
3281 
3282  dHistMap_ShowerXErrorVsP[locIsBCALFlag]->Fill(locP, locXError);
3283  dHistMap_ShowerXErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locXError);
3284  dHistMap_ShowerXErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locXError);
3285 
3286  dHistMap_ShowerYErrorVsP[locIsBCALFlag]->Fill(locP, locYError);
3287  dHistMap_ShowerYErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locYError);
3288  dHistMap_ShowerYErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locYError);
3289 
3290  dHistMap_ShowerZErrorVsP[locIsBCALFlag]->Fill(locP, locZError);
3291  dHistMap_ShowerZErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locZError);
3292  dHistMap_ShowerZErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locZError);
3293 
3294  dHistMap_ShowerTErrorVsP[locIsBCALFlag]->Fill(locP, locTError);
3295  dHistMap_ShowerTErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locTError);
3296  dHistMap_ShowerTErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locTError);
3297  }
3298  Unlock_Action();
3299  }
3300 
3301  return true;
3302 }
3303 
3305 {
3306  string locHistName;
3307 
3308  bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
3309 
3310  //CREATE THE HISTOGRAMS
3311  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
3312  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
3313  {
3315 
3316  //2D Summary
3317  locHistName = "NumHighLevelObjects";
3318  if(gDirectory->Get(locHistName.c_str()) != NULL) //already created by another thread, or directory name is duplicate (e.g. two identical steps)
3319  dHist_NumHighLevelObjects = static_cast<TH2D*>(gDirectory->Get(locHistName.c_str()));
3320  else
3321  {
3322  dHist_NumHighLevelObjects = new TH2D(locHistName.c_str(), ";;# Objects / Event", 13, 0.5, 13.5, dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3323  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(1, "DRFTime");
3324  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(2, "DSCHit");
3325  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(3, "DTOFPoint");
3326  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(4, "DBCALShower");
3327  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(5, "DFCALShower");
3328  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(6, "DTimeBasedTrack");
3329  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(7, "TrackSCMatches");
3330  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(8, "TrackTOFMatches");
3331  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(9, "TrackBCALMatches");
3332  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(10, "TrackFCALMatches");
3333  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(11, "DBeamPhoton");
3334  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(12, "DChargedTrack");
3335  dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(13, "DNeutralShower");
3336  }
3337 
3338  //Charged
3339  locHistName = "NumChargedTracks";
3340  dHist_NumChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3341  locHistName = "NumPosChargedTracks";
3342  dHist_NumPosChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3343  locHistName = "NumNegChargedTracks";
3344  dHist_NumNegChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3345 
3346  //TBT
3347  locHistName = "NumTimeBasedTracks";
3348  dHist_NumTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3349  locHistName = "NumPosTimeBasedTracks";
3350  dHist_NumPosTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3351  locHistName = "NumNegTimeBasedTracks";
3352  dHist_NumNegTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3353 
3354  if(!locIsRESTEvent)
3355  {
3356  //WBT
3357  locHistName = "NumWireBasedTracks";
3358  dHist_NumWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3359  locHistName = "NumPosWireBasedTracks";
3360  dHist_NumPosWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3361  locHistName = "NumNegWireBasedTracks";
3362  dHist_NumNegWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3363 
3364  //Track Candidates
3365  locHistName = "NumTrackCandidates";
3366  dHist_NumTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3367  locHistName = "NumPosTrackCandidates";
3368  dHist_NumPosTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3369  locHistName = "NumNegTrackCandidates";
3370  dHist_NumNegTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3371 
3372  //CDC Track Candidates
3373  locHistName = "NumPosTrackCandidates_CDC";
3374  dHist_NumPosTrackCandidates_CDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} CDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3375  locHistName = "NumNegTrackCandidates_CDC";
3376  dHist_NumNegTrackCandidates_CDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} CDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3377 
3378  //FDC Track Candidates
3379  locHistName = "NumPosTrackCandidates_FDC";
3380  dHist_NumPosTrackCandidates_FDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} FDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3381  locHistName = "NumNegTrackCandidates_FDC";
3382  dHist_NumNegTrackCandidates_FDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} FDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3383  }
3384 
3385  //Beam Photons
3386  locHistName = "NumBeamPhotons";
3387  dHist_NumBeamPhotons = GetOrCreate_Histogram<TH1D>(locHistName, ";# DBeamPhoton", dMaxNumBeamPhotons + 1, -0.5, (float)dMaxNumBeamPhotons + 0.5);
3388 
3389  //Showers / Neutrals / TOF / SC
3390  locHistName = "NumFCALShowers";
3391  dHist_NumFCALShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DFCALShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3392  locHistName = "NumBCALShowers";
3393  dHist_NumBCALShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DBCALShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3394  locHistName = "NumNeutralShowers";
3395  dHist_NumNeutralShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DNeutralShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3396  locHistName = "NumTOFPoints";
3397  dHist_NumTOFPoints = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTOFPoint", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3398  locHistName = "NumSCHits";
3399  dHist_NumSCHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DSCHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3400 
3401  if(!locIsRESTEvent)
3402  {
3403  locHistName = "NumTAGMHits";
3404  dHist_NumTAGMHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTAGMHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3405  locHistName = "NumTAGHHits";
3406  dHist_NumTAGHHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTAGHHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3407  }
3408 
3409  //Matches
3410  locHistName = "NumTrackBCALMatches";
3411  dHist_NumTrackBCALMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-BCAL Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3412  locHistName = "NumTrackFCALMatches";
3413  dHist_NumTrackFCALMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-FCAL Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3414  locHistName = "NumTrackTOFMatches";
3415  dHist_NumTrackTOFMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-TOF Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3416  locHistName = "NumTrackSCMatches";
3417  dHist_NumTrackSCMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-SC Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3418 
3419  if(!locIsRESTEvent)
3420  {
3421  //Hits
3422  locHistName = "NumCDCHits";
3423  dHist_NumCDCHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DCDCHit", dMaxNumCDCHits + 1, -0.5, (float)dMaxNumCDCHits + 0.5);
3424  locHistName = "NumFDCWireHits";
3425  dHist_NumFDCWireHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# Wire DFDCHit", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
3426  locHistName = "NumFDCCathodeHits";
3427  dHist_NumFDCCathodeHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# Cathode DFDCHit", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
3428  locHistName = "NumFDCPseudoHits";
3429  dHist_NumFDCPseudoHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# FDC Pseudo Hits", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
3430 
3431  locHistName = "NumTOFHits";
3432  dHist_NumTOFHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DTOFHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
3433  locHistName = "NumBCALHits";
3434  dHist_NumBCALHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DBCALHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
3435  locHistName = "NumFCALHits";
3436  dHist_NumFCALHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DFCALHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
3437 
3438  locHistName = "NumRFSignals";
3439  dHist_NumRFSignals = GetOrCreate_Histogram<TH1I>(locHistName, ";# DRFDigiTime + # DRFTDCDigiTime", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3440  }
3441 
3442  //Return to the base directory
3444  }
3445  japp->RootUnLock(); //RELEASE ROOT LOCK!!
3446 }
3447 
3448 bool DHistogramAction_NumReconstructedObjects::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3449 {
3451  return true; //else double-counting!
3452 
3453  bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
3454 
3455  vector<const DTrackTimeBased*> locTrackTimeBasedVector;
3456  locEventLoop->Get(locTrackTimeBasedVector);
3457 
3458  vector<const DBeamPhoton*> locBeamPhotons;
3459  locEventLoop->Get(locBeamPhotons);
3460 
3461  vector<const DFCALShower*> locFCALShowers;
3462  locEventLoop->Get(locFCALShowers);
3463 
3464  vector<const DChargedTrack*> locChargedTracks;
3465  locEventLoop->Get(locChargedTracks);
3466 
3467  vector<const DBCALShower*> locBCALShowers;
3468  locEventLoop->Get(locBCALShowers);
3469 
3470  vector<const DNeutralShower*> locNeutralShowers;
3471  locEventLoop->Get(locNeutralShowers);
3472 
3473  vector<const DTOFPoint*> locTOFPoints;
3474  locEventLoop->Get(locTOFPoints);
3475 
3476  vector<const DSCHit*> locSCHits;
3477  locEventLoop->Get(locSCHits);
3478 
3479  vector<const DRFTime*> locRFTimes;
3480  locEventLoop->Get(locRFTimes);
3481 
3482  const DDetectorMatches* locDetectorMatches = NULL;
3483  locEventLoop->GetSingle(locDetectorMatches);
3484 
3485  //if not REST
3486  vector<const DTrackWireBased*> locTrackWireBasedVector;
3487  vector<const DTrackCandidate*> locTrackCandidates;
3488  vector<const DTrackCandidate*> locTrackCandidates_CDC;
3489  vector<const DTrackCandidate*> locTrackCandidates_FDC;
3490  vector<const DCDCHit*> locCDCHits;
3491  vector<const DFDCHit*> locFDCHits;
3492  vector<const DTOFHit*> locTOFHits;
3493  vector<const DBCALHit*> locBCALHits;
3494  vector<const DFCALHit*> locFCALHits;
3495  vector<const DTAGMHit*> locTAGMHits;
3496  vector<const DTAGHHit*> locTAGHHits;
3497  vector<const DFDCPseudo*> locFDCPseudoHits;
3498  vector<const DRFDigiTime*> locRFDigiTimes;
3499  vector<const DRFTDCDigiTime*> locRFTDCDigiTimes;
3500 
3501  size_t locNumFDCWireHits = 0, locNumFDCCathodeHits = 0;
3502  if(!locIsRESTEvent)
3503  {
3504  locEventLoop->Get(locTrackWireBasedVector);
3505  locEventLoop->Get(locTrackCandidates);
3506  locEventLoop->Get(locTrackCandidates_CDC, "CDC");
3507  locEventLoop->Get(locTrackCandidates_FDC, "FDCCathodes");
3508  locEventLoop->Get(locCDCHits);
3509  locEventLoop->Get(locFDCHits);
3510  locEventLoop->Get(locFDCPseudoHits);
3511  locEventLoop->Get(locTOFHits);
3512  locEventLoop->Get(locBCALHits);
3513  locEventLoop->Get(locFCALHits);
3514  locEventLoop->Get(locTAGHHits);
3515  locEventLoop->Get(locTAGMHits);
3516  locEventLoop->Get(locRFDigiTimes);
3517  locEventLoop->Get(locRFTDCDigiTimes);
3518 
3519  for(size_t loc_i = 0; loc_i < locFDCHits.size(); ++loc_i)
3520  {
3521  if(locFDCHits[loc_i]->type == DFDCHit::AnodeWire)
3522  ++locNumFDCWireHits;
3523  else
3524  ++locNumFDCCathodeHits;
3525  }
3526  }
3527 
3528  //FILL HISTOGRAMS
3529  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3530  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3531  Lock_Action(); //ACQUIRE ROOT LOCK!!
3532  {
3533  //# High-Level Objects
3534  dHist_NumHighLevelObjects->Fill(1, (Double_t)locRFTimes.size());
3535  dHist_NumHighLevelObjects->Fill(2, (Double_t)locSCHits.size());
3536  dHist_NumHighLevelObjects->Fill(3, (Double_t)locTOFPoints.size());
3537  dHist_NumHighLevelObjects->Fill(4, (Double_t)locBCALShowers.size());
3538  dHist_NumHighLevelObjects->Fill(5, (Double_t)locFCALShowers.size());
3539  dHist_NumHighLevelObjects->Fill(6, (Double_t)locTrackTimeBasedVector.size());
3540  dHist_NumHighLevelObjects->Fill(7, (Double_t)locDetectorMatches->Get_NumTrackSCMatches());
3541  dHist_NumHighLevelObjects->Fill(8, (Double_t)locDetectorMatches->Get_NumTrackTOFMatches());
3542  dHist_NumHighLevelObjects->Fill(9, (Double_t)locDetectorMatches->Get_NumTrackBCALMatches());
3543  dHist_NumHighLevelObjects->Fill(10, (Double_t)locDetectorMatches->Get_NumTrackFCALMatches());
3544  dHist_NumHighLevelObjects->Fill(11, (Double_t)locBeamPhotons.size());
3545  dHist_NumHighLevelObjects->Fill(12, (Double_t)locChargedTracks.size());
3546  dHist_NumHighLevelObjects->Fill(13, (Double_t)locNeutralShowers.size());
3547 
3548  //Charged
3549  unsigned int locNumPos = 0, locNumNeg = 0;
3550  for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
3551  {
3552  if(ParticleCharge(locChargedTracks[loc_i]->Get_BestFOM()->PID()) > 0)
3553  ++locNumPos;
3554  else
3555  ++locNumNeg;
3556  }
3557  dHist_NumChargedTracks->Fill(locChargedTracks.size());
3558  dHist_NumPosChargedTracks->Fill(locNumPos);
3559  dHist_NumNegChargedTracks->Fill(locNumNeg);
3560 
3561  //TBT
3562  locNumPos = 0; locNumNeg = 0;
3563  for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
3564  {
3565  if(ParticleCharge(locTrackTimeBasedVector[loc_i]->PID()) > 0)
3566  ++locNumPos;
3567  else
3568  ++locNumNeg;
3569  }
3570  dHist_NumTimeBasedTracks->Fill(locTrackTimeBasedVector.size());
3571  dHist_NumPosTimeBasedTracks->Fill(locNumPos);
3572  dHist_NumNegTimeBasedTracks->Fill(locNumNeg);
3573 
3574  if(!locIsRESTEvent)
3575  {
3576  //WBT
3577  locNumPos = 0; locNumNeg = 0;
3578  for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
3579  {
3580  if(ParticleCharge(locTrackWireBasedVector[loc_i]->PID()) > 0)
3581  ++locNumPos;
3582  else
3583  ++locNumNeg;
3584  }
3585  dHist_NumWireBasedTracks->Fill(locTrackWireBasedVector.size());
3586  dHist_NumPosWireBasedTracks->Fill(locNumPos);
3587  dHist_NumNegWireBasedTracks->Fill(locNumNeg);
3588 
3589  //Candidates
3590  locNumPos = 0; locNumNeg = 0;
3591  for(size_t loc_i = 0; loc_i < locTrackCandidates.size(); ++loc_i)
3592  {
3593  if(locTrackCandidates[loc_i]->charge() > 0.0)
3594  ++locNumPos;
3595  else
3596  ++locNumNeg;
3597  }
3598  dHist_NumTrackCandidates->Fill(locTrackCandidates.size());
3599  dHist_NumPosTrackCandidates->Fill(locNumPos);
3600  dHist_NumNegTrackCandidates->Fill(locNumNeg);
3601 
3602  //CDC Candidates
3603  locNumPos = 0; locNumNeg = 0;
3604  for(size_t loc_i = 0; loc_i < locTrackCandidates_CDC.size(); ++loc_i)
3605  {
3606  if(locTrackCandidates_CDC[loc_i]->charge() > 0.0)
3607  ++locNumPos;
3608  else
3609  ++locNumNeg;
3610  }
3611  dHist_NumPosTrackCandidates_CDC->Fill(locNumPos);
3612  dHist_NumNegTrackCandidates_CDC->Fill(locNumNeg);
3613 
3614  //FDC Candidates
3615  locNumPos = 0; locNumNeg = 0;
3616  for(size_t loc_i = 0; loc_i < locTrackCandidates_FDC.size(); ++loc_i)
3617  {
3618  if(locTrackCandidates_FDC[loc_i]->charge() > 0.0)
3619  ++locNumPos;
3620  else
3621  ++locNumNeg;
3622  }
3623  dHist_NumPosTrackCandidates_FDC->Fill(locNumPos);
3624  dHist_NumNegTrackCandidates_FDC->Fill(locNumNeg);
3625  }
3626 
3627  //Beam Photons
3628  dHist_NumBeamPhotons->Fill((Double_t)locBeamPhotons.size());
3629 
3630  //Showers
3631  dHist_NumFCALShowers->Fill((Double_t)locFCALShowers.size());
3632  dHist_NumBCALShowers->Fill((Double_t)locBCALShowers.size());
3633  dHist_NumNeutralShowers->Fill((Double_t)locNeutralShowers.size());
3634 
3635  //TOF & SC
3636  dHist_NumTOFPoints->Fill((Double_t)locTOFPoints.size());
3637  dHist_NumSCHits->Fill((Double_t)locSCHits.size());
3638 
3639  //TAGGER
3640  if(!locIsRESTEvent)
3641  {
3642  dHist_NumTAGMHits->Fill((Double_t)locTAGMHits.size());
3643  dHist_NumTAGHHits->Fill((Double_t)locTAGHHits.size());
3644  }
3645 
3646  //Matches
3647  dHist_NumTrackBCALMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackBCALMatches());
3648  dHist_NumTrackFCALMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackFCALMatches());
3649  dHist_NumTrackTOFMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackTOFMatches());
3650  dHist_NumTrackSCMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackSCMatches());
3651 
3652  //Hits
3653  if(!locIsRESTEvent)
3654  {
3655  dHist_NumCDCHits->Fill((Double_t)locCDCHits.size());
3656  dHist_NumFDCWireHits->Fill((Double_t)locNumFDCWireHits);
3657  dHist_NumFDCCathodeHits->Fill((Double_t)locNumFDCCathodeHits);
3658  dHist_NumFDCPseudoHits->Fill((Double_t)locFDCPseudoHits.size());
3659  dHist_NumTOFHits->Fill((Double_t)locTOFHits.size());
3660  dHist_NumBCALHits->Fill((Double_t)locBCALHits.size());
3661  dHist_NumFCALHits->Fill((Double_t)locFCALHits.size());
3662  dHist_NumRFSignals->Fill((Double_t)(locRFDigiTimes.size() + locRFTDCDigiTimes.size()));
3663  }
3664  }
3665  Unlock_Action(); //RELEASE ROOT LOCK!!
3666 
3667  return true;
3668 }
3669 
3671 {
3672  string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
3673  if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
3674  gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
3675  if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
3676  gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
3677 
3678  //CREATE THE HISTOGRAMS
3679  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
3680  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
3681  {
3682  //So: Default tag is "", User can set it to something else
3683  //In here, if tag is "", get from gparms, if not, leave it alone
3684  //If gparms value does not exist, set it to (and use) "PreSelect"
3685  if(dTrackSelectionTag == "NotATag")
3686  dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
3687  if(dShowerSelectionTag == "NotATag")
3688  dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
3689 
3691 
3692  string locHistName("NumReconstructedParticles");
3693  if(gDirectory->Get(locHistName.c_str()) != NULL) //already created by another thread, or directory name is duplicate (e.g. two identical steps)
3694  dHist_NumReconstructedParticles = static_cast<TH2D*>(gDirectory->Get(locHistName.c_str()));
3695  else
3696  {
3697  dHist_NumReconstructedParticles = new TH2D("NumReconstructedParticles", ";Particle Type;Num Particles / Event", 5 + dFinalStatePIDs.size(), -0.5, 4.5 + dFinalStatePIDs.size(), dMaxNumTracks + 1, -0.5, (float)dMaxNumTracks + 0.5);
3698  dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(1, "# Total");
3699  dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(2, "# q != 0");
3700  dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(3, "# q = 0");
3701  dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(4, "# q = +");
3702  dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(5, "# q = -");
3703  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3704  {
3705  string locLabelName = string("# ") + string(ParticleName_ROOT(dFinalStatePIDs[loc_i]));
3706  dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(6 + loc_i, locLabelName.c_str());
3707  }
3708  }
3709 
3710  locHistName = "NumGoodReconstructedParticles";
3711  if(gDirectory->Get(locHistName.c_str()) != NULL) //already created by another thread, or directory name is duplicate (e.g. two identical steps)
3712  dHist_NumGoodReconstructedParticles = static_cast<TH2D*>(gDirectory->Get(locHistName.c_str()));
3713  else
3714  {
3715  dHist_NumGoodReconstructedParticles = new TH2D("NumGoodReconstructedParticles", ";Particle Type;Num Particles / Event", 5 + dFinalStatePIDs.size(), -0.5, 4.5 + dFinalStatePIDs.size(), dMaxNumTracks + 1, -0.5, (float)dMaxNumTracks + 0.5);
3716  dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(1, "# Total");
3717  dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(2, "# q != 0");
3718  dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(3, "# q = 0");
3719  dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(4, "# q = +");
3720  dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(5, "# q = -");
3721  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3722  {
3723  string locLabelName = string("# ") + string(ParticleName_ROOT(dFinalStatePIDs[loc_i]));
3724  dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(6 + loc_i, locLabelName.c_str());
3725  }
3726  }
3727 
3728  //Return to the base directory
3730  }
3731  japp->RootUnLock(); //RELEASE ROOT LOCK!!
3732 }
3733 
3734 bool DHistogramAction_TrackMultiplicity::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3735 {
3737  return true; //else double-counting!
3738 
3739  vector<const DChargedTrack*> locChargedTracks;
3740  locEventLoop->Get(locChargedTracks);
3741 
3742  vector<const DChargedTrack*> locGoodChargedTracks;
3743  locEventLoop->Get(locGoodChargedTracks, dTrackSelectionTag.c_str());
3744 
3745  // get #tracks by PID/q type
3746  size_t locNumPositiveTracks = 0, locNumNegativeTracks = 0;
3747  map<Particle_t, size_t> locNumTracksByPID;
3748  for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
3749  {
3750  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
3751  Particle_t locPID = locChargedTrackHypothesis->PID();
3752 
3753  if(locChargedTrackHypothesis->charge() > 0.0)
3754  ++locNumPositiveTracks;
3755  else
3756  ++locNumNegativeTracks;
3757 
3758  if(locNumTracksByPID.find(locPID) != locNumTracksByPID.end())
3759  ++locNumTracksByPID[locPID];
3760  else
3761  locNumTracksByPID[locPID] = 1;
3762  }
3763 
3764  // get # good tracks by PID/q type
3765  size_t locNumGoodPositiveTracks = 0, locNumGoodNegativeTracks = 0;
3766  map<Particle_t, size_t> locNumGoodTracksByPID;
3767  for(size_t loc_i = 0; loc_i < locGoodChargedTracks.size(); ++loc_i)
3768  {
3769  const DChargedTrackHypothesis* locChargedTrackHypothesis = locGoodChargedTracks[loc_i]->Get_BestFOM();
3770  Particle_t locPID = locChargedTrackHypothesis->PID();
3771 
3772  double locPIDFOM = locChargedTrackHypothesis->Get_FOM();
3773 
3774  if(locChargedTrackHypothesis->charge() > 0.0)
3775  ++locNumGoodPositiveTracks;
3776  else
3777  ++locNumGoodNegativeTracks;
3778 
3779  if(locPIDFOM < dMinPIDFOM)
3780  continue;
3781 
3782  if(locNumGoodTracksByPID.find(locPID) != locNumGoodTracksByPID.end())
3783  ++locNumGoodTracksByPID[locPID];
3784  else
3785  locNumGoodTracksByPID[locPID] = 1;
3786  }
3787 
3788  vector<const DNeutralParticle*> locNeutralParticles;
3789  locEventLoop->Get(locNeutralParticles);
3790 
3791  vector<const DNeutralParticle*> locGoodNeutralParticles;
3792  locEventLoop->Get(locGoodNeutralParticles, dShowerSelectionTag.c_str());
3793 
3794  // neutrals by pid
3795  for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
3796  {
3797  const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
3798  if(locNeutralParticleHypothesis->Get_FOM() < dMinPIDFOM)
3799  continue;
3800 
3801  Particle_t locPID = locNeutralParticleHypothesis->PID();
3802  if(locNumTracksByPID.find(locPID) != locNumTracksByPID.end())
3803  ++locNumTracksByPID[locPID];
3804  else
3805  locNumTracksByPID[locPID] = 1;
3806  }
3807 
3808  // good neutrals
3809  for(size_t loc_i = 0; loc_i < locGoodNeutralParticles.size(); ++loc_i)
3810  {
3811  const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locGoodNeutralParticles[loc_i]->Get_BestFOM();
3812  if(locNeutralParticleHypothesis->Get_FOM() < dMinPIDFOM)
3813  continue;
3814 
3815  Particle_t locPID = locNeutralParticleHypothesis->PID();
3816  if(locNumGoodTracksByPID.find(locPID) != locNumGoodTracksByPID.end())
3817  ++locNumGoodTracksByPID[locPID];
3818  else
3819  locNumGoodTracksByPID[locPID] = 1;
3820  }
3821 
3822  size_t locNumGoodTracks = locNumGoodPositiveTracks + locNumGoodNegativeTracks;
3823 
3824  //FILL HISTOGRAMS
3825  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3826  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3827  Lock_Action();
3828  {
3829  dHist_NumReconstructedParticles->Fill(0.0, (Double_t)(locChargedTracks.size() + locNeutralParticles.size()));
3830  dHist_NumReconstructedParticles->Fill(1.0, (Double_t)locChargedTracks.size());
3831  dHist_NumReconstructedParticles->Fill(2.0, (Double_t)locNeutralParticles.size());
3832  dHist_NumReconstructedParticles->Fill(3.0, (Double_t)locNumPositiveTracks);
3833  dHist_NumReconstructedParticles->Fill(4.0, (Double_t)locNumNegativeTracks);
3834  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3835  dHist_NumReconstructedParticles->Fill(5.0 + (Double_t)loc_i, (Double_t)locNumTracksByPID[dFinalStatePIDs[loc_i]]);
3836 
3837  dHist_NumGoodReconstructedParticles->Fill(0.0, (Double_t)(locNumGoodTracks + locGoodNeutralParticles.size()));
3838  dHist_NumGoodReconstructedParticles->Fill(1.0, (Double_t)locNumGoodTracks);
3839  dHist_NumGoodReconstructedParticles->Fill(2.0, (Double_t)locGoodNeutralParticles.size());
3840  dHist_NumGoodReconstructedParticles->Fill(3.0, (Double_t)locNumGoodPositiveTracks);
3841  dHist_NumGoodReconstructedParticles->Fill(4.0, (Double_t)locNumGoodNegativeTracks);
3842  for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3843  dHist_NumGoodReconstructedParticles->Fill(5.0 + (Double_t)loc_i, (Double_t)locNumGoodTracksByPID[dFinalStatePIDs[loc_i]]);
3844  }
3845  Unlock_Action();
3846 
3847  return true;
3848 }
3849 
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
map< DetectorSystem_t, map< int, TH2I * > > dHistMap_BetaVsP
TDirectoryFile * ChangeTo_BaseDirectory(void)
double getEnergy() const
Definition: DFCALShower.h:156
unsigned int PredictSCSector(const DReferenceTrajectory *rt, DVector3 *locOutputProjPos=nullptr, bool *locProjBarrelRegion=nullptr, double *locMinDPhi=nullptr) const
void Get_FDCPlanes(unsigned int locBitPattern, set< int > &locFDCPlanes) const
void Get_CDCRings(unsigned int locBitPattern, set< int > &locCDCRings) const
int dHorizontalBar
Definition: DTOFPoint.h:38
map< Particle_t, TH2I * > dHistMap_TrackPyErrorVsTheta
map< bool, TH2I * > dHistMap_TOFPointTrackDeltaYVsVerticalPaddle
DRFTime_factory * locRFTimeFactory
static char * ParticleName_ROOT(Particle_t p)
Definition: particleType.h:851
const DTOFPaddleHit * Get_ClosestTOFPaddleHit_Horizontal(const DReferenceTrajectory *locReferenceTrajectory, const vector< const DTOFPaddleHit * > &locTOFPaddleHits, double locInputStartTime, double &locBestDeltaY, double &locBestDistance) const
DResourcePool< DParticleComboStep > dResourcePool_ParticleComboStep
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
bool Get_DistanceToNearestTrack(const DBCALShower *locBCALShower, double &locDistance) const
double z(void) const
map< pair< int, bool >, TH2I * > dHistMap_SCEnergyVsTheta
DResourcePool< DParticleCombo > dResourcePool_ParticleCombo
TVector3 DVector3
Definition: DVector3.h:14
Definition: GlueX.h:16
unsigned int dKinFitNDF
Definition: DVertex.h:31
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
map< pair< int, bool >, TH2I * > dHistMap_BCALShowerTrackDepthVsP
char string[256]
map< bool, TH2I * > dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit
axes Fill(100, 100)
double getTime() const
Definition: DFCALShower.h:161
map< Particle_t, TH1I * > dHistMap_NumPhotons_DIRC
TDirectoryFile * CreateAndChangeTo_Directory(TDirectoryFile *locBaseDirectory, string locDirName, string locDirTitle)
#define y
map< DetectorSystem_t, map< bool, TH2I * > > dHistMap_PVsTheta_HasHit
DResourcePool< DKinFitConstraint_Mass > dResourcePool_MassConstraint
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
const DVector3 & position(void) const
Definition: GlueX.h:24
const DTrackTimeBased * Get_TrackTimeBased(void) const
map< bool, TH2I * > dHistMap_TOFPointTrackDeltaXVsVerticalPaddle
TDirectoryFile * CreateAndChangeTo_ActionDirectory(void)
void Read_MemoryUsage(double &vm_usage, double &resident_set)
unsigned int potential_cdc_hits_on_track
map< DetectorSystem_t, map< Particle_t, TH2I * > > dHistMap_DeltadEdXVsP
size_t Get_NumTrackFCALMatches(void) const
size_t Get_NumTrackSCMatches(void) const
DResourcePool< DKinematicData > dResourcePool_KinematicData
static char * ParticleType(Particle_t p)
Definition: particleType.h:142
DetectorSystem_t
Definition: GlueX.h:15
unsigned int dFDCPlanes
double measuredBeta(void) const
double Calc_PropagatedRFTime(const DKinematicData *locKinematicData, const DEventRFBunch *locEventRFBunch) const
static int ParticleCharge(Particle_t p)
#define X(str)
Definition: hddm-c.cpp:83
Definition: GlueX.h:17
Definition: GlueX.h:19
shared_ptr< const DTOFHitMatchParams > Get_TOFHitMatchParams(void) const
map< int, TH2I * > dHistMap_PVsTheta_GoodWireBased_BadTimeBased
JApplication * japp
bool PredictFCALHit(const DReferenceTrajectory *rt, unsigned int &row, unsigned int &col, DVector3 *intersection=nullptr) const
DetectorSystem_t dDetectorSystem
map< Particle_t, TH2I * > dHistMap_Ldiff_kpiVsP_DIRC
unsigned int dFDCPlanes
map< Particle_t, TH2I * > dHistMap_ThetaCVsP_DIRC
double time(void) const
map< DetectorSystem_t, map< bool, TH2I * > > dHistMap_PhiVsTheta_NoHit
map< DetectorSystem_t, map< int, TH2I * > > dHistMap_dEdXVsP
bool Get_BestSCMatchParams(const DTrackingData *locTrack, const DDetectorMatches *locDetectorMatches, shared_ptr< const DSCHitMatchParams > &locBestMatchParams) const
map< bool, TH2I * > dHistMap_TOFPaddleTrackYVsVerticalPaddle_HasHit
map< bool, TH2I * > dHistMap_TOFPointTrackDeltaXVsHorizontalPaddle
bool PredictTOFPaddles(const DReferenceTrajectory *rt, unsigned int &hbar, unsigned int &vbar, DVector3 *intersection=nullptr) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
static double E3[100]
map< pair< int, bool >, TH2I * > dHistMap_FCALShowerTrackDepthVsP
void Initialize(JEventLoop *locEventLoop)
void Fill_Hists(JEventLoop *locEventLoop, bool locUseTruePIDFlag)
Definition: GlueX.h:20
Definition: GlueX.h:18
map< const JObject *, map< DKinFitPullType, double > > dKinFitPulls
Definition: DVertex.h:33
map< bool, TH2I * > dHistMap_TOFPointTrackDeltaYVsHorizontalPaddle
map< DetectorSystem_t, map< bool, TH2I * > > dHistMap_PVsTheta_NoHit
double GetMostProbabledEdx_DC(double p, double mass, double dx, bool locIsCDCFlag) const
Definition: DParticleID.cc:657
map< pair< int, bool >, TH1I * > dHistMap_BCALShowerEnergy
DGeometry * GetDGeometry(unsigned int run_number)
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
Definition: GlueX.h:22
shared_ptr< TMatrixFSym > dCovarianceMatrix
map< DetectorSystem_t, map< Particle_t, TH2I * > > dHistMap_DeltaTVsP
const DTOFPaddleHit * Get_ClosestTOFPaddleHit_Vertical(const DReferenceTrajectory *locReferenceTrajectory, const vector< const DTOFPaddleHit * > &locTOFPaddleHits, double locInputStartTime, double &locBestDeltaX, double &locBestDistance) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
DResourcePool< DNeutralParticleHypothesis > dResourcePool_NeutralParticleHypothesis
void Initialize(JEventLoop *locEventLoop)
Definition: DCutActions.h:512
DResourcePool< DChargedTrackHypothesis > dResourcePool_ChargedTrackHypothesis
shared_ptr< const DDIRCMatchParams > Get_DIRCMatchParams(void) const
DResourcePool< DSourceCombo > dResourcePool_SourceCombo
bool Get_IsMatchedToHit(const DTrackingData *locTrack) const
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
unsigned int potential_fdc_hits_on_track
bool Get_ClosestToTrack(const DReferenceTrajectory *rt, const vector< const DBCALShower * > &locBCALShowers, bool locCutFlag, double &locStartTime, shared_ptr< const DBCALShowerMatchParams > &locBestMatchParams, double *locStartTimeVariance=nullptr, DVector3 *locBestProjPos=nullptr, DVector3 *locBestProjMom=nullptr) const
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
DLorentzVector dSpacetimeVertex
Definition: DVertex.h:28
DResourcePool< DEventRFBunch > dResourcePool_EventRFBunch
shared_ptr< const DSCHitMatchParams > Get_SCHitMatchParams(void) const
const char * SystemName(DetectorSystem_t sys)
Definition: GlueX.h:38
DResourcePool< DKinFitConstraint_Spacetime > dResourcePool_SpacetimeConstraint
DResourcePool< DChargedTrackHypothesis::DTimingInfo > dResourcePool_ChargedHypoTimingInfo
void Initialize(JEventLoop *locEventLoop)
map< Particle_t, TH2I * > dHistMap_TrackPzErrorVsTheta
DLorentzVector lorentzMomentum(void) const
map< pair< int, bool >, TH1I * > dHistMap_FCALShowerTrackDepth
map< pair< int, bool >, TH1I * > dHistMap_FCALShowerEnergy
DResourcePool< DKinFitChainStep > dResourcePool_KinFitChainStep
double sqrt(double)
shared_ptr< const DFCALShowerMatchParams > Get_FCALShowerMatchParams(void) const
DResourcePool< DKinFitParticle > dResourcePool_KinFitParticle
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo=NULL)
unsigned int dCDCRings
map< pair< int, bool >, TH1I * > dHistMap_BCALShowerTrackDepth
int dVerticalBar
Definition: DTOFPoint.h:39
DResourcePool< DChargedTrackHypothesis::DTrackingInfo > dResourcePool_ChargedHypoTrackingInfo
bool Get_CalledPriorWithComboFlag(void) const
unsigned int dCDCRings
bool Get_Extrapolations(const DKinematicData *locTrack, map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > &extrapolations) const
DResourcePool< vector< const DSourceCombo * > > dResourcePool_SourceComboVector
size_t Get_NumTrackBCALMatches(void) const
void Initialize(JEventLoop *locEventLoop)
size_t Get_NumTrackTOFMatches(void) const
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
map< pair< int, bool >, TH2I * > dHistMap_SCPhiVsTheta
DResourcePool< DKinematicData::DKinematicInfo > dResourcePool_KinematicInfo
const DVector3 & momentum(void) const
bool Cut_TrackHitPattern(const DParticleID *locParticleID, const DKinematicData *locTrack) const
Definition: DCutActions.cc:809
map< int, TH2I * > dHistMap_PVsTheta_GoodWireBased_GoodTimeBased
map< DetectorSystem_t, map< bool, TH2I * > > dHistMap_PhiVsTheta_HasHit
size_t Get_NumObjectsAllThreads(void) const
map< Particle_t, TH2I * > dHistMap_Ldiff_pkVsP_DIRC
DetectorSystem_t dTimeSource
Definition: DEventRFBunch.h:28
void Initialize(JEventLoop *locEventLoop)
map< Particle_t, map< DKinFitPullType, TH1I * > > dHistMap_KinFitPulls
map< DetectorSystem_t, map< Particle_t, TH2I * > > dHistMap_DeltaBetaVsP
DResourcePool< DKinFitChain > dResourcePool_KinFitChain
DResourcePool< DBeamPhoton > dResourcePool_BeamPhotons
DResourcePool< DKinFitResults > dResourcePool_KinFitResults
shared_ptr< const DBCALShowerMatchParams > Get_BCALShowerMatchParams(void) const
double dKinFitChiSq
Definition: DVertex.h:32
map< Particle_t, TH2I * > dHistMap_TrackPxErrorVsTheta
DResourcePool< TMatrixFSym > dResourcePool_TMatrixFSym
shared_ptr< const TMatrixFSym > errorMatrix(void) const
DResourcePool< DKinFitConstraint_Vertex > dResourcePool_VertexConstraint
unsigned int dNumParticleVotes
Definition: DEventRFBunch.h:32
map< bool, TH2I * > dHistMap_TOFPaddleHorizontalPaddleVsTrackX_HasHit
void GetScintMPdEandSigma(double p, double M, double x, double &most_prob_dE, double &sigma_dE) const
Definition: DParticleID.cc:622
double Step_TimeToNearInputTime(double locTimeToStep, double locTimeToStepTo) const
DResourcePool< DNeutralParticleHypothesis::DTimingInfo > dResourcePool_NeutralHypoTimingInfo
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
bool PredictBCALWedge(const DReferenceTrajectory *rt, unsigned int &module, unsigned int &sector, DVector3 *intersection=nullptr) const
map< bool, TH2I * > dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
Particle_t PID(void) const
bool Get_IsMatchedToDetector(const DTrackingData *locTrack, DetectorSystem_t locDetectorSystem) const
void Fill_MatchingHists(JEventLoop *locEventLoop, bool locIsTimeBased)
DResourcePool< DKinFitConstraint_P4 > dResourcePool_P4Constraint
double mass(void) const
void Initialize(JEventLoop *locEventLoop)
Particle_t
Definition: particleType.h:12