File: | plugins/Analysis/fcal_charged/DEventProcessor_fcal_charged.cc |
Location: | line 518, column 6 |
Description: | Value stored to 'x9' is never read |
1 | // $Id$ |
2 | // |
3 | // File: DEventProcessor_fcal_charged.cc |
4 | // Modified version of BCAL_Eff |
5 | // |
6 | |
7 | #include "DEventProcessor_fcal_charged.h" |
8 | |
9 | #include <TLorentzVector.h> |
10 | #include "TMath.h" |
11 | |
12 | #include "DANA/DApplication.h" |
13 | #include "BCAL/DBCALShower.h" |
14 | #include "BCAL/DBCALTruthShower.h" |
15 | #include "BCAL/DBCALCluster.h" |
16 | #include "BCAL/DBCALPoint.h" |
17 | #include "BCAL/DBCALHit.h" |
18 | #include "FCAL/DFCALShower.h" |
19 | #include "FCAL/DFCALCluster.h" |
20 | #include "TRACKING/DMCThrown.h" |
21 | #include "ANALYSIS/DAnalysisUtilities.h" |
22 | //#include "TRACKING/DTrackFinder.h" |
23 | #include "GlueX.h" |
24 | |
25 | #include "BCAL/DBCALGeometry.h" |
26 | #include "FCAL/DFCALGeometry.h" |
27 | |
28 | #include "TOF/DTOFPoint.h" |
29 | #include "DVector3.h" |
30 | #include <vector> |
31 | #include <deque> |
32 | #include <string> |
33 | #include <iostream> |
34 | #include <algorithm> |
35 | #include <stdio.h> |
36 | #include <stdlib.h> |
37 | |
38 | static TH1I* h1_stats = NULL__null; |
39 | static TH1I* h1_deltaX = NULL__null; |
40 | static TH1I* h1_deltaY = NULL__null; |
41 | static TH1I* h1_deltaX_tof = NULL__null; |
42 | static TH1I* h1_deltaY_tof = NULL__null; |
43 | static TH1I* h1_Efcal = NULL__null; |
44 | static TH1I* h1_Ebcal = NULL__null; |
45 | static TH1I* h1_tfcal = NULL__null; |
46 | static TH1I* h1_intOverPeak = NULL__null; |
47 | |
48 | static TH1I* h1_N9 = NULL__null; |
49 | static TH1I* h1_E9 = NULL__null; |
50 | static TH1I* h1_t9 = NULL__null; |
51 | static TH1I* h1_t9sigma = NULL__null; |
52 | static TH1I* h1_N1 = NULL__null; |
53 | static TH1I* h1_E1 = NULL__null; |
54 | static TH1I* h1_t1 = NULL__null; |
55 | static TH1I* h1_t1sigma = NULL__null; |
56 | |
57 | static TH2I* h2_dEdx_vsp = NULL__null; |
58 | static TH2I* h2_dEdx9_vsp = NULL__null; |
59 | static TH2I* h2_YvsX9 = NULL__null; |
60 | static TH2I* h2_YvsXcheck = NULL__null; |
61 | static TH2I* h2_YvsX1 = NULL__null; |
62 | static TH2I* h2_E9vsp = NULL__null; |
63 | static TH2I* h2_E1vsp = NULL__null; |
64 | static TH2I* h2_E2vsp = NULL__null; |
65 | static TH2I* h2_E1vspPos = NULL__null; |
66 | static TH2I* h2_dEdxvsE9 = NULL__null; |
67 | static TH2I* h2_intOverPeak_vsEfcal = NULL__null; |
68 | static TH2I* h2_E1_vsintOverPeak = NULL__null; |
69 | static TH2I* h2_E1peak_vsp = NULL__null; |
70 | static TH2I* h2_E1vsRlocal = NULL__null; |
71 | static TH2I* h2_E9_vsTheta = NULL__null; |
72 | static TH2I* h2_E9_vsPhi = NULL__null; |
73 | static TH2I* h2_E1_vsTheta = NULL__null; |
74 | static TH2I* h2_E1_vsPhi = NULL__null; |
75 | static TH2D* h2_E9_vsThetaPhiw = NULL__null; |
76 | static TH2D* h2_E1_vsThetaPhiw = NULL__null; |
77 | static TH2D* h2_E1_vsThetaPhi = NULL__null; |
78 | |
79 | |
80 | |
81 | |
82 | // Routine used to create our DEventProcessor |
83 | |
84 | extern "C" |
85 | { |
86 | void InitPlugin(JApplication *locApplication) |
87 | { |
88 | InitJANAPlugin(locApplication); |
89 | locApplication->AddProcessor(new DEventProcessor_fcal_charged()); //register this plugin |
90 | } |
91 | } // "C" |
92 | |
93 | //------------------ |
94 | // init |
95 | //------------------ |
96 | jerror_t DEventProcessor_fcal_charged::init(void) |
97 | { |
98 | // First thread to get here makes all histograms. If one pointer is |
99 | // already not NULL, assume all histograms are defined and return now |
100 | if(h1_deltaX != NULL__null){ |
101 | printf ("TEST>> DEventProcessor_fcal_charged::init - Other threads continue\n"); |
102 | return NOERROR; |
103 | } |
104 | |
105 | |
106 | // ... create historgrams or trees ... |
107 | |
108 | // TDirectory *dir = new TDirectoryFile("FCAL","FCAL"); |
109 | // dir->cd(); |
110 | // create root folder for bcal and cd to it, store main dir |
111 | TDirectory *main = gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory ())); |
112 | gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory ()))->mkdir("fcal_charged")->cd(); |
113 | |
114 | |
115 | // double pi0_mass = 0.1349766; |
116 | int const nbins=100; |
117 | |
118 | h1_stats = new TH1I("h1_stats", "Run Statistics", 20,0,20); |
119 | |
120 | h1_deltaX = new TH1I("h1_deltaX", "Fcal X - Track X", nbins,-25,25); |
121 | h1_deltaX->SetXTitle("Fcal - Track X (cm)"); |
122 | h1_deltaX->SetYTitle("counts"); |
123 | h1_deltaY = new TH1I("h1_deltaY", "Fcal Y - Track X", nbins,-25,25); |
124 | h1_deltaY->SetXTitle("Fcal - Track Y (cm)"); |
125 | h1_deltaY->SetYTitle("counts"); |
126 | |
127 | h1_deltaX_tof = new TH1I("h1_deltaX_tof", "TOF X - Track X E1", nbins,-25,25); |
128 | h1_deltaX_tof->SetXTitle("TOF X - Track X (cm)"); |
129 | h1_deltaX_tof->SetYTitle("counts"); |
130 | h1_deltaY_tof = new TH1I("h1_deltaY_tof", "TOF Y - Track Y E1", nbins,-25,25); |
131 | h1_deltaY_tof->SetXTitle("TOF Y - Track Y (cm)"); |
132 | h1_deltaY_tof->SetYTitle("counts"); |
133 | |
134 | h1_Ebcal = new TH1I("h1_Ebcal", "Sum of point energy", nbins,0,1); |
135 | h1_Ebcal->SetXTitle("Bcal point energy (GeV)"); |
136 | h1_Ebcal->SetYTitle("counts"); |
137 | |
138 | h1_Efcal = new TH1I("h1_Efcal", "Hit energy", nbins,0,4); |
139 | h1_Efcal->SetXTitle("Fcal hit energy (GeV)"); |
140 | h1_Efcal->SetYTitle("counts"); |
141 | h1_tfcal = new TH1I("h1_tfcal", "Hit time", 250,-25,225); |
142 | h1_tfcal->SetXTitle("Fcal hit time (ns)"); |
143 | h1_tfcal->SetYTitle("counts"); |
144 | h1_intOverPeak = new TH1I("h1_intOverPeak", "Hit intOverPeak",nbins,0,25); |
145 | h1_intOverPeak->SetXTitle("Fcal hit intOverPeak"); |
146 | h1_intOverPeak->SetYTitle("counts"); |
147 | |
148 | h2_dEdx_vsp = new TH2I("h2_dEdx_vsp", "Track dEdx vs p",nbins,0,4,nbins,0,10); |
149 | h2_dEdx_vsp->SetXTitle("p (GeV/c)"); |
150 | h2_dEdx_vsp->SetYTitle("dEdx (keV/cm)"); |
151 | h2_dEdx9_vsp = new TH2I("h2_dEdx9_vsp", "Track dEdx9 vs p",nbins,0,4,nbins,0,10); |
152 | h2_dEdx9_vsp->SetXTitle("p (GeV/c)"); |
153 | h2_dEdx9_vsp->SetYTitle("dEdx (keV/cm)"); |
154 | |
155 | h1_N9 = new TH1I("h1_N9", "Hit N9",25,0,25); |
156 | h1_N9->SetXTitle("Number of Hits N9"); |
157 | h1_N9->SetYTitle("counts"); |
158 | h1_E9 = new TH1I("h1_E9", "Energy E9",nbins,0,2); |
159 | h1_E9->SetXTitle("Energy E9 (GeV)"); |
160 | h1_E9->SetYTitle("counts"); |
161 | h1_t9 = new TH1I("h1_t9", "Time t9",250,-25,225); |
162 | h1_t9->SetXTitle("Time t9 (ns)"); |
163 | h1_t9->SetYTitle("counts"); |
164 | h1_t9sigma = new TH1I("h1_t9sigma", "Time t9sigma",nbins,0,10); |
165 | h1_t9sigma->SetXTitle("Time spread t9sigma (ns)"); |
166 | h1_t9sigma->SetYTitle("counts"); |
167 | |
168 | h2_YvsX9 = new TH2I("h2_YvsX9", "Hit position Y vs X, E9",240,-120,120,240,-120,120); |
169 | h2_YvsX9->SetXTitle("X (cm)"); |
170 | h2_YvsX9->SetYTitle("Y (cm)"); |
171 | h2_YvsXcheck = new TH2I("h2_YvsXcheck", "Delta Hit Y vs X Checkered",nbins,-5,5,nbins,-5,5); |
172 | h2_YvsXcheck->SetXTitle("X (cm)"); |
173 | h2_YvsXcheck->SetYTitle("Y (cm)"); |
174 | h2_E1vsRlocal = new TH2I("h2_E1vsRlocal", "E1 vs Rtrk rel to Block center (cm)",nbins,0,5,nbins,0,4); |
175 | h2_E1vsRlocal->SetXTitle("Rlocal (cm)"); |
176 | h2_E1vsRlocal->SetYTitle("E1 (GeV)"); |
177 | h2_E9vsp = new TH2I("h2_E9vsp", "E9 vs p",nbins,0,4,nbins,0,4); |
178 | h2_E9vsp->SetXTitle("p (GeV)"); |
179 | h2_E9vsp->SetYTitle("E9 (GeV)"); |
180 | h2_dEdxvsE9 = new TH2I("h2_dEdxvsE9", "dEdx vs E9 energy",nbins,0,4,nbins,0,4); |
181 | h2_dEdxvsE9->SetXTitle("E9 (GeV)"); |
182 | h2_dEdxvsE9->SetYTitle("dEdx (keV/cm)"); |
183 | |
184 | h1_N1 = new TH1I("h1_N1", "Hit N1",25,0,25); |
185 | h1_N1->SetXTitle("Number of Hits N1"); |
186 | h1_N1->SetYTitle("counts"); |
187 | h1_E1 = new TH1I("h1_E1", "Energy E1",nbins,0,2); |
188 | h1_E1->SetXTitle("Energy E1 (GeV)"); |
189 | h1_E1->SetYTitle("counts"); |
190 | h1_t1 = new TH1I("h1_t1", "Time t1",250,-25,225); |
191 | h1_t1->SetXTitle("Time t1 (ns)"); |
192 | h1_t1->SetYTitle("counts"); |
193 | h1_t1sigma = new TH1I("h1_t1sigma", "Time t1sigma",nbins,0,10); |
194 | h1_t1sigma->SetXTitle("Time spread t1sigma (ns)"); |
195 | h1_t1sigma->SetYTitle("counts"); |
196 | |
197 | h2_YvsX1 = new TH2I("h2_YvsX1", "Hit position Y vs X, E1",240,-120,120,240,-120,120); |
198 | h2_YvsX1->SetXTitle("X (cm)"); |
199 | h2_YvsX1->SetYTitle("Y (cm)"); |
200 | h2_E1vsp = new TH2I("h2_E1vsp", "E1 vs p",nbins,0,4,nbins,0,4); |
201 | h2_E1vsp->SetXTitle("p (GeV)"); |
202 | h2_E1vsp->SetYTitle("E1 (GeV)"); |
203 | h2_E2vsp = new TH2I("h2_E2vsp", "E2 vs p",nbins,0,4,nbins,0,4); |
204 | h2_E2vsp->SetXTitle("p (GeV)"); |
205 | h2_E2vsp->SetYTitle("E2 (GeV)"); |
206 | h2_E1vspPos = new TH2I("h2_E1vspPos", "E1 vs p Q>0",nbins,0,4,nbins,0,4); |
207 | h2_E1vspPos->SetXTitle("p (GeV)"); |
208 | h2_E1vspPos->SetYTitle("E1 (GeV)"); |
209 | h2_E1peak_vsp = new TH2I("h2_E1peak_vsp", "E1peak*6 vs p",nbins,0,4,nbins,0,4); |
210 | h2_E1peak_vsp->SetXTitle("p (GeV)"); |
211 | h2_E1peak_vsp->SetYTitle("E1 peak*6(GeV)"); |
212 | h2_intOverPeak_vsEfcal = new TH2I("h2_intOverPeak_vsEfcal", "intOverPeak vs Efcal",nbins,0,1,nbins,0,25); |
213 | h2_intOverPeak_vsEfcal->SetXTitle("Efcal (GeV)"); |
214 | h2_intOverPeak_vsEfcal->SetYTitle("IntOverPeak"); |
215 | h2_E1_vsintOverPeak = new TH2I("h2_E1_vsintOverPeak", "E1 vs intOverPeak",nbins,0,25,nbins,0,4); |
216 | h2_E1_vsintOverPeak->SetXTitle("IntOverPeak"); |
217 | h2_E1_vsintOverPeak->SetYTitle("E1 (GeV)"); |
218 | |
219 | h2_E9_vsTheta = new TH2I("h2_E9_vsTheta", "E9 vs Theta",90,0,30,nbins,0,4); |
220 | h2_E9_vsTheta->SetXTitle("Theta (deg)"); |
221 | h2_E9_vsTheta->SetYTitle("E9 (GeV)"); |
222 | h2_E1_vsTheta = new TH2I("h2_E1_vsTheta", "E1 vs Theta",90,0,30,nbins,0,4); |
223 | h2_E1_vsTheta->SetXTitle("Theta (deg)"); |
224 | h2_E1_vsTheta->SetYTitle("E1 (GeV)"); |
225 | h2_E9_vsPhi = new TH2I("h2_E9_vsPhi", "E9 vs Phi",90,-180,180,nbins,0,4); |
226 | h2_E9_vsPhi->SetXTitle("Phi (deg)"); |
227 | h2_E9_vsPhi->SetYTitle("E9 (GeV)"); |
228 | h2_E1_vsPhi = new TH2I("h2_E1_vsPhi", "E1 vs Phi",90,-180,180,nbins,0,4); |
229 | h2_E1_vsPhi->SetXTitle("Phi (deg)"); |
230 | h2_E1_vsPhi->SetYTitle("E1 (GeV)"); |
231 | |
232 | h2_E9_vsThetaPhiw = new TH2D("h2_E9_vsThetaPhiw", "E9 vs ThetaPhiw",90,-180,180,90,0,30); |
233 | h2_E9_vsThetaPhiw->SetXTitle("Phi (deg)"); |
234 | h2_E9_vsThetaPhiw->SetYTitle("Theta (deg)"); |
235 | h2_E1_vsThetaPhi = new TH2D("h2_E1_vsThetaPhi", "Unity vs ThetaPhi",90,-180,180,90,0,30); |
236 | h2_E1_vsThetaPhi->SetXTitle("Phi (deg)"); |
237 | h2_E1_vsThetaPhi->SetYTitle("Theta (deg)"); |
238 | h2_E1_vsThetaPhiw = new TH2D("h2_E1_vsThetaPhiw", "E1 vs ThetaPhiw",90,-180,180,90,0,30); |
239 | h2_E1_vsThetaPhiw->SetXTitle("Phi (deg)"); |
240 | h2_E1_vsThetaPhiw->SetYTitle("Theta (deg)"); |
241 | |
242 | // back to main dir |
243 | printf ("TEST>> DEventProcessor_fcal_charged::init - First thread created histograms\n"); |
244 | main->cd(); |
245 | |
246 | return NOERROR; |
247 | } |
248 | |
249 | |
250 | //------------------ |
251 | // brun |
252 | //------------------ |
253 | jerror_t DEventProcessor_fcal_charged::brun(jana::JEventLoop* locEventLoop, int locRunNumber) |
254 | { |
255 | // This is called whenever the run number changes |
256 | //BCAL_Neutrals->Fill(); |
257 | //cout << " run number = " << RunNumber << endl; |
258 | /* |
259 | //Optional: Retrieve REST writer for writing out skims |
260 | locEventLoop->GetSingle(dEventWriterREST); |
261 | */ |
262 | |
263 | //vector<const DTrackFinder *> finders; |
264 | //locEventLoop->Get(finders); |
265 | //finder = const_cast<DTrackFinder*>(finders[0]); |
266 | |
267 | /* |
268 | //Recommeded: Create output ROOT TTrees (nothing is done if already created) |
269 | locEventLoop->GetSingle(dEventWriterROOT); |
270 | dEventWriterROOT->Create_DataTrees(locEventLoop); |
271 | */ |
272 | |
273 | return NOERROR; |
274 | } |
275 | |
276 | //------------------ |
277 | // evnt |
278 | //------------------ |
279 | |
280 | |
281 | jerror_t DEventProcessor_fcal_charged::evnt(jana::JEventLoop* locEventLoop, int locEventNumber) |
282 | { |
283 | |
284 | // This is called for every event. Use of common resources like writing |
285 | // to a file or filling a histogram should be mutex protected. Using |
286 | // locEventLoop->Get(...) to get reconstructed objects (and thereby activating the |
287 | // reconstruction algorithm) should be done outside of any mutex lock |
288 | // since multiple threads may call this method at the same time. |
289 | // |
290 | // Here's an example: |
291 | // |
292 | // vector<const MyDataClass*> mydataclasses; |
293 | // locEventLoop->Get(mydataclasses); |
294 | // |
295 | // japp->RootWriteLock(); |
296 | // ... fill historgrams or trees ... |
297 | // japp->RootUnLock(); |
298 | |
299 | // DOCUMENTATION: |
300 | // ANALYSIS library: https://halldweb1.jlab.org/wiki/index.php/GlueX_Analysis_Software |
301 | |
302 | vector<const DFCALShower*> locFCALShowers; |
303 | vector<const DBCALPoint*> bcalpoints; |
304 | vector<const DTOFPoint*> tofpoints; |
305 | vector<const DFCALHit*> fcalhits; |
306 | vector<const DFCALCluster*> locFCALClusters; |
307 | vector<const DVertex*> kinfitVertex; |
308 | //const DDetectorMatches* locDetectorMatches = NULL; |
309 | //locEventLoop->GetSingle(locDetectorMatches); |
310 | locEventLoop->Get(locFCALShowers); |
311 | locEventLoop->Get(bcalpoints);; |
312 | locEventLoop->Get(tofpoints); |
313 | locEventLoop->Get(fcalhits); |
314 | locEventLoop->Get(locFCALClusters); |
315 | locEventLoop->Get(kinfitVertex); |
316 | |
317 | vector<const DChargedTrack*> locChargedTrack; |
318 | locEventLoop->Get(locChargedTrack); |
319 | |
320 | const DFCALGeometry *fcalgeom=NULL__null; |
321 | locEventLoop->GetSingle(fcalgeom); |
322 | if (fcalgeom==NULL__null){ |
323 | jerr << "FCAL geometry not found!" << endl; |
324 | return RESOURCE_UNAVAILABLE; |
325 | } |
326 | |
327 | DVector3 trkpos(0.0,0.0,0.0); |
328 | DVector3 proj_mom(0.0,0.0,0.0); |
329 | DVector3 trkpos_tof(0.0,0.0,0.0); |
330 | DVector3 proj_mom_tof(0.0,0.0,0.0); |
331 | Double_t zfcal=638; |
332 | Double_t ztof=618.8; |
333 | DVector3 fcal_origin(0.0,0.0,zfcal); |
334 | DVector3 fcal_normal(0.0,0.0,1.0); |
335 | DVector3 tof_origin(0.0,0.0,ztof); |
336 | DVector3 tof_normal(0.0,0.0,1.0); |
337 | |
338 | /* Double_t Lfcal=45.; |
339 | Double_t zfcal_back=zfcal+Lfcal; |
340 | DVector3 trkpos_fcal_back(0.0,0.0,0.0);; |
341 | DVector3 proj_mom_back(0.0,0.0,0.0); |
342 | DVector3 fcal_origin_back(0.0,0.0,zfcal_back); |
343 | DVector3 fcal_normal_back(0.0,0.0,1.0);*/ |
344 | |
345 | if (locEventNumber%100 == 0) printf ("EventNumber=%d\n",locEventNumber); |
346 | |
347 | // FILL HISTOGRAMS |
348 | // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock |
349 | japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK |
350 | |
351 | double p; |
352 | double dEdx; |
353 | double pmin=1; |
354 | |
355 | |
356 | for (unsigned int i=0; i < locChargedTrack.size() ; ++i){ |
357 | const DChargedTrack *ctrk = locChargedTrack[i]; |
358 | const DChargedTrackHypothesis *bestctrack = ctrk->Get_BestFOM(); |
359 | vector<const DTrackTimeBased*> locTrackTimeBased; |
360 | bestctrack->Get(locTrackTimeBased); // locTrackTimeBased[0] contains the best FOM track |
361 | |
362 | double FOM = TMath::Prob(locTrackTimeBased[0]->chisq, locTrackTimeBased[0]->Ndof); |
363 | |
364 | // get intersection with fcal front face and backface; |
365 | vector<DTrackFitter::Extrapolation_t>extrapolations=locTrackTimeBased[0]->extrapolations.at(SYS_FCAL); |
366 | if (extrapolations.size()>0){ |
367 | trkpos=extrapolations[0].position; |
368 | h1_stats->Fill(0); |
369 | |
370 | vector<DTrackFitter::Extrapolation_t>tof_extrapolations=locTrackTimeBased[0]->extrapolations.at(SYS_TOF); |
371 | if (tof_extrapolations.size()==0) continue; |
372 | h1_stats->Fill(1); |
373 | trkpos_tof=tof_extrapolations[0].position; |
374 | |
375 | // sum up all energy and Bcal and keep only if likely BCAL trigger |
376 | |
377 | float bcal_energy = 0; |
378 | for (unsigned int j=0; j < bcalpoints.size(); ++j) { |
379 | bcal_energy += bcalpoints[j]->E(); |
380 | } |
381 | if (bcal_energy < 0.15) { |
382 | continue; // require that bcal be sufficient for trigger |
383 | } |
384 | h1_Ebcal->Fill(bcal_energy); |
385 | h1_stats->Fill(2); |
386 | |
387 | |
388 | double q = locTrackTimeBased[0]->charge(); |
389 | p = locTrackTimeBased[0]->momentum().Mag(); |
390 | double trkmass = locTrackTimeBased[0]->mass(); |
391 | |
392 | double radius = sqrt(trkpos.X()*trkpos.X() + trkpos.Y()*trkpos.Y()); // distance of track from beamline |
393 | dEdx = (locTrackTimeBased[0]->dNumHitsUsedFordEdx_CDC >= locTrackTimeBased[0]->dNumHitsUsedFordEdx_FDC) ? locTrackTimeBased[0]->ddEdx_CDC_amp : locTrackTimeBased[0]->ddEdx_FDC; |
394 | dEdx *= 1e6; // convert to keV |
395 | if(!(trkmass < 0.15 && FOM>0.01 && radius>20)) { |
396 | continue; // select tracks of interest |
397 | } |
398 | h1_stats->Fill(3); |
399 | |
400 | if (! (p>pmin)) { |
401 | continue; // require minimum momentum to pion analyzis |
402 | } |
403 | h1_stats->Fill (4); |
404 | |
405 | h2_dEdx_vsp->Fill(p,dEdx); |
406 | // use position at the back of fcal // asume zero field out here |
407 | // double xfcal_back = trkpos.X() + proj_mom.X()*Lfcal/proj_mom.Z(); |
408 | // double yfcal_back = trkpos.Y() + proj_mom.Y()*Lfcal/proj_mom.Z(); |
409 | double trkposX = trkpos.X(); |
410 | double trkposY = trkpos.Y(); |
411 | int trkrow = fcalgeom->row((float)trkposY); |
412 | int trkcol = fcalgeom->column((float)trkposX); |
413 | double dX_tof = 10000; |
414 | double dY_tof = 10000; |
415 | double dR_tof = 10000; |
416 | |
417 | // get tof matches |
418 | |
419 | for (unsigned int j=0; j < tofpoints.size(); ++j) { |
420 | double tofx = 10000; |
421 | double tofy = 10000; |
422 | if (tofpoints[j]->Is_XPositionWellDefined() && tofpoints[j]->Is_YPositionWellDefined()) { |
423 | DVector3 pos_tof = tofpoints[j]->pos; |
424 | tofx = pos_tof.X(); |
425 | tofy = pos_tof.Y(); |
426 | printf ("Event=%d tofx=%f, tofy=%f, trkx=%f, trky=%f \n",locEventNumber,tofx,tofy,trkposX,trkposY); |
427 | double dX_tof1 = tofx - trkpos_tof.X(); |
428 | double dY_tof1 = tofy - trkpos_tof.Y(); |
429 | if (sqrt(dX_tof1*dX_tof1 + dY_tof1*dY_tof1) < dR_tof) { |
430 | dX_tof = dX_tof1; |
431 | dY_tof = dY_tof1; |
432 | dR_tof = sqrt(dX_tof*dX_tof + dY_tof*dY_tof); |
433 | } |
434 | } |
435 | } |
436 | |
437 | printf ("\n1 Event=%d dX_tof=%f, DY_tof=%f, trkx=%f, trky=%f \n",locEventNumber,dX_tof,dY_tof,trkposX,trkposY); |
438 | if ( !(abs(dX_tof)<10 && abs(dY_tof)<10) ) continue; |
439 | |
440 | h1_stats->Fill(5); |
441 | printf ("2 Event=%d dX_tof=%f, DY_tof=%f, trkx=%f, trky=%f \n",locEventNumber,dX_tof,dY_tof,trkposX,trkposY); |
442 | printf ("Event=%d trkmass=%f radius=%f p=%f dEdx=%g,trkrow=%d trkcol=%d x=%f y=%f z=%f\n",locEventNumber,trkmass,radius,p,dEdx,trkrow,trkcol,trkposX,trkposY,trkpos.Z()); |
443 | |
444 | // get components of p at FCAL |
445 | // double theta = proj_mom.Theta()*TVector3::RAD2DEG; |
446 | double theta = proj_mom.Theta()*180./M_PI3.14159265358979323846; |
447 | double phi = proj_mom.Phi()*180./M_PI3.14159265358979323846; |
448 | printf ("Event=%d p=%f,theta=%f, phi=%f\n",locEventNumber,p,theta,phi); |
449 | |
450 | |
451 | // loop over fcal hits |
452 | |
453 | double E9=0; // energy, x, y of 9 blocks surrounding track |
454 | double E9peak=0; |
455 | double x9=0; |
456 | double y9=0; |
457 | double t9=0; |
458 | double t9sq=0; |
459 | double t9sigma=0; |
460 | int N9=0; |
461 | int Delta_block=2; // =1 for 3x3, =2 for 5x5 |
462 | //int row_E1=-1000; |
463 | //int col_E1=-1000; |
464 | //int drow_E1=1000; |
465 | //int dcol_E1=1000; |
466 | double dX_E1=-1000; |
467 | double dY_E1=-1000; |
468 | |
469 | int row_offset = 0; // offset actual row to look for randoms |
470 | int col_offset = 0; |
471 | trkrow += row_offset; |
472 | trkcol += col_offset; |
473 | |
474 | for (unsigned int j=0; j < fcalhits.size(); ++j) { |
475 | int row = fcalhits[j]->row; |
476 | int col = fcalhits[j]->column; |
477 | double x = fcalhits[j]->x; |
478 | double y = fcalhits[j]->y; |
479 | double Efcal = fcalhits[j]->E; |
480 | double tfcal= fcalhits[j]->t; |
481 | double intOverPeak = fcalhits[j]->intOverPeak; |
482 | // printf ("Event=%d row=%d col=%d x=%f y=%f Efcal=%f tfcal=%f intOverPeak=%f\n",locEventNumber,row,col,x,y,Efcal,tfcal,intOverPeak); |
483 | |
484 | // fill histograms |
485 | int drow = abs(row - trkrow); |
486 | int dcol = abs(col - trkcol); |
487 | |
488 | h1_deltaX->Fill(x - trkposX); |
489 | h1_deltaY->Fill(y - trkposY); |
490 | h1_Efcal->Fill(Efcal); |
491 | h1_tfcal->Fill(tfcal); |
492 | h1_intOverPeak->Fill(intOverPeak); |
493 | h2_intOverPeak_vsEfcal->Fill(Efcal,intOverPeak); |
494 | |
495 | // select hits |
496 | if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50) && (intOverPeak>2.5 && intOverPeak<9)) { |
497 | // if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50)) { |
498 | E9 += Efcal; |
499 | E9peak += Efcal*6/intOverPeak; // factor of 6 so that E9peak ~ E9 |
500 | x9 += x; |
501 | y9 += y; |
502 | t9 += tfcal; |
503 | t9sq += tfcal*tfcal; |
504 | N9 += 1; |
505 | |
506 | //save for later |
507 | //row_E1 = row; |
508 | //col_E1 = col; |
509 | //drow_E1 = drow; |
510 | //dcol_E1 = dcol; |
511 | dX_E1 = x - trkposX; |
512 | dY_E1 = y - trkposY; |
513 | printf ("Event=%d, trkposX=%f, trkposY=%f, dX_E1=%f, dY_E1=%f\n",locEventNumber,trkposX,trkposY,dX_E1,dY_E1); |
514 | } |
515 | |
516 | } // end loop over fcal hits |
517 | |
518 | x9 = N9>0? x9/N9 : 0; |
Value stored to 'x9' is never read | |
519 | y9 = N9>0? y9/N9 : 0; |
520 | t9 = N9>0? t9/N9 : 0; |
521 | t9sigma = N9>0? sqrt(t9sq/N9 - t9*t9): 0; |
522 | // printf ("\nEvent=%d N9=%d E9=%f x9=%f y9=%f t9=%f t9sigma=%f\n",locEventNumber,N9,E9,x9,y9,t9,t9sigma); |
523 | |
524 | |
525 | h1_stats->Fill(6); |
526 | if (E9 > 0.8 && theta<4) continue; |
527 | h1_stats->Fill(7); |
528 | if (N9>0) { |
529 | h1_N9->Fill(N9); |
530 | h1_E9->Fill(E9); |
531 | h1_t9->Fill(t9); |
532 | h1_t9sigma->Fill(t9sigma); |
533 | h2_YvsX9->Fill(trkposX,trkposY); |
534 | h2_dEdx9_vsp->Fill(p,dEdx); |
535 | h2_E9vsp->Fill(p,E9); |
536 | h2_dEdxvsE9->Fill(E9,dEdx); |
537 | h1_stats->Fill(8); |
538 | h2_E9_vsThetaPhiw->Fill(phi,theta,E9); |
539 | h2_E9_vsTheta->Fill(theta,E9); |
540 | h2_E9_vsPhi->Fill(phi,E9); |
541 | } |
542 | if (N9==1) { |
543 | h1_N1->Fill(N9); |
544 | h1_E1->Fill(E9); |
545 | h1_t1->Fill(t9); |
546 | h1_t1sigma->Fill(t9sigma); |
547 | h2_YvsX1->Fill(trkposX,trkposY); |
548 | h2_E1vsp->Fill(p,E9); |
549 | if (q > 0) h2_E1vspPos->Fill(p,E9); |
550 | h2_E1peak_vsp->Fill(p,E9peak); |
551 | h2_E1_vsintOverPeak->Fill(E9*6/E9peak,E9); // note that this only works because a single cell fires |
552 | h1_deltaX_tof->Fill(dX_tof); |
553 | h1_deltaY_tof->Fill(dY_tof); |
554 | h1_stats->Fill(9); |
555 | h2_E1_vsThetaPhi->Fill(phi,theta); |
556 | h2_E1_vsThetaPhiw->Fill(phi,theta,E9); |
557 | h2_E1_vsTheta->Fill(theta,E9); |
558 | h2_E1_vsPhi->Fill(phi,E9); |
559 | double Rlocal = sqrt(dX_E1*dX_E1 + dY_E1*dY_E1); |
560 | h2_E1vsRlocal->Fill(Rlocal,E9); |
561 | // if (((row_E1%2==0 && col_E1%2==0) || (row_E1%2==1 && col_E1%2==1)) ) h2_YvsXcheck->Fill(dX_E1,dY_E1); |
562 | h2_YvsXcheck->Fill(dX_E1,dY_E1); |
563 | } |
564 | if (N9==1 || N9==2) { |
565 | h2_E2vsp->Fill(p,E9); |
566 | h1_stats->Fill(10); |
567 | } |
568 | |
569 | |
570 | |
571 | } // end track intersection if statement |
572 | } |
573 | |
574 | japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK |
575 | |
576 | |
577 | /* |
578 | //Optional: Save event to output REST file. Use this to create skims. |
579 | dEventWriterREST->Write_RESTEvent(locEventLoop, "FCAL_Shower"); //string is part of output file name |
580 | */ |
581 | |
582 | return NOERROR; |
583 | } |
584 | |
585 | //------------------ |
586 | // erun |
587 | //------------------ |
588 | jerror_t DEventProcessor_fcal_charged::erun(void) |
589 | { |
590 | // This is called whenever the run number changes, before it is |
591 | // changed to give you a chance to clean up before processing |
592 | // events from the next run number. |
593 | |
594 | |
595 | return NOERROR; |
596 | } |
597 | |
598 | //------------------ |
599 | // fini |
600 | //------------------ |
601 | jerror_t DEventProcessor_fcal_charged::fini(void) |
602 | { |
603 | // Called before program exit after event processing is finished. |
604 | |
605 | return NOERROR; |
606 | } |
607 |