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