UpMuTrigger_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: UpMuTrigger
3 // Module Type: filter
4 // File: UpMuTrigger_module.cc
5 //
6 ////////////////////////////////////////////////////////////////////////
7 
13 #include "fhiclcpp/ParameterSet.h"
21 #include "TTimeStamp.h"
22 #include "TFile.h"
23 #include "TTree.h"
24 
36 
37 #include <TGraphErrors.h>
38 #include <TPaveText.h>
39 #include <TCanvas.h>
40 #include <TF1.h>
41 #include <TH1F.h>
42 #include <TApplication.h>
43 #include <TNtuple.h>
44 #include <TMath.h>
45 #include <string>
46 #include <vector>
47 #include <algorithm>
48 #include <numeric>
49 #include <TFile.h>
50 
51 #include <boost/math/special_functions/gamma.hpp>
52 #include <boost/math/special_functions.hpp>
53 
54 namespace novaddt {
55  class UpMuTrigger;
56 }
57 
59 public:
60  explicit UpMuTrigger(fhicl::ParameterSet const & p);
61  virtual ~UpMuTrigger();
62  void endJob() override;
63  void beginJob() override;
64  virtual bool filter(art::Event & e) override;
65 
66 private:
67 /*
68  art::ServiceHandle<art::TFileService> tfs;
69  TNtuple* fupmu;
70  TH1F *fChi2;
71  TH1F *fLLR;
72 */
73 
74  unsigned _prescale;
80 
83 
85 
86  double c0x, c0y, p0, cw, pw; // initialized later based on detector
87  double _TrackLen;
88  unsigned _TrackHitsXY;
89  unsigned _TrackHitsX;
90  unsigned _TrackHitsY;
91  int _dX;
92  int _dY;
93  int _dZ;
94  double _Chi2;
95  double _LLR;
96  double _R2X;
97  double _R2Y;
98  double _MinSlope;
99  double _MaxSlope;
100 
101  //Float_t Chi2_plot;
102  //Float_t LLR_plot;
103 
104  double GetT(double tdc, double dist_to_readout);
105  double GetRes(int ADC);
106  double GetRes2(int ADC);
107  void GetXYZ(DAQHit const& daqhit, TVector3 start, TVector3 end, double *xyzd);
108 
109  void LinFit(const std::vector<double>& x, const std::vector<double>& y, double *fitpar);
110  void LinFit(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& ye, double *fitpar);
111  void LinFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error,
112  double& slope, double& chi2, double& P_up, double& P_dn);
113 
114  double GetPigtail(DAQHit const& daqhit) const;
115  double getDCMOff(uint16_t cell, uint8_t view);
116 };
117 
119  : _prescale (p.get<unsigned>("prescale"))
120  , _slicelabel (p.get<std::string>("slice_label"))
121  , _sliceinstance (p.get<std::string>("slice_instance"))
122  , _tracklabel (p.get<std::string>("track_label"))
123  , _trackinstance (p.get<std::string>("track_instance"))
124  , _trackToSlice (p.get<std::string>("TrackToSliceInstanceLabel"))
125  , _detector (p.get<std::string>("detector"))
126  , _containedbit (p.get<bool>("containedbit"))
127 
128  , _TrackLen (p.get<double>("TrackLen"))
129  , _TrackHitsXY (p.get<unsigned>("TrackHitsXY"))
130  , _TrackHitsX (p.get<unsigned>("TrackHitsX"))
131  , _TrackHitsY (p.get<unsigned>("TrackHitsY"))
132  , _dX (p.get<int>("dX"))
133  , _dY (p.get<int>("dY"))
134  , _dZ (p.get<int>("dZ"))
135  , _Chi2 (p.get<double>("Chi2"))
136  , _LLR (p.get<double>("LLR"))
137  , _R2X (p.get<double>("R2X"))
138  , _R2Y (p.get<double>("R2Y"))
139  , _MinSlope (p.get<double>("MinSlope"))
140  , _MaxSlope (p.get<double>("MaxSlope"))
141 {
142  assert(_detector == "ndos" || _detector == "fd");
143 
144  //initialize cell and plane locations
145  if (_detector == "ndos") {
146  c0x = -119.5;
147  c0y = -199.5;
148  p0 = 4.57;
149  cw = 3.97;
150  pw = 6.654;
151  } else if (_detector == "fd") {
152  c0x = -759.5;
153  c0y = -759.5;
154  p0 = 4.57;
155  cw = 3.97;
156  pw = 6.654;
157  }
158 
159  _trigger_counts=0;
160  _after_prescale=0;
161 
162  produces<std::vector<TriggerDecision>>();
163 }
164 
166 {
167  // Clean up dynamic memory and other resources here.
168 }
169 
171 {
172  std::unique_ptr<std::vector<TriggerDecision>>
173  trigger_decisions(new std::vector<TriggerDecision>);
174  bool result = false;
175 
177  e.getByLabel(_tracklabel, tracks);
178 
179  // get the hit list for each track. Time cosmusing.
181  assert(fohl.isValid());
182  assert(fohl.size() == tracks->size());
183 
184  for(size_t i_track = 0; i_track < tracks->size(); ++i_track){
185  novaddt::Track3D track = tracks->at(i_track);
186  // get the hit list for this track
187  art::Ptr<novaddt::HitList> this_hit_list = fohl.at(i_track);
188 
189  // use these to define the trigger window
190  // safest not to assume this is time sorted
191  uint64_t min_tdc = this_hit_list->front().TDC().val;
192  uint64_t max_tdc = this_hit_list->back().TDC().val;
193  for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
194  if (this_hit_list->at(i_hit).TDC().val < min_tdc) min_tdc = this_hit_list->at(i_hit).TDC().val;
195  if (this_hit_list->at(i_hit).TDC().val > max_tdc) max_tdc = this_hit_list->at(i_hit).TDC().val;
196  } // end of loop on hits
197 
198  novaddt::TDC absolute_t0 = min_tdc;
199 
200  int StartX = track.Start().X();
201  int StartY = track.Start().Y();
202  int StartZ = track.Start().Z();
203 
204  int EndX = track.End().X();
205  int EndY = track.End().Y();
206  int EndZ = track.End().Z();
207  if(StartY > EndY){
208  std::swap(StartX, EndX);
209  std::swap(StartY, EndY);
210  std::swap(StartZ, EndZ);
211  }
212 
213  TVector3 length;
214  length.SetXYZ(cw*(track.End().X() - track.Start().X()),
215  cw*(track.End().Y() - track.Start().Y()),
216  pw*(track.End().Z() - track.Start().Z())
217  );
218 
219  double tr_length = length.Mag();
220  double exp_trav_time = tr_length/29.97; // c=29.97 cm/ns; 1/64E6=15.625 ns
221 
222  if (tr_length < _TrackLen) continue; // Select tracks longer than _TrackLen.
223  if (TMath::Abs(EndX - StartX) < _dX) continue; // Select tracks longer than _dX cells in X.
224  if (TMath::Abs(EndY - StartY) < _dY) continue; // Select tracks longer than _dY cells in Y.
225  if (TMath::Abs(EndZ - StartZ) < _dZ) continue; // Select tracks longer than _dZ planes in Z.
226  if (!track.Is3D()) continue; // Select only 3D tracks
227 
228  int T0 = 0;
229  std::vector<double> x_hit;
230  std::vector<double> y_hit;
231  std::vector<double> zy_hit;
232  std::vector<double> zx_hit;
233 
234  // Get T0 that corresponds to the minimum cell number in Y-view
235  for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
236  uint16_t iC = this_hit_list->at(i_hit).Cell().val;
237  uint16_t iP = this_hit_list->at(i_hit).Plane().val;
238  uint8_t view = this_hit_list->at(i_hit).View().val;
239  int iTDC = static_cast<int>(this_hit_list->at(i_hit).TDC().val - absolute_t0.val);
240 
241  if(view == daqchannelmap::X_VIEW){
242  x_hit.push_back(iC);
243  zx_hit.push_back(iP);
244  }
245 
246  if(view == daqchannelmap::Y_VIEW){
247  if(iC == StartY)
248  T0 = iTDC;
249  y_hit.push_back(iC);
250  zy_hit.push_back(iP);
251  }
252 
253  }
254 
255  // Linear fit of the hits in XZ/YZ views. Drop tracks that has bad linear fit
256 
257 
258  double fitpar[3];
259  LinFit(x_hit, zx_hit, fitpar);
260  double r2x = fitpar[2];
261  LinFit(y_hit, zy_hit, fitpar);
262  double r2y = fitpar[2];
263 
264  if(r2x < _R2X) continue;
265  if(r2y < _R2Y) continue;
266 
267  // Vectors of expected and measured times for TGraph
268  std::vector<double> eT_vecxy;
269  std::vector<double> mT_vecxy;
270  std::vector<double> mT_vecx;
271  std::vector<double> mT_vecy;
272 
273  // Time errors
274  std::vector<double> mT_vecxy_e;
275 
276  for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
277  uint16_t iC = this_hit_list->at(i_hit).Cell().val;
278  uint16_t iP = this_hit_list->at(i_hit).Plane().val;
279  int hit_adc = this_hit_list->at(i_hit).ADC().val;
280 
281  if (hit_adc < 50) continue;
282 
283  // Calculate time of the hit
284  DAQHit const& dqhit = this_hit_list->at(i_hit);
285  double pos[4];
286  TVector3 start = track.Start();
287  TVector3 end = track.End();
288  GetXYZ(dqhit, start, end, pos);
289  double pigtail = GetPigtail(dqhit);
290  int iTDC = static_cast<int>(this_hit_list->at(i_hit).TDC().val - absolute_t0.val);
291  double meas_time = static_cast<double>(iTDC - T0);
292  meas_time = meas_time + 0.01*static_cast<double>((this_hit_list->at(i_hit).TDC().frac));
293 
294  double tns = GetT(meas_time, pos[3] + pigtail);
295  // Subtract DCM offsects
296  double DCM_off = getDCMOff(iC, this_hit_list->at(i_hit).View().val);
297  tns = tns - DCM_off;
298 
299  if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW){
300  if(iP > TMath::Max(StartZ, EndZ)) continue;
301  if(iP < TMath::Min(StartZ, EndZ)) continue;
302 
303  if(iC > TMath::Max(StartX, EndX)) continue;
304  if(iC < TMath::Min(StartX, EndX)) continue;
305  }
306  else{
307  if(iP > TMath::Max(StartZ, EndZ)) continue;
308  if(iP < TMath::Min(StartZ, EndZ)) continue;
309 
310  if(iC > TMath::Max(StartY, EndY)) continue;
311  if(iC < TMath::Min(StartY, EndY)) continue;
312  }
313 
314  double exp_time = -999;
315  std::string h_view = "y";
316 
317  if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW){
318  h_view = "x";
319  exp_time = exp_trav_time*(iC-StartX)/(EndX-StartX);
320  }
321  else
322  exp_time = exp_trav_time*(iC-StartY)/(EndY-StartY);
323 
324  if(false && tr_length > 1200. && tns<-1000){
325 
326  std::cout << ", i_track: " << std::setw(3) << i_track
327  << ", i_hit: " << std::setw(3) << i_hit
328  << ", iP: " << std::setw(3) << iP
329  << ", iC: " << std::setw(3) << iC
330  << ", T0: " << std::setw(3) << T0
331  << ", pigtail: " << std::setw(3) << pigtail
332  << ", DCM_off: " << std::setw(3) << DCM_off
333  << ", ADC: " << std::setw(5) << (int)this_hit_list->at(i_hit).ADC().val
334  << ", TDC: " << (uint64_t)this_hit_list->at(i_hit).TDC().val
335  << ", meas_time: " << meas_time
336  << ", tns: " << tns
337  << ", exp_time: " << exp_time
338  << std::endl;
339  }
340 
341  assert(pos[3]>=0);
342 
343  eT_vecxy.push_back(exp_time);
344  mT_vecxy.push_back(tns);
345  mT_vecxy_e.push_back(GetRes(hit_adc)); // mT_vecxy_e.push_back(GetRes2(hit_adc));
346  if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW)
347  mT_vecx.push_back(tns);
348  if( this_hit_list->at(i_hit).View().val == daqchannelmap::Y_VIEW)
349  mT_vecy.push_back(tns);
350 
351  } // end of loop on hits
352 
353  unsigned nxhit = mT_vecx.size();
354  unsigned nyhit = mT_vecy.size();
355 
356  // Remove tracks with less than fTrackHits hits in X+Y, X, Y views
357  if(nxhit < _TrackHitsX) continue;
358  if(nyhit < _TrackHitsY) continue;
359  if((nxhit+nyhit) < _TrackHitsXY) continue;
360 
361  double slope_xy, chi_xy, P_up_xy, P_dn_xy;
362  LinFitLLR(eT_vecxy, mT_vecxy, mT_vecxy_e, slope_xy, chi_xy, P_up_xy, P_dn_xy);
363  double LLR = log(P_up_xy/P_dn_xy);
364 
365 
366  // Remove tracks with LLR smaller than input in the config fcl file
367  if(LLR < _LLR) continue;
368  // Remove tracks with Chi2 fit values higher than input in the config fcl file
369  if(chi_xy > _Chi2) continue;
370 
371  // Remove tracks with Slope fit values lower/higher than input in the config fcl file
372  if(slope_xy < _MinSlope) continue;
373  if(slope_xy > _MaxSlope) continue;
374 
375  if(false){
376  std::cout<< std::setiosflags(std::ios::fixed)
377  << "Run: " << std::setprecision(4) << std::setw(8) << std::left << e.id().run()
378  << ", Event: " << std::setprecision(4) << std::setw(8) << std::left << e.id().event()
379  << ", LLR: " << std::setprecision(4) << std::setw(5) << std::left << LLR
380  << ", Chi2: " << std::setprecision(4) << std::setw(5) << std::left << chi_xy
381  << ", slope: " << std::setprecision(4) << std::setw(5) << std::left << slope_xy
382  << ", tr_length: " << std::setprecision(3) << std::setw(3) << std::left << tr_length
383  << ", dX: " << std::setprecision(4) << std::setw(4) << std::left << TMath::Abs(EndX - StartX)
384  << ", dY: " << std::setprecision(4) << std::setw(4) << std::left << TMath::Abs(EndY - StartY)
385  << ", dZ: " << std::setprecision(4) << std::setw(4) << std::left << TMath::Abs(EndZ - StartZ)
386  << ", 3d: " << std::setprecision(1) << std::setw(1) << std::left << track.Is3D()
387  << ", nxhit: " << std::setprecision(4) << std::setw(3) << std::left << nxhit
388  << ", nyhit: " << std::setprecision(4) << std::setw(3) << std::left << nyhit
389  << ", nxyhit: " << std::setprecision(4) << std::setw(3) << std::left << nxhit+nyhit
390  << ", r2x: " << std::setprecision(4) << std::setw(5) << std::left << r2x
391  << ", r2y: " << std::setprecision(4) << std::setw(5) << std::left << r2y
392  << std::endl;
393  }
394 
395 
396  _trigger_counts++;
397  //fChi2->Fill(chi_xy);
398  //fLLR->Fill(LLR);
399 
400  //float track_entries[2]=
401  // { (float)LLR,
402  // (float)chi_xy };
403  //fupmu->Fill(track_entries);
404 
406  _after_prescale++;
407 
408  if(_containedbit)
409  trigger_decisions->emplace_back
410  (min_tdc, max_tdc - min_tdc, daqdataformats::trigID::TRIG_ID_DATA_CONTAINED, _prescale);
411  else
412  trigger_decisions->emplace_back
413  (min_tdc, max_tdc - min_tdc, daqdataformats::trigID::TRIG_ID_DATA_UPMU, _prescale);
414 
415  result = true;
416  break; // don't need to loop any more on this slice
417  } // end if (prescale)
418 
419  } //end loop over tracks
420 
421  e.put(std::move(trigger_decisions));
422  return result;
423 
424 } // end Filter()
425 //----------------------------------------------------------------------------------
427 {
428 
429  //fupmu = tfs->make<TNtuple>("upmu_ntuple", "Track Ntuple","LLR:Chi2");
430  //fChi2 = tfs->make<TH1F>("chi_xy", "Chi2; Chi2; ", 200,-1,2);
431  //fLLR = tfs->make<TH1F>("LLR", "LLR ; LLR ; ", 200,0,80);
432  }
433 //----------------------------------------------------------------------------------
435 {
436  std::cout << "=== novaddt::UpmuTrigger endJob" << std::endl;
437 }
438 //----------------------------------------------------------------------------------
439 void novaddt::UpMuTrigger::GetXYZ(DAQHit const& daqhit, TVector3 start, TVector3 end, double *xyzd) {
440  UInt_t plane = daqhit.Plane().val;
441  UInt_t cell = daqhit.Cell().val;
442 
443  // The parameters are taken from the fit of plane vs Z
444  double x = c0x + cw*cell;
445  double y = x;
446  double z = p0 + pw*plane;
447  double d = -999.9;
448 
449  double x0 = c0x + cw*start.X(); double x1 = c0x + cw*end.X();
450  double y0 = c0y + cw*start.Y(); double y1 = c0y + cw*end.Y();
451  double z0 = p0 + pw*start.Z(); double z1 = p0 + pw*end.Z();
452 
453  if(daqhit.View().val == daqchannelmap::X_VIEW){
454  y = (y1 - y0)/(z1 - z0);
455  y = y*(z - z0);
456  y = y + y0;
457  d = 800 - y;
458  }
459  if(daqhit.View().val == daqchannelmap::Y_VIEW){
460  x = (x1 - x0)/(z1 - z0);
461  x = x*(z - z0);
462  x = x + x0;
463  d = 800 - x;
464  }
465 
466  if(d<0) d=0;
467  xyzd[0] = x;
468  xyzd[1] = y;
469  xyzd[2] = z;
470  xyzd[3] = d;
471 }
472 
473 double novaddt::UpMuTrigger::GetT(double tdc, double dist_to_readout) {
474  // This needs to be calibrated and put in a DB table
475  double speedOfFiberTransport = 15.3; // cm/ns, "first principles" calib.
476  double TDC_to_ns = 15.625; // 1/64E6=15.625 ns
477 
478  // Differs from c/n due to zigzag
479  // paths down fiber. But, this
480  // value may be way off (10%? 30%?).
481  return (tdc*TDC_to_ns - dist_to_readout/speedOfFiberTransport);
482 }
483 
484 
485 
487 
488  //The fit function from Evan: f(x) = p0/(p1+x^p2) + p3;
489  double p0 = 70267.4;
490  double p1 = 677.305;
491  double p2 = 1.48702;
492  double p3 = 8.46997;
493 
494  double res = p0/(p1 + pow( static_cast<double>(ADC), p2) ) + p3;
495  return res;
496 
497 }
498 
500 {
501  //The fit function from Evan: f(x) = p0/(p1+x^p2) + p3;
502 
503  //1.273e+06
504  double p0 = 1273000;
505 
506  //1.781e+04
507  double p1 = 17810;
508 
509  //
510  double p2 = 2.09;
511  double p3 = 8.246;
512 
513  double res = p0/(p1 + pow( static_cast<double>(ADC), p2) ) + p3;
514  return res;
515 }
516 
517 double novaddt::UpMuTrigger::GetPigtail(DAQHit const& daqhit) const{
518  if (_detector == "fd") {
519  UInt_t cellid = daqhit.Cell().val;
520  //NOTE: Do we think 32 cells per module can ever change
521  //This number can be found in DAQChannelMap but did not want to add that
522  //dependancy to geometry unless needed.
523  const int kCellsPerModule = 32;
524  cellid = cellid % kCellsPerModule;
525  // In vertical planes, high cell numbers have longer fibres.
526  // In horizontal planes it's the opposite, so flip it round.
527  if(daqhit.View().val == daqchannelmap::X_VIEW)
528  cellid = kCellsPerModule-cellid-1;
529 
530  // This should really never happen, but just to be safe...
531  if(cellid >= kCellsPerModule) return 100;
532  // Email from Tom Chase 2011-04-29
533  // NB: this isn't just a ~3.8cm change per cell. Every 8 cells something
534  // different happens.
535  const double kPigtails[kCellsPerModule] = {
536  34.5738, 38.4379, 42.3020, 46.1660, 50.0301, 53.8942, 57.7582, 61.6223,
537  64.7504, 68.6144, 72.4785, 76.3426, 80.2067, 84.0707, 87.9348, 91.0790,
538  95.3301, 99.1941, 103.058, 106.922, 110.786, 114.650, 118.514, 122.379,
539  125.507, 129.371, 133.235, 137.099, 140.963, 144.827, 148.691, 150.751
540  };
541  return kPigtails[cellid];
542  }
543  else
544  return 0; //to be updated as necessary
545 }
546 
547 
548 void novaddt::UpMuTrigger::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
549  // http://en.wikipedia.org/wiki/Simple_linear_regression
550  const auto n = x_val.size();
551  const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n; // <x>
552  const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n; // <y>
553  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n; // <xx>
554  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n; // <yy>
555  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n; // <xy>
556 
557  const auto b = (xy - x*y)/(xx - x*x); // slope
558  const auto a = y - b*x; // intercept
559  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
560  fitpar[0] = a;
561  fitpar[1] = b;
562  fitpar[2] = r*r;
563 
564 }
565 
566 void novaddt::UpMuTrigger::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, const std::vector<double>& y_err, double *fitpar){
567  // http://en.wikipedia.org/wiki/Simple_linear_regression
568  int n = x_val.size();
569  double x = 0;
570  double y = 0;
571  double xx = 0;
572  double yy = 0;
573  double xy = 0;
574  double ee = 0;
575 
576  for ( int i=0; i<n; i++ ){
577 
578  x = x + x_val[i]/y_err[i]/y_err[i];
579  y = y + y_val[i]/y_err[i]/y_err[i];
580 
581  xx = xx + x_val[i]*x_val[i]/y_err[i]/y_err[i];
582  yy = yy + y_val[i]*y_val[i]/y_err[i]/y_err[i];
583  xy = xy + x_val[i]*y_val[i]/y_err[i]/y_err[i];
584  ee = ee + 1./y_err[i]/y_err[i];
585  }
586 
587  const auto b = (xy*ee - x*y)/(xx*ee - x*x); // slope
588  const auto a = (xy - b*xx)/x; // intercept
589  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
590  fitpar[0] = a;
591  fitpar[1] = b;
592  fitpar[2] = r*r;
593 
594 
595 }
596 
597 void novaddt::UpMuTrigger::LinFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error,
598  double& slope, double& chi2, double& P_up, double& P_dn) {
599 
600  int n = x_hit.size();
601 
602  double fitpar[3];
603  LinFit(x_hit, y_hit, y_error, fitpar);
604  double a = fitpar[0];
605  double b = fitpar[1];
606 
607 
608  int totAllowOutlier = static_cast<int>(0.1*n);
609  const double sigma_cut = 5;
610  int iOutlier = 0;
611  std::vector<double> x_filt, y_filt, ye_filt, x_filt_e;
612  for (int i = 0; i < n; i++) {
613  double y_fit = a + b*x_hit[i];
614  if ( fabs(y_hit[i] - y_fit) < sigma_cut*y_error[i] || iOutlier>totAllowOutlier)
615  { // > 5 * sigma. Accept 10% of 5sigma outliers. These numbers should be optimized.
616  x_filt.push_back(x_hit[i]);
617  y_filt.push_back(y_hit[i]);
618  ye_filt.push_back(y_error[i]);
619  x_filt_e.push_back(0);
620  }
621  else
622  iOutlier++;
623 
624  }
625 
626  LinFit(x_filt, y_filt, ye_filt, fitpar);
627  a = fitpar[0];
628  b = fitpar[1];
629  n = x_filt.size();
630  if (n < 5){
631  slope = 0;
632  chi2 = 999;
633  P_dn = 1e-30;
634  P_dn = 1e-30;
635  return;
636  }
637 
638 
639  slope = fitpar[1];
640 
641  // chi2
642  chi2=0.0;
643  for (int i=0; i<n; i++) {
644  double y_expected = a + b*x_filt.at(i);
645  chi2+=pow((y_filt.at(i)-y_expected) / ye_filt.at(i), 2.0);
646  }
647 
648  chi2/=static_cast<double>(n-2); // NDF = N points - 2 parameters in the fit
649 
650 
651  // Calculate up/down intercepts
652  double one_over_ee = 0.0;
653  double x_over_ee = 0.0;
654  double y_over_ee = 0.0;
655 
656  for (int i=0; i<n; i++) {
657  double e = ye_filt.at(i);
658  one_over_ee+=1.0/e/e;
659  x_over_ee+=x_filt.at(i)/e/e;
660  y_over_ee+=y_filt.at(i)/e/e;
661  }
662 
663  // up/down intercepts defined below
664  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
665  double down_inter = (y_over_ee+x_over_ee)/one_over_ee;
666 
667 
668  // Calculate up/down chi2
669  double up_chi2=0.0, down_chi2=0.0;
670  for (int i=0;i<n;i++) {
671 
672  double e = ye_filt.at(i);
673  double up_expected = up_inter + x_filt.at(i);
674  double down_expected = down_inter - x_filt.at(i);
675  up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
676  down_chi2 += pow((y_filt.at(i)-down_expected) / e, 2.0);
677  }
678 
679  // if statement is here for safety to prevent gamma_q failing
680 
681  double prob_up = boost::math::gamma_q(static_cast<double>(n-2)/2.0,up_chi2/2.0);
682  double prob_down = boost::math::gamma_q(static_cast<double>(n-2)/2.0,down_chi2/2.0);
683 
684  // Use Root TMath if you prefer. Root and boost libraries give identical results.
685  // double prob_up = TMath::Prob(up_chi2, n-2);
686  // double prob_down = TMath::Prob(down_chi2, n-2);
687 
688  if (prob_up<1e-30) prob_up = 1e-30;
689  if (prob_down<1e-30) prob_down = 1e-30;
690 
691  P_up = prob_up;
692  P_dn = prob_down;
693 
694  if(false){
695 
696  TCanvas *can = new TCanvas("can", "can", 800, 600);
697  // Debugging part. Will produce the time distribution fits using Root fitter and analytical fuction
698  TGraphErrors *eTmT_grx = new TGraphErrors(x_filt.size(), &x_filt[0], &y_filt[0], &x_filt_e[0], &ye_filt[0]);
699  TF1 *SlopeFit = new TF1("SlopeFit", "pol1", -10, 150);
700  TF1 *UpFit = new TF1("UpFit", "pol1", -10, 150);
701  TF1 *UpFit1 = new TF1("UpFit1", "pol1", -10, 150);
702  TF1 *DnFit2 = new TF1("DnFit2", "pol1", -10, 150);
703  TF1 *DnFit = new TF1("DnFit", "pol1", -10, 150);
704  UpFit->SetLineColor(4);
705  DnFit->SetLineColor(4);
706  DnFit->SetLineStyle(3);
707  UpFit->SetParameter(0, up_inter);
708  UpFit->SetParameter(1, 1);
709 
710  DnFit2->SetParameter(0, fitpar[0]);
711  DnFit2->SetParameter(1, fitpar[1]);
712  DnFit2->SetLineColor(kGreen);
713  DnFit2->SetLineStyle(3);
714 
715 
716  DnFit->SetParameter(0, down_inter);
717  DnFit->SetParameter(1, -1);
718 
719 
720 
721  SlopeFit->SetParameter(0, a);
722  SlopeFit->SetParameter(1, b);
723  SlopeFit->SetLineColor(3);
724  SlopeFit->SetLineStyle(3);
725  can->Draw();
726  can->cd();
727 
728  eTmT_grx->GetXaxis()->SetRangeUser(-10, 100);
729  eTmT_grx->Draw("AP");
730  SlopeFit->Draw("same");
731  // UpFit->Draw("same");
732  DnFit2->Draw("same");
733  DnFit->Draw("same");
734 
735  UpFit1->SetLineColor(kRed);
736  UpFit1->SetLineWidth(1);
737  // UpFit1->SetLineStyle(4);
738 
739  eTmT_grx->Fit("UpFit1", "BR", "", -10, 150);
740 
741  TPaveText *pt = new TPaveText(.75, .9, 0.9, 1.0, "brNDC");
742  pt->SetBorderSize(1);
743  pt->AddText(Form("P_up = %.1e", P_up ));
744  pt->AddText(Form("P_dn = %.1e", P_dn ));
745  pt->AddText(Form("LLR = %.1f", log(P_up/P_dn) ));
746  pt->AddText(Form("Chi2 = %.2f", chi2 ));
747  pt->AddText(Form("NDF = %d", n-2 ));
748 
749 
750  pt->Draw();
751  can->Update();
752  can->SaveAs(Form("plots/Chi2_%f.pdf", chi2));
753  // can->Delete();
754  }
755 
756 
757 
758 
759 
760 }
761 
762 double novaddt::UpMuTrigger::getDCMOff(uint16_t cell, uint8_t view){
763 
764  // Cell number starts from 1
765 
766  if (_detector == "ndos") {
767  if (view == daqchannelmap::Y_VIEW)
768  if (cell > 31 && cell < 64) return -26;
769  else if (cell > 63) return -40;
770  else return 0;
771  else // X_VIEW
772  if (cell > 7 && cell < 32) return -10;
773  else if (cell > 31) return 5;
774  else return 0;
775  }
776  else if (_detector == "fd")
777  return (cell/64)*40.;
778  else return 0;
779 }
780 
double GetT(double tdc, double dist_to_readout)
Double_t xx
Definition: macro.C:12
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
value_type val
Definition: BaseProducts.h:34
Float_t y1[n_points_granero]
Definition: compare.C:5
double getDCMOff(uint16_t cell, uint8_t view)
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
value_type val
Definition: BaseProducts.h:109
constexpr T pow(T x)
Definition: pow.h:75
Definition: event.h:19
const double T0
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: EventID.h:98
double chi2()
TVector3 const & End() const
Definition: Track3D.h:46
void LinFitLLR(std::vector< double > &x_hit, std::vector< double > &y_hit, std::vector< double > &y_error, double &slope, double &chi2, double &P_up, double &P_dn)
Identifier for the Y measuring view of the detector (side)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
length
Definition: demo0.py:21
const double a
novaddt::View const & View() const
Definition: DAQHit.h:72
Identifier for the X measuring view of the detector (top)
Float_t d
Definition: plot.C:236
TVector3 const & Start() const
Definition: Track3D.h:45
value_type val
Definition: BaseProducts.h:84
bool const & Is3D() const
Definition: Track3D.h:43
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
EventNumber_t event() const
Definition: EventID.h:116
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
UpMuTrigger(fhicl::ParameterSet const &p)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
virtual bool filter(art::Event &e) override
const hit & b
Definition: hits.cxx:21
value_type val
Definition: BaseProducts.h:137
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
Float_t e
Definition: plot.C:35
double GetPigtail(DAQHit const &daqhit) const
Definition: fwd.h:28
void GetXYZ(DAQHit const &daqhit, TVector3 start, TVector3 end, double *xyzd)
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_q.hpp:12
EventID id() const
Definition: Event.h:56