Bug Summary

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