Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DEventProcessor_acceptance_hists.cc
Go to the documentation of this file.
1 // $Id: DEventProcessor_acceptance_hists.cc 1816 2006-06-06 14:38:18Z davidl $
2 //
3 // File: DEventProcessor_acceptance_hists.cc
4 // Created: Sun Apr 24 06:45:21 EDT 2005
5 // Creator: davidl (on Darwin Harriet.local 7.8.0 powerpc)
6 //
7 
8 #include <iostream>
9 #include <map>
10 #include <cmath>
11 using namespace std;
12 
13 #include <TThread.h>
14 
15 #include <JANA/JEventLoop.h>
16 
18 
19 #include <DANA/DApplication.h>
20 #include <TRACKING/DMCThrown.h>
21 #include <TRACKING/DMCTrackHit.h>
22 #include <FDC/DFDCHit.h>
23 #include <FDC/DFDCGeometry.h>
24 #include <CDC/DCDCTrackHit.h>
25 
26 #define MIN_CDC_HITS 8
27 #define MIN_FDC_HITS 8
28 #define MIN_CDC_FDC_HITS 8
29 #define MIN_BCAL_HITS 4
30 #define MIN_FCAL_HITS 4
31 #define MIN_TOF_HITS 1
32 #define MIN_UPV_HITS 4
33 
34 // The executable should define the ROOTfile global variable. It will
35 // be automatically linked when dlopen is called.
36 //extern TFile *ROOTfile;
37 
38 // Routine used to create our DEventProcessor
39 extern "C"{
40 void InitPlugin(JApplication *app){
41  InitJANAPlugin(app);
42  app->AddProcessor(new DEventProcessor_acceptance_hists());
43 }
44 } // "C"
45 
46 
47 //------------------
48 // DEventProcessor_acceptance_hists
49 //------------------
51 {
52 }
53 
54 //------------------
55 // ~DEventProcessor_acceptance_hists
56 //------------------
58 {
59 }
60 
61 //------------------
62 // init
63 //------------------
65 {
66  // open ROOT file (if needed)
67  //if(ROOTfile != NULL) ROOTfile->cd();
68 
69  // Create ACCEPTANCE directory
70  TDirectory *dir = new TDirectoryFile("ACCEPTANCE","ACCEPTANCE");
71  dir->cd();
72 
73  // Create histograms
74  int N_p_bins = 100;
75  float p_lo = 0.0;
76  float p_hi = 12.0;
77  int N_theta_bins = 400;
78  float theta_lo = 0.0;
79  float theta_hi = M_PI*57.3;
80  CDC = new TH2F("CDC","CDC acceptance",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
81  FDC = new TH2F("FDC","FDC acceptance",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
82  CDC_FDC = new TH2F("CDC_FDC","CDC_FDC acceptance",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
83  BCAL = new TH2F("BCAL","BCAL acceptance",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
84  FCAL = new TH2F("FCAL","FCAL acceptance",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
85  TOF = new TH2F("TOF","TOF acceptance",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
86  thrown_charged = new TH2F("thrown_charged","thrown charged particles",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
87  thrown_photon = new TH2F("thrown_photon","thrown photons",N_p_bins, p_lo, p_hi, N_theta_bins, theta_lo, theta_hi);
88 
89  FDC_anode_hits_per_event = new TH1D("FDC_anode_hits_per_event","FDC anode hits/event", 201, -0.5, 200.5);
90  FDC_anode_hits_per_layer = new TH1D("FDC_anode_hits_per_layer","FDC anode hits/layer", 24, 0.5, 24.5);
91  FDC_anode_hits_per_wire = new TH1D("FDC_anode_hits_per_wire","FDC anode hits/wire", 96, 0.5, 96.5);
92 
93  CDC_nhits_vs_pthrown = new TH1D("CDC_nhits_vs_pthrown","Number of CDC hits per event vs. thrown momentum", 40, 0.0, 9.0);
94  FDC_nhits_vs_pthrown = new TH1D("FDC_nhits_vs_pthrown","Number of FDC anode hits per event vs. thrown momentum", 40, 0.0, 9.0);
95  pthrown = new TH1D("pthrown","thrown momentum", 40, 0.0, 9.0);
96 
97  // Go back up to the parent directory
98  dir->cd("../");
99 
100  return NOERROR;
101 }
102 
103 //------------------
104 // evnt
105 //------------------
106 jerror_t DEventProcessor_acceptance_hists::evnt(JEventLoop *loop, uint64_t eventnumber)
107 {
108  vector<const DMCThrown*> mcthrowns;
109  vector<const DMCTrackHit*> mctrackhits;
110  loop->Get(mcthrowns);
111  loop->Get(mctrackhits);
112 
113  // Count CDC hits
114  double Ncdc_anode = 0.0;
115  vector<const DCDCTrackHit*> cdctrackhits;
116  loop->Get(cdctrackhits);
117  Ncdc_anode = (double)cdctrackhits.size();
118 
119  // Count FDC anode hits
120  double Nfdc_anode = 0.0;
121  vector<const DFDCHit*> fdchits;
122  loop->Get(fdchits);
123 
124  // FILL HISTOGRAMS
125  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
126  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
127  {
128  for(unsigned int i=0; i<fdchits.size(); i++){
129  if(fdchits[i]->type==0){
130  Nfdc_anode+=1.0;
131 
132  FDC_anode_hits_per_layer->Fill(DFDCGeometry::gLayer(fdchits[i]));
133  FDC_anode_hits_per_wire->Fill(fdchits[i]->element);
134  }
135  }
136  FDC_anode_hits_per_event->Fill(Nfdc_anode);
137 
138  // Loop over thrown tracks
139  for(unsigned int i=0;i<mcthrowns.size();i++){
140  const DMCThrown *mcthrown = mcthrowns[i];
141 
142  if(mcthrown->charge() != 0.0){
143  thrown_charged->Fill(mcthrown->momentum().Mag(), mcthrown->momentum().Theta()*57.3);
144  }else if(mcthrown->type == 1){
145  thrown_photon->Fill(mcthrown->momentum().Mag(), mcthrown->momentum().Theta()*57.3);
146  }
147  }
148 
149  // Loop over track hits
150  map<int,int> cdchits;
151  map<int,int> fdchitmap;
152  map<int,int> bcalhits;
153  for(unsigned int i=0;i<mctrackhits.size();i++){
154  const DMCTrackHit *mctrackhit = mctrackhits[i];
155 
156  switch(mctrackhit->system){
157  case SYS_CDC:
158  if(mctrackhit->primary)cdchits[mctrackhit->track]++;
159  break;
160  case SYS_FDC:
161  if(mctrackhit->primary)fdchitmap[mctrackhit->track]++;
162  break;
163  case SYS_BCAL:
164  bcalhits[mctrackhit->track]++;
165  break;
166  default:
167  break;
168  }
169  }
170 
171  // Simple 1-D histos for CDC and FDC as a function of thrown momentum
172  if(mcthrowns.size()==1){
173  double p = mcthrowns[0]->momentum().Mag();
174  CDC_nhits_vs_pthrown->Fill(p, Ncdc_anode);
175  FDC_nhits_vs_pthrown->Fill(p, Nfdc_anode);
176  pthrown->Fill(p);
177  }
178 
179 
180  // NOTE: In the sections that follow we have to assume that
181  // the track number of the hit corresponds to the position of
182  // the DMCThrown object in the list of DMCThrown objects.
183  // This should be a good assumption, but I don't know that
184  // it is (or always will be) guaranteed.
185 
186  // Loop over tracks in the CDC
187  map<int,int>::iterator iter;
188  for(iter=cdchits.begin(); iter!=cdchits.end(); iter++){
189 
190  // Find thrown parameters for this track (if any)
191  if(iter->first<=0 || iter->first>(int)mcthrowns.size())continue;
192  const DMCThrown *mcthrown = mcthrowns[iter->first-1];
193 
194  if(iter->second >= MIN_CDC_HITS)CDC->Fill(mcthrown->momentum().Mag(), mcthrown->momentum().Theta()*57.3);
195  }
196 
197  // Loop over tracks in the FDC
198  for(iter=fdchitmap.begin(); iter!=fdchitmap.end(); iter++){
199 
200  // Find thrown parameters for this track (if any)
201  if(iter->first<=0 || iter->first>(int)mcthrowns.size())continue;
202  const DMCThrown *mcthrown = mcthrowns[iter->first-1];
203 
204  if(iter->second >= MIN_FDC_HITS)FDC->Fill(mcthrown->momentum().Mag(), mcthrown->momentum().Theta()*57.3);
205 
206  // Fill CDC + FDC histo
207  if(cdchits.find(iter->first) != cdchits.end()){
208  int cdc_fdc_hits = cdchits.find(iter->first)->second + iter->second;
209  if(cdc_fdc_hits >= MIN_CDC_FDC_HITS)
210  CDC_FDC->Fill(mcthrown->momentum().Mag(), mcthrown->momentum().Theta()*57.3);
211  }
212  }
213  }
214  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
215 
216  return NOERROR;
217 }
218 
219 //------------------
220 // erun
221 //------------------
223 {
224  // Since we are modifying histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
225  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
226  {
227  CDC->Divide(thrown_charged);
228  FDC->Divide(thrown_charged);
229  CDC_FDC->Divide(thrown_charged);
230  }
231  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
232 
233  return NOERROR;
234 }
235 
236 //------------------
237 // fini
238 //------------------
240 {
241 
242  return NOERROR;
243 }
jerror_t fini(void)
Invoked via DEventProcessor virtual method.
Definition: GlueX.h:17
Definition: GlueX.h:19
Definition: CDC.h:14
JApplication * japp
DetectorSystem_t system
particle type
Definition: DMCTrackHit.h:25
#define MIN_CDC_FDC_HITS
jerror_t init(void)
Invoked via DEventProcessor virtual method.
static int gLayer(const DFDCHit *h)
DFDCGeometry::gLayer(): Compute the global layer (detection layer 1-24) number based on module and la...
Definition: DFDCGeometry.h:59
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via DEventProcessor virtual method.
InitPlugin_t InitPlugin
jerror_t erun(void)
Invoked via DEventProcessor virtual method.
Definition: GlueX.h:18
double charge(void) const
Definition: FCAL.h:14
int primary
primary track=1 not primary track=0
Definition: DMCTrackHit.h:23
const DVector3 & momentum(void) const
TDirectory * dir
Definition: bcal_hist_eff.C:31
int type
GEANT particle ID.
Definition: DMCThrown.h:20
Definition: TOF.h:14
int track
Track number.
Definition: DMCTrackHit.h:21