Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventProcessor_photoneff_hists.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: DEventProcessor_photoneff_hists.cc
4 // Created: Thu Feb 12 09:43:13 EST 2009
5 // Creator: davidl (on Darwin harriet.jlab.org 9.6.0 i386)
6 //
7 
8 #include <iostream>
9 #include <cmath>
10 using namespace std;
11 
12 #include <TThread.h>
13 
14 #include <TROOT.h>
15 
17 
18 #include <JANA/JApplication.h>
19 #include <JANA/JEventLoop.h>
20 
21 #include <DANA/DApplication.h>
22 #include <DVector2.h>
23 #include <particleType.h>
24 #include <FCAL/DFCALHit.h>
25 #include <BCAL/DBCALMCResponse.h>
26 //#include <FCAL/DFCALMCResponse.h>
27 
28 // Routine used to create our DEventProcessor
29 extern "C"{
30 void InitPlugin(JApplication *app){
31  InitJANAPlugin(app);
32  app->AddProcessor(new DEventProcessor_photoneff_hists());
33 }
34 } // "C"
35 
36 
37 //------------------
38 // DEventProcessor_photoneff_hists
39 //------------------
41 {
42  phtn_ptr = &phtn;
43 
44 
45  pthread_mutex_init(&mutex, NULL);
46 
47 }
48 
49 //------------------
50 // ~DEventProcessor_photoneff_hists
51 //------------------
53 {
54 
55 }
56 
57 //------------------
58 // init
59 //------------------
61 {
62  // Create TRACKING directory
63  TDirectory *dir = (TDirectory*)gROOT->FindObject("CALORIMETRY");
64  if(!dir)dir = new TDirectoryFile("CALORIMETRY","CALORIMETRY");
65  dir->cd();
66 
67  // Create Trees
68  phtneff = new TTree("phtneff","Photon Reconstruction Efficiency");
69  phtneff->Branch("E","photon",&phtn_ptr);
70 
71  dir->cd("../");
72 
73  // Create ROOT file
74  RootFile = new TFile("photoneff.root","RECREATE");
75 
76  JParameterManager *parms = app->GetJParameterManager();
77 
78  DEBUG = 1;
79 
80  parms->SetDefaultParameter("TRKEFF:DEBUG", DEBUG);
81 
82  return NOERROR;
83 }
84 
85 //------------------
86 // brun
87 //------------------
88 jerror_t DEventProcessor_photoneff_hists::brun(JEventLoop *loop, int32_t runnumber)
89 {
90 
91  return NOERROR;
92 }
93 
94 //------------------
95 // erun
96 //------------------
98 {
99  return NOERROR;
100 }
101 
102 //------------------
103 // fini
104 //------------------
106 {
107  return NOERROR;
108 }
109 
110 //------------------
111 // evnt
112 //------------------
113 jerror_t DEventProcessor_photoneff_hists::evnt(JEventLoop *loop, uint64_t eventnumber)
114 {
115  vector<const DFCALHit*> fcalhits;
116  vector<const DBCALMCResponse*> bcalhits;
117  vector<const DPhoton*> photons;
118  vector<const DMCThrown*> mcthrowns;
119  vector<const DMCTrajectoryPoint*> mctrajpoints;
120 
121  loop->Get(fcalhits);
122  loop->Get(bcalhits);
123  loop->Get(photons);
124  loop->Get(mcthrowns);
125  loop->Get(mctrajpoints);
126 
127  // Get hit list for all throwns
128  for(unsigned int i=0; i<mcthrowns.size(); i++){
129  const DMCThrown *mcthrown = mcthrowns[i];
130 
131  // if this is a charged track, then skip it
132  if(fabs(mcthrowns[i]->charge())!=0.0)continue;
133 
134  // Momentum of thrown particle
135  DVector3 pthrown = mcthrown->momentum();
136 
137  bool locIsRecon = isReconstructable(mcthrown, loop);
138 
139  // Although we are only filling objects local to this plugin, TTree::Fill() periodically writes to file: Global ROOT lock
140  japp->RootWriteLock(); //ACQUIRE ROOT LOCK
141 
142  phtn.pthrown = pthrown;
143  // Initialize with the "not found" values
144  phtn.pfit.SetXYZ(0,0,0);
145  phtn.chisq=1.0E20;
146  phtn.Ndof=-1;
147  phtn.delta_E_over_E=1.0E20;
148  phtn.delta_theta=1.0E20;
149  phtn.delta_phi=1.0E20;
150  phtn.isreconstructable = locIsRecon;
151  phtn.Nbcal = 0;
152  phtn.Nfcal = 0;
153  phtn.event = eventnumber;
154 
155  double fom_best = 1.0E8;
156 
157  // Loop over found/fit photons
158  for(unsigned int j=0; j<photons.size(); j++){
159  const DPhoton *photon = photons[j];
160 
161  // Copy momentum vector to convenient local variables
162  DVector3 pfit = photon->momentum();
163 
164  // Calculate residuals from momentum parameters from DParticle
165  double delta_E_over_E = (pfit.Mag()-pthrown.Mag())/pthrown.Mag();
166  double delta_theta = (pfit.Theta() - pthrown.Theta())*1000.0;
167  double delta_phi = (pfit.Phi() - pthrown.Phi())*1000.0;
168 
169  // Formulate a figure of merit to decide if this fit photon is closer to
170  // the thrown photon than the best one we found so far.
171  double fom = fabs(delta_E_over_E);
172  if(fom<fom_best){
173  fom = fom_best;
174 
175  phtn.pfit = pfit;
176  //phtn.chisq = photon->chisq;
177  //phtn.Ndof = photon->Ndof;
178  phtn.delta_E_over_E = delta_E_over_E;
179  phtn.delta_theta = delta_theta;
180  phtn.delta_phi = delta_phi;
181  phtn.Nbcal = bcalhits.size(); // Need hits actually used in DPhoton!
182  phtn.Nfcal = fcalhits.size(); // Need hits actually used in DPhoton!
183 
184  }
185  }
186 
187  phtneff->Fill();
188 
189  japp->RootUnLock(); //RELEASE ROOT LOCK
190  }
191 
192  return NOERROR;
193 }
194 
195 //------------------
196 // isReconstructable
197 //------------------
198 bool DEventProcessor_photoneff_hists::isReconstructable(const DMCThrown *mcthrown, JEventLoop *loop)
199 {
200  /// For photon reconstruction, we just check if the total energy
201  /// in ether the BCAL or FCAL is greater than 50% of the thrown
202  /// value to determine if it is reconstructible.
203  ///
204  /// This is not entirely accurate since a photon could convert
205  /// upstream and spray that much energy into a calorimeter
206  /// and the photon would then not be reconstructible.
207  ///
208  /// On the other hand, if you are really looking for detection
209  /// efficiency, then you'll want that to be included in the
210  /// calculation so maybe this is a good thing ....?
211 
212  vector<const DBCALMCResponse*> bcalhits;
213  vector<const DFCALHit*> fcalhits;
214 
215  loop->Get(bcalhits);
216  loop->Get(fcalhits);
217 
218  double Ebcal = 0.0;
219  for(unsigned int i=0; i<bcalhits.size(); i++)Ebcal += (bcalhits[i])->E;
220  if(Ebcal>0.5*mcthrown->energy())return true;
221 
222  double Efcal = 0.0;
223  for(unsigned int i=0; i<fcalhits.size(); i++)Efcal += (fcalhits[i])->E;
224  if(Efcal>0.5*mcthrown->energy())return true;
225 
226  return false;
227 }
228 
Definition: photon.h:15
TVector3 DVector3
Definition: DVector3.h:14
double energy(void) const
jerror_t brun(JEventLoop *loop, int32_t runnumber)
JApplication * japp
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
bool isReconstructable(const DMCThrown *mcthrown, JEventLoop *loop)
InitPlugin_t InitPlugin
jerror_t init(void)
Invoked via DEventProcessor virtual method.
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
const DVector3 & momentum(void) const
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
TDirectory * dir
Definition: bcal_hist_eff.C:31
#define DEBUG
Definition: tw_corr.C:41