libraries/ANALYSIS/DHistogramActions

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