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*, double> > locHorizontalTOFPaddleTrackDistanceMap; //double: delta-y
1145 map<const DKinematicData*, pair<const DTOFPaddleHit*, double> > locVerticalTOFPaddleTrackDistanceMap; //double: delta-x
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;
1160 double locStartTime = locParticleID->Calc_PropagatedRFTime(locKinematicData, locEventRFBunch);
1161 const DTOFPaddleHit* locClosestTOFPaddleHit_Vertical = locParticleID->Get_ClosestTOFPaddleHit_Vertical(locReferenceTrajectory, locTOFPaddleHits, locStartTime, locBestDeltaX);
1162 const DTOFPaddleHit* locClosestTOFPaddleHit_Horizontal = locParticleID->Get_ClosestTOFPaddleHit_Horizontal(locReferenceTrajectory, locTOFPaddleHits, locStartTime, locBestDeltaY);
1163 if(locClosestTOFPaddleHit_Vertical != NULL__null)
1164 locVerticalTOFPaddleTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPaddleHit*, double>(locClosestTOFPaddleHit_Vertical, locBestDeltaX);
1165 if(locClosestTOFPaddleHit_Horizontal != NULL__null)
1166 locHorizontalTOFPaddleTrackDistanceMap[locTrackIterator->second] = pair<const DTOFPaddleHit*, double>(locClosestTOFPaddleHit_Horizontal, locBestDeltaY);
1167 }
1168
1169 //TRACK / SC CLOSEST MATCHES
1170 map<const DKinematicData*, double> locSCTrackDistanceMap;
1171 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1172 {
1173 double locBestDeltaPhi = 999.0;
1174 if(locParticleID->Get_ClosestToTrack_SC(locTrackIterator->second, locSCHits, locBestDeltaPhi) != NULL__null)
1175 locSCTrackDistanceMap[locTrackIterator->second] = locBestDeltaPhi;
1176 }
1177
1178 //PROJECTED HIT POSITIONS
1179 map<const DKinematicData*, pair<int, bool> > dProjectedSCPaddleMap; //pair: paddle, hit-barrel-flag (false if bend/nose)
1180 map<const DKinematicData*, pair<int, int> > dProjectedTOF2DPaddlesMap; //pair: vertical, horizontal
1181 map<const DKinematicData*, pair<float, float> > dProjectedTOFXYMap; //pair: x, y
1182 map<const DKinematicData*, pair<int, int> > dProjectedFCALRowColumnMap; //pair: column, row
1183 map<const DKinematicData*, pair<float, float> > dProjectedFCALXYMap; //pair: x, y
1184 map<const DKinematicData*, pair<float, int> > dProjectedBCALModuleSectorMap; //pair: z, module
1185 map<const DKinematicData*, pair<float, float> > dProjectedBCALPhiZMap; //pair: z, phi
1186 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1187 {
1188 const DKinematicData* locTrack = locTrackIterator->second;
1189 const DReferenceTrajectory* locReferenceTrajectory = NULL__null;
1190 if(locIsTimeBased)
1191 locReferenceTrajectory = (static_cast<const DTrackTimeBased*>(locTrack))->rt;
1192 else
1193 locReferenceTrajectory = (static_cast<const DTrackWireBased*>(locTrack))->rt;
1194 if(locReferenceTrajectory == NULL__null)
1195 break; //e.g. REST data: no trajectory
1196
1197 //SC
1198 DVector3 locSCIntersection;
1199 bool locProjBarrelFlag = false;
1200 unsigned int locProjectedSCPaddle = locParticleID->PredictSCSector(locReferenceTrajectory, 999.9, &locSCIntersection, &locProjBarrelFlag);
1201 if(locProjectedSCPaddle != 0)
1202 dProjectedSCPaddleMap[locTrack] = pair<int, bool>(locProjectedSCPaddle, locProjBarrelFlag);
1203
1204 //TOF
1205 DVector3 locTOFIntersection;
1206 unsigned int locHorizontalBar = 0, locVerticalBar = 0;
1207 if(locParticleID->PredictTOFPaddles(locReferenceTrajectory, locHorizontalBar, locVerticalBar, &locTOFIntersection))
1208 {
1209 dProjectedTOF2DPaddlesMap[locTrack] = pair<int, int>(locVerticalBar, locHorizontalBar);
1210 dProjectedTOFXYMap[locTrack] = pair<float, float>(locTOFIntersection.X(), locTOFIntersection.Y());
1211 }
1212
1213 //FCAL
1214 DVector3 locFCALIntersection;
1215 unsigned int locRow = 0, locColumn = 0;
1216 if(locParticleID->PredictFCALHit(locReferenceTrajectory, locRow, locColumn, &locFCALIntersection))
1217 {
1218 dProjectedFCALRowColumnMap[locTrack] = pair<int, int>(locColumn, locRow);
1219 dProjectedFCALXYMap[locTrack] = pair<float, float>(locFCALIntersection.X(), locFCALIntersection.Y());
1220 }
1221
1222 //BCAL
1223 DVector3 locBCALIntersection;
1224 unsigned int locModule = 0, locSector = 0;
1225 if(locParticleID->PredictBCALWedge(locReferenceTrajectory, locModule, locSector, &locBCALIntersection))
1226 {
1227 dProjectedBCALModuleSectorMap[locTrack] = pair<float, int>(locBCALIntersection.Z(), locModule);
1228 dProjectedBCALPhiZMap[locTrack] = pair<float, float>(locBCALIntersection.Z(), locBCALIntersection.Phi()*180.0/TMath::Pi());
1229 }
1230 }
1231
1232 //FILL HISTOGRAMS
1233 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1234 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1235 Lock_Action(); //ACQUIRE ROOT LOCK!!
1236 {
1237 /********************************************************** MATCHING DISTANCE **********************************************************/
1238
1239 //BCAL
1240 map<const DKinematicData*, pair<double, double> >::iterator locBCALIterator = locBCALTrackDistanceMap.begin();
1241 for(; locBCALIterator != locBCALTrackDistanceMap.end(); ++locBCALIterator)
1242 {
1243 const DKinematicData* locTrack = locBCALIterator->first;
1244 dHistMap_BCALDeltaPhiVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locBCALIterator->second.first*180.0/TMath::Pi());
1245 dHistMap_BCALDeltaZVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locBCALIterator->second.second);
1246 }
1247
1248 //FCAL
1249 map<const DKinematicData*, double>::iterator locFCALIterator = locFCALTrackDistanceMap.begin();
1250 for(; locFCALIterator != locFCALTrackDistanceMap.end(); ++locFCALIterator)
1251 {
1252 const DKinematicData* locTrack = locFCALIterator->first;
1253 dHistMap_FCALTrackDistanceVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locFCALIterator->second);
1254 dHistMap_FCALTrackDistanceVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locFCALIterator->second);
1255 }
1256
1257 //TOF Paddle
1258 //Horizontal
1259 map<const DKinematicData*, pair<const DTOFPaddleHit*, double> >::iterator locTOFPaddleIterator = locHorizontalTOFPaddleTrackDistanceMap.begin();
1260 for(; locTOFPaddleIterator != locHorizontalTOFPaddleTrackDistanceMap.end(); ++locTOFPaddleIterator)
1261 {
1262 double locDistance = locTOFPaddleIterator->second.second;
1263 dHistMap_TOFPaddleTrackDeltaY[locIsTimeBased]->Fill(locDistance);
1264 }
1265 //Vertical
1266 locTOFPaddleIterator = locVerticalTOFPaddleTrackDistanceMap.begin();
1267 for(; locTOFPaddleIterator != locVerticalTOFPaddleTrackDistanceMap.end(); ++locTOFPaddleIterator)
1268 {
1269 double locDistance = locTOFPaddleIterator->second.second;
1270 dHistMap_TOFPaddleTrackDeltaX[locIsTimeBased]->Fill(locDistance);
1271 }
1272
1273 //TOF Point
1274 map<const DKinematicData*, pair<const DTOFPoint*, pair<double, double> > >::iterator locTOFPointIterator = locTOFPointTrackDistanceMap.begin();
1275 for(; locTOFPointIterator != locTOFPointTrackDistanceMap.end(); ++locTOFPointIterator)
1276 {
1277 const DKinematicData* locTrack = locTOFPointIterator->first;
1278 const DTOFPoint* locTOFPoint = locTOFPointIterator->second.first;
1279 double locDeltaX = locTOFPointIterator->second.second.first;
1280 double locDeltaY = locTOFPointIterator->second.second.second;
1281
1282 double locDistance = sqrt(locDeltaX*locDeltaX + locDeltaY*locDeltaY);
1283 if((fabs(locDeltaX) < 500.0) && (fabs(locDeltaY) < 500.0)) //else position not well-defined
1284 {
1285 dHistMap_TOFPointTrackDistanceVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDistance);
1286 dHistMap_TOFPointTrackDistanceVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locDistance);
1287 if((locTOFPoint->dHorizontalBar != 0) && (locTOFPoint->dVerticalBar != 0))
1288 dHistMap_TOFPointTrackDistance_BothPlanes[locIsTimeBased]->Fill(locDistance);
1289 else
1290 dHistMap_TOFPointTrackDistance_OnePlane[locIsTimeBased]->Fill(locDistance);
1291 }
1292
1293 dHistMap_TOFPointTrackDeltaXVsHorizontalPaddle[locIsTimeBased]->Fill(locTOFPoint->dHorizontalBar, locDeltaX);
1294 dHistMap_TOFPointTrackDeltaXVsVerticalPaddle[locIsTimeBased]->Fill(locTOFPoint->dVerticalBar, locDeltaX);
1295
1296 dHistMap_TOFPointTrackDeltaYVsHorizontalPaddle[locIsTimeBased]->Fill(locTOFPoint->dHorizontalBar, locDeltaY);
1297 dHistMap_TOFPointTrackDeltaYVsVerticalPaddle[locIsTimeBased]->Fill(locTOFPoint->dVerticalBar, locDeltaY);
1298 }
1299
1300 //SC
1301 map<const DKinematicData*, double>::iterator locSCIterator = locSCTrackDistanceMap.begin();
1302 if(locSCHits.size() <= 4) //don't fill if every paddle fired!
1303 {
1304 for(; locSCIterator != locSCTrackDistanceMap.end(); ++locSCIterator)
1305 {
1306 const DKinematicData* locTrack = locSCIterator->first;
1307 double locDeltaPhi = locSCIterator->second*180.0/TMath::Pi();
1308 dHistMap_SCTrackDeltaPhiVsP[locIsTimeBased]->Fill(locTrack->momentum().Mag(), locDeltaPhi);
1309 dHistMap_SCTrackDeltaPhiVsTheta[locIsTimeBased]->Fill(locTrack->momentum().Theta()*180.0/TMath::Pi(), locDeltaPhi);
1310 }
1311 }
1312
1313 /********************************************************* MATCHING EFFICINECY *********************************************************/
1314
1315 //Does-it-match, by detector
1316 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1317 {
1318 const DKinematicData* locTrack = locTrackIterator->second;
1319 double locTheta = locTrack->momentum().Theta()*180.0/TMath::Pi();
1320 double locPhi = locTrack->momentum().Phi()*180.0/TMath::Pi();
1321 double locP = locTrack->momentum().Mag();
1322
1323 //BCAL
1324 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_START))
1325 {
1326 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_BCAL))
1327 {
1328 dHistMap_PVsTheta_HasHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locP);
1329 dHistMap_PhiVsTheta_HasHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1330 if(dProjectedBCALModuleSectorMap.find(locTrack) != dProjectedBCALModuleSectorMap.end())
1331 {
1332 pair<float, float>& locPositionPair = dProjectedBCALPhiZMap[locTrack];
1333 dHistMap_TrackBCALPhiVsZ_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1334 pair<float, int>& locElementPair = dProjectedBCALModuleSectorMap[locTrack];
1335 dHistMap_TrackBCALModuleVsZ_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1336 }
1337 }
1338 else
1339 {
1340 dHistMap_PVsTheta_NoHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locP);
1341 dHistMap_PhiVsTheta_NoHit[SYS_BCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1342 if(dProjectedBCALModuleSectorMap.find(locTrack) != dProjectedBCALModuleSectorMap.end())
1343 {
1344 pair<float, float>& locPositionPair = dProjectedBCALPhiZMap[locTrack];
1345 dHistMap_TrackBCALPhiVsZ_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1346 pair<float, int>& locElementPair = dProjectedBCALModuleSectorMap[locTrack];
1347 dHistMap_TrackBCALModuleVsZ_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1348 }
1349 }
1350 }
1351
1352 //FCAL
1353 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1354 {
1355 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1356 {
1357 dHistMap_PVsTheta_HasHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locP);
1358 dHistMap_PhiVsTheta_HasHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1359 if(dProjectedFCALRowColumnMap.find(locTrack) != dProjectedFCALRowColumnMap.end())
1360 {
1361 pair<float, float>& locPositionPair = dProjectedFCALXYMap[locTrack];
1362 dHistMap_TrackFCALYVsX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1363 pair<int, int>& locElementPair = dProjectedFCALRowColumnMap[locTrack];
1364 dHistMap_TrackFCALRowVsColumn_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1365 }
1366 }
1367 else
1368 {
1369 dHistMap_PVsTheta_NoHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locP);
1370 dHistMap_PhiVsTheta_NoHit[SYS_FCAL][locIsTimeBased]->Fill(locTheta, locPhi);
1371 if(dProjectedFCALRowColumnMap.find(locTrack) != dProjectedFCALRowColumnMap.end())
1372 {
1373 pair<float, float>& locPositionPair = dProjectedFCALXYMap[locTrack];
1374 dHistMap_TrackFCALYVsX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1375 pair<int, int>& locElementPair = dProjectedFCALRowColumnMap[locTrack];
1376 dHistMap_TrackFCALRowVsColumn_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1377 }
1378 }
1379 }
1380
1381 //TOF Paddle
1382 if((dProjectedTOFXYMap.find(locTrack) != dProjectedTOFXYMap.end()) && locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1383 {
1384 pair<float, float>& locPositionPair = dProjectedTOFXYMap[locTrack];
1385 pair<int, int>& locPaddlePair = dProjectedTOF2DPaddlesMap[locTrack]; //vertical, horizontal
1386
1387 //Horizontal
1388 if(locHorizontalTOFPaddleTrackDistanceMap.find(locTrack) != locHorizontalTOFPaddleTrackDistanceMap.end())
1389 {
1390 double locDistance = locHorizontalTOFPaddleTrackDistanceMap[locTrack].second;
1391 if(locDistance <= dMinTOFPaddleMatchDistance) //match
1392 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1393 else //no match
1394 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1395 }
1396 else // no match
1397 dHistMap_TOFPaddleHorizontalPaddleVsTrackX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPaddlePair.second);
1398
1399 //Vertical
1400 if(locVerticalTOFPaddleTrackDistanceMap.find(locTrack) != locVerticalTOFPaddleTrackDistanceMap.end())
1401 {
1402 double locDistance = locVerticalTOFPaddleTrackDistanceMap[locTrack].second;
1403 if(locDistance <= dMinTOFPaddleMatchDistance) //match
1404 dHistMap_TOFPaddleTrackYVsVerticalPaddle_HasHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1405 else //no match
1406 dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1407 }
1408 else // no match
1409 dHistMap_TOFPaddleTrackYVsVerticalPaddle_NoHit[locIsTimeBased]->Fill(locPaddlePair.first, locPositionPair.second);
1410 }
1411
1412 //TOF Point
1413 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL))
1414 {
1415 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1416 {
1417 dHistMap_PVsTheta_HasHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locP);
1418 dHistMap_PhiVsTheta_HasHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locPhi);
1419 if(dProjectedTOFXYMap.find(locTrack) != dProjectedTOFXYMap.end())
1420 {
1421 pair<float, float>& locPositionPair = dProjectedTOFXYMap[locTrack];
1422 dHistMap_TrackTOFYVsX_HasHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1423 pair<int, int>& locElementPair = dProjectedTOF2DPaddlesMap[locTrack];
1424 dHistMap_TrackTOF2DPaddles_HasHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1425 }
1426 }
1427 else
1428 {
1429 dHistMap_PVsTheta_NoHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locP);
1430 dHistMap_PhiVsTheta_NoHit[SYS_TOF][locIsTimeBased]->Fill(locTheta, locPhi);
1431 if(dProjectedTOFXYMap.find(locTrack) != dProjectedTOFXYMap.end())
1432 {
1433 pair<float, float>& locPositionPair = dProjectedTOFXYMap[locTrack];
1434 dHistMap_TrackTOFYVsX_NoHit[locIsTimeBased]->Fill(locPositionPair.first, locPositionPair.second);
1435 pair<int, int>& locElementPair = dProjectedTOF2DPaddlesMap[locTrack];
1436 dHistMap_TrackTOF2DPaddles_NoHit[locIsTimeBased]->Fill(locElementPair.first, locElementPair.second);
1437 }
1438 }
1439 }
1440
1441 //SC
1442 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_BCAL) || locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_FCAL) || locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_TOF))
1443 {
1444 if(locDetectorMatches->Get_IsMatchedToDetector(locTrack, SYS_START))
1445 {
1446 dHistMap_PVsTheta_HasHit[SYS_START][locIsTimeBased]->Fill(locTheta, locP);
1447 dHistMap_PhiVsTheta_HasHit[SYS_START][locIsTimeBased]->Fill(locTheta, locPhi);
1448 if(dProjectedSCPaddleMap.find(locTrack) != dProjectedSCPaddleMap.end())
1449 {
1450 dHistMap_SCPaddleVsTheta_HasHit[locIsTimeBased]->Fill(locTheta, dProjectedSCPaddleMap[locTrack].first);
1451 if(dProjectedSCPaddleMap[locTrack].second)
1452 dHistMap_SCPaddle_BarrelRegion_HasHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1453 else
1454 dHistMap_SCPaddle_NoseRegion_HasHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1455 }
1456 }
1457 else
1458 {
1459 dHistMap_PVsTheta_NoHit[SYS_START][locIsTimeBased]->Fill(locTheta, locP);
1460 dHistMap_PhiVsTheta_NoHit[SYS_START][locIsTimeBased]->Fill(locTheta, locPhi);
1461 if(dProjectedSCPaddleMap.find(locTrack) != dProjectedSCPaddleMap.end())
1462 {
1463 dHistMap_SCPaddleVsTheta_NoHit[locIsTimeBased]->Fill(locTheta, dProjectedSCPaddleMap[locTrack].first);
1464 if(dProjectedSCPaddleMap[locTrack].second)
1465 dHistMap_SCPaddle_BarrelRegion_NoHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1466 else
1467 dHistMap_SCPaddle_NoseRegion_NoHit[locIsTimeBased]->Fill(dProjectedSCPaddleMap[locTrack].first);
1468 }
1469 }
1470 }
1471 }
1472 //Is-Matched to Something
1473 for(locTrackIterator = locBestTrackMap.begin(); locTrackIterator != locBestTrackMap.end(); ++locTrackIterator)
1474 {
1475 const DKinematicData* locTrack = locTrackIterator->second;
1476 double locTheta = locTrack->momentum().Theta()*180.0/TMath::Pi();
1477 double locP = locTrack->momentum().Mag();
1478 if(locDetectorMatches->Get_IsMatchedToHit(locTrack))
1479 dHistMap_TrackPVsTheta_HitMatch[locIsTimeBased]->Fill(locTheta, locP);
1480 else
1481 dHistMap_TrackPVsTheta_NoHitMatch[locIsTimeBased]->Fill(locTheta, locP);
1482 }
1483 }
1484 Unlock_Action(); //RELEASE ROOT LOCK!!
1485}
1486
1487void DHistogramAction_DetectorPID::Initialize(JEventLoop* locEventLoop)
1488{
1489 //Create any histograms/trees/etc. within a ROOT lock.
1490 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
1491
1492 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
1493 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
1494
1495 string locHistName, locHistTitle, locParticleROOTName;
1496
1497 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
1498 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
1499 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
1500 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
1501 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
1502
1503 //CREATE THE HISTOGRAMS
1504 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
1505 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
1506 {
1507 //So: Default tag is "", User can set it to something else
1508 //In here, if tag is "", get from gparms, if not, leave it alone
1509 //If gparms value does not exist, set it to (and use) "PreSelect"
1510 if(dTrackSelectionTag == "NotATag")
1511 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
1512 if(dShowerSelectionTag == "NotATag")
1513 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
1514
1515 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
1516 //If another thread has already created the folder, it just changes to it.
1517 CreateAndChangeTo_ActionDirectory();
1518
1519 //q = 0
1520 locParticleROOTName = ParticleName_ROOT(Gamma);
1521 //BCAL
1522 CreateAndChangeTo_Directory("BCAL", "BCAL");
1523
1524 locHistName = "BetaVsP_q0";
1525 locHistTitle = "BCAL q^{0};Shower Energy (GeV);#beta";
1526 dHistMap_BetaVsP[SYS_BCAL][0] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1527
1528 locHistName = "DeltaTVsShowerE_Photon";
1529 locHistTitle = string("BCAL q^{0}") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat_{BCAL - RF}");
1530 dHistMap_DeltaTVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1531
1532/*
1533 //Uncomment when ready!
1534 locHistName = "TimePullVsShowerE_Photon";
1535 locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
1536 dHistMap_TimePullVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1537
1538 locHistName = "TimeFOMVsShowerE_Photon";
1539 locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
1540 dHistMap_TimeFOMVsP[SYS_BCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1541*/
1542 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1543
1544 //FCAL
1545 CreateAndChangeTo_Directory("FCAL", "FCAL");
1546
1547 locHistName = "BetaVsP_q0";
1548 locHistTitle = "FCAL q^{0};Shower Energy (GeV);#beta";
1549 dHistMap_BetaVsP[SYS_FCAL][0] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1550
1551 locHistName = "DeltaTVsShowerE_Photon";
1552 locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat_{FCAL - RF}");
1553 dHistMap_DeltaTVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1554
1555/*
1556 //Uncomment when ready!
1557 locHistName = "TimePullVsShowerE_Photon";
1558 locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");
1559 dHistMap_TimePullVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1560
1561 locHistName = "TimeFOMVsShowerE_Photon";
1562 locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");
1563 dHistMap_TimeFOMVsP[SYS_FCAL][Gamma] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1564*/
1565 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1566
1567 //q +/-
1568 for(int locCharge = -1; locCharge <= 1; locCharge += 2)
1569 {
1570 string locParticleName = (locCharge == -1) ? "q-" : "q+";
1571 string locParticleROOTName = (locCharge == -1) ? "q^{-}" : "q^{+}";
1572
1573 //SC
1574 CreateAndChangeTo_Directory("SC", "SC");
1575
1576 locHistName = string("dEdXVsP_") + locParticleName;
1577 locHistTitle = locParticleROOTName + ";p (GeV/c);SC dE/dX (MeV/cm)";
1578 dHistMap_dEdXVsP[SYS_START][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1579
1580 locHistName = string("BetaVsP_") + locParticleName;
1581 locHistTitle = string("SC ") + locParticleROOTName + string(";p (GeV/c);#beta");
1582 dHistMap_BetaVsP[SYS_START][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1583
1584 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1585
1586 //TOF
1587 CreateAndChangeTo_Directory("TOF", "TOF");
1588
1589 locHistName = string("dEdXVsP_") + locParticleName;
1590 locHistTitle = locParticleROOTName + ";p (GeV/c);TOF dE/dX (MeV/cm)";
1591 dHistMap_dEdXVsP[SYS_TOF][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1592
1593 locHistName = string("BetaVsP_") + locParticleName;
1594 locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);#beta");
1595 dHistMap_BetaVsP[SYS_TOF][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1596
1597 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1598
1599 //BCAL
1600 CreateAndChangeTo_Directory("BCAL", "BCAL");
1601
1602 locHistName = string("BetaVsP_") + locParticleName;
1603 locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);#beta");
1604 dHistMap_BetaVsP[SYS_BCAL][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1605
1606 locHistName = string("EOverPVsP_") + locParticleName;
1607 locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);E_{Shower}/p_{Track} (c);");
1608 dHistMap_BCALEOverPVsP[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1609
1610 locHistName = string("EOverPVsTheta_") + locParticleName;
1611 locHistTitle = string("BCAL ") + locParticleROOTName + string(";#theta#circ;E_{Shower}/p_{Track} (c);");
1612 dHistMap_BCALEOverPVsTheta[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALThetaBins, dMinBCALTheta, dMaxBCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1613
1614 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1615
1616 //FCAL
1617 CreateAndChangeTo_Directory("FCAL", "FCAL");
1618
1619 locHistName = string("BetaVsP_") + locParticleName;
1620 locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#beta");
1621 dHistMap_BetaVsP[SYS_FCAL][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);
1622
1623 locHistName = string("EOverPVsP_") + locParticleName;
1624 locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);E_{Shower}/p_{Track} (c);");
1625 dHistMap_FCALEOverPVsP[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1626
1627 locHistName = string("EOverPVsTheta_") + locParticleName;
1628 locHistTitle = string("FCAL ") + locParticleROOTName + string(";#theta#circ;E_{Shower}/p_{Track} (c);");
1629 dHistMap_FCALEOverPVsTheta[locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DFCALThetaBins, dMinFCALTheta, dMaxFCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);
1630
1631 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1632
1633 //CDC
1634 CreateAndChangeTo_Directory("CDC", "CDC");
1635
1636 locHistName = string("dEdXVsP_") + locParticleName;
1637 locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx (keV/cm)");
1638 dHistMap_dEdXVsP[SYS_CDC][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1639
1640 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1641
1642 //FDC
1643 CreateAndChangeTo_Directory("FDC", "FDC");
1644
1645 locHistName = string("dEdXVsP_") + locParticleName;
1646 locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC dE/dx (keV/cm)");
1647 dHistMap_dEdXVsP[SYS_FDC][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);
1648
1649 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1650 }
1651
1652 //delta's by PID
1653 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
1654 {
1655 Particle_t locPID = dFinalStatePIDs[loc_i];
1656 string locParticleName = ParticleType(locPID);
1657 string locParticleROOTName = ParticleName_ROOT(locPID);
1658
1659 //SC
1660 CreateAndChangeTo_Directory("SC", "SC");
1661
1662 locHistName = string("DeltadEdXVsP_") + locParticleName;
1663 locHistTitle = locParticleROOTName + " Candidates;p (GeV/c);SC #Delta(dE/dX) (MeV/cm)";
1664 dHistMap_DeltadEdXVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1665
1666 locHistName = string("DeltaBetaVsP_") + locParticleName;
1667 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1668 dHistMap_DeltaBetaVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1669
1670 locHistName = string("DeltaTVsP_") + locParticleName;
1671 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{SC - RF}");
1672 dHistMap_DeltaTVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1673
1674/*
1675 //Uncomment when ready!
1676 locHistName = string("TimePullVsP_") + locParticleName;
1677 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1678 dHistMap_TimePullVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1679
1680 locHistName = string("TimeFOMVsP_") + locParticleName;
1681 locHistTitle = string("SC ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1682 dHistMap_TimeFOMVsP[SYS_START][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1683*/
1684 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1685
1686 //TOF
1687 CreateAndChangeTo_Directory("TOF", "TOF");
1688
1689 locHistName = string("DeltadEdXVsP_") + locParticleName;
1690 locHistTitle = locParticleROOTName + " Candidates;p (GeV/c);TOF #Delta(dE/dX) (MeV/cm)";
1691 dHistMap_DeltadEdXVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1692
1693 locHistName = string("DeltaBetaVsP_") + locParticleName;
1694 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1695 dHistMap_DeltaBetaVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1696
1697 locHistName = string("DeltaTVsP_") + locParticleName;
1698 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{TOF - RF}");
1699 dHistMap_DeltaTVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1700
1701/*
1702 //Uncomment when ready!
1703 locHistName = string("TimePullVsP_") + locParticleName;
1704 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1705 dHistMap_TimePullVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1706
1707 locHistName = string("TimeFOMVsP_") + locParticleName;
1708 locHistTitle = string("TOF ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1709 dHistMap_TimeFOMVsP[SYS_TOF][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1710*/
1711 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1712
1713 //BCAL
1714 CreateAndChangeTo_Directory("BCAL", "BCAL");
1715
1716 locHistName = string("DeltaBetaVsP_") + locParticleName;
1717 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1718 dHistMap_DeltaBetaVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1719
1720 locHistName = string("DeltaTVsP_") + locParticleName;
1721 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{BCAL - RF}");
1722 dHistMap_DeltaTVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1723
1724/*
1725 //Uncomment when ready!
1726 locHistName = string("TimePullVsP_") + locParticleName;
1727 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1728 dHistMap_TimePullVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1729
1730 locHistName = string("TimeFOMVsP_") + locParticleName;
1731 locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1732 dHistMap_TimeFOMVsP[SYS_BCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1733*/
1734 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1735
1736 //FCAL
1737 CreateAndChangeTo_Directory("FCAL", "FCAL");
1738
1739 locHistName = string("DeltaBetaVsP_") + locParticleName;
1740 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Delta#beta");
1741 dHistMap_DeltaBetaVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
1742
1743 locHistName = string("DeltaTVsP_") + locParticleName;
1744 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{FCAL - RF}");
1745 dHistMap_DeltaTVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
1746
1747/*
1748 //Uncomment when ready!
1749 locHistName = string("TimePullVsP_") + locParticleName;
1750 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat/#sigma_{#Deltat}");
1751 dHistMap_TimePullVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1752
1753 locHistName = string("TimeFOMVsP_") + locParticleName;
1754 locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;p (GeV/c);Timing PID Confidence Level");
1755 dHistMap_TimeFOMVsP[SYS_FCAL][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1756*/
1757 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1758
1759 //CDC
1760 CreateAndChangeTo_Directory("CDC", "CDC");
1761
1762 locHistName = string("DeltadEdXVsP_") + locParticleName;
1763 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) (keV/cm)");
1764 dHistMap_DeltadEdXVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1765
1766/*
1767 //Uncomment when ready!
1768 locHistName = string("dEdXPullVsP_") + locParticleName;
1769 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
1770 dHistMap_dEdXPullVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1771
1772 locHistName = string("dEdXFOMVsP_") + locParticleName;
1773 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx PID Confidence Level");
1774 dHistMap_dEdXFOMVsP[SYS_CDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1775*/
1776 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1777
1778 //FDC
1779 CreateAndChangeTo_Directory("FDC", "FDC");
1780
1781 locHistName = string("DeltadEdXVsP_") + locParticleName;
1782 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX) (keV/cm)");
1783 dHistMap_DeltadEdXVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);
1784
1785/*
1786 //Uncomment when ready!
1787 locHistName = string("dEdXPullVsP_") + locParticleName;
1788 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");
1789 dHistMap_dEdXPullVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);
1790
1791 locHistName = string("dEdXFOMVsP_") + locParticleName;
1792 locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC dE/dx PID Confidence Level");
1793 dHistMap_dEdXFOMVsP[SYS_FDC][locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);
1794*/
1795 gDirectory(TDirectory::CurrentDirectory())->cd("..");
1796 }
1797
1798 //Return to the base directory
1799 ChangeTo_BaseDirectory();
1800 }
1801 japp->RootUnLock(); //RELEASE ROOT LOCK!!
1802}
1803
1804bool DHistogramAction_DetectorPID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
1805{
1806 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
1807
1808 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
1809 if(Get_NumPreviousParticleCombos() != 0)
1810 return true;
1811
1812 vector<const DChargedTrack*> locChargedTracks;
1813 locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
1814
1815 vector<const DNeutralParticle*> locNeutralParticles;
1816 locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
1817
1818 const DDetectorMatches* locDetectorMatches = NULL__null;
1819 locEventLoop->GetSingle(locDetectorMatches);
1820
1821 const DEventRFBunch* locEventRFBunch = NULL__null;
1822 locEventLoop->GetSingle(locEventRFBunch);
1823
1824 const DParticleID* locParticleID = NULL__null;
1825 locEventLoop->GetSingle(locParticleID);
1826
1827 //FILL HISTOGRAMS
1828 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
1829 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
1830 Lock_Action(); //ACQUIRE ROOT LOCK!!
1831 {
1832 if(locEventRFBunch->dTimeSource != SYS_NULL) //only histogram beta for neutrals if the t0 is well known
1833 {
1834 for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
1835 {
1836 //doesn't matter which hypothesis you use for beta: t0 is from DEventVertex time
1837 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
1838 double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
1839 const DNeutralShower* locNeutralShower = locNeutralParticles[loc_i]->dNeutralShower;
1840 double locShowerEnergy = locNeutralShower->dEnergy;
1841
1842 double locDeltaT = locNeutralParticleHypothesis->time() - locEventRFBunch->dTime;
1843 if(locNeutralShower->dDetectorSystem == SYS_BCAL)
1844 {
1845 dHistMap_BetaVsP[SYS_BCAL][0]->Fill(locShowerEnergy, locBeta_Timing);
1846 dHistMap_DeltaTVsP[SYS_BCAL][Gamma]->Fill(locShowerEnergy, locDeltaT);
1847 }
1848 else
1849 {
1850 dHistMap_BetaVsP[SYS_FCAL][0]->Fill(locShowerEnergy, locBeta_Timing);
1851 dHistMap_DeltaTVsP[SYS_FCAL][Gamma]->Fill(locShowerEnergy, locDeltaT);
1852 }
1853 }
1854 }
1855
1856 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
1857 {
1858 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestTrackingFOM();
1859 int locCharge = ParticleCharge(locChargedTrackHypothesis->PID());
1860 if(dHistMap_dEdXVsP[SYS_START].find(locCharge) == dHistMap_dEdXVsP[SYS_START].end())
1861 continue;
1862
1863 double locStartTime = locParticleID->Calc_PropagatedRFTime(locChargedTrackHypothesis, locEventRFBunch);
1864
1865 const DTrackTimeBased* locTrackTimeBased = NULL__null;
1866 locChargedTrackHypothesis->GetSingle(locTrackTimeBased);
1867
1868 Particle_t locPID = locChargedTrackHypothesis->PID();
1869 double locP = locTrackTimeBased->momentum().Mag();
1870 double locTheta = locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi();
1871
1872 //if RF time is indeterminate, start time will be NaN
1873 const DBCALShowerMatchParams* locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
1874 const DFCALShowerMatchParams* locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
1875 const DTOFHitMatchParams* locTOFHitMatchParams = locChargedTrackHypothesis->Get_TOFHitMatchParams();
1876 const DSCHitMatchParams* locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
1877
1878 if(locSCHitMatchParams != NULL__null)
1879 {
1880 dHistMap_dEdXVsP[SYS_START][locCharge]->Fill(locP, locSCHitMatchParams->dEdx*1.0E3);
1881 if((locEventRFBunch->dTimeSource != SYS_START) && (locEventRFBunch->dNumParticleVotes > 1))
1882 {
1883 //If SC was used for RF time, don't compute delta-beta
1884 double locBeta_Timing = locSCHitMatchParams->dPathLength/(29.9792458*(locSCHitMatchParams->dHitTime - locChargedTrackHypothesis->t0()));
1885 dHistMap_BetaVsP[SYS_START][locCharge]->Fill(locP, locBeta_Timing);
1886 if(dHistMap_DeltaBetaVsP[SYS_START].find(locPID) != dHistMap_DeltaBetaVsP[SYS_START].end())
1887 {
1888 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1889 dHistMap_DeltaBetaVsP[SYS_START][locPID]->Fill(locP, locDeltaBeta);
1890 double locDeltaT = locSCHitMatchParams->dHitTime - locSCHitMatchParams->dFlightTime - locStartTime;
1891 dHistMap_DeltaTVsP[SYS_START][locPID]->Fill(locP, locDeltaT);
1892 }
1893 }
1894 if(dHistMap_DeltadEdXVsP[SYS_START].find(locPID) != dHistMap_DeltadEdXVsP[SYS_START].end())
1895 {
1896 double locdx = locSCHitMatchParams->dHitEnergy/locSCHitMatchParams->dEdx;
1897 double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
1898 locParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
1899 dHistMap_DeltadEdXVsP[SYS_START][locPID]->Fill(locP, (locSCHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
1900 }
1901 }
1902 if(locTOFHitMatchParams != NULL__null)
1903 {
1904 dHistMap_dEdXVsP[SYS_TOF][locCharge]->Fill(locP, locTOFHitMatchParams->dEdx*1.0E3);
1905 if(locEventRFBunch->dNumParticleVotes > 1)
1906 {
1907 double locBeta_Timing = locTOFHitMatchParams->dPathLength/(29.9792458*(locTOFHitMatchParams->dHitTime - locChargedTrackHypothesis->t0()));
1908 dHistMap_BetaVsP[SYS_TOF][locCharge]->Fill(locP, locBeta_Timing);
1909 if(dHistMap_DeltaBetaVsP[SYS_TOF].find(locPID) != dHistMap_DeltaBetaVsP[SYS_TOF].end())
1910 {
1911 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1912 dHistMap_DeltaBetaVsP[SYS_TOF][locPID]->Fill(locP, locDeltaBeta);
1913 double locDeltaT = locTOFHitMatchParams->dHitTime - locTOFHitMatchParams->dFlightTime - locStartTime;
1914 dHistMap_DeltaTVsP[SYS_TOF][locPID]->Fill(locP, locDeltaT);
1915 }
1916 }
1917 if(dHistMap_DeltadEdXVsP[SYS_TOF].find(locPID) != dHistMap_DeltadEdXVsP[SYS_TOF].end())
1918 {
1919 double locdx = locTOFHitMatchParams->dHitEnergy/locTOFHitMatchParams->dEdx;
1920 double locProbabledEdx = 0.0, locSigmadEdx = 0.0;
1921 locParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);
1922 dHistMap_DeltadEdXVsP[SYS_TOF][locPID]->Fill(locP, (locTOFHitMatchParams->dEdx - locProbabledEdx)*1.0E3);
1923 }
1924 }
1925 if(locBCALShowerMatchParams != NULL__null)
1926 {
1927 const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
1928 double locEOverP = locBCALShower->E/locP;
1929 dHistMap_BCALEOverPVsP[locCharge]->Fill(locP, locEOverP);
1930 dHistMap_BCALEOverPVsTheta[locCharge]->Fill(locTheta, locEOverP);
1931 if(locEventRFBunch->dNumParticleVotes > 1)
1932 {
1933 double locBeta_Timing = locBCALShowerMatchParams->dPathLength/(29.9792458*(locBCALShower->t - locChargedTrackHypothesis->t0()));
1934 dHistMap_BetaVsP[SYS_BCAL][locCharge]->Fill(locP, locBeta_Timing);
1935 if(dHistMap_DeltaBetaVsP[SYS_BCAL].find(locPID) != dHistMap_DeltaBetaVsP[SYS_BCAL].end())
1936 {
1937 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1938 dHistMap_DeltaBetaVsP[SYS_BCAL][locPID]->Fill(locP, locDeltaBeta);
1939 double locDeltaT = locBCALShower->t - locBCALShowerMatchParams->dFlightTime - locStartTime;
1940 dHistMap_DeltaTVsP[SYS_BCAL][locPID]->Fill(locP, locDeltaT);
1941 }
1942 }
1943 }
1944 if(locFCALShowerMatchParams != NULL__null)
1945 {
1946 const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
1947 double locEOverP = locFCALShower->getEnergy()/locP;
1948 dHistMap_FCALEOverPVsP[locCharge]->Fill(locP, locEOverP);
1949 dHistMap_FCALEOverPVsTheta[locCharge]->Fill(locTheta, locEOverP);
1950 if(locEventRFBunch->dNumParticleVotes > 1)
1951 {
1952 double locBeta_Timing = locFCALShowerMatchParams->dPathLength/(29.9792458*(locFCALShower->getTime() - locChargedTrackHypothesis->t0()));
1953 dHistMap_BetaVsP[SYS_FCAL][locCharge]->Fill(locP, locBeta_Timing);
1954 if(dHistMap_DeltaBetaVsP[SYS_FCAL].find(locPID) != dHistMap_DeltaBetaVsP[SYS_FCAL].end())
1955 {
1956 double locDeltaBeta = locChargedTrackHypothesis->lorentzMomentum().Beta() - locBeta_Timing;
1957 dHistMap_DeltaBetaVsP[SYS_FCAL][locPID]->Fill(locP, locDeltaBeta);
1958 double locDeltaT = locFCALShower->getTime() - locFCALShowerMatchParams->dFlightTime - locStartTime;
1959 dHistMap_DeltaTVsP[SYS_FCAL][locPID]->Fill(locP, locDeltaT);
1960 }
1961 }
1962 }
1963
1964 if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC > 0)
1965 {
1966 dHistMap_dEdXVsP[SYS_CDC][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_CDC*1.0E6);
1967 if(dHistMap_DeltadEdXVsP[SYS_CDC].find(locPID) != dHistMap_DeltadEdXVsP[SYS_CDC].end())
1968 {
1969 double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_CDC, true);
1970 dHistMap_DeltadEdXVsP[SYS_CDC][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_CDC - locProbabledEdx)*1.0E6);
1971 }
1972 }
1973 if(locTrackTimeBased->dNumHitsUsedFordEdx_FDC > 0)
1974 {
1975 dHistMap_dEdXVsP[SYS_FDC][locCharge]->Fill(locP, locTrackTimeBased->ddEdx_FDC*1.0E6);
1976 if(dHistMap_DeltadEdXVsP[SYS_FDC].find(locPID) != dHistMap_DeltadEdXVsP[SYS_FDC].end())
1977 {
1978 double locProbabledEdx = locParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_FDC, false);
1979 dHistMap_DeltadEdXVsP[SYS_FDC][locPID]->Fill(locP, (locTrackTimeBased->ddEdx_FDC - locProbabledEdx)*1.0E6);
1980 }
1981 }
1982 }
1983 }
1984 Unlock_Action(); //RELEASE ROOT LOCK!!
1985
1986 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
1987}
1988
1989void DHistogramAction_Neutrals::Initialize(JEventLoop* locEventLoop)
1990{
1991 //Create any histograms/trees/etc. within a ROOT lock.
1992 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
1993
1994 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
1995 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
1996
1997 string locHistName;
1998
1999 DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
2000 DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
2001 double locTargetZCenter = 0.0;
2002 locGeometry->GetTargetZ(locTargetZCenter);
2003
2004 //CREATE THE HISTOGRAMS
2005 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2006 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2007 {
2008 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
2009 //If another thread has already created the folder, it just changes to it.
2010 CreateAndChangeTo_ActionDirectory();
2011
2012 if(dTargetCenter.Z() < -9.8E9)
2013 dTargetCenter.SetXYZ(0.0, 0.0, locTargetZCenter);
2014
2015 //BCAL
2016 locHistName = "BCALTrackDOCA";
2017 dHist_BCALTrackDOCA = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower Distance to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
2018 locHistName = "BCALTrackDeltaPhi";
2019 dHist_BCALTrackDeltaPhi = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #Delta#phi#circ to Nearest Track", dNumDeltaPhiBins, dMinDeltaPhi, dMaxDeltaPhi);
2020 locHistName = "BCALTrackDeltaZ";
2021 dHist_BCALTrackDeltaZ = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Shower #DeltaZ to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
2022 locHistName = "BCALNeutralShowerEnergy";
2023 dHist_BCALNeutralShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Neutral Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxBCALP);
2024 locHistName = "BCALNeutralShowerDeltaT";
2025 dHist_BCALNeutralShowerDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, ";BCAL Neutral Shower #Deltat (Propagated - RF) (ns)", dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
2026 locHistName = "BCALNeutralShowerDeltaTVsE";
2027 dHist_BCALNeutralShowerDeltaTVsE = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Neutral Shower Energy (GeV);BCAL Neutral Shower #Deltat (ns)", dNum2DShowerEnergyBins, dMinShowerEnergy, dMaxBCALP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
2028 locHistName = "BCALNeutralShowerDeltaTVsZ";
2029 dHist_BCALNeutralShowerDeltaTVsZ = GetOrCreate_Histogram<TH2I>(locHistName, ";BCAL Neutral Shower Z (cm);BCAL Neutral Shower #Deltat (ns)", dNum2DBCALZBins, 0.0, 450.0, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
2030
2031 //FCAL
2032 locHistName = "FCALTrackDOCA";
2033 dHist_FCALTrackDOCA = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Shower Distance to Nearest Track (cm)", dNumTrackDOCABins, dMinTrackDOCA, dMaxTrackDOCA);
2034 locHistName = "FCALNeutralShowerEnergy";
2035 dHist_FCALNeutralShowerEnergy = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Neutral Shower Energy (GeV)", dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2036 locHistName = "FCALNeutralShowerDeltaT";
2037 dHist_FCALNeutralShowerDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, ";FCAL Neutral Shower #Deltat (Propagated - RF) (ns)", dNumDeltaTBins, dMinDeltaT, dMaxDeltaT);
2038 locHistName = "FCALNeutralShowerDeltaTVsE";
2039 dHist_FCALNeutralShowerDeltaTVsE = GetOrCreate_Histogram<TH2I>(locHistName, ";FCAL Neutral Shower Energy (GeV);FCAL Neutral Shower #Deltat (ns)", dNum2DShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);
2040
2041 //Return to the base directory
2042 ChangeTo_BaseDirectory();
2043 }
2044 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2045}
2046
2047bool DHistogramAction_Neutrals::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2048{
2049 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
2050
2051 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
2052 if(Get_NumPreviousParticleCombos() != 0)
2053 return true;
2054
2055 vector<const DNeutralShower*> locNeutralShowers;
2056 locEventLoop->Get(locNeutralShowers);
2057
2058 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
2059 locEventLoop->Get(locTrackTimeBasedVector);
2060
2061 const DDetectorMatches* locDetectorMatches = NULL__null;
2062 locEventLoop->GetSingle(locDetectorMatches);
2063
2064 vector<const DEventRFBunch*> locEventRFBunches;
2065 locEventLoop->Get(locEventRFBunches);
2066 double locStartTime = locEventRFBunches.empty() ? 0.0 : locEventRFBunches[0]->dTime;
2067
2068 //FILL HISTOGRAMS
2069 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2070 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2071 Lock_Action(); //ACQUIRE ROOT LOCK!!
2072 {
2073 for(size_t loc_i = 0; loc_i < locNeutralShowers.size(); ++loc_i)
2074 {
2075 //assume is photon
2076 double locPathLength = (locNeutralShowers[loc_i]->dSpacetimeVertex.Vect() - dTargetCenter).Mag();
2077 double locDeltaT = locNeutralShowers[loc_i]->dSpacetimeVertex.T() - locPathLength/29.9792458 - locStartTime;
2078
2079 if(locNeutralShowers[loc_i]->dDetectorSystem == SYS_FCAL)
2080 {
2081 const DFCALShower* locFCALShower = NULL__null;
2082 locNeutralShowers[loc_i]->GetSingle(locFCALShower);
2083
2084 double locDistance = 9.9E9;
2085 if(locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistance))
2086 dHist_FCALTrackDOCA->Fill(locDistance);
2087
2088 dHist_FCALNeutralShowerEnergy->Fill(locNeutralShowers[loc_i]->dEnergy);
2089 dHist_FCALNeutralShowerDeltaT->Fill(locDeltaT);
2090 dHist_FCALNeutralShowerDeltaTVsE->Fill(locNeutralShowers[loc_i]->dEnergy, locDeltaT);
2091 }
2092 else
2093 {
2094 const DBCALShower* locBCALShower = NULL__null;
2095 locNeutralShowers[loc_i]->GetSingle(locBCALShower);
2096
2097 double locDistance = 9.9E9, locDeltaPhi = 9.9E9, locDeltaZ = 9.9E9;
2098 if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDistance))
2099 dHist_BCALTrackDOCA->Fill(locDistance);
2100 if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDeltaPhi, locDeltaZ))
2101 {
2102 dHist_BCALTrackDeltaPhi->Fill(180.0*locDeltaPhi/TMath::Pi());
2103 dHist_BCALTrackDeltaZ->Fill(locDeltaZ);
2104 }
2105
2106 dHist_BCALNeutralShowerEnergy->Fill(locNeutralShowers[loc_i]->dEnergy);
2107 dHist_BCALNeutralShowerDeltaT->Fill(locDeltaT);
2108 dHist_BCALNeutralShowerDeltaTVsE->Fill(locNeutralShowers[loc_i]->dEnergy, locDeltaT);
2109 dHist_BCALNeutralShowerDeltaTVsZ->Fill(locNeutralShowers[loc_i]->dSpacetimeVertex.Z(), locDeltaT);
2110 }
2111 }
2112 }
2113 Unlock_Action(); //RELEASE ROOT LOCK!!
2114
2115 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2116}
2117
2118
2119void DHistogramAction_DetectorMatchParams::Initialize(JEventLoop* locEventLoop)
2120{
2121 //Create any histograms/trees/etc. within a ROOT lock.
2122 //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time.
2123
2124 //When creating a reaction-independent action, only modify member variables within a ROOT lock.
2125 //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously.
2126
2127 string locHistName, locHistTitle;
2128
2129 vector<const DMCThrown*> locMCThrowns;
2130 locEventLoop->Get(locMCThrowns);
2131
2132 DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
2133 DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
2134 double locTargetZCenter = 0.0;
2135 locGeometry->GetTargetZ(locTargetZCenter);
2136
2137 string locTrackSelectionTag = "NotATag";
2138 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2139 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2140
2141 //CREATE THE HISTOGRAMS
2142 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2143 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2144 {
2145 //So: Default tag is "", User can set it to something else
2146 //In here, if tag is "", get from gparms, if not, leave it alone
2147 //If gparms value does not exist, set it to (and use) "PreSelect"
2148 if(dTrackSelectionTag == "NotATag")
2149 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2150
2151 //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
2152 //If another thread has already created the folder, it just changes to it.
2153 CreateAndChangeTo_ActionDirectory();
2154
2155 if(dTargetCenterZ < -9.8E9)
2156 dTargetCenterZ = locTargetZCenter; //only set if not already set
2157
2158 //Track Matched to Hit
2159 for(int locTruePIDFlag = 0; locTruePIDFlag < 2; ++locTruePIDFlag)
2160 {
2161 if(locMCThrowns.empty() && (locTruePIDFlag == 1))
2162 continue; //not a simulated event: don't histogram thrown info!
2163
2164 string locDirName = (locTruePIDFlag == 1) ? "TruePID" : "ReconstructedPID";
2165 CreateAndChangeTo_Directory(locDirName.c_str(), locDirName.c_str());
2166
2167 //By PID
2168 for(int loc_i = -2; loc_i < int(dTrackingPIDs.size()); ++loc_i) //-2 = q-, -1 = q+
2169 {
2170 string locParticleName, locParticleROOTName;
2171 int locPID = loc_i;
2172 if(loc_i == -2)
2173 {
2174 locParticleName = "q-";
2175 locParticleROOTName = "#it{q}^{-}";
2176 }
2177 else if(loc_i == -1)
2178 {
2179 locParticleName = "q+";
2180 locParticleROOTName = "#it{q}^{+}";
2181 }
2182 else
2183 {
2184 locParticleName = ParticleType(dTrackingPIDs[loc_i]);
2185 locParticleROOTName = ParticleName_ROOT(dTrackingPIDs[loc_i]);
2186 locPID = int(dTrackingPIDs[loc_i]);
2187 }
2188 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2189 pair<int, bool> locPIDPair(locPID, bool(locTruePIDFlag));
2190
2191 //BCAL
2192 locHistName = "BCALShowerEnergy";
2193 locHistTitle = locParticleROOTName + ";BCAL Shower Energy (GeV)";
2194 dHistMap_BCALShowerEnergy[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerEnergyBins, dMinShowerEnergy, dMaxBCALP);
2195
2196 locHistName = "BCALShowerTrackDepth";
2197 locHistTitle = locParticleROOTName + ";BCAL Shower Track Depth (cm)";
2198 dHistMap_BCALShowerTrackDepth[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2199
2200 locHistName = "BCALShowerTrackDepthVsP";
2201 locHistTitle = locParticleROOTName + ";p (GeV/c);BCAL Shower Track Depth (cm)";
2202 dHistMap_BCALShowerTrackDepthVsP[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2203
2204 //FCAL
2205 locHistName = "FCALShowerEnergy";
2206 locHistTitle = locParticleROOTName + ";FCAL Shower Energy (GeV)";
2207 dHistMap_FCALShowerEnergy[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerEnergyBins, dMinShowerEnergy, dMaxShowerEnergy);
2208
2209 locHistName = "FCALShowerTrackDepth";
2210 locHistTitle = locParticleROOTName + ";FCAL Shower Track Depth (cm)";
2211 dHistMap_FCALShowerTrackDepth[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2212
2213 locHistName = "FCALShowerTrackDepthVsP";
2214 locHistTitle = locParticleROOTName + ";p (GeV/c);FCAL Shower Track Depth (cm)";
2215 dHistMap_FCALShowerTrackDepthVsP[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumShowerDepthBins, dMinShowerDepth, dMaxShowerDepth);
2216
2217 //SC
2218 locHistName = "SCEnergyVsTheta";
2219 locHistTitle = locParticleROOTName + ";#theta#circ;SC Point Energy (MeV)";
2220 dHistMap_SCEnergyVsTheta[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DHitEnergyBins, dMinHitEnergy, dMaxHitEnergy);
2221
2222 locHistName = "SCPhiVsTheta";
2223 locHistTitle = locParticleROOTName + ";#theta#circ;#phi#circ";
2224 dHistMap_SCPhiVsTheta[locPIDPair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
2225
2226 gDirectory(TDirectory::CurrentDirectory())->cd(".."); //end of PID
2227 }
2228 gDirectory(TDirectory::CurrentDirectory())->cd(".."); //end of true/recon PID
2229 }
2230
2231 //Return to the base directory
2232 ChangeTo_BaseDirectory();
2233 }
2234 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2235}
2236
2237bool DHistogramAction_DetectorMatchParams::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2238{
2239 //Expect locParticleCombo to be NULL since this is a reaction-independent action.
2240
2241 //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms)
2242 if(Get_NumPreviousParticleCombos() != 0)
2243 return true;
2244
2245 vector<const DMCThrown*> locMCThrowns;
2246 locEventLoop->Get(locMCThrowns);
2247
2248 Fill_Hists(locEventLoop, false);
2249 if(!locMCThrowns.empty())
2250 Fill_Hists(locEventLoop, true);
2251
2252 return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
2253}
2254
2255void DHistogramAction_DetectorMatchParams::Fill_Hists(JEventLoop* locEventLoop, bool locUseTruePIDFlag)
2256{
2257 vector<const DChargedTrack*> locChargedTracks;
2258 locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
2259
2260 vector<const DMCThrownMatching*> locMCThrownMatchingVector;
2261 locEventLoop->Get(locMCThrownMatchingVector);
2262
2263 const DEventRFBunch* locEventRFBunch = NULL__null;
2264 locEventLoop->GetSingle(locEventRFBunch);
2265
2266 //FILL HISTOGRAMS
2267 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2268 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2269 Lock_Action(); //ACQUIRE ROOT LOCK!!
2270 {
2271 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2272 {
2273 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2274
2275 if(locUseTruePIDFlag && (!locMCThrownMatchingVector.empty()))
2276 {
2277 double locMatchFOM = 0.0;
2278 const DMCThrown* locMCThrown = locMCThrownMatchingVector[0]->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);
2279 if(locMCThrown == NULL__null)
2280 continue;
2281 //OK, have the thrown. Now, grab the best charged track hypothesis to get the best matching
2282 locChargedTrackHypothesis = locMCThrownMatchingVector[0]->Get_MatchingChargedHypothesis(locMCThrown, locMatchFOM);
2283 }
2284
2285 pair<int, bool> locPIDPair(int(locChargedTrackHypothesis->PID()), locUseTruePIDFlag);
2286 bool locDisregardPIDFlag = (dHistMap_BCALShowerEnergy.find(locPIDPair) == dHistMap_BCALShowerEnergy.end());
2287 int locQIndex = (locChargedTrackHypothesis->charge() > 0.0) ? -1 : -2;
2288 pair<int, bool> locChargePair(locQIndex, locUseTruePIDFlag);
2289
2290 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2291 const DFCALShowerMatchParams* locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
2292 const DSCHitMatchParams* locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
2293 const DBCALShowerMatchParams* locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
2294
2295 //BCAL
2296 if(locBCALShowerMatchParams != NULL__null)
2297 {
2298 const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
2299 dHistMap_BCALShowerEnergy[locChargePair]->Fill(locBCALShower->E);
2300 dHistMap_BCALShowerTrackDepth[locChargePair]->Fill(locBCALShowerMatchParams->dx);
2301 dHistMap_BCALShowerTrackDepthVsP[locChargePair]->Fill(locMomentum.Mag(), locBCALShowerMatchParams->dx);
2302
2303 if(!locDisregardPIDFlag)
2304 {
2305 dHistMap_BCALShowerEnergy[locPIDPair]->Fill(locBCALShower->E);
2306 dHistMap_BCALShowerTrackDepth[locPIDPair]->Fill(locBCALShowerMatchParams->dx);
2307 dHistMap_BCALShowerTrackDepthVsP[locPIDPair]->Fill(locMomentum.Mag(), locBCALShowerMatchParams->dx);
2308 }
2309 }
2310
2311 //FCAL
2312 if(locFCALShowerMatchParams != NULL__null)
2313 {
2314 const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
2315 dHistMap_FCALShowerEnergy[locChargePair]->Fill(locFCALShower->getEnergy());
2316 dHistMap_FCALShowerTrackDepth[locChargePair]->Fill(locFCALShowerMatchParams->dx);
2317 dHistMap_FCALShowerTrackDepthVsP[locChargePair]->Fill(locMomentum.Mag(), locFCALShowerMatchParams->dx);
2318
2319 if(!locDisregardPIDFlag)
2320 {
2321 dHistMap_FCALShowerEnergy[locPIDPair]->Fill(locFCALShower->getEnergy());
2322 dHistMap_FCALShowerTrackDepth[locPIDPair]->Fill(locFCALShowerMatchParams->dx);
2323 dHistMap_FCALShowerTrackDepthVsP[locPIDPair]->Fill(locMomentum.Mag(), locFCALShowerMatchParams->dx);
2324 }
2325 }
2326
2327 //SC
2328 if(locSCHitMatchParams != NULL__null)
2329 {
2330 dHistMap_SCEnergyVsTheta[locChargePair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locSCHitMatchParams->dHitEnergy*1.0E3);
2331 dHistMap_SCPhiVsTheta[locChargePair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locMomentum.Phi()*180.0/TMath::Pi());
2332
2333 if(!locDisregardPIDFlag)
2334 {
2335 dHistMap_SCEnergyVsTheta[locPIDPair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locSCHitMatchParams->dHitEnergy*1.0E3);
2336 dHistMap_SCPhiVsTheta[locPIDPair]->Fill(locMomentum.Theta()*180.0/TMath::Pi(), locMomentum.Phi()*180.0/TMath::Pi());
2337 }
2338 }
2339 }
2340 }
2341 Unlock_Action(); //RELEASE ROOT LOCK!!
2342}
2343
2344void DHistogramAction_EventVertex::Initialize(JEventLoop* locEventLoop)
2345{
2346 string locHistName, locHistTitle;
2347
2348 string locTrackSelectionTag = "NotATag";
2349 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2350 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2351
2352 //CREATE THE HISTOGRAMS
2353 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2354 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2355 {
2356 //So: Default tag is "", User can set it to something else
2357 //In here, if tag is "", get from gparms, if not, leave it alone
2358 //If gparms value does not exist, set it to (and use) "PreSelect"
2359 if(dTrackSelectionTag == "NotATag")
2360 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2361
2362 CreateAndChangeTo_ActionDirectory();
2363
2364 // Event RF Bunch Time
2365 locHistName = "RFTrackDeltaT";
2366 locHistTitle = ";#Deltat_{RF - Track} (ns)";
2367 dRFTrackDeltaT = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumRFTBins, -3.0, 3.0);
2368
2369 // ALL EVENTS
2370 CreateAndChangeTo_Directory("AllEvents", "AllEvents");
2371
2372 // Event Vertex-Z
2373 locHistName = "EventVertexZ";
2374 locHistTitle = ";Event Vertex-Z (cm)";
2375 dEventVertexZ_AllEvents = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2376
2377 // Event Vertex-Y Vs Vertex-X
2378 locHistName = "EventVertexYVsX";
2379 locHistTitle = ";Event Vertex-X (cm);Event Vertex-Y (cm)";
2380 dEventVertexYVsX_AllEvents = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2381
2382 // Event Vertex-T
2383 locHistName = "EventVertexT";
2384 locHistTitle = ";Event Vertex Time (ns)";
2385 dEventVertexT_AllEvents = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2386
2387 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2388
2389
2390 // 2+ Good Tracks
2391 CreateAndChangeTo_Directory("2+GoodTracks", "2+GoodTracks");
2392
2393 // Event Vertex-Z
2394 locHistName = "EventVertexZ";
2395 locHistTitle = ";Event Vertex-Z (cm)";
2396 dEventVertexZ_2OrMoreGoodTracks = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2397
2398 // Event Vertex-Y Vs Vertex-X
2399 locHistName = "EventVertexYVsX";
2400 locHistTitle = ";Event Vertex-X (cm);Event Vertex-Y (cm)";
2401 dEventVertexYVsX_2OrMoreGoodTracks = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2402
2403 // Event Vertex-T
2404 locHistName = "EventVertexT";
2405 locHistTitle = ";Event Vertex Time (ns)";
2406 dEventVertexT_2OrMoreGoodTracks = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2407
2408 // Confidence Level
2409 locHistName = "ConfidenceLevel";
2410 dHist_KinFitConfidenceLevel = GetOrCreate_Histogram<TH1I>(locHistName, "Event Vertex Kinematic Fit;Confidence Level;# Events", dNumConfidenceLevelBins, 0.0, 1.0);
2411
2412 //final particle pulls
2413 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2414 {
2415 Particle_t locPID = dFinalStatePIDs[loc_i];
2416 string locParticleDirName = string("Pulls_") + string(ParticleType(locPID));
2417 string locParticleROOTName = ParticleName_ROOT(locPID);
2418 CreateAndChangeTo_Directory(locParticleDirName, locParticleDirName);
2419
2420 //Px Pull
2421 locHistName = "Pull_Px";
2422 locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{x} Pull;# Events");
2423 dHistMap_KinFitPulls[locPID][d_PxPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2424
2425 //Py Pull
2426 locHistName = "Pull_Py";
2427 locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{y} Pull;# Events");
2428 dHistMap_KinFitPulls[locPID][d_PyPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2429
2430 //Pz Pull
2431 locHistName = "Pull_Pz";
2432 locHistTitle = locParticleROOTName + string(", Vertex Fit;p_{z} Pull;# Events");
2433 dHistMap_KinFitPulls[locPID][d_PzPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2434
2435 //Xx Pull
2436 locHistName = "Pull_Xx";
2437 locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{x} Pull;# Events");
2438 dHistMap_KinFitPulls[locPID][d_XxPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2439
2440 //Xy Pull
2441 locHistName = "Pull_Xy";
2442 locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{y} Pull;# Events");
2443 dHistMap_KinFitPulls[locPID][d_XyPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2444
2445 //Xz Pull
2446 locHistName = "Pull_Xz";
2447 locHistTitle = locParticleROOTName + string(", Vertex Fit;x_{z} Pull;# Events");
2448 dHistMap_KinFitPulls[locPID][d_XzPull] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);
2449
2450 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2451 }
2452
2453 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2454
2455 //Return to the base directory
2456 ChangeTo_BaseDirectory();
2457 }
2458 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2459}
2460
2461bool DHistogramAction_EventVertex::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2462{
2463 if(Get_NumPreviousParticleCombos() != 0)
2464 return true; //else double-counting!
2465
2466 const DVertex* locVertex = NULL__null;
2467 locEventLoop->GetSingle(locVertex);
2468
2469 vector<const DChargedTrack*> locChargedTracks;
2470 locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());
2471
2472 const DDetectorMatches* locDetectorMatches = NULL__null;
2473 locEventLoop->GetSingle(locDetectorMatches);
2474
2475 const DParticleID* locParticleID = NULL__null;
2476 locEventLoop->GetSingle(locParticleID);
2477
2478 const DEventRFBunch* locEventRFBunch = NULL__null;
2479 locEventLoop->GetSingle(locEventRFBunch);
2480
2481 //Make sure that brun() is called (to get rf period) before using.
2482 //Cannot call JEventLoop->Get() because object may be in datastream (REST), bypassing factory brun() call.
2483 //Must do here rather than in Initialize() function because this object is shared by all threads (which each have their own factory)
2484 DRFTime_factory* locRFTimeFactory = static_cast<DRFTime_factory*>(locEventLoop->GetFactory("DRFTime"));
2485 if(!locRFTimeFactory->brun_was_called())
2486 {
2487 locRFTimeFactory->brun(locEventLoop, locEventLoop->GetJEvent().GetRunNumber());
2488 locRFTimeFactory->Set_brun_called();
2489 }
2490
2491 //Get time-based tracks: use best PID FOM
2492 //Note that these may not be the PIDs that were used in the fit!!!
2493 //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
2494 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
2495 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2496 {
2497 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2498 const DTrackTimeBased* locTrackTimeBased = NULL__null;
2499 locChargedTrackHypothesis->GetSingle(locTrackTimeBased);
2500 if(locTrackTimeBased != NULL__null)
2501 locTrackTimeBasedVector.push_back(locTrackTimeBased);
2502 }
2503
2504 //Event Vertex
2505 //FILL HISTOGRAMS
2506 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2507 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2508 Lock_Action();
2509 {
2510 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
2511 {
2512 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
2513 double locPropagatedRFTime = locParticleID->Calc_PropagatedRFTime(locChargedTrackHypothesis, locEventRFBunch);
2514 double locShiftedRFTime = locRFTimeFactory->Step_TimeToNearInputTime(locPropagatedRFTime, locChargedTrackHypothesis->time());
2515 double locDeltaT = locShiftedRFTime - locChargedTrackHypothesis->time();
2516 dRFTrackDeltaT->Fill(locDeltaT);
2517 }
2518 dEventVertexZ_AllEvents->Fill(locVertex->dSpacetimeVertex.Z());
2519 dEventVertexYVsX_AllEvents->Fill(locVertex->dSpacetimeVertex.X(), locVertex->dSpacetimeVertex.Y());
2520 dEventVertexT_AllEvents->Fill(locVertex->dSpacetimeVertex.T());
2521
2522 if(locChargedTracks.size() >= 2)
2523 {
2524 dEventVertexZ_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.Z());
2525 dEventVertexYVsX_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.X(), locVertex->dSpacetimeVertex.Y());
2526 dEventVertexT_2OrMoreGoodTracks->Fill(locVertex->dSpacetimeVertex.T());
2527 }
2528 }
2529 Unlock_Action();
2530
2531 if(locVertex->dKinFitNDF == 0)
2532 return true; //kin fit not performed or didn't converge: no results to histogram
2533
2534 double locConfidenceLevel = TMath::Prob(locVertex->dKinFitChiSq, locVertex->dKinFitNDF);
2535
2536 //FILL HISTOGRAMS
2537 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2538 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2539 Lock_Action();
2540 {
2541 dHist_KinFitConfidenceLevel->Fill(locConfidenceLevel);
2542
2543 //pulls
2544 if(locConfidenceLevel > dPullHistConfidenceLevelCut)
2545 {
2546 for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
2547 {
2548 const DKinematicData* locKinematicData = static_cast<const DKinematicData*>(locTrackTimeBasedVector[loc_i]);
2549 Particle_t locPID = locKinematicData->PID();
2550 if(dHistMap_KinFitPulls.find(locPID) == dHistMap_KinFitPulls.end())
2551 continue; //PID not histogrammed
2552
2553 map<const JObject*, map<DKinFitPullType, double> >::const_iterator locParticleIterator = locVertex->dKinFitPulls.find(locKinematicData);
2554 if(locParticleIterator == locVertex->dKinFitPulls.end())
2555 continue;
2556
2557 const map<DKinFitPullType, double>& locPullMap = locParticleIterator->second;
2558 map<DKinFitPullType, double>::const_iterator locPullIterator = locPullMap.begin();
2559 for(; locPullIterator != locPullMap.end(); ++locPullIterator)
2560 dHistMap_KinFitPulls[locPID][locPullIterator->first]->Fill(locPullIterator->second);
2561 }
2562 }
2563 }
2564 Unlock_Action();
2565
2566 return true;
2567}
2568
2569void DHistogramAction_DetectedParticleKinematics::Initialize(JEventLoop* locEventLoop)
2570{
2571 string locHistName, locHistTitle, locParticleName, locParticleROOTName;
2572 Particle_t locPID;
2573
2574 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
2575 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2576 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2577 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
2578 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
2579
2580 //CREATE THE HISTOGRAMS
2581 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2582 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2583 {
2584 //So: Default tag is "", User can set it to something else
2585 //In here, if tag is "", get from gparms, if not, leave it alone
2586 //If gparms value does not exist, set it to (and use) "PreSelect"
2587 if(dTrackSelectionTag == "NotATag")
2588 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2589 if(dShowerSelectionTag == "NotATag")
2590 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
2591
2592 CreateAndChangeTo_ActionDirectory();
2593
2594 // Beam Particle
2595 locPID = Gamma;
2596 locParticleName = string("Beam_") + ParticleType(locPID);
2597 locParticleROOTName = ParticleName_ROOT(locPID);
2598 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2599 locHistName = "Momentum";
2600 locHistTitle = string("Beam ") + locParticleROOTName + string(";p (GeV/c)");
2601 dBeamParticle_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBeamEBins, dMinP, dMaxBeamE);
2602 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2603
2604 //PID
2605 CreateAndChangeTo_Directory("PID", "PID");
2606 {
2607 //beta vs p
2608 locHistName = "BetaVsP_Q+";
2609 locHistTitle = "q^{+};p (GeV/c);#beta";
2610 dHistMap_QBetaVsP[1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2611
2612 locHistName = "BetaVsP_Q-";
2613 locHistTitle = "q^{-};p (GeV/c);#beta";
2614 dHistMap_QBetaVsP[-1] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2615 }
2616 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2617
2618 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2619 {
2620 locPID = dFinalStatePIDs[loc_i];
2621 locParticleName = ParticleType(locPID);
2622 locParticleROOTName = ParticleName_ROOT(locPID);
2623 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2624
2625 // Momentum
2626 locHistName = "Momentum";
2627 locHistTitle = locParticleROOTName + string(";p (GeV/c)");
2628 dHistMap_P[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);
2629
2630 // Theta
2631 locHistName = "Theta";
2632 locHistTitle = locParticleROOTName + string(";#theta#circ");
2633 dHistMap_Theta[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);
2634
2635 // Phi
2636 locHistName = "Phi";
2637 locHistTitle = locParticleROOTName + string(";#phi#circ");
2638 dHistMap_Phi[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);
2639
2640 // P Vs Theta
2641 locHistName = "PVsTheta";
2642 locHistTitle = locParticleROOTName + string(";#theta#circ;p (GeV/c)");
2643 dHistMap_PVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);
2644
2645 // Phi Vs Theta
2646 locHistName = "PhiVsTheta";
2647 locHistTitle = locParticleROOTName + string(";#theta#circ;#phi#circ");
2648 dHistMap_PhiVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);
2649
2650 //beta vs p
2651 locHistName = "BetaVsP";
2652 locHistTitle = locParticleROOTName + string(";p (GeV/c);#beta");
2653 dHistMap_BetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);
2654
2655 //delta-beta vs p
2656 locHistName = "DeltaBetaVsP";
2657 locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#beta");
2658 dHistMap_DeltaBetaVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);
2659
2660 // Vertex-Z
2661 locHistName = "VertexZ";
2662 locHistTitle = locParticleROOTName + string(";Vertex-Z (cm)");
2663 dHistMap_VertexZ[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);
2664
2665 // Vertex-Y Vs Vertex-X
2666 locHistName = "VertexYVsX";
2667 locHistTitle = locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");
2668 dHistMap_VertexYVsX[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);
2669
2670 // Vertex-T
2671 locHistName = "VertexT";
2672 locHistTitle = locParticleROOTName + string(";Vertex-T (ns)");
2673 dHistMap_VertexT[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumTBins, dMinT, dMaxT);
2674
2675 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2676 }
2677
2678 //Return to the base directory
2679 ChangeTo_BaseDirectory();
2680 }
2681 japp->RootUnLock(); //RELEASE ROOT LOCK!!
2682}
2683
2684bool DHistogramAction_DetectedParticleKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
2685{
2686 if(Get_NumPreviousParticleCombos() != 0)
2687 return true; //else double-counting!
2688
2689 vector<const DBeamPhoton*> locBeamPhotons;
2690 locEventLoop->Get(locBeamPhotons);
2691
2692 //FILL HISTOGRAMS
2693 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2694 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2695 Lock_Action();
2696 {
2697 for(size_t loc_i = 0; loc_i < locBeamPhotons.size(); ++loc_i)
2698 dBeamParticle_P->Fill(locBeamPhotons[loc_i]->energy());
2699 }
2700 Unlock_Action();
2701
2702 vector<const DChargedTrack*> locPreSelectChargedTracks;
2703 locEventLoop->Get(locPreSelectChargedTracks, dTrackSelectionTag.c_str());
2704
2705 for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
2706 {
2707 const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestTrackingFOM();
2708 int locCharge = ParticleCharge(locChargedTrackHypothesis->PID());
2709
2710 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2711 double locP = locMomentum.Mag();
2712 double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
2713
2714 if(dHistMap_QBetaVsP.find(locCharge) == dHistMap_QBetaVsP.end())
2715 continue;
2716
2717 //FILL HISTOGRAMS
2718 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2719 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2720 Lock_Action();
2721 {
2722 //Extremely inefficient, I know ...
2723 dHistMap_QBetaVsP[locCharge]->Fill(locP, locBeta_Timing);
2724 }
2725 Unlock_Action();
2726 }
2727
2728 for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
2729 {
2730 const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestFOM();
2731
2732 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
2733 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2734 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2735 double locP = locMomentum.Mag();
2736 double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();
2737 double locDeltaBeta = locChargedTrackHypothesis->deltaBeta();
2738
2739 Particle_t locPID = (locChargedTrackHypothesis->dFOM < dMinPIDFOM) ? Unknown : locChargedTrackHypothesis->PID();
2740 if(dHistMap_P.find(locPID) == dHistMap_P.end())
2741 continue; //not interested in histogramming
2742
2743 //FILL HISTOGRAMS
2744 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2745 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2746 Lock_Action();
2747 {
2748 dHistMap_P[locPID]->Fill(locP);
2749 dHistMap_Phi[locPID]->Fill(locPhi);
2750 dHistMap_Theta[locPID]->Fill(locTheta);
2751 dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
2752 dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
2753 dHistMap_BetaVsP[locPID]->Fill(locP, locBeta_Timing);
2754 dHistMap_DeltaBetaVsP[locPID]->Fill(locP, locDeltaBeta);
2755 dHistMap_VertexZ[locPID]->Fill(locChargedTrackHypothesis->position().Z());
2756 dHistMap_VertexYVsX[locPID]->Fill(locChargedTrackHypothesis->position().X(), locChargedTrackHypothesis->position().Y());
2757 dHistMap_VertexT[locPID]->Fill(locChargedTrackHypothesis->time());
2758 }
2759 Unlock_Action();
2760 }
2761
2762 vector<const DNeutralParticle*> locNeutralParticles;
2763 locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
2764
2765 for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
2766 {
2767 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(Gamma);
2768 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
2769 continue;
2770
2771 Particle_t locPID = locNeutralParticleHypothesis->PID();
2772 if(dHistMap_P.find(locPID) == dHistMap_P.end())
2773 continue; //e.g. a decaying particle, or not interested in histogramming
2774
2775 DVector3 locMomentum = locNeutralParticleHypothesis->momentum();
2776 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
2777 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
2778 double locP = locMomentum.Mag();
2779
2780 double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();
2781 double locDeltaBeta = locNeutralParticleHypothesis->deltaBeta();
2782
2783 //FILL HISTOGRAMS
2784 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
2785 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
2786 Lock_Action();
2787 {
2788 dHistMap_P[locPID]->Fill(locP);
2789 dHistMap_Phi[locPID]->Fill(locPhi);
2790 dHistMap_Theta[locPID]->Fill(locTheta);
2791 dHistMap_PVsTheta[locPID]->Fill(locTheta, locP);
2792 dHistMap_PhiVsTheta[locPID]->Fill(locTheta, locPhi);
2793 dHistMap_BetaVsP[locPID]->Fill(locP, locBeta_Timing);
2794 dHistMap_DeltaBetaVsP[locPID]->Fill(locP, locDeltaBeta);
2795 dHistMap_VertexZ[locPID]->Fill(locNeutralParticleHypothesis->position().Z());
2796 dHistMap_VertexYVsX[locPID]->Fill(locNeutralParticleHypothesis->position().X(), locNeutralParticleHypothesis->position().Y());
2797 dHistMap_VertexT[locPID]->Fill(locNeutralParticleHypothesis->time());
2798 }
2799 Unlock_Action();
2800 }
2801 return true;
2802}
2803
2804void DHistogramAction_TrackShowerErrors::Initialize(JEventLoop* locEventLoop)
2805{
2806 string locHistName, locHistTitle, locParticleName, locParticleROOTName;
2807 Particle_t locPID;
2808
2809 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
2810 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
2811 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
2812 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
2813 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
2814
2815 //CREATE THE HISTOGRAMS
2816 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
2817 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
2818 {
2819 //So: Default tag is "", User can set it to something else
2820 //In here, if tag is "", get from gparms, if not, leave it alone
2821 //If gparms value does not exist, set it to (and use) "PreSelect"
2822 if(dTrackSelectionTag == "NotATag")
2823 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
2824 if(dShowerSelectionTag == "NotATag")
2825 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
2826
2827 CreateAndChangeTo_ActionDirectory();
2828
2829 //TRACKS
2830 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
2831 {
2832 locPID = dFinalStatePIDs[loc_i];
2833 locParticleName = ParticleType(locPID);
2834 locParticleROOTName = ParticleName_ROOT(locPID);
2835 CreateAndChangeTo_Directory(locParticleName, locParticleName);
2836
2837 // Px
2838 locHistName = "PxErrorVsP";
2839 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{p_{x}}");
2840 dHistMap_TrackPxErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
2841
2842 locHistName = "PxErrorVsTheta";
2843 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{p_{x}}");
2844 dHistMap_TrackPxErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
2845
2846 locHistName = "PxErrorVsPhi";
2847 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{p_{x}}");
2848 dHistMap_TrackPxErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
2849
2850 // Py
2851 locHistName = "PyErrorVsP";
2852 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{p_{y}}");
2853 dHistMap_TrackPyErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
2854
2855 locHistName = "PyErrorVsTheta";
2856 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{p_{y}}");
2857 dHistMap_TrackPyErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
2858
2859 locHistName = "PyErrorVsPhi";
2860 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{p_{y}}");
2861 dHistMap_TrackPyErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DPxyErrorBins, 0.0, dMaxPxyError);
2862
2863 // Pz
2864 locHistName = "PzErrorVsP";
2865 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{p_{z}}");
2866 dHistMap_TrackPzErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPzErrorBins, 0.0, dMaxPzError);
2867
2868 locHistName = "PzErrorVsTheta";
2869 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{p_{z}}");
2870 dHistMap_TrackPzErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPzErrorBins, 0.0, dMaxPzError);
2871
2872 locHistName = "PzErrorVsPhi";
2873 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{p_{z}}");
2874 dHistMap_TrackPzErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DPzErrorBins, 0.0, dMaxPzError);
2875
2876 // X
2877 locHistName = "XErrorVsP";
2878 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{x}");
2879 dHistMap_TrackXErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
2880
2881 locHistName = "XErrorVsTheta";
2882 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{x}");
2883 dHistMap_TrackXErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
2884
2885 locHistName = "XErrorVsPhi";
2886 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{x}");
2887 dHistMap_TrackXErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
2888
2889 // Y
2890 locHistName = "YErrorVsP";
2891 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{y}");
2892 dHistMap_TrackYErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
2893
2894 locHistName = "YErrorVsTheta";
2895 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{y}");
2896 dHistMap_TrackYErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
2897
2898 locHistName = "YErrorVsPhi";
2899 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{y}");
2900 dHistMap_TrackYErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
2901
2902 // Z
2903 locHistName = "ZErrorVsP";
2904 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{z}");
2905 dHistMap_TrackZErrorVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DZErrorBins, 0.0, dMaxZError);
2906
2907 locHistName = "ZErrorVsTheta";
2908 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{z}");
2909 dHistMap_TrackZErrorVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DZErrorBins, 0.0, dMaxZError);
2910
2911 locHistName = "ZErrorVsPhi";
2912 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{z}");
2913 dHistMap_TrackZErrorVsPhi[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DZErrorBins, 0.0, dMaxZError);
2914
2915 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2916 }
2917
2918 //SHOWERS
2919 for(bool locIsBCALFlag : {false, true})
2920 {
2921 string locDirName = locIsBCALFlag ? "Photon_BCAL" : "Photon_FCAL";
2922 locParticleROOTName = ParticleName_ROOT(Gamma);
2923 CreateAndChangeTo_Directory(locDirName, locDirName);
2924
2925 double locMaxP = locIsBCALFlag ? dMaxPBCAL : dMaxP;
2926 double locMinTheta = locIsBCALFlag ? dMinThetaBCAL : dMinTheta;
2927 double locMaxTheta = locIsBCALFlag ? dMaxTheta : dMaxThetaFCAL;
2928
2929 // E
2930 locHistName = "EErrorVsP";
2931 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{E}");
2932 dHistMap_ShowerEErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DEErrorBins, 0.0, dMaxEError);
2933
2934 locHistName = "EErrorVsTheta";
2935 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{E}");
2936 dHistMap_ShowerEErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DEErrorBins, 0.0, dMaxEError);
2937
2938 locHistName = "EErrorVsPhi";
2939 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{E}");
2940 dHistMap_ShowerEErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DEErrorBins, 0.0, dMaxEError);
2941
2942 // X
2943 locHistName = "XErrorVsP";
2944 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{x}");
2945 dHistMap_ShowerXErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
2946
2947 locHistName = "XErrorVsTheta";
2948 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{x}");
2949 dHistMap_ShowerXErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
2950
2951 locHistName = "XErrorVsPhi";
2952 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{x}");
2953 dHistMap_ShowerXErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
2954
2955 // Y
2956 locHistName = "YErrorVsP";
2957 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{y}");
2958 dHistMap_ShowerYErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DXYErrorBins, 0.0, dMaxXYError);
2959
2960 locHistName = "YErrorVsTheta";
2961 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{y}");
2962 dHistMap_ShowerYErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DXYErrorBins, 0.0, dMaxXYError);
2963
2964 locHistName = "YErrorVsPhi";
2965 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{y}");
2966 dHistMap_ShowerYErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DXYErrorBins, 0.0, dMaxXYError);
2967
2968 // Z
2969 locHistName = "ZErrorVsP";
2970 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{z}");
2971 dHistMap_ShowerZErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DZErrorBins, 0.0, dMaxZError);
2972
2973 locHistName = "ZErrorVsTheta";
2974 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{z}");
2975 dHistMap_ShowerZErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DZErrorBins, 0.0, dMaxZError);
2976
2977 locHistName = "ZErrorVsPhi";
2978 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{z}");
2979 dHistMap_ShowerZErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DZErrorBins, 0.0, dMaxZError);
2980
2981 // E
2982 locHistName = "TErrorVsP";
2983 locHistTitle = locParticleROOTName + string(";p (GeV/c);#sigma_{E}");
2984 dHistMap_ShowerTErrorVsP[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, locMaxP, dNum2DTErrorBins, 0.0, dMaxTError);
2985
2986 locHistName = "TErrorVsTheta";
2987 locHistTitle = locParticleROOTName + string(";#theta#circ;#sigma_{E}");
2988 dHistMap_ShowerTErrorVsTheta[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, locMinTheta, locMaxTheta, dNum2DTErrorBins, 0.0, dMaxTError);
2989
2990 locHistName = "TErrorVsPhi";
2991 locHistTitle = locParticleROOTName + string(";#phi#circ;#sigma_{E}");
2992 dHistMap_ShowerTErrorVsPhi[locIsBCALFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DTErrorBins, 0.0, dMaxTError);
2993
2994 gDirectory(TDirectory::CurrentDirectory())->cd("..");
2995 }
2996
2997 //Return to the base directory
2998 ChangeTo_BaseDirectory();
2999 }
3000 japp->RootUnLock(); //RELEASE ROOT LOCK!!
3001}
3002
3003bool DHistogramAction_TrackShowerErrors::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3004{
3005 if(Get_NumPreviousParticleCombos() != 0)
3006 return true; //else double-counting!
3007
3008 vector<const DChargedTrack*> locPreSelectChargedTracks;
3009 locEventLoop->Get(locPreSelectChargedTracks, dTrackSelectionTag.c_str());
3010
3011 for(size_t loc_i = 0; loc_i < locPreSelectChargedTracks.size(); ++loc_i)
3012 {
3013 const DChargedTrackHypothesis* locChargedTrackHypothesis = locPreSelectChargedTracks[loc_i]->Get_BestFOM();
3014
3015 Particle_t locPID = (locChargedTrackHypothesis->dFOM < dMinPIDFOM) ? Unknown : locChargedTrackHypothesis->PID();
3016 if(dHistMap_TrackPxErrorVsP.find(locPID) == dHistMap_TrackPxErrorVsP.end())
3017 continue; //not interested in histogramming
3018
3019 DVector3 locMomentum = locChargedTrackHypothesis->momentum();
3020 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
3021 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
3022 double locP = locMomentum.Mag();
3023
3024 const DMatrixDSym& locCovarianceMatrix = locChargedTrackHypothesis->errorMatrix();
3025 double locPxError = sqrt(locCovarianceMatrix(0, 0));
3026 double locPyError = sqrt(locCovarianceMatrix(1, 1));
3027 double locPzError = sqrt(locCovarianceMatrix(2, 2));
3028 double locXError = sqrt(locCovarianceMatrix(3, 3));
3029 double locYError = sqrt(locCovarianceMatrix(4, 4));
3030 double locZError = sqrt(locCovarianceMatrix(5, 5));
3031
3032 //FILL HISTOGRAMS
3033 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3034 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3035 Lock_Action();
3036 {
3037 dHistMap_TrackPxErrorVsP[locPID]->Fill(locP, locPxError);
3038 dHistMap_TrackPxErrorVsTheta[locPID]->Fill(locTheta, locPxError);
3039 dHistMap_TrackPxErrorVsPhi[locPID]->Fill(locPhi, locPxError);
3040
3041 dHistMap_TrackPyErrorVsP[locPID]->Fill(locP, locPyError);
3042 dHistMap_TrackPyErrorVsTheta[locPID]->Fill(locTheta, locPyError);
3043 dHistMap_TrackPyErrorVsPhi[locPID]->Fill(locPhi, locPyError);
3044
3045 dHistMap_TrackPzErrorVsP[locPID]->Fill(locP, locPzError);
3046 dHistMap_TrackPzErrorVsTheta[locPID]->Fill(locTheta, locPzError);
3047 dHistMap_TrackPzErrorVsPhi[locPID]->Fill(locPhi, locPzError);
3048
3049 dHistMap_TrackXErrorVsP[locPID]->Fill(locP, locXError);
3050 dHistMap_TrackXErrorVsTheta[locPID]->Fill(locTheta, locXError);
3051 dHistMap_TrackXErrorVsPhi[locPID]->Fill(locPhi, locXError);
3052
3053 dHistMap_TrackYErrorVsP[locPID]->Fill(locP, locYError);
3054 dHistMap_TrackYErrorVsTheta[locPID]->Fill(locTheta, locYError);
3055 dHistMap_TrackYErrorVsPhi[locPID]->Fill(locPhi, locYError);
3056
3057 dHistMap_TrackZErrorVsP[locPID]->Fill(locP, locZError);
3058 dHistMap_TrackZErrorVsTheta[locPID]->Fill(locTheta, locZError);
3059 dHistMap_TrackZErrorVsPhi[locPID]->Fill(locPhi, locZError);
3060 }
3061 Unlock_Action();
3062 }
3063
3064 vector<const DNeutralParticle*> locNeutralParticles;
3065 locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());
3066
3067 for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
3068 {
3069 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(Gamma);
3070 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
3071 continue;
3072
3073 DVector3 locMomentum = locNeutralParticleHypothesis->momentum();
3074 double locPhi = locMomentum.Phi()*180.0/TMath::Pi();
3075 double locTheta = locMomentum.Theta()*180.0/TMath::Pi();
3076 double locP = locMomentum.Mag();
3077
3078 const DNeutralShower* locNeutralShower = locNeutralParticles[loc_i]->dNeutralShower;
3079 const DMatrixDSym& locCovarianceMatrix = locNeutralShower->dCovarianceMatrix;
3080 bool locIsBCALFlag = (locNeutralShower->dDetectorSystem == SYS_BCAL);
3081 double locEError = sqrt(locCovarianceMatrix(0, 0));
3082 double locXError = sqrt(locCovarianceMatrix(1, 1));
3083 double locYError = sqrt(locCovarianceMatrix(2, 2));
3084 double locZError = sqrt(locCovarianceMatrix(3, 3));
3085 double locTError = sqrt(locCovarianceMatrix(4, 4));
3086
3087 //FILL HISTOGRAMS
3088 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3089 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3090 Lock_Action();
3091 {
3092 dHistMap_ShowerEErrorVsP[locIsBCALFlag]->Fill(locP, locEError);
3093 dHistMap_ShowerEErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locEError);
3094 dHistMap_ShowerEErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locEError);
3095
3096 dHistMap_ShowerXErrorVsP[locIsBCALFlag]->Fill(locP, locXError);
3097 dHistMap_ShowerXErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locXError);
3098 dHistMap_ShowerXErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locXError);
3099
3100 dHistMap_ShowerYErrorVsP[locIsBCALFlag]->Fill(locP, locYError);
3101 dHistMap_ShowerYErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locYError);
3102 dHistMap_ShowerYErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locYError);
3103
3104 dHistMap_ShowerZErrorVsP[locIsBCALFlag]->Fill(locP, locZError);
3105 dHistMap_ShowerZErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locZError);
3106 dHistMap_ShowerZErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locZError);
3107
3108 dHistMap_ShowerTErrorVsP[locIsBCALFlag]->Fill(locP, locTError);
3109 dHistMap_ShowerTErrorVsTheta[locIsBCALFlag]->Fill(locTheta, locTError);
3110 dHistMap_ShowerTErrorVsPhi[locIsBCALFlag]->Fill(locPhi, locTError);
3111 }
3112 Unlock_Action();
3113 }
3114
3115 return true;
3116}
3117
3118void DHistogramAction_NumReconstructedObjects::Initialize(JEventLoop* locEventLoop)
3119{
3120 string locHistName;
3121
3122 bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
3123
3124 //CREATE THE HISTOGRAMS
3125 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
3126 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
3127 {
3128 CreateAndChangeTo_ActionDirectory();
3129
3130 //2D Summary
3131 locHistName = "NumHighLevelObjects";
3132 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)
3133 dHist_NumHighLevelObjects = static_cast<TH2D*>(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()));
3134 else
3135 {
3136 dHist_NumHighLevelObjects = new TH2D(locHistName.c_str(), ";;# Objects / Event", 13, 0.5, 13.5, dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3137 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(1, "DRFTime");
3138 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(2, "DSCHit");
3139 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(3, "DTOFPoint");
3140 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(4, "DBCALShower");
3141 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(5, "DFCALShower");
3142 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(6, "DTimeBasedTrack");
3143 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(7, "TrackSCMatches");
3144 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(8, "TrackTOFMatches");
3145 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(9, "TrackBCALMatches");
3146 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(10, "TrackFCALMatches");
3147 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(11, "DBeamPhoton");
3148 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(12, "DChargedTrack");
3149 dHist_NumHighLevelObjects->GetXaxis()->SetBinLabel(13, "DNeutralShower");
3150 }
3151
3152 //Charged
3153 locHistName = "NumChargedTracks";
3154 dHist_NumChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3155 locHistName = "NumPosChargedTracks";
3156 dHist_NumPosChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3157 locHistName = "NumNegChargedTracks";
3158 dHist_NumNegChargedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} DChargedTrack", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3159
3160 //TBT
3161 locHistName = "NumTimeBasedTracks";
3162 dHist_NumTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3163 locHistName = "NumPosTimeBasedTracks";
3164 dHist_NumPosTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3165 locHistName = "NumNegTimeBasedTracks";
3166 dHist_NumNegTimeBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Time-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3167
3168 if(!locIsRESTEvent)
3169 {
3170 //WBT
3171 locHistName = "NumWireBasedTracks";
3172 dHist_NumWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3173 locHistName = "NumPosWireBasedTracks";
3174 dHist_NumPosWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3175 locHistName = "NumNegWireBasedTracks";
3176 dHist_NumNegWireBasedTracks = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Wire-Based Tracks", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3177
3178 //Track Candidates
3179 locHistName = "NumTrackCandidates";
3180 dHist_NumTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3181 locHistName = "NumPosTrackCandidates";
3182 dHist_NumPosTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3183 locHistName = "NumNegTrackCandidates";
3184 dHist_NumNegTrackCandidates = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3185
3186 //CDC Track Candidates
3187 locHistName = "NumPosTrackCandidates_CDC";
3188 dHist_NumPosTrackCandidates_CDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} CDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3189 locHistName = "NumNegTrackCandidates_CDC";
3190 dHist_NumNegTrackCandidates_CDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} CDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3191
3192 //FDC Track Candidates
3193 locHistName = "NumPosTrackCandidates_FDC";
3194 dHist_NumPosTrackCandidates_FDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{+} FDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3195 locHistName = "NumNegTrackCandidates_FDC";
3196 dHist_NumNegTrackCandidates_FDC = GetOrCreate_Histogram<TH1D>(locHistName, ";# #it{q}^{-} FDC Track Candidates", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3197 }
3198
3199 //Beam Photons
3200 locHistName = "NumBeamPhotons";
3201 dHist_NumBeamPhotons = GetOrCreate_Histogram<TH1D>(locHistName, ";# DBeamPhoton", dMaxNumBeamPhotons + 1, -0.5, (float)dMaxNumBeamPhotons + 0.5);
3202
3203 //Showers / Neutrals / TOF / SC
3204 locHistName = "NumFCALShowers";
3205 dHist_NumFCALShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DFCALShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3206 locHistName = "NumBCALShowers";
3207 dHist_NumBCALShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DBCALShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3208 locHistName = "NumNeutralShowers";
3209 dHist_NumNeutralShowers = GetOrCreate_Histogram<TH1D>(locHistName, ";# DNeutralShower", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3210 locHistName = "NumTOFPoints";
3211 dHist_NumTOFPoints = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTOFPoint", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3212 locHistName = "NumSCHits";
3213 dHist_NumSCHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DSCHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3214
3215 if(!locIsRESTEvent)
3216 {
3217 locHistName = "NumTAGMHits";
3218 dHist_NumTAGMHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTAGMHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3219 locHistName = "NumTAGHHits";
3220 dHist_NumTAGHHits = GetOrCreate_Histogram<TH1D>(locHistName, ";# DTAGHHit", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3221 }
3222
3223 //Matches
3224 locHistName = "NumTrackBCALMatches";
3225 dHist_NumTrackBCALMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-BCAL Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3226 locHistName = "NumTrackFCALMatches";
3227 dHist_NumTrackFCALMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-FCAL Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3228 locHistName = "NumTrackTOFMatches";
3229 dHist_NumTrackTOFMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-TOF Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3230 locHistName = "NumTrackSCMatches";
3231 dHist_NumTrackSCMatches = GetOrCreate_Histogram<TH1D>(locHistName, ";# Track-SC Matches", dMaxNumMatchObjects + 1, -0.5, (float)dMaxNumMatchObjects + 0.5);
3232
3233 if(!locIsRESTEvent)
3234 {
3235 //Hits
3236 locHistName = "NumCDCHits";
3237 dHist_NumCDCHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DCDCHit", dMaxNumCDCHits + 1, -0.5, (float)dMaxNumCDCHits + 0.5);
3238 locHistName = "NumFDCWireHits";
3239 dHist_NumFDCWireHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# Wire DFDCHit", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
3240 locHistName = "NumFDCCathodeHits";
3241 dHist_NumFDCCathodeHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# Cathode DFDCHit", dMaxNumFDCHits/2 + 1, -0.5, (float)dMaxNumFDCHits - 0.5 + 2.0);
3242
3243 locHistName = "NumTOFHits";
3244 dHist_NumTOFHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DTOFHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
3245 locHistName = "NumBCALHits";
3246 dHist_NumBCALHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DBCALHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
3247 locHistName = "NumFCALHits";
3248 dHist_NumFCALHits = GetOrCreate_Histogram<TH1I>(locHistName, ";# DFCALHit", dMaxNumTOFCalorimeterHits + 1, -0.5, (float)dMaxNumTOFCalorimeterHits + 0.5);
3249
3250 locHistName = "NumRFSignals";
3251 dHist_NumRFSignals = GetOrCreate_Histogram<TH1I>(locHistName, ";# DRFDigiTime + # DRFTDCDigiTime", dMaxNumObjects + 1, -0.5, (float)dMaxNumObjects + 0.5);
3252 }
3253
3254 //Return to the base directory
3255 ChangeTo_BaseDirectory();
3256 }
3257 japp->RootUnLock(); //RELEASE ROOT LOCK!!
3258}
3259
3260bool DHistogramAction_NumReconstructedObjects::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3261{
3262 if(Get_NumPreviousParticleCombos() != 0)
3263 return true; //else double-counting!
3264
3265 bool locIsRESTEvent = locEventLoop->GetJEvent().GetStatusBit(kSTATUS_REST);
3266
3267 vector<const DTrackTimeBased*> locTrackTimeBasedVector;
3268 locEventLoop->Get(locTrackTimeBasedVector);
3269
3270 vector<const DBeamPhoton*> locBeamPhotons;
3271 locEventLoop->Get(locBeamPhotons);
3272
3273 vector<const DFCALShower*> locFCALShowers;
3274 locEventLoop->Get(locFCALShowers);
3275
3276 vector<const DChargedTrack*> locChargedTracks;
3277 locEventLoop->Get(locChargedTracks);
3278
3279 vector<const DBCALShower*> locBCALShowers;
3280 locEventLoop->Get(locBCALShowers);
3281
3282 vector<const DNeutralShower*> locNeutralShowers;
3283 locEventLoop->Get(locNeutralShowers);
3284
3285 vector<const DTOFPoint*> locTOFPoints;
3286 locEventLoop->Get(locTOFPoints);
3287
3288 vector<const DSCHit*> locSCHits;
3289 locEventLoop->Get(locSCHits);
3290
3291 vector<const DRFTime*> locRFTimes;
3292 locEventLoop->Get(locRFTimes);
3293
3294 const DDetectorMatches* locDetectorMatches = NULL__null;
3295 locEventLoop->GetSingle(locDetectorMatches);
3296
3297 //if not REST
3298 vector<const DTrackWireBased*> locTrackWireBasedVector;
3299 vector<const DTrackCandidate*> locTrackCandidates;
3300 vector<const DTrackCandidate*> locTrackCandidates_CDC;
3301 vector<const DTrackCandidate*> locTrackCandidates_FDC;
3302 vector<const DCDCHit*> locCDCHits;
3303 vector<const DFDCHit*> locFDCHits;
3304 vector<const DTOFHit*> locTOFHits;
3305 vector<const DBCALHit*> locBCALHits;
3306 vector<const DFCALHit*> locFCALHits;
3307 vector<const DTAGMHit*> locTAGMHits;
3308 vector<const DTAGHHit*> locTAGHHits;
3309 vector<const DRFDigiTime*> locRFDigiTimes;
3310 vector<const DRFTDCDigiTime*> locRFTDCDigiTimes;
3311
3312 size_t locNumFDCWireHits = 0, locNumFDCCathodeHits = 0;
3313 if(!locIsRESTEvent)
3314 {
3315 locEventLoop->Get(locTrackWireBasedVector);
3316 locEventLoop->Get(locTrackCandidates);
3317 locEventLoop->Get(locTrackCandidates_CDC, "CDC");
3318 locEventLoop->Get(locTrackCandidates_FDC, "FDCCathodes");
3319 locEventLoop->Get(locCDCHits);
3320 locEventLoop->Get(locFDCHits);
3321 locEventLoop->Get(locTOFHits);
3322 locEventLoop->Get(locBCALHits);
3323 locEventLoop->Get(locFCALHits);
3324 locEventLoop->Get(locTAGHHits);
3325 locEventLoop->Get(locTAGMHits);
3326 locEventLoop->Get(locRFDigiTimes);
3327 locEventLoop->Get(locRFTDCDigiTimes);
3328
3329 for(size_t loc_i = 0; loc_i < locFDCHits.size(); ++loc_i)
3330 {
3331 if(locFDCHits[loc_i]->type == DFDCHit::AnodeWire)
3332 ++locNumFDCWireHits;
3333 else
3334 ++locNumFDCCathodeHits;
3335 }
3336 }
3337
3338 //FILL HISTOGRAMS
3339 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3340 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3341 Lock_Action(); //ACQUIRE ROOT LOCK!!
3342 {
3343 //# High-Level Objects
3344 dHist_NumHighLevelObjects->Fill(1, (Double_t)locRFTimes.size());
3345 dHist_NumHighLevelObjects->Fill(2, (Double_t)locSCHits.size());
3346 dHist_NumHighLevelObjects->Fill(3, (Double_t)locTOFPoints.size());
3347 dHist_NumHighLevelObjects->Fill(4, (Double_t)locBCALShowers.size());
3348 dHist_NumHighLevelObjects->Fill(5, (Double_t)locFCALShowers.size());
3349 dHist_NumHighLevelObjects->Fill(6, (Double_t)locTrackTimeBasedVector.size());
3350 dHist_NumHighLevelObjects->Fill(7, (Double_t)locDetectorMatches->Get_NumTrackSCMatches());
3351 dHist_NumHighLevelObjects->Fill(8, (Double_t)locDetectorMatches->Get_NumTrackTOFMatches());
3352 dHist_NumHighLevelObjects->Fill(9, (Double_t)locDetectorMatches->Get_NumTrackBCALMatches());
3353 dHist_NumHighLevelObjects->Fill(10, (Double_t)locDetectorMatches->Get_NumTrackFCALMatches());
3354 dHist_NumHighLevelObjects->Fill(11, (Double_t)locBeamPhotons.size());
3355 dHist_NumHighLevelObjects->Fill(12, (Double_t)locChargedTracks.size());
3356 dHist_NumHighLevelObjects->Fill(13, (Double_t)locNeutralShowers.size());
3357
3358 //Charged
3359 unsigned int locNumPos = 0, locNumNeg = 0;
3360 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
3361 {
3362 if(ParticleCharge(locChargedTracks[loc_i]->Get_BestFOM()->PID()) > 0)
3363 ++locNumPos;
3364 else
3365 ++locNumNeg;
3366 }
3367 dHist_NumChargedTracks->Fill(locChargedTracks.size());
3368 dHist_NumPosChargedTracks->Fill(locNumPos);
3369 dHist_NumNegChargedTracks->Fill(locNumNeg);
3370
3371 //TBT
3372 locNumPos = 0; locNumNeg = 0;
3373 for(size_t loc_i = 0; loc_i < locTrackTimeBasedVector.size(); ++loc_i)
3374 {
3375 if(ParticleCharge(locTrackTimeBasedVector[loc_i]->PID()) > 0)
3376 ++locNumPos;
3377 else
3378 ++locNumNeg;
3379 }
3380 dHist_NumTimeBasedTracks->Fill(locTrackTimeBasedVector.size());
3381 dHist_NumPosTimeBasedTracks->Fill(locNumPos);
3382 dHist_NumNegTimeBasedTracks->Fill(locNumNeg);
3383
3384 if(!locIsRESTEvent)
3385 {
3386 //WBT
3387 locNumPos = 0; locNumNeg = 0;
3388 for(size_t loc_i = 0; loc_i < locTrackWireBasedVector.size(); ++loc_i)
3389 {
3390 if(ParticleCharge(locTrackWireBasedVector[loc_i]->PID()) > 0)
3391 ++locNumPos;
3392 else
3393 ++locNumNeg;
3394 }
3395 dHist_NumWireBasedTracks->Fill(locTrackWireBasedVector.size());
3396 dHist_NumPosWireBasedTracks->Fill(locNumPos);
3397 dHist_NumNegWireBasedTracks->Fill(locNumNeg);
3398
3399 //Candidates
3400 locNumPos = 0; locNumNeg = 0;
3401 for(size_t loc_i = 0; loc_i < locTrackCandidates.size(); ++loc_i)
3402 {
3403 if(locTrackCandidates[loc_i]->charge() > 0.0)
3404 ++locNumPos;
3405 else
3406 ++locNumNeg;
3407 }
3408 dHist_NumTrackCandidates->Fill(locTrackCandidates.size());
3409 dHist_NumPosTrackCandidates->Fill(locNumPos);
3410 dHist_NumNegTrackCandidates->Fill(locNumNeg);
3411
3412 //CDC Candidates
3413 locNumPos = 0; locNumNeg = 0;
3414 for(size_t loc_i = 0; loc_i < locTrackCandidates_CDC.size(); ++loc_i)
3415 {
3416 if(locTrackCandidates_CDC[loc_i]->charge() > 0.0)
3417 ++locNumPos;
3418 else
3419 ++locNumNeg;
3420 }
3421 dHist_NumPosTrackCandidates_CDC->Fill(locNumPos);
3422 dHist_NumNegTrackCandidates_CDC->Fill(locNumNeg);
3423
3424 //FDC Candidates
3425 locNumPos = 0; locNumNeg = 0;
3426 for(size_t loc_i = 0; loc_i < locTrackCandidates_FDC.size(); ++loc_i)
3427 {
3428 if(locTrackCandidates_FDC[loc_i]->charge() > 0.0)
3429 ++locNumPos;
3430 else
3431 ++locNumNeg;
3432 }
3433 dHist_NumPosTrackCandidates_FDC->Fill(locNumPos);
3434 dHist_NumNegTrackCandidates_FDC->Fill(locNumNeg);
3435 }
3436
3437 //Beam Photons
3438 dHist_NumBeamPhotons->Fill((Double_t)locBeamPhotons.size());
3439
3440 //Showers
3441 dHist_NumFCALShowers->Fill((Double_t)locFCALShowers.size());
3442 dHist_NumBCALShowers->Fill((Double_t)locBCALShowers.size());
3443 dHist_NumNeutralShowers->Fill((Double_t)locNeutralShowers.size());
3444
3445 //TOF & SC
3446 dHist_NumTOFPoints->Fill((Double_t)locTOFPoints.size());
3447 dHist_NumSCHits->Fill((Double_t)locSCHits.size());
3448
3449 //TAGGER
3450 if(!locIsRESTEvent)
3451 {
3452 dHist_NumTAGMHits->Fill((Double_t)locTAGMHits.size());
3453 dHist_NumTAGHHits->Fill((Double_t)locTAGHHits.size());
3454 }
3455
3456 //Matches
3457 dHist_NumTrackBCALMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackBCALMatches());
3458 dHist_NumTrackFCALMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackFCALMatches());
3459 dHist_NumTrackTOFMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackTOFMatches());
3460 dHist_NumTrackSCMatches->Fill((Double_t)locDetectorMatches->Get_NumTrackSCMatches());
3461
3462 //Hits
3463 if(!locIsRESTEvent)
3464 {
3465 dHist_NumCDCHits->Fill((Double_t)locCDCHits.size());
3466 dHist_NumFDCWireHits->Fill((Double_t)locNumFDCWireHits);
3467 dHist_NumFDCCathodeHits->Fill((Double_t)locNumFDCCathodeHits);
3468 dHist_NumTOFHits->Fill((Double_t)locTOFHits.size());
3469 dHist_NumBCALHits->Fill((Double_t)locBCALHits.size());
3470 dHist_NumFCALHits->Fill((Double_t)locFCALHits.size());
3471 dHist_NumRFSignals->Fill((Double_t)(locRFDigiTimes.size() + locRFTDCDigiTimes.size()));
3472 }
3473 }
3474 Unlock_Action(); //RELEASE ROOT LOCK!!
3475
3476 return true;
3477}
3478
3479void DHistogramAction_TrackMultiplicity::Initialize(JEventLoop* locEventLoop)
3480{
3481 string locTrackSelectionTag = "NotATag", locShowerSelectionTag = "NotATag";
3482 if(gPARMS->Exists("COMBO:TRACK_SELECT_TAG"))
3483 gPARMS->GetParameter("COMBO:TRACK_SELECT_TAG", locTrackSelectionTag);
3484 if(gPARMS->Exists("COMBO:SHOWER_SELECT_TAG"))
3485 gPARMS->GetParameter("COMBO:SHOWER_SELECT_TAG", locShowerSelectionTag);
3486
3487 //CREATE THE HISTOGRAMS
3488 //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
3489 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
3490 {
3491 //So: Default tag is "", User can set it to something else
3492 //In here, if tag is "", get from gparms, if not, leave it alone
3493 //If gparms value does not exist, set it to (and use) "PreSelect"
3494 if(dTrackSelectionTag == "NotATag")
3495 dTrackSelectionTag = (locTrackSelectionTag == "NotATag") ? "PreSelect" : locTrackSelectionTag;
3496 if(dShowerSelectionTag == "NotATag")
3497 dShowerSelectionTag = (locShowerSelectionTag == "NotATag") ? "PreSelect" : locShowerSelectionTag;
3498
3499 CreateAndChangeTo_ActionDirectory();
3500
3501 string locHistName("NumReconstructedParticles");
3502 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)
3503 dHist_NumReconstructedParticles = static_cast<TH2D*>(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()));
3504 else
3505 {
3506 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);
3507 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(1, "# Total");
3508 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(2, "# q != 0");
3509 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(3, "# q = 0");
3510 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(4, "# q = +");
3511 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(5, "# q = -");
3512 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3513 {
3514 string locLabelName = string("# ") + string(ParticleName_ROOT(dFinalStatePIDs[loc_i]));
3515 dHist_NumReconstructedParticles->GetXaxis()->SetBinLabel(6 + loc_i, locLabelName.c_str());
3516 }
3517 }
3518
3519 locHistName = "NumGoodReconstructedParticles";
3520 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)
3521 dHist_NumGoodReconstructedParticles = static_cast<TH2D*>(gDirectory(TDirectory::CurrentDirectory())->Get(locHistName.c_str()));
3522 else
3523 {
3524 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);
3525 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(1, "# Total");
3526 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(2, "# q != 0");
3527 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(3, "# q = 0");
3528 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(4, "# q = +");
3529 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(5, "# q = -");
3530 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3531 {
3532 string locLabelName = string("# ") + string(ParticleName_ROOT(dFinalStatePIDs[loc_i]));
3533 dHist_NumGoodReconstructedParticles->GetXaxis()->SetBinLabel(6 + loc_i, locLabelName.c_str());
3534 }
3535 }
3536
3537 //Return to the base directory
3538 ChangeTo_BaseDirectory();
3539 }
3540 japp->RootUnLock(); //RELEASE ROOT LOCK!!
3541}
3542
3543bool DHistogramAction_TrackMultiplicity::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
3544{
3545 if(Get_NumPreviousParticleCombos() != 0)
3546 return true; //else double-counting!
3547
3548 vector<const DChargedTrack*> locChargedTracks;
3549 locEventLoop->Get(locChargedTracks);
3550
3551 vector<const DChargedTrack*> locGoodChargedTracks;
3552 locEventLoop->Get(locGoodChargedTracks, dTrackSelectionTag.c_str());
3553
3554 // get #tracks by PID/q type
3555 size_t locNumPositiveTracks = 0, locNumNegativeTracks = 0;
3556 map<Particle_t, size_t> locNumTracksByPID;
3557 for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i)
3558 {
3559 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM();
3560 Particle_t locPID = locChargedTrackHypothesis->PID();
3561
3562 if(locChargedTrackHypothesis->charge() > 0.0)
3563 ++locNumPositiveTracks;
3564 else
3565 ++locNumNegativeTracks;
3566
3567 if(locNumTracksByPID.find(locPID) != locNumTracksByPID.end())
3568 ++locNumTracksByPID[locPID];
3569 else
3570 locNumTracksByPID[locPID] = 1;
3571 }
3572
3573 // get # good tracks by PID/q type
3574 size_t locNumGoodPositiveTracks = 0, locNumGoodNegativeTracks = 0;
3575 map<Particle_t, size_t> locNumGoodTracksByPID;
3576 for(size_t loc_i = 0; loc_i < locGoodChargedTracks.size(); ++loc_i)
3577 {
3578 const DChargedTrackHypothesis* locChargedTrackHypothesis = locGoodChargedTracks[loc_i]->Get_BestFOM();
3579 Particle_t locPID = locChargedTrackHypothesis->PID();
3580
3581 double locPIDFOM = locChargedTrackHypothesis->dFOM;
3582
3583 if(locChargedTrackHypothesis->charge() > 0.0)
3584 ++locNumGoodPositiveTracks;
3585 else
3586 ++locNumGoodNegativeTracks;
3587
3588 if(locPIDFOM < dMinPIDFOM)
3589 continue;
3590
3591 if(locNumGoodTracksByPID.find(locPID) != locNumGoodTracksByPID.end())
3592 ++locNumGoodTracksByPID[locPID];
3593 else
3594 locNumGoodTracksByPID[locPID] = 1;
3595 }
3596
3597 vector<const DNeutralParticle*> locNeutralParticles;
3598 locEventLoop->Get(locNeutralParticles);
3599
3600 vector<const DNeutralParticle*> locGoodNeutralParticles;
3601 locEventLoop->Get(locGoodNeutralParticles, dShowerSelectionTag.c_str());
3602
3603 // neutrals by pid
3604 for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i)
3605 {
3606 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM();
3607 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
3608 continue;
3609
3610 Particle_t locPID = locNeutralParticleHypothesis->PID();
3611 if(locNumTracksByPID.find(locPID) != locNumTracksByPID.end())
3612 ++locNumTracksByPID[locPID];
3613 else
3614 locNumTracksByPID[locPID] = 1;
3615 }
3616
3617 // good neutrals
3618 for(size_t loc_i = 0; loc_i < locGoodNeutralParticles.size(); ++loc_i)
3619 {
3620 const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locGoodNeutralParticles[loc_i]->Get_BestFOM();
3621 if(locNeutralParticleHypothesis->dFOM < dMinPIDFOM)
3622 continue;
3623
3624 Particle_t locPID = locNeutralParticleHypothesis->PID();
3625 if(locNumGoodTracksByPID.find(locPID) != locNumGoodTracksByPID.end())
3626 ++locNumGoodTracksByPID[locPID];
3627 else
3628 locNumGoodTracksByPID[locPID] = 1;
3629 }
3630
3631 size_t locNumGoodTracks = locNumGoodPositiveTracks + locNumGoodNegativeTracks;
3632
3633 //FILL HISTOGRAMS
3634 //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
3635 //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
3636 Lock_Action();
3637 {
3638 dHist_NumReconstructedParticles->Fill(0.0, (Double_t)(locChargedTracks.size() + locNeutralParticles.size()));
3639 dHist_NumReconstructedParticles->Fill(1.0, (Double_t)locChargedTracks.size());
3640 dHist_NumReconstructedParticles->Fill(2.0, (Double_t)locNeutralParticles.size());
3641 dHist_NumReconstructedParticles->Fill(3.0, (Double_t)locNumPositiveTracks);
3642 dHist_NumReconstructedParticles->Fill(4.0, (Double_t)locNumNegativeTracks);
3643 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3644 dHist_NumReconstructedParticles->Fill(5.0 + (Double_t)loc_i, (Double_t)locNumTracksByPID[dFinalStatePIDs[loc_i]]);
3645
3646 dHist_NumGoodReconstructedParticles->Fill(0.0, (Double_t)(locNumGoodTracks + locGoodNeutralParticles.size()));
3647 dHist_NumGoodReconstructedParticles->Fill(1.0, (Double_t)locNumGoodTracks);
3648 dHist_NumGoodReconstructedParticles->Fill(2.0, (Double_t)locGoodNeutralParticles.size());
3649 dHist_NumGoodReconstructedParticles->Fill(3.0, (Double_t)locNumGoodPositiveTracks);
3650 dHist_NumGoodReconstructedParticles->Fill(4.0, (Double_t)locNumGoodNegativeTracks);
3651 for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i)
3652 dHist_NumGoodReconstructedParticles->Fill(5.0 + (Double_t)loc_i, (Double_t)locNumGoodTracksByPID[dFinalStatePIDs[loc_i]]);
3653 }
3654 Unlock_Action();
3655
3656 return true;
3657}
3658