Bug Summary

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