Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_BCAL_SiPM_saturation.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_BCAL_SiPM_saturation.cc
4 // Modified file from BCAL_neutron_discriminator.cc ES 5/10/2018
5 // Created: Thu Apr 5 16:36:00 EDT 2018
6 // Creator: dalton (on Linux gluon119.jlab.org 2.6.32-642.3.1.el6.x86_64 x86_64)
7 //
8 
10 #include "HistogramTools.h"
11 #include "TRACKING/DMCThrown.h"
12 // #include "BCAL/DBCALHit.h"
13 // #include "BCAL/DBCALTDCHit.h"
14 // #include "BCAL/DBCALCluster.h"
15 // #include "BCAL/DBCALDigiHit.h"
16 #include "BCAL/DBCALPoint.h"
17 // #include "BCAL/DBCALUnifiedHit.h"
18 // #include "BCAL/DBCALGeometry.h"
19 #include "BCAL/DBCALShower.h"
20 #include "DANA/DStatusBits.h"
21 // #include "PID/DChargedTrack.h"
22 // #include "PID/DEventRFBunch.h"
23 // #include "PID/DDetectorMatches.h"
24 #include "PID/DNeutralShower.h"
25 // #include "PID/DVertex.h"
26 // #include "TRACKING/DTrackTimeBased.h"
27 // #include "TRIGGER/DL1Trigger.h"
28 using namespace jana;
29 
30 
31 // Routine used to create our JEventProcessor
32 #include <JANA/JApplication.h>
33 #include <JANA/JFactory.h>
34 extern "C"{
35 void InitPlugin(JApplication *app){
36  InitJANAPlugin(app);
37  app->AddProcessor(new JEventProcessor_BCAL_SiPM_saturation());
38 }
39 } // "C"
40 
41 
42 //------------------
43 // JEventProcessor_BCAL_SiPM_saturation (Constructor)
44 //------------------
46 {
47  VERBOSE = 0;
48  if(gPARMS){
49  gPARMS->SetDefaultParameter("BCAL_SiPM_saturation:VERBOSE", VERBOSE, "Verbosity level");
50  }
51 }
52 
53 //------------------
54 // ~JEventProcessor_BCAL_SiPM_saturation (Destructor)
55 //------------------
57 {
58 
59 }
60 
61 //------------------
62 // init
63 //------------------
65 {
66  // This is called once at program startup.
67 
68  return NOERROR;
69 }
70 
71 //------------------
72 // brun
73 //------------------
74 jerror_t JEventProcessor_BCAL_SiPM_saturation::brun(JEventLoop *eventLoop, int32_t runnumber)
75 {
76  // This is called whenever the run number changes
77 
78  attenuation_parameters.clear();
79  eventLoop->GetCalib("/BCAL/attenuation_parameters", attenuation_parameters);
80  /*int channel = 0;
81  for (int module=1; module<=48; module++) {
82  for (int layer=1; layer<=4; layer++) {
83  for (int sector=1; sector<=4; sector++) {
84  printf("%2i %2i %2i %12.4f %12.4f %12.4f\n",
85  module,layer,sector,
86  attenuation_parameters[channel][0],
87  attenuation_parameters[channel][1],
88  attenuation_parameters[channel][2]);
89  }
90  channel++;
91  }
92  }*/
93 
94  return NOERROR;
95 }
96 
97 //------------------
98 // evnt
99 //------------------
100 jerror_t JEventProcessor_BCAL_SiPM_saturation::evnt(JEventLoop *loop, uint64_t eventnumber)
101 {
102 
103  // If not Physics event call BCAL showers to initialize parameters
104  if (!loop->GetJEvent().GetStatusBit(kSTATUS_PHYSICS_EVENT)) {
105  vector<const DBCALShower*> BCALShowers;
106  loop->Get(BCALShowers);
107  }
108 
109  // Check to see if this is a physics event
110  if (loop->GetJEvent().GetStatusBit(kSTATUS_PHYSICS_EVENT)) {
111 
112  // Check if thrown information is available
113  vector<const DMCThrown*> MCThrowns;
114  loop->Get(MCThrowns);
115  unsigned int NumThrown=MCThrowns.size();
116  double Ethrown=0;
117  double thetathrown=0;
118  if (NumThrown == 1) {
119  const DMCThrown* locMCThrown = MCThrowns[0];
120  Ethrown = locMCThrown->energy();
121  DVector3 pgen = locMCThrown->momentum();
122  //pgen.Print();
123  thetathrown = 180./3.14159*pgen.Theta();
124  }
125 
126  if (NumThrown==0 || thetathrown > 12) { // Either data or ignore MC photons hitting end of bcal
127 
128  // Get vector of neutral showers in event
129  vector<const DNeutralShower*> NeutralShowers;
130  loop->Get(NeutralShowers);
131 
132  unsigned int NumShowers=NeutralShowers.size();
133  float Eshower = 0;
134 
135  // loop over neutral showers
136  for (unsigned int i = 0; i < NumShowers; i++){
137  const DNeutralShower* locNeutralShower = NeutralShowers[i];
138  // Must be BCAL shower
139  DetectorSystem_t locDetector = locNeutralShower->dDetectorSystem;
140  if (locDetector != SYS_BCAL) continue;
141  // Get shower properties
142  vector<const DBCALShower*> BCALShowers;
143  locNeutralShower->Get(BCALShowers);
144 
145  // Should be only one BCAL shower for each neutral shower
146  const DBCALShower* locBCALShower= BCALShowers[0];
147  Eshower = locBCALShower->E;
148 
149  // for MC select showers that are greater than 50% of thrown energy
150  if (NumThrown == 1 && Eshower < Ethrown/2.)
151  continue;
152 
153  /*float E_preshower = locBCALShower->E_preshower;
154  float z = locBCALShower->z;
155  float x = locBCALShower->x;
156  float y = locBCALShower->y;
157  float R = sqrt(x*x+y*y);
158  float sigLong = locBCALShower->sigLong; // rho
159  float sigTrans = locBCALShower->sigTrans; // phi
160  float sigTheta = locBCALShower->sigTheta;*/
161 
162  // cout << " Shower i=" << i << " E_shower=" << Eshower << " E_preshower=" << E_preshower << " x=" << x << " y=" << y << " z=" << z << " R=" << R << " sigLong=" << sigLong << " sigTrans=" << sigTrans << " sigTheta=" << sigTheta << endl;
163 
164  Int_t nbins=100;
165 
166  // Fill histogram for showers
167  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Ethrown", Ethrown,
168  "BCAL SiPM Saturation; Thrown Energy (GeV)",4*nbins,0,10);
169 
170  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Eshower", Eshower,
171  "BCAL SiPM Saturation; Shower Energy (GeV)",10*nbins,0,10);
172  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Thrown Theta", thetathrown,
173  "BCAL SiPM Saturation; Thrown Theta (degrees)",4*90,0,90);
174 
175  // Fill 2D histograms
176  Fill2DHistogram ("BCAL_SiPM_saturation", "Hists2D", "Eshower_vs_Ethrown", Ethrown, Eshower,
177  "BCAL SiPM Saturation; Thrown Energy (GeV); Shower Energy (GeV)",
178  4*nbins,0,10,4*nbins,0,10);
179  Fill2DHistogram ("BCAL_SiPM_saturation", "Hists2D", "EDiff_vs_Ethrown", Ethrown, Eshower - Ethrown,
180  "BCAL SiPM Saturation; Thrown Energy (GeV); (Shower - Thrown) Energy (GeV)",
181  4*nbins,0,10,nbins,-0.5,0.5);
182  Fill2DHistogram ("BCAL_SiPM_saturation", "Hists2D", "EDiff/Ethrown_vs_Ethrown", Ethrown, Ethrown>0? (Eshower - Ethrown)/Ethrown : 0,
183  "BCAL SiPM Saturation; Thrown Energy (GeV); (EShower - EThrown)/Ethrown",
184  4*nbins,0,10,nbins,-0.2,0.2);
185 
186  // Get vector of points in this shower
187  vector<const DBCALPoint*> Points;
188  locNeutralShower->Get(Points);
189  uint Ncell = Points.size();
190 
191  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "NCell", Ncell,
192  "BCAL SiPM Saturation; Number of cells",nbins,0,100);
193 
194  for (unsigned int j = 0; j < Ncell; j++){
195  const DBCALPoint* locPoint = Points[j];
196  float t = locPoint->t();
197  float z = locPoint->z();
198  float Ept = locPoint->E();
199  int module= locPoint->module();
200  int layer = locPoint->layer();
201  int sector = locPoint->sector();
202  // cout << " j=" << j << " t=" << t << endl;
203 
204  // calculate point energy (code taken from DBCALPoint.cc)
205  float z_bcal_center = 212;
206  float fibLen = 390;
207  float zLocal = z - z_bcal_center;
208 
209 
210  float dUp = 0.5 * fibLen + zLocal;
211  float dDown = 0.5 * fibLen - zLocal;
212  if (dUp>fibLen) dUp=fibLen;
213  if (dUp<0) dUp=0;
214  if (dDown>fibLen) dDown=fibLen;
215  if (dDown<0) dDown=0;
216  int channel = (module-1)*16 + (layer-1)*4 + (sector-1);
217  // cout << " module=" << module << " layer=" << layer << " sector=" << sector << " channel=" << channel << endl;
218  // cout << " attenuation length 00=" << attenuation_parameters[channel][0] << endl;
219  float attenuation_length = attenuation_parameters[channel][0];
220  float attUp = exp( -dUp / attenuation_length );
221  float attDown = exp( -dDown / attenuation_length );
222 
223 
224  if (VERBOSE>=3) cout << " VERBOSE >=3" << " t=" << t << " z=" << z << endl;
225 
226 
227  // Fill 1D histograms
228 
229  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "layer", layer,
230  "BCAL SiPM Saturation; Layer Number",5,0,5);
231 
232  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Ept", Ept,
233  "BCAL SiPM Saturation; Point Energy (GeV)",10*nbins,0,10);
234 
235  // cout << " Point: Ept=" << Ept << endl;
236  float upHit=0;
237  float downHit=0;
238  float uppeak=0;
239  float downpeak=0;
240 
241  vector<const DBCALHit*> Hits;
242  locPoint->Get(Hits);
243  uint Nhits = Hits.size();
244 
245  for (unsigned int j = 0; j < Nhits; j++){
246  const DBCALHit* locHit = Hits[j];
247  float Ehit = locHit->E;
248  /*int module = locHit->module;
249  int layer = locHit->layer;
250  int sector = locHit->sector;*/
251  int end = locHit->end;
252  int pulse_peak = locHit->pulse_peak;
253 
254  // cout << " module=" << module << " layer=" << layer << " sector=" << sector << " pulse_peak/Ehit=" << pulse_peak/Ehit << endl;
255 
256  if (end == 0) {
257  upHit = Ehit;
258  uppeak = pulse_peak;
259  }
260  if (end == 1) {
261  downHit = Ehit;
262  downpeak = pulse_peak;
263  }
264 
265 
266  nbins=5100;
267 
268  // Fill 1D histograms
269  if (layer == 1) {
270  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit pulse_peak layer=1", pulse_peak,
271  "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
272 
273  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit integral layer=1", Ehit,
274  "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
275  }
276  else if (layer == 2) {
277  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit pulse_peak layer=2", pulse_peak,
278  "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
279 
280  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit integral layer=2", Ehit,
281  "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
282  }
283  else if (layer == 3) {
284  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit pulse_peak layer=3", pulse_peak,
285  "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
286 
287  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit integral layer=3", Ehit,
288  "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
289  }
290  else if (layer == 4) {
291  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit pulse_peak layer=4", pulse_peak,
292  "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
293 
294  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Hit integral layer=4", Ehit,
295  "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
296  }
297  else {
298  cout << " ***Illegal layer=" << layer << endl;
299  }
300 
301 
302  }
303 
304 
305 
306  // use these to correct the energy
307  float E_US = ( upHit / attUp );
308  float E_DS = ( downHit / attDown );
309  float Ecalc = ( E_US + E_DS ) / 2;
310 
311  /*cout << " i=" << i << " module=" << module << " layer=" << layer << " sector=" << sector
312  << " lambda=" << attenuation_length << " Eup=" << upHit << " Edown=" << downHit << " Point: Ept="
313  << Ept << " Ecalc=" << Ecalc << " Diff=" << Ept-Ecalc << endl;*/
314 
315  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Ecalc", Ecalc,
316  "BCAL SiPM Saturation; Calc Energy (GeV)",4*nbins,0,10);
317  Fill1DHistogram ("BCAL_SiPM_saturation", "Hists1D", "Ecalc-Ept", Ecalc-Ept,
318  "BCAL SiPM Saturation; Calc-Pt Energy (GeV)",nbins,-0.05,0.05);
319 
320  }
321 
322 
323 
324  }
325  }
326  }
327  return NOERROR;
328 }
329 
330 //------------------
331 // erun
332 //------------------
334 {
335  // This is called whenever the run number changes, before it is
336  // changed to give you a chance to clean up before processing
337  // events from the next run number.
338  return NOERROR;
339 }
340 
341 //------------------
342 // fini
343 //------------------
345 {
346  // Called before program exit after event processing is finished.
347  return NOERROR;
348 }
349 
float E
Definition: DBCALHit.h:30
int module() const
Definition: DBCALPoint.h:64
Int_t layer
TVector3 DVector3
Definition: DVector3.h:14
double energy(void) const
int sector() const
Definition: DBCALPoint.h:66
DetectorSystem_t
Definition: GlueX.h:15
int layer() const
Definition: DBCALPoint.h:65
Definition: GlueX.h:19
DetectorSystem_t dDetectorSystem
DBCALGeometry::End end
Definition: DBCALHit.h:28
float E() const
Definition: DBCALPoint.h:38
InitPlugin_t InitPlugin
const bool VERBOSE
void Fill1DHistogram(const char *plugin, const char *directoryName, const char *name, const double value, const char *title, int nBins, double xmin, double xmax, bool print=false)
float z() const
Definition: DBCALPoint.h:60
void Fill2DHistogram(const char *plugin, const char *directoryName, const char *name, const double valueX, const double valueY, const char *title, int nBinsX, double xmin, double xmax, int nBinsY, double ymin, double ymax, bool print=false)
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
const DVector3 & momentum(void) const
jerror_t fini(void)
Called after last event of last event source has been processed.
int pulse_peak
Definition: DBCALHit.h:29
float t() const
Definition: DBCALPoint.h:41
jerror_t init(void)
Called once at program start.