Bug Summary

File:plugins/Analysis/fcal_charged/DEventProcessor_fcal_charged.cc
Location:line 511, column 5
Description:Value stored to 'col_E1' is never read

Annotated Source Code

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
38DFCALGeometry fcalgeom;
39
40
41static TH1I* h1_stats = NULL__null;
42static TH1I* h1_deltaX = NULL__null;
43static TH1I* h1_deltaY = NULL__null;
44static TH1I* h1_deltaX_tof = NULL__null;
45static TH1I* h1_deltaY_tof = NULL__null;
46static TH1I* h1_Efcal = NULL__null;
47static TH1I* h1_Ebcal = NULL__null;
48static TH1I* h1_tfcal = NULL__null;
49static TH1I* h1_intOverPeak = NULL__null;
50
51static TH1I* h1_N9 = NULL__null;
52static TH1I* h1_E9 = NULL__null;
53static TH1I* h1_t9 = NULL__null;
54static TH1I* h1_t9sigma = NULL__null;
55static TH1I* h1_N1 = NULL__null;
56static TH1I* h1_E1 = NULL__null;
57static TH1I* h1_t1 = NULL__null;
58static TH1I* h1_t1sigma = NULL__null;
59
60static TH2I* h2_dEdx_vsp = NULL__null;
61static TH2I* h2_dEdx9_vsp = NULL__null;
62static TH2I* h2_YvsX9 = NULL__null;
63static TH2I* h2_YvsXcheck = NULL__null;
64static TH2I* h2_YvsX1 = NULL__null;
65static TH2I* h2_E9vsp = NULL__null;
66static TH2I* h2_E1vsp = NULL__null;
67static TH2I* h2_E2vsp = NULL__null;
68static TH2I* h2_E1vspPos = NULL__null;
69static TH2I* h2_dEdxvsE9 = NULL__null;
70static TH2I* h2_intOverPeak_vsEfcal = NULL__null;
71static TH2I* h2_E1_vsintOverPeak = NULL__null;
72static TH2I* h2_E1peak_vsp = NULL__null;
73static TH2I* h2_E1vsRlocal = NULL__null;
74static TH2I* h2_E9_vsTheta = NULL__null;
75static TH2I* h2_E9_vsPhi = NULL__null;
76static TH2I* h2_E1_vsTheta = NULL__null;
77static TH2I* h2_E1_vsPhi = NULL__null;
78static TH2D* h2_E9_vsThetaPhiw = NULL__null;
79static TH2D* h2_E1_vsThetaPhiw = NULL__null;
80static TH2D* h2_E1_vsThetaPhi = NULL__null;
81
82
83
84
85// Routine used to create our DEventProcessor
86
87extern "C"
88{
89 void InitPlugin(JApplication *locApplication)
90 {
91 InitJANAPlugin(locApplication);
92 locApplication->AddProcessor(new DEventProcessor_fcal_charged()); //register this plugin
93 }
94} // "C"
95
96//------------------
97// init
98//------------------
99jerror_t DEventProcessor_fcal_charged::init(void)
100{
101 // This is called once at program startup. If you are creating
102 // and filling historgrams in this plugin, you should lock the
103 // ROOT mutex like this:
104
105 japp->RootWriteLock();
106
107 // First thread to get here makes all histograms. If one pointer is
108 // already not NULL, assume all histograms are defined and return now
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 // ... create historgrams or trees ...
117
118 // TDirectory *dir = new TDirectoryFile("FCAL","FCAL");
119 // dir->cd();
120 // create root folder for bcal and cd to it, store main dir
121 TDirectory *main = gDirectory(TDirectory::CurrentDirectory());
122 gDirectory(TDirectory::CurrentDirectory())->mkdir("fcal_charged")->cd();
123
124
125 // double pi0_mass = 0.1349766;
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 // back to main dir
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// brun
265//------------------
266jerror_t DEventProcessor_fcal_charged::brun(jana::JEventLoop* locEventLoop, int locRunNumber)
267{
268 // This is called whenever the run number changes
269 //BCAL_Neutrals->Fill();
270 //cout << " run number = " << RunNumber << endl;
271 /*
272 //Optional: Retrieve REST writer for writing out skims
273 locEventLoop->GetSingle(dEventWriterREST);
274 */
275
276 //vector<const DTrackFinder *> finders;
277 //locEventLoop->Get(finders);
278 //finder = const_cast<DTrackFinder*>(finders[0]);
279
280 /*
281 //Recommeded: Create output ROOT TTrees (nothing is done if already created)
282 locEventLoop->GetSingle(dEventWriterROOT);
283 dEventWriterROOT->Create_DataTrees(locEventLoop);
284 */
285
286 return NOERROR;
287}
288
289//------------------
290// evnt
291//------------------
292
293
294jerror_t DEventProcessor_fcal_charged::evnt(jana::JEventLoop* locEventLoop, int locEventNumber)
295{
296
297 // This is called for every event. Use of common resources like writing
298 // to a file or filling a histogram should be mutex protected. Using
299 // locEventLoop->Get(...) to get reconstructed objects (and thereby activating the
300 // reconstruction algorithm) should be done outside of any mutex lock
301 // since multiple threads may call this method at the same time.
302 //
303 // Here's an example:
304 //
305 // vector<const MyDataClass*> mydataclasses;
306 // locEventLoop->Get(mydataclasses);
307 //
308 // japp->RootWriteLock();
309 // ... fill historgrams or trees ...
310 // japp->RootUnLock();
311
312 // DOCUMENTATION:
313 // ANALYSIS library: https://halldweb1.jlab.org/wiki/index.php/GlueX_Analysis_Software
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 //const DDetectorMatches* locDetectorMatches = NULL;
322 //locEventLoop->GetSingle(locDetectorMatches);
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 /* Double_t Lfcal=45.;
345 Double_t zfcal_back=zfcal+Lfcal;
346 DVector3 trkpos_fcal_back(0.0,0.0,0.0);;
347 DVector3 proj_mom_back(0.0,0.0,0.0);
348 DVector3 fcal_origin_back(0.0,0.0,zfcal_back);
349 DVector3 fcal_normal_back(0.0,0.0,1.0);*/
350
351 if (locEventNumber%100 == 0) printf ("EventNumber=%d\n",locEventNumber);
352
353 japp->RootWriteLock(); // lock before loop so that histograms can be filled.
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); // locTrackTimeBased[0] contains the best FOM track
365
366 double FOM = TMath::Prob(locTrackTimeBased[0]->chisq, locTrackTimeBased[0]->Ndof);
367
368 // get intersection with fcal front face and backface;
369
370 if (locTrackTimeBased[0]->rt->GetIntersectionWithPlane(fcal_origin,fcal_normal,trkpos,proj_mom,NULL__null,NULL__null,NULL__null,SYS_FCAL)==NOERROR){
371 // if (locTrackTimeBased[0]->rt->GetIntersectionWithPlane(fcal_origin_back,fcal_normal_back,trkpos_fcal_back,proj_mom_back,NULL,NULL,NULL)!=NOERROR)
372 // continue; // get track projection to back of fcal
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; // get track projection to tof plane
376 }
377 h1_stats->Fill(1);
378
379 // sum up all energy and Bcal and keep only if likely BCAL trigger
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; // require that bcal be sufficient for trigger
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()); // distance of track from beamline
397 dEdx = locTrackTimeBased[0]->dEdx()*1e6; // convert to keV
398 if(!(trkmass < 0.15 && FOM>0.01 && radius>20)) {
399 continue; // select tracks of interest
400 }
401 h1_stats->Fill(3);
402
403 if (! (p>pmin)) {
404 continue; // require minimum momentum to pion analyzis
405 }
406 h1_stats->Fill (4);
407
408 h2_dEdx_vsp->Fill(p,dEdx);
409 // use position at the back of fcal // asume zero field out here
410 // double xfcal_back = trkpos.X() + proj_mom.X()*Lfcal/proj_mom.Z();
411 // double yfcal_back = trkpos.Y() + proj_mom.Y()*Lfcal/proj_mom.Z();
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 // get tof matches
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 // get components of p at FCAL
448 // double theta = proj_mom.Theta()*TVector3::RAD2DEG;
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 // loop over fcal hits
455
456 double E9=0; // energy, x, y of 9 blocks surrounding track
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; // =1 for 3x3, =2 for 5x5
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; // offset actual row to look for randoms
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 // 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);
486
487 // fill histograms
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 // select hits
499 if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50) && (intOverPeak>2.5 && intOverPeak<9)) {
500 // if (drow<=Delta_block && dcol<=Delta_block && (tfcal>-15 && tfcal<50)) {
501 E9 += Efcal;
502 E9peak += Efcal*6/intOverPeak; // factor of 6 so that E9peak ~ E9
503 x9 += x;
504 y9 += y;
505 t9 += tfcal;
506 t9sq += tfcal*tfcal;
507 N9 += 1;
508
509 //save for later
510 row_E1 = row;
511 col_E1 = col;
Value stored to 'col_E1' is never read
512 drow_E1 = drow;
513 dcol_E1 = dcol;
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 } // end loop over fcal hits
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 // printf ("\nEvent=%d N9=%d E9=%f x9=%f y9=%f t9=%f t9sigma=%f\n",locEventNumber,N9,E9,x9,y9,t9,t9sigma);
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); // note that this only works because a single cell fires
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 // if (((row_E1%2==0 && col_E1%2==0) || (row_E1%2==1 && col_E1%2==1)) ) h2_YvsXcheck->Fill(dX_E1,dY_E1);
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 } // end track intersection if statement
575 }
576
577
578 //UnlockState();
579 japp->RootUnLock();
580
581
582 /*
583 //Optional: Save event to output REST file. Use this to create skims.
584 dEventWriterREST->Write_RESTEvent(locEventLoop, "FCAL_Shower"); //string is part of output file name
585 */
586
587 return NOERROR;
588}
589
590//------------------
591// erun
592//------------------
593jerror_t DEventProcessor_fcal_charged::erun(void)
594{
595 // This is called whenever the run number changes, before it is
596 // changed to give you a chance to clean up before processing
597 // events from the next run number.
598
599 japp->RootWriteLock();
600
601
602 japp->RootUnLock();
603
604
605 return NOERROR;
606}
607
608//------------------
609// fini
610//------------------
611jerror_t DEventProcessor_fcal_charged::fini(void)
612{
613 // Called before program exit after event processing is finished.
614
615 japp->RootWriteLock();
616
617
618 japp->RootUnLock();
619 return NOERROR;
620}
621