Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
p2pi_trees/DCustomAction_p2pi_unusedHists.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DCustomAction_p2pi_unusedHists.cc
4 // Created: Thu Jan 22 08:06:18 EST 2015
5 // Creator: jrsteven (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64)
6 //
7 
9 
10 void DCustomAction_p2pi_unusedHists::Initialize(JEventLoop* locEventLoop)
11 {
12 
13  // get PID algos
14  const DParticleID* locParticleID = NULL;
15  locEventLoop->GetSingle(locParticleID);
16  dParticleID = locParticleID;
17 
18  vector<const DFCALGeometry*> locFCALGeometry;
19  locEventLoop->Get(locFCALGeometry);
20  dFCALGeometry = locFCALGeometry[0];
21 
22  //CREATE THE HISTOGRAMS
23  //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock
24  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
25  {
26  //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action.
27  //If another thread has already created the folder, it just changes to it.
29 
30  dNShowerBCAL_FCAL = GetOrCreate_Histogram<TH2I>("NShowerBCAL_FCAL", "Number of non-matched showers; FCAL; BCAL", 10, 0, 10, 10, 0, 10);
31  dHistMap_FCALShowerDeltaR_P = GetOrCreate_Histogram<TH2I>("FCALShowerDeltaR_P", "Unused FCAL showers #DeltaR vs momentum; momentum; #DeltaR", 120, 0., 6., 200, 0, 50);
32  dHistMap_FCALShowerDeltaR_Theta = GetOrCreate_Histogram<TH2I>("FCALShowerDeltaR_Theta", "Unused FCAL showers #DeltaR vs #theta; #theta; #DeltaR", 150, 0, 15, 200, 0, 50);
33  dHistMap_FCALShowerDeltaD_P = GetOrCreate_Histogram<TH2I>("FCALShowerDeltaD_P", "Unused FCAL showers #DeltaD vs momentum; momentum; #DeltaR", 120, 0., 6., 200, 0, 50);
34  dHistMap_FCALShowerDeltaD_Theta = GetOrCreate_Histogram<TH2I>("FCALShowerDeltaD_Theta", "Unused FCAL showers #DeltaD vs #theta; #theta; #DeltaD", 150, 0, 15, 200, 0, 50);
35  dHistMap_FCALShowerDeltaD_DeltaT = GetOrCreate_Histogram<TH2I>("FCALShowerDeltaD_DeltaT", "Unused FCAL showers #DeltaD vs #Deltat; #Deltat; #DeltaD", 200, -100., 100., 200, 0, 50);
36 
37  dHistMap_BCALShowerDeltaPhi_DeltaZ = GetOrCreate_Histogram<TH2I>("BCALShowerDeltaPhi_DeltaZ", "Unused BCAL showers #Delta#phi vs #DeltaZ; #DeltaZ; #Delta#phi", 200, -50, 50, 360, -180, 180);
38  dHistMap_BCALShowerDeltaPhi_P = GetOrCreate_Histogram<TH2I>("BCALShowerDeltaPhi_P", "Unused BCAL showers #Delta#phi vs momentum; momentum; #Delta#phi", 120, 0., 6., 360, -180, 180);
39  dHistMap_BCALShowerDeltaPhi_Theta = GetOrCreate_Histogram<TH2I>("BCALShowerDeltaPhi_Theta", "Unused BCAL showers #Delta#phi vs #theta; #theta; #Delta#phi", 150, 0, 150, 360, -180, 180);
40  dHistMap_BCALShowerDeltaPhi_DeltaT = GetOrCreate_Histogram<TH2I>("BCALShowerDeltaPhi_DeltaT", "Unused BCAL showers #Delta#phi vs #Deltat; #Deltat; #Delta$phi", 200, -100., 100., 360, -180, 180);
41  dHistMap_BCALShowerDeltaZ_DeltaT = GetOrCreate_Histogram<TH2I>("BCALShowerDeltaZ_DeltaT", "Unused FCAL showers #DeltaZ vs #Deltat; #Deltat; #DeltaZ", 200, -100., 100., 200, -50, 50);
42 
43  int nbins = 100;
44  h1_deltaX = GetOrCreate_Histogram<TH1I>("h1_deltaX", "Fcal X - Track X", nbins,-25,25);
45  h1_deltaY = GetOrCreate_Histogram<TH1I>("h1_deltaY", "Fcal Y - Track Y", nbins,-25,25);
46  h1_Efcal = GetOrCreate_Histogram<TH1I>("h1_Efcal", "Hit energy", nbins,0,4);
47  h1_tfcal = GetOrCreate_Histogram<TH1I>("h1_tfcal", "Hit time", 250,-25,225);
48 
49  h1_N9 = GetOrCreate_Histogram<TH1I>("h1_N9", "Hit N9", 25, 0, 25);
50  h1_E9 = GetOrCreate_Histogram<TH1I>("h1_E9", "Energy E9 (GeV)", nbins, 0, 2);
51  h1_t9 = GetOrCreate_Histogram<TH1I>("h1_t9", "Time t9 (ns)", 250, -25, 225);
52  h1_t9sigma = GetOrCreate_Histogram<TH1I>("h1_t9sigma", "Time t9sigma",nbins,0,10);
53 
54  h1_N1 = GetOrCreate_Histogram<TH1I>("h1_N1", "Hit N1", 25, 0, 25);
55  h1_E1 = GetOrCreate_Histogram<TH1I>("h1_E1", "Energy E1 (GeV)", nbins, 0, 2);
56  h1_t1 = GetOrCreate_Histogram<TH1I>("h1_t1", "Time t1 (ns)", 250, -25, 225);
57  h1_t1sigma = GetOrCreate_Histogram<TH1I>("h1_t1sigma", "Time t1sigma",nbins,0,10);
58 
59  h2_YvsX9 = GetOrCreate_Histogram<TH2I>("h2_YvsX9", "Hit position Y vs X, E9",240,-120,120,240,-120,120);
60  h2_dEdx9_vsp= GetOrCreate_Histogram<TH2I>("h2_dEdx9_vsp", "Track dEdx9 vs p",nbins,0,4,nbins,0,10);
61  h2_E9vsp= GetOrCreate_Histogram<TH2I>("h2_E9vsp", "E9 vs p",nbins,0,4,nbins,0,4);
62  h2_dEdxvsE9= GetOrCreate_Histogram<TH2I>("h2_dEdxvsE9", "dEdx vs E9 energy",nbins,0,4,nbins,0,4);
63  h2_E9_vsTheta= GetOrCreate_Histogram<TH2I>("h2_E9_vsTheta", "E9 vs Theta",90,0,30,nbins,0,4);
64  h2_E9_vsPhi= GetOrCreate_Histogram<TH2I>("h2_E9_vsPhi", "E9 vs Phi",90,-180,180,nbins,0,4);
65 
66  h2_YvsX1 = GetOrCreate_Histogram<TH2I>("h2_YvsX1", "Hit position Y vs X, E1",240,-120,120,240,-120,120);
67  h2_dEdx1_vsp = GetOrCreate_Histogram<TH2I>("h2_dEdx1_vsp", "Track dEdx1 vs p",nbins,0,4,nbins,0,10);
68  h2_E1vsp = GetOrCreate_Histogram<TH2I>("h2_E1vsp", "E1 vs p",nbins,0,4,nbins,0,4);
69  h2_dEdxvsE1 = GetOrCreate_Histogram<TH2I>("h2_dEdxvsE1", "dEdx vs E1 energy",nbins,0,4,nbins,0,4);
70  h2_E1_vsTheta = GetOrCreate_Histogram<TH2I>("h2_E1_vsTheta", "E1 vs Theta",90,0,30,nbins,0,4);
71  h2_E1_vsPhi = GetOrCreate_Histogram<TH2I>("h2_E1_vsPhi", "E1 vs Phi",90,-180,180,nbins,0,4);
72  h2_E1vsRlocal = GetOrCreate_Histogram<TH2I>("h2_E1vsRlocal", "E1 vs Rtrk rel to Block center (cm)",nbins,0,5,nbins,0,4);
73  h2_YvsXcheck = GetOrCreate_Histogram<TH2I>("h2_E1vsRlocal", "E1 vs Rtrk rel to Block center (cm)",nbins,0,5,nbins,0,4);
74 
75  TString locHistName;
76  TString locHistTitle;
77  TString match[2] = {"Match",""};
78  map<int, TString> charge;
79  charge[-1] = "Neg"; charge[+1] = "Pos";
80 
81  for(int locDummy = 0; locDummy < 2; locDummy++) {
82  bool locMatch = (locDummy == 0);
83 
84  for(int locCharge = -1; locCharge < 2; locCharge+=2) {
85  locHistName = "TrackNhits_Theta" + match[locDummy] + charge[locCharge];
86  locHistTitle = "Number of hits on track vs #theta: " + match[locDummy] + charge[locCharge];
87  locHistTitle += "; #theta; # DOF";
88  dHistMap_TrackNhits_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90, 0., 90., 35, 0, 35);
89 
90  locHistName = "TrackNhitsCDC_Theta" + match[locDummy] + charge[locCharge];
91  locHistTitle = "Number of CDC ring hits on track vs #theta: " + match[locDummy] + charge[locCharge];
92  locHistTitle += "; #theta; # CDC ring hits";
93  dHistMap_TrackNhitsCDC_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90, 0., 90., 35, 0, 35);
94 
95  locHistName = "TrackNhitsFDC_Theta" + match[locDummy] + charge[locCharge];
96  locHistTitle = "Number of FDC plane hits on track vs #theta: " + match[locDummy] + charge[locCharge];
97  locHistTitle += "; #theta; # FDC plane hits";
98  dHistMap_TrackNhitsFDC_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90, 0., 90., 35, 0, 35);
99 
100  locHistName = "TrackNhitsFDCVsCDC_Theta13_16" + match[locDummy] + charge[locCharge];
101  locHistTitle = "Number of FDC plane hits vs CDC ring hits (13<#theta<16): " + match[locDummy] + charge[locCharge];
102  locHistTitle += "; # CDC ring hits; # FDC plane hits";
103  dHistMap_TrackNhitsFDCVsCDC_Theta13_16[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 35, 0, 35, 35, 0, 35);
104 
105  locHistName = "TrackNhitsFDCVsCDC_Theta16_20" + match[locDummy] + charge[locCharge];
106  locHistTitle = "Number of FDC plane hits vs CDC ring hits (16<#theta<20): " + match[locDummy] + charge[locCharge];
107  locHistTitle += "; # CDC ring hits; # FDC plane hits";
108  dHistMap_TrackNhitsFDCVsCDC_Theta16_20[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 35, 0, 35, 35, 0, 35);
109 
110  locHistName = "TrackChiSq_Theta" + match[locDummy] + charge[locCharge];
111  locHistTitle = "FOM for track #chi^{2} vs #theta: " + match[locDummy] + charge[locCharge];
112  locHistTitle += "; #theta; #chi^{2}";
113  dHistMap_TrackChiSq_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90, 0., 90., 500, 0., 100.);
114 
115  locHistName = "TrackFOM_Theta" + match[locDummy] + charge[locCharge];
116  locHistTitle = "FOM for track fit vs #theta: " + match[locDummy] + charge[locCharge];
117  locHistTitle += "; #theta; FOM";
118  dHistMap_TrackFOM_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90, 0., 90., 1000, 0., 1.);
119 
120  locHistName = "TrackP_Theta" + match[locDummy] + charge[locCharge];
121  locHistTitle = "P vs #theta: " + match[locDummy] + charge[locCharge];
122  locHistTitle += "; #theta; momentum";
123  dHistMap_TrackP_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90, 0., 90., 120, 0., 3.);
124 
125  locHistName = "TrackPOCAXY" + match[locDummy] + charge[locCharge];
126  locHistTitle = "POCA to beamline Y vs X: " + match[locDummy] + charge[locCharge];
127  locHistTitle += "; POCA X; POCAY";
128  dHistMap_TrackPOCAXY[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 200, -20., 20., 200, -20., 20.);
129 
130  locHistName = "TrackPOCAZ" + match[locDummy] + charge[locCharge];
131  locHistTitle = "POCA to beamline Z: " + match[locDummy] + charge[locCharge];
132  locHistTitle += "; POCA Z";
133  dHistMap_TrackPOCAZ[locMatch][locCharge] = GetOrCreate_Histogram<TH1I>(locHistName.Data(), locHistTitle.Data(), 500, 0., 250.);
134 
135  locHistName = "TrackCDCHitRadius_NHits" + match[locDummy] + charge[locCharge];
136  locHistTitle = "CDC Hit radius vs # hits: " + match[locDummy] + charge[locCharge];
137  locHistTitle += "; # hits; CDC Hit Radius";
138  dHistMap_TrackCDCHitRadius_Nhits[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 50, 0, 50, 50, 0, 50);
139 
140  locHistName = "HighNDFTrackCDCHitRadius_PT" + match[locDummy] + charge[locCharge];
141  locHistTitle = "High NDF tracks CDC Hit radius vs P_{T}: " + match[locDummy] + charge[locCharge];
142  locHistTitle += "; P_{T} (GeV); CDC Hit Radius";
143  dHistMap_HighNDFTrackCDCHitRadius_PT[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 120, 0., 3., 50, 0, 50);
144 
145  locHistName = "LowNDFTrackCDCHitRadius_PT" + match[locDummy] + charge[locCharge];
146  locHistTitle = "Low NDF tracks CDC Hit radius vs P_{T}: " + match[locDummy] + charge[locCharge];
147  locHistTitle += "; P_{T} (GeV); CDC Hit Radius";
148  dHistMap_LowNDFTrackCDCHitRadius_PT[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 120, 0., 3., 50, 0, 50);
149 
150  locHistName = "LowNDFTrackP_VertZ" + match[locDummy] + charge[locCharge];
151  locHistTitle = "Low NDF tracks P vs vertex Z: " + match[locDummy] + charge[locCharge];
152  locHistTitle += "; Vertex Z (cm); momentum";
153  dHistMap_LowNDFTrackP_VertZ[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 500, 0., 250., 120, 0., 3.);
154 
155  locHistName = "LowNDFTrackPT_Phi" + match[locDummy] + charge[locCharge];
156  locHistTitle = "Low NDF tracks P_{T} vs #phi: " + match[locDummy] + charge[locCharge];
157  locHistTitle += "; #phi; transverse momentum";
158  dHistMap_LowNDFTrackPT_Phi[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 360., -180., 180., 120, 0., 3.);
159 
160  locHistName = "TrackMCHitFraction_Theta" + match[locDummy] + charge[locCharge];
161  locHistTitle = "Track MC Hit Fraction vs #theta: " + match[locDummy] + charge[locCharge];
162  locHistTitle += "; #theta; MC Hit Fraction";
163  dHistMap_TrackMCHitFraction_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90., 0., 90., 100, 0., 1.);
164 
165  locHistName = "TrackMCMomRes_Theta" + match[locDummy] + charge[locCharge];
166  locHistTitle = "Track MC Momentum Resolution vs #theta: " + match[locDummy] + charge[locCharge];
167  locHistTitle += "; #theta; MC Momentum Resolution #Delta p/p";
168  dHistMap_TrackMCMomRes_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90., 0., 90., 100, -0.4, 0.4);
169 
170  locHistName = "TrackMCThetaRes_Theta" + match[locDummy] + charge[locCharge];
171  locHistTitle = "Track MC Theta Resolution vs #theta: " + match[locDummy] + charge[locCharge];
172  locHistTitle += "; #theta; MC Theta Resolution #Delta#theta (degrees)";
173  dHistMap_TrackMCThetaRes_Theta[locMatch][locCharge] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 90., 0., 90., 100, -5, 5);
174  }
175 
176  vector<DetectorSystem_t> locDetectorSystems;
177  locDetectorSystems.push_back(SYS_FCAL); locDetectorSystems.push_back(SYS_BCAL);
178  for(size_t loc_i = 0; loc_i < locDetectorSystems.size(); ++loc_i)
179  {
180  DetectorSystem_t locSystem = locDetectorSystems[loc_i];
181  TString locSystemName = (locSystem == SYS_FCAL) ? "FCAL" : "BCAL";
182  double thetaMax = (locSystem == SYS_FCAL) ? 15. : 150.;
183 
184  locHistName = "ShowerEnergy_Theta" + match[locDummy] + locSystemName;
185  locHistTitle = "Shower Energy vs #theta: " + match[locDummy] + locSystemName;
186  locHistTitle += ";#theta; Energy";
187  dHistMap_ShowerEnergy_Theta[locMatch][locSystem] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 100., 0., 5.);
188 
189  locHistName = "ShowerPhi_Theta" + match[locDummy] + locSystemName;
190  locHistTitle = "Shower #phi vs #theta: " + match[locDummy] + locSystemName;
191  locHistTitle += ";#theta; #phi";
192  dHistMap_ShowerPhi_Theta[locMatch][locSystem] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 360., -180., 180.);
193 
194  locHistName = "ShowerNClusters" + match[locDummy] + locSystemName;
195  locHistTitle = "Number of clusters in shower: " + match[locDummy] + locSystemName;
196  dHistMap_ShowerNclusters[locMatch][locSystem] = GetOrCreate_Histogram<TH1I>(locHistName.Data(), locHistTitle.Data(), 15, 0, 15);
197 
198  locHistName = "ShowerNHits" + match[locDummy] + locSystemName;
199  locHistTitle = "Number of hits in shower: " + match[locDummy] + locSystemName;
200  dHistMap_ShowerNhits[locMatch][locSystem] = GetOrCreate_Histogram<TH1I>(locHistName.Data(), locHistTitle.Data(), 15, 0, 15);
201 
202  locHistName = "ShowerMaxEnergy_NHits" + match[locDummy] + locSystemName;
203  locHistTitle = "E_{max}/E_{tot} vs number of hits in shower: " + match[locDummy] + locSystemName;
204  dHistMap_ShowerMaxEnergy_Nhits[locMatch][locSystem] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 15, 0, 15, 500, 0., 1.);
205 
206  locHistName = "ShowerDeltaT_NHits" + match[locDummy] + locSystemName;
207  locHistTitle = "t_{shower} - t_{#gamma} vs number of hits in shower: " + match[locDummy] + locSystemName;
208  dHistMap_ShowerDeltaT_Nhits[locMatch][locSystem] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 15, 0, 15, 200, -100., 100.);
209 
210  locHistName = "ShowerDeltaT_E" + match[locDummy] + locSystemName;
211  locHistTitle = "t_{shower} - t_{#gamma} vs shower energy: " + match[locDummy] + locSystemName;
212  dHistMap_ShowerDeltaT_E[locMatch][locSystem] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 100., 0., 5., 200, -100., 100.);
213 
214  locHistName = "ShowerE_Theta" + match[locDummy] + locSystemName;
215  locHistTitle = "t_{shower} - t_{#gamma} vs shower #theta: " + match[locDummy] + locSystemName;
216  dHistMap_ShowerE_Theta[locMatch][locSystem] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 120., 0., 6.);
217 
218  if(locSystem == SYS_BCAL){
219  locHistName = "ShowerNcell" + match[locDummy] + locSystemName;
220  locHistTitle = "Number of points in shower: " + match[locDummy] + locSystemName;
221  dHistMap_BCALShowerNcell[locMatch] = GetOrCreate_Histogram<TH1I>(locHistName.Data(), locHistTitle.Data(), 15, 0, 15);
222 
223  locHistName = "Layer1Energy_Theta" + match[locDummy] + locSystemName;
224  locHistTitle = "Shower Layer 1 Energy vs #theta: " + match[locDummy] + locSystemName;
225  locHistTitle += ";#theta; Energy Layer 1";
226  dHistMap_Layer1Energy_Theta[locMatch] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 500., 0., 1.);
227 
228  locHistName = "Layer2Energy_Theta" + match[locDummy] + locSystemName;
229  locHistTitle = "Shower Layer 2 Energy vs #theta: " + match[locDummy] + locSystemName;
230  locHistTitle += ";#theta; Energy Layer 2";
231  dHistMap_Layer2Energy_Theta[locMatch] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 500., 0., 1.);
232 
233  locHistName = "Layer3Energy_Theta" + match[locDummy] + locSystemName;
234  locHistTitle = "Shower Layer 3 Energy vs #theta: " + match[locDummy] + locSystemName;
235  locHistTitle += ";#theta; Energy Layer 3";
236  dHistMap_Layer3Energy_Theta[locMatch] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 500., 0., 1.);
237 
238  locHistName = "Layer4Energy_Theta" + match[locDummy] + locSystemName;
239  locHistTitle = "Shower Layer 4 Energy vs #theta: " + match[locDummy] + locSystemName;
240  locHistTitle += ";#theta; Energy Layer 4";
241  dHistMap_Layer4Energy_Theta[locMatch] = GetOrCreate_Histogram<TH2I>(locHistName.Data(), locHistTitle.Data(), 150, 0., thetaMax, 500., 0., 1.);
242  }
243  }
244  }
245  //Return to the base directory
247  }
248  japp->RootUnLock(); //RELEASE ROOT LOCK!!
249 }
250 
251 bool DCustomAction_p2pi_unusedHists::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
252 {
253 
254  // should only have one reaction step
255  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(0);
256 
257  // get beam photon energy and final state particles
258  const DKinematicData* locBeamPhoton = NULL;
259  deque<const DKinematicData*> locParticles;
260  if(!Get_UseKinFitResultsFlag()) { //measured
261  locBeamPhoton = locParticleComboStep->Get_InitialParticle_Measured();
262  locParticleComboStep->Get_FinalParticles_Measured(locParticles);
263  }
264  else {
265  locBeamPhoton = locParticleComboStep->Get_InitialParticle();
266  locParticleComboStep->Get_FinalParticles(locParticles);
267  }
268  double locBeamPhotonTime = locBeamPhoton->time();
269 
270  // detector matches for charged track -to- shower matching
271  const DDetectorMatches* locDetectorMatches = NULL;
272  locEventLoop->GetSingle(locDetectorMatches);
273 
274  // thrown particles for track resolution
275  vector<const DMCThrownMatching*> locMCThrownMatchingVector;
276  locEventLoop->Get(locMCThrownMatchingVector);
277  vector<const DMCThrown*> locMCThrowns;
278  locEventLoop->Get(locMCThrowns);
279 
280  // get all charged tracks and showers for comparison to combo (no "pre-select" cut here)
281  vector<const DChargedTrack*> locChargedTracks;
282  locEventLoop->Get(locChargedTracks); //, "PreSelect");
283  vector<const DNeutralShower*> locNeutralShowers;
284  locEventLoop->Get(locNeutralShowers); //, "PreSelect");
285  vector<const DFCALShower*> locFCALShowers;
286  locEventLoop->Get(locFCALShowers);
287  vector<const DBCALShower*> locBCALShowers;
288  locEventLoop->Get(locBCALShowers);
289 
290  vector<const DChargedTrack*> locUnusedChargedTracks;
291  vector<const DNeutralShower*> locUnusedNeutralShowers;
292 
293  // loop over charged tracks to make list of unused
294  for(size_t loc_j = 0; loc_j < locChargedTracks.size(); ++loc_j) {
295  int nMatched = 0;
296 
297  double locMatchFOM;
298  const DMCThrown* locMCThrown = NULL;
299  if(!locMCThrowns.empty())
300  locMCThrown = locMCThrownMatchingVector[0]->Get_MatchingMCThrown(locChargedTracks[loc_j], locMatchFOM);
301 
302  // add tracks not in combo to vector of unused tracks
303  for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) {
304 
305  const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i));
306  if(locChargedTrack == NULL) continue; // should never happen
307 
308  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_BestFOM();
309  if(locChargedTracks[loc_j]->Get_BestFOM()->candidateid == locChargedTrackHypothesis->candidateid) {
310  nMatched++;
311  FillTrack(locEventLoop, locChargedTracks[loc_j], true, locMCThrown);
312  }
313  }
314 
315  // plot properties of unused charged tracks
316  if(nMatched==0) {
317  locUnusedChargedTracks.push_back(locChargedTracks[loc_j]);
318  FillTrack(locEventLoop, locChargedTracks[loc_j], false, locMCThrown);
319  }
320  }
321 
322  // loop over BCAL and FCAL showers to fill histograms for matched tracks
323  for(size_t loc_j = 0; loc_j < locBCALShowers.size(); ++loc_j) {
324 
325  // add showers matched to tracks in combo
326  for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) {
327 
328  const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i));
329  if(locChargedTrack == NULL) continue; // should never happen
330 
331  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_BestFOM();
332  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
333  const DReferenceTrajectory* rt = locTrackTimeBased->rt;
334  if(rt == NULL) continue;
335 
336  const DBCALShower *locBCALShower = locBCALShowers[loc_j];
337  DVector3 bcal_pos(locBCALShower->x, locBCALShower->y, locBCALShower->z);
338 
339  double locFlightTime = 9.9E9, locPathLength = 9.9E9, locFlightTimeVariance = 9.9E9;
340  rt->DistToRTwithTime(bcal_pos, &locPathLength, &locFlightTime, &locFlightTimeVariance, SYS_BCAL);
341 
342  DVector3 proj_pos = rt->GetLastDOCAPoint();
343  if(proj_pos.Perp() < 65.0)
344  continue; // not inside BCAL!
345 
346  double dz = bcal_pos.z() - proj_pos.z();
347  double dphi = bcal_pos.Phi() - proj_pos.Phi();
348  while(dphi > TMath::Pi())
349  dphi -= 2*TMath::Pi();
350  while(dphi < -1.*TMath::Pi())
351  dphi += 2*TMath::Pi();
352 
353  double locDeltaT = locBCALShowers[loc_j]->t - locBeamPhotonTime;
354 
355  //FILL HISTOGRAMS
356  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
357  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
358  Lock_Action(); //ACQUIRE ROOT LOCK!!
359  {
360  dHistMap_BCALShowerDeltaPhi_DeltaZ->Fill(dz, dphi*180./TMath::Pi());
361  dHistMap_BCALShowerDeltaPhi_P->Fill(locTrackTimeBased->momentum().Mag(), dphi*180./TMath::Pi());
362  dHistMap_BCALShowerDeltaPhi_Theta->Fill(locTrackTimeBased->momentum().Theta()*180./TMath::Pi(), dphi*180./TMath::Pi());
363  dHistMap_BCALShowerDeltaPhi_DeltaT->Fill(locDeltaT, dphi*180./TMath::Pi());
364  dHistMap_BCALShowerDeltaZ_DeltaT->Fill(locDeltaT, dz);
365  }
366  Unlock_Action(); //RELEASE ROOT LOCK!!
367 
368  DBCALShowerMatchParams locBCALShowerMatchParams;
369  bool foundBCAL = dParticleID->Get_BestBCALMatchParams(locTrackTimeBased, locDetectorMatches, locBCALShowerMatchParams);
370  if(foundBCAL){
371  if(locBCALShowerMatchParams.dBCALShower == locBCALShower) {
372  DNeutralShower* locNeutralShower = new DNeutralShower();
373  locNeutralShower->dDetectorSystem = SYS_BCAL;
374  locNeutralShower->dEnergy = locBCALShowers[loc_j]->E;
375  locNeutralShower->dSpacetimeVertex.SetXYZT(locBCALShowers[loc_j]->x, locBCALShowers[loc_j]->y, locBCALShowers[loc_j]->z, locBCALShowers[loc_j]->t);
376  locNeutralShower->AddAssociatedObject(locBCALShowers[loc_j]);
377  double locFlightTime = locBCALShowerMatchParams.dFlightTime;
378  FillShower(locNeutralShower, true, locBeamPhotonTime, locFlightTime);
379  delete locNeutralShower;
380  }
381  }
382  }
383  }
384 
385  for(size_t loc_j = 0; loc_j < locFCALShowers.size(); ++loc_j) {
386 
387  // add showers matched to tracks in combo
388  for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) {
389 
390  const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i));
391  if(locChargedTrack == NULL) continue; // should never happen
392 
393  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_BestFOM();
394  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
395  const DReferenceTrajectory* rt = locTrackTimeBased->rt;
396  if(rt == NULL) continue;
397 
398  const DFCALShower *locFCALShower = locFCALShowers[loc_j];
399  DVector3 fcal_pos = locFCALShower->getPosition();
400  DVector3 norm(0.0, 0.0, 1.0); //normal vector for the FCAL plane
401  DVector3 proj_pos, proj_mom;
402  double locPathLength = 9.9E9, locFlightTime = 9.9E9, locFlightTimeVariance = 9.9e9;
403  if(rt->GetIntersectionWithPlane(fcal_pos, norm, proj_pos, proj_mom, &locPathLength, &locFlightTime, &locFlightTimeVariance, SYS_FCAL) != NOERROR)
404  continue;
405 
406  double dd = (fcal_pos - proj_pos).Mag();
407  double dr = (fcal_pos - proj_pos).Perp();
408  double dphi = fcal_pos.Phi() - proj_pos.Phi();
409  while(dphi > TMath::Pi())
410  dphi -= 2*TMath::Pi();
411  while(dphi < -1.*TMath::Pi())
412  dphi += 2*TMath::Pi();
413 
414  double locDeltaT = locFCALShowers[loc_j]->getTime() - locBeamPhotonTime;
415 
416  //FILL HISTOGRAMS
417  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
418  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
419  Lock_Action(); //ACQUIRE ROOT LOCK!!
420  {
421  dHistMap_FCALShowerDeltaR_P->Fill(locTrackTimeBased->momentum().Mag(), dr);
422  dHistMap_FCALShowerDeltaR_Theta->Fill(locTrackTimeBased->momentum().Theta()*180./TMath::Pi(), dr);
423  dHistMap_FCALShowerDeltaD_P->Fill(locTrackTimeBased->momentum().Mag(), dd);
424  dHistMap_FCALShowerDeltaD_Theta->Fill(locTrackTimeBased->momentum().Theta()*180./TMath::Pi(), dd);
425  dHistMap_FCALShowerDeltaD_DeltaT->Fill(locDeltaT, dd);
426  }
427  Unlock_Action(); //RELEASE ROOT LOCK!!
428 
429  DFCALShowerMatchParams locFCALShowerMatchParams;
430  bool foundFCAL = dParticleID->Get_BestFCALMatchParams(locTrackTimeBased, locDetectorMatches, locFCALShowerMatchParams);
431  if(foundFCAL){
432  if(locFCALShowerMatchParams.dFCALShower == locFCALShower) {
433  DNeutralShower* locNeutralShower = new DNeutralShower();
434  locNeutralShower->dDetectorSystem = SYS_FCAL;
435  locNeutralShower->dEnergy = locFCALShowers[loc_j]->getEnergy();
436  locNeutralShower->dSpacetimeVertex.SetVect(locFCALShowers[loc_j]->getPosition());
437  locNeutralShower->dSpacetimeVertex.SetT(locFCALShowers[loc_j]->getTime());
438  locNeutralShower->AddAssociatedObject(locFCALShowers[loc_j]);
439  double locFlightTime = locFCALShowerMatchParams.dFlightTime;
440  FillShower(locNeutralShower, true, locBeamPhotonTime, locFlightTime);
441  delete locNeutralShower;
442  }
443  }
444  }
445  }
446 
447  int nUnmatchedFCAL = 0;
448  int nUnmatchedBCAL = 0;
449 
450  /////////////////////////////////////////////////////////////////////////////////////////
451  // loop over neutral showers which should all be unused (unless bad matching criteria) //
452  /////////////////////////////////////////////////////////////////////////////////////////
453  for(size_t loc_j = 0; loc_j < locNeutralShowers.size(); ++loc_j) {
454  locUnusedNeutralShowers.push_back(locNeutralShowers[loc_j]);
455  double locFlightTime = locNeutralShowers[loc_j]->dSpacetimeVertex.Vect().Mag()/SPEED_OF_LIGHT;
456  FillShower(locNeutralShowers[loc_j], false, locBeamPhotonTime, locFlightTime);
457 
458  if(locNeutralShowers[loc_j]->dDetectorSystem == SYS_FCAL) {
459  nUnmatchedFCAL++;
460  }
461  if(locNeutralShowers[loc_j]->dDetectorSystem == SYS_BCAL) {
462  nUnmatchedBCAL++;
463  }
464  }
465 
466  //FILL HISTOGRAMS
467  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
468  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
469  Lock_Action(); //ACQUIRE ROOT LOCK!!
470  {
471  dNShowerBCAL_FCAL->Fill(nUnmatchedFCAL, nUnmatchedBCAL);
472  }
473  Unlock_Action(); //RELEASE ROOT LOCK!!
474 
475  return true; //return false if you want to use this action to apply a cut (and it fails the cut!)
476 }
477 
478 
479 void DCustomAction_p2pi_unusedHists::FillTrack(JEventLoop* locEventLoop, const DChargedTrack* locChargedTrack, bool locMatch, const DMCThrown* locMCThrown)
480 {
481 
482  const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_BestFOM();
483  int locCharge = locChargedTrackHypothesis->charge();
484 
485  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
486 
487  double nHits = locTrackTimeBased->Ndof + 5.;
488  double locTheta = locTrackTimeBased->momentum().Theta()*180/TMath::Pi();
489 
490  set<int> locCDCRings;
491  dParticleID->Get_CDCRings(locTrackTimeBased->dCDCRings, locCDCRings);
492 
493  set<int> locFDCPlanes;
494  dParticleID->Get_FDCPlanes(locTrackTimeBased->dFDCPlanes, locFDCPlanes);
495 
496  double locHitFraction = locTrackTimeBased->dNumHitsMatchedToThrown/nHits;
497  double locMomentumRes = -999;
498  double locThetaRes = -999;
499  if(locMCThrown) {
500  locMomentumRes = (locTrackTimeBased->momentum().Mag() - locMCThrown->momentum().Mag())/locMCThrown->momentum().Mag();
501  locThetaRes = (locTrackTimeBased->momentum().Theta() - locMCThrown->momentum().Theta())*180./TMath::Pi();
502  }
503 
504  //FILL HISTOGRAMS
505  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
506  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
507  Lock_Action(); //ACQUIRE ROOT LOCK!!
508  {
509  dHistMap_TrackNhits_Theta[locMatch][locCharge]->Fill(locTheta, nHits);
510  dHistMap_TrackNhitsCDC_Theta[locMatch][locCharge]->Fill(locTheta, locCDCRings.size());
511  dHistMap_TrackNhitsFDC_Theta[locMatch][locCharge]->Fill(locTheta, locFDCPlanes.size());
512  if(locTheta>13. && locTheta<16.)
513  dHistMap_TrackNhitsFDCVsCDC_Theta13_16[locMatch][locCharge]->Fill(locCDCRings.size(), locFDCPlanes.size());
514  if(locTheta>16. && locTheta<20.)
515  dHistMap_TrackNhitsFDCVsCDC_Theta16_20[locMatch][locCharge]->Fill(locCDCRings.size(), locFDCPlanes.size());
516  dHistMap_TrackChiSq_Theta[locMatch][locCharge]->Fill(locTheta, locTrackTimeBased->chisq);
517  dHistMap_TrackFOM_Theta[locMatch][locCharge]->Fill(locTheta, locTrackTimeBased->FOM);
518  dHistMap_TrackP_Theta[locMatch][locCharge]->Fill(locTheta, locTrackTimeBased->momentum().Mag());
519  dHistMap_TrackPOCAXY[locMatch][locCharge]->Fill(locTrackTimeBased->position().X(), locTrackTimeBased->position().Y());
520  dHistMap_TrackPOCAZ[locMatch][locCharge]->Fill(locTrackTimeBased->position().Z());
521 
522  dHistMap_TrackMCHitFraction_Theta[locMatch][locCharge]->Fill(locTrackTimeBased->momentum().Theta()*180.0/TMath::Pi(), locHitFraction);
523  if(locMCThrown) {
524  dHistMap_TrackMCMomRes_Theta[locMatch][locCharge]->Fill(locTheta, locMomentumRes);
525  dHistMap_TrackMCThetaRes_Theta[locMatch][locCharge]->Fill(locTheta, locThetaRes);
526  }
527 
528  for(set<int>::iterator locIterator = locFDCPlanes.begin(); locIterator != locFDCPlanes.end(); ++locIterator) {
529  //dHistMap_TrackFDCPlaneVs_Nhits->Fill(nHits, *locIterator);
530  }
531 
532  // understand low Nhits tracks
533  if(locTrackTimeBased->momentum().Theta()*180./TMath::Pi() > 30.){
534  for(set<int>::iterator locIterator = locCDCRings.begin(); locIterator != locCDCRings.end(); ++locIterator) {
535  dHistMap_TrackCDCHitRadius_Nhits[locMatch][locCharge]->Fill(nHits, *locIterator);
536  if(nHits < 20) dHistMap_LowNDFTrackCDCHitRadius_PT[locMatch][locCharge]->Fill(locTrackTimeBased->momentum().Perp(), *locIterator);
537  else dHistMap_HighNDFTrackCDCHitRadius_PT[locMatch][locCharge]->Fill(locTrackTimeBased->momentum().Perp(), *locIterator);
538  }
539 
540  if(nHits < 20){
541  dHistMap_LowNDFTrackP_VertZ[locMatch][locCharge]->Fill(locTrackTimeBased->position().Z(), locTrackTimeBased->momentum().Mag());
542  dHistMap_LowNDFTrackPT_Phi[locMatch][locCharge]->Fill(locTrackTimeBased->momentum().Phi()*180/TMath::Pi(), locTrackTimeBased->momentum().Perp());
543  }
544  }
545  }
546  Unlock_Action(); //RELEASE ROOT LOCK!!
547 
548  if(!locMatch) return;
549 
550  // do FCAL energy sums to compare with Elton's E9 and E25 in fcal_charged
551  Double_t zfcal=625;
552  DVector3 fcal_origin(0.0,0.0,zfcal);
553  DVector3 fcal_normal(0.0,0.0,1.0);
554  DVector3 trkpos(0.0,0.0,0.0);
555  DVector3 proj_mom(0.0,0.0,0.0);
556  double theta = locTrackTimeBased->momentum().Theta()*180./TMath::Pi();
557  double phi = locTrackTimeBased->momentum().Phi()*180./TMath::Pi();
558  double p = locTrackTimeBased->momentum().Mag();
559  double dEdx = locTrackTimeBased->dEdx()*1e6;
560 
561  vector<const DFCALHit*> fcalhits;
562  locEventLoop->Get(fcalhits);
563  if(fcalhits.empty()) return;
564 
565  if (locTrackTimeBased->rt->GetIntersectionWithPlane(fcal_origin,fcal_normal,trkpos,proj_mom,NULL,NULL,NULL,SYS_FCAL)==NOERROR){
566  double trkposX = trkpos.X();
567  double trkposY = trkpos.Y();
568  int trkrow = dFCALGeometry->row((float)trkposY);
569  int trkcol = dFCALGeometry->column((float)trkposX);
570 
571  // loop over fcal hits
572 
573  double E9=0; // energy, x, y of 9 blocks surrounding track
574  double E9peak=0;
575  double x9=0;
576  double y9=0;
577  double t9=0;
578  double t9sq=0;
579  double t9sigma=0;
580  int N9=0;
581  int Delta_block=1; // =1 for 3x3, =2 for 5x5
582  double dX_E1=-1000;
583  double dY_E1=-1000;
584 
585  int row_offset = 0; // offset actual row to look for randoms
586  int col_offset = 0;
587  trkrow += row_offset;
588  trkcol += col_offset;
589 
590  for (unsigned int j=0; j < fcalhits.size(); ++j) {
591  int row = fcalhits[j]->row;
592  int col = fcalhits[j]->column;
593  double x = fcalhits[j]->x;
594  double y = fcalhits[j]->y;
595  double Efcal = fcalhits[j]->E;
596  double tfcal= fcalhits[j]->t;
597  double intOverPeak = fcalhits[j]->intOverPeak;
598 
599  // fill histograms
600  int drow = abs(row - trkrow);
601  int dcol = abs(col - trkcol);
602 
603  h1_deltaX->Fill(x - trkposX);
604  h1_deltaY->Fill(y - trkposY);
605  h1_Efcal->Fill(Efcal);
606  h1_tfcal->Fill(tfcal);
607 
608  // select hits
609  if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50) && (intOverPeak>2.5 && intOverPeak<9)) {
610  E9 += Efcal;
611  E9peak += Efcal*6/intOverPeak; // factor of 6 so that E9peak ~ E9
612  x9 += x;
613  y9 += y;
614  t9 += tfcal;
615  t9sq += tfcal*tfcal;
616  N9 += 1;
617 
618  dX_E1 = x - trkposX;
619  dY_E1 = y - trkposY;
620  }
621 
622  } // end loop over fcal hits
623 
624 // x9 = N9>0? x9/N9 : 0;
625 // y9 = N9>0? y9/N9 : 0;
626  t9 = N9>0? t9/N9 : 0;
627  t9sigma = N9>0? sqrt(t9sq/N9 - t9*t9): 0;
628 
629  //FILL HISTOGRAMS
630  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
631  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
632  Lock_Action(); //ACQUIRE ROOT LOCK!!
633  {
634  if (N9>0) {
635  h1_N9->Fill(N9);
636  h1_E9->Fill(E9);
637  h1_t9->Fill(t9);
638  h1_t9sigma->Fill(t9sigma);
639  h2_YvsX9->Fill(trkposX,trkposY);
640  h2_dEdx9_vsp->Fill(p,dEdx);
641  h2_E9vsp->Fill(p,E9);
642  h2_dEdxvsE9->Fill(E9,dEdx);
643  h2_E9_vsTheta->Fill(theta,E9);
644  h2_E9_vsPhi->Fill(phi,E9);
645  }
646  if (N9==1) {
647  h1_N1->Fill(N9);
648  h1_E1->Fill(E9);
649  h1_t1->Fill(t9);
650  h1_t1sigma->Fill(t9sigma);
651  h2_YvsX1->Fill(trkposX,trkposY);
652  h2_dEdx1_vsp->Fill(p,dEdx);
653  h2_E1vsp->Fill(p,E9);
654  h2_dEdxvsE1->Fill(E9,dEdx);
655  h2_E1_vsTheta->Fill(theta,E9);
656  h2_E1_vsPhi->Fill(phi,E9);
657  double Rlocal = sqrt(dX_E1*dX_E1 + dY_E1*dY_E1);
658  h2_E1vsRlocal->Fill(Rlocal,E9);
659  h2_YvsXcheck->Fill(dX_E1,dY_E1);
660  }
661  }
662  Unlock_Action(); //RELEASE ROOT LOCK!!
663  }
664 
665  return;
666 }
667 
668 
669 void DCustomAction_p2pi_unusedHists::FillShower(const DNeutralShower* locNeutralShower, bool locMatch, double locBeamPhotonTime, double locFlightTime)
670 {
671 
672  int nClusters = 0;
673  int nHits = 0;
674 
675  double layerE[4] = {0., 0., 0., 0.};
676 
677  double locEnergy = locNeutralShower->dEnergy;
678  DVector3 locPosition = locNeutralShower->dSpacetimeVertex.Vect();
679  double locDeltaT = locNeutralShower->dSpacetimeVertex.T() - locFlightTime - locBeamPhotonTime;
680 
681  double locEnergyCluster = 0.;
682  double locMaxEnergyCluster = 0.;
683 
684  DetectorSystem_t locSystem = locNeutralShower->dDetectorSystem;
685  if(locSystem == SYS_FCAL) {
686 
687  const DFCALShower* locFCALShower = NULL;
688  locNeutralShower->GetSingleT(locFCALShower);
689 
690  vector<const DFCALCluster*> locFCALClusters;
691  locFCALShower->Get(locFCALClusters);
692  nClusters = locFCALClusters.size();
693 
694  // get hits in FCAL shower
695  for(unsigned int i=0; i<locFCALClusters.size(); i++){
696  const vector<DFCALCluster::DFCALClusterHit_t> locFCALHits = locFCALClusters[i]->GetHits();
697 
698  locEnergyCluster = locFCALClusters[i]->getEnergy();
699  locMaxEnergyCluster = locFCALClusters[i]->getEmax();
700 
701  for(unsigned int j=0; j<locFCALHits.size(); j++){
702 // const DFCALCluster::DFCALClusterHit_t hit = locFCALHits[j];
703  nHits++;
704  }
705  }
706  }
707  if(locSystem == SYS_BCAL) {
708 
709  const DBCALShower* locBCALShower = NULL;
710  locNeutralShower->GetSingleT(locBCALShower);
711 
712  //FILL HISTOGRAMS
713  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
714  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
715  Lock_Action(); //ACQUIRE ROOT LOCK!!
716  {
717  dHistMap_BCALShowerNcell[locMatch]->Fill(locBCALShower->N_cell);
718  }
719  Unlock_Action(); //RELEASE ROOT LOCK!!
720 
721  vector<const DBCALCluster*> locBCALClusters;
722  locBCALShower->Get(locBCALClusters);
723  nClusters = locBCALClusters.size();
724 
725  // get points in BCAL shower
726  for(unsigned int i=0; i<locBCALClusters.size(); i++){
727  vector<const DBCALPoint*> locBCALPoints;
728  locBCALClusters[i]->Get(locBCALPoints);
729 
730  locEnergyCluster = locBCALClusters[i]->E();
731  if(locBCALPoints.size() == 1) locMaxEnergyCluster = locBCALClusters[i]->E();
732 
733  for(unsigned int j=0; j<locBCALPoints.size(); j++){
734  const DBCALPoint *point = locBCALPoints[j];
735  nHits++;
736 
737  if(point->E() > locMaxEnergyCluster)
738  locMaxEnergyCluster = point->E();
739 
740  layerE[point->layer()-1] += point->E();
741  }
742  }
743  }
744 
745  //FILL HISTOGRAMS
746  //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock
747  //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex
748  Lock_Action(); //ACQUIRE ROOT LOCK!!
749  {
750  dHistMap_ShowerEnergy_Theta[locMatch][locSystem]->Fill(locPosition.Theta()*180./TMath::Pi(), locEnergy);
751  dHistMap_ShowerPhi_Theta[locMatch][locSystem]->Fill(locPosition.Theta()*180./TMath::Pi(), locPosition.Phi()*180./TMath::Pi());
752  dHistMap_ShowerNclusters[locMatch][locSystem]->Fill(nClusters);
753  dHistMap_ShowerNhits[locMatch][locSystem]->Fill(nHits);
754 
755  dHistMap_ShowerMaxEnergy_Nhits[locMatch][locSystem]->Fill(nHits, locMaxEnergyCluster/locEnergyCluster);
756  dHistMap_ShowerDeltaT_Nhits[locMatch][locSystem]->Fill(nHits, locDeltaT);
757  dHistMap_ShowerDeltaT_E[locMatch][locSystem]->Fill(locEnergy, locDeltaT);
758  if(locDeltaT > -5. && locDeltaT < 5. && locSystem == SYS_FCAL)
759  dHistMap_ShowerE_Theta[locMatch][locSystem]->Fill(locPosition.Theta()*180./TMath::Pi(), locEnergy);
760  if(locDeltaT > -5. && locDeltaT < 5. && locSystem == SYS_BCAL)
761  dHistMap_ShowerE_Theta[locMatch][locSystem]->Fill(locPosition.Theta()*180./TMath::Pi(), locEnergy);
762 
763  if(locSystem == SYS_BCAL) {
764  dHistMap_Layer1Energy_Theta[locMatch]->Fill(locPosition.Theta()*180./TMath::Pi(), layerE[0]/locEnergy);
765  dHistMap_Layer2Energy_Theta[locMatch]->Fill(locPosition.Theta()*180./TMath::Pi(), layerE[1]/locEnergy);
766  dHistMap_Layer3Energy_Theta[locMatch]->Fill(locPosition.Theta()*180./TMath::Pi(), layerE[2]/locEnergy);
767  dHistMap_Layer4Energy_Theta[locMatch]->Fill(locPosition.Theta()*180./TMath::Pi(), layerE[3]/locEnergy);
768  }
769 
770  }
771  Unlock_Action(); //RELEASE ROOT LOCK!!
772 
773  return;
774 }
vector< const DKinematicData * > Get_FinalParticles_Measured(void) const
TDirectoryFile * ChangeTo_BaseDirectory(void)
void Get_FDCPlanes(unsigned int locBitPattern, set< int > &locFDCPlanes) const
void Get_CDCRings(unsigned int locBitPattern, set< int > &locCDCRings) const
map< bool, map< DetectorSystem_t, TH2I * > > dHistMap_ShowerEnergy_Theta
map< bool, map< int, TH2I * > > dHistMap_TrackNhitsFDC_Theta
#define SPEED_OF_LIGHT
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
map< bool, map< DetectorSystem_t, TH1I * > > dHistMap_ShowerNclusters
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
#define y
DVector3 GetLastDOCAPoint(void) const
map< bool, map< DetectorSystem_t, TH2I * > > dHistMap_ShowerDeltaT_Nhits
map< bool, map< DetectorSystem_t, TH2I * > > dHistMap_ShowerMaxEnergy_Nhits
const DTrackTimeBased * Get_TrackTimeBased(void) const
TDirectoryFile * CreateAndChangeTo_ActionDirectory(void)
DetectorSystem_t
Definition: GlueX.h:15
const DKinematicData * Get_InitialParticle_Measured(void) const
map< bool, map< int, TH2I * > > dHistMap_TrackCDCHitRadius_Nhits
int layer() const
Definition: DBCALPoint.h:65
bool Get_UseKinFitResultsFlag(void) const
Definition: GlueX.h:19
map< bool, map< DetectorSystem_t, TH1I * > > dHistMap_ShowerNhits
JApplication * japp
DetectorSystem_t dDetectorSystem
DLorentzVector dSpacetimeVertex
map< bool, map< int, TH2I * > > dHistMap_TrackMCThetaRes_Theta
map< bool, map< int, TH2I * > > dHistMap_LowNDFTrackCDCHitRadius_PT
double time(void) const
float E() const
Definition: DBCALPoint.h:38
map< bool, map< int, TH2I * > > dHistMap_TrackMCMomRes_Theta
map< bool, map< int, TH2I * > > dHistMap_TrackNhitsCDC_Theta
bool Get_BestFCALMatchParams(const DTrackingData *locTrack, const DDetectorMatches *locDetectorMatches, shared_ptr< const DFCALShowerMatchParams > &locBestMatchParams) const
double DistToRTwithTime(DVector3 hit, double *s=NULL, double *t=NULL, double *var_t=NULL, DetectorSystem_t detector=SYS_NULL) const
bool Get_BestBCALMatchParams(const DTrackingData *locTrack, const DDetectorMatches *locDetectorMatches, shared_ptr< const DBCALShowerMatchParams > &locBestMatchParams) const
void FillShower(const DNeutralShower *locNeutralShower, bool locMatch, double locBeamPhotonTime, double locFlightTime)
Definition: GlueX.h:22
map< bool, map< int, TH2I * > > dHistMap_HighNDFTrackCDCHitRadius_PT
vector< const DKinematicData * > Get_FinalParticles(void) const
map< bool, map< int, TH2I * > > dHistMap_TrackNhitsFDCVsCDC_Theta16_20
map< bool, map< int, TH2I * > > dHistMap_TrackNhitsFDCVsCDC_Theta13_16
void FillTrack(JEventLoop *locEventLoop, const DChargedTrack *locChargedTrack, bool locMatch, const DMCThrown *locMCThrown)
const DFCALShower * dFCALShower
map< bool, map< int, TH2I * > > dHistMap_TrackChiSq_Theta
int Ndof
Number of degrees of freedom in the fit.
double charge(void) const
map< bool, map< int, TH2I * > > dHistMap_TrackMCHitFraction_Theta
int column(int channel) const
Definition: DFCALGeometry.h:65
double sqrt(double)
const JObject * Get_FinalParticle_SourceObject(size_t locFinalParticleIndex) const
map< bool, map< int, TH2I * > > dHistMap_TrackNhits_Theta
const DKinematicData * Get_InitialParticle(void) const
const DVector3 & momentum(void) const
map< bool, map< DetectorSystem_t, TH2I * > > dHistMap_ShowerE_Theta
const DChargedTrackHypothesis * Get_BestFOM(void) const
Definition: DChargedTrack.h:69
int row(int channel) const
Definition: DFCALGeometry.h:64
map< bool, map< int, TH2I * > > dHistMap_LowNDFTrackP_VertZ
const DParticleComboStep * Get_ParticleComboStep(size_t locStepIndex) const
map< bool, map< DetectorSystem_t, TH2I * > > dHistMap_ShowerPhi_Theta
DVector3 getPosition() const
Definition: DFCALShower.h:151
map< bool, map< DetectorSystem_t, TH2I * > > dHistMap_ShowerDeltaT_E
jerror_t GetIntersectionWithPlane(const DVector3 &origin, const DVector3 &norm, DVector3 &pos, double *s=NULL, double *t=NULL, double *var_t=NULL, DetectorSystem_t detector=SYS_NULL) const
const DBCALShower * dBCALShower
map< bool, map< int, TH2I * > > dHistMap_LowNDFTrackPT_Phi