Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_PS_online.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_PS_online.cc
4 // Created: Fri Nov 9 11:58:09 EST 2012
5 // Creator: wolin (on Linux stan.jlab.org 2.6.32-279.11.1.el6.x86_64 x86_64)
6 
7 
8 #include <stdint.h>
9 #include <vector>
10 
12 #include <JANA/JApplication.h>
13 
14 using namespace std;
15 using namespace jana;
16 
20 
21 #include <TDirectory.h>
22 #include <TH1.h>
23 #include <TH2.h>
24 #include <TProfile.h>
25 
26 // number of PS columns (tiles) per arm
28 const int Narms = DPSGeometry::NUM_ARMS; // 2
29 
30 const int NmultBins = 10; //number of bins for multiplicity histograms
31 
32 double counts_cut = 500.0;
33 // root hist pointers
34 static TH1I *ps_num_events;
35 static TH1I *hHit_NHits;
36 static TH1I *hHit_Arm;
37 static TH2I *hHit_NHitsVsArm;
38 static TH1I *hHit_Occupancy[Narms];
39 static TH1I *hHit_Energy[Narms];
40 static TH2I *hHit_EnergyVsColumn[Narms];
41 static TH1I *hHit_Integral[Narms];
43 //static TH1I *hHit_Npix[Narms];
44 static TH1I *hHit_Time[Narms];
45 static TH2I *hHit_TimeVsColumn[Narms];
46 
47 static TH1I *hDigiHit_NHits;
48 static TH1I *hDigiHit_Arm;
49 static TH2I *hDigiHit_NHitsVsArm;
51 static TH1I *hDigiHit_Pedestal[Narms];
52 static TProfile *hDigiHit_PedestalVsColumn[Narms];
54 static TH1I *hDigiHit_Occupancy[Narms];
55 static TH1I *hDigiHit_RawPeak[Narms];
63 static TH1I *hDigiHit_PulseTime[Narms];
64 static TH1I *hDigiHit_Time[Narms];
67 
68 static TH1I *hDigiHit_NHits_cut;
69 static TH1I *hDigiHit_Arm_cut;
72 static TH1I *hDigiHit_Time_cut[Narms];
75 //----------------------------------------------------------------------------------
76 
77 // Routine used to create our JEventProcessor
78 extern "C"{
79  void InitPlugin(JApplication *app){
80  InitJANAPlugin(app);
81  app->AddProcessor(new JEventProcessor_PS_online());
82  }
83 }
84 
85 
86 //----------------------------------------------------------------------------------
87 
88 
90 }
91 
92 
93 //----------------------------------------------------------------------------------
94 
95 
97 }
98 
99 
100 //----------------------------------------------------------------------------------
101 
103 
104  // create root folder for ps and cd to it, store main dir
105  TDirectory *mainDir = gDirectory;
106  TDirectory *psDir = gDirectory->mkdir("PS");
107  psDir->cd();
108 
109  // book hists
110  ps_num_events = new TH1I("ps_num_events","PS Number of events",1,0.5,1.5);
111  // fADC250 hit-level hists (after calibration)
112  TDirectory *hitDir = gDirectory->mkdir("Hit"); hitDir->cd();
113  hHit_NHits = new TH1I("Hit_NHits","PS hit multiplicity;hits;events",NmultBins,0.5,0.5+NmultBins);
114  hHit_Arm = new TH1I("Hit_Arm","PS arm;arm;hits",Narms,-0.5,-0.5+Narms);
115  hHit_NHitsVsArm = new TH2I("Hit_NHitsVsArm","PS hit multiplicity vs. arm;arm;hits",Narms,-0.5,-0.5+Narms,NmultBins,0.5,0.5+NmultBins);
116  TString arm_str[] = {"Left","Right"};
117  for (int i = 0; i < Narms; i++) {
118  gDirectory->mkdir(arm_str[i]+"Arm")->cd();
119  TString strN = "_" + arm_str[i] + "Arm";
120  TString strT = ", " + arm_str[i] + " arm";
121  hHit_Occupancy[i] = new TH1I("Hit_Occupancy"+strN,"PS occupancy"+strT+";column (tile);hits / column",Ncols,0.5,0.5+Ncols);
122  hHit_Energy[i] = new TH1I("Hit_Energy"+strN,"PS energy"+strT+";energy [GeV];hits / column",120,0.5,6.5);
123  hHit_EnergyVsColumn[i] = new TH2I("Hit_EnergyVsColumn"+strN,"PS energy vs. column"+strT+";column (tile);energy [GeV]",Ncols,0.5,0.5+Ncols,120,0.5,6.5);
124  hHit_Integral[i] = new TH1I("Hit_Integral"+strN,"PS fADC pulse integral"+strT+";pulse integral;hits",1000,0.0,30000.0);
125  hHit_IntegralVsColumn[i] = new TH2I("Hit_IntegralVsColumn"+strN,"PS fADC pulse integral vs. column"+strT+";column (tile);pulse integral",Ncols,0.5,0.5+Ncols,1000,0.0,30000.0);
126  //hHit_Npix[i] = new TH1I("Hit_Npix"+strN,"PS fADC number of pixels"+strT+";pixels;hits",200,0,100);
127  hHit_Time[i] = new TH1I("Hit_Time"+strN,"PS fADC time"+strT+";time [ns];hits / 400 ps",1000,-200.0,200.0);
128  hHit_TimeVsColumn[i] = new TH2I("Hit_TimeVsColumn"+strN,"PS fADC time vs. column"+strT+";column (tile);time [ns]",Ncols,0.5,0.5+Ncols,1000,-200.0,200.0);
129  hitDir->cd();
130  }
131  // fADC250 digihit-level hists
132  psDir->cd();
133  TDirectory *digihitDir = gDirectory->mkdir("DigiHit"); digihitDir->cd();
134  hDigiHit_NHits = new TH1I("DigiHit_NHits","PS fADC hit multiplicity;raw hits;events",NmultBins,0.5,0.5+NmultBins);
135  hDigiHit_Arm = new TH1I("DigiHit_Arm","PS arm;arm;raw hits",Narms,-0.5,-0.5+Narms);
136  hDigiHit_NHitsVsArm = new TH2I("DigiHit_NHitsVsArm","PS hit multiplicity vs. arm;arm;raw hits",Narms,-0.5,-0.5+Narms,NmultBins,0.5,0.5+NmultBins);
137  hDigiHit_NHits_cut = new TH1I("DigiHit_NHits_cut","PS fADC hit multiplicity (> 500 ADC integral counts);raw hits;events",NmultBins,0.5,0.5+NmultBins);
138  hDigiHit_Arm_cut = new TH1I("DigiHit_Arm_cut","PS arm (> 500 ADC integral counts);arm;raw hits",Narms,-0.5,-0.5+Narms);
139  hDigiHit_NHitsVsArm_cut = new TH2I("DigiHit_NHitsVsArm_cut","PS hit multiplicity vs. arm (> 500 ADC integral counts);arm;raw hits",Narms,-0.5,-0.5+Narms,NmultBins,0.5,0.5+NmultBins);
140  for (int i = 0; i < Narms; i++) {
141  gDirectory->mkdir(arm_str[i]+"Arm")->cd();
142  TString strN = "_" + arm_str[i] + "Arm";
143  TString strT = ", " + arm_str[i] + " arm";
144  hDigiHit_NSamplesPedestal[i] = new TH1I("DigiHit_NSamplesPedestal"+strN,"PS fADC pedestal samples"+strT+";pedestal samples;raw hits",50,-0.5,49.5);
145  hDigiHit_Pedestal[i] = new TH1I("DigiHit_Pedestal"+strN,"PS fADC pedestals"+strT+";pedestal [fADC counts];raw hits",200,0.0,200.0);
146  hDigiHit_PedestalVsColumn[i] = new TProfile("DigiHit_PedestalVsColumn"+strN,"PS pedestal vs. column"+strT+";column (tile);average pedestal [fADC counts]",Ncols,0.5,0.5+Ncols,"s");
147  hDigiHit_QualityFactor[i] = new TH1I("DigiHit_QualityFactor"+strN,"PS fADC quality factor"+strT+";quality factor;raw hits",4,-0.5,3.5);
148  hDigiHit_Occupancy[i] = new TH1I("DigiHit_Occupancy"+strN,"PS fADC hit occupancy"+strT+";column (tile);raw hits / column",Ncols,0.5,0.5+Ncols);
149  hDigiHit_RawPeak[i] = new TH1I("DigiHit_RawPeak"+strN,"PS fADC pulse peak (raw)"+strT+";pulse peak (raw);raw hits",410,0.0,4100.0);
150  hDigiHit_RawPeakVsColumn[i] = new TH2I("DigiHit_RawPeakVsColumn"+strN,"PS fADC pulse peak (raw) vs. column"+strT+";column (tile);pulse peak (raw)",Ncols,0.5,0.5+Ncols,410,0.0,4100.0);
151  hDigiHit_RawIntegral[i] = new TH1I("DigiHit_RawIntegral"+strN,"PS fADC pulse integral (raw)"+strT+";pulse integral (raw);raw hits",1000,0.0,30000.0);
152  hDigiHit_RawIntegralVsColumn[i] = new TH2I("DigiHit_RawIntegralVsColumn"+strN,"PS fADC pulse integral (raw) vs. column"+strT+";column (tile);pulse integral (raw)",Ncols,0.5,0.5+Ncols,1000,0.0,30000.0);
153  hDigiHit_NSamplesIntegral[i] = new TH1I("DigiHit_NSamplesIntegral"+strN,"PS fADC integral samples"+strT+";integral samples;raw hits",60,-0.5,59.5);
154  hDigiHit_PeakVsColumn[i] = new TH2I("DigiHit_PeakVsColumn"+strN,"PS fADC pulse peak vs. column"+strT+";column (tile);pulse peak",Ncols,0.5,0.5+Ncols,410,0.0,4100.0);
155  hDigiHit_IntegralVsPeak[i] = new TH2I("DigiHit_IntegralVsPeak"+strN,"PS fADC pulse integral vs. peak"+strT+";pulse peak;pulse integral",410,0.0,4100.0,1000,0.0,30000.0);
156  hDigiHit_IntegralVsColumn[i] = new TH2I("DigiHit_IntegralVsColumn"+strN,"PS fADC pulse integral vs. column"+strT+";column (tile);pulse integral",Ncols,0.5,0.5+Ncols,1000,0.0,30000.0);
157  hDigiHit_PulseTime[i] = new TH1I("DigiHit_PulseTime"+strN,"PS fADC pulse time"+strT+";pulse time [62.5 ps];raw hits",2000,0.0,6500.0);
158  hDigiHit_Time[i] = new TH1I("DigiHit_Time"+strN,"PS fADC pulse time"+strT+";pulse time [ns];raw hits / 400 ps",1000,0.0,400.0);
159  hDigiHit_TimeVsColumn[i] = new TH2I("DigiHit_TimeVsColumn"+strN,"PS fADC pulse time vs. column"+strT+";column (tile);pulse time [ns]",Ncols,0.5,0.5+Ncols,1000,0.0,400.0);
160  hDigiHit_TimeVsIntegral[i] = new TH2I("DigiHit_TimeVsIntegral"+strN,"PS fADC pulse time vs. integral"+strT+";pulse integral;pulse time [ns]",1000,0.0,30000.0,1000,0.0,400.0);
161  // digihit-level hists after cut on ADC integral counts
162  hDigiHit_Occupancy_cut[i] = new TH1I("DigiHit_Occupancy_cut"+strN,"PS fADC hit occupancy (> 500 ADC integral counts)"+strT+";column (tile);raw hits / column",Ncols,0.5,0.5+Ncols);
163  hDigiHit_Time_cut[i] = new TH1I("DigiHit_Time_cut"+strN,"PS fADC pulse time (> 500 ADC integral counts)"+strT+";pulse time [ns];raw hits / 400 ps",1000,0.0,400.0);
164  hDigiHit_TimeVsColumn_cut[i] = new TH2I("DigiHit_TimeVsColumn_cut"+strN,"PS fADC pulse time vs. column (> 500 ADC integral counts)"+strT+";column (tile);pulse time [ns]",Ncols,0.5,0.5+Ncols,1000,0.0,400.0);
165  hDigiHit_TimeVsQF_cut[i] = new TH2I("DigiHit_TimeVsQF_cut"+strN,"PS fADC pulse time vs. quality factor (> 500 ADC integral counts)"+strT+";fADC quality factor;pulse time [ns]",4,-0.5,3.5,1000,0.0,400.0);
166  digihitDir->cd();
167  }
168  // back to main dir
169  mainDir->cd();
170 
171  return NOERROR;
172 }
173 
174 
175 //----------------------------------------------------------------------------------
176 
177 
178 jerror_t JEventProcessor_PS_online::brun(JEventLoop *eventLoop, int32_t runnumber) {
179  // This is called whenever the run number changes
180 
181  // extract the PS geometry
182  vector<const DPSGeometry*> psGeomVect;
183  eventLoop->Get(psGeomVect);
184  if (psGeomVect.size() == 0) return OBJECT_NOT_AVAILABLE;
185  const DPSGeometry& psGeom = *(psGeomVect[0]);
186 
187  // Get photon energy bin lows for variable-width energy binning
188  double Elows[Narms][Ncols+1];
189  double cols[Ncols+1];
190  for (int i = 0; i < Ncols; i++) {
191  Elows[0][i] = psGeom.getElow(0,i+1);
192  Elows[1][i] = psGeom.getElow(1,i+1);
193  cols[i] = 0.5 + i;
194  }
195  // Add the upper limits
196  Elows[0][Ncols] = psGeom.getEhigh(0,Ncols);
197  Elows[1][Ncols] = psGeom.getEhigh(1,Ncols);
198  cols[Ncols] = 0.5 + Ncols;
199 
200  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
201  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
202 
203  // Set variable-width energy bins if histogram is empty
204  for (int i = 0; i < Narms; i++) {
205  if (hHit_Energy[i]->GetEntries() == 0) hHit_Energy[i]->SetBins(Ncols,Elows[i]);
206  if (hHit_EnergyVsColumn[i]->GetEntries() == 0) hHit_EnergyVsColumn[i]->SetBins(Ncols,cols,Ncols,Elows[i]);
207  }
208 
209  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
210 
211  return NOERROR;
212 }
213 
214 
215 //----------------------------------------------------------------------------------
216 
217 
218 jerror_t JEventProcessor_PS_online::evnt(JEventLoop *eventLoop, uint64_t eventnumber) {
219  // This is called for every event. Use of common resources like writing
220  // to a file or filling a histogram should be mutex protected. Using
221  // loop-Get(...) to get reconstructed objects (and thereby activating the
222  // reconstruction algorithm) should be done outside of any mutex lock
223  // since multiple threads may call this method at the same time.
224 
225  // get data for PS
226  vector<const DPSDigiHit*> digihits;
227  eventLoop->Get(digihits);
228  vector<const DPSHit*> hits;
229  eventLoop->Get(hits);
230 
231  // FILL HISTOGRAMS
232  // Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
233  japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
234 
235  if (digihits.size() > 0) ps_num_events->Fill(1);
236 
237  // Fill digihit hists
238  int NDigiHits[] = {0,0};
239  int NDigiHits_cut[] = {0,0};
240  hDigiHit_NHits->Fill(digihits.size());
241  for (const auto& hit : digihits) {
242  int arm = hit->arm;
243  double ped = (double)hit->pedestal/hit->nsamples_pedestal;
244  hDigiHit_NSamplesPedestal[arm]->Fill(hit->nsamples_pedestal);
245  hDigiHit_Pedestal[arm]->Fill(ped);
246  hDigiHit_RawPeak[arm]->Fill(hit->pulse_peak);
247  if (ped == 0.0 || hit->pulse_peak == 0) continue;
248  hDigiHit_PedestalVsColumn[arm]->Fill(hit->column,ped);
249  NDigiHits[arm]++;
250  hDigiHit_Arm->Fill(arm);
251  hDigiHit_Occupancy[arm]->Fill(hit->column);
252  hDigiHit_RawPeakVsColumn[arm]->Fill(hit->column,hit->pulse_peak);
253  hDigiHit_RawIntegral[arm]->Fill(hit->pulse_integral);
254  hDigiHit_RawIntegralVsColumn[arm]->Fill(hit->column,hit->pulse_integral);
255  hDigiHit_NSamplesIntegral[arm]->Fill(hit->nsamples_integral);
256  hDigiHit_PeakVsColumn[arm]->Fill(hit->column,hit->pulse_peak-ped);
257  double pI = hit->pulse_integral-hit->nsamples_integral*ped;
258  hDigiHit_IntegralVsColumn[arm]->Fill(hit->column,pI);
259  hDigiHit_IntegralVsPeak[arm]->Fill(hit->pulse_peak-ped,pI);
260  hDigiHit_PulseTime[arm]->Fill(hit->pulse_time);
261  double t_ns = 0.0625*hit->pulse_time;
262  hDigiHit_Time[arm]->Fill(t_ns);
263  hDigiHit_TimeVsColumn[arm]->Fill(hit->column,t_ns);
264  hDigiHit_TimeVsIntegral[arm]->Fill(pI,t_ns);
265  hDigiHit_QualityFactor[arm]->Fill(hit->QF);
266  if (pI > counts_cut) {
267  NDigiHits_cut[arm]++;
268  hDigiHit_Arm_cut->Fill(arm);
269  hDigiHit_Occupancy_cut[arm]->Fill(hit->column);
270  hDigiHit_Time_cut[arm]->Fill(t_ns);
271  hDigiHit_TimeVsColumn_cut[arm]->Fill(hit->column,t_ns);
272  hDigiHit_TimeVsQF_cut[arm]->Fill(hit->QF,t_ns);
273  }
274  }
275  hDigiHit_NHitsVsArm->Fill(0.,NDigiHits[0]); hDigiHit_NHitsVsArm->Fill(1.,NDigiHits[1]);
276  hDigiHit_NHits_cut->Fill(NDigiHits_cut[0]+NDigiHits_cut[1]);
277  hDigiHit_NHitsVsArm_cut->Fill(0.,NDigiHits_cut[0]); hDigiHit_NHitsVsArm_cut->Fill(1.,NDigiHits_cut[1]);
278 
279  // Fill calibrated-hit hists
280  int NHits[] = {0,0};
281  hHit_NHits->Fill(hits.size());
282  for (const auto& hit : hits) {
283  int arm = hit->arm;
284  NHits[arm]++;
285  hHit_Arm->Fill(arm);
286  hHit_Occupancy[arm]->Fill(hit->column);
287  hHit_Energy[arm]->Fill(hit->E);
288  hHit_EnergyVsColumn[arm]->Fill(hit->column,hit->E);
289  hHit_Integral[arm]->Fill(hit->integral);
290  hHit_IntegralVsColumn[arm]->Fill(hit->column,hit->integral);
291  //hHit_Npix[arm]->Fill(hit->npix_fadc);
292  hHit_Time[arm]->Fill(hit->t);
293  hHit_TimeVsColumn[arm]->Fill(hit->column,hit->t);
294  }
295  hHit_NHitsVsArm->Fill(0.,NHits[0]); hHit_NHitsVsArm->Fill(1.,NHits[1]);
296 
297  japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK
298 
299  return NOERROR;
300 }
301 
302 
303 //----------------------------------------------------------------------------------
304 
305 
307  // This is called whenever the run number changes, before it is
308  // changed to give you a chance to clean up before processing
309  // events from the next run number.
310  return NOERROR;
311 }
312 
313 
314 //----------------------------------------------------------------------------------
315 
316 
318  // Called before program exit after event processing is finished.
319  return NOERROR;
320 }
321 
322 
323 //----------------------------------------------------------------------------------
324 //----------------------------------------------------------------------------------
static TH1I * ps_num_events
static TH1I * hDigiHit_NHits_cut
double counts_cut
static TH2I * hDigiHit_RawPeakVsColumn[Narms]
static TH1I * hDigiHit_Arm_cut
const int Ncols
static const int NUM_ARMS
Definition: DPSGeometry.h:22
static TH1I * hDigiHit_Time[Narms]
static TH1I * hDigiHit_Pedestal[Narms]
static TH2I * hDigiHit_PeakVsColumn[Narms]
static TH1I * hDigiHit_Arm
static TH2I * hDigiHit_IntegralVsColumn[Narms]
double getEhigh(int arm, int column) const
Definition: DPSGeometry.cc:48
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
const int NmultBins
static TH2I * hDigiHit_NHitsVsArm
TDirectory * mainDir
Definition: p2k_hists.C:2
static TH1I * hHit_NHits
static TH2I * hDigiHit_TimeVsColumn_cut[Narms]
static TH2I * hDigiHit_RawIntegralVsColumn[Narms]
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
static TH2I * hDigiHit_TimeVsIntegral[Narms]
JApplication * japp
static TH1I * hDigiHit_RawIntegral[Narms]
static TH1I * hDigiHit_RawPeak[Narms]
static TH1I * hDigiHit_NSamplesPedestal[Narms]
static TH2I * hDigiHit_TimeVsQF_cut[Narms]
static TH1I * hDigiHit_Time_cut[Narms]
static TH2I * hHit_NHitsVsArm
InitPlugin_t InitPlugin
static TProfile * hDigiHit_PedestalVsColumn[Narms]
static TH2I * hDigiHit_TimeVsColumn[Narms]
static TH2I * hHit_IntegralVsColumn[Narms]
static TH1I * hDigiHit_NSamplesIntegral[Narms]
static TH1I * hHit_Energy[Narms]
static TH1I * hDigiHit_Occupancy_cut[Narms]
static TH1I * hHit_Occupancy[Narms]
const int Narms
static TH1I * hHit_Time[Narms]
static TH1I * hHit_Arm
static TH1I * hDigiHit_Occupancy[Narms]
static TH2I * hHit_TimeVsColumn[Narms]
static TH1I * hDigiHit_QualityFactor[Narms]
jerror_t init(void)
Called once at program start.
static TH2I * hHit_EnergyVsColumn[Narms]
static TH2I * hDigiHit_IntegralVsPeak[Narms]
static TH1I * hDigiHit_PulseTime[Narms]
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
jerror_t fini(void)
Called after last event of last event source has been processed.
static const int NUM_FINE_COLUMNS
Definition: DPSGeometry.h:26
static TH1I * hHit_Integral[Narms]
static TH2I * hDigiHit_NHitsVsArm_cut
static TH1I * hDigiHit_NHits
double getElow(int arm, int column) const
Definition: DPSGeometry.cc:40