UpMuAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: UpMuAna
3 // Module Type: analyzer
4 // File: UpMuAna_module.cc
5 //
6 ////////////////////////////////////////////////////////////////////////
7 
12 
19 #include "fhiclcpp/ParameterSet.h"
20 
21 #include "RawData/RawTrigger.h"
22 #include "Simulation/FLSHit.h"
23 #include "Simulation/FLSHitList.h"
29 #include "TrackFit/LocatedChan.h"
30 #include "RecoBase/CellHit.h"
31 #include <TF1.h>
32 #include <TH1F.h>
33 #include <TCanvas.h>
34 #include <TApplication.h>
35 #include <TGraphErrors.h>
36 #include <TNtuple.h>
37 #include <TPaveText.h>
38 #include <TMath.h>
39 #include <algorithm>
40 #include <vector>
41 #include <numeric>
42 #include <boost/math/special_functions/gamma.hpp>
43 #include <boost/math/special_functions.hpp>
44 #include "Simulation/FLSHit.h"
45 #include "Simulation/FLSHitList.h"
46 #include "Simulation/Particle.h"
49 #include "Utilities/NTree.h"
50 #include "MCCheater/BackTracker.h"
52 
53 
54 
55 namespace novaddt {
56  class UpMuAna;
57 }
58 
60 public:
61  explicit UpMuAna(fhicl::ParameterSet const & p);
62  virtual ~UpMuAna();
63 
64  void analyze(art::Event const & e) override;
65 
66  void beginJob() override;
67 
68 private:
69  bool isMC;
73  unsigned int _nSlices = 0;
74  unsigned int _nTracks = 0;
75  const float c0 = -759.5; // x or y coordinates of cell=0 [cm]
76  const float p0 = 4.57; // z coordinates of plane=0 [cm]
77  const float cw = 3.97; // Cell width [cm]
78  const float pw = 6.654; // Plane width [cm]
79 
80  std::string _trackModuleLabel; ///< label of module making the Tracks
81  std::string _trackInstanceLabel; ///< instance label making the Tracks
83  std::string _sliceModuleLabel; ///< label of module making the HitList
84  std::string _sliceInstanceLabel; ///< instance label making the HitList
87 
88  float _TrackLen = 0;
89  unsigned _TrackHitsXY = 0;
90  unsigned _TrackHitsX = 0;
91  unsigned _TrackHitsY = 0;
92  int _dX = 2;
93  int _dY = 5;
94  int _dZ = 2;
95  float _Chi2 = 999.0;
96  float _LLR = -999.0;
97  float _R2X = 0.0;
98  float _R2Y = 0.0;
99  float _MinSlope = -999.0;
100  float _MaxSlope = 999.0;
101  float TrueLength = -999.;
102  float TrueEDep = -999.;
103 
104 
105 
106  void GetXYZD(DAQHit const& daqhit, double w, double *xyzd);
107 
108  float GetZ(int plane);
109  float GetT(double TDC, double dist_to_readout);
110  float GetRes(int ADC);
111 
112  float GetX(int cell);
113  void GetXYZ(DAQHit const& daqhit, TVector3 start, TVector3 end, double *xyzd);
114  void LinFit(const std::vector<double>& x, const std::vector<double>& y, float *fitpar);
115  void LinFit(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& ye, float *fitpar);
116 
117  void LinFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error,
118  double& slope, double& chi2, double& P_up, double& P_dn);
119 
120 
121  float RootFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error, float slope);
122 
123  double GetPigtail(DAQHit const& daqhit) const;
124 
125  void populateFLSmap(std::map<geo::OfflineChan, const sim::FLSHit*> &pcpflsmap,
126  std::map<geo::OfflineChan, unsigned short> &pcpnummap,
127  art::Event const &e);
128 
129  std::map<std::string, TNtuple*> ntup_;
130  TH1F *dTdist;
131  TCanvas *can;
132 
134 };
135 
136 
138  : EDAnalyzer(p)
139  , isMC (p.get<bool> ("isMC"))
140  , UseTimeOffsets (p.get<bool> ("UseTimeOffsets"))
141  , _trackModuleLabel (p.get<std::string >("TrackModuleLabel" ))
142  , _trackInstanceLabel(p.get<std::string >("TrackInstanceLabel", ""))
143  , _trackToSliceInstanceLabel(p.get<std::string >("TrackToSliceInstanceLabel", ""))
144  , _sliceModuleLabel (p.get<std::string >("SliceModuleLabel" ))
145  , _sliceInstanceLabel(p.get<std::string >("SliceInstanceLabel", ""))
146  , _TrackLen (p.get<float>("TrackLen"))
147  , _TrackHitsXY (p.get<unsigned>("TrackHitsXY"))
148  , _TrackHitsX (p.get<unsigned>("TrackHitsX"))
149  , _TrackHitsY (p.get<unsigned>("TrackHitsY"))
150  , _dX (p.get<int>("dX"))
151  , _dY (p.get<int>("dY"))
152  , _dZ (p.get<int>("dZ"))
153  , _Chi2 (p.get<float>("Chi2"))
154  , _LLR (p.get<float>("LLR"))
155  , _R2X (p.get<float>("R2X"))
156  , _R2Y (p.get<float>("R2Y"))
157  , _MinSlope (p.get<float>("MinSlope"))
158  , _MaxSlope (p.get<float>("MaxSlope"))
159 {
160  if (isMC) {
161  _FLSHitListLabel = p.get<std::string>("FLSHitListLabel");
162  _RawDigitLabel = p.get<std::string>("RawDigitLabel");
163  }
164 }
165 
167 {
168 }
169 
171 {
172 
173  ntup_["particle_ntuple"] = tfs_->make<TNtuple>
174  ("particle_ntuple", "Upward going particles",
175  "Run:Event:ParticleID:PDGid:InitKE:endKE:MuEneDep:MuInitKE:Mother:Process");
176 
177  ntup_["hit_ntuple"] = tfs_->make<TNtuple>
178  ("hit_ntuple", "Per hit variables",
179  "Run:SubRun:Event:TrackID:iHit:View:eT:mT:mT_uncor:compMT:truthT:PDG:TDC:ADC:Plane:Cell:DCM:X:Y:Z:D:TrackLen");
180 
181  ntup_["track_ntuple"] = tfs_->make<TNtuple>
182  ("track_ntuple", "Slope NTuple",
183  "Run:SubRun:Event:StartX:StartY:StartZ:EndX:EndY:EndZ:dTX:dTY:TrackLen:TrueLenght:TrueEDep:TrackID:R2X:R2Y:nXhits:nYhits:SlopeX:SlopeY:SlopeXY:Chi2X:Chi2Y:Chi2XY:ProbUpX:ProbDnX:ProbUpY:ProbDnY:ProbUpXY:ProbDnXY:ProbUpRoot:ProbDnRoot");
184 
185  can = tfs_->make<TCanvas> ("can", "can", 800, 600);
186 
187 }
188 
190 {
191 
192  std::map<geo::OfflineChan, const sim::FLSHit*> pcpflsmap;
193  std::map<geo::OfflineChan, unsigned short> pcpnummap;
194  if (isMC) populateFLSmap(pcpflsmap, pcpnummap, e);
195 
196 
199 
201  e.getByLabel("track", "", tracks);
202 
203  std::map<unsigned short, std::set<unsigned short>> planesAndCells;
204 
205  // get the hit list for each track
207  // get the tracks for each hit list
209 
210  art::FindManyP<novaddt::Track3D> find_tracks(slices, e, the_slices);
211  assert(find_tracks.isValid());
212 
213  // Loop through the slices
214  for(size_t i_slice = 0; i_slice < slices->size(); ++i_slice){
215 
216  novaddt::HitList hit_list = slices->at(i_slice);
217  novaddt::TDC absolute_t0 = hit_list.front().TDC().val;
218 
219 
220  // find the tracks for this hit list
221  std::vector< art::Ptr<novaddt::Track3D> > tracks = find_tracks.at(i_slice);
222  if(tracks.size() == 0) continue;
223  // get the hit list for each track
225  assert(fohl.isValid());
226  assert(fohl.size() == tracks.size());
227 
228 
229 
230  size_t keep_track = 0;
231  size_t keep_hlist = 0;
232  // int StartZ = track->Start().Z();
233  // int EndZ = track->End().Z();
234 
235  for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
236  art::Ptr<novaddt::HitList> this_hit_list = fohl.at(i_track);
237  art::Ptr<novaddt::Track3D> track = tracks[i_track];
238  if (!track->Is3D()) continue;
239 
240  if(this_hit_list->size() > keep_hlist){
241  keep_track = i_track;
242  keep_hlist = this_hit_list->size();
243  }
244  }
245 
246 
247  for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
248 
249  if(keep_track != i_track) continue;
250 
251  art::Ptr<novaddt::Track3D> track = tracks[i_track];
252  // get the hit list for this track
253  art::Ptr<novaddt::HitList> this_hit_list = fohl.at(i_track);
254 
255  int StartX = track->Start().X();
256  int StartY = track->Start().Y();
257  int StartZ = track->Start().Z();
258 
259  int EndX = track->End().X();
260  int EndY = track->End().Y();
261  int EndZ = track->End().Z();
262  if(StartY > EndY){
263  std::swap(StartX, EndX);
264  std::swap(StartY, EndY);
265  std::swap(StartZ, EndZ);
266  }
267 
268  TVector3 length;
269  length.SetXYZ(cw*(track->End().X() - track->Start().X()),
270  cw*(track->End().Y() - track->Start().Y()),
271  pw*(track->End().Z() - track->Start().Z())
272  );
273 
274  float tr_length = length.Mag();
275  float exp_trav_time = tr_length/29.97; // c=29.97 cm/ns; 1/64E6=15.625 ns
276 
277  if (tr_length < _TrackLen) continue; // Let's select tracks 20 m long.
278  if (!track->Is3D()) continue;
279  if (TMath::Abs(EndX-StartX) < _dX) continue; // Select tracks longer than _dX cells in X.
280  if (TMath::Abs(EndY-StartY) < _dY) continue; // Select tracks longer than _dY cells in Y.
281  if (TMath::Abs(EndZ-StartZ) < _dZ) continue; // Select tracks longer than _dZ planes in Z.
282 
283 
284  std::cout.precision(4);
285  if(e.id().event() == 0){
286  std::cout << "\t3D-track[" << i_track
287  << "]\t view: " << track->View()
288  << "\thits: " << this_hit_list->size()
289  << "\tX: " << StartX << " - " << EndX
290  << "\tY: " << StartY << " - " << EndY
291  << "\tZ: " << StartZ << " - " << EndZ
292  << "\tLen: " << tr_length/100
293  << "\tExp time: " << exp_trav_time << "\t3D: " << track->Is3D()
294  << std::endl;
295  }
296 
297 
298  int T0 = 0;
299  int T0X = 0;
300  int T0Y = 0;
301  int T1X = 0;
302  int T1Y = 0;
303 
304  std::vector<double> x_hit;
305  std::vector<double> y_hit;
306  std::vector<double> zy_hit;
307  std::vector<double> zx_hit;
308 
309  // Get T0 that corresponds to the minimum cell number in Y-view
310  for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
311  uint16_t iC = this_hit_list->at(i_hit).Cell().val;
312  uint16_t iP = this_hit_list->at(i_hit).Plane().val;
313  uint8_t view = this_hit_list->at(i_hit).View().val;
314  int iTDC = int(this_hit_list->at(i_hit).TDC().val - absolute_t0);
315 
316  if(view == daqchannelmap::X_VIEW){
317  if(iC == StartX)
318  T0X = iTDC;
319  if(iC == EndX)
320  T1X = iTDC;
321  x_hit.push_back(iC);
322  zx_hit.push_back(iP);
323  }
324 
325  if(view == daqchannelmap::Y_VIEW){
326  if(iC == StartY)
327  T0Y = iTDC;
328  if(iC == EndY)
329  T1Y = iTDC;
330  y_hit.push_back(iC);
331  zy_hit.push_back(iP);
332  }
333 
334  T0 = T0Y;
335  }
336 
337 
338  // Linear fit of the hits in XZ/YZ views. Drop tracks that has bad linear fit
339  float fitpar[3];
340  LinFit(x_hit, zx_hit, fitpar);
341  float r2x = fitpar[2];
342  LinFit(y_hit, zy_hit, fitpar);
343  float r2y = fitpar[2];
344 
345  if(r2x < _R2X) continue;
346  if(r2y < _R2Y) continue;
347 
348  int dTX = (T0X - T1X);
349  int dTY = (T0Y - T1Y);
350  // Vectors of expected and measured times for TGraph
351  std::vector<Double_t> eT_vecx;
352  std::vector<Double_t> mT_vecx;
353  std::vector<Double_t> eT_vecy;
354  std::vector<Double_t> mT_vecy;
355  std::vector<Double_t> eT_vecxy;
356  std::vector<Double_t> mT_vecxy;
357 
358  // errors
359  std::vector<Double_t> mT_vecx_e;
360  std::vector<Double_t> mT_vecy_e;
361  std::vector<Double_t> mT_vecxy_e;
362  std::vector<Double_t> eT_vecx_e;
363  std::vector<Double_t> eT_vecy_e;
364  std::vector<Double_t> eT_vecxy_e;
365 
366  for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
367  uint16_t iC = this_hit_list->at(i_hit).Cell().val;
368  uint16_t iP = this_hit_list->at(i_hit).Plane().val;
369  int iTDC = int(this_hit_list->at(i_hit).TDC().val - absolute_t0);
370  double meas_time = (double)(iTDC - T0);
371  meas_time = meas_time + 0.01*static_cast<double>((this_hit_list->at(i_hit).TDC().frac));
372  int hit_adc = (int)this_hit_list->at(i_hit).ADC().val;
373 
374  if (hit_adc < 50) continue;
375 
376  // Calculate time of the hit
377  DAQHit const& dqhit = this_hit_list->at(i_hit);
378  double pos[4];
379  TVector3 start = track->Start();
380  TVector3 end = track->End();
381  GetXYZ(dqhit, start, end, pos);
382  double pigtail = GetPigtail(dqhit);
383  double tns = GetT(meas_time, pos[3] + pigtail);
384 
385  float compMT = 0;
386  double truthT = 0;
387  int PDG = -999;
388 
389  if (isMC) {
390  geo::OfflineChan chan(iP, iC);
391  const sim::FLSHit *fls = pcpflsmap[chan];
392  const short flag = pcpnummap[chan];
393  if (fls && flag==1) {// only 1 entry for this channel
394  compMT = GetT(meas_time + (int)(T0 - absolute_t0), pos[3]);
395  truthT = 0.5*(fls->GetEntryT() + fls->GetExitT());
396  PDG = fls->GetPDG();
397 
398  }
399  }
400 
401 
402 
403 
404 
405  double tns2 = tns;
406  // Subtract DCM offsects
407  double DCM_off = 0.;
408  if(UseTimeOffsets)
409  DCM_off = (iC/64)*40.;
410 
411  tns = tns - DCM_off;
412 
413 
414  if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW){
415  if(iP > TMath::Max(StartZ, EndZ)) continue;
416  if(iP < TMath::Min(StartZ, EndZ)) continue;
417 
418  if(iC > TMath::Max(StartX, EndX)) continue;
419  if(iC < TMath::Min(StartX, EndX)) continue;
420  }
421  else{
422  if(iP > TMath::Max(StartZ, EndZ)) continue;
423  if(iP < TMath::Min(StartZ, EndZ)) continue;
424 
425  if(iC > TMath::Max(StartY, EndY)) continue;
426  if(iC < TMath::Min(StartY, EndY)) continue;
427  }
428 
429  float exp_time = -999;
430  std::string h_view = "y";
431 
432  if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW){
433  h_view = "x";
434  exp_time = exp_trav_time*(iC-StartX)/(EndX-StartX);
435  }
436  else
437  exp_time = exp_trav_time*(iC-StartY)/(EndY-StartY);
438 
439  const bool debug = false;
440  if(debug){
441  std::cout<< "hit[" << i_hit
442  << "]: \tp: " << this_hit_list->at(i_hit).Plane().val
443  << ",\tc: " << this_hit_list->at(i_hit).Cell().val
444  << ",\tDCM: " << this_hit_list->at(i_hit).getDCM().val
445  << ",\t view: " << h_view
446  << ",\tTDC: " << this_hit_list->at(i_hit).TDC().val
447  << ",\tT0: " << T0
448  << ",\texp_t: " << exp_time
449  << ",\tmt: " << tns
450  << ",\tmt2: " << tns2
451  << ",\tXP: " << pos[0]
452  << ",\tYP: " << pos[1]
453  << std::endl;
454  }
455 
456  assert(pos[3]>0);
457 
458  if(this_hit_list->at(i_hit).View().val == daqchannelmap::Y_VIEW){
459  eT_vecy.push_back(exp_time);
460  mT_vecy.push_back(tns);
461  eT_vecy_e.push_back(0);
462  mT_vecy_e.push_back(GetRes(hit_adc));
463  }
464 
465  if(this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW){
466  eT_vecx.push_back(exp_time);
467  mT_vecx.push_back(tns);
468  eT_vecx_e.push_back(0);
469  mT_vecx_e.push_back(GetRes(hit_adc));
470  }
471 
472  eT_vecxy.push_back(exp_time);
473  mT_vecxy.push_back(tns);
474  eT_vecxy_e.push_back(0);
475  mT_vecxy_e.push_back(GetRes(hit_adc));
476 
477 
478  float ntuple_entries[22] = {
479  (float)e.id().run(),(float)e.id().subRun(), (float)e.id().event(), (float)_nTracks,
480  (float)i_hit, (float)(this_hit_list->at(i_hit).View().val == daqchannelmap::Y_VIEW), exp_time,
481  (float)tns, (float)tns2, compMT,
482  (float) truthT, (float)PDG, (float)meas_time, (float)this_hit_list->at(i_hit).ADC().val,
483  (float)this_hit_list->at(i_hit).Plane().val, (float)this_hit_list->at(i_hit).Cell().val, (float)this_hit_list->at(i_hit).getDCM().val,
484  (float)pos[0], (float)pos[1], (float)pos[2],
485  (float)pos[3], (float)tr_length
486  };
487  ntup_["hit_ntuple"]->Fill(ntuple_entries);
488 
489  } // end of loop on hits
490 
491 
492  _nTracks++;
493 
494  unsigned int nxhit = eT_vecx.size();
495  unsigned int nyhit = eT_vecy.size();
496 
497  if(nxhit < _TrackHitsX) continue;
498  if(nyhit < _TrackHitsY) continue;
499  if((nxhit+nyhit) < _TrackHitsXY) continue;
500 
501  // Calculate the likelihood of upward/downward going hypothesis
502  float P_up_root = 0;
503  float P_down_root = 0;
504  P_up_root = RootFitLLR(eT_vecxy, mT_vecxy, mT_vecxy_e, 1);
505  P_down_root = RootFitLLR(eT_vecxy, mT_vecxy, mT_vecxy_e, -1);
506 
507  if(P_down_root<1e-30) P_down_root = 1e-30;
508  if(P_up_root<1e-30) P_up_root = 1e-30;
509 
510  double slope_x, chi_x, P_up_x, P_dn_x;
511  double slope_y, chi_y, P_up_y, P_dn_y;
512  double slope_xy, chi_xy, P_up_xy, P_dn_xy;
513 
514  LinFitLLR(eT_vecx, mT_vecx, mT_vecx_e, slope_x, chi_x, P_up_x, P_dn_x);
515  LinFitLLR(eT_vecy, mT_vecy, mT_vecy_e, slope_y, chi_y, P_up_y, P_dn_y);
516  LinFitLLR(eT_vecxy, mT_vecxy, mT_vecxy_e, slope_xy, chi_xy, P_up_xy, P_dn_xy);
517 
518  if(slope_xy < _MinSlope) continue;
519  if(slope_xy > _MaxSlope) continue;
520 
521  // Remove tracks with LLR smaller than input in the config fcl file
522  double LLR = log(P_up_xy/P_dn_xy);
523  if(LLR < _LLR) continue;
524  // Remove tracks with Chi2 fit values higher than input in the config fcl file
525  if(chi_xy > _Chi2) continue;
526 
527 
528  float tr_ntuple[33] = {(float)e.id().run(),(float)e.id().subRun(), (float)e.id().event(),
529  (float)StartX, (float)StartY, (float)StartZ, (float)EndX, (float)EndY, (float)EndZ, (float)dTX, (float)dTY,
530  (float)tr_length, TrueLength, TrueEDep, (float)_nTracks,
531  (float)r2x, (float)r2y,
532  (float)nxhit, (float)nyhit,
533  (float)slope_x, (float)slope_y, (float)slope_xy,
534  (float)chi_x, (float)chi_y, (float)chi_xy,
535  (float)P_up_x, (float)P_dn_x,
536  (float)P_up_y, (float)P_dn_y,
537  (float)P_up_xy, (float)P_dn_xy,
538  (float)P_up_root, (float)P_down_root};
539 
540  ntup_["track_ntuple"]->Fill(tr_ntuple);
541 
542  } // end of loop on tracks
543 
544  } // end of loop on slices
545 
546 }
547 
548 void novaddt::UpMuAna::GetXYZ(DAQHit const& daqhit, TVector3 start, TVector3 end, double *xyzd) {
549  UInt_t plane = daqhit.Plane().val;
550  UInt_t cell = daqhit.Cell().val;
551 
552  // The parameters are taken from the fit of plane vs Z
553  float x = c0 + cw*cell;
554  float y = x;
555  float z = p0 + pw*plane;
556  float d = -999.9;
557 
558  float x0 = c0 + cw*start.X(); float x1 = c0 + cw*end.X();
559  float y0 = c0 + cw*start.Y(); float y1 = c0 + cw*end.Y();
560  float z0 = p0 + pw*start.Z(); float z1 = p0 + pw*end.Z();
561 
562  if(daqhit.View().val == daqchannelmap::X_VIEW){
563  y = (y1 - y0)/(z1 - z0);
564  y = y*(z - z0);
565  y = y + y0;
566  d = 800 - y;
567  }
568  if(daqhit.View().val == daqchannelmap::Y_VIEW){
569  x = (x1 - x0)/(z1 - z0);
570  x = x*(z - z0);
571  x = x + x0;
572  d = 800 - x;
573  }
574 
575  if(d<0) d=0;
576  xyzd[0] = x;
577  xyzd[1] = y;
578  xyzd[2] = z;
579  xyzd[3] = d;
580 }
581 
582 float novaddt::UpMuAna::GetT(double TDC, double dist_to_readout) {
583  // This needs to be calibrated and put in a DB table
584  float speedOfFiberTransport = 15.3; // cm/ns, "first principles" calib.
585  float TDC_to_ns = 15.625; // 1/64E6=15.625 ns
586 
587  // Differs from c/n due to zigzag
588  // paths down fiber. But, this
589  // value may be way off (10%? 30%?).
590  return (TDC*TDC_to_ns - dist_to_readout/speedOfFiberTransport);
591 }
592 
593 
595 
596  //The fit function from Evan: f(x) = p0/(p1+x^p2) + p3;
597  float p0 = 70267.4;
598  float p1 = 677.305;
599  float p2 = 1.48702;
600  float p3 = 8.46997;
601 
602  float res = p0/(p1 + pow( (float)ADC, p2) ) + p3;
603  return res;
604 
605 }
606 
607 
608 
609 double novaddt::UpMuAna::GetPigtail(DAQHit const& daqhit) const{
610  UInt_t cellid = daqhit.Cell().val;
611  //NOTE: Do we think 32 cells per module can ever change
612  //This number can be found in DAQChannelMap but did not want to add that
613  //dependancy to geometry unless needed.
614  const int kCellsPerModule = 32;
615  cellid = cellid % kCellsPerModule;
616  // In vertical planes, high cell numbers have longer fibres.
617  // In horizontal planes it's the opposite, so flip it round.
618  if(daqhit.View().val == daqchannelmap::X_VIEW)
619  cellid = kCellsPerModule-cellid-1;
620 
621  // This should really never happen, but just to be safe...
622  if(cellid >= kCellsPerModule) return 100;
623  // Email from Tom Chase 2011-04-29
624  // NB: this isn't just a ~3.8cm change per cell. Every 8 cells something
625  // different happens.
626  const double kPigtails[kCellsPerModule] = {
627  34.5738, 38.4379, 42.3020, 46.1660, 50.0301, 53.8942, 57.7582, 61.6223,
628  64.7504, 68.6144, 72.4785, 76.3426, 80.2067, 84.0707, 87.9348, 91.0790,
629  95.3301, 99.1941, 103.058, 106.922, 110.786, 114.650, 118.514, 122.379,
630  125.507, 129.371, 133.235, 137.099, 140.963, 144.827, 148.691, 150.751
631  };
632  return kPigtails[cellid];
633  }
634 
635 
636 void novaddt::UpMuAna::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, float *fitpar){
637  // http://en.wikipedia.org/wiki/Simple_linear_regression
638  const auto n = x_val.size();
639  const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n; // <x>
640  const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n; // <y>
641  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n; // <xx>
642  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n; // <yy>
643  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n; // <xy>
644 
645  const auto b = (xy - x*y)/(xx - x*x); // slope
646  const auto a = y - b*x; // intercept
647  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
648  fitpar[0] = a;
649  fitpar[1] = b;
650  fitpar[2] = r*r;
651 
652 }
653 
654 
655 void novaddt::UpMuAna::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, const std::vector<double>& y_err, float *fitpar){
656  // http://en.wikipedia.org/wiki/Simple_linear_regression
657  int n = x_val.size();
658  double x = 0;
659  double y = 0;
660  double xx = 0;
661  double yy = 0;
662  double xy = 0;
663  double ee = 0;
664 
665  for ( int i=0; i<n; i++ ){
666 
667  x = x + x_val[i]/y_err[i]/y_err[i];
668  y = y + y_val[i]/y_err[i]/y_err[i];
669 
670  xx = xx + x_val[i]*x_val[i]/y_err[i]/y_err[i];
671  yy = yy + y_val[i]*y_val[i]/y_err[i]/y_err[i];
672  xy = xy + x_val[i]*y_val[i]/y_err[i]/y_err[i];
673  ee = ee + 1./y_err[i]/y_err[i];
674  }
675 
676  const auto b = (xy*ee - x*y)/(xx*ee - x*x); // slope
677  const auto a = (xy - b*xx)/x; // intercept
678  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
679  fitpar[0] = a;
680  fitpar[1] = b;
681  fitpar[2] = r*r;
682 
683 
684 }
685 
686 
687 
688 void novaddt::UpMuAna::LinFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error,
689  double& slope, double& chi2, double& P_up, double& P_dn) {
690 
691  int n = x_hit.size();
692 
693  float fitpar[3];
694  LinFit(x_hit, y_hit, y_error, fitpar);
695  float a = fitpar[0];
696  float b = fitpar[1];
697 
698 
699  int totAllowOutlier = (int)(0.1*n);
700  int iOutlier = 0;
701  std::vector<double> x_filt, y_filt, ye_filt, x_filt_e;
702  for (int i = 0; i < n; i++) {
703  float y_fit = a + b*x_hit[i];
704  if ( fabs(y_hit[i] - y_fit) < 5*y_error[i] || iOutlier>totAllowOutlier)
705  { // > 5 * sigma. Accept 10% of 5sigma outliers. These numbers should be optimized.
706  x_filt.push_back(x_hit[i]);
707  y_filt.push_back(y_hit[i]);
708  ye_filt.push_back(y_error[i]);
709  x_filt_e.push_back(0);
710  }
711  else
712  iOutlier++;
713 
714  }
715 
716  LinFit(x_filt, y_filt, ye_filt, fitpar);
717  a = fitpar[0];
718  b = fitpar[1];
719  n = x_filt.size();
720  if (n < 5){
721  slope = 0;
722  chi2 = 999;
723  P_dn = 1e-30;
724  P_dn = 1e-30;
725  return;
726  }
727 
728 
729  slope = fitpar[1];
730 
731  // chi2
732  chi2=0.0;
733  for (int i=0; i<n; i++) {
734  double y_expected = a + b*x_filt.at(i);
735  chi2+=pow((y_filt.at(i)-y_expected) / ye_filt.at(i), 2.0);
736  }
737 
738  chi2/=(double)(n-2); // NDF = N points - 2 parameters in the fit
739 
740 
741  // Calculate up/down intercepts
742  double one_over_ee = 0.0;
743  double x_over_ee = 0.0;
744  double y_over_ee = 0.0;
745 
746  for (int i=0; i<n; i++) {
747  double e = ye_filt.at(i);
748  one_over_ee+=1.0/e/e;
749  x_over_ee+=x_filt.at(i)/e/e;
750  y_over_ee+=y_filt.at(i)/e/e;
751  }
752 
753  // up/down intercepts defined below
754  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
755  double down_inter = (y_over_ee+x_over_ee)/one_over_ee;
756 
757 
758  // Calculate up/down chi2
759  double up_chi2=0.0, down_chi2=0.0;
760  for (int i=0;i<n;i++) {
761 
762  double e = ye_filt.at(i);
763  double up_expected = up_inter + x_filt.at(i);
764  double down_expected = down_inter - x_filt.at(i);
765  up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
766  down_chi2 += pow((y_filt.at(i)-down_expected) / e, 2.0);
767  }
768 
769  // if statement is here for safety to prevent gamma_q failing
770 
771  double prob_up = boost::math::gamma_q((double)(n-2)/2.0,up_chi2/2.0);
772  double prob_down = boost::math::gamma_q((double)(n-2)/2.0,down_chi2/2.0);
773 
774  // Use Root TMath if you prefer. Root and boost libraries give identical results.
775  // double prob_up = TMath::Prob(up_chi2, n-2);
776  // double prob_down = TMath::Prob(down_chi2, n-2);
777 
778  if (prob_up<1e-30) prob_up = 1e-30;
779  if (prob_down<1e-30) prob_down = 1e-30;
780 
781  P_up = prob_up;
782  P_dn = prob_down;
783 
784 
785  if(false){
786  // Debugging part. Will produce the time distribution fits using Root fitter and analytical fuction
787  TGraphErrors *eTmT_grx = new TGraphErrors(x_filt.size(), &x_filt[0], &y_filt[0], &x_filt_e[0], &ye_filt[0]);
788  TF1 *SlopeFit = new TF1("SlopeFit", "pol1", -10, 150);
789  TF1 *UpFit = new TF1("UpFit", "pol1", -10, 150);
790  TF1 *UpFit1 = new TF1("UpFit1", "pol1", -10, 150);
791  TF1 *DnFit2 = new TF1("DnFit2", "pol1", -10, 150);
792  TF1 *DnFit = new TF1("DnFit", "pol1", -10, 150);
793  UpFit->SetLineColor(4);
794  DnFit->SetLineColor(4);
795  DnFit->SetLineStyle(3);
796  UpFit->SetParameter(0, up_inter);
797  UpFit->SetParameter(1, 1);
798 
799  DnFit2->SetParameter(0, fitpar[0]);
800  DnFit2->SetParameter(1, fitpar[1]);
801  DnFit2->SetLineColor(kGreen);
802  DnFit2->SetLineStyle(3);
803 
804 
805  DnFit->SetParameter(0, down_inter);
806  DnFit->SetParameter(1, -1);
807 
808 
809 
810  SlopeFit->SetParameter(0, a);
811  SlopeFit->SetParameter(1, b);
812  SlopeFit->SetLineColor(3);
813  SlopeFit->SetLineStyle(3);
814  can->Draw();
815  can->cd();
816 
817  eTmT_grx->GetXaxis()->SetRangeUser(-10, 100);
818  eTmT_grx->Draw("AP");
819  SlopeFit->Draw("same");
820  // UpFit->Draw("same");
821  DnFit2->Draw("same");
822  DnFit->Draw("same");
823 
824  UpFit1->SetLineColor(kRed);
825  UpFit1->SetLineWidth(1);
826  // UpFit1->SetLineStyle(4);
827 
828  eTmT_grx->Fit("UpFit1", "BR", "", -10, 150);
829 
830  TPaveText *pt = new TPaveText(.75, .9, 0.9, 1.0, "brNDC");
831  pt->SetBorderSize(1);
832  pt->AddText(Form("P_up = %.1e", P_up ));
833  pt->AddText(Form("P_dn = %.1e", P_dn ));
834  pt->AddText(Form("LLR = %.1f", log(P_up/P_dn) ));
835  pt->AddText(Form("Chi2 = %.2f", chi2 ));
836  pt->AddText(Form("NDF = %d", n-2 ));
837 
838  pt->Draw();
839  can->Update();
840  can->SaveAs(Form("plots/Chi2_%f.pdf", chi2));
841  }
842 
843 }
844 
845 
846 float novaddt::UpMuAna::RootFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error, float slope){
847 
848  std::vector<double> x_error(x_hit.size(), 0.0);
849  TGraphErrors *gr = new TGraphErrors(x_hit.size(), &x_hit[0], &y_hit[0], &x_error[0], &y_error[0]);
850 
851  TF1 *SlopeFit = new TF1("SlopeFit", "pol1", -10, 150);
852  SlopeFit->FixParameter(1, slope);
853  gr->Fit("SlopeFit", "BRQ");
854 
855  int tot_out = (int)(0.05*gr->GetN());
856  int iout = 0;
857 
858  for(int i=0; i<gr->GetN(); i++){
859  float x = gr->GetX()[i];
860  float y = gr->GetY()[i];
861  float ye = gr->GetEY()[i];
862  float p0 = SlopeFit->GetParameter(0);
863  float p1 = SlopeFit->GetParameter(1);
864  float eTime = p0 + p1*x;
865  float sigma = TMath::Abs((eTime - y)/ye);
866 
867  const float sigma_cut = 5;
868  if(sigma>sigma_cut && iout<tot_out){
869  gr->RemovePoint(i);
870  iout++;
871  i--;
872  }
873  }
874 
875  gr->Fit("SlopeFit", "BRQ");
876  float fit_prob = SlopeFit->GetProb();
877 
878  SlopeFit->Delete();
879  gr->Delete();
880 
881  return fit_prob;
882 }
883 
884 
885 // Make a mapping from a channel to a FLSHit for later comparison.
886 // This method only works well for events with a single track.
888 (std::map<geo::OfflineChan, const sim::FLSHit*> &pcpflsmap,
889  std::map<geo::OfflineChan, unsigned short> &pcpnummap,
890  art::Event const &e)
891 {
892  // Get the event time to remove offset in later comparison.
895  e.getByLabel(_FLSHitListLabel, flshitlist);
896 
897  // if (flshitlist->empty()) std::cout << "FLSHitList is empty!" << std::endl;
898 
899  float muEneDep = 0;
900  float muEn = 0;
901  float minY = 999.;
902  float totLength = 0.;
903 
904  for (unsigned int i=0; i<flshitlist->size(); ++i) {
905  art::Ptr<sim::FLSHitList> hlist(flshitlist, i);
906  const std::vector<sim::FLSHit>& hits = hlist->fHits;
907 
908  for (unsigned int j=0; j<hits.size(); ++j) {
909  if (hits[j].GetTrackID() < 0) continue; // this is incoming particle
910 
911  unsigned int planeid = hits[j].GetPlaneID();
912  unsigned int cellid = hits[j].GetCellID();
913 
914  if(TMath::Abs(hits[j].GetPDG()) == 13){
915  muEneDep = muEneDep + hits[j].GetTotalEnergyLoss();
916  totLength = totLength + hits[j].GetTotalPathLength();
917  if(muEn < hits[j].GetEntryEnergy()) muEn = hits[j].GetEntryEnergy();
918  if(minY>hits[j].GetYAverage()) minY = hits[j].GetYAverage();
919  }
920 
921  geo::OfflineChan chan(planeid, cellid);
922  pcpflsmap[chan]= &hits[j];
923  pcpnummap[chan]++;
924  }
925  }
926 
927  TrueLength = totLength;
928  TrueEDep = muEneDep;
929 
930 
932  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
933  // loop thru the simulated particles
934  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
935  const sim::Particle* p = (*i).second;
936  int pdg = p->PdgCode();
937  int trID = p->TrackId();
938  int mother = p->Mother();
939  int process = -999;
940  float initKE = p->Momentum().E() - p->Momentum().M();
941  float muInitKE = muEn;
942  float endKE = p->EndMomentum().E() - p->Momentum().M();
943 
944  if(p->Process() == "Primary") process = 0; // primary muon or neutrino
945  if(p->Process() == "muMinusCaptureAtRest") process = 1; // mu- decay
946  if(p->Process() == "Decay") process = 2; // mu+ decay
947 
948  float ntuple_entries[10] = {
949  (float)e.id().run(), (float)e.id().event(), (float)trID,
950  (float)pdg, initKE, endKE, muEneDep, muInitKE, (float)mother, (float)process
951  };
952  ntup_["particle_ntuple"]->Fill(ntuple_entries);
953 
954  }
955 
956 }
957 
void beginJob() override
double GetPigtail(DAQHit const &daqhit) const
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
std::vector< sim::FLSHit > fHits
Definition: FLSHitList.h:21
enum BeamMode kRed
double minY
Definition: plot_hist.C:9
Double_t xx
Definition: macro.C:12
void LinFit(const std::vector< double > &x, const std::vector< double > &y, float *fitpar)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
float GetRes(int ADC)
Float_t x1[n_points_granero]
Definition: compare.C:5
float GetZ(int plane)
float GetExitT() const
Definition: FLSHit.h:57
int Mother() const
Definition: MCParticle.h:212
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
value_type val
Definition: BaseProducts.h:109
constexpr T pow(T x)
Definition: pow.h:75
list_type::const_iterator const_iterator
float GetX(int cell)
std::string _sliceInstanceLabel
instance label making the HitList
Definition: event.h:19
const double T0
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: EventID.h:98
std::string Process() const
Definition: MCParticle.h:214
double chi2()
std::map< std::string, TNtuple * > ntup_
TVector3 const & End() const
Definition: Track3D.h:46
int TrackId() const
Definition: MCParticle.h:209
Identifier for the Y measuring view of the detector (side)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
UpMuAna(fhicl::ParameterSet const &p)
int GetPDG() const
PDG.
Definition: FLSHit.h:43
void analyze(art::Event const &e) override
void hits()
Definition: readHits.C:15
length
Definition: demo0.py:21
TDC()=default
void GetXYZD(DAQHit const &daqhit, double w, double *xyzd)
unsigned int _nTracks
std::string _trackModuleLabel
label of module making the Tracks
const double a
novaddt::View const & View() const
Definition: DAQHit.h:72
T get(std::string const &key) const
Definition: ParameterSet.h:231
Identifier for the X measuring view of the detector (top)
Float_t d
Definition: plot.C:236
std::string _FLSHitListLabel
TVector3 const & Start() const
Definition: Track3D.h:45
std::string _sliceModuleLabel
label of module making the HitList
const double j
Definition: BetheBloch.cxx:29
value_type val
Definition: BaseProducts.h:84
int const & View() const
Definition: Track3D.h:44
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)
bool const & Is3D() const
Definition: Track3D.h:43
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
z
Definition: test.py:28
void GetXYZ(DAQHit const &daqhit, TVector3 start, TVector3 end, double *xyzd)
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
std::string _trackInstanceLabel
instance label making the Tracks
std::string _instance
void populateFLSmap(std::map< geo::OfflineChan, const sim::FLSHit * > &pcpflsmap, std::map< geo::OfflineChan, unsigned short > &pcpnummap, art::Event const &e)
EventNumber_t event() const
Definition: EventID.h:116
float GetEntryT() const
Definition: FLSHit.h:51
SubRunNumber_t subRun() const
Definition: EventID.h:110
std::string _trackToSliceInstanceLabel
T * make(ARGS...args) const
unsigned int _nSlices
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
float GetT(double TDC, double dist_to_readout)
A (plane, cell) pair.
Definition: OfflineChan.h:17
std::string _RawDigitLabel
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
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
enum BeamMode kGreen
Float_t w
Definition: plot.C:20
art::ServiceHandle< art::TFileService > tfs_
const TLorentzVector & EndMomentum() const
Definition: MCParticle.h:239
Definition: fwd.h:28
std::string _hitslabel
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_q.hpp:12
EventID id() const
Definition: Event.h:56
float RootFitLLR(std::vector< double > &x_hit, std::vector< double > &y_hit, std::vector< double > &y_error, float slope)
enum BeamMode string