1 | |
2 | |
3 | |
4 | |
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 | |
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 | DFCALGeometry fcalgeom; |
39 | |
40 | |
41 | static TH1I* h1_stats = NULL__null; |
42 | static TH1I* h1_deltaX = NULL__null; |
43 | static TH1I* h1_deltaY = NULL__null; |
44 | static TH1I* h1_deltaX_tof = NULL__null; |
45 | static TH1I* h1_deltaY_tof = NULL__null; |
46 | static TH1I* h1_Efcal = NULL__null; |
47 | static TH1I* h1_Ebcal = NULL__null; |
48 | static TH1I* h1_tfcal = NULL__null; |
49 | static TH1I* h1_intOverPeak = NULL__null; |
50 | |
51 | static TH1I* h1_N9 = NULL__null; |
52 | static TH1I* h1_E9 = NULL__null; |
53 | static TH1I* h1_t9 = NULL__null; |
54 | static TH1I* h1_t9sigma = NULL__null; |
55 | static TH1I* h1_N1 = NULL__null; |
56 | static TH1I* h1_E1 = NULL__null; |
57 | static TH1I* h1_t1 = NULL__null; |
58 | static TH1I* h1_t1sigma = NULL__null; |
59 | |
60 | static TH2I* h2_dEdx_vsp = NULL__null; |
61 | static TH2I* h2_dEdx9_vsp = NULL__null; |
62 | static TH2I* h2_YvsX9 = NULL__null; |
63 | static TH2I* h2_YvsXcheck = NULL__null; |
64 | static TH2I* h2_YvsX1 = NULL__null; |
65 | static TH2I* h2_E9vsp = NULL__null; |
66 | static TH2I* h2_E1vsp = NULL__null; |
67 | static TH2I* h2_E2vsp = NULL__null; |
68 | static TH2I* h2_E1vspPos = NULL__null; |
69 | static TH2I* h2_dEdxvsE9 = NULL__null; |
70 | static TH2I* h2_intOverPeak_vsEfcal = NULL__null; |
71 | static TH2I* h2_E1_vsintOverPeak = NULL__null; |
72 | static TH2I* h2_E1peak_vsp = NULL__null; |
73 | static TH2I* h2_E1vsRlocal = NULL__null; |
74 | static TH2I* h2_E9_vsTheta = NULL__null; |
75 | static TH2I* h2_E9_vsPhi = NULL__null; |
76 | static TH2I* h2_E1_vsTheta = NULL__null; |
77 | static TH2I* h2_E1_vsPhi = NULL__null; |
78 | static TH2D* h2_E9_vsThetaPhiw = NULL__null; |
79 | static TH2D* h2_E1_vsThetaPhiw = NULL__null; |
80 | static TH2D* h2_E1_vsThetaPhi = NULL__null; |
81 | |
82 | |
83 | |
84 | |
85 | |
86 | |
87 | extern "C" |
88 | { |
89 | void InitPlugin(JApplication *locApplication) |
90 | { |
91 | InitJANAPlugin(locApplication); |
92 | locApplication->AddProcessor(new DEventProcessor_fcal_charged()); |
93 | } |
94 | } |
95 | |
96 | |
97 | |
98 | |
99 | jerror_t DEventProcessor_fcal_charged::init(void) |
100 | { |
101 | |
102 | |
103 | |
104 | |
105 | japp->RootWriteLock(); |
106 | |
107 | |
108 | |
109 | if(h1_deltaX != NULL__null){ |
110 | printf ("TEST>> DEventProcessor_fcal_charged::init - Other threads continue\n"); |
111 | japp->RootUnLock(); |
112 | return NOERROR; |
113 | } |
114 | |
115 | |
116 | |
117 | |
118 | |
119 | |
120 | |
121 | TDirectory *main = gDirectory(TDirectory::CurrentDirectory()); |
122 | gDirectory(TDirectory::CurrentDirectory())->mkdir("fcal_charged")->cd(); |
123 | |
124 | |
125 | |
126 | int const nbins=100; |
127 | |
128 | h1_stats = new TH1I("h1_stats", "Run Statistics", 20,0,20); |
129 | |
130 | h1_deltaX = new TH1I("h1_deltaX", "Fcal X - Track X", nbins,-25,25); |
131 | h1_deltaX->SetXTitle("Fcal - Track X (cm)"); |
132 | h1_deltaX->SetYTitle("counts"); |
133 | h1_deltaY = new TH1I("h1_deltaY", "Fcal Y - Track X", nbins,-25,25); |
134 | h1_deltaY->SetXTitle("Fcal - Track Y (cm)"); |
135 | h1_deltaY->SetYTitle("counts"); |
136 | |
137 | h1_deltaX_tof = new TH1I("h1_deltaX_tof", "TOF X - Track X E1", nbins,-25,25); |
138 | h1_deltaX_tof->SetXTitle("TOF X - Track X (cm)"); |
139 | h1_deltaX_tof->SetYTitle("counts"); |
140 | h1_deltaY_tof = new TH1I("h1_deltaY_tof", "TOF Y - Track Y E1", nbins,-25,25); |
141 | h1_deltaY_tof->SetXTitle("TOF Y - Track Y (cm)"); |
142 | h1_deltaY_tof->SetYTitle("counts"); |
143 | |
144 | h1_Ebcal = new TH1I("h1_Ebcal", "Sum of point energy", nbins,0,1); |
145 | h1_Ebcal->SetXTitle("Bcal point energy (GeV)"); |
146 | h1_Ebcal->SetYTitle("counts"); |
147 | |
148 | h1_Efcal = new TH1I("h1_Efcal", "Hit energy", nbins,0,4); |
149 | h1_Efcal->SetXTitle("Fcal hit energy (GeV)"); |
150 | h1_Efcal->SetYTitle("counts"); |
151 | h1_tfcal = new TH1I("h1_tfcal", "Hit time", 250,-25,225); |
152 | h1_tfcal->SetXTitle("Fcal hit time (ns)"); |
153 | h1_tfcal->SetYTitle("counts"); |
154 | h1_intOverPeak = new TH1I("h1_intOverPeak", "Hit intOverPeak",nbins,0,25); |
155 | h1_intOverPeak->SetXTitle("Fcal hit intOverPeak"); |
156 | h1_intOverPeak->SetYTitle("counts"); |
157 | |
158 | h2_dEdx_vsp = new TH2I("h2_dEdx_vsp", "Track dEdx vs p",nbins,0,4,nbins,0,10); |
159 | h2_dEdx_vsp->SetXTitle("p (GeV/c)"); |
160 | h2_dEdx_vsp->SetYTitle("dEdx (keV/cm)"); |
161 | h2_dEdx9_vsp = new TH2I("h2_dEdx9_vsp", "Track dEdx9 vs p",nbins,0,4,nbins,0,10); |
162 | h2_dEdx9_vsp->SetXTitle("p (GeV/c)"); |
163 | h2_dEdx9_vsp->SetYTitle("dEdx (keV/cm)"); |
164 | |
165 | h1_N9 = new TH1I("h1_N9", "Hit N9",25,0,25); |
166 | h1_N9->SetXTitle("Number of Hits N9"); |
167 | h1_N9->SetYTitle("counts"); |
168 | h1_E9 = new TH1I("h1_E9", "Energy E9",nbins,0,2); |
169 | h1_E9->SetXTitle("Energy E9 (GeV)"); |
170 | h1_E9->SetYTitle("counts"); |
171 | h1_t9 = new TH1I("h1_t9", "Time t9",250,-25,225); |
172 | h1_t9->SetXTitle("Time t9 (ns)"); |
173 | h1_t9->SetYTitle("counts"); |
174 | h1_t9sigma = new TH1I("h1_t9sigma", "Time t9sigma",nbins,0,10); |
175 | h1_t9sigma->SetXTitle("Time spread t9sigma (ns)"); |
176 | h1_t9sigma->SetYTitle("counts"); |
177 | |
178 | h2_YvsX9 = new TH2I("h2_YvsX9", "Hit position Y vs X, E9",240,-120,120,240,-120,120); |
179 | h2_YvsX9->SetXTitle("X (cm)"); |
180 | h2_YvsX9->SetYTitle("Y (cm)"); |
181 | h2_YvsXcheck = new TH2I("h2_YvsXcheck", "Delta Hit Y vs X Checkered",nbins,-5,5,nbins,-5,5); |
182 | h2_YvsXcheck->SetXTitle("X (cm)"); |
183 | h2_YvsXcheck->SetYTitle("Y (cm)"); |
184 | h2_E1vsRlocal = new TH2I("h2_E1vsRlocal", "E1 vs Rtrk rel to Block center (cm)",nbins,0,5,nbins,0,4); |
185 | h2_E1vsRlocal->SetXTitle("Rlocal (cm)"); |
186 | h2_E1vsRlocal->SetYTitle("E1 (GeV)"); |
187 | h2_E9vsp = new TH2I("h2_E9vsp", "E9 vs p",nbins,0,4,nbins,0,4); |
188 | h2_E9vsp->SetXTitle("p (GeV)"); |
189 | h2_E9vsp->SetYTitle("E9 (GeV)"); |
190 | h2_dEdxvsE9 = new TH2I("h2_dEdxvsE9", "dEdx vs E9 energy",nbins,0,4,nbins,0,4); |
191 | h2_dEdxvsE9->SetXTitle("E9 (GeV)"); |
192 | h2_dEdxvsE9->SetYTitle("dEdx (keV/cm)"); |
193 | |
194 | h1_N1 = new TH1I("h1_N1", "Hit N1",25,0,25); |
195 | h1_N1->SetXTitle("Number of Hits N1"); |
196 | h1_N1->SetYTitle("counts"); |
197 | h1_E1 = new TH1I("h1_E1", "Energy E1",nbins,0,2); |
198 | h1_E1->SetXTitle("Energy E1 (GeV)"); |
199 | h1_E1->SetYTitle("counts"); |
200 | h1_t1 = new TH1I("h1_t1", "Time t1",250,-25,225); |
201 | h1_t1->SetXTitle("Time t1 (ns)"); |
202 | h1_t1->SetYTitle("counts"); |
203 | h1_t1sigma = new TH1I("h1_t1sigma", "Time t1sigma",nbins,0,10); |
204 | h1_t1sigma->SetXTitle("Time spread t1sigma (ns)"); |
205 | h1_t1sigma->SetYTitle("counts"); |
206 | |
207 | h2_YvsX1 = new TH2I("h2_YvsX1", "Hit position Y vs X, E1",240,-120,120,240,-120,120); |
208 | h2_YvsX1->SetXTitle("X (cm)"); |
209 | h2_YvsX1->SetYTitle("Y (cm)"); |
210 | h2_E1vsp = new TH2I("h2_E1vsp", "E1 vs p",nbins,0,4,nbins,0,4); |
211 | h2_E1vsp->SetXTitle("p (GeV)"); |
212 | h2_E1vsp->SetYTitle("E1 (GeV)"); |
213 | h2_E2vsp = new TH2I("h2_E2vsp", "E2 vs p",nbins,0,4,nbins,0,4); |
214 | h2_E2vsp->SetXTitle("p (GeV)"); |
215 | h2_E2vsp->SetYTitle("E2 (GeV)"); |
216 | h2_E1vspPos = new TH2I("h2_E1vspPos", "E1 vs p Q>0",nbins,0,4,nbins,0,4); |
217 | h2_E1vspPos->SetXTitle("p (GeV)"); |
218 | h2_E1vspPos->SetYTitle("E1 (GeV)"); |
219 | h2_E1peak_vsp = new TH2I("h2_E1peak_vsp", "E1peak*6 vs p",nbins,0,4,nbins,0,4); |
220 | h2_E1peak_vsp->SetXTitle("p (GeV)"); |
221 | h2_E1peak_vsp->SetYTitle("E1 peak*6(GeV)"); |
222 | h2_intOverPeak_vsEfcal = new TH2I("h2_intOverPeak_vsEfcal", "intOverPeak vs Efcal",nbins,0,1,nbins,0,25); |
223 | h2_intOverPeak_vsEfcal->SetXTitle("Efcal (GeV)"); |
224 | h2_intOverPeak_vsEfcal->SetYTitle("IntOverPeak"); |
225 | h2_E1_vsintOverPeak = new TH2I("h2_E1_vsintOverPeak", "E1 vs intOverPeak",nbins,0,25,nbins,0,4); |
226 | h2_E1_vsintOverPeak->SetXTitle("IntOverPeak"); |
227 | h2_E1_vsintOverPeak->SetYTitle("E1 (GeV)"); |
228 | |
229 | h2_E9_vsTheta = new TH2I("h2_E9_vsTheta", "E9 vs Theta",90,0,30,nbins,0,4); |
230 | h2_E9_vsTheta->SetXTitle("Theta (deg)"); |
231 | h2_E9_vsTheta->SetYTitle("E9 (GeV)"); |
232 | h2_E1_vsTheta = new TH2I("h2_E1_vsTheta", "E1 vs Theta",90,0,30,nbins,0,4); |
233 | h2_E1_vsTheta->SetXTitle("Theta (deg)"); |
234 | h2_E1_vsTheta->SetYTitle("E1 (GeV)"); |
235 | h2_E9_vsPhi = new TH2I("h2_E9_vsPhi", "E9 vs Phi",90,-180,180,nbins,0,4); |
236 | h2_E9_vsPhi->SetXTitle("Phi (deg)"); |
237 | h2_E9_vsPhi->SetYTitle("E9 (GeV)"); |
238 | h2_E1_vsPhi = new TH2I("h2_E1_vsPhi", "E1 vs Phi",90,-180,180,nbins,0,4); |
239 | h2_E1_vsPhi->SetXTitle("Phi (deg)"); |
240 | h2_E1_vsPhi->SetYTitle("E1 (GeV)"); |
241 | |
242 | h2_E9_vsThetaPhiw = new TH2D("h2_E9_vsThetaPhiw", "E9 vs ThetaPhiw",90,-180,180,90,0,30); |
243 | h2_E9_vsThetaPhiw->SetXTitle("Phi (deg)"); |
244 | h2_E9_vsThetaPhiw->SetYTitle("Theta (deg)"); |
245 | h2_E1_vsThetaPhi = new TH2D("h2_E1_vsThetaPhi", "Unity vs ThetaPhi",90,-180,180,90,0,30); |
246 | h2_E1_vsThetaPhi->SetXTitle("Phi (deg)"); |
247 | h2_E1_vsThetaPhi->SetYTitle("Theta (deg)"); |
248 | h2_E1_vsThetaPhiw = new TH2D("h2_E1_vsThetaPhiw", "E1 vs ThetaPhiw",90,-180,180,90,0,30); |
249 | h2_E1_vsThetaPhiw->SetXTitle("Phi (deg)"); |
250 | h2_E1_vsThetaPhiw->SetYTitle("Theta (deg)"); |
251 | |
252 | |
253 | printf ("TEST>> DEventProcessor_fcal_charged::init - First thread created histograms\n"); |
254 | main->cd(); |
255 | |
256 | japp->RootUnLock(); |
257 | |
258 | |
259 | return NOERROR; |
260 | } |
261 | |
262 | |
263 | |
264 | |
265 | |
266 | jerror_t DEventProcessor_fcal_charged::brun(jana::JEventLoop* locEventLoop, int locRunNumber) |
267 | { |
268 | |
269 | |
270 | |
271 | |
272 | |
273 | |
274 | |
275 | |
276 | |
277 | |
278 | |
279 | |
280 | |
281 | |
282 | |
283 | |
284 | |
285 | |
286 | return NOERROR; |
287 | } |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | jerror_t DEventProcessor_fcal_charged::evnt(jana::JEventLoop* locEventLoop, int locEventNumber) |
295 | { |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | |
302 | |
303 | |
304 | |
305 | |
306 | |
307 | |
308 | |
309 | |
310 | |
311 | |
312 | |
313 | |
314 | |
315 | vector<const DFCALShower*> locFCALShowers; |
316 | vector<const DBCALPoint*> bcalpoints; |
317 | vector<const DTOFPoint*> tofpoints; |
318 | vector<const DFCALHit*> fcalhits; |
319 | vector<const DFCALCluster*> locFCALClusters; |
320 | vector<const DVertex*> kinfitVertex; |
321 | |
322 | |
323 | locEventLoop->Get(locFCALShowers); |
324 | locEventLoop->Get(bcalpoints);; |
325 | locEventLoop->Get(tofpoints); |
326 | locEventLoop->Get(fcalhits); |
327 | locEventLoop->Get(locFCALClusters); |
328 | locEventLoop->Get(kinfitVertex); |
329 | |
330 | vector<const DChargedTrack*> locChargedTrack; |
331 | locEventLoop->Get(locChargedTrack); |
332 | |
333 | DVector3 trkpos(0.0,0.0,0.0); |
334 | DVector3 proj_mom(0.0,0.0,0.0); |
335 | DVector3 trkpos_tof(0.0,0.0,0.0); |
336 | DVector3 proj_mom_tof(0.0,0.0,0.0); |
337 | Double_t zfcal=638; |
338 | Double_t ztof=618.8; |
339 | DVector3 fcal_origin(0.0,0.0,zfcal); |
340 | DVector3 fcal_normal(0.0,0.0,1.0); |
341 | DVector3 tof_origin(0.0,0.0,ztof); |
342 | DVector3 tof_normal(0.0,0.0,1.0); |
343 | |
344 | |
345 | |
346 | |
347 | |
348 | |
349 | |
350 | |
351 | if (locEventNumber%100 == 0) printf ("EventNumber=%d\n",locEventNumber); |
352 | |
353 | japp->RootWriteLock(); |
354 | |
355 | double p; |
356 | double dEdx; |
357 | double pmin=1; |
358 | |
359 | |
360 | for (unsigned int i=0; i < locChargedTrack.size() ; ++i){ |
361 | const DChargedTrack *ctrk = locChargedTrack[i]; |
362 | const DChargedTrackHypothesis *bestctrack = ctrk->Get_BestFOM(); |
363 | vector<const DTrackTimeBased*> locTrackTimeBased; |
364 | bestctrack->Get(locTrackTimeBased); |
365 | |
366 | double FOM = TMath::Prob(locTrackTimeBased[0]->chisq, locTrackTimeBased[0]->Ndof); |
367 | |
368 | |
369 | |
370 | if (locTrackTimeBased[0]->rt->GetIntersectionWithPlane(fcal_origin,fcal_normal,trkpos,proj_mom,NULL__null,NULL__null,NULL__null,SYS_FCAL)==NOERROR){ |
371 | |
372 | |
373 | h1_stats->Fill(0); |
374 | if (locTrackTimeBased[0]->rt->GetIntersectionWithPlane(tof_origin,tof_normal,trkpos_tof,proj_mom_tof,NULL__null,NULL__null,NULL__null,SYS_TOF)!=NOERROR) { |
375 | continue; |
376 | } |
377 | h1_stats->Fill(1); |
378 | |
379 | |
380 | |
381 | float bcal_energy = 0; |
382 | for (unsigned int j=0; j < bcalpoints.size(); ++j) { |
383 | bcal_energy += bcalpoints[j]->E(); |
384 | } |
385 | if (bcal_energy < 0.15) { |
386 | continue; |
387 | } |
388 | h1_Ebcal->Fill(bcal_energy); |
389 | h1_stats->Fill(2); |
390 | |
391 | |
392 | double q = locTrackTimeBased[0]->rt->q; |
393 | p = locTrackTimeBased[0]->momentum().Mag(); |
394 | double trkmass = locTrackTimeBased[0]->mass(); |
395 | |
396 | double radius = sqrt(trkpos.X()*trkpos.X() + trkpos.Y()*trkpos.Y()); |
397 | dEdx = locTrackTimeBased[0]->dEdx()*1e6; |
398 | if(!(trkmass < 0.15 && FOM>0.01 && radius>20)) { |
399 | continue; |
400 | } |
401 | h1_stats->Fill(3); |
402 | |
403 | if (! (p>pmin)) { |
404 | continue; |
405 | } |
406 | h1_stats->Fill (4); |
407 | |
408 | h2_dEdx_vsp->Fill(p,dEdx); |
409 | |
410 | |
411 | |
412 | double trkposX = trkpos.X(); |
413 | double trkposY = trkpos.Y(); |
414 | int trkrow = fcalgeom.row((float)trkposY); |
415 | int trkcol = fcalgeom.column((float)trkposX); |
416 | double dX_tof = 10000; |
417 | double dY_tof = 10000; |
418 | double dR_tof = 10000; |
419 | |
420 | |
421 | |
422 | for (unsigned int j=0; j < tofpoints.size(); ++j) { |
423 | double tofx = 10000; |
424 | double tofy = 10000; |
425 | if (tofpoints[j]->Is_XPositionWellDefined() && tofpoints[j]->Is_YPositionWellDefined()) { |
426 | DVector3 pos_tof = tofpoints[j]->pos; |
427 | tofx = pos_tof.X(); |
428 | tofy = pos_tof.Y(); |
429 | printf ("Event=%d tofx=%f, tofy=%f, trkx=%f, trky=%f \n",locEventNumber,tofx,tofy,trkposX,trkposY); |
430 | double dX_tof1 = tofx - trkpos_tof.X(); |
431 | double dY_tof1 = tofy - trkpos_tof.Y(); |
432 | if (sqrt(dX_tof1*dX_tof1 + dY_tof1*dY_tof1) < dR_tof) { |
433 | dX_tof = dX_tof1; |
434 | dY_tof = dY_tof1; |
435 | dR_tof = sqrt(dX_tof*dX_tof + dY_tof*dY_tof); |
436 | } |
437 | } |
438 | } |
439 | |
440 | printf ("\n1 Event=%d dX_tof=%f, DY_tof=%f, trkx=%f, trky=%f \n",locEventNumber,dX_tof,dY_tof,trkposX,trkposY); |
441 | if ( !(abs(dX_tof)<10 && abs(dY_tof)<10) ) continue; |
442 | |
443 | h1_stats->Fill(5); |
444 | printf ("2 Event=%d dX_tof=%f, DY_tof=%f, trkx=%f, trky=%f \n",locEventNumber,dX_tof,dY_tof,trkposX,trkposY); |
445 | 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()); |
446 | |
447 | |
448 | |
449 | double theta = proj_mom.Theta()*180./M_PI3.14159265358979323846; |
450 | double phi = proj_mom.Phi()*180./M_PI3.14159265358979323846; |
451 | printf ("Event=%d p=%f,theta=%f, phi=%f\n",locEventNumber,p,theta,phi); |
452 | |
453 | |
454 | |
455 | |
456 | double E9=0; |
457 | double E9peak=0; |
458 | double x9=0; |
459 | double y9=0; |
460 | double t9=0; |
461 | double t9sq=0; |
462 | double t9sigma=0; |
463 | int N9=0; |
464 | int Delta_block=2; |
465 | int row_E1=-1000; |
466 | int col_E1=-1000; |
467 | int drow_E1=1000; |
468 | int dcol_E1=1000; |
469 | double dX_E1=-1000; |
470 | double dY_E1=-1000; |
471 | |
472 | int row_offset = 0; |
473 | int col_offset = 0; |
474 | trkrow += row_offset; |
475 | trkcol += col_offset; |
476 | |
477 | for (unsigned int j=0; j < fcalhits.size(); ++j) { |
478 | int row = fcalhits[j]->row; |
479 | int col = fcalhits[j]->column; |
480 | double x = fcalhits[j]->x; |
481 | double y = fcalhits[j]->y; |
482 | double Efcal = fcalhits[j]->E; |
483 | double tfcal= fcalhits[j]->t; |
484 | double intOverPeak = fcalhits[j]->intOverPeak; |
485 | |
486 | |
487 | |
488 | int drow = abs(row - trkrow); |
489 | int dcol = abs(col - trkcol); |
490 | |
491 | h1_deltaX->Fill(x - trkposX); |
492 | h1_deltaY->Fill(y - trkposY); |
493 | h1_Efcal->Fill(Efcal); |
494 | h1_tfcal->Fill(tfcal); |
495 | h1_intOverPeak->Fill(intOverPeak); |
496 | h2_intOverPeak_vsEfcal->Fill(Efcal,intOverPeak); |
497 | |
498 | |
499 | if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50) && (intOverPeak>2.5 && intOverPeak<9)) { |
500 | |
501 | E9 += Efcal; |
502 | E9peak += Efcal*6/intOverPeak; |
503 | x9 += x; |
504 | y9 += y; |
505 | t9 += tfcal; |
506 | t9sq += tfcal*tfcal; |
507 | N9 += 1; |
508 | |
509 | |
510 | row_E1 = row; |
511 | col_E1 = col; |
512 | drow_E1 = drow; |
513 | dcol_E1 = dcol; |
| Value stored to 'dcol_E1' is never read |
514 | dX_E1 = x - trkposX; |
515 | dY_E1 = y - trkposY; |
516 | printf ("Event=%d, trkposX=%f, trkposY=%f, dX_E1=%f, dY_E1=%f\n",locEventNumber,trkposX,trkposY,dX_E1,dY_E1); |
517 | } |
518 | |
519 | } |
520 | |
521 | x9 = N9>0? x9/N9 : 0; |
522 | y9 = N9>0? y9/N9 : 0; |
523 | t9 = N9>0? t9/N9 : 0; |
524 | t9sigma = N9>0? sqrt(t9sq/N9 - t9*t9): 0; |
525 | |
526 | |
527 | |
528 | h1_stats->Fill(6); |
529 | if (E9 > 0.8 && theta<4) continue; |
530 | h1_stats->Fill(7); |
531 | if (N9>0) { |
532 | h1_N9->Fill(N9); |
533 | h1_E9->Fill(E9); |
534 | h1_t9->Fill(t9); |
535 | h1_t9sigma->Fill(t9sigma); |
536 | h2_YvsX9->Fill(trkposX,trkposY); |
537 | h2_dEdx9_vsp->Fill(p,dEdx); |
538 | h2_E9vsp->Fill(p,E9); |
539 | h2_dEdxvsE9->Fill(E9,dEdx); |
540 | h1_stats->Fill(8); |
541 | h2_E9_vsThetaPhiw->Fill(phi,theta,E9); |
542 | h2_E9_vsTheta->Fill(theta,E9); |
543 | h2_E9_vsPhi->Fill(phi,E9); |
544 | } |
545 | if (N9==1) { |
546 | h1_N1->Fill(N9); |
547 | h1_E1->Fill(E9); |
548 | h1_t1->Fill(t9); |
549 | h1_t1sigma->Fill(t9sigma); |
550 | h2_YvsX1->Fill(trkposX,trkposY); |
551 | h2_E1vsp->Fill(p,E9); |
552 | if (q > 0) h2_E1vspPos->Fill(p,E9); |
553 | h2_E1peak_vsp->Fill(p,E9peak); |
554 | h2_E1_vsintOverPeak->Fill(E9*6/E9peak,E9); |
555 | h1_deltaX_tof->Fill(dX_tof); |
556 | h1_deltaY_tof->Fill(dY_tof); |
557 | h1_stats->Fill(9); |
558 | h2_E1_vsThetaPhi->Fill(phi,theta); |
559 | h2_E1_vsThetaPhiw->Fill(phi,theta,E9); |
560 | h2_E1_vsTheta->Fill(theta,E9); |
561 | h2_E1_vsPhi->Fill(phi,E9); |
562 | double Rlocal = sqrt(dX_E1*dX_E1 + dY_E1*dY_E1); |
563 | h2_E1vsRlocal->Fill(Rlocal,E9); |
564 | |
565 | h2_YvsXcheck->Fill(dX_E1,dY_E1); |
566 | } |
567 | if (N9==1 || N9==2) { |
568 | h2_E2vsp->Fill(p,E9); |
569 | h1_stats->Fill(10); |
570 | } |
571 | |
572 | |
573 | |
574 | } |
575 | } |
576 | |
577 | |
578 | |
579 | japp->RootUnLock(); |
580 | |
581 | |
582 | |
583 | |
584 | |
585 | |
586 | |
587 | return NOERROR; |
588 | } |
589 | |
590 | |
591 | |
592 | |
593 | jerror_t DEventProcessor_fcal_charged::erun(void) |
594 | { |
595 | |
596 | |
597 | |
598 | |
599 | japp->RootWriteLock(); |
600 | |
601 | |
602 | japp->RootUnLock(); |
603 | |
604 | |
605 | return NOERROR; |
606 | } |
607 | |
608 | |
609 | |
610 | |
611 | jerror_t DEventProcessor_fcal_charged::fini(void) |
612 | { |
613 | |
614 | |
615 | japp->RootWriteLock(); |
616 | |
617 | |
618 | japp->RootUnLock(); |
619 | return NOERROR; |
620 | } |
621 | |