Bug Summary

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

Annotated Source Code

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