Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_BCAL_TDC_Timing.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_BCAL_TDC_Timing.cc
4 // Created: Tue Jul 28 10:55:56 EDT 2015
5 // Creator: mstaib (on Linux egbert 2.6.32-504.30.3.el6.x86_64 x86_64)
6 //
7 
8 /*************************************
9  This plugin is designed to calibrate the TDC times for the BCAL.
10  These calibrations include the time-walk for each TDC channel, effective velocities and relative offsets between ends of the detector.
11  Creation of constants from histograms is controlled by ROOT scripts in the FitScripts directory.
12  *************************************/
13 
15 #include "HistogramTools.h" // To make my life easier
16 #include "BCAL/DBCALHit.h"
17 #include "BCAL/DBCALTDCHit.h"
18 #include "BCAL/DBCALCluster.h"
19 #include "BCAL/DBCALDigiHit.h"
20 #include "BCAL/DBCALPoint.h"
21 #include "BCAL/DBCALUnifiedHit.h"
22 #include "BCAL/DBCALGeometry.h"
23 #include "DANA/DStatusBits.h"
24 #include "PID/DChargedTrack.h"
25 #include "PID/DEventRFBunch.h"
26 #include "PID/DDetectorMatches.h"
27 #include "PID/DNeutralShower.h"
28 #include "PID/DVertex.h"
30 #include "TRIGGER/DL1Trigger.h"
31 
32 using namespace jana;
33 
34 #include "DANA/DApplication.h"
35 // Routine used to create our JEventProcessor
36 #include <JANA/JApplication.h>
37 #include <JANA/JFactory.h>
38 extern "C"{
39 void InitPlugin(JApplication *app){
40  InitJANAPlugin(app);
41  app->AddProcessor(new JEventProcessor_BCAL_TDC_Timing());
42 }
43 } // "C"
44 
45 
46 //------------------
47 // JEventProcessor_BCAL_TDC_Timing (Constructor)
48 //------------------
50 {
51  VERBOSE = 0;
52  VERBOSEHISTOGRAMS = 0;
53 
54  if(gPARMS){
55  gPARMS->SetDefaultParameter("BCAL_TDC_Timing:VERBOSE", VERBOSE, "Verbosity level");
56  gPARMS->SetDefaultParameter("BCAL_TDC_Timing:VERBOSEHISTOGRAMS", VERBOSEHISTOGRAMS, "Create more histograms (default 0 for monitoring)");
57  }
58 
59 }
60 
61 //------------------
62 // ~JEventProcessor_BCAL_TDC_Timing (Destructor)
63 //------------------
65 {
66 
67 }
68 
69 //------------------
70 // init
71 //------------------
73 {
74  // This is called once at program startup on a single thread.
75 
76  return NOERROR;
77 }
78 
79 //------------------
80 // brun
81 //------------------
82 jerror_t JEventProcessor_BCAL_TDC_Timing::brun(JEventLoop *loop, int32_t runnumber)
83 {
84  // This is called whenever the run number changes
85  DApplication* app = dynamic_cast<DApplication*>(loop->GetJApplication());
86  DGeometry* geom = app->GetDGeometry(runnumber);
87  geom->GetTargetZ(Z_TARGET);
88 
89  // load BCAL geometry
90  vector<const DBCALGeometry *> BCALGeomVec;
91  loop->Get(BCALGeomVec);
92  if(BCALGeomVec.size() == 0)
93  throw JException("Could not load DBCALGeometry object!");
94  auto locBCALGeom = BCALGeomVec[0];
95 
96  printf("locBCALGeom->GetBCAL_center()=%f\nZ_TARGET=%f\n",
97  locBCALGeom->GetBCAL_center(), Z_TARGET);
98 
99  // //////
100  // // THIS NEEDS TO CHANGE IF THERE IS A NEW TABLE
101  // //////
102  // //get timewalk corrections from CCDB
103  // JCalibration *jcalib = loop->GetJCalibration();
104  // //these tables hold: module layer sector end c0 c1 c2 c3
105  // vector<vector<float> > tdc_timewalk_table;
106  // jcalib->Get("BCAL/timewalk_tdc",tdc_timewalk_table);
107 
108  // for (vector<vector<float> >::const_iterator iter = tdc_timewalk_table.begin();
109  // iter != tdc_timewalk_table.end();
110  // ++iter) {
111  // if (iter->size() != 8) {
112  // cout << "DBCALUnifiedHit_factory: Wrong number of values in timewalk_tdc table (should be 8)" << endl;
113  // continue;
114  // }
115  // //be really careful about float->int conversions
116  // int module = (int)((*iter)[0]+0.5);
117  // int layer = (int)((*iter)[1]+0.5);
118  // int sector = (int)((*iter)[2]+0.5);
119  // int endi = (int)((*iter)[3]+0.5);
120  // DBCALGeometry::End end = (endi==0) ? DBCALGeometry::kUpstream : DBCALGeometry::kDownstream;
121  // float c0 = (*iter)[4];
122  // float c1 = (*iter)[5];
123  // float c2 = (*iter)[6];
124  // float a_thresh = (*iter)[7];
125  // int cellId = locBCALGeom->cellId(module, layer, sector);
126  // readout_channel channel(cellId,end);
127  // tdc_timewalk_map[channel] = timewalk_coefficients(c0,c1,c2,a_thresh);
128  // }
129 
130  /// Read in initial calibration constants and write to root file for use in later calibration
131  vector<double> raw_channel_global_offset;
132  //if(print_messages) jout << "In BCAL_TDC_Timing, loading constants..." << endl;
133  if(loop->GetCalib("/BCAL/channel_global_offset", raw_channel_global_offset))
134  jout << "Error loading /BCAL/channel_global_offset !" << endl;
135 
136  japp->RootFillLock(this);
137  //TH1D *CCDB_raw_channel_global_offset = new TH1D("CCDB_raw_channel_global_offset","Offsets at time of running;channel;offset",768,0.5,768.5);
138  int counter = 1;
139  Fill1DWeightedHistogram("BCAL_Global_Offsets", "Target Time", "CCDB_raw_channel_global_offset",
140  769, 1,
141  "Offsets at time of running;CCDB Index;CCDB timing offset [ns]",
142  769, 0.5, 769.5);
143  for (vector<double>::iterator iter = raw_channel_global_offset.begin(); iter != raw_channel_global_offset.end(); ++iter) {
144  //CCDB_raw_channel_global_offset->SetBinContent(counter,*iter);
145  Fill1DWeightedHistogram("BCAL_Global_Offsets", "Target Time", "CCDB_raw_channel_global_offset",
146  counter, *iter,
147  "Offsets at time of running;CCDB Index;CCDB timing offset [ns]",
148  769, 0.5, 769.5);
149  counter++;
150  }
151  japp->RootFillUnLock(this);
152 
153  return NOERROR;
154 }
155 
156 //------------------
157 // evnt
158 //------------------
159 jerror_t JEventProcessor_BCAL_TDC_Timing::evnt(JEventLoop *loop, uint64_t eventnumber)
160 {
161  // First check that this is not a font panel trigger or no trigger
162  bool goodtrigger=1;
163 
164  const DL1Trigger *trig = NULL;
165  try {
166  loop->GetSingle(trig);
167  } catch (...) {}
168  if (trig) {
169  if (trig->fp_trig_mask){
170  goodtrigger=0;
171  }
172  } else {
173  // HDDM files are from simulation, so keep them even though they have no trigger
174  bool locIsHDDMEvent = loop->GetJEvent().GetStatusBit(kSTATUS_HDDM);
175  if (!locIsHDDMEvent) goodtrigger=0;
176  }
177 
178  // load BCAL geometry
179  vector<const DBCALGeometry *> BCALGeomVec;
180  loop->Get(BCALGeomVec);
181  if(BCALGeomVec.size() == 0)
182  throw JException("Could not load DBCALGeometry object!");
183  auto locBCALGeom = BCALGeomVec[0];
184 
185  vector<const DTrackFitter *> fitters;
186  loop->Get(fitters);
187 
188  if(fitters.size()<1){
189  _DBG_<<"Unable to get a DTrackFinder object!"<<endl;
190  return RESOURCE_UNAVAILABLE;
191  }
192 
193  auto fitter = fitters[0];
194 
195 
196 
197  // calculate total BCAL energy in order to catch BCAL LED events
198  vector<const DBCALHit *> bcal_hits;
199  loop->Get(bcal_hits);
200  double total_bcal_energy = 0.;
201  for(unsigned int i=0; i<bcal_hits.size(); i++) {
202  total_bcal_energy += bcal_hits[i]->E;
203  }
204  if (total_bcal_energy > 12.) goodtrigger=0;
205 
206  if (!goodtrigger) {
207  return NOERROR;
208  }
209 
210  vector<const DBCALUnifiedHit *> bcalUnifiedHitVector;
211  loop->Get(bcalUnifiedHitVector);
212 
213  /**********************************************
214  _____ _ __ __ _ _
215  |_ _(_)_ __ ___ ___ \ \ / /_ _| | | __
216  | | | | '_ ` _ \ / _ \____\ \ /\ / / _` | | |/ /
217  | | | | | | | | | __/_____\ V V / (_| | | <
218  |_| |_|_| |_| |_|\___| \_/\_/ \__,_|_|_|\_\
219  ********************************************/
220  // The very first thing to do is correct for the timewalk. If this calibration is screwed up, the calibrations that follow will be wrong...
221  // The following plots can be used for the calibration or to check the calibration if it has already been done
222 
223  double MIN_TDIFF = -10.0, MAX_TDIFF = 10.0, MAX_TDIFF_WIDE = 30.0;
224  const int ndtbins = 100;
225  double dtbins[ndtbins+1];
226  for (int i=0; i<=ndtbins; i++) {
227  dtbins[i] = MIN_TDIFF + (MAX_TDIFF-MIN_TDIFF)/ndtbins*i;
228  }
229  const int ndtbinswide = 200;
230  double dtbinswide[ndtbinswide+1];
231  for (int i=0; i<=ndtbinswide; i++) {
232  dtbinswide[i] = MIN_TDIFF + (MAX_TDIFF_WIDE-MIN_TDIFF)/ndtbinswide*i;
233  }
234  const int npeakbins = 100;
235  double peakbins[npeakbins+1] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,20,
236  21,23,24,26,28,30,32,34,36,39,42,44,48,51,54,58,62,66,71,76,
237  81,86,92,98,105,112,120,128,137,146,156,167,178,190,203,217,
238  231,247,264,281,300,321,343,366,390,417,445,475,507,541,578,
239  617,659,703,750,801,855,913,974,1040,1110,1185,1265,1351,1442,
240  1539,1643,1754,1872,1998,2133,2277,2430,2594,2769,2956,3155,
241  3368,3595,3837,4096};
242 
243  for (unsigned int i = 0; i < bcalUnifiedHitVector.size(); i++){
244  //int the_cell = (bcalUnifiedHitVector[i]->module - 1) * 16 + (bcalUnifiedHitVector[i]->layer - 1) * 4 + bcalUnifiedHitVector[i]->sector;
245  // There is one less layer of TDCs so the numbering relects this
246  int the_tdc_cell = (bcalUnifiedHitVector[i]->module - 1) * 12 + (bcalUnifiedHitVector[i]->layer - 1) * 4 + bcalUnifiedHitVector[i]->sector;
247  //int cellId = locBCALGeom->cellId(bcalUnifiedHitVector[i]->module, bcalUnifiedHitVector[i]->layer, bcalUnifiedHitVector[i]->sector);
248  // Get the underlying associated objects
249  const DBCALHit * thisADCHit;
250  const DBCALTDCHit * thisTDCHit;
251  bcalUnifiedHitVector[i]->GetSingle(thisADCHit);
252  bcalUnifiedHitVector[i]->GetSingle(thisTDCHit);
253 
254  // From the ADC hit we can extract the pulse peak information...
255  int pulse_peak = thisADCHit->pulse_peak;
256 
257  // The raw information from the DBCALHit and DBCALTDCHit is not corrected for timewalk yet, so we can always plot the before and after.
258  if (thisADCHit != NULL && thisTDCHit != NULL){
259  char name[200];
260  sprintf(name, "M%02iL%iS%i", bcalUnifiedHitVector[i]->module, bcalUnifiedHitVector[i]->layer, bcalUnifiedHitVector[i]->sector);
261  if (bcalUnifiedHitVector[i]->end == 0){
262  Fill2DHistogram ("BCAL_TDC_Timing", "Upstream_TimewalkVsPeak", name,
263  pulse_peak, thisTDCHit->t - thisADCHit->t,
264  "Timewalk; Pulse Peak [ADC Counts]; t_{TDC} - t_{ADC} [ns]",
265  npeakbins, peakbins, ndtbinswide, dtbinswide);
266  Fill2DHistogram ("BCAL_TDC_Timing", "Timewalk_All", "Upstream_Channel_Deltat",
267  the_tdc_cell, thisTDCHit->t - thisADCHit->t,
268  "BCAL Upstream t_{TDC}-t_{ADC}; cellID; t_{TDC} - t_{ADC} [ns] ",
269  576, 0.5, 576.5, ndtbins, MIN_TDIFF, MAX_TDIFF);
270  }
271 
272  else{
273  Fill2DHistogram ("BCAL_TDC_Timing", "Downstream_TimewalkVsPeak", name,
274  pulse_peak, thisTDCHit->t - thisADCHit->t,
275  "Timewalk; Pulse Peak [ADC Counts]; t_{TDC} - t_{ADC} [ns]",
276  npeakbins, peakbins, ndtbinswide, dtbinswide);
277  Fill2DHistogram ("BCAL_TDC_Timing", "Timewalk_All", "Downstream_Channel_Deltat",
278  the_tdc_cell, thisTDCHit->t - thisADCHit->t,
279  "BCAL Upstream t_{TDC}-t_{ADC}; cellID; t_{TDC} - t_{ADC} [ns] ",
280  576, 0.5, 576.5, ndtbins, MIN_TDIFF, MAX_TDIFF);
281  }
282  }
283  // Next look directly at the DBCALUnifiedHit to get the corrected times and plot those seperately
284  if (bcalUnifiedHitVector[i]->has_TDC_hit){
285  double correctedTDCTime = bcalUnifiedHitVector[i]->t_TDC;
286  char name[200];
287  sprintf(name, "M%02iL%iS%i", bcalUnifiedHitVector[i]->module, bcalUnifiedHitVector[i]->layer, bcalUnifiedHitVector[i]->sector);
288  if (bcalUnifiedHitVector[i]->end == 0){
289  Fill2DHistogram ("BCAL_TDC_Timing", "Upstream_TimewalkVsPeak_Corrected", name,
290  pulse_peak, correctedTDCTime - bcalUnifiedHitVector[i]->t_ADC,
291  "Timewalk; Pulse Peak [ADC Counts]; t_{TDC} - t_{ADC} [ns]",
292  npeakbins, peakbins, ndtbins ,dtbins);
293  Fill2DHistogram ("BCAL_TDC_Timing", "Timewalk_All", "Upstream_Channel_Deltat_TimewalkCorrected",
294  the_tdc_cell, correctedTDCTime - bcalUnifiedHitVector[i]->t_ADC,
295  "BCAL Upstream t_{TDC}-t_{ADC} corrected; cellID; t_{TDC} - t_{ADC} [ns] ",
296  576, 0.5, 576.5, ndtbins, MIN_TDIFF, MAX_TDIFF);
297  }
298  else{
299  Fill2DHistogram ("BCAL_TDC_Timing", "Downstream_TimewalkVsPeak_Corrected", name,
300  pulse_peak, correctedTDCTime - bcalUnifiedHitVector[i]->t_ADC,
301  "Timewalk; Pulse Peak [ADC Counts]; t_{TDC} - t_{ADC} [ns]",
302  npeakbins, peakbins, ndtbins ,dtbins);
303  Fill2DHistogram ("BCAL_TDC_Timing", "Timewalk_All", "Downstream_Channel_Deltat_TimewalkCorrected",
304  the_tdc_cell, correctedTDCTime - bcalUnifiedHitVector[i]->t_ADC,
305  "BCAL Downstream t_{TDC}-t_{ADC} corrected; cellID; t_{TDC} - t_{ADC} [ns] ",
306  576, 0.5, 576.5, ndtbins, MIN_TDIFF, MAX_TDIFF);
307  }
308  }
309  }
310 
311 
312  //
313  // Attenuation Length
314  //
315  // vector<const DBCALPoint *> BCALPointVector;
316  // loop->Get(BCALPointVector);
317 
318  // for (unsigned int i = 0; i < BCALPointVector.size(); i++){
319  // const DBCALPoint *thisPoint = BCALPointVector[i];
320  // vector<const DBCALDigiHit*> digihits;
321  // thisPoint->Get(digihits);
322  // if (digihits.size()!=2) {
323  // printf("Warning: BCAL_attenlength_gainratio: event %llu: wrong number of BCALDigiHit objects found %i\n",
324  // (long long unsigned int)eventnumber,(int)digihits.size());
325  // continue;
326  // }
327  // if (digihits[0]->end==digihits[1]->end) {
328  // printf("Warning: BCAL_attenlength_gainratio: event %llu: two hits in same end of point\n",(long long unsigned int)eventnumber);
329  // continue;
330  // }
331  // float integralUS, integralDS;
332  // // end 0=upstream, 1=downstream
333  // if (digihits[0]->end==0) {
334  // integralUS = digihits[0]->pulse_integral - ((float)digihits[0]->nsamples_integral*(float)digihits[0]->pedestal)/
335  // (float)digihits[0]->nsamples_pedestal;
336  // integralDS = digihits[1]->pulse_integral - ((float)digihits[1]->nsamples_integral*(float)digihits[1]->pedestal)/
337  // (float)digihits[1]->nsamples_pedestal;
338  // } else {
339  // integralDS = digihits[0]->pulse_integral - ((float)digihits[0]->nsamples_integral*(float)digihits[0]->pedestal)/
340  // (float)digihits[0]->nsamples_pedestal;
341  // integralUS = digihits[1]->pulse_integral - ((float)digihits[1]->nsamples_integral*(float)digihits[1]->pedestal)/
342  // (float)digihits[1]->nsamples_pedestal;
343  // }
344  // float intratio = (float)integralUS/(float)integralDS;
345  // float logintratio = log(intratio);
346 
347  // float zminlocal = -250;
348  // float zmaxlocal = 250;
349  // char name[200], title[200];
350  // sprintf(title,"Attenuation;Z_{Track} (cm);log of integral ratio US/DS");
351  // Fill2DHistogram ("BCAL_atten_gain", "logintratiovsZtrack", "AllPoints",
352  // thisPoint->z, logintratio, title,
353  // 250, zminlocal, zmaxlocal, 250, -3, 3);
354  // sprintf(title,"Attenuation (M%i,L%i,S%i);Z_{Track} (cm);log of integral ratio US/DS",
355  // thisPoint->module(), thisPoint->layer(), thisPoint->sector());
356  // Fill2DHistogram ("BCAL_atten_gain", "logintratiovsZtrack", channame,
357  // thisPoint->z, logintratio, title,
358  // 250, zminlocal, zmaxlocal, 250, -3, 3);
359  // }
360 
361 
362 
363 
364  /*************************************************
365  _________ _____ _______ _
366  /_ __/ _ \/ ___/ /_ __(_)_ _ (_)__ ___ _
367  / / / // / /__ / / / / ' \/ / _ \/ _ `/
368  /_/ /____/\___/ /_/ /_/_/_/_/_/_//_/\_, /
369  /___/
370  **************************************************/
371 
372  // Now we will use the track matching in order to work out the time offsets and effective velocities
373  // We just need to grab the charged particles in order to get their detector matches
374  // Note that this method is only really applicable for runs with the solenoid on...
375 
376  // We need the RF bunch for the event in order to check the global alignemtn of the timing
377  const DEventRFBunch *thisRFBunch = NULL;
378  loop->GetSingle(thisRFBunch);
379 
380  vector <const DChargedTrack *> chargedTrackVector;
381  loop->Get(chargedTrackVector);
382 
383  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 1, "Success profile;Step", 16, -0.5, 15.5);
384  for (unsigned int iTrack = 0; iTrack < chargedTrackVector.size(); iTrack++){
385  // Pick out the best charged track hypothesis for this charged track based only on the Tracking FOM
386  // const DChargedTrackHypothesis* bestHypothesis = chargedTrackVector[iTrack]->Get_BestTrackingFOM();
387  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 2, "Success profile;Step", 16, -0.5, 15.5);
388  // get charge and choose pion hypothesis as most likely
389  int charge = chargedTrackVector[iTrack]->Get_Charge();
390  char q[2];
391  const DChargedTrackHypothesis* bestHypothesis;
392  if (charge>0) {
393  bestHypothesis = chargedTrackVector[iTrack]->Get_Hypothesis(PiPlus);
394  sprintf(q,"+");
395  } else {
396  bestHypothesis = chargedTrackVector[iTrack]->Get_Hypothesis(PiMinus);
397  sprintf(q,"-");
398  }
399  if (bestHypothesis == NULL) continue;
400  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 3, "Success profile;Step", 16, -0.5, 15.5);
401 
402  // Now from this hypothesis we can get the detector matches to the BCAL
403  auto bcalMatch = bestHypothesis->Get_BCALShowerMatchParams();
404  auto scMatch = bestHypothesis->Get_SCHitMatchParams(); // Needed for quality cut later
405  DVector3 position = bestHypothesis->position();
406  //DVector3 momentum = bestHypothesis->momentum();
407  //float Z_track = position.z();
408  //float_track = momentum.Mag();
409 
410  if (bcalMatch == NULL) continue;
411  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 4, "Success profile;Step", 16, -0.5, 15.5);
412  if (scMatch == NULL) continue;
413  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 5, "Success profile;Step", 16, -0.5, 15.5);
414 
415  double dDeltaZToShower = bcalMatch->dDeltaZToShower;
416  double dDeltaPhiToShower = bcalMatch->dDeltaPhiToShower;
417 
418  Fill2DHistogram("BCAL_Global_Offsets", "Showers", "BCAL Match",
419  dDeltaZToShower, dDeltaPhiToShower, "BCAL Match;#Delta Z [cm]; #Delta#phi [rad]",
420  200, -25, 25, 200, -0.1, 0.1);
421 
422  const DTrackTimeBased *timeBasedTrack = bestHypothesis->Get_TrackTimeBased();
423  if (timeBasedTrack->FOM < 0.0027) continue; // 3-sigma cut on tracking FOM
424  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 6, "Success profile;Step", 16, -0.5, 15.5);
425  if (timeBasedTrack->Ndof < 10) continue; // CDC: 5 params in fit, 10 dof => [15 hits]; FDC [10 hits]
426 
427  // Use CDC dEdx to help reject protons
428  double dEdx=1e6*timeBasedTrack->ddEdx_CDC_amp;
429  double P_track=timeBasedTrack->momentum().Mag();
430  bool dEdx_pion = 0;
431  if (dEdx<2.5) dEdx_pion = 1;
432 
433  // Get the shower from the match
434  const DBCALShower *thisShower = bcalMatch->dBCALShower;
435 
436  // Fill histograms based on the shower
437  char name[200], title[200];
438  DVector3 proj_pos;
439  double flightTime=0.;
440  //double innerpathLength, innerflightTime;
441  double shower_x = thisShower->x;
442  double shower_y = thisShower->y;
443 
444  double t_shower = thisShower->t;
445  double E_shower = thisShower->E;
446  double Z_shower = thisShower->z;
447  // UNPROJECTION CODE
448  // double shower_deltaz = thisShower->z-Z_TARGET;
449  // double straightPathLength = sqrt(r_shower*r_shower + shower_deltaz*shower_deltaz);
450  // double modulepath = (r_shower-locBCALGeom->GetBCAL_inner_rad())/r_shower*straightPathLength; // project time back to shower location
451  // double t_module = modulepath/SPEED_OF_LIGHT; // project time back to shower location
452  // t_shower += t_module;
453 
454  // printf("shower (%5.1f,%5.1f,%5.1f) r=%5.1f t=%5.1f E=%5.3f rho=%5.1f module rho=%5.1f t=%5.1f new t=%5.1f\n",
455  // shower_x,shower_y,shower_deltaz,r_shower,thisShower->t,E_shower,straightPathLength,modulepath,t_module,t_shower);
456  //int res1 = rt->GetIntersectionWithRadius(r_shower,proj_pos, &pathLength, &flightTime);
457  //int res2 = rt->GetIntersectionWithRadius(locBCALGeom->GetBCAL_inner_rad(),proj_pos, &innerpathLength, &innerflightTime);
458  //if (res1==NOERROR && res2==NOERROR) {
459  DVector3 bcalpos(shower_x,shower_y,Z_shower);
460  double R=bcalpos.Perp();
461  double pathLength=0.;
462  DVector3 proj_mom;
463  vector<DTrackFitter::Extrapolation_t>extrapolations=timeBasedTrack->extrapolations.at(SYS_BCAL);
464  if (fitter->ExtrapolateToRadius(R,extrapolations,proj_pos,proj_mom,
465  flightTime,pathLength)){
466 
467  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 7, "Success profile;Step", 16, -0.5, 15.5);
468  if (thisRFBunch->dNumParticleVotes >= 2){ // Require good RF bunch and this track match the SC
469  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 8, "Success profile;Step", 16, -0.5, 15.5);
470  // We have the flight time to our BCAL point, so we can get the target time
471  double targetCenterTime = t_shower - flightTime - ((timeBasedTrack->position()).Z() - Z_TARGET) / SPEED_OF_LIGHT;
472  sprintf(name, "AllShowers_q%s", q);
473  sprintf(title, "Charged shower; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]");
474  Fill2DHistogram("BCAL_Global_Offsets", "Showers", name,
475  E_shower, targetCenterTime - thisRFBunch->dTime, title,
476  1000, 0.0, 5.0, 200, -10, 10);
477  sprintf(name, "dEdxVsP_q%s", q);
478  sprintf(title, "CDC dE/dx vs P; P [GeV]; dE/dx [keV/cm]");
479  Fill2DHistogram("BCAL_Global_Offsets", "Showers_PID", name,
480  P_track, dEdx, title,
481  200, 0.0, 5.0, 200, 0.0, 5.0);
482  sprintf(name, "deltaTVsdEdx_q%s", q);
483  sprintf(title, "PID; dE/dx [keV/cm]; t_{Target} - t_{RF} [ns]");
484  Fill2DHistogram("BCAL_Global_Offsets", "Showers_PID", name,
485  dEdx, targetCenterTime - thisRFBunch->dTime, title,
486  200, 0.0, 5.0, 200, -10, 10);
487  sprintf(name, "EVsP_q%s", q);
488  sprintf(title, "PID; P [GeV]; E/P");
489  Fill2DHistogram("BCAL_Global_Offsets", "Showers_PID", name,
490  P_track, E_shower, title,
491  200, 0.0, 5.0, 200, 0, 5.0);
492  sprintf(name, "EoverPVsP_q%s", q);
493  sprintf(title, "PID; P [GeV]; E [GeV]");
494  Fill2DHistogram("BCAL_Global_Offsets", "Showers_PID", name,
495  P_track, E_shower/P_track, title,
496  200, 0.0, 5.0, 200, 0, 2);
497  if (dEdx_pion) {
498  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 9, "Success profile;Step", 16, -0.5, 15.5);
499  sprintf(name, "PionShowers_q%s", q);
500  sprintf(title, "Pion showers; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]");
501  Fill2DHistogram("BCAL_Global_Offsets", "Showers", name,
502  E_shower, targetCenterTime - thisRFBunch->dTime, title,
503  1000, 0.0, 5.0, 200, -10, 10);
504  sprintf(name, "PionShowersVsP_q%s", q);
505  sprintf(title, "Pion showers; P [GeV]; t_{Target} - t_{RF} [ns]");
506  Fill2DHistogram("BCAL_Global_Offsets", "Showers", name,
507  P_track, targetCenterTime - thisRFBunch->dTime, title,
508  1000, 0.0, 5.0, 200, -10, 10);
509  sprintf(name, "PionShowersVsZ_q%s", q);
510  sprintf(title, "Pion showers; Z [cm]; t_{Target} - t_{RF} [ns]");
511  Fill2DHistogram("BCAL_Global_Offsets", "Showers", name,
512  Z_shower, targetCenterTime - thisRFBunch->dTime, title,
513  880, 0.0, 420.0, 200, -10, 10);
514  }
515  }
516  }
517 
518 
519  // Get the points from the shower
520  vector <const DBCALPoint*> pointVector;
521  thisShower->Get(pointVector);
522 
523  int N_points = pointVector.size();
524  // N_points vs E_shower
525  sprintf(name, "NpointVsEshower_q%s", q);
526  sprintf(title, "PID; E_{shower} [GeV]; N_{point}");
527  Fill2DHistogram("BCAL_Global_Offsets", "Points", name,
528  E_shower, N_points, title,
529  500, 0.0, 5.0, 50, 0, 50);
530 
531  // Loop over the points within the cluster
532  for (unsigned int iPoint = 0; iPoint < pointVector.size(); iPoint++){
533  const DBCALPoint *thisPoint = pointVector[iPoint];
534  //if (thisPoint->E() < 0.05) continue; // The timing is known not to be great for very low energy, so only use our best info
535  double rpoint = thisPoint->r();
536  float E_point = thisPoint->E();
537  double Z_point = thisPoint->z();
538 
539  float zminhall = 0;
540  float zmaxhall = 450;
541  float zminlocal = -250;
542  float zmaxlocal = 250;
543  Fill2DHistogram ("BCAL_Global_Offsets", "Z Position", "AllPointsVsShower",
544  Z_shower, Z_point + Z_TARGET,
545  "Z_{Point} vs Z_{Shower};Z_{Shower} [cm];Z_{Point} [cm]",
546  225, zminhall, zmaxhall, 225, zminhall, zmaxhall);
547  if (fitter->ExtrapolateToRadius(rpoint,extrapolations,proj_pos,
548  proj_mom,
549  flightTime,pathLength)){
550 
551  // Now proj_pos contains the projected position of the track at this particular point within the BCAL
552  // We can plot the difference of the projected position and the BCAL position as a function of the channel
553  char channame[255], layername[255], chargename[255];
554  sprintf(channame, "M%02iL%iS%i", thisPoint->module(), thisPoint->layer(), thisPoint->sector());
555  sprintf(layername, "AllLayer%i", thisPoint->layer());
556  sprintf(chargename, "All_q%s", q);
557  // These results are in slightly different coordinate systems. We want one where the center of the BCAL is z=0
558  double trackHitZ = proj_pos.z();
559  double localTrackHitZ = proj_pos.z() - locBCALGeom->GetBCAL_center();
560  double BCALHitZ = thisPoint->z() + Z_TARGET;
561  //double localBCALHitZ = thisPoint->z() + Z_TARGET - locBCALGeom->GetBCAL_center();
562  double deltaZ = trackHitZ-BCALHitZ;
563  double Deltat = thisPoint->t_US() - thisPoint->t_DS();
564 
565  Fill2DHistogram ("BCAL_TDC_Offsets", "ZvsDeltat", "AllPoints",
566  Deltat, trackHitZ,
567  "Z_{Track} vs #Delta t;#Delta t = t_{US}-t_{DS};Z_{Track} [cm]",
568  480, -30, 30, 250, zminhall, zmaxhall); // simulation has 16 values in each Deltat=1
569  Fill2DHistogram ("BCAL_TDC_Offsets", "ZvsDeltat", layername,
570  Deltat, trackHitZ,
571  "Z_{Track} vs #Delta t;#Delta t = t_{US}-t_{DS};Z_{Track} [cm]",
572  480, -30, 30, 250, zminhall, zmaxhall);
573  Fill2DHistogram ("BCAL_TDC_Offsets", "ZvsDeltat", chargename,
574  Deltat, trackHitZ,
575  "Z_{Track} vs #Delta t;#Delta t = t_{US}-t_{DS};Z_{Track} [cm]",
576  480, -30, 30, 250, zminhall, zmaxhall);
577  sprintf(title, "%s Z_{Track} vs #Delta t;#Delta t = t_{US}-t_{DS};Z_{Track} [cm]", channame);
578  Fill2DHistogram ("BCAL_TDC_Offsets", "ZvsDeltat", channame,
579  Deltat, trackHitZ, title,
580  480, -30, 30, 250, zminhall, zmaxhall);
581 
582  double trackTheta = 180/3.14159265358*atan2(timeBasedTrack->pperp(),timeBasedTrack->pz());
583  Fill2DHistogram ("BCAL_TDC_Offsets", "DeltatvsTheta", layername,
584  trackTheta, Deltat,
585  "#Delta t vs #theta_{Track};#theta_{Track} (deg);#Delta t = t_{US}-t_{DS} (ns)",
586  360,0,180, 480, -30, 30);
587 
588 
589  int the_cell = (thisPoint->module() - 1) * 16 + (thisPoint->layer() - 1) * 4 + thisPoint->sector();
590  float Deltat_Zcorr = Deltat - (trackHitZ-212)/8.1;
591  sprintf(title, "#Delta t (Hit) corrected for Z;#Delta t - Z_{Track}/v_{eff}");
592  Fill1DHistogram ("BCAL_Global_Offsets", "Deltat", "AllPoints",
593  Deltat_Zcorr, title, 70, -10, 14);
594  Fill2DHistogram ("BCAL_Global_Offsets", "Deltat", "VsCell",
595  the_cell, Deltat_Zcorr,
596  "#Delta t (Hit) corrected for Z;#Delta t - Z_{Track}/v_{eff}",
597  768, 0.5, 768.5, 70, -10, 14);
598 
599  Fill2DHistogram ("BCAL_TDC_Offsets", "Delta Z", "AllPoints",
600  trackHitZ, deltaZ,
601  "#Delta Z vs Z_{Track};Z_{Track} [cm];#Delta Z = Z_{Track} - Z_{Point}",
602  250, zminhall, zmaxhall, 100, -50, 50);
603  Fill2DHistogram ("BCAL_TDC_Offsets", "Delta Z", layername,
604  trackHitZ, deltaZ,
605  "#Delta Z vs Z_{Track};Z_{Track} [cm];#Delta Z = Z_{Track} - Z_{Point}",
606  250, zminhall, zmaxhall, 100, -50, 50);
607  Fill2DHistogram ("BCAL_TDC_Offsets", "Delta Z", chargename,
608  trackHitZ, deltaZ,
609  "#Delta Z vs Z_{Track};Z_{Track} [cm];#Delta Z = Z_{Track} - Z_{Point}",
610  250, zminhall, zmaxhall, 100, -50, 50);
611  Fill2DHistogram ("BCAL_TDC_Offsets", "Delta Z", channame,
612  trackHitZ, deltaZ,
613  "#Delta Z vs Z_{Track};Z_{Track} [cm];#Delta Z = Z_{Track} - Z_{Point}",
614  250, zminhall, zmaxhall, 100, -50, 50);
615 
616  Fill2DHistogram ("BCAL_TDC_Offsets", "Z Position", "AllPoints",
617  trackHitZ, BCALHitZ,
618  "Z_{point} Vs. Z_{Track}; Z_{Track} [cm]; Z_{Point} [cm]",
619  500, zminhall, zmaxhall, 500, zminhall, zmaxhall);
620  Fill2DHistogram ("BCAL_TDC_Offsets", "Z Position", layername,
621  trackHitZ, BCALHitZ,
622  "Z_{point} Vs. Z_{Track}; Z_{Track} [cm]; Z_{Point} [cm]",
623  500, zminhall, zmaxhall, 500, zminhall, zmaxhall);
624  Fill2DHistogram ("BCAL_TDC_Offsets", "Z Position", chargename,
625  trackHitZ, BCALHitZ,
626  "Z_{point} Vs. Z_{Track}; Z_{Track} [cm]; Z_{Point} [cm]",
627  500, zminhall, zmaxhall, 500, zminhall, zmaxhall);
628  Fill2DHistogram ("BCAL_TDC_Offsets", "Z Position", channame,
629  trackHitZ, BCALHitZ,
630  "Z_{point} Vs. Z_{Track}; Z_{Track} [cm]; Z_{Point} [cm]",
631  500, zminhall, zmaxhall, 500, zminhall, zmaxhall);
632 
633  // Get the unifiedhits
634  vector <const DBCALUnifiedHit*> unifiedhitVector;
635  thisPoint->Get(unifiedhitVector);
636  int up = unifiedhitVector[0]->end; // up=1 if unifiedhitVector[0]->end = 1 (is downstream)
637  int down = unifiedhitVector[1]->end; // down=1 if unifiedhitVector[1]->end = 1 (is downstream)
638  const DBCALUnifiedHit *thisUnifiedhitup = unifiedhitVector[up];
639  const DBCALUnifiedHit *thisUnifiedhitdown = unifiedhitVector[down];
640  float t_up = thisUnifiedhitup->t;
641  float t_ADC_up = thisUnifiedhitup->t_ADC;
642  float t_TDC_up = thisUnifiedhitup->t_TDC;
643  float t_down = thisUnifiedhitdown->t;
644  float t_ADC_down = thisUnifiedhitdown->t_ADC;
645  float t_TDC_down = thisUnifiedhitdown->t_TDC;
646  char type[10];
647  sprintf(type,"Mixed");
648  if (t_up == t_ADC_up && t_down == t_ADC_down) sprintf(type,"ADC");
649  if (t_up == t_TDC_up && t_down == t_TDC_down) sprintf(type,"TDC");
650 
651  const DBCALHit * thisADCHit_up;
652  thisUnifiedhitup->GetSingle(thisADCHit_up);
653  const DBCALHit * thisADCHit_down;
654  thisUnifiedhitdown->GetSingle(thisADCHit_down);
655 
656  // Get raw times
657  float Deltat_raw = thisADCHit_up->t_raw - thisADCHit_down->t_raw;
658  float Deltat_raw_Zcorr = Deltat_raw - (trackHitZ-212)/8.1;
659  Fill2DHistogram("BCAL_Global_Offsets", "Deltat_raw", "VsCell",
660  the_cell, Deltat_raw_Zcorr,
661  "#Delta t (Hit) corrected for Z;#Delta t_{raw} - Z_{Track}/v_{eff}",
662  768, 0.5, 768.5, 70, -10, 14);
663  sprintf(title, "#Delta t (Hit) corrected for Z;#Delta t_{raw} - Z_{Track}/v_{eff}");
664  Fill1DHistogram ("BCAL_Global_Offsets", "Deltat_raw", "AllPoints",
665  Deltat_raw_Zcorr, title, 70, -10, 14);
666 
667  // Attenuation Length
668  vector<const DBCALDigiHit*> digihits;
669  thisPoint->Get(digihits);
670  if (digihits.size()!=2) {
671  printf("Warning: BCAL_attenlength_gainratio: event %llu: wrong number of BCALDigiHit objects found %i\n",
672  (long long unsigned int)eventnumber,(int)digihits.size());
673  continue;
674  }
675  if (digihits[0]->end==digihits[1]->end) {
676  printf("Warning: BCAL_attenlength_gainratio: event %llu: two hits in same end of point\n",(long long unsigned int)eventnumber);
677  continue;
678  }
679  float integralUS, integralDS;
680  // end 0=upstream, 1=downstream
681  if (digihits[0]->end==0) {
682  integralUS = digihits[0]->pulse_integral - ((float)digihits[0]->nsamples_integral*(float)digihits[0]->pedestal)/
683  (float)digihits[0]->nsamples_pedestal;
684  integralDS = digihits[1]->pulse_integral - ((float)digihits[1]->nsamples_integral*(float)digihits[1]->pedestal)/
685  (float)digihits[1]->nsamples_pedestal;
686  } else {
687  integralDS = digihits[0]->pulse_integral - ((float)digihits[0]->nsamples_integral*(float)digihits[0]->pedestal)/
688  (float)digihits[0]->nsamples_pedestal;
689  integralUS = digihits[1]->pulse_integral - ((float)digihits[1]->nsamples_integral*(float)digihits[1]->pedestal)/
690  (float)digihits[1]->nsamples_pedestal;
691  }
692  float intratio = (float)integralUS/(float)integralDS;
693  float logintratio = log(intratio);
694 
695  if (VERBOSEHISTOGRAMS) {
696  sprintf(title,"Attenuation;Z_{Track} (cm);log of integral ratio US/DS");
697  Fill2DHistogram ("BCAL_atten_gain", "logintratiovsZtrack", "AllPoints",
698  localTrackHitZ, logintratio, title,
699  250, zminlocal, zmaxlocal, 250, -3, 3);
700  sprintf(title,"Attenuation (M%i,L%i,S%i);Z_{Track} (cm);log of integral ratio US/DS",
701  thisPoint->module(), thisPoint->layer(), thisPoint->sector());
702  Fill2DHistogram ("BCAL_atten_gain", "logintratiovsZtrack", channame,
703  localTrackHitZ, logintratio, title,
704  250, zminlocal, zmaxlocal, 250, -3, 3);
705  }
706 
707  // Now fill some histograms that are useful for aligning the BCAL with the rest of the detector systems
708  if (thisRFBunch->dNumParticleVotes >= 2 && scMatch != NULL && dEdx_pion){ // Require good RF bunch and this track match the SC
709  // Get the time of the BCAL point
710  double pointTime = thisPoint->t();
711  // We have the flight time to our BCAL point, so we can get the target time
712  double vertexTime = (timeBasedTrack->position().Z() - Z_TARGET) / SPEED_OF_LIGHT;; // time for beam to go from center of target to vertex
713  double targetCenterTime = pointTime - flightTime - vertexTime;
714 
715  // Now we just plot the difference in from the RF Time to get out the correction
716  if (E_point > 0.05) { // The timing is known not to be great for very low energy, so only use our best info
717  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell",
718  the_cell, targetCenterTime - thisRFBunch->dTime,
719  "Charged shower points; CCDB Index; t_{Target} - t_{RF} [ns]",
720  768, 0.5, 768.5, 200, -10, 10);
721  sprintf(name , "deltaTVsCell_q%s", q);
722  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", name,
723  the_cell, targetCenterTime - thisRFBunch->dTime,
724  "Charged shower points; CCDB Index; t_{Target} - t_{RF} [ns]",
725  768, 0.5, 768.5, 200, -10, 10);
726  sprintf(name , "deltaTVsCell_q%s_Eweight", q);
727  Fill2DWeightedHistogram("BCAL_Global_Offsets", "Target Time", name,
728  the_cell, targetCenterTime - thisRFBunch->dTime, E_point,
729  "Charged shower points (E weighted); CCDB Index; t_{Target} - t_{RF} [ns]",
730  768, 0.5, 768.5, 200, -10, 10);
731  sprintf(name , "deltaTVsCell_q%s_E2weight", q);
732  Fill2DWeightedHistogram("BCAL_Global_Offsets", "Target Time", name,
733  the_cell, targetCenterTime - thisRFBunch->dTime, E_point*E_point,
734  "Charged shower points (E^{2} weighted); CCDB Index; t_{Target} - t_{RF} [ns]",
735  768, 0.5, 768.5, 200, -10, 10);
736  sprintf(name , "deltaTVsLayer_q%s", q); // for simulation separate by layer
737  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", name,
738  thisPoint->layer(), targetCenterTime - thisRFBunch->dTime,
739  "Charged shower points; CCDB Index; t_{Target} - t_{RF} [ns]",
740  4, 0.5, 4.5, 200, -10, 10);
741  }
742 
743  float pulse_peak_max = max(thisADCHit_up->pulse_peak,thisADCHit_down->pulse_peak);
744  float pulse_peak_min = min(thisADCHit_up->pulse_peak,thisADCHit_down->pulse_peak);
745 
746  double fibLen = locBCALGeom->GetBCAL_length();
747  double c_effective = locBCALGeom->GetBCAL_c_effective();
748  double BCALtrackHitZ = trackHitZ - (locBCALGeom->GetBCAL_center() - fibLen/2); // position wrt BCAL front edge
749  double barproptime_up = BCALtrackHitZ / c_effective;
750  double barproptime_down = (fibLen-BCALtrackHitZ) / c_effective;
751  double hitup_TargetCenterTime = t_up - barproptime_up - flightTime - vertexTime;
752  double hitdown_TargetCenterTime = t_down - barproptime_down - flightTime - vertexTime;
753  double hittimediff = t_down - t_up - 2*localTrackHitZ/c_effective;
754 
755  int the_cell = (thisPoint->module() - 1) * 16 + (thisPoint->layer() - 1) * 4 + thisPoint->sector();
756  //int channel = end*768 + the_cell;
757  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", "hitDeltaTVsChannel",
758  the_cell, hitup_TargetCenterTime - thisRFBunch->dTime,
759  "Charged shower hit; CCDB Index; t_{Target} - t_{RF} [ns]",
760  1536, 0.5, 1536.5, 200, -10, 10);
761  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", "hitDeltaTVsChannel",
762  the_cell+768, hitdown_TargetCenterTime - thisRFBunch->dTime,
763  "Charged shower hit; CCDB Index; t_{Target} - t_{RF} [ns]",
764  1536, 0.5, 1536.5, 200, -10, 10);
765  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", "hittimediff",
766  the_cell, hittimediff,
767  "Charged shower hit; CCDB Index; t_{Target} - t_{RF} [ns]",
768  768, 0.5, 768.5, 200, -10, 10);
769 
770  sprintf(title, "Charged shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
771  sprintf(name, "AllHits_%s_q%s",type,q);
772  Fill2DHistogram("BCAL_Global_Offsets", "Hits_deltaTVsE", name,
773  E_point, targetCenterTime - thisRFBunch->dTime, title,
774  1000, 0.0, 2.0, 200, -10, 10);
775 
776  sprintf(title, "Charged shower points; peak [counts]; t_{Target} - t_{RF} [ns]");
777  sprintf(name, "AllHits_%s_q%s",type,q);
778  Fill2DHistogram("BCAL_Global_Offsets", "Hits_deltaTVsPPmax", name,
779  pulse_peak_max, targetCenterTime - thisRFBunch->dTime, title,
780  1000, 0.0, 4000, 200, -10, 10);
781  sprintf(name, "AllHits_%s_q%s",type,q);
782  Fill2DHistogram("BCAL_Global_Offsets", "Hits_deltaTVsPPmin", name,
783  pulse_peak_min, targetCenterTime - thisRFBunch->dTime, title,
784  1000, 0.0, 4000, 200, -10, 10);
785 
786  int layer = thisPoint->layer();
787  // E_point vs E_shower
788  sprintf(name, "EpointVsEshower_q%s", q);
789  sprintf(title, "PID; E_{shower} [GeV]; E_{point} [GeV]");
790  Fill2DHistogram("BCAL_Global_Offsets", "Points", name,
791  E_shower, E_point, title,
792  1000, 0.0, 5.0, 1000, 0, 2);
793  sprintf(name, "EpointVsEshower_Layer%i_q%s", layer, q);
794  sprintf(title, "PID; E_{shower} [GeV]; E_{point} [GeV]");
795  Fill2DHistogram("BCAL_Global_Offsets", "Points", name,
796  E_shower, E_point, title,
797  1000, 0.0, 5.0, 1000, 0, 2);
798  // deltaT vs E_point
799  sprintf(name, "AllPoints_q%s", q);
800  sprintf(title, "Charged shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
801  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
802  E_point, targetCenterTime - thisRFBunch->dTime, title,
803  1000, 0.0, 2.0, 200, -10, 10);
804  sprintf(name, "Layer%i_q%s", layer, q);
805  sprintf(title, "Charged shower points, layer %i; E_{point} [GeV]; t_{Target} - t_{RF} [ns]", layer);
806  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
807  E_point, targetCenterTime - thisRFBunch->dTime, title,
808  1000, 0.0, 2.0, 200, -10, 10);
809  sprintf(name, "AllPoints_%s_q%s", type, q);
810  sprintf(title, "Charged shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
811  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
812  E_point, targetCenterTime - thisRFBunch->dTime, title,
813  1000, 0.0, 2.0, 200, -10, 10);
814  sprintf(name, "Layer%i_%s_q%s", layer, type, q);
815  sprintf(title, "Charged shower points, layer %i; E_{point} [GeV]; t_{Target} - t_{RF} [ns]", layer);
816  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
817  E_point, targetCenterTime - thisRFBunch->dTime, title,
818  1000, 0.0, 2.0, 200, -10, 10);
819  // deltaT vs E_shower
820  sprintf(name, "AllPoints_q%s", q);
821  sprintf(title, "Charged shower points; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]");
822  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsShowerEnergy", name,
823  E_shower, targetCenterTime - thisRFBunch->dTime, title,
824  1000, 0.0, 2.0, 200, -10, 10);
825  sprintf(name, "Layer%i_q%s", layer, q);
826  sprintf(title, "Charged shower points, layer %i; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]", layer);
827  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsShowerEnergy", name,
828  E_shower, targetCenterTime - thisRFBunch->dTime, title,
829  1000, 0.0, 2.0, 200, -10, 10);
830  }
831  }
832  }
833  }
834 
835  /*************************************************
836  ___ _
837  |\| _ _|_ __ _ | | o __ o __ (_|
838  | |(/_|_| |_ | (_| | | | ||| | | |__|
839 
840  **************************************************/
841 
842 
843  // double pathLength, flightTime;
844  // //double r_shower = sqrt(shower_x*shower_x+shower_y*shower_y);
845  // double t_shower = shower_t;
846  // double E_shower = shower_E;
847  char name[200], title[200];
848  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 10, "Success profile;Step", 16, -0.5, 15.5);
849  if (thisRFBunch->dNumParticleVotes >= 2){ // Require good RF bunch
850  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 11, "Success profile;Step", 16, -0.5, 15.5);
851  vector<const DVertex*> locVertex;
852  loop->Get(locVertex);
853  // *** get unmatched BCAL showers from neutral showers
854  // *** not restrictive enough so do the matching locally
855  //vector <const DNeutralShower *> neutralShower;
856  // loop->Get(neutralShower);
857  // for (unsigned int ishower = 0; ishower < neutralShower.size(); ishower++){
858  // vector <const DBCALShower*> bcalshowervector;
859  // if (neutralShower[ishower]->dDetectorSystem == SYS_BCAL) {
860  // neutralShower[ishower]->Get(bcalshowervector);
861  // }
862  vector <const DBCALShower *> locBCALShowers;
863  loop->Get(locBCALShowers);
864  vector<const DTrackTimeBased*> locTrackTimeBased;
865  loop->Get(locTrackTimeBased);
866  for (unsigned int ibcalshower = 0; ibcalshower < locBCALShowers.size(); ibcalshower++){
867  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 12, "Success profile;Step", 16, -0.5, 15.5);
868  const DBCALShower *bcalshower = locBCALShowers[ibcalshower];
869  double x = bcalshower->x;
870  double y = bcalshower->y;
871  double z = bcalshower->z;
872  DVector3 showerpos(x,y,z);
873  double R_shower = showerpos.Perp();
874  // *** Remove matched BCAL showers
875  bool matched=0;
876  for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i) {
877  DVector3 trackpos(0.0,0.0,0.0);
878  vector<DTrackFitter::Extrapolation_t>extrapolations=locTrackTimeBased[i]->extrapolations.at(SYS_BCAL);
879  if (fitter->ExtrapolateToRadius(R_shower,extrapolations,trackpos)){
880 
881  double dPhi = 180./3.14159265358*(trackpos.Phi()-showerpos.Phi());
882  double dZ = (trackpos.Z() - z);
883  sprintf(name, "Matching");
884  sprintf(title, "Shower-Track position difference;dZ [cm]; d#phi [degrees]");
885  Fill2DHistogram("BCAL_Global_Offsets", "Showers_PID", name,
886  dZ, dPhi, title,
887  200, -60.0, 60.0, 200, -30, 30);
888  // analysis shows 40 and 15 are better
889  if (TMath::Abs(dZ < 40.0) && TMath::Abs(dPhi) < 15) matched=1;
890  }
891  }
892  if (matched) continue;
893  Fill1DHistogram("BCAL_Global_Offsets", "Debug", "Success", 13, "Success profile;Step", 16, -0.5, 15.5);
894 
895  float vertexX = locVertex[0]->dSpacetimeVertex.X();
896  float vertexY = locVertex[0]->dSpacetimeVertex.Y();
897  float vertexZ = locVertex[0]->dSpacetimeVertex.Z();
898  float xdiff = bcalshower->x - vertexX;
899  float ydiff = bcalshower->y - vertexY;
900  float zdiff = bcalshower->z - vertexZ;
901  float pathLength = sqrt(xdiff*xdiff + ydiff*ydiff + zdiff*zdiff);
902  // UNPROJECTION CODE
903  // since time is projected to BCAL inner radius, shorten pathlength to match
904  //pathLength *= locBCALGeom->GetBCAL_inner_rad()/R_shower;
905 
906  float flightTime = pathLength/SPEED_OF_LIGHT;
907  float vertexTime = (vertexZ - Z_TARGET) / SPEED_OF_LIGHT; // time for beam to go from center of target to vertex
908  double targetCenterTime = bcalshower->t - flightTime - vertexTime;
909  double deltaTime = targetCenterTime - thisRFBunch->dTime;
910  double E_shower = bcalshower->E;
911  //printf("shower (%5.1f,%5.1f,%5.1f) vertex (%5.1f,%5.1f,%5.1f) length (%5.1f,%5.1f,%5.1f) time sh %5.1f flight %5.1f %5.1f target %5.1f rf %5.1f\n",
912  // bcalshower->x,bcalshower->y,bcalshower->z,vertexX,vertexY,vertexZ,
913  // xdiff,ydiff,zdiff, bcalshower->t,flightTime,vertexTime,targetCenterTime,thisRFBunch->dTime);
914  sprintf(name, "AllShowers_q0");
915  sprintf(title, "Neutral Showers; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]");
916  Fill2DHistogram("BCAL_Global_Offsets", "Showers", name,
917  E_shower, deltaTime, title,
918  1000, 0.0, 5.0, 200, -10, 10);
919  sprintf(name, "AllShowersVsZ_q0");
920  sprintf(title, "Neutral showers; Z [cm]; t_{Target} - t_{RF} [ns]");
921  Fill2DHistogram("BCAL_Global_Offsets", "Showers", name,
922  bcalshower->z, deltaTime, title,
923  880, 0.0, 440.0, 200, -10, 10);
924 
925  // // Get the points from the shower
926  vector <const DBCALPoint*> pointVector;
927  bcalshower->Get(pointVector);
928  int N_points = pointVector.size();
929  // N_points vs E_shower
930  sprintf(name, "NpointVsEshower_q0");
931  sprintf(title, "PID; E_{shower} [GeV]; N_{point}");
932  Fill2DHistogram("BCAL_Global_Offsets", "Points", name,
933  E_shower, N_points, title,
934  500, 0.0, 5.0, 50, 0, 50);
935  // Loop over the points within the cluster
936  for (unsigned int iPoint = 0; iPoint < pointVector.size(); iPoint++){
937  const DBCALPoint *thisPoint = pointVector[iPoint];
938  float E_point = thisPoint->E();
939  // Now we just plot the difference in from the RF Time to get out the correction
940  pathLength = thisPoint->rho(); // This is an approximation: the photon doesn't come from the target center
941  flightTime = pathLength/SPEED_OF_LIGHT;
942  vertexTime = 0; // no vertex time since we are using target center from point
943  targetCenterTime = thisPoint->t() - flightTime - vertexTime;
944  deltaTime = targetCenterTime - thisRFBunch->dTime;
945 
946  if (E_point > 0.05) { // The timing is known not to be great for very low energy, so only use our best info
947  int the_cell = (thisPoint->module() - 1) * 16 + (thisPoint->layer() - 1) * 4 + thisPoint->sector();
948  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell_q0",
949  the_cell, deltaTime,
950  "Neutral shower points; CCDB Index; t_{Target} - t_{RF} [ns]",
951  768, 0.5, 768.5, 200, -10, 10);
952  Fill2DHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell_q0",
953  the_cell, deltaTime,
954  "Neutral shower points; CCDB Index; t_{Target} - t_{RF} [ns]",
955  768, 0.5, 768.5, 200, -10, 10);
956  Fill2DWeightedHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell_q0_Eweight",
957  the_cell, deltaTime, E_point,
958  "Neutral shower points (E weighted); CCDB Index; t_{Target} - t_{RF} [ns]",
959  768, 0.5, 768.5, 200, -10, 10);
960  Fill2DWeightedHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell_q0_E2weight",
961  the_cell, deltaTime, E_point*E_point,
962  "Neutral shower points (E^{2} weighted); CCDB Index; t_{Target} - t_{RF} [ns]",
963  768, 0.5, 768.5, 200, -10, 10);
964  }
965  int layer = thisPoint->layer();
966  // deltaT vs E_point
967  sprintf(name, "AllPoints_q0");
968  sprintf(title, "Neutral shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
969  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
970  E_point, deltaTime, title,
971  1000, 0.0, 2.0, 200, -10, 10);
972  sprintf(name, "Layer%i_q0", layer);
973  sprintf(title, "Neutral shower points, layer %i; E_{point} [GeV]; t_{Target} - t_{RF} [ns]", layer);
974  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
975  E_point, deltaTime, title,
976  1000, 0.0, 2.0, 200, -10, 10);
977  // sprintf(name, "AllPoints_%s_q0", type);
978  // sprintf(title, "Neutral shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
979  // Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
980  // E_point, deltaTime, title,
981  // 1000, 0.0, 2.0, 200, -10, 10);
982  // sprintf(name, "Layer%i_%s_q0", layer, type);
983  // sprintf(title, "Neutral shower points, layer %i; E_{point} [GeV]; t_{Target} - t_{RF} [ns]", layer);
984  // Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsEnergy", name,
985  // E_point, deltaTime, title,
986  // 1000, 0.0, 2.0, 200, -10, 10);
987  // deltaT vs E_shower
988  sprintf(name, "AllPoints_q0");
989  sprintf(title, "Neutral shower points; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]");
990  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsShowerEnergy", name,
991  E_shower, deltaTime, title,
992  1000, 0.0, 2.0, 200, -10, 10);
993  sprintf(name, "Layer%i_q0", layer);
994  sprintf(title, "Neutral shower points, layer %i; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]", layer);
995  Fill2DHistogram("BCAL_Global_Offsets", "Points_deltaTVsShowerEnergy", name,
996  E_shower, deltaTime, title,
997  1000, 0.0, 2.0, 200, -10, 10);
998 
999  // Now simulate the projection to the inner radius
1000 
1001  pathLength = thisPoint->rho()*locBCALGeom->GetBCAL_inner_rad()/R_shower;
1002  flightTime = pathLength/SPEED_OF_LIGHT;
1003  vertexTime = 0; // no vertex time since we are using target center from point
1004  targetCenterTime = thisPoint->tInnerRadius() - flightTime - vertexTime;
1005  float altDeltaTime = targetCenterTime - thisRFBunch->dTime;
1006 
1007  // altDeltaT vs E_point
1008  sprintf(name, "AllPoints_q0");
1009  sprintf(title, "Neutral shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
1010  Fill2DHistogram("BCAL_Global_Offsets", "Points_altDeltaTVsEnergy", name,
1011  E_point, altDeltaTime, title,
1012  1000, 0.0, 2.0, 200, -10, 10);
1013  sprintf(name, "Layer%i_q0", layer);
1014  sprintf(title, "Neutral shower points, layer %i; E_{point} [GeV]; t_{Target} - t_{RF} [ns]", layer);
1015  Fill2DHistogram("BCAL_Global_Offsets", "Points_altDeltaTVsEnergy", name,
1016  E_point, altDeltaTime, title,
1017  1000, 0.0, 2.0, 200, -10, 10);
1018  // sprintf(name, "AllPoints_%s_q0", type);
1019  // sprintf(title, "Neutral shower points; E_{point} [GeV]; t_{Target} - t_{RF} [ns]");
1020  // Fill2DHistogram("BCAL_Global_Offsets", "Points_altDeltaTVsEnergy", name,
1021  // E_point, altDeltaTime, title,
1022  // 1000, 0.0, 2.0, 200, -10, 10);
1023  // sprintf(name, "Layer%i_%s_q0", layer, type);
1024  // sprintf(title, "Neutral shower points, layer %i; E_{point} [GeV]; t_{Target} - t_{RF} [ns]", layer);
1025  // Fill2DHistogram("BCAL_Global_Offsets", "Points_altDeltaTVsEnergy", name,
1026  // E_point, altDeltaTime, title,
1027  // 1000, 0.0, 2.0, 200, -10, 10);
1028  // altDeltaT vs E_shower
1029  sprintf(name, "AllPoints_q0");
1030  sprintf(title, "Neutral shower points; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]");
1031  Fill2DHistogram("BCAL_Global_Offsets", "Points_altDeltaTVsShowerEnergy", name,
1032  E_shower, altDeltaTime, title,
1033  1000, 0.0, 2.0, 200, -10, 10);
1034  sprintf(name, "Layer%i_q0", layer);
1035  sprintf(title, "Neutral shower points, layer %i; E_{shower} [GeV]; t_{Target} - t_{RF} [ns]", layer);
1036  Fill2DHistogram("BCAL_Global_Offsets", "Points_altDeltaTVsShowerEnergy", name,
1037  E_shower, altDeltaTime, title,
1038  1000, 0.0, 2.0, 200, -10, 10);
1039  }
1040  }
1041  }
1042 
1043  return NOERROR;
1044 }
1045 
1046 //------------------
1047 // erun
1048 //------------------
1050 {
1051  // This is called whenever the run number changes, before it is
1052  // changed to give you a chance to clean up before processing
1053  // events from the next run number.
1054  return NOERROR;
1055 }
1056 
1057 //------------------
1058 // fini
1059 //------------------
1061 {
1062  // Called before program exit after event processing is finished.
1063  SortDirectories(); // Sort the histogram directories by name
1064  return NOERROR;
1065 }
1066 
float rho() const
Definition: DBCALPoint.h:50
bool ExtrapolateToRadius(double R, const vector< Extrapolation_t > &extraps, DVector3 &pos, DVector3 &mom, double &t, double &s) const
#define SPEED_OF_LIGHT
int module() const
Definition: DBCALPoint.h:64
float t_TDC
Time of TDC hit that is closes t to the ADC time.
float tInnerRadius() const
Definition: DBCALPoint.cc:143
Int_t layer
TVector3 DVector3
Definition: DVector3.h:14
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
uint32_t fp_trig_mask
Definition: DL1Trigger.h:19
double pz(void) const
sprintf(text,"Post KinFit Cut")
#define y
void Fill1DWeightedHistogram(const char *plugin, const char *directoryName, const char *name, const double value, const double weight, const char *title, int nBins, double xmin, double xmax, bool print=false)
int sector() const
Definition: DBCALPoint.h:66
const DVector3 & position(void) const
const DTrackTimeBased * Get_TrackTimeBased(void) const
trig[33-1]
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
jerror_t init(void)
Called once at program start.
int layer() const
Definition: DBCALPoint.h:65
Definition: GlueX.h:19
JApplication * japp
void Fill2DWeightedHistogram(const char *plugin, const char *directoryName, const char *name, const double valueX, const double valueY, const double weight, const char *title, int nBinsX, double xmin, double xmax, int nBinsY, double ymin, double ymax, bool print=false)
double counter
Definition: FitGains.C:151
void SortDirectories()
float E() const
Definition: DBCALPoint.h:38
InitPlugin_t InitPlugin
const bool VERBOSE
DGeometry * GetDGeometry(unsigned int run_number)
#define _DBG_
Definition: HDEVIO.h:12
float t_US() const
Return the time of US Hit.
Definition: DBCALPoint.h:42
float t_ADC
Time from fADC.
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
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
int Ndof
Number of degrees of freedom in the fit.
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 erun(void)
Called everytime run number changes, provided brun has been called.
float t
Definition: DBCALHit.h:31
shared_ptr< const DSCHitMatchParams > Get_SCHitMatchParams(void) const
#define Z_TARGET
Definition: DRiemannFit.cc:12
double sqrt(double)
const DVector3 & momentum(void) const
int pulse_peak
Definition: DBCALHit.h:29
const DTrackFitter * fitter
shared_ptr< const DBCALShowerMatchParams > Get_BCALShowerMatchParams(void) const
float t() const
Definition: DBCALPoint.h:41
float r() const
Definition: DBCALPoint.h:62
map< DetectorSystem_t, vector< DTrackFitter::Extrapolation_t > > extrapolations
unsigned int dNumParticleVotes
Definition: DEventRFBunch.h:32
printf("string=%s", string)
double pperp(void) const
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
float t_DS() const
Return the time of DS Hit.
Definition: DBCALPoint.h:43
jerror_t fini(void)
Called after last event of last event source has been processed.
float t_raw
Uncalibrated time in ns.
Definition: DBCALHit.h:32
float t
Unified time, obtained from ADC and/or TDC and used for further analysis.