Bug Summary

File:plugins/Analysis/fcal_charged/DEventProcessor_fcal_charged.cc
Location:line 519, column 6
Description:Value stored to 'y9' 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
38static TH1I* h1_stats = NULL__null;
39static TH1I* h1_deltaX = NULL__null;
40static TH1I* h1_deltaY = NULL__null;
41static TH1I* h1_deltaX_tof = NULL__null;
42static TH1I* h1_deltaY_tof = NULL__null;
43static TH1I* h1_Efcal = NULL__null;
44static TH1I* h1_Ebcal = NULL__null;
45static TH1I* h1_tfcal = NULL__null;
46static TH1I* h1_intOverPeak = NULL__null;
47
48static TH1I* h1_N9 = NULL__null;
49static TH1I* h1_E9 = NULL__null;
50static TH1I* h1_t9 = NULL__null;
51static TH1I* h1_t9sigma = NULL__null;
52static TH1I* h1_N1 = NULL__null;
53static TH1I* h1_E1 = NULL__null;
54static TH1I* h1_t1 = NULL__null;
55static TH1I* h1_t1sigma = NULL__null;
56
57static TH2I* h2_dEdx_vsp = NULL__null;
58static TH2I* h2_dEdx9_vsp = NULL__null;
59static TH2I* h2_YvsX9 = NULL__null;
60static TH2I* h2_YvsXcheck = NULL__null;
61static TH2I* h2_YvsX1 = NULL__null;
62static TH2I* h2_E9vsp = NULL__null;
63static TH2I* h2_E1vsp = NULL__null;
64static TH2I* h2_E2vsp = NULL__null;
65static TH2I* h2_E1vspPos = NULL__null;
66static TH2I* h2_dEdxvsE9 = NULL__null;
67static TH2I* h2_intOverPeak_vsEfcal = NULL__null;
68static TH2I* h2_E1_vsintOverPeak = NULL__null;
69static TH2I* h2_E1peak_vsp = NULL__null;
70static TH2I* h2_E1vsRlocal = NULL__null;
71static TH2I* h2_E9_vsTheta = NULL__null;
72static TH2I* h2_E9_vsPhi = NULL__null;
73static TH2I* h2_E1_vsTheta = NULL__null;
74static TH2I* h2_E1_vsPhi = NULL__null;
75static TH2D* h2_E9_vsThetaPhiw = NULL__null;
76static TH2D* h2_E1_vsThetaPhiw = NULL__null;
77static TH2D* h2_E1_vsThetaPhi = NULL__null;
78
79
80
81
82// Routine used to create our DEventProcessor
83
84extern "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//------------------
96jerror_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(TDirectory::CurrentDirectory());
112 gDirectory(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//------------------
253jerror_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
281jerror_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;
519 y9 = N9>0? y9/N9 : 0;
Value stored to 'y9' is never read
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//------------------
588jerror_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//------------------
601jerror_t DEventProcessor_fcal_charged::fini(void)
602{
603 // Called before program exit after event processing is finished.
604
605 return NOERROR;
606}
607