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