Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_BCAL_Eff.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_BCAL_Eff.cc Copied from Will McGinley
4 // Created: Fri Oct 10 16:41:18 EDT 2014
5 // Creator: wmcginle (on Linux ifarm1101 2.6.32-220.7.1.el6.x86_64 x86_64)
6 //
7 
8 #include "JEventProcessor_BCAL_Eff.h"
9 
10 #include <TLorentzVector.h>
11 #include "TMath.h"
12 
13 #include "DANA/DApplication.h"
14 #include "BCAL/DBCALShower.h"
15 #include "BCAL/DBCALTruthShower.h"
16 //#include "BCAL/DBCALCluster.h"
17 #include "BCAL/DBCALPoint.h"
18 #include "BCAL/DBCALHit.h"
19 #include "FCAL/DFCALShower.h"
20 #include "TRACKING/DMCThrown.h"
21 #include "ANALYSIS/DAnalysisUtilities.h"
22 //#include "TRACKING/DTrackFinder.h"
23 #include "TRIGGER/DTrigger.h"
24 
25 #include <vector>
26 #include <deque>
27 #include <string>
28 #include <iostream>
29 #include <algorithm>
30 #include <stdio.h>
31 #include <stdlib.h>
32 
33 
34 
35  static TH1I* h1_Num_matched_showers = NULL;
36  /*static TH1I* h1_E = NULL;
37  static TH1I* h1_E_raw = NULL;
38  static TH1I* h1_x = NULL;
39  static TH1I* h1_y = NULL;
40  static TH1I* h1_z = NULL;
41  static TH1I* h1_t = NULL;
42  static TH1I* h1_xErr = NULL;
43  static TH1I* h1_yErr = NULL;
44  static TH1I* h1_zErr = NULL;
45  static TH1I* h1_tErr = NULL;
46  static TH1I* h1_N_cell = NULL;*/
47 
48  static TH1I* h1trk_Num_matched_tracks = NULL;
49  static TH1I* h1trk_FOM = NULL;
50  static TH1I* h1trk_pmag = NULL;
51  static TH1I* h1trk_px = NULL;
52  static TH1I* h1trk_py = NULL;
53  static TH1I* h1trk_pz = NULL;
54  static TH1I* h1trk_energy = NULL;
55 
56  static TH2I* h2_Evsenergy = NULL;
57  static TH2I* h2_pmagvsenergy = NULL;
58 
59  //static TH1I* h1shower_Num_matched_showers = NULL;
60 
61 
62  static TH1I* h1pt_Num_points = NULL;
63  static TH1I* h1pt_module = NULL;
64  static TH1I* h1pt_layer = NULL;
65  static TH1I* h1pt_sector = NULL;
66  static TH1I* h1pt_cell_id = NULL;
67  static TH1I* h1pt_energy = NULL;
68  static TH1I* h1pt_energy_US = NULL;
69  static TH1I* h1pt_energy_DS = NULL;
70 
71  static TH1I* h1eff_layertot = NULL;
72  static TH1I* h1eff_layer = NULL;
73  static TH1F* h1eff_eff = NULL;
74 
75  static TH1I* h1eff_cellidtot = NULL;
76  static TH1I* h1eff_cellid = NULL;
77  static TH1F* h1eff_cellideff = NULL;
78 
79  static TH1I* h1eff2_layertot = NULL;
80  static TH1I* h1eff2_layer = NULL;
81  static TH1F* h1eff2_eff2 = NULL;
82 
83  static TH1I* h1eff2_cellidtot = NULL;
84  static TH1I* h1eff2_cellid = NULL;
85  static TH1F* h1eff2_cellideff2 = NULL;
86 
87  static TH1I* h1E_US_layer1 = NULL;
88  static TH1I* h1E_DS_layer1 = NULL;
89  static TH2I* h2E_USvsDS_layer1 = NULL;
90 
91  int event_count=0;
92 
93 
94 // Routine used to create our JEventProcessor
95 
96 extern "C"
97 {
98  void InitPlugin(JApplication *locApplication)
99  {
100  InitJANAPlugin(locApplication);
101  locApplication->AddProcessor(new JEventProcessor_BCAL_Eff()); //register this plugin
102  }
103 } // "C"
104 
105 //------------------
106 // init
107 //------------------
108 jerror_t JEventProcessor_BCAL_Eff::init(void)
109 {
110  // First thread to get here makes all histograms. If one pointer is
111  // already not NULL, assume all histograms are defined and return now
112  if(h1_Num_matched_showers != NULL){
113  return NOERROR;
114  }
115 
116  // ... create historgrams or trees ...
117 
118  // TDirectory *dir = new TDirectoryFile("BCAL","BCAL");
119  // dir->cd();
120  // create root folder for bcal and cd to it, store main dir
121  TDirectory *main = gDirectory;
122  gDirectory->mkdir("bcal_eff")->cd();
123 
124 
125  // double pi0_mass = 0.1349766;
126  int nbins=100;
127 
128  h1_Num_matched_showers = new TH1I("h1_Num_matched_showers", "Number of matched showers per event", 20,0,20);
129  h1_Num_matched_showers->SetXTitle("Number of showers");
130  h1_Num_matched_showers->SetYTitle("counts");
131 
132  /*h1_E = new TH1I("h1_E", "Matched energy per shower",nbins,0,2);
133  h1_E->SetXTitle("Energy per shower (GeV)");
134  h1_E->SetYTitle("counts");
135  h1_E_raw = new TH1I("h1_E_raw", "Matched raw energy per shower",nbins,0,2);
136  h1_E_raw->SetXTitle("Raw Energy per shower (GeV)");
137  h1_E_raw->SetYTitle("counts");
138  h1_x = new TH1I("h1_x", "x pos of shower",nbins,-100,100);
139  h1_x->SetXTitle("x (cm)");
140  h1_x->SetYTitle("counts");
141  h1_y = new TH1I("h1_y", "y pos of shower",nbins,-100,100);
142  h1_y->SetXTitle("y (cm)");
143  h1_y->SetYTitle("counts");
144  h1_z = new TH1I("h1_z", "z pos of shower",nbins,0,500);
145  h1_z->SetXTitle("z (cm)");
146  h1_z->SetYTitle("counts");
147  h1_t = new TH1I("h1_t", "time of shower",nbins,-100,100);
148  h1_t->SetXTitle("time of shower (ns)");
149  h1_t->SetYTitle("counts");
150  h1_xErr = new TH1I("h1_xErr", "x err of shower",nbins,0,20);
151  h1_xErr->SetXTitle("x err (cm)");
152  h1_xErr->SetYTitle("counts");
153  h1_yErr = new TH1I("h1_yErr", "y err of shower",nbins,0,20);
154  h1_yErr->SetXTitle("y err (cm)");
155  h1_yErr->SetYTitle("counts");
156  h1_zErr = new TH1I("h1_zErr", "z err of shower",nbins,0,20);
157  h1_zErr->SetXTitle("z err (cm)");
158  h1_zErr->SetYTitle("counts");
159  h1_tErr = new TH1I("h1_tErr", "t err of shower",nbins,0,5);
160  h1_tErr->SetXTitle("t err (ns)");
161  h1_tErr->SetYTitle("counts");
162  h1_N_cell = new TH1I("h1_N_cell", "N_cell of shower",20,0,20);
163  h1_N_cell->SetXTitle("N_cell");
164  h1_N_cell->SetYTitle("counts");*/
165 
166  h1trk_Num_matched_tracks = new TH1I("h1trk_Num_matched_tracks", "Number of matched tracks per event", 20,0,20);
167  h1trk_Num_matched_tracks->SetXTitle("Number of tracks");
168  h1trk_Num_matched_tracks->SetYTitle("counts");
169  h1trk_FOM = new TH1I("h1trk_FOM", "FOM for matched tracks", nbins,0,1);
170  h1trk_FOM->SetXTitle("FOM");
171  h1trk_FOM->SetYTitle("counts");
172  h1trk_pmag = new TH1I("h1trk_mag", "pmag for matched tracks", nbins,0,4);
173  h1trk_pmag->SetXTitle("pmag (GeV)");
174  h1trk_pmag->SetYTitle("counts");
175  h1trk_px = new TH1I("h1trk_px", "px for matched tracks", nbins,0,2);
176  h1trk_px->SetXTitle("px (GeV)");
177  h1trk_px->SetYTitle("counts");
178  h1trk_py = new TH1I("h1trk_py", "py for matched tracks", nbins,0,2);
179  h1trk_py->SetXTitle("Number of tracks");
180  h1trk_py->SetYTitle("counts");
181  h1trk_pz = new TH1I("h1trk_pz", "pz for matched tracks", nbins,0,4);
182  h1trk_pz->SetXTitle("pz (GeV)");
183  h1trk_pz->SetYTitle("counts");
184  h1trk_energy = new TH1I("h1trk_energy", "energy for matched tracks", nbins,0,4);
185  h1trk_energy->SetXTitle("energy (GeV)");
186  h1trk_energy->SetYTitle("counts");
187 
188  h2_Evsenergy = new TH2I("h2_Evsenergy", "E vs energy matched tracks", nbins,0,4,nbins,0,4);
189  h2_Evsenergy->SetXTitle("track energy (GeV)");
190  h2_Evsenergy->SetYTitle("Eshower (GeV)");
191  h2_pmagvsenergy = new TH2I("h2_pmagvsenergy", "pmag vs energy matched tracks", nbins,0,4,nbins,0,4);
192  h2_pmagvsenergy->SetXTitle("track energy (GeV)");
193  h2_pmagvsenergy->SetYTitle("track pmag (GeV)");
194 
195 
196  h1pt_Num_points = new TH1I("h1pt_Num_points", "pt points per Shower",20,0,20);
197  h1pt_Num_points->SetXTitle("Points per Shower");
198  h1pt_Num_points->SetYTitle("counts");
199  h1pt_module = new TH1I("h1pt_module", "pt module number",50,0,50);
200  h1pt_module->SetXTitle("Module");
201  h1pt_module->SetYTitle("counts");
202  h1pt_layer = new TH1I("h1pt_layer", "pt layer number",5,0,5);
203  h1pt_layer->SetXTitle("Layer");
204  h1pt_layer->SetYTitle("counts");
205  h1pt_sector = new TH1I("h1pt_sector", "pt sector number",5,0,5);
206  h1pt_sector->SetXTitle("Sector");
207  h1pt_sector->SetYTitle("counts");
208  h1pt_cell_id = new TH1I("h1pt_cell_id", "pt cell id number",800,0,800);
209  h1pt_cell_id->SetXTitle("Cell ID");
210  h1pt_cell_id->SetYTitle("counts");
211  h1pt_energy = new TH1I("h1pt_energy", "pt energy",nbins,0,1);
212  h1pt_energy->SetXTitle("Point Energy");
213  h1pt_energy->SetYTitle("counts");
214  h1pt_energy_US = new TH1I("h1pt_energy_US", "pt energy US",nbins,0,1);
215  h1pt_energy_US->SetXTitle("Point Energy US");
216  h1pt_energy_US->SetYTitle("counts");
217  h1pt_energy_DS = new TH1I("h1pt_energy_DS", "pt energy DS",nbins,0,1);
218  h1pt_energy_DS->SetXTitle("Point Energy DS");
219  h1pt_energy_DS->SetYTitle("counts");
220 
221  h1eff_layertot = new TH1I("h1eff_layertot", "eff layertot",5,0,5);
222  h1eff_layertot->SetXTitle("Layer");
223  h1eff_layertot->SetYTitle("counts");
224  //h1eff_layertot->Sumw2();
225  h1eff_layer = new TH1I("h1eff_layer", "eff layer hit",5,0,5);
226  h1eff_layer->SetXTitle("Layer hit");
227  h1eff_layer->SetYTitle("counts");
228  //h1eff_layer->Sumw2();
229  h1eff_eff = new TH1F("h1eff_eff", "efficiency",5,0,5);
230  h1eff_eff->SetXTitle("Layer");
231  h1eff_eff->SetYTitle("Efficiency");
232  h1eff_eff->Sumw2();
233 
234  h1eff_cellidtot = new TH1I("h1eff_cellidtot", "eff cellidtot",200,0,200);
235  h1eff_cellidtot->SetXTitle("Cellid");
236  h1eff_cellidtot->SetYTitle("counts");
237  //h1eff_cellidtot->Sumw2();
238  h1eff_cellid = new TH1I("h1eff_cellid", "eff cellid hit",200,0,200);
239  h1eff_cellid->SetXTitle("Cellid hit");
240  h1eff_cellid->SetYTitle("counts");
241  //h1eff_cellid->Sumw2();
242  h1eff_cellideff = new TH1F("h1eff_cellideff", "efficiency",200,0,200);
243  h1eff_cellideff->SetXTitle("Cellid");
244  h1eff_cellideff->SetYTitle("Efficiency");
245  h1eff_cellideff->Sumw2();
246 
247  h1eff2_layertot = new TH1I("h1eff2_layertot", "eff2 layertot",5,0,5);
248  h1eff2_layertot->SetXTitle("Layer");
249  h1eff2_layertot->SetYTitle("counts");
250  //h1eff2_layertot->Sumw2();
251  h1eff2_layer = new TH1I("h1eff2_layer", "eff2 layer hit",5,0,5);
252  h1eff2_layer->SetXTitle("Layer hit");
253  h1eff2_layer->SetYTitle("counts");
254  //h1eff2_layer->Sumw2();
255  h1eff2_eff2 = new TH1F("h1eff2_eff2", "eff2 efficiency",5,0,5);
256  h1eff2_eff2->SetXTitle("Layer");
257  h1eff2_eff2->SetYTitle("Efficiency");
258  h1eff2_eff2->Sumw2();
259 
260  h1eff2_cellidtot = new TH1I("h1eff2_cellidtot", "eff2 cellidtot",200,0,200);
261  h1eff2_cellidtot->SetXTitle("Cellid");
262  h1eff2_cellidtot->SetYTitle("counts");
263  //h1eff2_cellidtot->Sumw2();
264  h1eff2_cellid = new TH1I("h1eff2_cellid", "eff2 cellid hit",200,0,200);
265  h1eff2_cellid->SetXTitle("Cellid hit");
266  h1eff2_cellid->SetYTitle("counts");
267  //h1eff2_cellid->Sumw2();
268  h1eff2_cellideff2 = new TH1F("h1eff2_cellideff2", "eff2 efficiency",200,0,200);
269  h1eff2_cellideff2->SetXTitle("Cellid");
270  h1eff2_cellideff2->SetYTitle("Efficiency");
271  h1eff2_cellideff2->Sumw2();
272 
273  h1E_US_layer1 = new TH1I("h1E_US_layer1", "E US layer1",nbins,0,1);
274  h1E_US_layer1->SetXTitle("Upstream E (GeV)");
275  h1E_US_layer1->SetYTitle("Counts");
276  h1E_DS_layer1 = new TH1I("h1E_DS_layer1", "E DS layer1",nbins,0,1);
277  h1E_DS_layer1->SetXTitle("Downstream E (GeV)");
278  h1E_DS_layer1->SetYTitle("Counts");
279  h2E_USvsDS_layer1 = new TH2I("h2E_USvsDS_layer1", "E US vs DS ",nbins,0,1,nbins,0,1);
280  h2E_USvsDS_layer1->SetXTitle("Downstream E (GeV)");
281  h2E_USvsDS_layer1->SetYTitle("Upstream E (GeV)");
282 
283 
284 
285  // back to main dir
286  main->cd();
287 
288  return NOERROR;
289 }
290 
291 
292 //------------------
293 // brun
294 //------------------
295 jerror_t JEventProcessor_BCAL_Eff::brun(jana::JEventLoop* locEventLoop, int locRunNumber)
296 {
297  // This is called whenever the run number changes
298  Run_Number = locRunNumber;
299  //BCAL_Neutrals->Fill();
300  //cout << " run number = " << RunNumber << endl;
301  /*
302  //Optional: Retrieve REST writer for writing out skims
303  locEventLoop->GetSingle(dEventWriterREST);
304  */
305 
306  //vector<const DTrackFinder *> finders;
307  //locEventLoop->Get(finders);
308  //finder = const_cast<DTrackFinder*>(finders[0]);
309 
310  /*
311  //Recommeded: Create output ROOT TTrees (nothing is done if already created)
312  locEventLoop->GetSingle(dEventWriterROOT);
313  dEventWriterROOT->Create_DataTrees(locEventLoop);
314  */
315  /*
316  vector<const DTrackFitter *> fitters;
317  locEventLoop->Get(fitters);
318 
319  if(fitters.size()<1){
320  _DBG_<<"Unable to get a DTrackFinder object!"<<endl;
321  return RESOURCE_UNAVAILABLE;
322  }
323 
324  fitter = fitters[0];
325  */
326 
327  return NOERROR;
328 }
329 
330 //------------------
331 // evnt
332 //------------------
333 
334 
335 jerror_t JEventProcessor_BCAL_Eff::evnt(jana::JEventLoop* locEventLoop, uint64_t locEventNumber)
336 {
337 
338  // This is called for every event. Use of common resources like writing
339  // to a file or filling a histogram should be mutex protected. Using
340  // locEventLoop->Get(...) to get reconstructed objects (and thereby activating the
341  // reconstruction algorithm) should be done outside of any mutex lock
342  // since multiple threads may call this method at the same time.
343  //
344  // Here's an example:
345  //
346  // vector<const MyDataClass*> mydataclasses;
347  // locEventLoop->Get(mydataclasses);
348  //
349  // japp->RootWriteLock();
350  // ... fill historgrams or trees ...
351  // japp->RootUnLock();
352 
353  // DOCUMENTATION:
354  // ANALYSIS library: https://halldweb1.jlab.org/wiki/index.php/GlueX_Analysis_Software
355 
356  // select events with physics events, i.e., not LED and other front panel triggers
357  const DTrigger* locTrigger = NULL;
358  locEventLoop->GetSingle(locTrigger);
359  if(locTrigger->Get_L1FrontPanelTriggerBits() != 0)
360  return NOERROR;
361 
362  vector<const DTrackFitter *> fitters;
363  locEventLoop->Get(fitters);
364 
365  if(fitters.size()<1){
366  _DBG_<<"Unable to get a DTrackFinder object!"<<endl;
367  return RESOURCE_UNAVAILABLE;
368  }
369 
370  const DTrackFitter *fitter = fitters[0];
371 
372  vector<const DBCALShower*> locBCALShowers;
373  vector<const DBCALHit*> bcalhits;
374  vector<const DBCALPoint*> locBCALPoints;
375  vector<const DVertex*> kinfitVertex;
376  //const DDetectorMatches* locDetectorMatches = NULL;
377  //locEventLoop->GetSingle(locDetectorMatches);
378  locEventLoop->Get(locBCALShowers);
379  locEventLoop->Get(bcalhits);
380  locEventLoop->Get(locBCALPoints);
381  locEventLoop->Get(kinfitVertex);
382 
383  vector<const DTrackTimeBased*> locTrackTimeBased;
384  locEventLoop->Get(locTrackTimeBased);
385 
386  vector <const DBCALShower *> matchedShowers;
387  vector < const DBCALShower *> matchedShowersneg;
388  vector < const DBCALShower *> matchedShowerspos;
389  vector <const DTrackTimeBased* > matchedTracks;
390  DVector3 mypos(0.0,0.0,0.0);
391 
392  for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
393  vector<DTrackFitter::Extrapolation_t>extrapolations=locTrackTimeBased[i]->extrapolations.at(SYS_BCAL);
394  if (extrapolations.size()>0){
395  for (unsigned int j=0; j< locBCALShowers.size(); ++j){
396  double x = locBCALShowers[j]->x;
397  double y = locBCALShowers[j]->y;
398  double z = locBCALShowers[j]->z;
399  // double E = locBCALShowers[j]->E;
400  DVector3 pos_bcal(x,y,z);
401  double R = pos_bcal.Perp();
402  double phi = pos_bcal.Phi();
403  if (fitter->ExtrapolateToRadius(R,extrapolations,mypos)){
404 
405  // double q = locTrackTimeBased[i]->rt->q;
406  locTrackTimeBased[i]->momentum().Mag();
407  double trkmass = locTrackTimeBased[i]->mass();
408  //double dPhi = TMath::Abs(mypos.Phi()-phi);
409  double dPhi = mypos.Phi()-phi;
410  if (dPhi<-M_PI) dPhi+=2.*M_PI;
411  if (dPhi>M_PI) dPhi-=2.*M_PI;
412  double dZ = TMath::Abs(mypos.Z() - z);
413 
414  // save showers matched to tracks
415 
416  // assume program is run with plugin -PTRKFIT:MASS_HYPOTHESES_NEGATIVE=0.13957 -PTRKFIT:MASS_HYPOTHESES_POSITIVE=0.13957
417  // if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R) {
418  if(dZ < 30.0 && fabs(dPhi) < 0.18 && trkmass < 0.15) { // select pion hypothesis only
419  matchedShowers.push_back(locBCALShowers[j]);
420  matchedTracks.push_back(locTrackTimeBased[i]);
421  // printf ("Matched event=%d, i=%d, j=%d, p=%f, Ztrk=%f Zshr=%f, Phitrk=%f Phishw=%f\n",locEventNumber,i,j,p,mypos.Z(),z,mypos.Phi(),pos_bcal.Phi());
422  }
423  }
424  }
425  }
426  }
427 
428  // FILL HISTOGRAMS
429  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
430  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
431 
432  event_count++;
433 // if (event_count%100 == 0) printf ("Event count=%d, EventNumber=%d\n",event_count,locEventNumber);
434 
435 
436  // Get kinematic-fitted vertex
437 
438 #if 0 // disabling since it doesn't do anything other than cause compiler warnings 10/15/2015 DL
439  double kinfitVertexX, kinfitVertexY, kinfitVertexZ, kinfitVertexT;
440  for (unsigned int i = 0 ; i < kinfitVertex.size(); i++)
441  {
442  kinfitVertexX = kinfitVertex[i]->dSpacetimeVertex.X();
443  kinfitVertexY = kinfitVertex[i]->dSpacetimeVertex.Y();
444  kinfitVertexZ = kinfitVertex[i]->dSpacetimeVertex.Z();
445  kinfitVertexT = kinfitVertex[i]->dSpacetimeVertex.T();
446  // goodVertexZ->Fill(kinfitVertexZ);
447  }
448 #endif
449 
450  // now loop over matched showers
451 
452  Int_t numshowers_per_event = matchedShowers.size();
453  if (numshowers_per_event > 0) h1_Num_matched_showers->Fill(numshowers_per_event);
454  Int_t numtracks_per_event = matchedTracks.size();
455  if (numtracks_per_event > 0) h1trk_Num_matched_tracks->Fill(numtracks_per_event);
456  // if (numtracks_per_event > 0) printf ("Event=%d nshower=%d ntracks=%d\n\n",locEventNumber,numshowers_per_event,numtracks_per_event);
457 
458  for(int i=0; i<numshowers_per_event; i++){
459 
460  // fill histogram information for matched showers
461  Float_t E = matchedShowers[i]->E;
462  /*Float_t E_raw = matchedShowers[i]->E_raw;
463  Float_t x = matchedShowers[i]->x;
464  Float_t y = matchedShowers[i]->y;
465  Float_t z = matchedShowers[i]->z;
466  Float_t t = matchedShowers[i]->t;
467  Float_t xErr = matchedShowers[i]->xErr;
468  Float_t yErr = matchedShowers[i]->yErr;
469  Float_t zErr = matchedShowers[i]->zErr;
470  Float_t tErr = matchedShowers[i]->tErr;
471  Int_t N_cell = matchedShowers[i]->N_cell;
472 
473  h1_E->Fill(E);
474  h1_E_raw->Fill(E_raw);
475  h1_x->Fill(x);
476  h1_y->Fill(y);
477  h1_z->Fill(z);
478  h1_t->Fill(t);
479  h1_xErr->Fill(xErr);
480  h1_yErr->Fill(yErr);
481  h1_zErr->Fill(zErr);
482  h1_tErr->Fill(tErr);
483  h1_N_cell->Fill(N_cell);*/
484 
485  // fill histograms related to matched track
486 
487  double FOM = matchedTracks[i]->FOM;
488  double pmag = matchedTracks[i]->pmag();
489  double px = matchedTracks[i]->px();
490  double py = matchedTracks[i]->py();
491  double pz = matchedTracks[i]->pz();
492  double energy = matchedTracks[i]->energy();
493  h1trk_FOM->Fill(FOM);
494  h1trk_pmag->Fill(pmag);
495  h1trk_px->Fill(px);
496  h1trk_py->Fill(py);
497  h1trk_pz->Fill(pz);
498  h1trk_energy->Fill(energy);
499 
500  h2_Evsenergy->Fill(energy,E);
501  h2_pmagvsenergy->Fill(energy,pmag);
502 
503 
504  vector <int> matchedShowers_modules; // vector with list of modules in Showers
505 
506 
507  // get associated information in the shower
508 
509  vector<const DBCALPoint*> points;
510  matchedShowers[i]->Get(points);
511  matchedShowers_modules.clear();
512 
513  // loop over points in shower
514  Int_t numpoints_per_shower = points.size();
515  h1pt_Num_points->Fill(numpoints_per_shower);
516 
517  // initialize hit variables
518  Int_t pointlayer1=0;
519  Int_t pointlayer2=0;
520  Int_t pointlayer3=0;
521  Int_t pointlayer4=0;
522  Int_t module_ndx=0; // use one of the modules in shower as index
523  for (int jj=0; jj<numpoints_per_shower; jj++) {
524 
525  int module = points[jj]->module();
526  int layer = points[jj]->layer();
527  int sector = points[jj]->sector();
528  // int cell_id = DBCALGeometry::cellId(module, layer, sector);
529  int cell_id = (module-1)*16 + (layer-1)*4 + sector;
530  // cout << " cell ID = " << cell_id << " module = " << module1 << " layer = " << layer1 << " sector = " << sector1 << endl;
531  double energy = points[jj]->E();
532  double energy_US = points[jj]->E_US();
533  double energy_DS = points[jj]->E_DS();
534 
535  // fill histograms with point information
536  h1pt_module->Fill(module);
537  h1pt_layer->Fill(layer);
538  h1pt_sector->Fill(sector);
539  h1pt_cell_id->Fill(cell_id);
540  h1pt_energy->Fill(energy);
541  h1pt_energy_US->Fill(energy_US);
542  h1pt_energy_DS->Fill(energy_DS);
543 
544  if (layer == 1) pointlayer1 = 1;
545  if (layer == 2) pointlayer2 = 2;
546  if (layer == 3) pointlayer3 = 3;
547  if (layer == 4) pointlayer4 = 4;
548  // printf ("event=%d, pt jj=%d, module=%d, layer=%d, sector=%d, cell_id=%d, pointlayer1=%d, pointlayer2=%d, pointlayer3=%d, pointlayer4=%d\n",locEventNumber,jj,module,layer, sector,cell_id,pointlayer1,pointlayer2,pointlayer3,pointlayer4);
549 
550  // fill vector with module number if not yet in list
551  bool module_in_list = false;
552  for (unsigned int jjj=0; jjj<matchedShowers_modules.size(); jjj++) {
553  if (module == matchedShowers_modules[jjj]) module_in_list = true;
554  }
555  if (!module_in_list) matchedShowers_modules.push_back(module);
556  module_ndx = module;
557 
558  }
559 
560  // fill efficiency histograms, once per shower
561 
562  if (pointlayer2 + pointlayer3 + pointlayer4 > 0) {
563  h1eff_layertot->Fill(1);
564  h1eff_layer->Fill(pointlayer1);
565  int cellid = (module_ndx-1)*4 + 1;
566  h1eff_cellidtot->Fill(cellid);
567  if (pointlayer1>0) h1eff_cellid->Fill(cellid);
568  }
569 
570  if (pointlayer3 + pointlayer4 > 0) {
571  h1eff_layertot->Fill(2);
572  h1eff_layer->Fill(pointlayer2);
573  int cellid = (module_ndx-1)*4 + 2;
574  h1eff_cellidtot->Fill(cellid);
575  if (pointlayer2>0) h1eff_cellid->Fill(cellid);
576  }
577  if (pointlayer4 > 0) {
578  h1eff_layertot->Fill(3);
579  h1eff_layer->Fill(pointlayer3);
580  int cellid = (module_ndx-1)*4 + 3;
581  h1eff_cellidtot->Fill(cellid);
582  if (pointlayer3>0) h1eff_cellid->Fill(cellid);
583  }
584 
585  if (pointlayer1 + pointlayer2 + pointlayer3 > 3) {
586  h1eff_layertot->Fill(4);
587  h1eff_layer->Fill(pointlayer4);
588  int cellid = (module_ndx-1)*4 + 4;
589  h1eff_cellidtot->Fill(cellid);
590  if (pointlayer4>0) h1eff_cellid->Fill(cellid);
591  }
592 
593 
594 
595  int hitlayer1=0;
596  int hitlayer2=0;
597  int hitlayer3=0;
598  int hitlayer4=0;
599  Int_t hitmodule_ndx=0; // use one of the modules in shower as index
600  // loop over all points in BCAL
601  for (unsigned int jjj=0; jjj<bcalhits.size(); jjj++) {
602  int module = bcalhits[jjj]->module;
603 
604  // histogram all hits, whether or not in clustersssss
605 
606  if (bcalhits[jjj]->end == 0) h1E_US_layer1->Fill(bcalhits[jjj]->E);
607  if (bcalhits[jjj]->end == 1) h1E_DS_layer1->Fill(bcalhits[jjj]->E);
608 
609  // check to see if module is in cluster list
610  bool module_in_list = false;
611  for (unsigned int k=0; k<matchedShowers_modules.size(); k++) {
612  if (module == matchedShowers_modules[k]) module_in_list = true;
613  }
614  if (!module_in_list) continue;
615 
616  int layer = bcalhits[jjj]->layer;
617  // int sector = bcalhits[jjj]->sector;
618  // int end = bcalhits[jjj]->end;
619  // printf ("locPts event=%d, module=%d, layer=%d sector=%d end=%d\n",locEventNumber,module,layer,sector,end);
620 
621  if (layer == 1) hitlayer1 = 1;
622  if (layer == 2) hitlayer2 = 2;
623  if (layer == 3) hitlayer3 = 3;
624  if (layer == 4) hitlayer4 = 4;
625  hitmodule_ndx = module;
626 
627  }
628 
629  // fill efficiency histograms, once per shower
630 
631  if (hitlayer2 + hitlayer3 + hitlayer4 > 0) {
632  h1eff2_layertot->Fill(1);
633  h1eff2_layer->Fill(hitlayer1);;
634  int cellid = (hitmodule_ndx-1)*4 + 1;
635  h1eff2_cellidtot->Fill(cellid);
636  if (hitlayer1>0) h1eff2_cellid->Fill(cellid);
637  }
638 
639  if (hitlayer3 + hitlayer4 > 0) {
640  h1eff2_layertot->Fill(2);
641  h1eff2_layer->Fill(hitlayer2);;
642  int cellid = (hitmodule_ndx-1)*4 + 2;
643  h1eff2_cellidtot->Fill(cellid);
644  if (hitlayer2>0) h1eff2_cellid->Fill(cellid);
645  }
646  if (hitlayer4 > 0) {
647  h1eff2_layertot->Fill(3);
648  h1eff2_layer->Fill(hitlayer3);;
649  int cellid = (hitmodule_ndx-1)*4 + 3;
650  h1eff2_cellidtot->Fill(cellid);
651  if (hitlayer3>0) h1eff2_cellid->Fill(cellid);
652  }
653 
654  if (hitlayer1 + hitlayer2 + hitlayer3 > 3) {
655  h1eff2_layertot->Fill(4);
656  h1eff2_layer->Fill(hitlayer4);;
657  int cellid = (hitmodule_ndx-1)*4 + 4;
658  h1eff2_cellidtot->Fill(cellid);
659  if (hitlayer4>0) h1eff2_cellid->Fill(cellid);
660  }
661 
662 
663 
664  } // end loop over matched showers
665 
666  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
667 
668  /*
669  //Optional: Save event to output REST file. Use this to create skims.
670  dEventWriterREST->Write_RESTEvent(locEventLoop, "BCAL_Shower"); //string is part of output file name
671  */
672 
673  return NOERROR;
674 }
675 
676 //------------------
677 // erun
678 //------------------
679 jerror_t JEventProcessor_BCAL_Eff::erun(void)
680 {
681  // This is called whenever the run number changes, before it is
682  // changed to give you a chance to clean up before processing
683  // events from the next run number.
684 
685  // FILL HISTOGRAMS
686  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
687  //japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
688  // h1eff_eff->Divide(h1eff_layer,h1eff_layertot);
689  //japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
690 
691  return NOERROR;
692 }
693 
694 //------------------
695 // fini
696 //------------------
697 jerror_t JEventProcessor_BCAL_Eff::fini(void)
698 {
699  // Called before program exit after event processing is finished.
700 
701  h1eff_eff->Divide(h1eff_layer,h1eff_layertot,1,1,"B");
702  h1eff2_eff2->Divide(h1eff2_layer,h1eff2_layertot,1,1,"B");
703 
704  h1eff_cellideff->Divide(h1eff_cellid,h1eff_cellidtot,1,1,"B");
705  h1eff2_cellideff2->Divide(h1eff2_cellid,h1eff2_cellidtot,1,1,"B");
706 
707  return NOERROR;
708 }
709 
h1eff2_cellideff2
TH1F * h1eff2_cellideff2
Definition: bcal_hist_eff.C:40
DTrackFitter::ExtrapolateToRadius
bool ExtrapolateToRadius(double R, const vector< Extrapolation_t > &extraps, DVector3 &pos, DVector3 &mom, double &t, double &s) const
Definition: DTrackFitter.cc:473
DBCALShower.h
JEventProcessor_BCAL_Eff::brun
jerror_t brun(jana::JEventLoop *locEventLoop, int locRunNumber)
Called every time a new run number is detected.
Definition: JEventProcessor_BCAL_Eff.cc:295
JEventProcessor_BCAL_Eff::evnt
jerror_t evnt(jana::JEventLoop *locEventLoop, uint64_t locEventNumber)
Called every event.
Definition: JEventProcessor_BCAL_Eff.cc:335
DTrackFitter
The DTrackFitter class is a base class for different charged track fitting algorithms. It does not actually fit the track itself, but provides the interface and some common support features most algorthims will need to implement.
Definition: DTrackFitter.h:61
JEventProcessor_BCAL_Eff::init
jerror_t init(void)
Called once at program start.
Definition: JEventProcessor_BCAL_Eff.cc:108
layer
Int_t layer
Definition: Read_bcal_hadronic_eff2.C:50
DVector3
TVector3 DVector3
Definition: DVector3.h:14
h1eff2_layer
TH1F * h1eff2_layer
Definition: bcal_hist_eff.C:44
h2_pmagvsenergy
static TH2I * h2_pmagvsenergy
Definition: JEventProcessor_BCAL_Eff.cc:57
x
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
y
#define y
h1trk_FOM
static TH1I * h1trk_FOM
Definition: JEventProcessor_BCAL_Eff.cc:49
h1trk_pmag
static TH1I * h1trk_pmag
Definition: JEventProcessor_BCAL_Eff.cc:50
DApplication.h
DTrigger::Get_L1FrontPanelTriggerBits
uint32_t Get_L1FrontPanelTriggerBits(void) const
Definition: libraries/TRIGGER/DTrigger.h:14
h1pt_layer
static TH1I * h1pt_layer
Definition: JEventProcessor_BCAL_Eff.cc:64
JEventProcessor_BCAL_Eff::fini
jerror_t fini(void)
Called after last event of last event source has been processed.
Definition: JEventProcessor_BCAL_Eff.cc:697
h1eff2_eff2
TH1F * h1eff2_eff2
Definition: bcal_hist_eff.C:39
h1trk_pz
static TH1I * h1trk_pz
Definition: JEventProcessor_BCAL_Eff.cc:53
DBCALTruthShower.h
SYS_BCAL
Definition: GlueX.h:19
japp
JApplication * japp
h1trk_energy
static TH1I * h1trk_energy
Definition: JEventProcessor_BCAL_Eff.cc:54
h1eff_cellidtot
TH1F * h1eff_cellidtot
Definition: bcal_hist_eff.C:47
h1eff_layertot
TH1F * h1eff_layertot
Definition: bcal_hist_eff.C:43
event_count
int event_count
Definition: JEventProcessor_BCAL_point_calib.cc:97
DBCALPoint.h
h1pt_energy_DS
static TH1I * h1pt_energy_DS
Definition: JEventProcessor_BCAL_Eff.cc:69
InitPlugin
InitPlugin_t InitPlugin
h1pt_energy_US
static TH1I * h1pt_energy_US
Definition: JEventProcessor_BCAL_Eff.cc:68
h1eff2_layertot
TH1F * h1eff2_layertot
Definition: bcal_hist_eff.C:45
_DBG_
#define _DBG_
Definition: HDEVIO.h:12
h1eff_layer
TH1F * h1eff_layer
Definition: bcal_hist_eff.C:42
DBCALHit.h
h1eff_cellideff
TH1F * h1eff_cellideff
Definition: bcal_hist_eff.C:38
h1E_US_layer1
static TH1I * h1E_US_layer1
Definition: JEventProcessor_BCAL_Eff.cc:87
h1pt_module
static TH1I * h1pt_module
Definition: JEventProcessor_BCAL_Eff.cc:63
h1pt_energy
static TH1I * h1pt_energy
Definition: JEventProcessor_BCAL_Eff.cc:67
py
TH1D * py[NCHANNELS]
Definition: ST_Propagation_Time/macros/Resolution.C:72
JEventProcessor_BCAL_Eff::Run_Number
int Run_Number
Definition: JEventProcessor_BCAL_Eff.h:34
h1pt_sector
static TH1I * h1pt_sector
Definition: JEventProcessor_BCAL_Eff.cc:65
JEventProcessor_BCAL_Eff.h
h2_Evsenergy
static TH2I * h2_Evsenergy
Definition: JEventProcessor_BCAL_Eff.cc:56
DFCALShower.h
h1_Num_matched_showers
static TH1I * h1_Num_matched_showers
Definition: JEventProcessor_BCAL_Eff.cc:35
h1trk_py
static TH1I * h1trk_py
Definition: JEventProcessor_BCAL_Eff.cc:52
h1eff_eff
TH1F * h1eff_eff
Definition: bcal_hist_eff.C:37
h1eff_cellid
TH1F * h1eff_cellid
Definition: bcal_hist_eff.C:46
DAnalysisUtilities.h
h1trk_px
static TH1I * h1trk_px
Definition: JEventProcessor_BCAL_Eff.cc:51
fitter
const DTrackFitter * fitter
Definition: DEventProcessor_trk_profile.cc:14
h1eff2_cellid
TH1F * h1eff2_cellid
Definition: bcal_hist_eff.C:48
h1eff2_cellidtot
TH1F * h1eff2_cellidtot
Definition: bcal_hist_eff.C:49
JEventProcessor_BCAL_Eff::erun
jerror_t erun(void)
Called every time run number changes, provided brun has been called.
Definition: JEventProcessor_BCAL_Eff.cc:679
h1pt_cell_id
static TH1I * h1pt_cell_id
Definition: JEventProcessor_BCAL_Eff.cc:66
JEventProcessor_BCAL_Eff
Definition: JEventProcessor_BCAL_Eff.h:25
DMCThrown.h
h1trk_Num_matched_tracks
static TH1I * h1trk_Num_matched_tracks
Definition: JEventProcessor_BCAL_Eff.cc:48
h1pt_Num_points
static TH1I * h1pt_Num_points
Definition: JEventProcessor_BCAL_Eff.cc:62
h1E_DS_layer1
static TH1I * h1E_DS_layer1
Definition: JEventProcessor_BCAL_Eff.cc:88
main
int main(int argc, char *argv[])
Definition: gendoc.cc:6
h2E_USvsDS_layer1
static TH2I * h2E_USvsDS_layer1
Definition: JEventProcessor_BCAL_Eff.cc:89