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 | 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 | |
83 | |
84 | extern "C" |
85 | { |
86 | void InitPlugin(JApplication *locApplication) |
87 | { |
88 | InitJANAPlugin(locApplication); |
89 | locApplication->AddProcessor(new DEventProcessor_fcal_charged()); |
90 | } |
91 | } |
92 | |
93 | |
94 | |
95 | |
96 | jerror_t DEventProcessor_fcal_charged::init(void) |
97 | { |
98 | |
99 | |
100 | if(h1_deltaX != NULL__null){ |
101 | printf ("TEST>> DEventProcessor_fcal_charged::init - Other threads continue\n"); |
102 | return NOERROR; |
103 | } |
104 | |
105 | |
106 | |
107 | |
108 | |
109 | |
110 | |
111 | TDirectory *main = gDirectory(TDirectory::CurrentDirectory()); |
112 | gDirectory(TDirectory::CurrentDirectory())->mkdir("fcal_charged")->cd(); |
113 | |
114 | |
115 | |
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 | |
243 | printf ("TEST>> DEventProcessor_fcal_charged::init - First thread created histograms\n"); |
244 | main->cd(); |
245 | |
246 | return NOERROR; |
247 | } |
248 | |
249 | |
250 | |
251 | |
252 | |
253 | jerror_t DEventProcessor_fcal_charged::brun(jana::JEventLoop* locEventLoop, int locRunNumber) |
254 | { |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | |
264 | |
265 | |
266 | |
267 | |
268 | |
269 | |
270 | |
271 | |
272 | |
273 | return NOERROR; |
274 | } |
275 | |
276 | |
277 | |
278 | |
279 | |
280 | |
281 | jerror_t DEventProcessor_fcal_charged::evnt(jana::JEventLoop* locEventLoop, int locEventNumber) |
282 | { |
283 | |
284 | |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
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 | |
309 | |
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 | |
339 | |
340 | |
341 | |
342 | |
343 | |
344 | |
345 | if (locEventNumber%100 == 0) printf ("EventNumber=%d\n",locEventNumber); |
346 | |
347 | |
348 | |
349 | japp->RootFillLock(this); |
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); |
361 | |
362 | double FOM = TMath::Prob(locTrackTimeBased[0]->chisq, locTrackTimeBased[0]->Ndof); |
363 | |
364 | |
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 | |
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; |
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()); |
393 | dEdx = (locTrackTimeBased[0]->dNumHitsUsedFordEdx_CDC >= locTrackTimeBased[0]->dNumHitsUsedFordEdx_FDC) ? locTrackTimeBased[0]->ddEdx_CDC_amp : locTrackTimeBased[0]->ddEdx_FDC; |
394 | dEdx *= 1e6; |
395 | if(!(trkmass < 0.15 && FOM>0.01 && radius>20)) { |
396 | continue; |
397 | } |
398 | h1_stats->Fill(3); |
399 | |
400 | if (! (p>pmin)) { |
401 | continue; |
402 | } |
403 | h1_stats->Fill (4); |
404 | |
405 | h2_dEdx_vsp->Fill(p,dEdx); |
406 | |
407 | |
408 | |
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 | |
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 | |
445 | |
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 | |
452 | |
453 | double E9=0; |
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; |
462 | |
463 | |
464 | |
465 | |
466 | double dX_E1=-1000; |
467 | double dY_E1=-1000; |
468 | |
469 | int row_offset = 0; |
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 | |
483 | |
484 | |
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 | |
496 | if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50) && (intOverPeak>2.5 && intOverPeak<9)) { |
497 | |
498 | E9 += Efcal; |
499 | E9peak += Efcal*6/intOverPeak; |
500 | x9 += x; |
501 | y9 += y; |
502 | t9 += tfcal; |
503 | t9sq += tfcal*tfcal; |
504 | N9 += 1; |
505 | |
506 | |
507 | |
508 | |
509 | |
510 | |
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 | } |
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 | |
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); |
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 | |
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 | } |
572 | } |
573 | |
574 | japp->RootFillUnLock(this); |
575 | |
576 | |
577 | |
578 | |
579 | |
580 | |
581 | |
582 | return NOERROR; |
583 | } |
584 | |
585 | |
586 | |
587 | |
588 | jerror_t DEventProcessor_fcal_charged::erun(void) |
589 | { |
590 | |
591 | |
592 | |
593 | |
594 | |
595 | return NOERROR; |
596 | } |
597 | |
598 | |
599 | |
600 | |
601 | jerror_t DEventProcessor_fcal_charged::fini(void) |
602 | { |
603 | |
604 | |
605 | return NOERROR; |
606 | } |
607 | |