Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_BCAL_point_calib.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_BCAL_point_calib.cc
4 // Created: Mon Sep 26 09:38:23 EDT 2016
5 // Creator: gvasil (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64)
6 //
7 
9 
10 #include <TLorentzVector.h>
11 #include "TMath.h"
12 #include "TF1.h"
13 #include "TCanvas.h"
14 #include "TLatex.h"
15 
16 #include "DANA/DApplication.h"
17 #include "BCAL/DBCALShower.h"
18 #include "BCAL/DBCALTruthShower.h"
19 #include "BCAL/DBCALPoint.h"
20 #include "BCAL/DBCALHit.h"
21 #include "FCAL/DFCALShower.h"
22 #include "TRACKING/DMCThrown.h"
24 #include "BCAL/DBCALGeometry.h"
25 #include "TRIGGER/DTrigger.h"
26 
27 #include <vector>
28 #include "TGraphErrors.h"
29 #include <deque>
30 #include <string>
31 #include <iostream>
32 #include <algorithm>
33 #include <stdio.h>
34 #include <stdlib.h>
35 
36 using namespace jana;
37 using namespace std;
38 
39  // histograms for storing matched shower info
40  static TH1I* h1_Num_matched_showers = NULL;
41  static TH1I* h1_E = NULL;
42  static TH1I* h1_E_raw = NULL;
43  static TH1I* h1_x = NULL;
44  static TH1I* h1_y = NULL;
45  static TH1I* h1_z = NULL;
46  static TH1I* h1_t = NULL;
47  static TH1I* h1_N_cell = NULL;
48 
49  // histograms for storing matched track info
50  static TH1I* h1trk_Num_matched_tracks = NULL;
51  static TH1I* h1trk_FOM = NULL;
52  static TH1I* h1trk_pmag = NULL;
53  static TH1I* h1trk_px = NULL;
54  static TH1I* h1trk_py = NULL;
55  static TH1I* h1trk_pz = NULL;
56  static TH1I* h1trk_x = NULL;
57  static TH1I* h1trk_y = NULL;
58  static TH1I* h1trk_z = NULL;
59  static TH1I* h1trk_energy = NULL;
60 
61  static TH2I* h2_Evsenergy = NULL;
62  static TH2I* h2_pmagvsenergy = NULL;
63 
64  // histograms for storing point info
65  static TH1I* h1pt_Num_points = NULL;
66  static TH1I* h1pt_module = NULL;
67  static TH1I* h1pt_layer = NULL;
68  static TH1I* h1pt_sector = NULL;
69  static TH1I* h1pt_cell_id = NULL;
70  static TH1I* h1pt_energy = NULL;
71  static TH1I* h1pt_energy_US = NULL;
72  static TH1I* h1pt_energy_DS = NULL;
73  static TH1I* h1pt_z = NULL;
74  static TH1I* h1pt_r = NULL;
75  static TH1I* h1pt_r_size = NULL;
76  static TH1I* h1pt_phi = NULL;
77  static TH1I* h1pt_t = NULL;
78 
79  // histograms to check z_track and z_point correlation
80  static TH1I* h1_ztrack_minus_zpoint = NULL;
81  static TH1I* h1_abs_ztrack_minus_zpoint = NULL;
82  static TH2I* h2_ztrack_minus_zpoint_vs_ztrack = NULL;
83 
84  static TH2I* h2_zpoint_vs_ztrack = NULL;
85  static TH2I* h2_zpoint_vs_ztrack_thrown = NULL;
86 
87  // create a graph for each channel (cell)
88  static TGraph* h2_tgraph[768] = {NULL};
89  // graphs of thrown points for each channel (cell)
90  static TGraph* h2_tgraph_thrown[768] = {NULL};
91  // create a graph for each channel for the Delta t calibration (cell)
92  static TGraph* h2_tgraph_deltat[768] = {NULL};
93 
94  // histogram for the percentage of thrown points per channel
95  static TH1D* h1_thrown_per_channel = NULL;
96 
97  int event_count=0;
98 
99  // vectors for storing the z_track and z_point for each channel
100  vector< vector< double > > z_track_graph(768, vector< double >() );
101  vector< vector< double > > z_point_graph(768, vector< double >() );
102  vector< vector< double > > z_deltat_graph(768, vector< double >() );
103  // vectors for storing the thrown points for each channel (if cuts are applied)
104  vector< vector< double > > z_track_thrown(768, vector< double >() );
105  vector< vector< double > > z_point_thrown(768, vector< double >() );
106 
107 // Routine used to create our JEventProcessor
108 #include <JANA/JApplication.h>
109 #include <JANA/JFactory.h>
110 extern "C"{
111 void InitPlugin(JApplication *app){
112  InitJANAPlugin(app);
113  app->AddProcessor(new JEventProcessor_BCAL_point_calib());
114 }
115 } // "C"
116 
117 
118 //------------------
119 // JEventProcessor_BCAL_point_calib (Constructor)
120 //------------------
122 {
123 
124 }
125 
126 //------------------
127 // ~JEventProcessor_BCAL_point_calib (Destructor)
128 //------------------
130 {
131 
132 }
133 
134 //------------------
135 // init
136 //------------------
138 {
139  // This is called once at program startup.
140  // First thread to get here makes all histograms. If one pointer is
141  // already not NULL, assume all histograms are defined and return now
142 
143  // japp->RootFillLock(this);
144  DEBUG=false;
145  VERBOSE=false;
146  gPARMS->SetDefaultParameter("BCALPOINTCALIB:DEBUG", DEBUG, "Control the creation of extra histograms");
147  gPARMS->SetDefaultParameter("BCALPOINTCALIB:DEBUG", VERBOSE, "Control verbose output");
148 
149  if(h1_Num_matched_showers != NULL){
150  return NOERROR;
151  }
152 
153  // lock all root operations
154  japp->RootWriteLock();
155 
156  // TDirectory *dir = new TDirectoryFile("BCAL","BCAL");
157  // dir->cd();
158  // create root folder for bcal and cd to it, store main dir
159  TDirectory *main = gDirectory;
160  gDirectory->mkdir("bcal_point_calibs")->cd();
161 
162 
163  int nbins = 100;
164  int nbins2 = 500;
165 
166 
167  h1_Num_matched_showers = new TH1I("h1_Num_matched_showers", "Number of matched showers per event", 20,0,20);
168  h1_Num_matched_showers->SetXTitle("Number of showers");
169  h1_Num_matched_showers->SetYTitle("counts");
170 
171  h1trk_Num_matched_tracks = new TH1I("h1trk_Num_matched_tracks", "Number of matched tracks per event", 20,0,20);
172  h1trk_Num_matched_tracks->SetXTitle("Number of tracks");
173  h1trk_Num_matched_tracks->SetYTitle("counts");
174 
175  // thrown points per channel (if cuts are applied)
176  h1_thrown_per_channel = new TH1D("h1_thrown_per_channel","Percentage of thrown points per Channel",800,0,800);
177  h1_thrown_per_channel->SetXTitle("Channel #");
178  h1_thrown_per_channel->SetYTitle("Percentage (%)");
179  h1_thrown_per_channel->SetMarkerStyle(20);
180 
181  h1pt_Num_points = new TH1I("h1pt_Num_points", "pt points per shower",20,0,20);
182  h1pt_Num_points->SetXTitle("Points per Shower");
183  h1pt_Num_points->SetYTitle("counts");
184 
185  if(DEBUG) {
186 
187  // shower info
188  h1_E = new TH1I("h1_E", "Matched energy per shower",nbins,0,2);
189  h1_E->SetXTitle("Energy per shower (GeV)");
190  h1_E->SetYTitle("counts");
191  h1_E_raw = new TH1I("h1_E_raw", "Matched raw energy per shower",nbins,0,2);
192  h1_E_raw->SetXTitle("Raw Energy per shower (GeV)");
193  h1_E_raw->SetYTitle("counts");
194  h1_x = new TH1I("h1_x", "x pos of shower",nbins2,-100,100);
195  h1_x->SetXTitle("x (cm)");
196  h1_x->SetYTitle("counts");
197  h1_y = new TH1I("h1_y", "y pos of shower",nbins2,-100,100);
198  h1_y->SetXTitle("y (cm)");
199  h1_y->SetYTitle("counts");
200  h1_z = new TH1I("h1_z", "z pos of shower",nbins2,0,500);
201  h1_z->SetXTitle("z (cm)");
202  h1_z->SetYTitle("counts");
203  h1_t = new TH1I("h1_t", "time of shower",nbins2,-100,100);
204  h1_t->SetXTitle("time of shower (ns)");
205  h1_t->SetYTitle("counts");
206  h1_N_cell = new TH1I("h1_N_cell", "N_cell of shower",20,0,20);
207  h1_N_cell->SetXTitle("N_cell");
208  h1_N_cell->SetYTitle("counts");
209 
210  // track info
211  h1trk_FOM = new TH1I("h1trk_FOM", "FOM for matched tracks", nbins,0,1);
212  h1trk_FOM->SetXTitle("FOM");
213  h1trk_FOM->SetYTitle("counts");
214  h1trk_pmag = new TH1I("h1trk_mag", "pmag for matched tracks", nbins,0,4);
215  h1trk_pmag->SetXTitle("pmag (GeV)");
216  h1trk_pmag->SetYTitle("counts");
217  h1trk_px = new TH1I("h1trk_px", "px for matched tracks", nbins,0,2);
218  h1trk_px->SetXTitle("px (GeV)");
219  h1trk_px->SetYTitle("counts");
220  h1trk_py = new TH1I("h1trk_py", "py for matched tracks", nbins,0,2);
221  h1trk_py->SetXTitle("Number of tracks");
222  h1trk_py->SetYTitle("counts");
223  h1trk_pz = new TH1I("h1trk_pz", "pz for matched tracks", nbins,0,4);
224  h1trk_pz->SetXTitle("pz (GeV)");
225  h1trk_pz->SetYTitle("counts");
226  h1trk_x = new TH1I("h1trk_x", "x for matched tracks", nbins2,-30,30);
227  h1trk_x->SetXTitle("x (cm)");
228  h1trk_x->SetYTitle("counts");
229  h1trk_y = new TH1I("h1trk_y", "y for matched tracks", nbins2,-30,30);
230  h1trk_y->SetXTitle("y (cm)");
231  h1trk_y->SetYTitle("counts");
232  h1trk_z = new TH1I("h1trk_z", "z for matched tracks", nbins2,-50,250);
233  h1trk_z->SetXTitle("z (cm)");
234  h1trk_z->SetYTitle("counts");
235  h1trk_energy = new TH1I("h1trk_energy", "energy for matched tracks", nbins,0,4);
236  h1trk_energy->SetXTitle("energy (GeV)");
237  h1trk_energy->SetYTitle("counts");
238 
239  h2_Evsenergy = new TH2I("h2_Evsenergy", "E vs energy matched tracks", nbins,0,4,nbins,0,4);
240  h2_Evsenergy->SetXTitle("track energy (GeV)");
241  h2_Evsenergy->SetYTitle("Eshower (GeV)");
242  h2_pmagvsenergy = new TH2I("h2_pmagvsenergy", "pmag vs energy matched tracks", nbins,0,4,nbins,0,4);
243  h2_pmagvsenergy->SetXTitle("track energy (GeV)");
244  h2_pmagvsenergy->SetYTitle("track pmag (GeV)");
245 
246  // point info
247  h1pt_module = new TH1I("h1pt_module", "pt module number",50,0,50);
248  h1pt_module->SetXTitle("Module");
249  h1pt_module->SetYTitle("counts");
250  h1pt_layer = new TH1I("h1pt_layer", "pt layer number",5,0,5);
251  h1pt_layer->SetXTitle("Layer");
252  h1pt_layer->SetYTitle("counts");
253  h1pt_sector = new TH1I("h1pt_sector", "pt sector number",5,0,5);
254  h1pt_sector->SetXTitle("Sector");
255  h1pt_sector->SetYTitle("counts");
256  h1pt_cell_id = new TH1I("h1pt_cell_id", "pt cell id number",800,0,800);
257  h1pt_cell_id->SetXTitle("Cell ID");
258  h1pt_cell_id->SetYTitle("counts");
259  h1pt_energy = new TH1I("h1pt_energy", "pt energy",nbins,0,1);
260  h1pt_energy->SetXTitle("Point Energy");
261  h1pt_energy->SetYTitle("counts");
262  h1pt_energy_US = new TH1I("h1pt_energy_US", "pt energy US",nbins,0,1);
263  h1pt_energy_US->SetXTitle("Point Energy US");
264  h1pt_energy_US->SetYTitle("counts");
265  h1pt_energy_DS = new TH1I("h1pt_energy_DS", "pt energy DS",nbins,0,1);
266  h1pt_energy_DS->SetXTitle("Point Energy DS");
267  h1pt_energy_DS->SetYTitle("counts");
268  h1pt_z = new TH1I("h1pt_z", "z for point",nbins2,0,500);
269  h1pt_z->SetXTitle("z (cm)");
270  h1pt_z->SetYTitle("counts");
271  h1pt_r = new TH1I("h1pt_r", "r for point",nbins,60,90);
272  h1pt_r->SetXTitle("r (cm)");
273  h1pt_r->SetYTitle("counts");
274  h1pt_r_size = new TH1I("h1pt_r_size", "r size of point", 20, 0, 20);
275  h1pt_r_size->SetXTitle("r size (cm)");
276  h1pt_r_size->SetYTitle("counts");
277  h1pt_phi = new TH1I("h1pt_phi", "phi of point", 700, 0, 7);
278  h1pt_phi->SetXTitle("phi (rad)");
279  h1pt_phi->SetYTitle("counts");
280  h1pt_t = new TH1I("h1pt_t", "t for point",nbins2,-100,100);
281  h1pt_t->SetXTitle("t");
282  h1pt_t->SetYTitle("counts");
283 
284  // z_track and z_point correlation
285  h1_ztrack_minus_zpoint = new TH1I("h1_ztrack_minus_zpoint", "z_track - z_point", 400, -200, 200);
286  h1_ztrack_minus_zpoint->SetXTitle("z_track - z_point (cm)");
287  h1_ztrack_minus_zpoint->SetYTitle("counts");
288 
289  h1_abs_ztrack_minus_zpoint = new TH1I("h1_abs_ztrack_minus_zpoint", "|z_track - z_point|", 400, 0, 400);
290  h1_abs_ztrack_minus_zpoint->SetXTitle("|z_track - z_point| (cm)");
291  h1_abs_ztrack_minus_zpoint->SetYTitle("counts");
292  h2_ztrack_minus_zpoint_vs_ztrack = new TH2I("h2_ztrack_minus_zpoint_vs_ztrack", "z_track - z_point, vs z_track", 500, -100, 400, 400, -200, 200);
293  h2_ztrack_minus_zpoint_vs_ztrack->SetXTitle("z_track (cm)");
294  h2_ztrack_minus_zpoint_vs_ztrack->SetYTitle("z_track - z_point (cm)");
295  h2_zpoint_vs_ztrack = new TH2I("h2_zpoint_vs_ztrack", "z point vs z track for the BCAL",1000,0,500,1000,0,500);
296  h2_zpoint_vs_ztrack->SetXTitle("z of track (cm)");
297  h2_zpoint_vs_ztrack->SetYTitle("z of point (cm)");
298  h2_zpoint_vs_ztrack_thrown = new TH2I("h2_zpoint_vs_ztrack_thrown", "thrown z point vs z track for the BCAL",1000,0,500,1000,0,500);
299  h2_zpoint_vs_ztrack_thrown->SetXTitle("z of track (cm)");
300  h2_zpoint_vs_ztrack_thrown->SetYTitle("z of point (cm)");
301 
302  }
303 
304  // back to main dir
305  main->cd();
306 
307  // japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
308 
309  // unlock
310  japp->RootUnLock();
311 
312  return NOERROR;
313 
314 }
315 
316 //------------------
317 // brun
318 //------------------
319 jerror_t JEventProcessor_BCAL_point_calib::brun(JEventLoop *eventLoop, int32_t runnumber)
320 {
321  // This is called whenever the run number changes
322  Run_Number = runnumber;
323  //BCAL_Neutrals->Fill();
324  //cout << " run number = " << RunNumber << endl;
325 
326 
327  return NOERROR;
328 }
329 
330 //------------------
331 // evnt
332 //------------------
333 jerror_t JEventProcessor_BCAL_point_calib::evnt(JEventLoop *loop, uint64_t eventnumber)
334 {
335  // This is called for every event. Use of common resources like writing
336  // to a file or filling a histogram should be mutex protected. Using
337  // loop->Get(...) to get reconstructed objects (and thereby activating the
338  // reconstruction algorithm) should be done outside of any mutex lock
339  // since multiple threads may call this method at the same time.
340  // Here's an example:
341  //
342  // vector<const MyDataClass*> mydataclasses;
343  // loop->Get(mydataclasses);
344  //
345  // japp->RootFillLock(this);
346  // ... fill historgrams or trees ...
347  // japp->RootFillUnLock(this);
348 
349  // select events with physics events, i.e., not LED and other front panel triggers
350  const DTrigger* locTrigger = NULL;
351  loop->GetSingle(locTrigger);
352  if(locTrigger->Get_L1FrontPanelTriggerBits() != 0)
353  return NOERROR;
354 
355  // load BCAL geometry
356  vector<const DBCALGeometry *> BCALGeomVec;
357  loop->Get(BCALGeomVec);
358  if(BCALGeomVec.size() == 0)
359  throw JException("Could not load DBCALGeometry object!");
360  auto locBCALGeom = BCALGeomVec[0];
361 
362  vector<const DTrackFitter *> fitters;
363  loop->Get(fitters);
364 
365  if(fitters.size()<1){
366  _DBG_<<"Unable to get a DTrackFinder object!"<<endl;
367  return RESOURCE_UNAVAILABLE;
368  }
369 
370  auto fitter = fitters[0];
371 
372  vector<const DBCALShower*> locBCALShowers;
373  vector<const DBCALHit*> bcalhits;
374  vector<const DBCALPoint*> locBCALPoints;
375 
376  loop->Get(locBCALShowers);
377  loop->Get(bcalhits);
378  loop->Get(locBCALPoints);
379 
380  vector<const DTrackTimeBased*> locTrackTimeBased;
381  loop->Get(locTrackTimeBased);
382 
383  vector <const DBCALShower *> matchedShowers;
384  vector < const DBCALShower *> matchedShowersneg;
385  vector < const DBCALShower *> matchedShowerspos;
386  vector <const DTrackTimeBased* > matchedTracks;
387  DVector3 mypos(0.0,0.0,0.0);
388 
389  map< const DBCALShower*, vector<const DBCALPoint*> > matchedShowerPoints_cache;
390  map< const DBCALPoint *, const DBCALUnifiedHit *> upstreamHit_cache;
391  map< const DBCALPoint *, const DBCALUnifiedHit *> downstreamHit_cache;
392 
393  // the following two loops loop through a) the showers in each event and b) the tracks in each event and determine the z and phi of both. If the quantities dZ and dPhi are "small" enough, the relevant track and shower are appended to the end of the matchedTracks and matchedShowers vectors respectively
394 
395  for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
396  vector<DTrackFitter::Extrapolation_t>extrapolations=locTrackTimeBased[i]->extrapolations.at(SYS_BCAL);
397  if (extrapolations.size()==0) continue;
398 
399  for (unsigned int j=0; j< locBCALShowers.size(); ++j){
400 
401  double x = locBCALShowers[j]->x;
402  double y = locBCALShowers[j]->y;
403  double z = locBCALShowers[j]->z;
404 
405  DVector3 pos_bcal(x,y,z);
406  double R = pos_bcal.Perp();
407  double phi = pos_bcal.Phi();
408 
409  if (fitter->ExtrapolateToRadius(R,extrapolations,mypos)){
410  double dPhi = mypos.Phi()-phi;
411  if (dPhi<-M_PI) dPhi+=2.*M_PI;
412  if (dPhi>M_PI) dPhi-=2.*M_PI;
413 
414  double dZ = TMath::Abs(mypos.Z() - z);
415 
416  // save showers matched to tracks
417  if(dZ < 30.0 && fabs(dPhi) < 0.18) {
418  matchedShowers.push_back(locBCALShowers[j]);
419  matchedTracks.push_back(locTrackTimeBased[i]);
420  }
421  }
422  } // end of shower loop
423  } // end of track loop
424 
425 
426 
427  event_count++;
428  if (VERBOSE && (event_count%100 == 0)) printf ("Event count=%d, EventNumber=%lu\n",event_count,eventnumber);
429 
430  // vectors to store the intersection of the matched tracks with the middle of each BCAL layer
431  DVector3 mypos_1(0.0,0.0,0.0); // middle of Layer 1
432  DVector3 mypos_2(0.0,0.0,0.0); // middle of Layer 2
433  DVector3 mypos_3(0.0,0.0,0.0); // middle of Layer 3
434  DVector3 mypos_4(0.0,0.0,0.0); // middle of Layer 4
435 
436  // get the inner radius of each layer (and the outer radius of layer 4)
437  const float*bcal_radii = locBCALGeom->GetBCAL_radii();
438 
439  // the radii at the middle of each layer
440  double R1 = 0.5*(bcal_radii[0] + bcal_radii[1]);
441  double R2 = 0.5*(bcal_radii[1] + bcal_radii[2]);
442  double R3 = 0.5*(bcal_radii[2] + bcal_radii[3]);
443  double R4 = 0.5*(bcal_radii[3] + bcal_radii[4]);
444 
445  // loop over matched showers
446  Int_t numshowers_per_event = matchedShowers.size();
447  if (numshowers_per_event > 0)
448  h1_Num_matched_showers->Fill(numshowers_per_event);
449  Int_t numtracks_per_event = matchedTracks.size();
450  if (numtracks_per_event > 0)
451  h1trk_Num_matched_tracks->Fill(numtracks_per_event);
452 
453  // cache points for each matched points to speed up the locked loop below
454  for(int i=0; i<numshowers_per_event; i++) {
455  vector<const DBCALPoint*> points;
456  matchedShowers[i]->Get(points);
457  matchedShowerPoints_cache[ matchedShowers[i] ] = points;
458 
459  vector<const DBCALUnifiedHit*> hits;
460  for(unsigned int k=0; k<points.size(); k++) {
461  points[k]->Get(hits);
462 
463  if(hits.size() != 2)
464  continue;
465 
466  if(hits[0]->end == DBCALGeometry::kUpstream) {
467  upstreamHit_cache[ points[k] ] = hits[0];
468  downstreamHit_cache[ points[k] ] = hits[1];
469  } else {
470  upstreamHit_cache[ points[k] ] = hits[1];
471  downstreamHit_cache[ points[k] ] = hits[0];
472  }
473  }
474  }
475 
476  // FILL HISTOGRAMS
477  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
478  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
479 
480  for(int i=0; i<numshowers_per_event; i++){
481  if(DEBUG) {
482  // fill histogram information for matched showers
483  Float_t E = matchedShowers[i]->E;
484  Float_t E_raw = matchedShowers[i]->E_raw;
485  Float_t x = matchedShowers[i]->x;
486  Float_t y = matchedShowers[i]->y;
487  Float_t z = matchedShowers[i]->z;
488  Float_t t = matchedShowers[i]->t;
489  Int_t N_cell = matchedShowers[i]->N_cell;
490 
491  h1_E->Fill(E);
492  h1_E_raw->Fill(E_raw);
493  h1_x->Fill(x);
494  h1_y->Fill(y);
495  h1_z->Fill(z);
496  h1_t->Fill(t);
497  h1_N_cell->Fill(N_cell);
498 
499  // fill histograms related to matched track
500  double FOM = matchedTracks[i]->FOM;
501  double pmag = matchedTracks[i]->pmag();
502  double px = matchedTracks[i]->px();
503  double py = matchedTracks[i]->py();
504  double pz = matchedTracks[i]->pz();
505  double x_track = matchedTracks[i]->x();
506  double y_track = matchedTracks[i]->y();
507  double z_track = matchedTracks[i]->z();
508 
509  double energy = matchedTracks[i]->energy();
510  h1trk_FOM->Fill(FOM);
511  h1trk_pmag->Fill(pmag);
512  h1trk_px->Fill(px);
513  h1trk_py->Fill(py);
514  h1trk_pz->Fill(pz);
515  h1trk_x->Fill(x_track);
516  h1trk_y->Fill(y_track);
517  h1trk_z->Fill(z_track);
518  h1trk_energy->Fill(energy);
519 
520  h2_Evsenergy->Fill(energy,E);
521  h2_pmagvsenergy->Fill(energy,pmag);
522  }
523 
524  // get mypos for given radius R
525  vector<DTrackFitter::Extrapolation_t>extrapolations=matchedTracks[i]->extrapolations.at(SYS_BCAL);
526  fitter->ExtrapolateToRadius(R1,extrapolations,mypos_1); // middle of layer 1
527  fitter->ExtrapolateToRadius(R2,extrapolations,mypos_2); // middle of layer 2
528  fitter->ExtrapolateToRadius(R3,extrapolations,mypos_3); // middle of layer 3
529  fitter->ExtrapolateToRadius(R4,extrapolations,mypos_4); // middle of layer 4
530 
531  // create an array with the values of z_of_track for each layer
532  double mypos_z_array[4] = {};
533  mypos_z_array[0] = mypos_1.Z(); // Layer 1
534  mypos_z_array[1] = mypos_2.Z(); // Layer 2
535  mypos_z_array[2] = mypos_3.Z(); // Layer 3
536  mypos_z_array[3] = mypos_4.Z(); // Layer 4
537 
538  double target_center = 65.0;
539 
540  // get associated information in the shower
541 
542  //vector<const DBCALPoint*> points;
543  //matchedShowers[i]->Get(points);
544  vector<const DBCALPoint*> &points = matchedShowerPoints_cache[ matchedShowers[i] ];
545 
546  // loop over points in shower
547  Int_t numpoints_per_shower = points.size();
548  h1pt_Num_points->Fill(numpoints_per_shower);
549 
550  for (int jj=0; jj<numpoints_per_shower; jj++) {
551  int module = points[jj]->module();
552  int layer = points[jj]->layer();
553  int sector = points[jj]->sector();
554  int cell_id = (module-1)*16 + (layer-1)*4 + sector;
555  double energy = points[jj]->E();
556  double energy_US = points[jj]->E_US();
557  double energy_DS = points[jj]->E_DS();
558  double z_point_wrt_target_center = points[jj]->z();
559  double z_point = z_point_wrt_target_center + target_center;
560  double r_point = points[jj]->r();
561  double phi_point = points[jj]->phi();
562  double t_point = points[jj]->t();
563 
564 
565  const DBCALUnifiedHit *up_hit = upstreamHit_cache[ points[jj] ];
566  const DBCALUnifiedHit *down_hit = downstreamHit_cache[ points[jj] ];
567  double deltat_point = up_hit->t_ADC - down_hit->t_ADC;
568 
569 
570  if(DEBUG) {
571  // fill histograms with point information
572  h1pt_module->Fill(module);
573  h1pt_layer->Fill(layer);
574  h1pt_sector->Fill(sector);
575  h1pt_cell_id->Fill(cell_id);
576  h1pt_energy->Fill(energy);
577  h1pt_energy_US->Fill(energy_US);
578  h1pt_energy_DS->Fill(energy_DS);
579  h1pt_z->Fill(z_point);
580  h1pt_r->Fill(r_point);
581  h1pt_phi->Fill(phi_point);
582  h1pt_t->Fill(t_point);
583  }
584 
585  // accept only positive z-coordinates
586  if(mypos_z_array[layer-1] > 0 && z_point > 0){
587  if(DEBUG) {
588  // fill histograms for z_track and z_point correlation
589  h1_ztrack_minus_zpoint->Fill(mypos_z_array[layer-1] - z_point);
590  h1_abs_ztrack_minus_zpoint->Fill(TMath::Abs(mypos_z_array[layer-1] - z_point));
591  h2_ztrack_minus_zpoint_vs_ztrack->Fill(mypos_z_array[layer-1], mypos_z_array[layer-1] - z_point);
592 
593  h2_zpoint_vs_ztrack->Fill(mypos_z_array[layer-1],z_point);
594  }
595 
596  // introduce |z_track - z_point| cuts here, if needed
597  if(TMath::Abs(mypos_z_array[layer-1] - z_point) < 400.0){
598  // fill z_track vs z_point graphs
599  z_track_graph[cell_id-1].push_back(mypos_z_array[layer-1]);
600  z_point_graph[cell_id-1].push_back(z_point);
601  z_deltat_graph[cell_id-1].push_back(deltat_point);
602  }
603  else
604  {
605  // fill graphs of rejected points, if cuts were applied
606  z_track_thrown[cell_id-1].push_back(mypos_z_array[layer-1]);
607  z_point_thrown[cell_id-1].push_back(z_point);
608  h2_zpoint_vs_ztrack_thrown->Fill(mypos_z_array[layer-1],z_point);
609  }
610  }
611 
612  } // end loop over points in shower
613 
614  } // end loop over matched showers
615 
616  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
617 
618  /*
619  //Optional: Save event to output REST file. Use this to create skims.
620  dEventWriterREST->Write_RESTEvent(loop, "BCAL_Shower"); //string is part of output file name
621  */
622 
623 
624  return NOERROR;
625 }
626 
627 //------------------
628 // erun
629 //------------------
631 {
632  // This is called whenever the run number changes, before it is
633  // changed to give you a chance to clean up before processing
634  // events from the next run number.
635  return NOERROR;
636 }
637 
638 //------------------
639 // fini
640 //------------------
642 {
643  // Called before program exit after event processing is finished.
644  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
645 
646  // lock all root operations
647  // japp->RootWriteLock();
648 
649  TDirectory *main = gDirectory;
650  TDirectory *locDirectory = (TDirectory*)gDirectory->FindObjectAny("bcal_point_calibs");
651 
652  if(locDirectory)
653  locDirectory->cd();
654 
655  int channels = 768;
656 
657  // array to store the percentage of thrown points per channel
658  double percent[768] = {0};
659  double number_of_points = 0.0;
660  double number_of_thrown_points = 0.0;
661 
662  // loop over channels (cells) and fill graphs
663  for(int m = 0; m < channels; ++m){
664  h2_tgraph[m] = new TGraph(z_track_graph[m].size(), &(z_track_graph[m][0]), &(z_point_graph[m][0]));
665  h2_tgraph[m]->SetTitle(Form("z of point vs z of track for channel %i", m+1));
666  h2_tgraph[m]->GetXaxis()->SetTitle("z of track (cm)");
667  h2_tgraph[m]->GetYaxis()->SetTitle("z of point (cm)");
668  h2_tgraph[m]->GetXaxis()->SetLimits(0.0, 500.0);
669  h2_tgraph[m]->GetYaxis()->SetRangeUser(0.0, 500.0);
670  h2_tgraph[m]->SetMarkerStyle(20);
671 
672  h2_tgraph[m]->Draw("A");
673  h2_tgraph[m]->Write(Form("h2_tgraph[%i]", m+1));
674 
675  h2_tgraph_deltat[m] = new TGraph(z_track_graph[m].size(), &(z_deltat_graph[m][0]), &(z_track_graph[m][0]));
676  h2_tgraph_deltat[m]->SetTitle(Form("Delta T of point vs z of track for channel %i", m+1));
677  h2_tgraph_deltat[m]->GetXaxis()->SetTitle("Delta T of point (cm)");
678  h2_tgraph_deltat[m]->GetYaxis()->SetTitle("z of track (cm)");
679  h2_tgraph_deltat[m]->GetXaxis()->SetLimits(-10.0, 10.0);
680  h2_tgraph_deltat[m]->GetYaxis()->SetRangeUser(0.0, 500.0);
681  h2_tgraph_deltat[m]->SetMarkerStyle(20);
682 
683  h2_tgraph_deltat[m]->Draw("A");
684  h2_tgraph_deltat[m]->Write(Form("h2_tgraph_deltat[%i]", m+1));
685 
686  // graphs of the thrown points, if cuts were applied
687  h2_tgraph_thrown[m] = new TGraph(z_track_thrown[m].size(), &(z_track_thrown[m][0]), &(z_point_thrown[m][0]));
688  h2_tgraph_thrown[m]->SetTitle(Form("Thrown z of point vs z of track for channel %i", m+1));
689  h2_tgraph_thrown[m]->GetXaxis()->SetTitle("z of track (cm)");
690  h2_tgraph_thrown[m]->GetYaxis()->SetTitle("z of point (cm)");
691  h2_tgraph_thrown[m]->GetXaxis()->SetLimits(0.0, 500.0);
692  h2_tgraph_thrown[m]->GetYaxis()->SetRangeUser(0.0, 500.0);
693  h2_tgraph_thrown[m]->SetMarkerStyle(20);
694  h2_tgraph_thrown[m]->Draw("A");
695  h2_tgraph_thrown[m]->Write(Form("h2_tgraph_thrown[%i]", m+1));
696 
697  // fill histos with percentages of thrown points per channel
698  number_of_points = h2_tgraph[m]->GetN();
699  number_of_thrown_points = h2_tgraph_thrown[m]->GetN();
700  percent[m] = 100 * (number_of_thrown_points / number_of_points);
701  if(percent[m] >= 0.0 && percent[m] <= 100.0){
702  h1_thrown_per_channel->Fill(m+1, percent[m]);
703  }
704  }
705 
706  main->cd();
707 
708  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
709 
710  // unlock
711  //japp->RootUnLock();
712 
713  return NOERROR;
714 }
715 
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
bool ExtrapolateToRadius(double R, const vector< Extrapolation_t > &extraps, DVector3 &pos, DVector3 &mom, double &t, double &s) const
vector< vector< double > > z_track_thrown(768, vector< double >())
static TH1I * h1trk_z
static TH1I * h1pt_sector
static TH1I * h1_ztrack_minus_zpoint
static TH1D * h1_thrown_per_channel
Int_t layer
static TH1I * h1pt_r_size
static TH1I * h1_y
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
static TH1I * h1_t
vector< vector< double > > z_track_graph(768, vector< double >())
#define y
static TH1I * h1pt_t
static TH1I * h1pt_energy_DS
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
static TH1I * h1_E
uint32_t Get_L1FrontPanelTriggerBits(void) const
vector< vector< double > > z_deltat_graph(768, vector< double >())
static TH1I * h1_x
static TH1I * h1pt_layer
static TH1I * h1pt_module
static TGraph * h2_tgraph_thrown[768]
static TH2I * h2_pmagvsenergy
static TH1I * h1pt_r
static TH2I * h2_zpoint_vs_ztrack_thrown
Definition: GlueX.h:19
JApplication * japp
TDirectory * locDirectory
static TH1I * h1pt_energy
static TH1I * h1trk_py
static TH2I * h2_Evsenergy
InitPlugin_t InitPlugin
const bool VERBOSE
static TH1I * h1pt_phi
static TH1I * h1pt_Num_points
vector< vector< double > > z_point_graph(768, vector< double >())
jerror_t init(void)
Called once at program start.
#define _DBG_
Definition: HDEVIO.h:12
float t_ADC
Time from fADC.
static TGraph * h2_tgraph[768]
TH1D * py[NCHANNELS]
static TH1I * h1pt_energy_US
static TH1I * h1trk_y
static TH1I * h1trk_Num_matched_tracks
static TH1I * h1pt_z
static TH1I * h1trk_energy
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
static TH2I * h2_zpoint_vs_ztrack
static TH1I * h1_N_cell
const DTrackFitter * fitter
jerror_t fini(void)
Called after last event of last event source has been processed.
static TH1I * h1_E_raw
static TH1I * h1trk_pmag
static TGraph * h2_tgraph_deltat[768]
static TH1I * h1trk_FOM
printf("string=%s", string)
static TH1I * h1_abs_ztrack_minus_zpoint
vector< vector< double > > z_point_thrown(768, vector< double >())
static TH2I * h2_ztrack_minus_zpoint_vs_ztrack
static TH1I * h1_z
#define DEBUG
Definition: tw_corr.C:41
int main(int argc, char *argv[])
Definition: gendoc.cc:6
static TH1I * h1trk_x
static TH1I * h1_Num_matched_showers
static TH1I * h1trk_pz
static TH1I * h1trk_px
static TH1I * h1pt_cell_id