Bug Summary

File:libraries/ANALYSIS/DHistogramActions_Independent.cc
Location:line 203, column 3
Description:Value stored to 'locBaseFactory' is never read

Annotated Source Code

1#include "ANALYSIS/DHistogramActions.h"
2
3void DHistogramAction_ObjectMemory::Initialize(JEventLoop* locEventLoop)
4{
5 string locHistName, locHistTitle;
6 vector<string> locBinLabels; //fill this
7
8 //ANALYSIS
9 dFactoryPoolBinMap["DParticleComboBlueprintStep"] = 1;
10 locBinLabels.push_back("DParticleComboBlueprintStep");
11
12 dFactoryPairsToTrack.push_back(pair<string, string>("DParticleComboBlueprint", ""));
13 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 2;
14 locBinLabels.push_back("DParticleComboBlueprint");
15
16 dFactoryPairsToTrack.push_back(pair<string, string>("DTrackTimeBased", "Combo"));
17 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 3;
18 locBinLabels.push_back("DTrackTimeBased_Combo");
19
20 dFactoryPairsToTrack.push_back(pair<string, string>("DEventRFBunch", "Combo"));
21 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 4;
22 locBinLabels.push_back("DEventRFBunch_Combo");
23
24 dFactoryPairsToTrack.push_back(pair<string, string>("DChargedTrackHypothesis", "Combo"));
25 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 5;
26 locBinLabels.push_back("DChargedTrackHypothesis_Combo");
27
28 dFactoryPairsToTrack.push_back(pair<string, string>("DNeutralParticleHypothesis", "Combo"));
29 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 6;
30 locBinLabels.push_back("DNeutralParticleHypothesis_Combo");
31
32 dFactoryPoolBinMap["DParticleComboStep_PreKinFit"] = 7;
33 locBinLabels.push_back("DParticleComboStep_PreKinFit");
34
35 dFactoryPoolBinMap["DKinematicData_ComboPreKinFit"] = 8;
36 locBinLabels.push_back("DKinematicData_ComboPreKinFit");
37
38 dFactoryPairsToTrack.push_back(pair<string, string>("DParticleCombo", "PreKinFit"));
39 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 9;
40 locBinLabels.push_back("DParticleCombo_PreKinFit");
41
42 dFactoryPoolBinMap["DKinFitParticle"] = 10;
43 locBinLabels.push_back("DKinFitParticle");
44
45 dFactoryPoolBinMap["DKinFitConstraint_Vertex"] = 11;
46 locBinLabels.push_back("DKinFitConstraint_Vertex");
47
48 dFactoryPoolBinMap["DKinFitConstraint_Spacetime"] = 12;
49 locBinLabels.push_back("DKinFitConstraint_Spacetime");
50
51 dFactoryPoolBinMap["DKinFitConstraint_P4"] = 13;
52 locBinLabels.push_back("DKinFitConstraint_P4");
53
54 dFactoryPoolBinMap["TMatrixDSym_KinFitter"] = 14;
55 locBinLabels.push_back("TMatrixDSym_KinFitter");
56
57 dFactoryPairsToTrack.push_back(pair<string, string>("DKinFitResults", ""));
58 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 15;
59 locBinLabels.push_back("DKinFitResults");
60
61 dFactoryPairsToTrack.push_back(pair<string, string>("DBeamPhoton", "KinFit"));
62 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 16;
63 locBinLabels.push_back("DBeamPhoton_KinFit");
64
65 dFactoryPairsToTrack.push_back(pair<string, string>("DChargedTrackHypothesis", "KinFit"));
66 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 17;
67 locBinLabels.push_back("DChargedTrackHypothesis_KinFit");
68
69 dFactoryPairsToTrack.push_back(pair<string, string>("DNeutralParticleHypothesis", "KinFit"));
70 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 18;
71 locBinLabels.push_back("DNeutralParticleHypothesis_KinFit");
72
73 dFactoryPoolBinMap["DParticleComboStep"] = 19;
74 locBinLabels.push_back("DParticleComboStep");
75
76 dFactoryPoolBinMap["DKinematicData_Combo"] = 20;
77 locBinLabels.push_back("DKinematicData_Combo");
78
79 dFactoryPairsToTrack.push_back(pair<string, string>("DParticleCombo", ""));
80 dFactoryPairBinMap[dFactoryPairsToTrack.back()] = 21;
81 locBinLabels.push_back("DParticleCombo");
82
83 //CREATE THE HISTOGRAMS
84 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
85 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
86 {
87 CreateAndChangeTo_ActionDirectory();
88
89 // Total Memory
90 locHistName = "TotalMemory";
91 locHistTitle = ";Event # ;Total Memory (MB)";
92 dHist_TotalMemory = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5);
93
94 // # Objects
95 locHistName = "NumObjects2D";
96 locHistTitle = "# Objects;Event #";
97 dHist_NumObjects = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5, locBinLabels.size(), 0.5, float(locBinLabels.size()) + 0.5);
98 for(size_t loc_i = 0; loc_i < locBinLabels.size(); ++loc_i)
99 dHist_NumObjects->GetYaxis()->SetBinLabel(1 + loc_i, locBinLabels[loc_i].c_str());
100
101 // Object Memory
102 locHistName = "Memory2D";
103 locHistTitle = "Memory (Bytes);Event #";
104 dHist_Memory = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5, locBinLabels.size(), 0.5, float(locBinLabels.size()) + 0.5);
105 for(size_t loc_i = 0; loc_i < locBinLabels.size(); ++loc_i)
106 dHist_Memory->GetYaxis()->SetBinLabel(1 + loc_i, locBinLabels[loc_i].c_str());
107
108 for(size_t loc_i = 0; loc_i < locBinLabels.size(); ++loc_i)
109 {
110 // # Objects
111 locHistName = string("NumObjects_") + locBinLabels[loc_i];
112 locHistTitle = locBinLabels[loc_i] + string(";Event # ;# Objects");
113 dHistMap_NumObjects[loc_i + 1] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5);
114
115 // # Objects
116 locHistName = string("Memory_") + locBinLabels[loc_i];
117 locHistTitle = locBinLabels[loc_i] + string(";Event # ;Memory (Bytes)");
118 dHistMap_Memory[loc_i + 1] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dMaxNumEvents, 0.5, float(dMaxNumEvents) + 0.5);
119 }
120
121 //Return to the base directory
122 ChangeTo_BaseDirectory();
123 }
124 japp->RootUnLock(); //RELEASE ROOT LOCK!!
125}
126
127bool DHistogramAction_ObjectMemory::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
128{
129 if(Get_NumPreviousParticleCombos() != 0)
130 return true; //else double-counting!
131
132 if(dEventCounter > dMaxNumEvents)
133 return true;
134
135 if(locParticleCombo != NULL__null)
136 return true; // Protect against infinite recursion (see below)
137
138 //THIS IS EXTREMELY DANGEROUS, AND SHOULD BE AVOIDED UNLESS YOU KNOW !!!!EXACTLY!!!! WHAT YOU ARE DOING
139 //And if you're doing this, you probably don't
140 //This will result in a infinite-recursion crash if it is called by DAnalysisResults_factory.
141 //This action should only be used directly in an event processsor
142 //This casuses the analysis to run, generating the objects needed for histogramming below.
143 vector<const DAnalysisResults*> locAnalysisResults;
144 locEventLoop->Get(locAnalysisResults);
145
146 //FACTORIES
147 //call get-n-rows first outside of lock, just to make sure
148 map<int, size_t> locNumObjectsMap; //int is bin
149 map<int, unsigned long long> locMemoryMap; //int is bin
150 double locTotalMemory = 0;
151 for(size_t loc_i = 0; loc_i < dFactoryPairsToTrack.size(); ++loc_i)
152 {
153 string locClassName = dFactoryPairsToTrack[loc_i].first;
154 JFactory_base* locFactory = locEventLoop->GetFactory(locClassName.c_str(), dFactoryPairsToTrack[loc_i].second.c_str());
155 unsigned long long locNumObjects = locFactory->GetNrows();
156 unsigned long long locDataClassSize = locFactory->GetDataClassSize();
157
158 unsigned long long locMemory = locDataClassSize*locNumObjects;
159 if(locClassName == "DChargedTrackHypothesis")
160 locMemory += locNumObjects*(7*7*8 + 5*5*8); //error matrices //8 = double
161 if((locClassName == "DNeutralParticleHypothesis") || (locClassName == "DBeamPhoton"))
162 locMemory += locNumObjects*(7*7*8); //error matrices //8 = double
163
164 int locBin = dFactoryPairBinMap[dFactoryPairsToTrack[loc_i]];
165 locNumObjectsMap[locBin] = locNumObjects;
166 locMemoryMap[locBin] = locMemory;
167 locTotalMemory += double(locMemory);
168 }
169
170 //RESOURCE POOLS
171 {
172 unsigned long long locMemory;
173 JFactory_base* locBaseFactory;
174 int locBin;
175
176 //DParticleComboBlueprintStep
177 locBin = dFactoryPoolBinMap["DParticleComboBlueprintStep"];
178 locBaseFactory = locEventLoop->GetFactory("DParticleComboBlueprint", "");
179 DParticleComboBlueprint_factory* locParticleComboBlueprintFactory = static_cast<DParticleComboBlueprint_factory*>(locBaseFactory);
180 locNumObjectsMap[locBin] = locParticleComboBlueprintFactory->Get_ParticleComboBlueprintStepPoolSize();
181 locMemory = sizeof(DParticleComboBlueprintStep)*locNumObjectsMap[locBin];
182 locMemoryMap[locBin] = locMemory;
183 locTotalMemory += double(locMemory);
184
185 //DParticleComboStep_PreKinFit
186 locBin = dFactoryPoolBinMap["DParticleComboStep_PreKinFit"];
187 locBaseFactory = locEventLoop->GetFactory("DParticleCombo", "PreKinFit");
188 DParticleCombo_factory_PreKinFit* locParticleComboFactory_PreKinFit = static_cast<DParticleCombo_factory_PreKinFit*>(locBaseFactory);
189 locNumObjectsMap[locBin] = locParticleComboFactory_PreKinFit->Get_ParticleComboStepPoolSize();
190 locMemory = sizeof(DParticleComboStep)*locNumObjectsMap[locBin];
191 locMemoryMap[locBin] = locMemory;
192 locTotalMemory += double(locMemory);
193
194 //DKinematicData_ComboPreKinFit
195 locBin = dFactoryPoolBinMap["DKinematicData_ComboPreKinFit"];
196 locNumObjectsMap[locBin] = locParticleComboFactory_PreKinFit->Get_KinematicDataPoolSize();
197 locMemory = sizeof(DKinematicData)*locNumObjectsMap[locBin];
198 locMemoryMap[locBin] = locMemory;
199 locTotalMemory += double(locMemory);
200
201 //DKinFitParticle
202 locBin = dFactoryPoolBinMap["DKinFitParticle"];
203 locBaseFactory = locEventLoop->GetFactory("DKinFitResults", "");
Value stored to 'locBaseFactory' is never read
204// DKinFitResults_factory* locKinFitResultsFactory = static_cast<DKinFitResults_factory*>(locBaseFactory);
205// locNumObjectsMap[locBin] = locKinFitResultsFactory->Get_KinFitParticlePoolSize();
206 locNumObjectsMap[locBin] = 1;
207 locMemory = sizeof(DKinFitParticle)*locNumObjectsMap[locBin];
208 locMemoryMap[locBin] = locMemory;
209 locTotalMemory += double(locMemory);
210
211 //DKinFitConstraint_Vertex
212 locBin = dFactoryPoolBinMap["DKinFitConstraint_Vertex"];
213// locNumObjectsMap[locBin] = locKinFitResultsFactory->Get_KinFitConstraintVertexPoolSize();
214 locNumObjectsMap[locBin] = 1;
215 locMemory = sizeof(DKinFitConstraint_Vertex)*locNumObjectsMap[locBin];
216 locMemoryMap[locBin] = locMemory;
217 locTotalMemory += double(locMemory);
218
219 //DKinFitConstraint_Spacetime
220 locBin = dFactoryPoolBinMap["DKinFitConstraint_Spacetime"];
221// locNumObjectsMap[locBin] = locKinFitResultsFactory->Get_KinFitConstraintSpacetimePoolSize();
222 locNumObjectsMap[locBin] = 1;
223 locMemory = sizeof(DKinFitConstraint_Spacetime)*locNumObjectsMap[locBin];
224 locMemoryMap[locBin] = locMemory;
225 locTotalMemory += double(locMemory);
226
227 //DKinFitConstraint_P4
228 locBin = dFactoryPoolBinMap["DKinFitConstraint_P4"];
229// locNumObjectsMap[locBin] = locKinFitResultsFactory->Get_KinFitConstraintP4PoolSize();
230 locNumObjectsMap[locBin] = 1;
231 locMemory = sizeof(DKinFitConstraint_P4)*locNumObjectsMap[locBin];
232 locMemoryMap[locBin] = locMemory;
233 locTotalMemory += double(locMemory);
234
235 //TMatrixDSym_KinFitter
236 locBin = dFactoryPoolBinMap["TMatrixDSym_KinFitter"];
237// locNumObjectsMap[locBin] = locKinFitResultsFactory->Get_MatrixDSymPoolSize() + locKinFitResultsFactory->Get_LargeMatrixDSymPoolSize();
238 locNumObjectsMap[locBin] = 1;
239// locMemory = (sizeof(TMatrixDSym) + 7*7*8)*locKinFitResultsFactory->Get_MatrixDSymPoolSize(); //assume 7x7 matrix of doubles (8)
240// locMemory += (sizeof(TMatrixDSym) + 30*30*8)*locKinFitResultsFactory->Get_LargeMatrixDSymPoolSize(); //assume 30x30 matrix of doubles (8)
241 locMemoryMap[locBin] = locMemory;
242 locTotalMemory += double(locMemory);
243
244 //DParticleComboStep
245 locBin = dFactoryPoolBinMap["DParticleComboStep"];
246 locBaseFactory = locEventLoop->GetFactory("DParticleCombo", "");
247 DParticleCombo_factory* locParticleComboFactory = static_cast<DParticleCombo_factory*>(locBaseFactory);
248 locNumObjectsMap[locBin] = locParticleComboFactory->Get_ParticleComboStepPoolSize();
249 locMemory = sizeof(DParticleComboStep)*locNumObjectsMap[locBin];
250 locMemoryMap[locBin] = locMemory;
251 locTotalMemory += double(locMemory);
252
253 //DKinematicData_Combo
254 locBin = dFactoryPoolBinMap["DKinematicData_Combo"];
255 locNumObjectsMap[locBin] = locParticleComboFactory->Get_KinematicDataPoolSize();
256 locMemory = sizeof(DKinematicData)*locNumObjectsMap[locBin];
257 locMemoryMap[locBin] = locMemory;
258 locTotalMemory += double(locMemory);
259 }
260 locTotalMemory /= (1024.0*1024.0); //convert to MB
261
262 map<int, size_t>::iterator locIterator = locNumObjectsMap.begin();
263
264 //FILL HISTOGRAMS
265 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
266 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
267 Lock_Action(); //ACQUIRE ROOT LOCK!!
268 {
269 ++dEventCounter;
270 if(dEventCounter <= dMaxNumEvents)
271 {
272 for(; locIterator != locNumObjectsMap.end(); ++locIterator)
273 {
274 int locObjectBin = locIterator->first;
275
276 dHistMap_NumObjects[locObjectBin]->SetBinContent(dEventCounter, locNumObjectsMap[locObjectBin]);
277 dHist_NumObjects->SetBinContent(dEventCounter, locObjectBin, locNumObjectsMap[locObjectBin]);
278
279 dHistMap_Memory[locObjectBin]->SetBinContent(dEventCounter, locMemoryMap[locObjectBin]);
280 dHist_Memory->SetBinContent(dEventCounter, locObjectBin, locMemoryMap[locObjectBin]);
281 }
282 dHist_TotalMemory->SetBinContent(dEventCounter, locTotalMemory);
283 }
284 }
285 Unlock_Action(); //RELEASE ROOT LOCK!!
286
287 return true;
288}
289
290void DHistogramAction_Reconstruction::Initialize(JEventLoop* locEventLoop)
291{
292 //Create any histograms/trees/etc. within a ROOT lock.
293 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
294
295 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
296 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
297
298 string locHistName, locHistTitle;
299
300 bool locIsRESTEvent = (string(locEventLoop->GetJEvent().GetJEventSource()->className()) == string("DEventSourceREST"));
301
302 vector<const DMCThrown*> locMCThrowns;
303 locEventLoop->Get(locMCThrowns);
304
305 //CREATE THE HISTOGRAMS
306 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
307 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
308 {
309 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
310 //If another thread has already created the folder, it just changes to it.
311 CreateAndChangeTo_ActionDirectory();
312
313 //FCAL
314 locHistName = "FCALShowerYVsX";
315 dHist_FCALShowerYVsX = GetOrCreate_Histogram<TH2I>(locHistName, ";FCAL Shower X (cm);FCAL Shower Y (cm)", dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
316 locHistName = "FCALShowerEnergy";
317 dHist_FCALShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
318
319 //BCAL
320 locHistName = "BCALShowerEnergy";
321 dHist_BCALShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxBCALP);
322 locHistName = "BCALShowerPhi";
323 dHist_BCALShowerPhi = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #phi#circ", dNumPhiBins, dMinPhi, dMaxPhi);
324 locHistName = "BCALShowerPhiVsZ";
325 dHist_BCALShowerPhiVsZ = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Shower Z (cm);BCAL Shower #phi#circ", dNum2DBCALZBins, 0.0, 450.0, dNum2DPhiBins, dMinPhi, dMaxPhi);
326
327 //TOF
328 locHistName = "TOFPointEnergy";
329 dHist_TOFPointEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";TOF Point Energy (MeV)", dNumHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
330 locHistName = "TOFPointYVsX";
331 dHist_TOFPointYVsX = GetOrCreate_Histogram<TH2I>(locHistName, ";TOF Point X (cm);TOF Point Y (cm)", dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
332
333 //SC
334 locHistName = "SCHitSector";
335 dHist_SCHitSector = GetOrCreate_Histogram<TH1I>(locHistName, ";SC Hit Sector", 30, 0.5, 30.5);
336 locHistName = "SCHitEnergy";
337 dHist_SCHitEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";SC Hit Energy (MeV)", dNumHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
338 locHistName = "SCHitEnergyVsSector";
339 dHist_SCHitEnergyVsSector = GetOrCreate_Histogram<TH2I>(locHistName, ";SC Hit Sector;SC Hit Energy (MeV)", 30, 0.5, 30.5, dNum2DHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
340
341 //TRACKING
342 CreateAndChangeTo_Directory("Tracking", "Tracking");
343 locHistName = "NumDCHitsPerTrack";
344 dHist_NumDCHitsPerTrack = GetOrCreate_Histogram<TH1I>(locHistName, ";# Track Hits", 50, 0.5, 50.5);
345 locHistName = "NumDCHitsPerTrackVsTheta";
346 dHist_NumDCHitsPerTrackVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, ";#theta#circ;# Track Hits", dNum2DThetaBins, dMinTheta, dMaxTheta, 46, 4.5, 50.5);
347 locHistName = "TrackingFOM_WireBased";
348 dHist_TrackingFOM_WireBased = GetOrCreate_Histogram<TH1I>(locHistName, ";Confidence Level", dNumFOMBins, 0.0, 1.0);
349 locHistName = "TrackingFOM";
350 dHist_TrackingFOM = GetOrCreate_Histogram<TH1I>(locHistName, ";Confidence Level", dNumFOMBins, 0.0, 1.0);
351 locHistName = "TrackingFOMVsP";
352 dHist_TrackingFOMVsP = GetOrCreate_Histogram<TH2I>(locHistName, ";p (GeV/c);Confidence Level", dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
353 locHistName = "TrackingFOMVsTheta";
354 dHist_TrackingFOMVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, ";#theta#circ;Confidence Level", dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DFOMBins, 0.0, 1.0);
355 locHistName = "TrackingFOMVsNumHits";
356 dHist_TrackingFOMVsNumHits = GetOrCreate_Histogram<TH2I>(locHistName, ";# Track Hits;Confidence Level", 46, 4.5, 50.5, dNum2DFOMBins, 0.0, 1.0);
357
358 if(!locIsRESTEvent)
359 {
360 locHistName = "CDCRingVsTheta_Candidates";
361 dHist_CDCRingVsTheta_Candidates = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Track Candidates;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
362 locHistName = "CDCRingVsTheta_WireBased";
363 dHist_CDCRingVsTheta_WireBased = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Wire-Based Tracks;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
364 }
365
366 locHistName = "CDCRingVsTheta_TimeBased";
367 dHist_CDCRingVsTheta_TimeBased = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Time-Based Tracks;#theta#circ;CDC Ring", dNum2DThetaBins, dMinTheta, dMaxTheta, 28, 0.5, 28.5);
368 locHistName = "CDCRingVsTheta_TimeBased_GoodTrackFOM";
369 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);
370
371 if(!locIsRESTEvent)
372 {
373 locHistName = "FDCPlaneVsTheta_Candidates";
374 dHist_FDCPlaneVsTheta_Candidates = GetOrCreate_Histogram<TH2I>(locHistName, "Hits on Track Candidates;p (GeV/c);FDC Plane", dNum2DThetaBins, dMinTheta, dMaxTheta, 24, 0.5, 24.5);
375 locHistName = "FDCPlaneVsTheta_WireBased";
376 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);
377 }
378
379 locHistName = "FDCPlaneVsTheta_TimeBased";
380 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);
381 locHistName = "FDCPlaneVsTheta_TimeBased_GoodTrackFOM";
382 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);
383
384 if(!locMCThrowns.empty())
385 {
386 locHistName = "MCMatchedHitsVsTheta";
387 dHist_MCMatchedHitsVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, "Fraction of Track Hits Matched to MC;#theta#circ;", dNum2DThetaBins, dMinTheta, dMaxTheta, 100, 0.0, 1.0);
388 locHistName = "MCMatchedHitsVsP";
389 dHist_MCMatchedHitsVsP = GetOrCreate_Histogram<TH2I>(locHistName, "Fraction of Track Hits Matched to MC;p (GeV/c);", dNum2DPBins, dMinP, dMaxP, 100, 0.0, 1.0);
390 }
391
392 for(int locCharge = -1; locCharge <= 1; locCharge += 2)
393 {
394 string locParticleROOTName = (locCharge == -1) ? "#it{q}^{-}" : "#it{q}^{+}";
395 string locParticleName = (locCharge == -1) ? "q-" : "q+";
396 CreateAndChangeTo_Directory(locParticleName, locParticleName);
397
398 if(!locIsRESTEvent)
399 {
400 // PVsTheta Track Candidates
401 locHistName = string("PVsTheta_Candidates_") + locParticleName;
402 locHistTitle = locParticleROOTName + string(" Track Candidates;#theta#circ;p (GeV/c)");
403 dHistMap_PVsTheta_Candidates[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
404
405 // PVsTheta Wire-Based Tracks
406 locHistName = string("PVsTheta_WireBased_") + locParticleName;
407 locHistTitle = locParticleROOTName + string(" Wire-Based Tracks;#theta#circ;p (GeV/c)");
408 dHistMap_PVsTheta_WireBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
409 }
410
411 // PVsTheta Time-Based Tracks
412 locHistName = string("PVsTheta_TimeBased_") + locParticleName;
413 locHistTitle = locParticleROOTName + string(" Time-Based Tracks;#theta#circ;p (GeV/c)");
414 dHistMap_PVsTheta_TimeBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
415
416 // PVsTheta Time-Based Tracks Good Track FOM
417 locHistName = string("PVsTheta_TimeBased_GoodTrackFOM_") + locParticleName;
418 locHistTitle = locParticleROOTName + string(" Time-Based Tracks, Good Tracking FOM;#theta#circ;p (GeV/c)");
419 dHistMap_PVsTheta_TimeBased_GoodTrackFOM[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
420
421 // PVsTheta Time-Based Tracks Low Track FOM
422 locHistName = string("PVsTheta_TimeBased_LowTrackFOM_") + locParticleName;
423 locHistTitle = locParticleROOTName + string(" Time-Based Tracks, Low Tracking FOM;#theta#circ;p (GeV/c)");
424 dHistMap_PVsTheta_TimeBased_LowTrackFOM[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
425
426 // PVsTheta Time-Based Tracks High Track FOM
427 locHistName = string("PVsTheta_TimeBased_HighTrackFOM_") + locParticleName;
428 locHistTitle = locParticleROOTName + string(" Time-Based Tracks, High Tracking FOM;#theta#circ;p (GeV/c)");
429 dHistMap_PVsTheta_TimeBased_HighTrackFOM[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
430
431 if(!locIsRESTEvent)
432 {
433 // PVsTheta: Good Wire-Based, Good Time-Based
434 locHistName = string("PVsTheta_GoodWireBased_GoodTimeBased_") + locParticleName;
435 locHistTitle = locParticleROOTName + string(" Good Wire-Based, Good Time-Based;#theta#circ;p (GeV/c)");
436 dHistMap_PVsTheta_GoodWireBased_GoodTimeBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
437
438 // PVsTheta: Good Wire-Based, Bad Time-Based
439 locHistName = string("PVsTheta_GoodWireBased_BadTimeBased_") + locParticleName;
440 locHistTitle = locParticleROOTName + string(" Good Wire-Based, Bad Time-Based;#theta#circ;p (GeV/c)");
441 dHistMap_PVsTheta_GoodWireBased_BadTimeBased[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
442 }
443
444 gDirectory(TDirectory::CurrentDirectory())->cd(".."); //end of charge
445 }
446 gDirectory(TDirectory::CurrentDirectory())->cd(".."); //End of "Tracking"
447
448 //Return to the base directory
449 ChangeTo_BaseDirectory();
450 }
451 japp->RootUnLock(); //RELEASE ROOT LOCK!!
452}
453
454bool DHistogramAction_Reconstruction::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
455{
456 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
457
458 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
459 if(Get_NumPreviousParticleCombos() != 0)
460 return true;
461
462 bool locIsRESTEvent = (string(locEventLoop->GetJEvent().GetJEventSource()->className()) == string("DEventSourceREST"));
463
464 vector<const DBCALShower*> locBCALShowers;
465 locEventLoop->Get(locBCALShowers);
466
467 vector<const DFCALShower*> locFCALShowers;
468 locEventLoop->Get(locFCALShowers);
469
470 vector<const DTOFPoint*> locTOFPoints;
471 locEventLoop->Get(locTOFPoints);
472
473 vector<const DSCHit*> locSCHits;
474 locEventLoop->Get(locSCHits);
475
476 const DDetectorMatches* locDetectorMatches = NULL__null;
477 locEventLoop->GetSingle(locDetectorMatches);
478
479 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
480 locEventLoop->Get(locTrackTimeBasedVector);
481
482 vector<const DMCThrownMatching*> locMCThrownMatchingVector;
483 locEventLoop->Get(locMCThrownMatchingVector);
484
485 vector<const DMCThrown*> locMCThrowns;
486 locEventLoop->Get(locMCThrowns, "FinalState");
487
488 const DParticleID* locParticleID = NULL__null;
489 locEventLoop->GetSingle(locParticleID);
490
491 const DDetectorMatches* locDetectorMatches_WireBased = NULL__null;
492 vector<const DTrackCandidate*> locTrackCandidates;
493 vector<const DTrackWireBased*> locTrackWireBasedVector;
494 if(!locIsRESTEvent)
495 {
496 vector<const DDetectorMatches*> locDetectorMatchesVector_WireBased;
497 locEventLoop->Get(locDetectorMatchesVector_WireBased, "WireBased");
498 if(!locDetectorMatchesVector_WireBased.empty())
499 locDetectorMatches_WireBased = locDetectorMatchesVector_WireBased[0];
500 locEventLoop->Get(locTrackCandidates);
501 locEventLoop->Get(locTrackWireBasedVector);
502 }
503
504 //select the best DTrackWireBased for each track: use best tracking FOM
505 map<JObject::oid_t, const DTrackWireBased*> locBestTrackWireBasedMap; //lowest tracking FOM for each candidate id
506 for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
507 {
508 JObject::oid_t locCandidateID = locTrackWireBasedVector[loc_i]->candidateid;
509 if(locBestTrackWireBasedMap.find(locCandidateID) == locBestTrackWireBasedMap.end())
510 locBestTrackWireBasedMap[locCandidateID] = locTrackWireBasedVector[loc_i];
511 else if(locTrackWireBasedVector[loc_i]->FOM > locBestTrackWireBasedMap[locCandidateID]->FOM)
512 locBestTrackWireBasedMap[locCandidateID] = locTrackWireBasedVector[loc_i];
513 }
514
515 //select the best DTrackTimeBased for each track: use best tracking FOM
516 //also, make map from WBT -> TBT (if not rest)
517 map<JObject::oid_t, const DTrackTimeBased*> locBestTrackTimeBasedMap; //lowest tracking FOM for each candidate id
518 map<const DTrackWireBased*, const DTrackTimeBased*> locWireToTimeBasedTrackMap;
519 for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
520 {
521 JObject::oid_t locCandidateID = locTrackTimeBasedVector[loc_i]->candidateid;
522 if(locBestTrackTimeBasedMap.find(locCandidateID) == locBestTrackTimeBasedMap.end())
523 locBestTrackTimeBasedMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
524 else if(locTrackTimeBasedVector[loc_i]->FOM > locBestTrackTimeBasedMap[locCandidateID]->FOM)
525 locBestTrackTimeBasedMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
526 if(locIsRESTEvent)
527 continue;
528 const DTrackWireBased* locTrackWireBased = NULL__null;
529 locTrackTimeBasedVector[loc_i]->GetSingle(locTrackWireBased);
530 locWireToTimeBasedTrackMap[locTrackWireBased] = locTrackTimeBasedVector[loc_i];
531 }
532
533 //FILL HISTOGRAMS
534 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
535 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
536 Lock_Action(); //ACQUIRE ROOT LOCK!!
537 {
538 for(size_t loc_i = 0; loc_i < locFCALShowers.size(); ++loc_i)
539 {
540 dHist_FCALShowerEnergy->Fill(locFCALShowers[loc_i]->getEnergy());
541 dHist_FCALShowerYVsX->Fill(locFCALShowers[loc_i]->getPosition().X(), locFCALShowers[loc_i]->getPosition().Y());
542 }
543
544 for(size_t loc_i = 0; loc_i < locBCALShowers.size(); ++loc_i)
545 {
546 dHist_BCALShowerEnergy->Fill(locBCALShowers[loc_i]->E);
547
548 DVector3 locBCALPosition(locBCALShowers[loc_i]->x, locBCALShowers[loc_i]->y, locBCALShowers[loc_i]->z);
549 double locBCALPhi = locBCALPosition.Phi()*180.0/TMath::Pi();
550 dHist_BCALShowerPhi->Fill(locBCALPhi);
551 dHist_BCALShowerPhiVsZ->Fill(locBCALPosition.Z(), locBCALPhi);
552 }
553
554 for(size_t loc_i = 0; loc_i < locTOFPoints.size(); ++loc_i)
555 {
556 dHist_TOFPointEnergy->Fill(locTOFPoints[loc_i]->dE*1.0E3);
557 dHist_TOFPointYVsX->Fill(locTOFPoints[loc_i]->pos.X(), locTOFPoints[loc_i]->pos.Y());
558 }
559
560 for(size_t loc_i = 0; loc_i < locSCHits.size(); ++loc_i)
561 {
562 dHist_SCHitSector->Fill(locSCHits[loc_i]->sector);
563 dHist_SCHitEnergy->Fill(locSCHits[loc_i]->dE*1.0E3);
564 dHist_SCHitEnergyVsSector->Fill(locSCHits[loc_i]->sector, locSCHits[loc_i]->dE*1.0E3);
565 }
566
567 for(size_t loc_i = 0; loc_i < locTrackCandidates.size(); ++loc_i)
568 {
569 int locCharge = (locTrackCandidates[loc_i]->charge() > 0.0) ? 1 : -1;
570 double locTheta = locTrackCandidates[loc_i]->momentum().Theta()*180.0/TMath::Pi();
571 double locP = locTrackCandidates[loc_i]->momentum().Mag();
572 dHistMap_PVsTheta_Candidates[locCharge]->Fill(locTheta, locP);
573
574 set<int> locCDCRings;
575 locParticleID->Get_CDCRings(locTrackCandidates[loc_i]->dCDCRings, locCDCRings);
576 for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator)
577 dHist_CDCRingVsTheta_Candidates->Fill(locTheta, *locIterator);
578
579 set<int> locFDCPlanes;
580 locParticleID->Get_FDCPlanes(locTrackCandidates[loc_i]->dFDCPlanes, locFDCPlanes);
581 for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator)
582 dHist_FDCPlaneVsTheta_Candidates->Fill(locTheta, *locIterator);
583 }
584
585 map<JObject::oid_t, const DTrackWireBased*>::iterator locWireBasedIterator = locBestTrackWireBasedMap.begin();
586 for(; locWireBasedIterator != locBestTrackWireBasedMap.end(); ++locWireBasedIterator)
587 {
588 const DTrackWireBased* locTrackWireBased = locWireBasedIterator->second;
589 int locCharge = (locTrackWireBased->charge() > 0.0) ? 1 : -1;
590 double locTheta = locTrackWireBased->momentum().Theta()*180.0/TMath::Pi();
591 double locP = locTrackWireBased->momentum().Mag();
592 dHistMap_PVsTheta_WireBased[locCharge]->Fill(locTheta, locP);
593
594 set<int> locCDCRings;
595 locParticleID->Get_CDCRings(locTrackWireBased->dCDCRings, locCDCRings);
596 for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator)
597 dHist_CDCRingVsTheta_WireBased->Fill(locTheta, *locIterator);
598
599 set<int> locFDCPlanes;
600 locParticleID->Get_FDCPlanes(locTrackWireBased->dFDCPlanes, locFDCPlanes);
601 for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator)
602 dHist_FDCPlaneVsTheta_WireBased->Fill(locTheta, *locIterator);
603
604 dHist_TrackingFOM_WireBased->Fill(locTrackWireBased->FOM);
605 }
606
607 map<JObject::oid_t, const DTrackTimeBased*>::iterator locTimeBasedIterator = locBestTrackTimeBasedMap.begin();
608 for(; locTimeBasedIterator != locBestTrackTimeBasedMap.end(); ++locTimeBasedIterator)
609 {
610 const DTrackTimeBased* locTrackTimeBased = locTimeBasedIterator->second;
611 int locCharge = (locTrackTimeBased->charge() > 0.0) ? 1 : -1;
612 double locTheta = locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi();
613 double locP = locTrackTimeBased->momentum().Mag();
614
615 dHistMap_PVsTheta_TimeBased[locCharge]->Fill(locTheta, locP);
616 dHist_NumDCHitsPerTrack->Fill(locTrackTimeBased->Ndof + 5);
617 dHist_NumDCHitsPerTrackVsTheta->Fill(locTheta, locTrackTimeBased->Ndof + 5);
618
619 dHist_TrackingFOM->Fill(locTrackTimeBased->FOM);
620 dHist_TrackingFOMVsTheta->Fill(locTheta, locTrackTimeBased->FOM);
621 dHist_TrackingFOMVsP->Fill(locP, locTrackTimeBased->FOM);
622 dHist_TrackingFOMVsNumHits->Fill(locTrackTimeBased->Ndof + 5, locTrackTimeBased->FOM);
623
624 set<int> locCDCRings;
625 locParticleID->Get_CDCRings(locTrackTimeBased->dCDCRings, locCDCRings);
626 for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator)
627 {
628 dHist_CDCRingVsTheta_TimeBased->Fill(locTheta, *locIterator);
629 if(locTrackTimeBased->FOM > dGoodTrackFOM)
630 dHist_CDCRingVsTheta_TimeBased_GoodTrackFOM->Fill(locTheta, *locIterator);
631 }
632
633 set<int> locFDCPlanes;
634 locParticleID->Get_FDCPlanes(locTrackTimeBased->dFDCPlanes, locFDCPlanes);
635 for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator)
636 {
637 dHist_FDCPlaneVsTheta_TimeBased->Fill(locTheta, *locIterator);
638 if(locTrackTimeBased->FOM > dGoodTrackFOM)
639 dHist_FDCPlaneVsTheta_TimeBased_GoodTrackFOM->Fill(locTheta, *locIterator);
640 }
641
642 if(locTrackTimeBased->FOM > dGoodTrackFOM)
643 dHistMap_PVsTheta_TimeBased_GoodTrackFOM[locCharge]->Fill(locTheta, locP);
644 else
645 dHistMap_PVsTheta_TimeBased_LowTrackFOM[locCharge]->Fill(locTheta, locP);
646 if(locTrackTimeBased->FOM > dHighTrackFOM)
647 dHistMap_PVsTheta_TimeBased_HighTrackFOM[locCharge]->Fill(locTheta, locP);
648 }
649
650 // If "Good" WBT, see if TBT is good
651 locWireBasedIterator = locBestTrackWireBasedMap.begin();
652 for(; locWireBasedIterator != locBestTrackWireBasedMap.end(); ++locWireBasedIterator)
653 {
654 if(locDetectorMatches_WireBased == NULL__null)
655 continue;
656 const DTrackWireBased* locTrackWireBased = locWireBasedIterator->second;
657 if(locTrackWireBased->FOM < dGoodTrackFOM)
658 continue; //no good
659 if(!locDetectorMatches_WireBased->Get_IsMatchedToHit(locTrackWireBased))
660 continue; //no good
661
662 int locCharge = (locTrackWireBased->charge() > 0.0) ? 1 : -1;
663 double locTheta = locTrackWireBased->momentum().Theta()*180.0/TMath::Pi();
664 double locP = locTrackWireBased->momentum().Mag();
665
666 map<const DTrackWireBased*, const DTrackTimeBased*>::iterator locReconIterator = locWireToTimeBasedTrackMap.find(locTrackWireBased);
667 if(locReconIterator == locWireToTimeBasedTrackMap.end())
668 {
669 dHistMap_PVsTheta_GoodWireBased_BadTimeBased[locCharge]->Fill(locTheta, locP);
670 continue; //no time-based
671 }
672
673 const DTrackTimeBased* locTrackTimeBased = locReconIterator->second;
674 if((locTrackTimeBased->FOM < dGoodTrackFOM) || (!locDetectorMatches->Get_IsMatchedToHit(locTrackTimeBased)))
675 dHistMap_PVsTheta_GoodWireBased_BadTimeBased[locCharge]->Fill(locTheta, locP);
676 else
677 dHistMap_PVsTheta_GoodWireBased_GoodTimeBased[locCharge]->Fill(locTheta, locP);
678 }
679
680 for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
681 {
682 if(fabs(locMCThrowns[loc_i]->charge()) < 0.9)
683 continue;
684
685 double locMatchFOM;
686 const DChargedTrackHypothesis* locChargedTrackHypothesis = locMCThrownMatchingVector[0]->Get_MatchingChargedHypothesis(locMCThrowns[loc_i], locMatchFOM);
687 if(locChargedTrackHypothesis == NULL__null)
688 continue;
689
690 const DTrackTimeBased* locTrackTimeBased = NULL__null;
691 locChargedTrackHypothesis->GetSingle(locTrackTimeBased);
692
693 double locHitFraction = 1.0*locTrackTimeBased->dNumHitsMatchedToThrown/(locTrackTimeBased->Ndof + 5);
694 dHist_MCMatchedHitsVsTheta->Fill(locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi(), locHitFraction);
695 dHist_MCMatchedHitsVsP->Fill(locTrackTimeBased->momentum().Mag(), locHitFraction);
696 }
697 }
698 Unlock_Action(); //RELEASE ROOT LOCK!!
699
700 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
701}
702
703void DHistogramAction_DetectorMatching::Initialize(JEventLoop* locEventLoop)
704{
705 //Create any histograms/trees/etc. within a ROOT lock.
706 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
707
708 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
709 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
710 string locHistName, locHistTitle;
711
712 bool locIsRESTEvent = (string(locEventLoop->GetJEvent().GetJEventSource()->className()) == string("DEventSourceREST"));
713
714 //CREATE THE HISTOGRAMS
715 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
716 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
717 {
718 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
719 //If another thread has already created the folder, it just changes to it.
720 CreateAndChangeTo_ActionDirectory();
721
722 //Loop over filling for time-based and wire-based tracks
723 for(unsigned int locDummy = 0; locDummy < 2; ++locDummy)
724 {
725 bool locIsTimeBased = (locDummy == 0);
726 if(locIsRESTEvent && (!locIsTimeBased))
727 continue;
728
729 string locDirectoryName = locIsTimeBased ? "TimeBased" : "WireBased";
730 CreateAndChangeTo_Directory(locDirectoryName.c_str(), locDirectoryName.c_str());
731 string locTrackString = locIsTimeBased ? "Time-Based Tracks" : "Wire-Based Tracks";
732
733 //Kinematics of has (no) hit
734 vector<DetectorSystem_t> locDetectorSystems;
735 locDetectorSystems.push_back(SYS_START); locDetectorSystems.push_back(SYS_BCAL);
736 locDetectorSystems.push_back(SYS_TOF); locDetectorSystems.push_back(SYS_FCAL);
737 for(size_t loc_i = 0; loc_i < locDetectorSystems.size(); ++loc_i)
738 {
739 DetectorSystem_t locSystem = locDetectorSystems[loc_i];
740
741 double locMaxTheta = ((locSystem == SYS_FCAL) || (locSystem == SYS_TOF)) ? 12.0 : dMaxTheta;
742 double locMaxP = (locSystem == SYS_BCAL) ? 3.0 : dMaxP;
743
744 string locSystemName = SystemName(locSystem);
745 if(locSystemName == "ST")
746 locSystemName = "SC";
747 string locDirName = locSystemName;
748 if(locSystemName == "TOF")
749 locDirName = "TOFPoint";
750
751 CreateAndChangeTo_Directory(locDirName, locDirName);
752
753 // PVsTheta Has Hit
754 locHistName = "PVsTheta_HasHit";
755 locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" Has Hit;#theta#circ;p (GeV/c)");
756 dHistMap_PVsTheta_HasHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPBins, dMinP, locMaxP);
757
758 // PVsTheta Has No Hit
759 locHistName = "PVsTheta_NoHit";
760 locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" No Hit;#theta#circ;p (GeV/c)");
761 dHistMap_PVsTheta_NoHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPBins, dMinP, locMaxP);
762
763 // PhiVsTheta Has Hit
764 locHistName = "PhiVsTheta_HasHit";
765 locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" Has Hit;#theta#circ;#phi#circ");
766 dHistMap_PhiVsTheta_HasHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
767
768 // PhiVsTheta Has No Hit
769 locHistName = "PhiVsTheta_NoHit";
770 locHistTitle = locTrackString + string(", Has Other Match, ") + locSystemName + string(" No Hit;#theta#circ;#phi#circ");
771 dHistMap_PhiVsTheta_NoHit[locSystem][locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, locMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
772
773 gDirectory(TDirectory::CurrentDirectory())->cd("..");
774 }
775
776 //SC
777 CreateAndChangeTo_Directory("SC", "SC");
778 locHistName = "SCPaddleVsTheta_HasHit";
779 locHistTitle = locTrackString + string(", Has Other Match, SC Has Hit;#theta#circ;Projected SC Paddle");
780 dHistMap_SCPaddleVsTheta_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, 30, 0.5, 30.5);
781
782 locHistName = "SCPaddleVsTheta_NoHit";
783 locHistTitle = locTrackString + string(", Has Other Match, SC No Hit;#theta#circ;Projected SC Paddle");
784 dHistMap_SCPaddleVsTheta_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, 30, 0.5, 30.5);
785
786 locHistName = "SCPaddle_BarrelRegion_HasHit";
787 locHistTitle = locTrackString + string(", Has Other Match, SC Barrel Region Has Hit;Projected SC Paddle");
788 dHistMap_SCPaddle_BarrelRegion_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
789
790 locHistName = "SCPaddle_BarrelRegion_NoHit";
791 locHistTitle = locTrackString + string(", Has Other Match, SC Barrel Region No Hit;Projected SC Paddle");
792 dHistMap_SCPaddle_BarrelRegion_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
793
794 locHistName = "SCPaddle_NoseRegion_HasHit";
795 locHistTitle = locTrackString + string(", Has Other Match, SC Front Region Has Hit;Projected SC Paddle");
796 dHistMap_SCPaddle_NoseRegion_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
797
798 locHistName = "SCPaddle_NoseRegion_NoHit";
799 locHistTitle = locTrackString + string(", Has Other Match, SC Front Region No Hit;Projected SC Paddle");
800 dHistMap_SCPaddle_NoseRegion_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 30, 0.5, 30.5);
801
802 locHistName = "SCTrackDeltaPhiVsP";
803 locHistTitle = locTrackString + string(";p (GeV/c);SC / Track #Delta#phi#circ");
804 dHistMap_SCTrackDeltaPhiVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaPhiBins, dSCMatchMinDeltaPhi, dSCMatchMaxDeltaPhi);
805
806 locHistName = "SCTrackDeltaPhiVsTheta";
807 locHistTitle = locTrackString + string(";#theta#circ;SC / Track #Delta#phi#circ");
808 dHistMap_SCTrackDeltaPhiVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaPhiBins, dSCMatchMinDeltaPhi, dSCMatchMaxDeltaPhi);
809 gDirectory(TDirectory::CurrentDirectory())->cd("..");
810
811 //TOFPaddle
812 CreateAndChangeTo_Directory("TOFPaddle", "TOFPaddle");
813 locHistName = "VerticalPaddleTrackDeltaX";
814 locHistTitle = locTrackString + string(";Vertical TOF Paddle / Track |#DeltaX| (cm)");
815 dHistMap_TOFPaddleTrackDeltaX[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins/2, 0.0, dMaxTrackMatchDOCA);
816
817 locHistName = "HorizontalPaddleTrackDeltaY";
818 locHistTitle = locTrackString + string(";Horizontal TOF Paddle / Track |#DeltaY| (cm)");
819 dHistMap_TOFPaddleTrackDeltaY[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins/2, 0.0, dMaxTrackMatchDOCA);
820
821 locHistName = "TrackYVsVerticalPaddle_HasHit";
822 locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle Has Hit;Projected Vertical Paddle;Projected TOF Hit Y (cm)");
823 dHistMap_TOFPaddleTrackYVsVerticalPaddle_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNumFCALTOFXYBins, -130.0, 130.0);
824
825 locHistName = "TrackYVsVerticalPaddle_NoHit";
826 locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle No Hit;Projected Vertical Paddle;Projected TOF Hit Y (cm)");
827 dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNumFCALTOFXYBins, -130.0, 130.0);
828
829 locHistName = "HorizontalPaddleVsTrackX_HasHit";
830 locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle Has Hit;Projected TOF Hit X (cm);Projected Horizontal Paddle");
831 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, 44, 0.5, 44.5);
832
833 locHistName = "HorizontalPaddleVsTrackX_NoHit";
834 locHistTitle = locTrackString + string(", Has Other Match, TOF Paddle No Hit;Projected TOF Hit X (cm);Projected Horizontal Paddle");
835 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, 44, 0.5, 44.5);
836 gDirectory(TDirectory::CurrentDirectory())->cd("..");
837
838 //TOFPoint
839 CreateAndChangeTo_Directory("TOFPoint", "TOFPoint");
840 locHistName = "TrackTOFYVsX_HasHit";
841 locHistTitle = locTrackString + string(", Has Other Match, TOF Has Hit;Projected TOF Hit X (cm);Projected TOF Hit Y (cm)");
842 dHistMap_TrackTOFYVsX_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
843
844 locHistName = "TrackTOFYVsX_NoHit";
845 locHistTitle = locTrackString + string(", Has Other Match, TOF No Hit;Projected TOF Hit X (cm);Projected TOF Hit Y (cm)");
846 dHistMap_TrackTOFYVsX_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
847
848 locHistName = "TrackTOF2DPaddles_HasHit";
849 locHistTitle = locTrackString + string(", Has Other Match, TOF Has Hit;Projected Vertical TOF Paddle;Projected Horizontal TOF Paddle");
850 dHistMap_TrackTOF2DPaddles_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, 44, 0.5, 44.5);
851
852 locHistName = "TrackTOF2DPaddles_NoHit";
853 locHistTitle = locTrackString + string(", Has Other Match, TOF No Hit;Projected Vertical TOF Paddle;Projected Horizontal TOF Paddle");
854 dHistMap_TrackTOF2DPaddles_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, 44, 0.5, 44.5);
855
856 locHistName = "TOFTrackDistanceVsP";
857 locHistTitle = locTrackString + string(";p (GeV/c);TOF / Track Distance (cm)");
858 dHistMap_TOFPointTrackDistanceVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
859
860 locHistName = "TOFTrackDistanceVsTheta";
861 locHistTitle = locTrackString + string(";#theta#circ;TOF / Track Distance (cm)");
862 dHistMap_TOFPointTrackDistanceVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, 20.0, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
863
864 locHistName = "TOFTrackDeltaXVsHorizontalPaddle";
865 locHistTitle = locTrackString + string(";TOF Horizontal Paddle;TOF / Track #DeltaX (cm)");
866 dHistMap_TOFPointTrackDeltaXVsHorizontalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
867
868 locHistName = "TOFTrackDeltaXVsVerticalPaddle";
869 locHistTitle = locTrackString + string(";TOF Vertical Paddle;TOF / Track #DeltaX (cm)");
870 dHistMap_TOFPointTrackDeltaXVsVerticalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
871
872 locHistName = "TOFTrackDeltaYVsHorizontalPaddle";
873 locHistTitle = locTrackString + string(";TOF Horizontal Paddle;TOF / Track #DeltaY (cm)");
874 dHistMap_TOFPointTrackDeltaYVsHorizontalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
875
876 locHistName = "TOFTrackDeltaYVsVerticalPaddle";
877 locHistTitle = locTrackString + string(";TOF Vertical Paddle;TOF / Track #DeltaY (cm)");
878 dHistMap_TOFPointTrackDeltaYVsVerticalPaddle[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 44, 0.5, 44.5, dNum2DTrackDOCABins, -1.0*dMaxTrackMatchDOCA, dMaxTrackMatchDOCA);
879
880 locHistName = "TOFTrackDistance_BothPlanes";
881 locHistTitle = locTrackString + string("TOF Hit in Both Planes;TOF / Track Distance (cm)");
882 dHistMap_TOFPointTrackDistance_BothPlanes[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
883
884 locHistName = "TOFTrackDistance_OnePlane";
885 locHistTitle = locTrackString + string("TOF Hit in One Plane;TOF / Track Distance (cm)");
886 dHistMap_TOFPointTrackDistance_OnePlane[locIsTimeBased] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
887 gDirectory(TDirectory::CurrentDirectory())->cd("..");
888
889 //FCAL
890 CreateAndChangeTo_Directory("FCAL", "FCAL");
891 locHistName = "TrackFCALYVsX_HasHit";
892 locHistTitle = locTrackString + string(", Has Other Match, FCAL Has Hit;Projected FCAL Hit X (cm);Projected FCAL Hit Y (cm)");
893 dHistMap_TrackFCALYVsX_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
894
895 locHistName = "TrackFCALYVsX_NoHit";
896 locHistTitle = locTrackString + string(", Has Other Match, FCAL No Hit;Projected FCAL Hit X (cm);Projected FCAL Hit Y (cm)");
897 dHistMap_TrackFCALYVsX_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFCALTOFXYBins, -130.0, 130.0, dNumFCALTOFXYBins, -130.0, 130.0);
898
899 locHistName = "TrackFCALRowVsColumn_HasHit";
900 locHistTitle = locTrackString + string(", Has Other Match, FCAL Has Hit;Projected FCAL Hit Column;Projected FCAL Hit Row");
901 dHistMap_TrackFCALRowVsColumn_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 59, -0.5, 58.5, 59, -0.5, 58.5);
902
903 locHistName = "TrackFCALRowVsColumn_NoHit";
904 locHistTitle = locTrackString + string(", Has Other Match, FCAL No Hit;Projected FCAL Hit Column;Projected FCAL Hit Row");
905 dHistMap_TrackFCALRowVsColumn_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, 59, -0.5, 58.5, 59, -0.5, 58.5);
906
907 locHistName = "FCALTrackDistanceVsP";
908 locHistTitle = locTrackString + string(";p (GeV/c);FCAL / Track Distance (cm)");
909 dHistMap_FCALTrackDistanceVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
910
911 locHistName = "FCALTrackDistanceVsTheta";
912 locHistTitle = locTrackString + string(";#theta#circ;FCAL / Track Distance (cm)");
913 dHistMap_FCALTrackDistanceVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, 20.0, dNum2DTrackDOCABins, dMinTrackDOCA, dMaxTrackMatchDOCA);
914 gDirectory(TDirectory::CurrentDirectory())->cd("..");
915
916 //BCAL
917 CreateAndChangeTo_Directory("BCAL", "BCAL");
918 locHistName = "TrackBCALModuleVsZ_HasHit";
919 locHistTitle = locTrackString + string(", Has Other Match, BCAL Has Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit Module");
920 dHistMap_TrackBCALModuleVsZ_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, 48, 0.5, 48.5);
921
922 locHistName = "TrackBCALModuleVsZ_NoHit";
923 locHistTitle = locTrackString + string(", Has Other Match, BCAL No Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit Module");
924 dHistMap_TrackBCALModuleVsZ_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, 48, 0.5, 48.5);
925
926 locHistName = "TrackBCALPhiVsZ_HasHit";
927 locHistTitle = locTrackString + string(", Has Other Match, BCAL Has Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit #phi#circ");
928 dHistMap_TrackBCALPhiVsZ_HasHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, dNum2DPhiBins, dMinPhi, dMaxPhi);
929
930 locHistName = "TrackBCALPhiVsZ_NoHit";
931 locHistTitle = locTrackString + string(", Has Other Match, BCAL No Hit;Projected BCAL Hit Z (cm);Projected BCAL Hit #phi#circ");
932 dHistMap_TrackBCALPhiVsZ_NoHit[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALZBins, 0.0, 450.0, dNum2DPhiBins, dMinPhi, dMaxPhi);
933
934 locHistName = "BCALDeltaPhiVsP";
935 locHistTitle = locTrackString + string(";p (GeV/c);BCAL / Track #Delta#phi#circ");
936 dHistMap_BCALDeltaPhiVsP[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, 4.0, dNum2DDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
937
938 locHistName = "BCALDeltaZVsTheta";
939 locHistTitle = locTrackString + string(";#theta#circ;BCAL / Track #Deltaz (cm)");
940 dHistMap_BCALDeltaZVsTheta[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DDeltaZBins, dMinDeltaZ, dMaxDeltaZ);
941 gDirectory(TDirectory::CurrentDirectory())->cd("..");
942
943 //TRACKING
944 locHistName = "PVsTheta_NoHitMatch";
945 locHistTitle = locTrackString + string(", No Hit Match;#theta#circ;p (GeV/c)");
946 dHistMap_TrackPVsTheta_NoHitMatch[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
947
948 locHistName = "PVsTheta_HitMatch";
949 locHistTitle = locTrackString + string(", Hit Match;#theta#circ;p (GeV/c)");
950 dHistMap_TrackPVsTheta_HitMatch[locIsTimeBased] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
951
952 gDirectory(TDirectory::CurrentDirectory())->cd("..");
953 }
954
955 //Return to the base directory
956 ChangeTo_BaseDirectory();
957 }
958 japp->RootUnLock(); //RELEASE ROOT LOCK!!
959}
960
961bool DHistogramAction_DetectorMatching::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
962{
963 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
964
965 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
966 if(Get_NumPreviousParticleCombos() != 0)
967 return true;
968
969 bool locIsRESTEvent = (string(locEventLoop->GetJEvent().GetJEventSource()->className()) == string("DEventSourceREST"));
970
971 Fill_MatchingHists(locEventLoop, true); //Time-based tracks
972 if(!locIsRESTEvent)
973 Fill_MatchingHists(locEventLoop, false); //Wire-based tracks
974
975 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
976}
977
978void DHistogramAction_DetectorMatching::Fill_MatchingHists(JEventLoop* locEventLoop, bool locIsTimeBased)
979{
980 const DParticleID* locParticleID = NULL__null;
981 locEventLoop->GetSingle(locParticleID);
982
983 //can't make this a class member: may cause race condition
984 DCutAction_TrackHitPattern locCutAction_TrackHitPattern(NULL__null, dMinHitRingsPerCDCSuperlayer, dMinHitPlanesPerFDCPackage);
985 locCutAction_TrackHitPattern.Initialize(locEventLoop);
986
987 //get the best tracks for each candidate id, based on good hit pattern & tracking FOM
988 map<JObject::oid_t, const DKinematicData*> locBestTrackMap; //lowest tracking FOM for each candidate id
989 if(locIsTimeBased)
990 {
991 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
992 locEventLoop->Get(locTrackTimeBasedVector);
993
994 //select the best DTrackTimeBased for each track: of tracks with good hit pattern, use best tracking FOM
995 for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
996 {
997 if(locTrackTimeBasedVector[loc_i]->FOM < dMinTrackingFOM)
998 continue;
999 if(!locCutAction_TrackHitPattern.Cut_TrackHitPattern(locTrackTimeBasedVector[loc_i]))
1000 continue;
1001 JObject::oid_t locCandidateID = locTrackTimeBasedVector[loc_i]->candidateid;
1002 if(locBestTrackMap.find(locCandidateID) == locBestTrackMap.end())
1003 locBestTrackMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
1004 else if(locTrackTimeBasedVector[loc_i]->FOM > (dynamic_cast<const DTrackTimeBased*>(locBestTrackMap[locCandidateID]))->FOM)
1005 locBestTrackMap[locCandidateID] = locTrackTimeBasedVector[loc_i];
1006 }
1007 }
1008 else
1009 {
1010 vector<const DTrackWireBased*> locTrackWireBasedVector;
1011 locEventLoop->Get(locTrackWireBasedVector);
1012
1013 //select the best DTrackWireBased for each track: of tracks with good hit pattern, use best tracking FOM
1014 for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
1015 {
1016 if(locTrackWireBasedVector[loc_i]->FOM < dMinTrackingFOM)
1017 continue;
1018 if(!locCutAction_TrackHitPattern.Cut_TrackHitPattern(locTrackWireBasedVector[loc_i]))
1019 continue;
1020 JObject::oid_t locCandidateID = locTrackWireBasedVector[loc_i]->candidateid;
1021 if(locBestTrackMap.find(locCandidateID) == locBestTrackMap.end())
1022 locBestTrackMap[locCandidateID] = locTrackWireBasedVector[loc_i];
1023 else if(locTrackWireBasedVector[loc_i]->FOM > (dynamic_cast<const DTrackWireBased*>(locBestTrackMap[locCandidateID]))->FOM)
1024 locBestTrackMap[locCandidateID] = locTrackWireBasedVector[loc_i];
1025 }
1026 }
1027
1028 vector<const DBCALShower*> locBCALShowers;
1029 locEventLoop->Get(locBCALShowers);
1030
1031 vector<const DFCALShower*> locFCALShowers;
1032 locEventLoop->Get(locFCALShowers);
1033
1034 vector<const DTOFPoint*> locTOFPoints;
1035 locEventLoop->Get(locTOFPoints);
1036
1037 vector<const DTOFPaddleHit*> locTOFPaddleHits;
1038 locEventLoop->Get(locTOFPaddleHits);
1039
1040 vector<const DSCHit*> locSCHits;
1041 locEventLoop->Get(locSCHits);
1042
1043 string locDetectorMatchesTag = locIsTimeBased ? "" : "WireBased";
1044 const DDetectorMatches* locDetectorMatches = NULL__null;
1045 locEventLoop->GetSingle(locDetectorMatches, locDetectorMatchesTag.c_str());
1046
1047 map<JObject::oid_t, const DKinematicData*>::iterator locTrackIterator;
1048
1049 //TRACK / BCAL CLOSEST MATCHES
1050 map<const DKinematicData*, pair<double, double> > locBCALTrackDistanceMap; //first double is delta-phi, second delta-z
1051 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1052 {
1053 double locBestMatchDeltaPhi = 0.0, locBestMatchDeltaZ = 0.0;
1054 if(locParticleID->Get_ClosestToTrack_BCAL(locTrackIterator->second, locBCALShowers, locBestMatchDeltaPhi, locBestMatchDeltaZ) != NULL__null)
1055 locBCALTrackDistanceMap[locTrackIterator->second] = pair<double, double>(locBestMatchDeltaPhi, locBestMatchDeltaZ);
1056 }
1057
1058 //TRACK / FCAL CLOSEST MATCHES
1059 map<const DKinematicData*, double> locFCALTrackDistanceMap;
1060 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1061 {
1062 double locBestDistance = 999.0;
1063 if(locParticleID->Get_ClosestToTrack_FCAL(locTrackIterator->second, locFCALShowers, locBestDistance) != NULL__null)
1064 locFCALTrackDistanceMap[locTrackIterator->second] = locBestDistance;
1065 }
1066
1067 //TRACK / TOF POINT CLOSEST MATCHES
1068 map<const DKinematicData*, pair<const DTOFPoint*, pair<double, double> > > locTOFPointTrackDistanceMap; //doubles: delta-x, delta-y
1069 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1070 {
1071 double locBestDeltaX, locBestDeltaY;
1072 const DTOFPoint* locClosestTOFPoint = locParticleID->Get_ClosestToTrack_TOFPoint(locTrackIterator->second, locTOFPoints, locBestDeltaX, locBestDeltaY);
1073 if(locClosestTOFPoint == NULL__null)
1074 continue;
1075 pair<double, double> locDeltas(locBestDeltaX, locBestDeltaY);
1076 locTOFPointTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPoint*, pair<double, double> >(locClosestTOFPoint, locDeltas);
1077 }
1078
1079 //TRACK / TOF PADDLE CLOSEST MATCHES
1080 map<const DKinematicData*, pair<const DTOFPaddleHit*, double> > locHorizontalTOFPaddleTrackDistanceMap; //double: delta-y
1081 map<const DKinematicData*, pair<const DTOFPaddleHit*, double> > locVerticalTOFPaddleTrackDistanceMap; //double: delta-x
1082 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1083 {
1084 double locBestDeltaX = 999.9, locBestDeltaY = 999.9;
1085 pair<const DTOFPaddleHit*, const DTOFPaddleHit*> locClosestTOFPaddleHits = locParticleID->Get_ClosestToTrack_TOFPaddles(locTrackIterator->second, locTOFPaddleHits, locBestDeltaX, locBestDeltaY);
1086 if(locClosestTOFPaddleHits.first != NULL__null)
1087 locVerticalTOFPaddleTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPaddleHit*, double>(locClosestTOFPaddleHits.first, locBestDeltaX);
1088 if(locClosestTOFPaddleHits.second != NULL__null)
1089 locHorizontalTOFPaddleTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPaddleHit*, double>(locClosestTOFPaddleHits.first, locBestDeltaY);
1090 }
1091
1092 //TRACK / SC CLOSEST MATCHES
1093 map<const DKinematicData*, double> locSCTrackDistanceMap;
1094 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1095 {
1096 double locBestDeltaPhi = 999.0;
1097 if(locParticleID->Get_ClosestToTrack_SC(locTrackIterator->second, locSCHits, locBestDeltaPhi) != NULL__null)
1098 locSCTrackDistanceMap[locTrackIterator->second] = locBestDeltaPhi;
1099 }
1100
1101 //PROJECTED HIT POSITIONS
1102 map<const DKinematicData*, pair<int, bool> > dProjectedSCPaddleMap; //pair: paddle, hit-barrel-flag (false if bend/nose)
1103 map<const DKinematicData*, pair<int, int> > dProjectedTOF2DPaddlesMap; //pair: vertical, horizontal
1104 map<const DKinematicData*, pair<float, float> > dProjectedTOFXYMap; //pair: x, y
1105 map<const DKinematicData*, pair<int, int> > dProjectedFCALRowColumnMap; //pair: column, row
1106 map<const DKinematicData*, pair<float, float> > dProjectedFCALXYMap; //pair: x, y
1107 map<const DKinematicData*, pair<float, int> > dProjectedBCALModuleSectorMap; //pair: z, module
1108 map<const DKinematicData*, pair<float, float> > dProjectedBCALPhiZMap; //pair: z, phi
1109 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1110 {
1111 const DKinematicData* locTrack = locTrackIterator->second;
1112 const DReferenceTrajectory* locReferenceTrajectory = NULL__null;
1113 if(locIsTimeBased)
1114 locReferenceTrajectory = (static_cast<const DTrackTimeBased*>(locTrack))->rt;
1115 else
1116 locReferenceTrajectory = (static_cast<const DTrackWireBased*>(locTrack))->rt;
1117 if(locReferenceTrajectory == NULL__null)
1118 break; //e.g. REST data: no trajectory
1119
1120 //SC
1121 DVector3 locSCIntersection;
1122 bool locProjBarrelFlag = false;
1123 unsigned int locProjectedSCPaddle = locParticleID->PredictSCSector(locReferenceTrajectory, 999.9, &locSCIntersection, &locProjBarrelFlag);
1124 if(locProjectedSCPaddle != 0)
1125 dProjectedSCPaddleMap[locTrack] = pair<int, bool>(locProjectedSCPaddle, locProjBarrelFlag);
1126
1127 //TOF
1128 DVector3 locTOFIntersection;
1129 unsigned int locHorizontalBar = 0, locVerticalBar = 0;
1130 if(locParticleID->PredictTOFPaddles(locReferenceTrajectory, locHorizontalBar, locVerticalBar, &locTOFIntersection))
1131 {
1132 dProjectedTOF2DPaddlesMap[locTrack] = pair<int, int>(locVerticalBar, locHorizontalBar);
1133 dProjectedTOFXYMap[locTrack] = pair<float, float>(locTOFIntersection.X(), locTOFIntersection.Y());
1134 }
1135
1136 //FCAL
1137 DVector3 locFCALIntersection;
1138 unsigned int locRow = 0, locColumn = 0;
1139 if(locParticleID->PredictFCALHit(locReferenceTrajectory, locRow, locColumn, &locFCALIntersection))
1140 {
1141 dProjectedFCALRowColumnMap[locTrack] = pair<int, int>(locColumn, locRow);
1142 dProjectedFCALXYMap[locTrack] = pair<float, float>(locFCALIntersection.X(), locFCALIntersection.Y());
1143 }
1144
1145 //BCAL
1146 DVector3 locBCALIntersection;
1147 unsigned int locModule = 0, locSector = 0;
1148 if(locParticleID->PredictBCALWedge(locReferenceTrajectory, locModule, locSector, &locBCALIntersection))
1149 {
1150 dProjectedBCALModuleSectorMap[locTrack] = pair<float, int>(locBCALIntersection.Z(), locModule);
1151 dProjectedBCALPhiZMap[locTrack] = pair<float, float>(locBCALIntersection.Z(), locBCALIntersection.Phi()*180.0/TMath::Pi());
1152 }
1153 }
1154
1155 //FILL HISTOGRAMS
1156 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1157 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1158 Lock_Action(); //ACQUIRE ROOT LOCK!!
1159 {
1160 /********************************************************** MATCHING DISTANCE **********************************************************/
1161
1162 //BCAL
1163 map<const DKinematicData*, pair<double, double> >::iterator locBCALIterator = locBCALTrackDistanceMap.begin();
1164 for(; locBCALIterator != locBCALTrackDistanceMap.end(); ++locBCALIterator)
1165 {
1166 const DKinematicData* locTrack = locBCALIterator->first;
1167 dHistMap_BCALDeltaPhiVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locBCALIterator->second.first*180.0/TMath::Pi());
1168 dHistMap_BCALDeltaZVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locBCALIterator->second.second);
1169 }
1170
1171 //FCAL
1172 map<const DKinematicData*, double>::iterator locFCALIterator = locFCALTrackDistanceMap.begin();
1173 for(; locFCALIterator != locFCALTrackDistanceMap.end(); ++locFCALIterator)
1174 {
1175 const DKinematicData* locTrack = locFCALIterator->first;
1176 dHistMap_FCALTrackDistanceVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locFCALIterator->second);
1177 dHistMap_FCALTrackDistanceVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locFCALIterator->second);
1178 }
1179
1180 //TOF Paddle
1181 //Horizontal
1182 map<const DKinematicData*, pair<const DTOFPaddleHit*, double> >::iterator locTOFPaddleIterator = locHorizontalTOFPaddleTrackDistanceMap.begin();
1183 for(; locTOFPaddleIterator != locHorizontalTOFPaddleTrackDistanceMap.end(); ++locTOFPaddleIterator)
1184 {
1185 double locDistance = locTOFPaddleIterator->second.second;
1186 dHistMap_TOFPaddleTrackDeltaY[locIsTimeBased]->Fill(locDistance);
1187 }
1188 //Vertical
1189 locTOFPaddleIterator = locVerticalTOFPaddleTrackDistanceMap.begin();
1190 for(; locTOFPaddleIterator != locVerticalTOFPaddleTrackDistanceMap.end(); ++locTOFPaddleIterator)
1191 {
1192 double locDistance = locTOFPaddleIterator->second.second;
1193 dHistMap_TOFPaddleTrackDeltaX[locIsTimeBased]->Fill(locDistance);
1194 }
1195
1196 //TOF Point
1197 map<const DKinematicData*, pair<const DTOFPoint*, pair<double, double> > >::iterator locTOFPointIterator = locTOFPointTrackDistanceMap.begin();
1198 for(; locTOFPointIterator != locTOFPointTrackDistanceMap.end(); ++locTOFPointIterator)
1199 {
1200 const DKinematicData* locTrack = locTOFPointIterator->first;
1201 const DTOFPoint* locTOFPoint = locTOFPointIterator->second.first;
1202 double locDeltaX = locTOFPointIterator->second.second.first;
1203 double locDeltaY = locTOFPointIterator->second.second.second;
1204
1205 double locDistance = sqrt(locDeltaX*locDeltaX + locDeltaY*locDeltaY);
1206 if((locDeltaX < 500.0) && (locDeltaY < 500.0)) //else position not well-defined
1207 {
1208 dHistMap_TOFPointTrackDistanceVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDistance);
1209 dHistMap_TOFPointTrackDistanceVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locDistance);
1210 if((locTOFPoint->dHorizontalBar != 0) && (locTOFPoint->dVerticalBar != 0))
1211 dHistMap_TOFPointTrackDistance_BothPlanes[locIsTimeBased]->Fill(locDistance);
1212 else
1213 dHistMap_TOFPointTrackDistance_OnePlane[locIsTimeBased]->Fill(locDistance);
1214 }
1215
1216 dHistMap_TOFPointTrackDeltaXVsHorizontalPaddle[locIsTimeBased]->Fill(locTOFPoint->dHorizontalBar, locDeltaX);
1217 dHistMap_TOFPointTrackDeltaXVsVerticalPaddle[locIsTimeBased]->Fill(locTOFPoint->dVerticalBar, locDeltaX);
1218
1219 dHistMap_TOFPointTrackDeltaYVsHorizontalPaddle[locIsTimeBased]->Fill(locTOFPoint->dHorizontalBar, locDeltaY);
1220 dHistMap_TOFPointTrackDeltaYVsVerticalPaddle[locIsTimeBased]->Fill(locTOFPoint->dVerticalBar, locDeltaY);
1221 }
1222
1223 //SC
1224 map<const DKinematicData*, double>::iterator locSCIterator = locSCTrackDistanceMap.begin();
1225 if(locSCHits.size() <= 4) //don't fill if every paddle fired!
1226 {
1227 for(; locSCIterator != locSCTrackDistanceMap.end(); ++locSCIterator)
1228 {
1229 const DKinematicData* locTrack = locSCIterator->first;
1230 double locDeltaPhi = locSCIterator->second*180.0/TMath::Pi();
1231 dHistMap_SCTrackDeltaPhiVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDeltaPhi);
1232 dHistMap_SCTrackDeltaPhiVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locDeltaPhi);
1233 }
1234 }
1235
1236 /********************************************************* MATCHING EFFICINECY *********************************************************/
1237
1238 //Does-it-match, by detector
1239 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1240 {
1241 const DKinematicData* locTrack = locTrackIterator->second;
1242 double locTheta = locTrack->momentum().Theta()*180.0/TMath::Pi();
1243 double locPhi = locTrack->momentum().Phi()*180.0/TMath::Pi();
1244 double locP = locTrack->momentum().Mag();
1245
1246 //BCAL
1247 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_START))
1248 {
1249 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_BCAL))
1250 {
1251 dHistMap_PVsTheta_HasHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locP);
1252 dHistMap_PhiVsTheta_HasHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1253 if(dProjectedBCALModuleSectorMap.find(locTrack) != dProjectedBCALModuleSectorMap.end())
1254 {
1255 pair<float, float>& locPositionPair = dProjectedBCALPhiZMap[locTrack];
1256 dHistMap_TrackBCALPhiVsZ_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1257 pair<float, int>& locElementPair = dProjectedBCALModuleSectorMap[locTrack];
1258 dHistMap_TrackBCALModuleVsZ_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1259 }
1260 }
1261 else
1262 {
1263 dHistMap_PVsTheta_NoHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locP);
1264 dHistMap_PhiVsTheta_NoHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1265 if(dProjectedBCALModuleSectorMap.find(locTrack) != dProjectedBCALModuleSectorMap.end())
1266 {
1267 pair<float, float>& locPositionPair = dProjectedBCALPhiZMap[locTrack];
1268 dHistMap_TrackBCALPhiVsZ_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1269 pair<float, int>& locElementPair = dProjectedBCALModuleSectorMap[locTrack];
1270 dHistMap_TrackBCALModuleVsZ_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1271 }
1272 }
1273 }
1274
1275 //FCAL
1276 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1277 {
1278 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1279 {
1280 dHistMap_PVsTheta_HasHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locP);
1281 dHistMap_PhiVsTheta_HasHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1282 if(dProjectedFCALRowColumnMap.find(locTrack) != dProjectedFCALRowColumnMap.end())
1283 {
1284 pair<float, float>& locPositionPair = dProjectedFCALXYMap[locTrack];
1285 dHistMap_TrackFCALYVsX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1286 pair<int, int>& locElementPair = dProjectedFCALRowColumnMap[locTrack];
1287 dHistMap_TrackFCALRowVsColumn_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1288 }
1289 }
1290 else
1291 {
1292 dHistMap_PVsTheta_NoHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locP);
1293 dHistMap_PhiVsTheta_NoHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1294 if(dProjectedFCALRowColumnMap.find(locTrack) != dProjectedFCALRowColumnMap.end())
1295 {
1296 pair<float, float>& locPositionPair = dProjectedFCALXYMap[locTrack];
1297 dHistMap_TrackFCALYVsX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1298 pair<int, int>& locElementPair = dProjectedFCALRowColumnMap[locTrack];
1299 dHistMap_TrackFCALRowVsColumn_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1300 }
1301 }
1302 }
1303
1304 //TOF Paddle
1305 if((dProjectedTOFXYMap.find(locTrack) != dProjectedTOFXYMap.end()) && locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1306 {
1307 pair<float, float>& locPositionPair = dProjectedTOFXYMap[locTrack];
1308 pair<int, int>& locPaddlePair = dProjectedTOF2DPaddlesMap[locTrack]; //vertical, horizontal
1309
1310 //Horizontal
1311 if(locHorizontalTOFPaddleTrackDistanceMap.find(locTrack) != locHorizontalTOFPaddleTrackDistanceMap.end())
1312 {
1313 double locDistance = locHorizontalTOFPaddleTrackDistanceMap[locTrack].second;
1314 if(locDistance <= dMinTOFPaddleMatchDistance) //match
1315 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1316 else //no match
1317 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1318 }
1319 else // no match
1320 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1321
1322 //Vertical
1323 if(locVerticalTOFPaddleTrackDistanceMap.find(locTrack) != locVerticalTOFPaddleTrackDistanceMap.end())
1324 {
1325 double locDistance = locVerticalTOFPaddleTrackDistanceMap[locTrack].second;
1326 if(locDistance <= dMinTOFPaddleMatchDistance) //match
1327 dHistMap_TOFPaddleTrackYVsVerticalPaddle_HasHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1328 else //no match
1329 dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1330 }
1331 else // no match
1332 dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1333 }
1334
1335 //TOF Point
1336 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1337 {
1338 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1339 {
1340 dHistMap_PVsTheta_HasHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locP);
1341 dHistMap_PhiVsTheta_HasHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locPhi);
1342 if(dProjectedTOFXYMap.find(locTrack) != dProjectedTOFXYMap.end())
1343 {
1344 pair<float, float>& locPositionPair = dProjectedTOFXYMap[locTrack];
1345 dHistMap_TrackTOFYVsX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1346 pair<int, int>& locElementPair = dProjectedTOF2DPaddlesMap[locTrack];
1347 dHistMap_TrackTOF2DPaddles_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1348 }
1349 }
1350 else
1351 {
1352 dHistMap_PVsTheta_NoHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locP);
1353 dHistMap_PhiVsTheta_NoHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locPhi);
1354 if(dProjectedTOFXYMap.find(locTrack) != dProjectedTOFXYMap.end())
1355 {
1356 pair<float, float>& locPositionPair = dProjectedTOFXYMap[locTrack];
1357 dHistMap_TrackTOFYVsX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1358 pair<int, int>& locElementPair = dProjectedTOF2DPaddlesMap[locTrack];
1359 dHistMap_TrackTOF2DPaddles_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1360 }
1361 }
1362 }
1363
1364 //SC
1365 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_BCAL) || locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL) || locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1366 {
1367 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_START))
1368 {
1369 dHistMap_PVsTheta_HasHit[SYS_START][locIsTimeBased]->Fill(locTheta, locP);
1370 dHistMap_PhiVsTheta_HasHit[SYS_START][locIsTimeBased]->Fill(locTheta, locPhi);
1371 if(dProjectedSCPaddleMap.find(locTrack) != dProjectedSCPaddleMap.end())
1372 {
1373 dHistMap_SCPaddleVsTheta_HasHit[locIsTimeBased]->Fill(locTheta, dProjectedSCPaddleMap[locTrack].first);
1374 if(dProjectedSCPaddleMap[locTrack].second)
1375 dHistMap_SCPaddle_BarrelRegion_HasHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1376 else
1377 dHistMap_SCPaddle_NoseRegion_HasHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1378 }
1379 }
1380 else
1381 {
1382 dHistMap_PVsTheta_NoHit[SYS_START][locIsTimeBased]->Fill(locTheta, locP);
1383 dHistMap_PhiVsTheta_NoHit[SYS_START][locIsTimeBased]->Fill(locTheta, locPhi);
1384 if(dProjectedSCPaddleMap.find(locTrack) != dProjectedSCPaddleMap.end())
1385 {
1386 dHistMap_SCPaddleVsTheta_NoHit[locIsTimeBased]->Fill(locTheta, dProjectedSCPaddleMap[locTrack].first);
1387 if(dProjectedSCPaddleMap[locTrack].second)
1388 dHistMap_SCPaddle_BarrelRegion_NoHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1389 else
1390 dHistMap_SCPaddle_NoseRegion_NoHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1391 }
1392 }
1393 }
1394 }
1395 //Is-Matched to Something
1396 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1397 {
1398 const DKinematicData* locTrack = locTrackIterator->second;
1399 double locTheta = locTrack->momentum().Theta()*180.0/TMath::Pi();
1400 double locP = locTrack->momentum().Mag();
1401 if(locDetectorMatches->Get_IsMatchedToHit(locTrack))
1402 dHistMap_TrackPVsTheta_HitMatch[locIsTimeBased]->Fill(locTheta, locP);
1403 else
1404 dHistMap_TrackPVsTheta_NoHitMatch[locIsTimeBased]->Fill(locTheta, locP);
1405 }
1406 }
1407 Unlock_Action(); //RELEASE ROOT LOCK!!
1408}
1409
1410void DHistogramAction_DetectorPID::Initialize(JEventLoop* locEventLoop)
1411{
1412 //Create any histograms/trees/etc. within a ROOT lock.
1413 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
1414
1415 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
1416 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
1417
1418 string locHistName, locHistTitle, locParticleROOTName;
1419
1420 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
1421 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
1422 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
1423 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
1424 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
1425
1426 //CREATE THE HISTOGRAMS
1427 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1428 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1429 {
1430 //So: Default tag is "", User can set it to something else
1431 //In here, if tag is "", get from gparms, if not, leave it alone
1432 //If gparms value does not exist, set it to (and use) "PreSelect"
1433 if(dTrackSelectionTag == "NotATag")
1434 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
1435 if(dShowerSelectionTag == "NotATag")
1436 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
1437
1438 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
1439 //If another thread has already created the folder, it just changes to it.
1440 CreateAndChangeTo_ActionDirectory();
1441
1442 //q = 0
1443 locParticleROOTName = ParticleName_ROOT(Gamma);
1444 //BCAL
1445 CreateAndChangeTo_Directory("BCAL", "BCAL");
1446
1447 locHistName = "BetaVsP_q0";
1448 locHistTitle = "BCAL q^{0};Shower Energy (GeV);#beta";
1449 dHistMap_BetaVsP[SYS_BCAL][0] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1450
1451 locHistName = "DeltaTVsShowerE_Photon";
1452 locHistTitle = string("BCAL q^{0}") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat_{BCAL - RF}");
1453 dHistMap_DeltaTVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1454
1455/*
1456 //Uncomment when ready!
1457 locHistName = "TimePullVsShowerE_Photon";
1458 locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
1459 dHistMap_TimePullVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1460
1461 locHistName = "TimeFOMVsShowerE_Photon";
1462 locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
1463 dHistMap_TimeFOMVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1464*/
1465 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1466
1467 //FCAL
1468 CreateAndChangeTo_Directory("FCAL", "FCAL");
1469
1470 locHistName = "BetaVsP_q0";
1471 locHistTitle = "FCAL q^{0};Shower Energy (GeV);#beta";
1472 dHistMap_BetaVsP[SYS_FCAL][0] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1473
1474 locHistName = "DeltaTVsShowerE_Photon";
1475 locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat_{FCAL - RF}");
1476 dHistMap_DeltaTVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1477
1478/*
1479 //Uncomment when ready!
1480 locHistName = "TimePullVsShowerE_Photon";
1481 locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
1482 dHistMap_TimePullVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1483
1484 locHistName = "TimeFOMVsShowerE_Photon";
1485 locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
1486 dHistMap_TimeFOMVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1487*/
1488 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1489
1490 //q +/-
1491 for(int locCharge = -1; locCharge <= 1; locCharge += 2)
1492 {
1493 string locParticleName = (locCharge == -1) ? "q-" : "q+";
1494 string locParticleROOTName = (locCharge == -1) ? "q^{-}" : "q^{+}";
1495
1496 //SC
1497 CreateAndChangeTo_Directory("SC", "SC");
1498
1499 locHistName = string("dEdXVsP_") + locParticleName;
1500 locHistTitle = locParticleROOTName + ";p (GeV/c);SC dE/dX (MeV/cm)";
1501 dHistMap_dEdXVsP[SYS_START][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1502
1503 locHistName = string("BetaVsP_") + locParticleName;
1504 locHistTitle = string("SC ") + locParticleROOTName + string(";p (GeV/c);#beta");
1505 dHistMap_BetaVsP[SYS_START][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1506
1507 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1508
1509 //TOF
1510 CreateAndChangeTo_Directory("TOF", "TOF");
1511
1512 locHistName = string("dEdXVsP_") + locParticleName;
1513 locHistTitle = locParticleROOTName + ";p (GeV/c);TOF dE/dX (MeV/cm)";
1514 dHistMap_dEdXVsP[SYS_TOF][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1515
1516 locHistName = string("BetaVsP_") + locParticleName;
1517 locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);#beta");
1518 dHistMap_BetaVsP[SYS_TOF][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1519
1520 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1521
1522 //BCAL
1523 CreateAndChangeTo_Directory("BCAL", "BCAL");
1524
1525 locHistName = string("BetaVsP_") + locParticleName;
1526 locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);#beta");
1527 dHistMap_BetaVsP[SYS_BCAL][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1528
1529 locHistName = string("EOverPVsP_") + locParticleName;
1530 locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);E_{Shower}/p_{Track} (c);");
1531 dHistMap_BCALEOverPVsP[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1532
1533 locHistName = string("EOverPVsTheta_") + locParticleName;
1534 locHistTitle = string("BCAL ") + locParticleROOTName + string(";#theta#circ;E_{Shower}/p_{Track} (c);");
1535 dHistMap_BCALEOverPVsTheta[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALThetaBins, dMinBCALTheta, dMaxBCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1536
1537 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1538
1539 //FCAL
1540 CreateAndChangeTo_Directory("FCAL", "FCAL");
1541
1542 locHistName = string("BetaVsP_") + locParticleName;
1543 locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#beta");
1544 dHistMap_BetaVsP[SYS_FCAL][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1545
1546 locHistName = string("EOverPVsP_") + locParticleName;
1547 locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);E_{Shower}/p_{Track} (c);");
1548 dHistMap_FCALEOverPVsP[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1549
1550 locHistName = string("EOverPVsTheta_") + locParticleName;
1551 locHistTitle = string("FCAL ") + locParticleROOTName + string(";#theta#circ;E_{Shower}/p_{Track} (c);");
1552 dHistMap_FCALEOverPVsTheta[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DFCALThetaBins, dMinFCALTheta, dMaxFCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1553
1554 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1555
1556 //CDC
1557 CreateAndChangeTo_Directory("CDC", "CDC");
1558
1559 locHistName = string("dEdXVsP_") + locParticleName;
1560 locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx (keV/cm)");
1561 dHistMap_dEdXVsP[SYS_CDC][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1562
1563 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1564
1565 //FDC
1566 CreateAndChangeTo_Directory("FDC", "FDC");
1567
1568 locHistName = string("dEdXVsP_") + locParticleName;
1569 locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC dE/dx (keV/cm)");
1570 dHistMap_dEdXVsP[SYS_FDC][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1571
1572 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1573 }
1574
1575 //delta's by PID
1576 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
1577 {
1578 Particle_t locPID = dFinalStatePIDs[loc_i];
1579 string locParticleName = ParticleType(locPID);
1580 string locParticleROOTName = ParticleName_ROOT(locPID);
1581
1582 //SC
1583 CreateAndChangeTo_Directory("SC", "SC");
1584
1585 locHistName = string("DeltadEdXVsP_") + locParticleName;
1586 locHistTitle = locParticleROOTName + " Candidates;p (GeV/c);SC #Delta(dE/dX) (MeV/cm)";
1587 dHistMap_DeltadEdXVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1588
1589 locHistName = string("DeltaBetaVsP_") + locParticleName;
1590 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1591 dHistMap_DeltaBetaVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1592
1593 locHistName = string("DeltaTVsP_") + locParticleName;
1594 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{SC - RF}");
1595 dHistMap_DeltaTVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1596
1597/*
1598 //Uncomment when ready!
1599 locHistName = string("TimePullVsP_") + locParticleName;
1600 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1601 dHistMap_TimePullVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1602
1603 locHistName = string("TimeFOMVsP_") + locParticleName;
1604 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1605 dHistMap_TimeFOMVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1606*/
1607 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1608
1609 //TOF
1610 CreateAndChangeTo_Directory("TOF", "TOF");
1611
1612 locHistName = string("DeltadEdXVsP_") + locParticleName;
1613 locHistTitle = locParticleROOTName + " Candidates;p (GeV/c);TOF #Delta(dE/dX) (MeV/cm)";
1614 dHistMap_DeltadEdXVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1615
1616 locHistName = string("DeltaBetaVsP_") + locParticleName;
1617 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1618 dHistMap_DeltaBetaVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1619
1620 locHistName = string("DeltaTVsP_") + locParticleName;
1621 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{TOF - RF}");
1622 dHistMap_DeltaTVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1623
1624/*
1625 //Uncomment when ready!
1626 locHistName = string("TimePullVsP_") + locParticleName;
1627 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1628 dHistMap_TimePullVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1629
1630 locHistName = string("TimeFOMVsP_") + locParticleName;
1631 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1632 dHistMap_TimeFOMVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1633*/
1634 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1635
1636 //BCAL
1637 CreateAndChangeTo_Directory("BCAL", "BCAL");
1638
1639 locHistName = string("DeltaBetaVsP_") + locParticleName;
1640 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1641 dHistMap_DeltaBetaVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1642
1643 locHistName = string("DeltaTVsP_") + locParticleName;
1644 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{BCAL - RF}");
1645 dHistMap_DeltaTVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1646
1647/*
1648 //Uncomment when ready!
1649 locHistName = string("TimePullVsP_") + locParticleName;
1650 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1651 dHistMap_TimePullVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1652
1653 locHistName = string("TimeFOMVsP_") + locParticleName;
1654 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1655 dHistMap_TimeFOMVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1656*/
1657 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1658
1659 //FCAL
1660 CreateAndChangeTo_Directory("FCAL", "FCAL");
1661
1662 locHistName = string("DeltaBetaVsP_") + locParticleName;
1663 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1664 dHistMap_DeltaBetaVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1665
1666 locHistName = string("DeltaTVsP_") + locParticleName;
1667 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{FCAL - RF}");
1668 dHistMap_DeltaTVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1669
1670/*
1671 //Uncomment when ready!
1672 locHistName = string("TimePullVsP_") + locParticleName;
1673 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1674 dHistMap_TimePullVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1675
1676 locHistName = string("TimeFOMVsP_") + locParticleName;
1677 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1678 dHistMap_TimeFOMVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1679*/
1680 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1681
1682 //CDC
1683 CreateAndChangeTo_Directory("CDC", "CDC");
1684
1685 locHistName = string("DeltadEdXVsP_") + locParticleName;
1686 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) (keV/cm)");
1687 dHistMap_DeltadEdXVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1688
1689/*
1690 //Uncomment when ready!
1691 locHistName = string("dEdXPullVsP_") + locParticleName;
1692 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
1693 dHistMap_dEdXPullVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1694
1695 locHistName = string("dEdXFOMVsP_") + locParticleName;
1696 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx PID Confidence Level");
1697 dHistMap_dEdXFOMVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1698*/
1699 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1700
1701 //FDC
1702 CreateAndChangeTo_Directory("FDC", "FDC");
1703
1704 locHistName = string("DeltadEdXVsP_") + locParticleName;
1705 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX) (keV/cm)");
1706 dHistMap_DeltadEdXVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1707
1708/*
1709 //Uncomment when ready!
1710 locHistName = string("dEdXPullVsP_") + locParticleName;
1711 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
1712 dHistMap_dEdXPullVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1713
1714 locHistName = string("dEdXFOMVsP_") + locParticleName;
1715 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC dE/dx PID Confidence Level");
1716 dHistMap_dEdXFOMVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1717*/
1718 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1719 }
1720
1721 //Return to the base directory
1722 ChangeTo_BaseDirectory();
1723 }
1724 japp->RootUnLock(); //RELEASE ROOT LOCK!!
1725}
1726
1727bool DHistogramAction_DetectorPID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1728{
1729 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
1730
1731 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
1732 if(Get_NumPreviousParticleCombos() != 0)
1733 return true;
1734
1735 vector<const DChargedTrack*> locChargedTracks;
1736 locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
1737
1738 vector<const DNeutralParticle*> locNeutralParticles;
1739 locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
1740
1741 const DDetectorMatches* locDetectorMatches = NULL__null;
1742 locEventLoop->GetSingle(locDetectorMatches);
1743
1744 const DEventRFBunch* locEventRFBunch = NULL__null;
1745 locEventLoop->GetSingle(locEventRFBunch);
1746
1747 const DParticleID* locParticleID = NULL__null;
1748 locEventLoop->GetSingle(locParticleID);
1749
1750 //FILL HISTOGRAMS
1751 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1752 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1753 Lock_Action(); //ACQUIRE ROOT LOCK!!
1754 {
1755 if(locEventRFBunch->dTimeSource != SYS_NULL) //only histogram beta for neutrals if the t0 is well known
1756 {
1757 for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
1758 {
1759 //doesn't matter which hypothesis you use for beta: t0 is from DEventVertex time
1760 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
1761 double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
1762 const DNeutralShower* locNeutralShower = locNeutralParticles[loc_i]->dNeutralShower;
1763 double locShowerEnergy = locNeutralShower->dEnergy;
1764
1765 double locDeltaT = locNeutralParticleHypothesis->time() - locEventRFBunch->dTime;
1766 if(locNeutralShower->dDetectorSystem == SYS_BCAL)
1767 {
1768 dHistMap_BetaVsP[SYS_BCAL][0]->Fill(locShowerEnergy, locBeta_Timing);
1769 dHistMap_DeltaTVsP[SYS_BCAL][Gamma]->Fill(locShowerEnergy, locDeltaT);
1770 }
1771 else
1772 {
1773 dHistMap_BetaVsP[SYS_FCAL][0]->Fill(locShowerEnergy, locBeta_Timing);
1774 dHistMap_DeltaTVsP[SYS_FCAL][Gamma]->Fill(locShowerEnergy, locDeltaT);
1775 }
1776 }
1777 }
1778
1779 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
1780 {
1781 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestTrackingFOM();
1782 int locCharge = ParticleCharge(locChargedTrackHypothesis->PID());
1783 if(dHistMap_dEdXVsP[SYS_START].find(locCharge) == dHistMap_dEdXVsP[SYS_START].end())
1784 continue;
1785
1786 double locStartTime = locParticleID->Calc_PropagatedRFTime(locChargedTrackHypothesis, locEventRFBunch);
1787
1788 const DTrackTimeBased* locTrackTimeBased = NULL__null;
1789 locChargedTrackHypothesis->GetSingle(locTrackTimeBased);
1790
1791 Particle_t locPID = locChargedTrackHypothesis->PID();
1792 double locP = locTrackTimeBased->momentum().Mag();
1793 double locTheta = locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi();
1794
1795 //if RF time is indeterminate, start time will be NaN
1796 const DBCALShowerMatchParams* locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
1797 const DFCALShowerMatchParams* locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
1798 const DTOFHitMatchParams* locTOFHitMatchParams = locChargedTrackHypothesis->Get_TOFHitMatchParams();
1799 const DSCHitMatchParams* locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
1800
1801 if(locSCHitMatchParams != NULL__null)
1802 {
1803 dHistMap_dEdXVsP[SYS_START][locCharge]->Fill(locP, locSCHitMatchParams->dEdx*1.0E3);
1804 if((locEventRFBunch->dTimeSource != SYS_START) && (locEventRFBunch->dNumParticleVotes > 1))
1805 {
1806 //If SC was used for RF time, don't compute delta-beta
1807 double locBeta_Timing = locSCHitMatchParams->dPathLength/(29.9792458*(locSCHitMatchParams->dHitTime - locChargedTrackHypothesis->t0()));
1808 dHistMap_BetaVsP[SYS_START][locCharge]->Fill(locP, locBeta_Timing);
1809 if(dHistMap_DeltaBetaVsP[SYS_START].find(locPID) != dHistMap_DeltaBetaVsP[SYS_START].end())
1810 {
1811 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1812 dHistMap_DeltaBetaVsP[SYS_START][locPID]->Fill(locP, locDeltaBeta);
1813 double locDeltaT = locSCHitMatchParams->dHitTime - locSCHitMatchParams->dFlightTime - locStartTime;
1814 dHistMap_DeltaTVsP[SYS_START][locPID]->Fill(locP, locDeltaT);
1815 }
1816 }
1817 if(dHistMap_DeltadEdXVsP[SYS_START].find(locPID) != dHistMap_DeltadEdXVsP[SYS_START].end())
1818 {
1819 double locdx = locSCHitMatchParams->dHitEnergy/locSCHitMatchParams->dEdx;
1820 double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
1821 locParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
1822 dHistMap_DeltadEdXVsP[SYS_START][locPID]->Fill(locP, (locSCHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
1823 }
1824 }
1825 if(locTOFHitMatchParams != NULL__null)
1826 {
1827 dHistMap_dEdXVsP[SYS_TOF][locCharge]->Fill(locP, locTOFHitMatchParams->dEdx*1.0E3);
1828 if(locEventRFBunch->dNumParticleVotes > 1)
1829 {
1830 double locBeta_Timing = locTOFHitMatchParams->dPathLength/(29.9792458*(locTOFHitMatchParams->dHitTime - locChargedTrackHypothesis->t0()));
1831 dHistMap_BetaVsP[SYS_TOF][locCharge]->Fill(locP, locBeta_Timing);
1832 if(dHistMap_DeltaBetaVsP[SYS_TOF].find(locPID) != dHistMap_DeltaBetaVsP[SYS_TOF].end())
1833 {
1834 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1835 dHistMap_DeltaBetaVsP[SYS_TOF][locPID]->Fill(locP, locDeltaBeta);
1836 double locDeltaT = locTOFHitMatchParams->dHitTime - locTOFHitMatchParams->dFlightTime - locStartTime;
1837 dHistMap_DeltaTVsP[SYS_TOF][locPID]->Fill(locP, locDeltaT);
1838 }
1839 }
1840 if(dHistMap_DeltadEdXVsP[SYS_TOF].find(locPID) != dHistMap_DeltadEdXVsP[SYS_TOF].end())
1841 {
1842 double locdx = locTOFHitMatchParams->dHitEnergy/locTOFHitMatchParams->dEdx;
1843 double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
1844 locParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
1845 dHistMap_DeltadEdXVsP[SYS_TOF][locPID]->Fill(locP, (locTOFHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
1846 }
1847 }
1848 if(locBCALShowerMatchParams != NULL__null)
1849 {
1850 const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
1851 double locEOverP = locBCALShower->E/locP;
1852 dHistMap_BCALEOverPVsP[locCharge]->Fill(locP, locEOverP);
1853 dHistMap_BCALEOverPVsTheta[locCharge]->Fill(locTheta, locEOverP);
1854 if(locEventRFBunch->dNumParticleVotes > 1)
1855 {
1856 double locBeta_Timing = locBCALShowerMatchParams->dPathLength/(29.9792458*(locBCALShower->t - locChargedTrackHypothesis->t0()));
1857 dHistMap_BetaVsP[SYS_BCAL][locCharge]->Fill(locP, locBeta_Timing);
1858 if(dHistMap_DeltaBetaVsP[SYS_BCAL].find(locPID) != dHistMap_DeltaBetaVsP[SYS_BCAL].end())
1859 {
1860 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1861 dHistMap_DeltaBetaVsP[SYS_BCAL][locPID]->Fill(locP, locDeltaBeta);
1862 double locDeltaT = locBCALShower->t - locBCALShowerMatchParams->dFlightTime - locStartTime;
1863 dHistMap_DeltaTVsP[SYS_BCAL][locPID]->Fill(locP, locDeltaT);
1864 }
1865 }
1866 }
1867 if(locFCALShowerMatchParams != NULL__null)
1868 {
1869 const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
1870 double locEOverP = locFCALShower->getEnergy()/locP;
1871 dHistMap_FCALEOverPVsP[locCharge]->Fill(locP, locEOverP);
1872 dHistMap_FCALEOverPVsTheta[locCharge]->Fill(locTheta, locEOverP);
1873 if(locEventRFBunch->dNumParticleVotes > 1)
1874 {
1875 double locBeta_Timing = locFCALShowerMatchParams->dPathLength/(29.9792458*(locFCALShower->getTime() - locChargedTrackHypothesis->t0()));
1876 dHistMap_BetaVsP[SYS_FCAL][locCharge]->Fill(locP, locBeta_Timing);
1877 if(dHistMap_DeltaBetaVsP[SYS_FCAL].find(locPID) != dHistMap_DeltaBetaVsP[SYS_FCAL].end())
1878 {
1879 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1880 dHistMap_DeltaBetaVsP[SYS_FCAL][locPID]->Fill(locP, locDeltaBeta);
1881 double locDeltaT = locFCALShower->getTime() - locFCALShowerMatchParams->dFlightTime - locStartTime;
1882 dHistMap_DeltaTVsP[SYS_FCAL][locPID]->Fill(locP, locDeltaT);
1883 }
1884 }
1885 }
1886
1887 if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC > 0)
1888 {
1889 dHistMap_dEdXVsP[SYS_CDC][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_CDC*1.0E6);
1890 if(dHistMap_DeltadEdXVsP[SYS_CDC].find(locPID) != dHistMap_DeltadEdXVsP[SYS_CDC].end())
1891 {
1892 double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_CDC, true);
1893 dHistMap_DeltadEdXVsP[SYS_CDC][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_CDC - locProbabledEdx)*1.0E6);
1894 }
1895 }
1896 if(locTrackTimeBased->dNumHitsUsedFordEdx_FDC > 0)
1897 {
1898 dHistMap_dEdXVsP[SYS_FDC][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_FDC*1.0E6);
1899 if(dHistMap_DeltadEdXVsP[SYS_FDC].find(locPID) != dHistMap_DeltadEdXVsP[SYS_FDC].end())
1900 {
1901 double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_FDC, false);
1902 dHistMap_DeltadEdXVsP[SYS_FDC][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_FDC - locProbabledEdx)*1.0E6);
1903 }
1904 }
1905 }
1906 }
1907 Unlock_Action(); //RELEASE ROOT LOCK!!
1908
1909 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
1910}
1911
1912void DHistogramAction_Neutrals::Initialize(JEventLoop* locEventLoop)
1913{
1914 //Create any histograms/trees/etc. within a ROOT lock.
1915 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
1916
1917 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
1918 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
1919
1920 string locHistName;
1921
1922 DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
1923 DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
1924 double locTargetZCenter = 0.0;
1925 locGeometry->GetTargetZ(locTargetZCenter);
1926
1927 //CREATE THE HISTOGRAMS
1928 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1929 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1930 {
1931 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
1932 //If another thread has already created the folder, it just changes to it.
1933 CreateAndChangeTo_ActionDirectory();
1934
1935 if(dTargetCenter.Z() < -9.8E9)
1936 dTargetCenter.SetXYZ(0.0, 0.0, locTargetZCenter);
1937
1938 //BCAL
1939 locHistName = "BCALTrackDOCA";
1940 dHist_BCALTrackDOCA = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower Distance to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
1941 locHistName = "BCALTrackDeltaPhi";
1942 dHist_BCALTrackDeltaPhi = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #Delta#phi#circ to Nearest Track", dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
1943 locHistName = "BCALTrackDeltaZ";
1944 dHist_BCALTrackDeltaZ = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #DeltaZ to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
1945 locHistName = "BCALNeutralShowerEnergy";
1946 dHist_BCALNeutralShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Neutral Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxBCALP);
1947 locHistName = "BCALNeutralShowerDeltaT";
1948 dHist_BCALNeutralShowerDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Neutral Shower #Deltat (Propagated - RF) (ns)", dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1949 locHistName = "BCALNeutralShowerDeltaTVsE";
1950 dHist_BCALNeutralShowerDeltaTVsE = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Neutral Shower Energy (GeV);BCAL Neutral Shower #Deltat (ns)", dNum2DShowerEnergyBins, dMinShowerEnergy, dMaxBCALP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1951 locHistName = "BCALNeutralShowerDeltaTVsZ";
1952 dHist_BCALNeutralShowerDeltaTVsZ = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Neutral Shower Z (cm);BCAL Neutral Shower #Deltat (ns)", dNum2DBCALZBins, 0.0, 450.0, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1953
1954 //FCAL
1955 locHistName = "FCALTrackDOCA";
1956 dHist_FCALTrackDOCA = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Shower Distance to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
1957 locHistName = "FCALNeutralShowerEnergy";
1958 dHist_FCALNeutralShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Neutral Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
1959 locHistName = "FCALNeutralShowerDeltaT";
1960 dHist_FCALNeutralShowerDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Neutral Shower #Deltat (Propagated - RF) (ns)", dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
1961 locHistName = "FCALNeutralShowerDeltaTVsE";
1962 dHist_FCALNeutralShowerDeltaTVsE = GetOrCreate_Histogram<TH2I>(locHistName, ";FCAL Neutral Shower Energy (GeV);FCAL Neutral Shower #Deltat (ns)", dNum2DShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1963
1964 //Return to the base directory
1965 ChangeTo_BaseDirectory();
1966 }
1967 japp->RootUnLock(); //RELEASE ROOT LOCK!!
1968}
1969
1970bool DHistogramAction_Neutrals::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1971{
1972 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
1973
1974 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
1975 if(Get_NumPreviousParticleCombos() != 0)
1976 return true;
1977
1978 vector<const DNeutralShower*> locNeutralShowers;
1979 locEventLoop->Get(locNeutralShowers);
1980
1981 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
1982 locEventLoop->Get(locTrackTimeBasedVector);
1983
1984 const DDetectorMatches* locDetectorMatches = NULL__null;
1985 locEventLoop->GetSingle(locDetectorMatches);
1986
1987 vector<const DEventRFBunch*> locEventRFBunches;
1988 locEventLoop->Get(locEventRFBunches);
1989 double locStartTime = locEventRFBunches.empty() ? 0.0 : locEventRFBunches[0]->dTime;
1990
1991 //FILL HISTOGRAMS
1992 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1993 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1994 Lock_Action(); //ACQUIRE ROOT LOCK!!
1995 {
1996 for(size_t loc_i = 0; loc_i < locNeutralShowers.size(); ++loc_i)
1997 {
1998 //assume is photon
1999 double locPathLength = (locNeutralShowers[loc_i]->dSpacetimeVertex.Vect() - dTargetCenter).Mag();
2000 double locDeltaT = locNeutralShowers[loc_i]->dSpacetimeVertex.T() - locPathLength/29.9792458 - locStartTime;
2001
2002 if(locNeutralShowers[loc_i]->dDetectorSystem == SYS_FCAL)
2003 {
2004 const DFCALShower* locFCALShower = NULL__null;
2005 locNeutralShowers[loc_i]->GetSingle(locFCALShower);
2006
2007 double locDistance = 9.9E9;
2008 if(locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistance))
2009 dHist_FCALTrackDOCA->Fill(locDistance);
2010
2011 dHist_FCALNeutralShowerEnergy->Fill(locNeutralShowers[loc_i]->dEnergy);
2012 dHist_FCALNeutralShowerDeltaT->Fill(locDeltaT);
2013 dHist_FCALNeutralShowerDeltaTVsE->Fill(locNeutralShowers[loc_i]->dEnergy, locDeltaT);
2014 }
2015 else
2016 {
2017 const DBCALShower* locBCALShower = NULL__null;
2018 locNeutralShowers[loc_i]->GetSingle(locBCALShower);
2019
2020 double locDistance = 9.9E9, locDeltaPhi = 9.9E9, locDeltaZ = 9.9E9;
2021 if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDistance))
2022 dHist_BCALTrackDOCA->Fill(locDistance);
2023 if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDeltaPhi, locDeltaZ))
2024 {
2025 dHist_BCALTrackDeltaPhi->Fill(180.0*locDeltaPhi/TMath::Pi());
2026 dHist_BCALTrackDeltaZ->Fill(locDeltaZ);
2027 }
2028
2029 dHist_BCALNeutralShowerEnergy->Fill(locNeutralShowers[loc_i]->dEnergy);
2030 dHist_BCALNeutralShowerDeltaT->Fill(locDeltaT);
2031 dHist_BCALNeutralShowerDeltaTVsE->Fill(locNeutralShowers[loc_i]->dEnergy, locDeltaT);
2032 dHist_BCALNeutralShowerDeltaTVsZ->Fill(locNeutralShowers[loc_i]->dSpacetimeVertex.Z(), locDeltaT);
2033 }
2034 }
2035 }
2036 Unlock_Action(); //RELEASE ROOT LOCK!!
2037
2038 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2039}
2040
2041
2042void DHistogramAction_DetectorMatchParams::Initialize(JEventLoop* locEventLoop)
2043{
2044 //Create any histograms/trees/etc. within a ROOT lock.
2045 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
2046
2047 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
2048 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
2049
2050 string locHistName, locHistTitle;
2051
2052 vector<const DMCThrown*> locMCThrowns;
2053 locEventLoop->Get(locMCThrowns);
2054
2055 DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
2056 DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
2057 double locTargetZCenter = 0.0;
2058 locGeometry->GetTargetZ(locTargetZCenter);
2059
2060 string locTrackSelectionTag = "NotATag";
2061 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2062 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2063
2064 //CREATE THE HISTOGRAMS
2065 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2066 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2067 {
2068 //So: Default tag is "", User can set it to something else
2069 //In here, if tag is "", get from gparms, if not, leave it alone
2070 //If gparms value does not exist, set it to (and use) "PreSelect"
2071 if(dTrackSelectionTag == "NotATag")
2072 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2073
2074 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
2075 //If another thread has already created the folder, it just changes to it.
2076 CreateAndChangeTo_ActionDirectory();
2077
2078 if(dTargetCenterZ < -9.8E9)
2079 dTargetCenterZ = locTargetZCenter; //only set if not already set
2080
2081 //Track Matched to Hit
2082 for(int locTruePIDFlag = 0; locTruePIDFlag < 2; ++locTruePIDFlag)
2083 {
2084 if(locMCThrowns.empty() && (locTruePIDFlag == 1))
2085 continue; //not a simulated event: don't histogram thrown info!
2086
2087 string locDirName = (locTruePIDFlag == 1) ? "TruePID" : "ReconstructedPID";
2088 CreateAndChangeTo_Directory(locDirName.c_str(), locDirName.c_str());
2089
2090 //By PID
2091 for(int loc_i = -2; loc_i < int(dTrackingPIDs.size()); ++loc_i) //-2 = q-, -1 = q+
2092 {
2093 string locParticleName, locParticleROOTName;
2094 int locPID = loc_i;
2095 if(loc_i == -2)
2096 {
2097 locParticleName = "q-";
2098 locParticleROOTName = "#it{q}^{-}";
2099 }
2100 else if(loc_i == -1)
2101 {
2102 locParticleName = "q+";
2103 locParticleROOTName = "#it{q}^{+}";
2104 }
2105 else
2106 {
2107 locParticleName = ParticleType(dTrackingPIDs[loc_i]);
2108 locParticleROOTName = ParticleName_ROOT(dTrackingPIDs[loc_i]);
2109 locPID = int(dTrackingPIDs[loc_i]);
2110 }
2111 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2112 pair<int, bool> locPIDPair(locPID, bool(locTruePIDFlag));
2113
2114 //BCAL
2115 locHistName = "BCALShowerEnergy";
2116 locHistTitle = locParticleROOTName + ";BCAL Shower Energy (GeV)";
2117 dHistMap_BCALShowerEnergy[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerEnergyBins, dMinShowerEnergy, dMaxBCALP);
2118
2119 locHistName = "BCALShowerTrackDepth";
2120 locHistTitle = locParticleROOTName + ";BCAL Shower Track Depth (cm)";
2121 dHistMap_BCALShowerTrackDepth[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2122
2123 locHistName = "BCALShowerTrackDepthVsP";
2124 locHistTitle = locParticleROOTName + ";p (GeV/c);BCAL Shower Track Depth (cm)";
2125 dHistMap_BCALShowerTrackDepthVsP[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2126
2127 //FCAL
2128 locHistName = "FCALShowerEnergy";
2129 locHistTitle = locParticleROOTName + ";FCAL Shower Energy (GeV)";
2130 dHistMap_FCALShowerEnergy[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2131
2132 locHistName = "FCALShowerTrackDepth";
2133 locHistTitle = locParticleROOTName + ";FCAL Shower Track Depth (cm)";
2134 dHistMap_FCALShowerTrackDepth[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2135
2136 locHistName = "FCALShowerTrackDepthVsP";
2137 locHistTitle = locParticleROOTName + ";p (GeV/c);FCAL Shower Track Depth (cm)";
2138 dHistMap_FCALShowerTrackDepthVsP[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2139
2140 //SC
2141 locHistName = "SCEnergyVsTheta";
2142 locHistTitle = locParticleROOTName + ";#theta#circ;SC Point Energy (MeV)";
2143 dHistMap_SCEnergyVsTheta[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
2144
2145 locHistName = "SCPhiVsTheta";
2146 locHistTitle = locParticleROOTName + ";#theta#circ;#phi#circ";
2147 dHistMap_SCPhiVsTheta[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
2148
2149 gDirectory(TDirectory::CurrentDirectory())->cd(".."); //end of PID
2150 }
2151 gDirectory(TDirectory::CurrentDirectory())->cd(".."); //end of true/recon PID
2152 }
2153
2154 //Return to the base directory
2155 ChangeTo_BaseDirectory();
2156 }
2157 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2158}
2159
2160bool DHistogramAction_DetectorMatchParams::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2161{
2162 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
2163
2164 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
2165 if(Get_NumPreviousParticleCombos() != 0)
2166 return true;
2167
2168 vector<const DMCThrown*> locMCThrowns;
2169 locEventLoop->Get(locMCThrowns);
2170
2171 Fill_Hists(locEventLoop, false);
2172 if(!locMCThrowns.empty())
2173 Fill_Hists(locEventLoop, true);
2174
2175 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2176}
2177
2178void DHistogramAction_DetectorMatchParams::Fill_Hists(JEventLoop* locEventLoop, bool locUseTruePIDFlag)
2179{
2180 vector<const DChargedTrack*> locChargedTracks;
2181 locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
2182
2183 vector<const DMCThrownMatching*> locMCThrownMatchingVector;
2184 locEventLoop->Get(locMCThrownMatchingVector);
2185
2186 const DEventRFBunch* locEventRFBunch = NULL__null;
2187 locEventLoop->GetSingle(locEventRFBunch);
2188
2189 //FILL HISTOGRAMS
2190 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2191 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2192 Lock_Action(); //ACQUIRE ROOT LOCK!!
2193 {
2194 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2195 {
2196 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2197
2198 if(locUseTruePIDFlag && (!locMCThrownMatchingVector.empty()))
2199 {
2200 double locMatchFOM = 0.0;
2201 const DMCThrown* locMCThrown = locMCThrownMatchingVector[0]->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);
2202 if(locMCThrown == NULL__null)
2203 continue;
2204 //OK, have the thrown. Now, grab the best charged track hypothesis to get the best matching
2205 locChargedTrackHypothesis = locMCThrownMatchingVector[0]->Get_MatchingChargedHypothesis(locMCThrown, locMatchFOM);
2206 }
2207
2208 pair<int, bool> locPIDPair(int(locChargedTrackHypothesis->PID()), locUseTruePIDFlag);
2209 bool locDisregardPIDFlag = (dHistMap_BCALShowerEnergy.find(locPIDPair) == dHistMap_BCALShowerEnergy.end());
2210 int locQIndex = (locChargedTrackHypothesis->charge() > 0.0) ? -1 : -2;
2211 pair<int, bool> locChargePair(locQIndex, locUseTruePIDFlag);
2212
2213 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2214 const DFCALShowerMatchParams* locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
2215 const DSCHitMatchParams* locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
2216 const DBCALShowerMatchParams* locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
2217
2218 //BCAL
2219 if(locBCALShowerMatchParams != NULL__null)
2220 {
2221 const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
2222 dHistMap_BCALShowerEnergy[locChargePair]->Fill(locBCALShower->E);
2223 dHistMap_BCALShowerTrackDepth[locChargePair]->Fill(locBCALShowerMatchParams->dx);
2224 dHistMap_BCALShowerTrackDepthVsP[locChargePair]->Fill(locMomentum.Mag(), locBCALShowerMatchParams->dx);
2225
2226 if(!locDisregardPIDFlag)
2227 {
2228 dHistMap_BCALShowerEnergy[locPIDPair]->Fill(locBCALShower->E);
2229 dHistMap_BCALShowerTrackDepth[locPIDPair]->Fill(locBCALShowerMatchParams->dx);
2230 dHistMap_BCALShowerTrackDepthVsP[locPIDPair]->Fill(locMomentum.Mag(), locBCALShowerMatchParams->dx);
2231 }
2232 }
2233
2234 //FCAL
2235 if(locFCALShowerMatchParams != NULL__null)
2236 {
2237 const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
2238 dHistMap_FCALShowerEnergy[locChargePair]->Fill(locFCALShower->getEnergy());
2239 dHistMap_FCALShowerTrackDepth[locChargePair]->Fill(locFCALShowerMatchParams->dx);
2240 dHistMap_FCALShowerTrackDepthVsP[locChargePair]->Fill(locMomentum.Mag(), locFCALShowerMatchParams->dx);
2241
2242 if(!locDisregardPIDFlag)
2243 {
2244 dHistMap_FCALShowerEnergy[locPIDPair]->Fill(locFCALShower->getEnergy());
2245 dHistMap_FCALShowerTrackDepth[locPIDPair]->Fill(locFCALShowerMatchParams->dx);
2246 dHistMap_FCALShowerTrackDepthVsP[locPIDPair]->Fill(locMomentum.Mag(), locFCALShowerMatchParams->dx);
2247 }
2248 }
2249
2250 //SC
2251 if(locSCHitMatchParams != NULL__null)
2252 {
2253 dHistMap_SCEnergyVsTheta[locChargePair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locSCHitMatchParams->dHitEnergy*1.0E3);
2254 dHistMap_SCPhiVsTheta[locChargePair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locMomentum.Phi()*180.0/TMath::Pi());
2255
2256 if(!locDisregardPIDFlag)
2257 {
2258 dHistMap_SCEnergyVsTheta[locPIDPair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locSCHitMatchParams->dHitEnergy*1.0E3);
2259 dHistMap_SCPhiVsTheta[locPIDPair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locMomentum.Phi()*180.0/TMath::Pi());
2260 }
2261 }
2262 }
2263 }
2264 Unlock_Action(); //RELEASE ROOT LOCK!!
2265}
2266
2267void DHistogramAction_EventVertex::Initialize(JEventLoop* locEventLoop)
2268{
2269 string locHistName, locHistTitle;
2270
2271 string locTrackSelectionTag = "NotATag";
2272 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2273 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2274
2275 //CREATE THE HISTOGRAMS
2276 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2277 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2278 {
2279 //So: Default tag is "", User can set it to something else
2280 //In here, if tag is "", get from gparms, if not, leave it alone
2281 //If gparms value does not exist, set it to (and use) "PreSelect"
2282 if(dTrackSelectionTag == "NotATag")
2283 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2284
2285 CreateAndChangeTo_ActionDirectory();
2286
2287 // Event RF Bunch Time
2288 locHistName = "RFTrackDeltaT";
2289 locHistTitle = ";#Deltat_{RF - Track} (ns)";
2290 dRFTrackDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumRFTBins, -3.0, 3.0);
2291
2292 // ALL EVENTS
2293 CreateAndChangeTo_Directory("AllEvents", "AllEvents");
2294
2295 // Event Vertex-Z
2296 locHistName = "EventVertexZ";
2297 locHistTitle = ";Event Vertex-Z (cm)";
2298 dEventVertexZ_AllEvents = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2299
2300 // Event Vertex-Y Vs Vertex-X
2301 locHistName = "EventVertexYVsX";
2302 locHistTitle = ";Event Vertex-X (cm);Event Vertex-Y (cm)";
2303 dEventVertexYVsX_AllEvents = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2304
2305 // Event Vertex-T
2306 locHistName = "EventVertexT";
2307 locHistTitle = ";Event Vertex Time (ns)";
2308 dEventVertexT_AllEvents = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2309
2310 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2311
2312
2313 // 2+ Good Tracks
2314 CreateAndChangeTo_Directory("2+GoodTracks", "2+GoodTracks");
2315
2316 // Event Vertex-Z
2317 locHistName = "EventVertexZ";
2318 locHistTitle = ";Event Vertex-Z (cm)";
2319 dEventVertexZ_2OrMoreGoodTracks = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2320
2321 // Event Vertex-Y Vs Vertex-X
2322 locHistName = "EventVertexYVsX";
2323 locHistTitle = ";Event Vertex-X (cm);Event Vertex-Y (cm)";
2324 dEventVertexYVsX_2OrMoreGoodTracks = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2325
2326 // Event Vertex-T
2327 locHistName = "EventVertexT";
2328 locHistTitle = ";Event Vertex Time (ns)";
2329 dEventVertexT_2OrMoreGoodTracks = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2330
2331 // Confidence Level
2332 locHistName = "ConfidenceLevel";
2333 dHist_KinFitConfidenceLevel = GetOrCreate_Histogram<TH1I>(locHistName, "Event Vertex Kinematic Fit;Confidence Level;# Events", dNumConfidenceLevelBins, 0.0, 1.0);
2334
2335 //final particle pulls
2336 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2337 {
2338 Particle_t locPID = dFinalStatePIDs[loc_i];
2339 string locParticleDirName = string("Pulls_") + string(ParticleType(locPID));
2340 string locParticleROOTName = ParticleName_ROOT(locPID);
2341 CreateAndChangeTo_Directory(locParticleDirName, locParticleDirName);
2342
2343 //Px Pull
2344 locHistName = "Pull_Px";
2345 locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{x} Pull;# Events");
2346 dHistMap_KinFitPulls[locPID][d_PxPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2347
2348 //Py Pull
2349 locHistName = "Pull_Py";
2350 locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{y} Pull;# Events");
2351 dHistMap_KinFitPulls[locPID][d_PyPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2352
2353 //Pz Pull
2354 locHistName = "Pull_Pz";
2355 locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{z} Pull;# Events");
2356 dHistMap_KinFitPulls[locPID][d_PzPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2357
2358 //Xx Pull
2359 locHistName = "Pull_Xx";
2360 locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{x} Pull;# Events");
2361 dHistMap_KinFitPulls[locPID][d_XxPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2362
2363 //Xy Pull
2364 locHistName = "Pull_Xy";
2365 locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{y} Pull;# Events");
2366 dHistMap_KinFitPulls[locPID][d_XyPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2367
2368 //Xz Pull
2369 locHistName = "Pull_Xz";
2370 locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{z} Pull;# Events");
2371 dHistMap_KinFitPulls[locPID][d_XzPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2372
2373 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2374 }
2375
2376 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2377
2378 //Return to the base directory
2379 ChangeTo_BaseDirectory();
2380 }
2381 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2382}
2383
2384bool DHistogramAction_EventVertex::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2385{
2386 if(Get_NumPreviousParticleCombos() != 0)
2387 return true; //else double-counting!
2388
2389 const DVertex* locVertex = NULL__null;
2390 locEventLoop->GetSingle(locVertex);
2391
2392 vector<const DChargedTrack*> locChargedTracks;
2393 locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
2394
2395 const DDetectorMatches* locDetectorMatches = NULL__null;
2396 locEventLoop->GetSingle(locDetectorMatches);
2397
2398 const DParticleID* locParticleID = NULL__null;
2399 locEventLoop->GetSingle(locParticleID);
2400
2401 const DEventRFBunch* locEventRFBunch = NULL__null;
2402 locEventLoop->GetSingle(locEventRFBunch);
2403
2404 //Make sure that brun() is called (to get rf period) before using.
2405 //Cannot call JEventLoop->Get() because object may be in datastream (REST), bypassing factory brun() call.
2406 //Must do here rather than in Initialize() function because this object is shared by all threads (which each have their own factory)
2407 DRFTime_factory* locRFTimeFactory = static_cast<DRFTime_factory*>(locEventLoop->GetFactory("DRFTime"));
2408 if(!locRFTimeFactory->brun_was_called())
2409 {
2410 locRFTimeFactory->brun(locEventLoop, locEventLoop->GetJEvent().GetRunNumber());
2411 locRFTimeFactory->Set_brun_called();
2412 }
2413
2414 //Get time-based tracks: use best PID FOM
2415 //Note that these may not be the PIDs that were used in the fit!!!
2416 //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
2417 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
2418 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2419 {
2420 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2421 const DTrackTimeBased* locTrackTimeBased = NULL__null;
2422 locChargedTrackHypothesis->GetSingle(locTrackTimeBased);
2423 if(locTrackTimeBased != NULL__null)
2424 locTrackTimeBasedVector.push_back(locTrackTimeBased);
2425 }
2426
2427 //Event Vertex
2428 //FILL HISTOGRAMS
2429 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2430 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2431 Lock_Action();
2432 {
2433 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2434 {
2435 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2436 double locPropagatedRFTime = locParticleID->Calc_PropagatedRFTime(locChargedTrackHypothesis, locEventRFBunch);
2437 double locShiftedRFTime = locRFTimeFactory->Step_TimeToNearInputTime(locPropagatedRFTime, locChargedTrackHypothesis->time());
2438 double locDeltaT = locShiftedRFTime - locChargedTrackHypothesis->time();
2439 dRFTrackDeltaT->Fill(locDeltaT);
2440 }
2441 dEventVertexZ_AllEvents->Fill(locVertex->dSpacetimeVertex.Z());
2442 dEventVertexYVsX_AllEvents->Fill(locVertex->dSpacetimeVertex.X(), locVertex->dSpacetimeVertex.Y());
2443 dEventVertexT_AllEvents->Fill(locVertex->dSpacetimeVertex.T());
2444
2445 if(locChargedTracks.size() >= 2)
2446 {
2447 dEventVertexZ_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.Z());
2448 dEventVertexYVsX_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.X(), locVertex->dSpacetimeVertex.Y());
2449 dEventVertexT_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.T());
2450 }
2451 }
2452 Unlock_Action();
2453
2454 if(locVertex->dKinFitNDF == 0)
2455 return true; //kin fit not performed or didn't converge: no results to histogram
2456
2457 double locConfidenceLevel = TMath::Prob(locVertex->dKinFitChiSq, locVertex->dKinFitNDF);
2458
2459 //FILL HISTOGRAMS
2460 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2461 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2462 Lock_Action();
2463 {
2464 dHist_KinFitConfidenceLevel->Fill(locConfidenceLevel);
2465
2466 //pulls
2467 if(locConfidenceLevel > dPullHistConfidenceLevelCut)
2468 {
2469 for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
2470 {
2471 const DKinematicData* locKinematicData = static_cast<const DKinematicData*>(locTrackTimeBasedVector[loc_i]);
2472 Particle_t locPID = locKinematicData->PID();
2473 if(dHistMap_KinFitPulls.find(locPID) == dHistMap_KinFitPulls.end())
2474 continue; //PID not histogrammed
2475
2476 map<const JObject*, map<DKinFitPullType, double> >::const_iterator locParticleIterator = locVertex->dKinFitPulls.find(locKinematicData);
2477 if(locParticleIterator == locVertex->dKinFitPulls.end())
2478 continue;
2479
2480 const map<DKinFitPullType, double>& locPullMap = locParticleIterator->second;
2481 map<DKinFitPullType, double>::const_iterator locPullIterator = locPullMap.begin();
2482 for(; locPullIterator != locPullMap.end(); ++locPullIterator)
2483 dHistMap_KinFitPulls[locPID][locPullIterator->first]->Fill(locPullIterator->second);
2484 }
2485 }
2486 }
2487 Unlock_Action();
2488
2489 return true;
2490}
2491
2492void DHistogramAction_DetectedParticleKinematics::Initialize(JEventLoop* locEventLoop)
2493{
2494 string locHistName, locHistTitle, locParticleName, locParticleROOTName;
2495 Particle_t locPID;
2496
2497 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
2498 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2499 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2500 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
2501 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
2502
2503 //CREATE THE HISTOGRAMS
2504 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2505 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2506 {
2507 //So: Default tag is "", User can set it to something else
2508 //In here, if tag is "", get from gparms, if not, leave it alone
2509 //If gparms value does not exist, set it to (and use) "PreSelect"
2510 if(dTrackSelectionTag == "NotATag")
2511 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2512 if(dShowerSelectionTag == "NotATag")
2513 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
2514
2515 CreateAndChangeTo_ActionDirectory();
2516
2517 // Beam Particle
2518 locPID = Gamma;
2519 locParticleName = string("Beam_") + ParticleType(locPID);
2520 locParticleROOTName = ParticleName_ROOT(locPID);
2521 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2522 locHistName = "Momentum";
2523 locHistTitle = string("Beam ") + locParticleROOTName + string(";p (GeV/c)");
2524 dBeamParticle_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBeamEBins, dMinP, dMaxBeamE);
2525 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2526
2527 //PID
2528 CreateAndChangeTo_Directory("PID", "PID");
2529 {
2530 //beta vs p
2531 locHistName = "BetaVsP_Q+";
2532 locHistTitle = "q^{+};p (GeV/c);#beta";
2533 dHistMap_QBetaVsP[1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2534
2535 locHistName = "BetaVsP_Q-";
2536 locHistTitle = "q^{-};p (GeV/c);#beta";
2537 dHistMap_QBetaVsP[-1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2538 }
2539 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2540
2541 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2542 {
2543 locPID = dFinalStatePIDs[loc_i];
2544 locParticleName = ParticleType(locPID);
2545 locParticleROOTName = ParticleName_ROOT(locPID);
2546 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2547
2548 // Momentum
2549 locHistName = "Momentum";
2550 locHistTitle = locParticleROOTName + string(";p (GeV/c)");
2551 dHistMap_P[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
2552
2553 // Theta
2554 locHistName = "Theta";
2555 locHistTitle = locParticleROOTName + string(";#theta#circ");
2556 dHistMap_Theta[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
2557
2558 // Phi
2559 locHistName = "Phi";
2560 locHistTitle = locParticleROOTName + string(";#phi#circ");
2561 dHistMap_Phi[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
2562
2563 // P Vs Theta
2564 locHistName = "PVsTheta";
2565 locHistTitle = locParticleROOTName + string(";#theta#circ;p (GeV/c)");
2566 dHistMap_PVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
2567
2568 // Phi Vs Theta
2569 locHistName = "PhiVsTheta";
2570 locHistTitle = locParticleROOTName + string(";#theta#circ;#phi#circ");
2571 dHistMap_PhiVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
2572
2573 //beta vs p
2574 locHistName = "BetaVsP";
2575 locHistTitle = locParticleROOTName + string(";p (GeV/c);#beta");
2576 dHistMap_BetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2577
2578 //delta-beta vs p
2579 locHistName = "DeltaBetaVsP";
2580 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#beta");
2581 dHistMap_DeltaBetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
2582
2583 // Vertex-Z
2584 locHistName = "VertexZ";
2585 locHistTitle = locParticleROOTName + string(";Vertex-Z (cm)");
2586 dHistMap_VertexZ[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2587
2588 // Vertex-Y Vs Vertex-X
2589 locHistName = "VertexYVsX";
2590 locHistTitle = locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
2591 dHistMap_VertexYVsX[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2592
2593 // Vertex-T
2594 locHistName = "VertexT";
2595 locHistTitle = locParticleROOTName + string(";Vertex-T (ns)");
2596 dHistMap_VertexT[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2597
2598 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2599 }
2600
2601 //Return to the base directory
2602 ChangeTo_BaseDirectory();
2603 }
2604 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2605}
2606
2607bool DHistogramAction_DetectedParticleKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2608{
2609 if(Get_NumPreviousParticleCombos() != 0)
2610 return true; //else double-counting!
2611
2612 vector<const DBeamPhoton*> locBeamPhotons;
2613 locEventLoop->Get(locBeamPhotons);
2614
2615 //FILL HISTOGRAMS
2616 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2617 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2618 Lock_Action();
2619 {
2620 for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
2621 dBeamParticle_P->Fill(locBeamPhotons[loc_i]->energy());
2622 }
2623 Unlock_Action();
2624
2625 vector<const DChargedTrack*> locPreSelectChargedTracks;
2626 locEventLoop->Get(locPreSelectChargedTracks, dTrackSelectionTag.c_str());
2627
2628 for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
2629 {
2630 const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestTrackingFOM();
2631 int locCharge = ParticleCharge(locChargedTrackHypothesis->PID());
2632
2633 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2634 double locP = locMomentum.Mag();
2635 double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
2636
2637 if(dHistMap_QBetaVsP.find(locCharge) == dHistMap_QBetaVsP.end())
2638 continue;
2639
2640 //FILL HISTOGRAMS
2641 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2642 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2643 Lock_Action();
2644 {
2645 //Extremely inefficient, I know ...
2646 dHistMap_QBetaVsP[locCharge]->Fill(locP, locBeta_Timing);
2647 }
2648 Unlock_Action();
2649 }
2650
2651 for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
2652 {
2653 const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestFOM();
2654
2655 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2656 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2657 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2658 double locP = locMomentum.Mag();
2659 double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
2660 double locDeltaBeta = locChargedTrackHypothesis->deltaBeta();
2661
2662 Particle_t locPID = (locChargedTrackHypothesis->dFOM < dMinPIDFOM) ? Unknown : locChargedTrackHypothesis->PID();
2663 if(dHistMap_P.find(locPID) == dHistMap_P.end())
2664 continue; //not interested in histogramming
2665
2666 //FILL HISTOGRAMS
2667 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2668 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2669 Lock_Action();
2670 {
2671 dHistMap_P[locPID]->Fill(locP);
2672 dHistMap_Phi[locPID]->Fill(locPhi);
2673 dHistMap_Theta[locPID]->Fill(locTheta);
2674 dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
2675 dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
2676 dHistMap_BetaVsP[locPID]->Fill(locP, locBeta_Timing);
2677 dHistMap_DeltaBetaVsP[locPID]->Fill(locP, locDeltaBeta);
2678 dHistMap_VertexZ[locPID]->Fill(locChargedTrackHypothesis->position().Z());
2679 dHistMap_VertexYVsX[locPID]->Fill(locChargedTrackHypothesis->position().X(), locChargedTrackHypothesis->position().Y());
2680 dHistMap_VertexT[locPID]->Fill(locChargedTrackHypothesis->time());
2681 }
2682 Unlock_Action();
2683 }
2684
2685 vector<const DNeutralParticle*> locNeutralParticles;
2686 locEventLoop->Get(locNeutralParticles);
2687
2688 vector<const DNeutralShower*> locNeutralShowers;
2689 locEventLoop->Get(locNeutralShowers, dShowerSelectionTag.c_str());
2690
2691 for(size_t loc_i = 0; loc_i < locNeutralShowers.size(); ++loc_i)
2692 {
2693 const DNeutralParticle* locNeutralParticle = NULL__null;
2694 for(size_t loc_j = 0; loc_j < locNeutralParticles.size(); ++loc_j)
2695 {
2696 if(locNeutralParticles[loc_j]->dNeutralShower != locNeutralShowers[loc_i])
2697 continue;
2698 locNeutralParticle = locNeutralParticles[loc_j];
2699 break;
2700 }
2701 if(locNeutralParticle == NULL__null)
2702 continue;
2703 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticle->Get_Hypothesis(Gamma);
2704 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
2705 continue;
2706
2707 Particle_t locPID = locNeutralParticleHypothesis->PID();
2708
2709 if(dHistMap_P.find(locPID) == dHistMap_P.end())
2710 continue; //e.g. a decaying particle, or not interested in histogramming
2711
2712 DVector3 locMomentum = locNeutralParticleHypothesis->momentum();
2713 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2714 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2715 double locP = locMomentum.Mag();
2716
2717 double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
2718 double locDeltaBeta = locNeutralParticleHypothesis->deltaBeta();
2719
2720 //FILL HISTOGRAMS
2721 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2722 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2723 Lock_Action();
2724 {
2725 dHistMap_P[locPID]->Fill(locP);
2726 dHistMap_Phi[locPID]->Fill(locPhi);
2727 dHistMap_Theta[locPID]->Fill(locTheta);
2728 dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
2729 dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
2730 dHistMap_BetaVsP[locPID]->Fill(locP, locBeta_Timing);
2731 dHistMap_DeltaBetaVsP[locPID]->Fill(locP, locDeltaBeta);
2732 dHistMap_VertexZ[locPID]->Fill(locNeutralParticleHypothesis->position().Z());
2733 dHistMap_VertexYVsX[locPID]->Fill(locNeutralParticleHypothesis->position().X(), locNeutralParticleHypothesis->position().Y());
2734 dHistMap_VertexT[locPID]->Fill(locNeutralParticleHypothesis->time());
2735 }
2736 Unlock_Action();
2737 }
2738 return true;
2739}
2740
2741void DHistogramAction_NumReconstructedObjects::Initialize(JEventLoop* locEventLoop)
2742{
2743 string locHistName;
2744
2745 bool locIsRESTEvent = (string(locEventLoop->GetJEvent().GetJEventSource()->className()) == string("DEventSourceREST"));
2746
2747 //CREATE THE HISTOGRAMS
2748 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2749 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2750 {
2751 CreateAndChangeTo_ActionDirectory();
2752
2753 //2D Summary
2754 locHistName = "NumHighLevelObjects";
2755 if(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()) != NULL__null) //already created by another thread, or directory name is duplicate (e.g. two identical steps)
2756 dHist_NumHighLevelObjects = static_cast<TH2D*>(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()));
2757 else
2758 {
2759 dHist_NumHighLevelObjects = new TH2D(locHistName.c_str(), ";;# Objects / Event", 13, 0.5, 13.5, dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2760 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(1, "DRFTime");
2761 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(2, "DSCHit");
2762 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(3, "DTOFPoint");
2763 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(4, "DBCALShower");
2764 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(5, "DFCALShower");
2765 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(6, "DTimeBasedTrack");
2766 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(7, "TrackSCMatches");
2767 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(8, "TrackTOFMatches");
2768 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(9, "TrackBCALMatches");
2769 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(10, "TrackFCALMatches");
2770 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(11, "DBeamPhoton");
2771 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(12, "DChargedTrack");
2772 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(13, "DNeutralShower");
2773 }
2774
2775 //Charged
2776 locHistName = "NumChargedTracks";
2777 dHist_NumChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2778 locHistName = "NumPosChargedTracks";
2779 dHist_NumPosChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2780 locHistName = "NumNegChargedTracks";
2781 dHist_NumNegChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2782
2783 //TBT
2784 locHistName = "NumTimeBasedTracks";
2785 dHist_NumTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2786 locHistName = "NumPosTimeBasedTracks";
2787 dHist_NumPosTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2788 locHistName = "NumNegTimeBasedTracks";
2789 dHist_NumNegTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2790
2791 if(!locIsRESTEvent)
2792 {
2793 //WBT
2794 locHistName = "NumWireBasedTracks";
2795 dHist_NumWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2796 locHistName = "NumPosWireBasedTracks";
2797 dHist_NumPosWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2798 locHistName = "NumNegWireBasedTracks";
2799 dHist_NumNegWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2800
2801 //Track Candidates
2802 locHistName = "NumTrackCandidates";
2803 dHist_NumTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2804 locHistName = "NumPosTrackCandidates";
2805 dHist_NumPosTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2806 locHistName = "NumNegTrackCandidates";
2807 dHist_NumNegTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2808
2809 //CDC Track Candidates
2810 locHistName = "NumPosTrackCandidates_CDC";
2811 dHist_NumPosTrackCandidates_CDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} CDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2812 locHistName = "NumNegTrackCandidates_CDC";
2813 dHist_NumNegTrackCandidates_CDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} CDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2814
2815 //FDC Track Candidates
2816 locHistName = "NumPosTrackCandidates_FDC";
2817 dHist_NumPosTrackCandidates_FDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} FDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2818 locHistName = "NumNegTrackCandidates_FDC";
2819 dHist_NumNegTrackCandidates_FDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} FDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2820 }
2821
2822 //Beam Photons
2823 locHistName = "NumBeamPhotons";
2824 dHist_NumBeamPhotons = GetOrCreate_Histogram<TH1D>(locHistName, ";# DBeamPhoton", dMaxNumBeamPhotons + 1, -0.5, (float)dMaxNumBeamPhotons + 0.5);
2825
2826 //Showers / Neutrals / TOF / SC
2827 locHistName = "NumFCALShowers";
2828 dHist_NumFCALShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DFCALShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2829 locHistName = "NumBCALShowers";
2830 dHist_NumBCALShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DBCALShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2831 locHistName = "NumNeutralShowers";
2832 dHist_NumNeutralShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DNeutralShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2833 locHistName = "NumTOFPoints";
2834 dHist_NumTOFPoints = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTOFPoint", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2835 locHistName = "NumSCHits";
2836 dHist_NumSCHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DSCHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2837
2838 if(!locIsRESTEvent)
2839 {
2840 locHistName = "NumTAGMHits";
2841 dHist_NumTAGMHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTAGMHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2842 locHistName = "NumTAGHHits";
2843 dHist_NumTAGHHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTAGHHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2844 }
2845
2846 //Matches
2847 locHistName = "NumTrackBCALMatches";
2848 dHist_NumTrackBCALMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-BCAL Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
2849 locHistName = "NumTrackFCALMatches";
2850 dHist_NumTrackFCALMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-FCAL Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
2851 locHistName = "NumTrackTOFMatches";
2852 dHist_NumTrackTOFMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-TOF Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
2853 locHistName = "NumTrackSCMatches";
2854 dHist_NumTrackSCMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-SC Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
2855
2856 if(!locIsRESTEvent)
2857 {
2858 //Hits
2859 locHistName = "NumCDCHits";
2860 dHist_NumCDCHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DCDCHit", dMaxNumCDCHits + 1, -0.5, (float)dMaxNumCDCHits + 0.5);
2861 locHistName = "NumFDCWireHits";
2862 dHist_NumFDCWireHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# Wire DFDCHit", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
2863 locHistName = "NumFDCCathodeHits";
2864 dHist_NumFDCCathodeHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# Cathode DFDCHit", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
2865
2866 locHistName = "NumTOFHits";
2867 dHist_NumTOFHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DTOFHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
2868 locHistName = "NumBCALHits";
2869 dHist_NumBCALHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DBCALHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
2870 locHistName = "NumFCALHits";
2871 dHist_NumFCALHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DFCALHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
2872
2873 locHistName = "NumRFSignals";
2874 dHist_NumRFSignals = GetOrCreate_Histogram<TH1I>(locHistName, ";# DRFDigiTime + # DRFTDCDigiTime", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
2875 }
2876
2877 //Return to the base directory
2878 ChangeTo_BaseDirectory();
2879 }
2880 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2881}
2882
2883bool DHistogramAction_NumReconstructedObjects::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2884{
2885 if(Get_NumPreviousParticleCombos() != 0)
2886 return true; //else double-counting!
2887
2888 bool locIsRESTEvent = (string(locEventLoop->GetJEvent().GetJEventSource()->className()) == string("DEventSourceREST"));
2889
2890 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
2891 locEventLoop->Get(locTrackTimeBasedVector);
2892
2893 vector<const DBeamPhoton*> locBeamPhotons;
2894 locEventLoop->Get(locBeamPhotons);
2895
2896 vector<const DFCALShower*> locFCALShowers;
2897 locEventLoop->Get(locFCALShowers);
2898
2899 vector<const DChargedTrack*> locChargedTracks;
2900 locEventLoop->Get(locChargedTracks);
2901
2902 vector<const DBCALShower*> locBCALShowers;
2903 locEventLoop->Get(locBCALShowers);
2904
2905 vector<const DNeutralShower*> locNeutralShowers;
2906 locEventLoop->Get(locNeutralShowers);
2907
2908 vector<const DTOFPoint*> locTOFPoints;
2909 locEventLoop->Get(locTOFPoints);
2910
2911 vector<const DSCHit*> locSCHits;
2912 locEventLoop->Get(locSCHits);
2913
2914 vector<const DRFTime*> locRFTimes;
2915 locEventLoop->Get(locRFTimes);
2916
2917 const DDetectorMatches* locDetectorMatches = NULL__null;
2918 locEventLoop->GetSingle(locDetectorMatches);
2919
2920 //if not REST
2921 vector<const DTrackWireBased*> locTrackWireBasedVector;
2922 vector<const DTrackCandidate*> locTrackCandidates;
2923 vector<const DTrackCandidate*> locTrackCandidates_CDC;
2924 vector<const DTrackCandidate*> locTrackCandidates_FDC;
2925 vector<const DCDCHit*> locCDCHits;
2926 vector<const DFDCHit*> locFDCHits;
2927 vector<const DTOFHit*> locTOFHits;
2928 vector<const DBCALHit*> locBCALHits;
2929 vector<const DFCALHit*> locFCALHits;
2930 vector<const DTAGMHit*> locTAGMHits;
2931 vector<const DTAGHHit*> locTAGHHits;
2932 vector<const DRFDigiTime*> locRFDigiTimes;
2933 vector<const DRFTDCDigiTime*> locRFTDCDigiTimes;
2934
2935 size_t locNumFDCWireHits = 0, locNumFDCCathodeHits = 0;
2936 if(!locIsRESTEvent)
2937 {
2938 locEventLoop->Get(locTrackWireBasedVector);
2939 locEventLoop->Get(locTrackCandidates);
2940 locEventLoop->Get(locTrackCandidates_CDC, "CDC");
2941 locEventLoop->Get(locTrackCandidates_FDC, "FDCCathodes");
2942 locEventLoop->Get(locCDCHits);
2943 locEventLoop->Get(locFDCHits);
2944 locEventLoop->Get(locTOFHits);
2945 locEventLoop->Get(locBCALHits);
2946 locEventLoop->Get(locFCALHits);
2947 locEventLoop->Get(locTAGHHits);
2948 locEventLoop->Get(locTAGMHits);
2949 locEventLoop->Get(locRFDigiTimes);
2950 locEventLoop->Get(locRFTDCDigiTimes);
2951
2952 for(size_t loc_i = 0; loc_i < locFDCHits.size(); ++loc_i)
2953 {
2954 if(locFDCHits[loc_i]->type == DFDCHit::AnodeWire)
2955 ++locNumFDCWireHits;
2956 else
2957 ++locNumFDCCathodeHits;
2958 }
2959 }
2960
2961 //FILL HISTOGRAMS
2962 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2963 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2964 Lock_Action(); //ACQUIRE ROOT LOCK!!
2965 {
2966 //# High-Level Objects
2967 dHist_NumHighLevelObjects->Fill(1, (Double_t)locRFTimes.size());
2968 dHist_NumHighLevelObjects->Fill(2, (Double_t)locSCHits.size());
2969 dHist_NumHighLevelObjects->Fill(3, (Double_t)locTOFPoints.size());
2970 dHist_NumHighLevelObjects->Fill(4, (Double_t)locBCALShowers.size());
2971 dHist_NumHighLevelObjects->Fill(5, (Double_t)locFCALShowers.size());
2972 dHist_NumHighLevelObjects->Fill(6, (Double_t)locTrackTimeBasedVector.size());
2973 dHist_NumHighLevelObjects->Fill(7, (Double_t)locDetectorMatches->Get_NumTrackSCMatches());
2974 dHist_NumHighLevelObjects->Fill(8, (Double_t)locDetectorMatches->Get_NumTrackTOFMatches());
2975 dHist_NumHighLevelObjects->Fill(9, (Double_t)locDetectorMatches->Get_NumTrackBCALMatches());
2976 dHist_NumHighLevelObjects->Fill(10, (Double_t)locDetectorMatches->Get_NumTrackFCALMatches());
2977 dHist_NumHighLevelObjects->Fill(11, (Double_t)locBeamPhotons.size());
2978 dHist_NumHighLevelObjects->Fill(12, (Double_t)locChargedTracks.size());
2979 dHist_NumHighLevelObjects->Fill(13, (Double_t)locNeutralShowers.size());
2980
2981 //Charged
2982 unsigned int locNumPos = 0, locNumNeg = 0;
2983 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2984 {
2985 if(ParticleCharge(locChargedTracks[loc_i]->Get_BestFOM()->PID()) > 0)
2986 ++locNumPos;
2987 else
2988 ++locNumNeg;
2989 }
2990 dHist_NumChargedTracks->Fill(locChargedTracks.size());
2991 dHist_NumPosChargedTracks->Fill(locNumPos);
2992 dHist_NumNegChargedTracks->Fill(locNumNeg);
2993
2994 //TBT
2995 locNumPos = 0; locNumNeg = 0;
2996 for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
2997 {
2998 if(ParticleCharge(locTrackTimeBasedVector[loc_i]->PID()) > 0)
2999 ++locNumPos;
3000 else
3001 ++locNumNeg;
3002 }
3003 dHist_NumTimeBasedTracks->Fill(locTrackTimeBasedVector.size());
3004 dHist_NumPosTimeBasedTracks->Fill(locNumPos);
3005 dHist_NumNegTimeBasedTracks->Fill(locNumNeg);
3006
3007 if(!locIsRESTEvent)
3008 {
3009 //WBT
3010 locNumPos = 0; locNumNeg = 0;
3011 for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
3012 {
3013 if(ParticleCharge(locTrackWireBasedVector[loc_i]->PID()) > 0)
3014 ++locNumPos;
3015 else
3016 ++locNumNeg;
3017 }
3018 dHist_NumWireBasedTracks->Fill(locTrackWireBasedVector.size());
3019 dHist_NumPosWireBasedTracks->Fill(locNumPos);
3020 dHist_NumNegWireBasedTracks->Fill(locNumNeg);
3021
3022 //Candidates
3023 locNumPos = 0; locNumNeg = 0;
3024 for(size_t loc_i = 0; loc_i < locTrackCandidates.size(); ++loc_i)
3025 {
3026 if(locTrackCandidates[loc_i]->charge() > 0.0)
3027 ++locNumPos;
3028 else
3029 ++locNumNeg;
3030 }
3031 dHist_NumTrackCandidates->Fill(locTrackCandidates.size());
3032 dHist_NumPosTrackCandidates->Fill(locNumPos);
3033 dHist_NumNegTrackCandidates->Fill(locNumNeg);
3034
3035 //CDC Candidates
3036 locNumPos = 0; locNumNeg = 0;
3037 for(size_t loc_i = 0; loc_i < locTrackCandidates_CDC.size(); ++loc_i)
3038 {
3039 if(locTrackCandidates_CDC[loc_i]->charge() > 0.0)
3040 ++locNumPos;
3041 else
3042 ++locNumNeg;
3043 }
3044 dHist_NumPosTrackCandidates_CDC->Fill(locNumPos);
3045 dHist_NumNegTrackCandidates_CDC->Fill(locNumNeg);
3046
3047 //FDC Candidates
3048 locNumPos = 0; locNumNeg = 0;
3049 for(size_t loc_i = 0; loc_i < locTrackCandidates_FDC.size(); ++loc_i)
3050 {
3051 if(locTrackCandidates_FDC[loc_i]->charge() > 0.0)
3052 ++locNumPos;
3053 else
3054 ++locNumNeg;
3055 }
3056 dHist_NumPosTrackCandidates_FDC->Fill(locNumPos);
3057 dHist_NumNegTrackCandidates_FDC->Fill(locNumNeg);
3058 }
3059
3060 //Beam Photons
3061 dHist_NumBeamPhotons->Fill((Double_t)locBeamPhotons.size());
3062
3063 //Showers
3064 dHist_NumFCALShowers->Fill((Double_t)locFCALShowers.size());
3065 dHist_NumBCALShowers->Fill((Double_t)locBCALShowers.size());
3066 dHist_NumNeutralShowers->Fill((Double_t)locNeutralShowers.size());
3067
3068 //TOF & SC
3069 dHist_NumTOFPoints->Fill((Double_t)locTOFPoints.size());
3070 dHist_NumSCHits->Fill((Double_t)locSCHits.size());
3071
3072 //TAGGER
3073 if(!locIsRESTEvent)
3074 {
3075 dHist_NumTAGMHits->Fill((Double_t)locTAGMHits.size());
3076 dHist_NumTAGHHits->Fill((Double_t)locTAGHHits.size());
3077 }
3078
3079 //Matches
3080 dHist_NumTrackBCALMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackBCALMatches());
3081 dHist_NumTrackFCALMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackFCALMatches());
3082 dHist_NumTrackTOFMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackTOFMatches());
3083 dHist_NumTrackSCMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackSCMatches());
3084
3085 //Hits
3086 if(!locIsRESTEvent)
3087 {
3088 dHist_NumCDCHits->Fill((Double_t)locCDCHits.size());
3089 dHist_NumFDCWireHits->Fill((Double_t)locNumFDCWireHits);
3090 dHist_NumFDCCathodeHits->Fill((Double_t)locNumFDCCathodeHits);
3091 dHist_NumTOFHits->Fill((Double_t)locTOFHits.size());
3092 dHist_NumBCALHits->Fill((Double_t)locBCALHits.size());
3093 dHist_NumFCALHits->Fill((Double_t)locFCALHits.size());
3094 dHist_NumRFSignals->Fill((Double_t)(locRFDigiTimes.size() + locRFTDCDigiTimes.size()));
3095 }
3096 }
3097 Unlock_Action(); //RELEASE ROOT LOCK!!
3098
3099 return true;
3100}
3101
3102void DHistogramAction_TrackMultiplicity::Initialize(JEventLoop* locEventLoop)
3103{
3104 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
3105 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
3106 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
3107 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
3108 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
3109
3110 //CREATE THE HISTOGRAMS
3111 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
3112 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
3113 {
3114 //So: Default tag is "", User can set it to something else
3115 //In here, if tag is "", get from gparms, if not, leave it alone
3116 //If gparms value does not exist, set it to (and use) "PreSelect"
3117 if(dTrackSelectionTag == "NotATag")
3118 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
3119 if(dShowerSelectionTag == "NotATag")
3120 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
3121
3122 CreateAndChangeTo_ActionDirectory();
3123
3124 string locHistName("NumReconstructedParticles");
3125 if(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()) != NULL__null) //already created by another thread, or directory name is duplicate (e.g. two identical steps)
3126 dHist_NumReconstructedParticles = static_cast<TH2D*>(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()));
3127 else
3128 {
3129 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);
3130 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(1, "# Total");
3131 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(2, "# q != 0");
3132 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(3, "# q = 0");
3133 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(4, "# q = +");
3134 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(5, "# q = -");
3135 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3136 {
3137 string locLabelName = string("# ") + string(ParticleName_ROOT(dFinalStatePIDs[loc_i]));
3138 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(6 + loc_i, locLabelName.c_str());
3139 }
3140 }
3141
3142 locHistName = "NumGoodReconstructedParticles";
3143 if(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()) != NULL__null) //already created by another thread, or directory name is duplicate (e.g. two identical steps)
3144 dHist_NumGoodReconstructedParticles = static_cast<TH2D*>(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()));
3145 else
3146 {
3147 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);
3148 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(1, "# Total");
3149 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(2, "# q != 0");
3150 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(3, "# q = 0");
3151 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(4, "# q = +");
3152 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(5, "# q = -");
3153 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3154 {
3155 string locLabelName = string("# ") + string(ParticleName_ROOT(dFinalStatePIDs[loc_i]));
3156 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(6 + loc_i, locLabelName.c_str());
3157 }
3158 }
3159
3160 //Return to the base directory
3161 ChangeTo_BaseDirectory();
3162 }
3163 japp->RootUnLock(); //RELEASE ROOT LOCK!!
3164}
3165
3166bool DHistogramAction_TrackMultiplicity::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3167{
3168 if(Get_NumPreviousParticleCombos() != 0)
3169 return true; //else double-counting!
3170
3171 vector<const DChargedTrack*> locChargedTracks;
3172 locEventLoop->Get(locChargedTracks);
3173
3174 vector<const DChargedTrack*> locGoodChargedTracks;
3175 locEventLoop->Get(locGoodChargedTracks, dTrackSelectionTag.c_str());
3176
3177 // get #tracks by PID/q type
3178 size_t locNumPositiveTracks = 0, locNumNegativeTracks = 0;
3179 map<Particle_t, size_t> locNumTracksByPID;
3180 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
3181 {
3182 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
3183 Particle_t locPID = locChargedTrackHypothesis->PID();
3184
3185 if(locChargedTrackHypothesis->charge() > 0.0)
3186 ++locNumPositiveTracks;
3187 else
3188 ++locNumNegativeTracks;
3189
3190 if(locNumTracksByPID.find(locPID) != locNumTracksByPID.end())
3191 ++locNumTracksByPID[locPID];
3192 else
3193 locNumTracksByPID[locPID] = 1;
3194 }
3195
3196 // get # good tracks by PID/q type
3197 size_t locNumGoodPositiveTracks = 0, locNumGoodNegativeTracks = 0;
3198 map<Particle_t, size_t> locNumGoodTracksByPID;
3199 for(size_t loc_i = 0; loc_i < locGoodChargedTracks.size(); ++loc_i)
3200 {
3201 const DChargedTrackHypothesis* locChargedTrackHypothesis = locGoodChargedTracks[loc_i]->Get_BestFOM();
3202 Particle_t locPID = locChargedTrackHypothesis->PID();
3203
3204 double locPIDFOM = locChargedTrackHypothesis->dFOM;
3205
3206 if(locChargedTrackHypothesis->charge() > 0.0)
3207 ++locNumGoodPositiveTracks;
3208 else
3209 ++locNumGoodNegativeTracks;
3210
3211 if(locPIDFOM < dMinPIDFOM)
3212 continue;
3213
3214 if(locNumGoodTracksByPID.find(locPID) != locNumGoodTracksByPID.end())
3215 ++locNumGoodTracksByPID[locPID];
3216 else
3217 locNumGoodTracksByPID[locPID] = 1;
3218 }
3219
3220 vector<const DNeutralParticle*> locNeutralParticles;
3221 locEventLoop->Get(locNeutralParticles);
3222
3223 vector<const DNeutralParticle*> locGoodNeutralParticles;
3224 locEventLoop->Get(locGoodNeutralParticles, dShowerSelectionTag.c_str());
3225
3226 // neutrals by pid
3227 for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
3228 {
3229 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
3230 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
3231 continue;
3232
3233 Particle_t locPID = locNeutralParticleHypothesis->PID();
3234 if(locNumTracksByPID.find(locPID) != locNumTracksByPID.end())
3235 ++locNumTracksByPID[locPID];
3236 else
3237 locNumTracksByPID[locPID] = 1;
3238 }
3239
3240 // good neutrals
3241 for(size_t loc_i = 0; loc_i < locGoodNeutralParticles.size(); ++loc_i)
3242 {
3243 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locGoodNeutralParticles[loc_i]->Get_BestFOM();
3244 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
3245 continue;
3246
3247 Particle_t locPID = locNeutralParticleHypothesis->PID();
3248 if(locNumGoodTracksByPID.find(locPID) != locNumGoodTracksByPID.end())
3249 ++locNumGoodTracksByPID[locPID];
3250 else
3251 locNumGoodTracksByPID[locPID] = 1;
3252 }
3253
3254 size_t locNumGoodTracks = locNumGoodPositiveTracks + locNumGoodNegativeTracks;
3255
3256 //FILL HISTOGRAMS
3257 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3258 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3259 Lock_Action();
3260 {
3261 dHist_NumReconstructedParticles->Fill(0.0, (Double_t)(locChargedTracks.size() + locNeutralParticles.size()));
3262 dHist_NumReconstructedParticles->Fill(1.0, (Double_t)locChargedTracks.size());
3263 dHist_NumReconstructedParticles->Fill(2.0, (Double_t)locNeutralParticles.size());
3264 dHist_NumReconstructedParticles->Fill(3.0, (Double_t)locNumPositiveTracks);
3265 dHist_NumReconstructedParticles->Fill(4.0, (Double_t)locNumNegativeTracks);
3266 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3267 dHist_NumReconstructedParticles->Fill(5.0 + (Double_t)loc_i, (Double_t)locNumTracksByPID[dFinalStatePIDs[loc_i]]);
3268
3269 dHist_NumGoodReconstructedParticles->Fill(0.0, (Double_t)(locNumGoodTracks + locGoodNeutralParticles.size()));
3270 dHist_NumGoodReconstructedParticles->Fill(1.0, (Double_t)locNumGoodTracks);
3271 dHist_NumGoodReconstructedParticles->Fill(2.0, (Double_t)locGoodNeutralParticles.size());
3272 dHist_NumGoodReconstructedParticles->Fill(3.0, (Double_t)locNumGoodPositiveTracks);
3273 dHist_NumGoodReconstructedParticles->Fill(4.0, (Double_t)locNumGoodNegativeTracks);
3274 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3275 dHist_NumGoodReconstructedParticles->Fill(5.0 + (Double_t)loc_i, (Double_t)locNumGoodTracksByPID[dFinalStatePIDs[loc_i]]);
3276 }
3277 Unlock_Action();
3278
3279 return true;
3280}
3281