UpMuAnalysis_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: UpMuAnalysis
3 // Module Type: analyzer
4 // File: UpMuAnalysis_module.cc
5 //
6 // Generated at Mon May 16 20:16:31 2016 by Aristeidis Tsaris using artmod
7 // from cetpkgsupport v1_08_07.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // Framework includes
20 #include "Utilities/func/MathUtil.h"
23 #include "fhiclcpp/ParameterSet.h"
25 
26 //NOVA includes
27 #include "RecoBase/Track.h"
28 #include "RecoBase/Cluster.h"
29 #include "RecoBase/RecoHit.h"
30 #include "RecoBase/Vertex.h"
31 #include <boost/math/special_functions.hpp> // for llr
33 #include "RawData/DAQHeader.h"
34 #include "RawData/RawTrigger.h"
35 #include "NovaTimingUtilities/TimingUtilities.h"
37 
38 // ROOT includes
39 #include "TVector3.h"
40 #include "TNtuple.h"
41 #include "TH1.h"
42 #include "TF1.h"
43 #include "TFormula.h"
44 #include "TGraph.h"
45 #include "TCanvas.h"
46 #include "TGraphErrors.h"
47 namespace upmuana {
48  class UpMuAnalysis;
49 
50  inline double getErr(double PE) {
51  return 165143/(1882.9+pow(PE,2.11447)) + 10.4321;
52  }
53 
54  inline double getDist(double x1, double y1, double z1,
55  double x2, double y2, double z2) {
56  return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
57  }
58 
59 
61  if (lhs.Y() == rhs.Y()) {
62  if (lhs.Z() == rhs.Z()) {
63  return lhs.T() < rhs.T();
64  }
65  else return lhs.Z() < rhs.Z();
66  }
67  else return lhs.Y() < rhs.Y();
68  }
69 
70  bool pairHitComp (std::pair<rb::CellHit, rb::RecoHit> lhp,
71  std::pair<rb::CellHit, rb::RecoHit> rhp) {
72  return recoHitComp(lhp.second, rhp.second);
73  }
74 
75 }
76 
78 public:
79  explicit UpMuAnalysis(fhicl::ParameterSet const & p);
80  // The destructor generated by the compiler is fine for classes
81  // without bare pointers or other resource use.
82 
83  // Plugins should not be copied or assigned.
84  UpMuAnalysis(UpMuAnalysis const &) = delete;
85  UpMuAnalysis(UpMuAnalysis &&) = delete;
86  UpMuAnalysis & operator = (UpMuAnalysis const &) = delete;
87  UpMuAnalysis & operator = (UpMuAnalysis &&) = delete;
88 
89  // Required functions.
90  void analyze(art::Event const & e) override;
91 
92  // Selected optional functions.
93  void beginJob() override;
94 
95 private:
96 
97  // Declare member data here.
101 
102  float containmentType(rb::Track);
103  double fMinY = -700.0, fMinX = -650.0, fMinZ = 100.0,
104  fMaxY = 400.0, fMaxX = 650.0, fMaxZ = 5800.0;
105 
106  double getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
107  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
108  std::pair<rb::CellHit,rb::RecoHit>)> ,
109  std::vector<rb::RecoHit> &outliers,
110  double &P_up, double &P_dn, double &chi2, double &slope);
111 
112  void LLR(std::vector<double>& eT,
113  std::vector<double>& mT,
114  std::vector<double>& mTErr, double& slope, double& chi2,
115  double& P_up, double& P_dn, std::vector<rb::RecoHit> &in,
116  std::vector<rb::RecoHit> &outliers);
117 
118  void LinFit(const std::vector<double>& x,
119  const std::vector<double>& y, double *fitpar);
120  void LinFit(const std::vector<double>& x,
121  const std::vector<double>& y,
122  const std::vector<double>& ye, double *fitpar);
123 
124  TVector3 AnglesToVector(double zen, double azi) const;
125 
128 
129  TNtuple* ntp_track;
130 
131 };
132 
133 
135  :
136  EDAnalyzer(p) // ,
137  // More initializers here.
138  , trackTag_(p.get<std::string>("trackInputTag"))
139  , fSliceModuleLabel (p.get<std::string>("SliceModuleLabel"))
140  , fSliceInstance (p.get<std::string>("SliceInstanceLabel"))
141  , ntp_track(nullptr)
142 {}
143 
145 {
146  // Implementation of optional member function here.
148  ntp_track = tfs->make<TNtuple>( "ntp_track", "Track Ntuple", "Run:SubRun:Event:SliceID:TrackID:Nhits:NRecohits:NOutliers:ProbUp:ProbDn:LLR:Chi2:Slope:LLRX:Chi2X:SlopeX:LLRY:Chi2Y:SlopeY:R2X:R2Y:StartX:StartY:StartZ:StartT:EndX:EndY:EndZ:EndT:TrackHitsX:TrackHitsY:Length:dirX:dirY:dirZ:eleAngle:totalE:containment:avgTX:avgTY:totalMSlices:SunZen:SunAzi:CosTheta:Azimuth:dotSun:zen_trk:azi_trk:event_time:chi2X_fit:chi2Y_fit:trackX_fromfit:trackY_fromfit:trackX_fromfit_slice:trackY_fromfit_slice");
149 
150 }
151 
153 {
154  // Get all slices
156  e.getByLabel( fSliceModuleLabel, slicecol);
158 
159  for(int i = 0; i < (int)slicecol->size(); i++){
160  art::Ptr<rb::Cluster> sli(slicecol, i);
161  slces.push_back(sli);
162  }
163 
164  // get all
166  e.getByLabel(trackTag_, tracks);
167  art::FindOneP<rb::Cluster> fmSlice(tracks, e, trackTag_);
168 
170  e.getByLabel("daq", raw_triggers); // this is need it
171  if (raw_triggers->size() != 1)
173  << __FILE__ << ":" << __LINE__ << "\n"
174  << "RawTrigger vector has incorrect size: " << raw_triggers->size()
175  << std::endl;
176 
177  std::vector<const rb::Cluster*> slices;
178  std::vector<float> slice_containment;
179 
180  art::Handle<rawdata::DAQHeader> raw_daqheader;
181  e.getByLabel("daq", raw_daqheader);
182 
183  double zen, azi;
184  unsigned long long event_time(raw_triggers->front().fTriggerTimingMarker_TimeStart);
185 
186  struct timespec ts;
188  convertNovaTimeToUnixTime(event_time, ts);
189 
190  fSunPos->GetSunPosition_FD(ts.tv_sec, zen, azi);
191 
192  //double zenR, aziR;
193  //fSunPos->GetRndmSunPosition_FD(ts, zenR, aziR);
194 
195  std::vector<unsigned> tracks_per_slice;
196  for (size_t i_track=0; i_track < tracks->size(); ++i_track){
197 
198  rb::Track theTrack = tracks->at(i_track);
199  float containment = containmentType(theTrack);
200 
201  float slice=-1;
202  bool found_slice = false;
203 
204  // Loop over Slices
205  const rb::Cluster *theSlice = &(*fmSlice.at(i_track));
206 
207  for (size_t i_slice=0; i_slice<slices.size(); ++i_slice) {
208  if (theSlice == slices.at(i_slice)){
209  slice = (float)i_slice; found_slice = true;
210  ++tracks_per_slice.at(i_slice);
211  //do the linear fit on the slice
212  if(containment != 4) slice_containment.at(i_slice) = 1;
213  break;
214  }
215  } //end loop slice
216 
217  if (!found_slice){
218  slice = (float)slices.size();
219  slices.push_back(theSlice);
220  tracks_per_slice.push_back(1);
221  if(containment!=4)
222  slice_containment.push_back(1);
223  else
224  slice_containment.push_back(4);
225  }
226 
227  // another loop to find slice id from vector slices
228  slice = -1;
229  for(size_t i_slice = 0; i_slice < slces.size(); i_slice++)
230  {
231  if(slces.at(i_slice).get() == theSlice)
232  {
233  slice = (float)i_slice;
234  break;
235  }
236  }
237 
238 
239  if (!theTrack.Is3D()) continue;
240 
241  art::PtrVector<rb::CellHit> trackHits = theTrack.AllCells();
242  std::set< std::pair<rb::CellHit,rb::RecoHit>,
243  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
244  std::pair<rb::CellHit,rb::RecoHit>)>
245  sortedTrackHits (pairHitComp);
246  std::set< std::pair<rb::CellHit,rb::RecoHit>,
247  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
248  std::pair<rb::CellHit,rb::RecoHit>)>
249  sortedTrackHitsX (pairHitComp);
250  std::set< std::pair<rb::CellHit,rb::RecoHit>,
251  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
252  std::pair<rb::CellHit,rb::RecoHit>)>
253  sortedTrackHitsY (pairHitComp);
254 
255  std::vector<double> x_hit;
256  std::vector<double> y_hit;
257  std::vector<double> zy_hit;
258  std::vector<double> zx_hit;
259 
260  // Loop over Hits
261  for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
262  art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
263  if (!theHit->GoodTiming()) continue;
264  rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
265  if (!theRecoHit.IsCalibrated()) continue;
266  sortedTrackHits.insert(std::make_pair(*theHit,theRecoHit));
267 
268  unsigned short iC = theHit->Cell();
269  unsigned short iP = theHit->Plane();
270  geo::View_t view = theHit->View();
271  if(view == geo::kY){
272  x_hit.push_back(iC);
273  zx_hit.push_back(iP);
274  sortedTrackHitsY.insert(std::make_pair(*theHit,theRecoHit));
275  }
276  else if(view == geo::kX){
277  y_hit.push_back(iC);
278  zy_hit.push_back(iP);
279  sortedTrackHitsX.insert(std::make_pair(*theHit,theRecoHit));
280  }
281  } // end loop hit
282 
283  double fitpar[3];
284  LinFit(x_hit, zx_hit, fitpar);
285  double r2x = fitpar[2];
286  LinFit(y_hit, zy_hit, fitpar);
287  double r2y = fitpar[2];
288  unsigned nxhit = x_hit.size();
289  unsigned nyhit = y_hit.size();
290 
291  //////////////////////////////////////////////////////////////////////////
292  //cp linear fit
293 
294  Double_t chi2X_fit = 1e6;
295  Double_t chi2Y_fit = 1e6;
296  const double wgt = 1;
297  Double_t trackX_fromfit = 0;
298  Double_t trackY_fromfit = 0;
299 
300  Double_t trackX_fromfit_slice = 0;
301  Double_t trackY_fromfit_slice = 0;
302 
303  double mvalX, cvalX, mvalY, cvalY;
304  TH1F *h1 = new TH1F("h1","xline",50,0,600);
305  TH1F *h2 = new TH1F("h2","yline",50,0,600);
306  // Loop over Slices
307  //for (size_t i = 0; i<slices.size(); ++i) {
308  std::vector<double> x_hit_fit;
309  std::vector<double> y_hit_fit;
310  std::vector<double> zy_hit_fit;
311  std::vector<double> zx_hit_fit;
312 
313  std::vector<double> weight_x;
314  std::vector<double> weight_y;
315 
316 
317  //const rb::Cluster* slice = theSlice;
318 
319  art::PtrVector<rb::CellHit> xcells = theSlice->XCells();
320  art::PtrVector<rb::CellHit> ycells = theSlice->YCells();
321 
322  for(size_t y_ind =0; y_ind < ycells.size(); y_ind++)
323  {
324  art::Ptr<rb::CellHit> theHity = ycells.at(y_ind);
325 
326  unsigned short iC = theHity->Cell();
327  unsigned short iP = theHity->Plane();
328  y_hit_fit.push_back(iC);
329  zy_hit_fit.push_back(iP);
330  weight_y.push_back(wgt);
331  trackY_fromfit += y_ind;
332  h1->Fill(y_hit_fit.at(y_ind));
333  }
334 
335  for(size_t x_ind = 0; x_ind < xcells.size(); x_ind++)
336  {
337  art::Ptr<rb::CellHit> theHitx = xcells.at(x_ind);
338 
339  unsigned short iC = theHitx->Cell();
340  unsigned short iP = theHitx->Plane();
341  x_hit_fit.push_back(iC);
342  zx_hit_fit.push_back(iP);
343  weight_x.push_back(wgt);
344  trackX_fromfit += x_ind;
345  h2->Fill(x_hit_fit.at(x_ind));
346  }
347  trackY_fromfit_slice += trackY_fromfit;
348  trackX_fromfit_slice += trackX_fromfit;
349  //}
350  if(x_hit_fit.size()>=2){
351  double chi2X_aux = util::LinFit(zx_hit_fit, x_hit_fit, weight_x, mvalX, cvalX);
352  if(chi2X_aux==chi2X_aux)chi2X_fit = chi2X_aux;
353  }
354  if(y_hit_fit.size()>=2){
355  double chi2Y_aux = util::LinFit(zy_hit_fit, y_hit_fit, weight_y, mvalY, cvalY);
356  if(chi2Y_aux==chi2Y_aux)chi2Y_fit = chi2Y_aux;
357  }
358 
359  ///////////////////////////////////////////////////////////
360  std::vector<rb::RecoHit> outliers, outliersX, outliersY;
361  double llr, P_up, P_dn, chi2, slope;
362  double llrX, P_upX, P_dnX, chi2X, slopeX;
363  double llrY, P_upY, P_dnY, chi2Y, slopeY;
364 
365  llr = getLLR(sortedTrackHits, outliers, P_up, P_dn, chi2, slope);
366  llrX = getLLR(sortedTrackHitsX, outliersX, P_upX, P_dnX, chi2X, slopeX);
367  llrY = getLLR(sortedTrackHitsY, outliersY, P_upY, P_dnY, chi2Y, slopeY);
368 
371  if (sortedTrackHits.size() >= 1) start = (sortedTrackHits.begin()->second);
372  else { start = rb::RecoHit(); end = rb::RecoHit(); }
373  if (sortedTrackHits.size() >= 2) end = (--sortedTrackHits.end())->second;
374  else end = start;
375  float dist = theTrack.TotalLength();
376 
377  // Average values
378  float tDirX=0, tDirY=0, tDirZ=0, tEleAngle, trackTotalE=0;
379  float avgTX=0, avgTY=0, tCosTheta, tAzimuth;
380 
381  for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
382  art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
383  rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
384 
385  TVector3 theDir = theTrack.InterpolateDir(theRecoHit.Z());
386 
387  tDirX += theDir.x()/trackHits.size();
388  tDirY += theDir.y()/trackHits.size();
389  tDirZ += theDir.z()/trackHits.size();
390 
391  if(theRecoHit.IsCalibrated()) trackTotalE += theRecoHit.GeV();
392 
393  if (theHit->View() == geo::kX)
394  avgTX += theRecoHit.T()/sortedTrackHitsX.size();
395  else if (theHit->View() == geo::kY)
396  avgTY += theRecoHit.T()/sortedTrackHitsY.size();
397  } // close the Hit loop
398 
399  tEleAngle = atan(tDirY/sqrt(tDirX*tDirX + tDirZ*tDirZ));
400  tCosTheta = fUtilities.CosTheta(theTrack.Dir().Y(), theTrack.Dir().Mag());
401  tAzimuth = fUtilities.Azimuth(theTrack.Dir().X(), theTrack.Dir().Z(), theTrack.Dir().Mag());
402 
403  TVector3 start_trk = theTrack.Start();
404  TVector3 end_trk = theTrack.Stop();
405  if(start_trk.Y() < end_trk.Y()) std::swap(start_trk, end_trk);
406  const TVector3 dir_trk = (start_trk-end_trk).Unit();
407  const double tzen_trk = acos(dir_trk.Y()) * 180 / M_PI;
408  double tazi_trk = atan2(-dir_trk.X(), dir_trk.Z()) * 180 / M_PI + 332.066111;
409  if(tazi_trk < 0) tazi_trk += 360;
410  if(tazi_trk > 360) tazi_trk -= 360;
411 
412  const double tdotSun = AnglesToVector(zen, azi).Dot(AnglesToVector(tzen_trk, tazi_trk));
413  //const double tdotSunR = AnglesToVector(zenR, aziR).Dot(AnglesToVector(tzen_trk, tazi_trk));
414 
415  float track_entries[55] =
416  {
417  (float)e.id().run(), (float)e.id().subRun(), (float)e.id().event(),
418  slice, (float)i_track, (float)trackHits.size(), (float)sortedTrackHits.size(),
419  (float)outliers.size(), (float)P_up, (float)P_dn, (float)llr,
420  (float)chi2, (float)slope, (float)llrX, (float)chi2X, (float)slopeX,
421  (float)llrY, (float)chi2Y, (float)slopeY, (float)r2x, (float)r2y,
422  start.X(), start.Y(), start.Z(), start.T(), end.X(), end.Y(),
423  end.Z(), end.T(), (float)nxhit,(float)nyhit, dist, tDirX, tDirY,
424  tDirZ, tEleAngle, trackTotalE, containment, avgTX, avgTY,
425  (float)raw_daqheader->TotalMicroSlices(), (float)zen, (float)azi,
426  tCosTheta, tAzimuth, (float)tdotSun, (float)tzen_trk, (float)tazi_trk,
427  (float)event_time, (float)chi2X_fit, (float)chi2Y_fit,(float)trackX_fromfit, (float)trackY_fromfit, (float)trackX_fromfit_slice, (float)trackY_fromfit_slice
428  };
429  ntp_track->Fill(track_entries);
430  } //end loop track
431 }
432 
433 // Creates unit vector pointing to a point on unit sphere with theta=zen, phi = azi
434 TVector3 upmuana::UpMuAnalysis::AnglesToVector(double zen, double azi) const
435 {
436  zen *= M_PI / 180;
437  azi *= M_PI / 180;
438  // Doesn't matter which axis is which
439  return TVector3(sin(zen)*cos(azi), sin(zen)*sin(azi), cos(zen));
440 }
441 
443 {
444  TVector3 start = theTrack.Start();
445  TVector3 stop = theTrack.Stop();
446 
447  if (start.y() > stop.y()) // assume upward-going, so switch
448  { TVector3 hold = start; start = stop; stop = hold; }
449 
450  if (start.y() < fMinY || start.x() < fMinX || start.x() > fMaxX ||
451  start.z() < fMinZ || start.z() > fMaxZ) { // entered from outside
452  if (stop.y() > fMaxY ||
453  stop.x() < fMinX ||
454  stop.x() > fMaxX ||
455  stop.z() < fMinZ ||
456  stop.z() > fMaxZ) { return 1; }//through-going
457  return 2; // stopped
458  } // started in detector
459  if (stop.y() > fMaxY || stop.x() < fMinX || stop.x() > fMaxX ||
460  stop.z() < fMinZ || stop.z() > fMaxZ)
461  { return 3; }// in-produced
462 
463  return 4; // fully contained
464 
465 }
466 
467 double upmuana::UpMuAnalysis::getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
468  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
469  std::pair<rb::CellHit,rb::RecoHit>)>
470  sortedHits,
471  std::vector<rb::RecoHit> &outliers,
472  double &P_up, double &P_dn, double &chi2,
473  double &slope)
474 {
475  if (sortedHits.size() < 2) {
476  P_up = 1e-30;
477  P_dn = 1e-30;
478  chi2 = 1e6;
479  slope = -1e6;
480  return 0;
481  }
482 
483  std::vector<double> measuredTimes;
484  std::vector<double> expectedTimes;
485  std::vector<double> wgts;
486 
487  // TGraph *g_plot = new TGraph (columns.size(), &columns[0], &columns[1]);
488  // g->SetMarkerStyle(21);
489  // g->Draw("alp");
490  // first hit
491  measuredTimes.push_back(0.0);
492  expectedTimes.push_back(0.0);
493  double err = getErr(sortedHits.begin()->first.PE());
494  wgts.push_back(1.0/err/err);
495 
496  double startX = sortedHits.begin()->second.X(),
497  startY = sortedHits.begin()->second.Y(),
498  startZ = sortedHits.begin()->second.Z(),
499  startT = sortedHits.begin()->second.T();
500 
501  std::vector<rb::RecoHit> in;
502 
503  for (auto i_hit = ++sortedHits.begin();
504  i_hit != sortedHits.end();
505  ++i_hit) {
506  in.push_back(i_hit->second);
507  measuredTimes.push_back(i_hit->second.T() - startT);
508  double dist = getDist(i_hit->second.X(), i_hit->second.Y(),
509  i_hit->second.Z(), startX, startY, startZ);
510  expectedTimes.push_back(dist/29.97);
511  err = getErr(i_hit->first.PE());
512  wgts.push_back(1.0/err/err);
513  }
514 
515  LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn,
516  in, outliers);
517 
518  return log(P_up/P_dn);
519 }
520 
521 void upmuana::UpMuAnalysis::LLR(std::vector<double>& eT,
522  std::vector<double>& mT,
523  std::vector<double>& wgts, double& slope,
524  double& chi2, double& P_up, double& P_dn,
525  std::vector<rb::RecoHit>& in,
526  std::vector<rb::RecoHit>& outliers)
527 {
528  // eT - x, mT - y
529  size_t n = mT.size();
530 
531  // y = a + bx
532  double a, b;
533  if (eT.size() >= 2)
534  chi2 = util::LinFit(eT, mT, wgts, b, a);
535  else {
536  chi2 = 1e6;
537  slope = -1e6;
538  P_up = 1e-30;
539  P_dn = 1e-30;
540  return;
541  }
542 
543  size_t totAllowOutlier = n/10;
544  size_t nOutliers = 0;
545  std::vector<double> x_filt, y_filt, w_filt;
546  for (size_t i=0; i < n; i++) {
547  double y_fit = a + b*eT.at(i);
548  if ( fabs(mT.at(i) - y_fit) < 5*(1.0/sqrt(wgts.at(i)))
549  || nOutliers>totAllowOutlier)
550  {
551  x_filt.push_back(eT.at(i));
552  y_filt.push_back(mT.at(i));
553  w_filt.push_back(wgts.at(i));
554  }
555  else {
556  ++nOutliers;
557  outliers.push_back(in[i]);
558  }
559  }
560 
561  if (x_filt.size() >= 2)
562  chi2 = util::LinFit(x_filt, y_filt, w_filt, b, a);
563  else {
564  chi2 = 1e6;
565  slope = -1e6;
566  P_up = 1e-30;
567  P_dn = 1e-30;
568  }
569  n = x_filt.size();
570  if (n < 5) {
571  slope = 0;
572  chi2 = 999;
573  P_up = 1e-30;
574  P_dn = 1;
575  return;
576  }
577 
578  slope = b;
579  chi2 /= (double)(n-2);
580 
581  double one_over_ee = 0.0,
582  x_over_ee = 0.0,
583  y_over_ee = 0.0;
584 
585  for (size_t i=0; i<n; ++i) {
586  one_over_ee += w_filt.at(i);
587  x_over_ee += x_filt.at(i)*w_filt.at(i);
588  y_over_ee += y_filt.at(i)*w_filt.at(i);
589  }
590 
591  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
592  double dn_inter = (y_over_ee+x_over_ee)/one_over_ee;
593 
594  double up_chi2 = 0.0, dn_chi2 = 0.0;
595  for (size_t i=0; i<n; ++i) {
596  double e = 1.0/sqrt(w_filt.at(i));
597  double up_expected = up_inter + x_filt.at(i);
598  double dn_expected = dn_inter - x_filt.at(i);
599  up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
600  dn_chi2 += pow((y_filt.at(i)-dn_expected) / e, 2.0);
601  }
602 
603  double prob_up = boost::math::gamma_q((double)(n-2)/2.0,up_chi2/2.0),
604  prob_dn = boost::math::gamma_q((double)(n-2)/2.0,dn_chi2/2.0);
605 
606  if (prob_up < 1e-30) prob_up = 1e-30;
607  if (prob_dn < 1e-30) prob_dn = 1e-30;
608 
609  P_up = prob_up;
610  P_dn = prob_dn;
611 
612 }
613 
614 
615 void upmuana::UpMuAnalysis::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
616 
617  const auto n = x_val.size();
618  const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n; // <x>
619  const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n; // <y>
620  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n; // <xx>
621  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n; // <yy>
622  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n; // <xy>
623 
624  const auto b = (xy - x*y)/(xx - x*x); // slope
625  const auto a = y - b*x; // intercept
626  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
627  fitpar[0] = a;
628  fitpar[1] = b;
629  fitpar[2] = r*r;
630 
631 }
632 
633 void upmuana::UpMuAnalysis::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, const std::vector<double>& y_err, double *fitpar){
634 
635  int n = x_val.size();
636  double x = 0;
637  double y = 0;
638  double xx = 0;
639  double yy = 0;
640  double xy = 0;
641  double ee = 0;
642 
643  for ( int i=0; i<n; i++ ){
644  x = x + x_val[i]/y_err[i]/y_err[i];
645  y = y + y_val[i]/y_err[i]/y_err[i];
646 
647  xx = xx + x_val[i]*x_val[i]/y_err[i]/y_err[i];
648  yy = yy + y_val[i]*y_val[i]/y_err[i]/y_err[i];
649  xy = xy + x_val[i]*y_val[i]/y_err[i]/y_err[i];
650  ee = ee + 1./y_err[i]/y_err[i];
651  }
652 
653  const auto b = (xy*ee - x*y)/(xx*ee - x*x); // slope
654  const auto a = (xy - b*xx)/x; // intercept
655  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
656  fitpar[0] = a;
657  fitpar[1] = b;
658  fitpar[2] = r*r;
659 
660 
661 }
662 
float T() const
Definition: RecoHit.h:39
Double_t xx
Definition: macro.C:12
const art::PtrVector< rb::CellHit > & XCells() const
Get all cells from the x-view.
Definition: Cluster.h:124
float containmentType(rb::Track)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
UpMuAnalysis & operator=(UpMuAnalysis const &)=delete
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
trk::CosmicTrackUtilities fUtilities
Float_t x1[n_points_granero]
Definition: compare.C:5
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
A collection of associated CellHits.
Definition: Cluster.h:47
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: EventID.h:98
virtual TVector3 Start() const
Definition: Prong.h:73
int TotalMicroSlices() const
Definition: DAQHeader.h:27
T acos(T number)
Definition: d0nt_math.hpp:54
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
double chi2()
double getErr(double PE)
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
double dist
Definition: runWimpSim.h:113
void LLR(std::vector< double > &eT, std::vector< double > &mT, std::vector< double > &mTErr, double &slope, double &chi2, double &P_up, double &P_dn, std::vector< rb::RecoHit > &in, std::vector< rb::RecoHit > &outliers)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
T atan(T number)
Definition: d0nt_math.hpp:66
Locate positions of celestial bodies.
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
UpMuAnalysis(fhicl::ParameterSet const &p)
unsigned short Cell() const
Definition: CellHit.h:40
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double getLLR(std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>, std::vector< rb::RecoHit > &outliers, double &P_up, double &P_dn, double &chi2, double &slope)
const double a
bool recoHitComp(rb::RecoHit lhs, rb::RecoHit rhs)
const art::PtrVector< rb::CellHit > & YCells() const
Get all cells from the x-view.
Definition: Cluster.h:126
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
float Azimuth(float const &dcosX, float const &dcosZ, float const &magnitude)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
float CosTheta(float const &dcosY, float const &magnitude)
reference at(size_type n)
Definition: PtrVector.h:365
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void GetSunPosition_FD(time_t time, double &zen, double &azi)
TVector3 Unit() const
Vertex location in position and time.
const ana::Var wgt
size_type size() const
Definition: PtrVector.h:308
TH1F * h2
Definition: plot.C:45
EventNumber_t event() const
Definition: EventID.h:116
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
TH1F * h1
bool pairHitComp(std::pair< rb::CellHit, rb::RecoHit > lhp, std::pair< rb::CellHit, rb::RecoHit > rhp)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
SubRunNumber_t subRun() const
Definition: EventID.h:110
T * make(ARGS...args) const
ifstream in
Definition: comparison.C:7
float GeV() const
Definition: RecoHit.cxx:69
T sin(T number)
Definition: d0nt_math.hpp:132
TVector3 AnglesToVector(double zen, double azi) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void analyze(art::Event const &e) override
float X() const
Definition: RecoHit.h:36
const hit & b
Definition: hits.cxx:21
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
virtual TVector3 InterpolateDir(double z) const
Definition: Track.cxx:348
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
float Y() const
Definition: RecoHit.h:37
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
bool GoodTiming() const
Definition: CellHit.h:48
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
double getDist(double x1, double y1, double z1, double x2, double y2, double z2)
Float_t e
Definition: plot.C:35
art::ServiceHandle< locator::CelestialLocator > fSunPos
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_q.hpp:12
size_t size() const
EventID id() const
Definition: Event.h:56
virtual bool Is3D() const
Definition: Prong.h:71
T atan2(T number)
Definition: d0nt_math.hpp:72