FmmTrackerValidation_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FmmTrackerValidation
3 // Module Type: analyzer
4 // File: FmmTrackerValidation_module.cc
5 //
6 // Generated at Wed Dec 18 15:41:32 2013 by Zukai Wang using artmod
7 // from cetpkgsupport v1_04_02.
8 ////////////////////////////////////////////////////////////////////////
9 #include "art/Framework/Core/EDAnalyzer.h"
10 #include "art/Framework/Core/ModuleMacros.h"
11 #include "art/Framework/Principal/Event.h"
12 #include "art/Framework/Principal/Handle.h"
13 #include "art/Framework/Principal/Run.h"
14 
15 #include "MCCheater/BackTracker.h"
16 #include "Calibrator/Calibrator.h"
17 #include "Simulation/Particle.h"
18 #include "Simulation/ParticleNavigator.h"
19 #include "art/Framework/Principal/SubRun.h"
20 #include "canvas/Utilities/InputTag.h"
21 #include "Utilities/func/MathUtil.h"
22 #include "fhiclcpp/ParameterSet.h"
23 #include "art/Framework/Services/Optional/TFileService.h"
24 
25 #include "Geometry/Geometry.h"
26 #include "GeometryObjects/CellGeo.h"
27 #include "GeometryObjects/Geo.h"
28 #include "GeometryObjects/PlaneGeo.h"
29 
30 #include "NovaDAQConventions/DAQConventions.h"
31 #include "RecoBase/CellHit.h"
32 #include "RecoBase/RecoHit.h"
33 #include "RecoBase/Cluster.h"
34 #include "RecoBase/Track.h"
35 
36 #include "FastMonopoleTriggers.h"
37 
38 #include "TVector3.h"
39 #include "TTree.h"
40 #include "TProfile.h"
41 #include "TRandom3.h"
42 
43 #include <algorithm>
44 #include <utility>
45 #include <cmath>
46 #include <iostream>
47 
48 namespace zcl {
49  class FmmTrackerValidation;
50 }
51 
52 class zcl::FmmTrackerValidation : public art::EDAnalyzer {
53  public:
54  explicit FmmTrackerValidation(fhicl::ParameterSet const & pset);
55  virtual ~FmmTrackerValidation();
56 
57  void analyze(const art::Event & evt);
58  void beginJob();
59  void endJob();
60 
61  private:
62  std::string fTrackerInput;
63  int fSaveImage;
64  const int ntrials = 100; //Number of MC trilas in determining the fitting error
65  const double X_err = 1.8; //System error used to calculate uncertainty in velocity recon.
66  const double Y_err = 1.8;
67  const double Z_err = 2.82;
68 
69  //Each entry is corresponding to 1 reconstructed track
70  TTree * RCV_Tree;
71  TRandom * genR = new TRandom3();
72  //Purely MC info from MCTruth
73  double beta_MC, Edep_MC, Edep_Birks_MC;
74  double xi_MC, xf_MC, yi_MC, yf_MC, zi_MC, zf_MC, ti_MC, tf_MC;
75  double x0, y0, z0, cosx, cosy, cosz, trlMC;
76  int nxhits_MC, nyhits_MC;
77 
78  //RC info
79  int triggered, trigger_penetrated, nxhits, nyhits, nhits_sat, nhits_true;
80  double trigger_nSatHits, trigger_meanADC, MIP_tot, ADC_tot, ADC_stdev, GeV_tot;
81  double tns_mean, tns_rms, tns_min, tns_max;
82  int max_gap, miss_plane, p_cross, p_overlap;
83  double Xi, Yi, Zi, Ti, Xf, Yf, Zf, Tf;//Notice the start and end point in the fmmtrack might be flipped.
84  double t_min, t_max, x_tmin, x_tmax, y_tmin, y_tmax, z_tmin, z_tmax;
85  double lxsqr, lysqr;
86  double vrc_sat, vSat_err, tSat_mean, tSat_rms;//velocity in cm/ns with saturated hits.
87  double vrc_nonsat, vnonSat_err, tnonSat_mean, tnonSat_rms;//velocity in cm/nswith non saturated hits.
88  double Chisqr_TS, Chisqr_TN;
89  double adc_W, GeV_W, r2xt, r2yt, r2GeVx, r2GeVy, time_gap_xz, time_gap_yz;
90  std::vector<int> hitxCell, hitxPlane, hityCell,hityPlane;
91  std::vector<double> hitxx, hitxz, hityy, hityz, hitxe, hitxt, hitye, hityt;
92 
93  double dist2(double dx, double dy);
94  double doca2(double x, double y, double x0, double y0, double dydx);
95 
96  //Linear Regression function
97  double vfit(std::vector<TVector3> & Positions,
98  std::vector<double> & Ts,
99  std::vector<double> & Terr,
100  double & v_err,
101  double & Tmean,
102  double & Trms,
103  double & dTsqr);
104  void LinFit(const std::vector<double>& x, const std::vector<double>& y, double *fitpar);
105 };
106 
107 zcl::FmmTrackerValidation::FmmTrackerValidation(fhicl::ParameterSet const & p) :
108  EDAnalyzer(p),
109  fTrackerInput (p.get< std::string > ("TrackerInput") ),
110  fSaveImage (p.get< int > ("SaveImage") )
111 {
112 }
113 
114 zcl::FmmTrackerValidation::~FmmTrackerValidation()
115 {
116  // Clean up dynamic memory and other resources here.
117 }
118 
119 void zcl::FmmTrackerValidation::beginJob()
120 {
121  art::ServiceHandle<art::TFileService> tfs;
122 
123  //Combining cheater info and the RC info to make the data structure better
124  //One RC track info per entry, if no tracks reconstructed in this event, fill in 1 entry anyway tagging the
125  //reconstruction failed
126  RCV_Tree = tfs->make<TTree>("RCV_Tree","Information from RC track info");
127 
128  RCV_Tree -> Branch("betaMC", &beta_MC, "betaMC/D");
129  RCV_Tree -> Branch("edepMC", &Edep_MC, "edepMC/D");
130  RCV_Tree -> Branch("edepbirksMC", &Edep_Birks_MC, "edepbirksMC/D");
131  RCV_Tree -> Branch("x0",&x0,"x0/D");//Entering point on surface of detector
132  RCV_Tree -> Branch("y0",&y0,"y0/D");//Entering point on surface of detector
133  RCV_Tree -> Branch("z0",&z0,"z0/D");//Entering point on surface of detector
134  RCV_Tree -> Branch("cosx",&cosx,"cosx/D");
135  RCV_Tree -> Branch("cosy",&cosy,"cosy/D");
136  RCV_Tree -> Branch("cosz",&cosz,"cosz/D");
137  RCV_Tree -> Branch("tiMC", &ti_MC, "tiMC/D");
138  RCV_Tree -> Branch("xiMC", &xi_MC, "xiMC/D");
139  RCV_Tree -> Branch("yiMC", &yi_MC, "yiMC/D");
140  RCV_Tree -> Branch("ziMC", &zi_MC, "ziMC/D");
141  RCV_Tree -> Branch("tfMC", &tf_MC, "tfMC/D");
142  RCV_Tree -> Branch("xfMC", &xf_MC, "xfMC/D");
143  RCV_Tree -> Branch("yfMC", &yf_MC, "yfMC/D");
144  RCV_Tree -> Branch("zfMC", &zf_MC, "zfMC/D");
145  RCV_Tree -> Branch("xhitsMC", &nxhits_MC, "xhitsMC/I");
146  RCV_Tree -> Branch("yhitsMC", &nyhits_MC, "yhitsMC/I");
147  RCV_Tree -> Branch("trlMC", &trlMC, "trlMC/D");
148  //Following is some trigger information.
149  RCV_Tree -> Branch("triggered", &triggered, "triggered/I");
150  RCV_Tree -> Branch("trigger_penetrated", &trigger_penetrated, "trigger_penetrated/I");
151  RCV_Tree -> Branch("trigger_nSatHits", &trigger_nSatHits, "trigger_nSatHits/D");
152  RCV_Tree -> Branch("trigger_meanADC", &trigger_meanADC, "trigger_meanADC/D");
153  //More recon information.
154  RCV_Tree -> Branch("nxhits", &nxhits, "nxhits/I");
155  RCV_Tree -> Branch("nyhits", &nyhits, "nyhits/I");
156  RCV_Tree -> Branch("nhits_true", &nhits_true, "nhits_true/I");
157  RCV_Tree -> Branch("nhits_sat", &nhits_sat, "nhits_sat/I");
158  RCV_Tree -> Branch("MIP_tot", &MIP_tot, "MIP_tot/D");
159  RCV_Tree -> Branch("ADC_tot", &ADC_tot, "ADC_tot/D");
160  RCV_Tree -> Branch("ADC_stdev", &ADC_stdev, "ADC_stdev/D");
161  RCV_Tree -> Branch("GeV_tot", &GeV_tot, "GeV_tot/D");
162  //Missing Planes and the max gap:
163  RCV_Tree -> Branch("max_gap", &max_gap, "max_gap/I");
164  RCV_Tree -> Branch("miss_plane", &miss_plane, "miss_plane/I");
165  RCV_Tree -> Branch("p_cross", &p_cross, "p_cross/I");
166  RCV_Tree -> Branch("p_overlap", &p_overlap, "p_overlap/I");
167  RCV_Tree -> Branch("Xi", &Xi, "Xi/D");
168  RCV_Tree -> Branch("Xf", &Xf, "Xf/D");
169  RCV_Tree -> Branch("Yi", &Yi, "Yi/D");
170  RCV_Tree -> Branch("Yf", &Yf, "Yf/D");
171  RCV_Tree -> Branch("Zi", &Zi, "Zi/D");
172  RCV_Tree -> Branch("Zf", &Zf, "Zf/D");
173  RCV_Tree -> Branch("Ti", &Ti, "Ti/D");
174  RCV_Tree -> Branch("Tf", &Tf, "Tf/D");
175  RCV_Tree -> Branch("tns_mean", &tns_mean, "tns_mean/D");
176  RCV_Tree -> Branch("tns_rms", &tns_rms, "tns_rms/D");
177  RCV_Tree -> Branch("tns_min", &tns_min, "tns_min/D");
178  RCV_Tree -> Branch("tns_max", &tns_max, "tns_max/D");
179  //This is double checking the velocity check:
180  RCV_Tree -> Branch("t_min", &t_min, "t_min/D");
181  RCV_Tree -> Branch("t_max", &t_max, "t_max/D");
182  RCV_Tree -> Branch("x_tmin", &x_tmin, "x_tmin/D");
183  RCV_Tree -> Branch("x_tmax", &x_tmax, "x_tmax/D");
184  RCV_Tree -> Branch("y_tmin", &y_tmin, "y_tmin/D");
185  RCV_Tree -> Branch("y_tmax", &y_tmax, "y_tmax/D");
186  RCV_Tree -> Branch("z_tmin", &z_tmin, "z_tmin/D");
187  RCV_Tree -> Branch("z_tmax", &z_tmax, "z_tmax/D");
188  //fitting parameters:
189  RCV_Tree -> Branch("lxsqr", &lxsqr, "lxsqr/D");
190  RCV_Tree -> Branch("lysqr", &lysqr, "lysqr/D");
191  RCV_Tree -> Branch("vrc_sat", &vrc_sat, "vrc_sat/D");
192  RCV_Tree -> Branch("tSat_mean", &tSat_mean,"tSat_mean/D");
193  RCV_Tree -> Branch("tSat_rms", &tSat_rms, "tSat_rms/D");
194  RCV_Tree -> Branch("vSat_err", &vSat_err, "vsat_err/D");
195  RCV_Tree -> Branch("vrc_nonsat", &vrc_nonsat, "vrc_nonsat/D");
196  RCV_Tree -> Branch("tnonSat_mean", &tnonSat_mean,"tnonSat_mean/D");
197  RCV_Tree -> Branch("tnonSat_rms", &tnonSat_rms, "tnonSat_rms/D");
198  RCV_Tree -> Branch("vnonSat_err", &vnonSat_err, "vnonSat_err/D");
199  RCV_Tree -> Branch("Chisqr_TS", &Chisqr_TS, "Chisqr_TS/D");
200  RCV_Tree -> Branch("Chisqr_TN", &Chisqr_TN, "Chisqr_TN/D");
201  RCV_Tree -> Branch("adc_W", &adc_W, "adc_W/D");
202  RCV_Tree -> Branch("GeV_W", &GeV_W, "GeV_W/D");
203  RCV_Tree->Branch("r2xt", &r2xt, "r2xt/D");//Regression coefficient for xt
204  RCV_Tree->Branch("r2yt", &r2yt, "r2yt/D");//Regression coefficient for yt
205  RCV_Tree->Branch("r2GeVx", &r2GeVx, "r2GeVx/D");//Regression coefficient for ex
206  RCV_Tree->Branch("r2GeVy", &r2GeVy, "r2GeVy/D");//Regression coefficient for ey
207  RCV_Tree->Branch("time_gap_xz", &time_gap_xz, "time_gap_xz/D");//Time gap for xz view, [0,1].
208  RCV_Tree->Branch("time_gap_yz", &time_gap_yz, "time_gap_yz/D");//Time gap for yz view, [0,1].
209  if(fSaveImage!=0){
210  //Some additional features added by Enhao to play with u-net segmentation.
211  RCV_Tree->Branch("hitxe", &hitxe);//vector with all x hit e
212  RCV_Tree->Branch("hitxt", &hitxt);//vector with all x hit t
213  RCV_Tree->Branch("hitye", &hitye);//vector with all y hit e
214  RCV_Tree->Branch("hityt", &hityt);//vector with all y hit t
215  RCV_Tree->Branch("hitxCell", &hitxCell);//hit cell number in xz view
216  RCV_Tree->Branch("hitxPlane", &hitxPlane);//hit plane number in xz view
217  RCV_Tree->Branch("hityCell", &hityCell);//hit cell number in yz view
218  RCV_Tree->Branch("hityPlane", &hityPlane);//hit cell number in yz view
219  }
220 }
221 
222 void zcl::FmmTrackerValidation::analyze(const art::Event & evt)
223 {
224 
225  art::ServiceHandle<geo::Geometry> geom;
226  art::ServiceHandle<cheat::BackTracker> bt;
227  art::ServiceHandle<calib::Calibrator> calib;
228 
229  art::Handle< std::vector<rb::Cluster> > slices;
230  evt.getByLabel("fmmslicer", slices);
231 
232  //Get MC info first:
233  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
234  nxhits_MC = 0;
235  nyhits_MC = 0;
236  trlMC = 0;
237  triggered = 0;
238  trigger_meanADC = 0;
239  trigger_nSatHits = 0;
240  trigger_penetrated = 0;
241  x0 = -9999;
242  y0 = -9999;
243  z0 = -9999;
244  cosx = -10;
245  cosy = -10;
246  cosz = -10;
247  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
248  nxhits_MC = 0;
249  nyhits_MC = 0;
250  trlMC = 0;
251 
252  const sim::Particle *p = (*i).second;
253  if (p->PdgCode()!=42 ) continue;
254  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId(), 42);
255  //const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId()); //This will include delta ray hits.
256  double momentum = p->P();
257  x0 = p->Vx();
258  y0 = p->Vy();
259  z0 = p->Vz();
260  cosx = (p->Px())/momentum;
261  cosy = (p->Py())/momentum;
262  cosz = (p->Pz())/momentum;
263 
264  double momentum2 = (p->P())*(p->P());
265  double mass2 = (p->Mass())*(p->Mass());
266  beta_MC = sqrt(momentum2 / (momentum2+mass2) );
267 
268  unsigned ntruehits_Tot = flshits.size();
269  if ( ntruehits_Tot==0 ) break; // There is only one monopole. If the monopole left no hits in detector... we can just fill MC infor and skip this event.
270  ti_MC = flshits[0].GetEntryT();
271  xi_MC = flshits[0].GetEntryX();
272  yi_MC = flshits[0].GetEntryY();
273  zi_MC = flshits[0].GetEntryZ();
274  tf_MC = flshits[0].GetExitT();
275  xf_MC = flshits[0].GetExitX();
276  yf_MC = flshits[0].GetExitY();
277  zf_MC = flshits[0].GetExitZ();
278  Edep_MC = flshits[0].GetEdep();
279  Edep_Birks_MC = flshits[0].GetEdepBirks();
280  trlMC = flshits[0].GetTotalPathLength();
281 
282  unsigned int planeID_i = flshits[0].GetPlaneID();
283  unsigned int cellID_i = flshits[0].GetCellID();
284  unsigned int planeID_f = planeID_i;
285  unsigned int cellID_f = cellID_i;
286 
287  if (planeID_i%2) ++nxhits_MC;
288  else ++nyhits_MC;
289 
290  for (unsigned h =1; h!= flshits.size(); ++h) {
291  double ti = flshits[h].GetEntryT();
292  double tf = flshits[h].GetExitT();
293  if (ti < ti_MC) {
294  ti_MC = ti;
295  xi_MC = flshits[h].GetEntryX();
296  yi_MC = flshits[h].GetEntryY();
297  zi_MC = flshits[h].GetEntryZ();
298  planeID_i = flshits[h].GetPlaneID();
299  cellID_i = flshits[h].GetCellID();
300  }
301  if (tf > tf_MC) {
302  tf_MC = tf;
303  xf_MC = flshits[h].GetExitX();
304  yf_MC = flshits[h].GetExitY();
305  zf_MC = flshits[h].GetExitZ();
306  planeID_f = flshits[h].GetPlaneID();
307  cellID_f = flshits[h].GetCellID();
308  }
309  if (flshits[h].GetPlaneID() % 2) ++nxhits_MC;
310  else ++nyhits_MC;
311 
312  Edep_MC += flshits[h].GetEdep();
313  Edep_Birks_MC += flshits[h].GetEdepBirks();
314  trlMC += flshits[h].GetTotalPathLength();
315  }
316 
317  //By now, you have got the local coordinate, you have to translate it to world:
318  const geo::CellGeo* startCell = geom->Plane(planeID_i)->Cell(cellID_i);
319  double local[3], world[3];
320  local[0] = xi_MC;
321  local[1] = yi_MC;
322  local[2] = zi_MC;
323  startCell->LocalToWorld(local,world);
324  xi_MC = world[0];
325  yi_MC = world[1];
326  zi_MC = world[2];
327 
328  const geo::CellGeo* stopCell = geom->Plane(planeID_f)->Cell(cellID_f);
329  local[0] = xf_MC;
330  local[1] = yf_MC;
331  local[2] = zf_MC;
332  stopCell->LocalToWorld(local,world);
333  xf_MC = world[0];
334  yf_MC = world[1];
335  zf_MC = world[2];
336  //1 monopole per event, so you can break out from the loop:
337  break;
338  }
339  //=====================================================================================
340  art::Handle< std::vector<rb::Track> > trks;
341  evt.getByLabel(fTrackerInput,trks);
342  nxhits = 0;
343  nyhits = 0;
344  nhits_sat = 0;
345  nhits_true = 0;
346  ADC_tot = 0;
347  ADC_stdev = 0;
348  //prepare vectors to save hit info
349  hitxe.clear();
350  hitxt.clear();
351  hitye.clear();
352  hityt.clear();
353  hitxx.clear();
354  hityy.clear();
355  hitxCell.clear();
356  hitxPlane.clear();
357  hityCell.clear();
358  hityPlane.clear();
359 
360  tns_mean = 0; //ADC Weighed average time
361  zcl::FastMonopoleTriggers trigger;
362  if(trigger.is_trigger_by_epoch1_fmmtrigger(slices, trigger_penetrated, trigger_nSatHits, trigger_meanADC )==false) {
363  RCV_Tree->Fill();
364  return;
365  }
366  triggered = 1;
367  //So events failed trigger or passed trigger without track get Filled once with MC info and zero RC info. Events passed trigger with tracks, will get the longest track filled RC info.
368 
369  bool filled = false;
370  //find the track with most true hits
371  int max_nhits_true = 0;
372  int trakIdx_for_max_nhits_true = -1;
373  for (unsigned trkIdx = 0; trkIdx != trks->size(); ++trkIdx ) {
374  nhits_true = 0;
375  const rb::Track& track = (*trks)[trkIdx];
376 
377  for (unsigned i = 0; i != track.NCell(); ++i) {
378  const art::Ptr<rb::CellHit> chit = track.Cell(i);
379  if(chit->IsMC()==true) ++nhits_true;
380  }
381  if(nhits_true>max_nhits_true){
382  max_nhits_true = nhits_true;
383  trakIdx_for_max_nhits_true=trkIdx;
384  }
385  }
386 
387  //take the entire 3D track:
388  if (trakIdx_for_max_nhits_true>-1) {
389  nhits_true = 0;
390  const rb::Track& track = (*trks)[trakIdx_for_max_nhits_true];
391  double ADC_sum = 0;
392  for (unsigned i = 0; i != track.NCell(); ++i) {
393  const art::Ptr<rb::CellHit> chit = track.Cell(i);
394  if(chit->IsMC()==true) ++nhits_true;
395  double tns = chit->TNS();
396  double adc = chit->ADC();
397  tns_mean += tns*adc;
398  ADC_sum += adc;
399  }
400  if (track.NCell()>=6){
401  //if(nhits_true*1.0/(nxhits_MC*1.0+nyhits_MC*1.0)<0.5) continue; // Oct 26 2017 Enhao, Feb 2018, this has some problem, for beta>0.1, some high velocity monopoles, they have a lot of hits didn't get included, hits from flash effect? We should use the longest track. Maybe there is better solution?
402  tns_mean /= ADC_sum;
403  //=============================================================
404  //Just determine the RMS of the time:
405  tns_rms = 0;
406  for (unsigned i = 0; i != track.NCell(); ++i) {
407  const art::Ptr<rb::CellHit> chit = track.Cell(i);
408  double tns = chit->TNS();
409  double adc = chit->ADC();
410  tns_rms = (tns-tns_mean)*adc*(tns-tns_mean)*adc;
411  }
412  tns_rms = sqrt(tns_rms/ADC_sum+250000/12.);//where does this 250000 come from?
413  //==============================================================
414  //Now, loop through the hits again to throw away the hits outside the 4 sigma range,
415  //Make sure they get updated!
416  tns_min = 99999999;
417  tns_max = -99999999;
418  nxhits = 0;
419  nyhits = 0;
420  nhits_sat = 0;
421  nhits_true = 0;
422  ADC_tot = 0;
423  std::vector<int> planes;
424 
425  int P_x0=9999;
426  int P_x1=0;
427  int P_y0=9999;
428  int P_y1=0;
429 
430  for (unsigned i = 0; i != track.NCell(); ++i) {
431  const art::Ptr<rb::CellHit> chit = track.Cell(i);
432  double tns = chit->TNS();
433  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
434  if (tns > tns_max) tns_max = tns;
435  if (tns < tns_min) tns_min = tns;
436  planes.emplace_back(chit->Plane());
437  ADC_tot += chit->ADC();
438 
439  if (chit->View()==geo::kX) {
440  ++nxhits;
441  if (chit->Plane() < P_x0) P_x0 = chit->Plane();
442  if (chit->Plane() > P_x1) P_x1 = chit->Plane();
443  }
444  else {
445  ++nyhits;
446  if (chit->Plane() < P_y0) P_y0 = chit->Plane();
447  if (chit->Plane() > P_y1) P_y1 = chit->Plane();
448  }
449  }
450 
451  if (planes.size()>=3 && P_x1-P_x0>=3 && P_y1-P_y0>=3) {
452 
453  //Now calculate the missing planes and the max gap:
454  std::sort (planes.begin(), planes.end());
455  max_gap = 0;
456  miss_plane = 0;
457  p_cross = planes.back()-planes.front()+1;
458  p_overlap = 1;
459  bool start_flag = false;
460  for (unsigned pidx =0; pidx+1<planes.size(); ++pidx) {
461  int dP = planes[pidx+1] - planes[pidx] - 1;
462  //Calculate the number of planes in the overlapped region:
463  if (!start_flag) {
464  if ( (planes[pidx+1] - planes[pidx]) % 2 ) start_flag = true;
465  }
466  else {
467  if ( (planes[pidx+1] - planes[pidx]) % 2 ) ++p_overlap;
468  }
469  if (dP > 0 ) {
470  miss_plane += dP;
471  if (dP > max_gap) max_gap = dP;
472  }
473  }
474  //========================================================
475  //Loop through the hits again: this time, determine dxdz and dydz
476  std::vector<double> x0s_x, x1s_x, y0s_y, y1s_y;
477  std::vector<double> x0s_z, x1s_z, y0s_z, y1s_z;
478  for (unsigned i = 0; i != track.NCell(); ++i) {
479  const art::Ptr<rb::CellHit> chit = track.Cell(i);
480  double tns = chit->TNS();
481  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
482 
483  double xyz[3];
484  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
485  if (chit->Plane() == P_x0 ) {
486  x0s_x.emplace_back(xyz[0]);
487  x0s_z.emplace_back(xyz[2]);
488  }
489  if (chit->Plane() == P_x1 ) {
490  x1s_x.emplace_back(xyz[0]);
491  x1s_z.emplace_back(xyz[2]);
492  }
493 
494  if (chit->Plane() == P_y0 ) {
495  y0s_y.emplace_back(xyz[1]);
496  y0s_z.emplace_back(xyz[2]);
497  }
498  if (chit->Plane() == P_y1 ) {
499  y1s_y.emplace_back(xyz[1]);
500  y1s_z.emplace_back(xyz[2]);
501  }
502  }
503  double z_x0=0;
504  double x_x0=0;
505  double z_y0=0;
506  double y_y0=0;
507  double z_x1=0;
508  double z_y1=0;
509  double dx = 0;
510  double dy = 0;
511  double dxdz = 0;
512  double dydz = 0;
513 
514  for (unsigned h0 = 0; h0 != x0s_x.size(); ++h0) {
515  for (unsigned h1 = 0; h1 != x1s_x.size(); ++h1) {
516  if (fabs(dx) < fabs(x1s_x[h1]-x0s_x[h0])) {
517  dx = x1s_x[h1]-x0s_x[h0];
518  dxdz = dx/(x1s_z[h1]-x0s_z[h0]);
519  z_x0 = x0s_z[h0];
520  x_x0 = x0s_x[h0];
521  z_x1 = x1s_z[h1];
522  }
523  }
524  }
525  for (unsigned h0 = 0; h0 != y0s_y.size(); ++h0) {
526  for (unsigned h1 = 0; h1 != y1s_y.size(); ++h1) {
527  if (fabs(dy) < fabs(y1s_y[h1]-y0s_y[h0])) {
528  dy = y1s_y[h1]-y0s_y[h0];
529  dydz = dy/(y1s_z[h1]-y0s_z[h0]);
530  z_y0 = y0s_z[h0];
531  y_y0 = y0s_y[h0];
532  z_y1 = y1s_z[h1];
533  }
534  }
535  }
536  //Now you can determine the two ends:
537  double zi = std::min(z_x0,z_y0);
538  double zf = std::max(z_x1,z_y1);
539  double xi = dxdz*(zi-z_x0)+x_x0;
540  double xf = dxdz*(zf-z_x0)+x_x0;
541  double yi = dydz*(zi-z_y0)+y_y0;
542  double yf = dydz*(zf-z_y0)+y_y0;
543  TVector3 Pi(xi,yi,zi); // Low Z end
544  TVector3 Pf(xf,yf,zf); // High Z end
545 
546  //Now, we know the direction of the track, can provide the W position,
547  //Loop through the hits attached to the track again!
548  MIP_tot = 0;
549  GeV_tot = 0;
550  double t_min = 9999999;
551  double t_max = -9999999;
552 
553  std::vector<double> ts_sat,terrs_sat;
554  std::vector<double> ts_unsat,terrs_unsat;
555  std::vector<TVector3> pos_sat, pos_unsat;
556 
557  lxsqr = 0;
558  lysqr = 0;
559 
560  //In the following loop, determine the ADC center from the geometry center:
561  TVector3 adc_Center(0,0,0);
562  TVector3 GeV_Center(0,0,0);
563  TVector3 geo_Center(0,0,0);
564  ADC_tot = 0;
565 
566  for (unsigned i = 0; i != track.NCell(); ++i) {
567  const art::Ptr<rb::CellHit> chit = track.Cell(i);
568  double tns = chit->TNS();
569  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
570  //Interpolate the W position:
571  double W;
572  double xyz[3];
573  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
574  if (chit->View()==geo::kX) {
575  W = dydz*(xyz[2]-z_y0)+y_y0;
576  lxsqr += doca2(xyz[2],xyz[0],z_x0,x_x0, dxdz);
577  }
578  else {
579  W = dxdz*(xyz[2]-z_x0)+x_x0;
580  lysqr += doca2(xyz[2],xyz[1],z_y0,y_y0, dydz);
581  }
582 
583  rb::RecoHit rhit(calib->MakeRecoHit(*chit,W));
584  if (!rhit.IsCalibrated()) continue;
585  float hitE = rhit.GeV();
586  //float hitZ = rhit.Z();
587  float hitT = rhit.T();
588  unsigned int hitCell = chit->Cell();
589  unsigned int hitPlane = chit->Plane();
590  if(chit->View() == geo :: kX){
591  float hitX = rhit.X();
592  hitxCell.push_back(hitCell);
593  hitxPlane.push_back(hitPlane);
594  hitxe.push_back(hitE);
595  hitxt.push_back(hitT);
596  hitxx.push_back(hitX);
597  }
598  else if(chit->View() == geo :: kY){
599  float hitY = rhit.Y();
600  hityCell.push_back(hitCell);
601  hityPlane.push_back(hitPlane);
602  hitye.push_back(hitE);
603  hityt.push_back(hitT);
604  hityy.push_back(hitY);
605  }//end y view
606 
607  MIP_tot += rhit.MIP();
608  GeV_tot += rhit.GeV();
609 
610  TVector3 pos(rhit.X(), rhit.Y(), rhit.Z());
611  double st = calib->GetTimeRes(*chit);
612 
613  adc_Center += (chit->ADC())*pos;
614  GeV_Center += rhit.GeV() * pos;
615  geo_Center += pos;
616  ADC_tot += chit->ADC();
617 
618  if(chit->IsMC()==true) ++nhits_true;
619 
620  if (chit->ADC(3) == 4095) {
621  pos_sat.emplace_back(pos-Pi);
622  ts_sat.emplace_back(rhit.T()-tns_min);
623  terrs_sat.emplace_back(st);
624  }
625 
626  else {
627  pos_unsat.emplace_back(pos-Pi);
628  ts_unsat.emplace_back(rhit.T()-tns_min);
629  terrs_unsat.emplace_back(st);
630  }
631 
632  if (rhit.T() < t_min ) {
633  x_tmin = rhit.X();
634  y_tmin = rhit.Y();
635  z_tmin = rhit.Z();
636  t_min = rhit.T();
637  }
638 
639  if (rhit.T() > t_max ) {
640  x_tmax = rhit.X();
641  y_tmax = rhit.Y();
642  z_tmax = rhit.Z();
643  t_max = rhit.T();
644  }
645  }
646 
647  double ADC_mean=ADC_tot/(nxhits+nyhits+0.0);
648  double ADC_sq_sum = 0;
649  for (unsigned i = 0; i != track.NCell(); ++i) {
650  const art::Ptr<rb::CellHit> chit = track.Cell(i);
651  double tns = chit->TNS();
652  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
653  ADC_sq_sum+=chit->ADC()*chit->ADC();
654  }
655  ADC_stdev = std::sqrt(ADC_sq_sum / (nxhits+nyhits+0.0) - ADC_mean*ADC_mean);
656 
657  adc_Center = adc_Center*(1/double(ADC_tot));
658  GeV_Center = GeV_Center*(1/double(GeV_tot));
659  geo_Center = geo_Center*(1/double(ts_sat.size()+ts_unsat.size()));
660 
661  adc_W = (adc_Center-geo_Center).Mag();
662  GeV_W = (GeV_Center - geo_Center).Mag();
663 
664  lxsqr /= double(nxhits);
665  lysqr /= double(nyhits);
666 
667  nhits_sat = pos_sat.size();
668 
669  vrc_sat = vfit(pos_sat, ts_sat, terrs_sat, vSat_err, tSat_mean, tSat_rms, Chisqr_TS);
670  vrc_nonsat = vfit(pos_unsat, ts_unsat, terrs_unsat, vnonSat_err, tnonSat_mean, tnonSat_rms, Chisqr_TN);
671 
672  if (nhits_sat<2) {
673  //Make them apparently invalid
674  vrc_sat = 60;
675  vSat_err = -1;
676  tSat_rms = -1;
677  tSat_mean = -1e6;
678  }
679  if (pos_unsat.size()<2) {
680  //Make them apparently invalid
681  vrc_nonsat = 60;
682  vnonSat_err = -1;
683  tnonSat_rms = -1;
684  tnonSat_mean = -1e6;
685  }
686 
687  double fitpar[3];
688  LinFit(hitxx, hitxt, fitpar);
689  r2xt = 2;
690  if( std::isfinite(fitpar[2]))
691  r2xt = fitpar[2];
692  LinFit(hityy, hityt, fitpar);
693  r2yt = 2;
694  if( std::isfinite(fitpar[2]))
695  r2yt = fitpar[2];
696  LinFit(hitxe, hitxx, fitpar);
697  r2GeVx = 2;
698  if( std::isfinite(fitpar[2]))
699  r2GeVx = fitpar[2];
700  LinFit(hitye, hityy, fitpar);
701  r2GeVy = 2;
702  if( std::isfinite(fitpar[2]))
703  r2GeVy = fitpar[2];
704 
705  //Just set the lowZ end to be the starting point:
706  Xi = xi;
707  Yi = yi;
708  Zi = zi;
709  Ti = t_min;
710  Xf = xf;
711  Yf = yf;
712  Zf = zf;
713  Tf = t_max;
714 
715  std::sort (hitxt.begin(), hitxt.end());
716  double max_time_gap = 0;
717  double time_cross = hitxt.back()-hitxt.front()+0.00000001;
718  for (unsigned idx =0; idx+1<hitxt.size(); ++idx) {
719  double dt = hitxt[idx+1] - hitxt[idx];
720  if (dt > max_time_gap) max_time_gap = dt;
721  }
722  time_gap_xz = max_time_gap / time_cross;
723  std::sort (hityt.begin(), hityt.end());
724  max_time_gap = 0;
725  time_cross = hityt.back()-hityt.front()+0.00000001;
726  for (unsigned idx =0; idx+1<hityt.size(); ++idx) {
727  double dt = hityt[idx+1] - hityt[idx];
728  if (dt > max_time_gap) max_time_gap = dt;
729  }
730  time_gap_yz = max_time_gap / time_cross;
731 
732  RCV_Tree->Fill();
733  filled = true;
734  }
735  }
736  }
737  if (!filled) RCV_Tree->Fill();
738 
739 } //end analyzing
740 
741 double zcl::FmmTrackerValidation::dist2(double dx, double dy) {
742  return dx*dx+dy*dy;
743 }
744 void zcl::FmmTrackerValidation::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
745  // http://en.wikipedia.org/wiki/Simple_linear_regression
746  const auto n = x_val.size();
747  const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n; // <x>
748  const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n; // <y>
749  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n;// <xx>
750  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n;// <yy>
751  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n;// <xy>
752 
753  const auto b = (xy - x*y)/(xx - x*x); // slope
754  const auto a = y - b*x; // intercept
755  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
756  fitpar[0] = a;
757  fitpar[1] = b;
758  fitpar[2] = r*r;
759 }
760 
761 double zcl::FmmTrackerValidation::doca2(double x, double y, double x0, double y0, double dydx)
762 {
763  if (fabs(x-x0) < 0.1 && fabs(y-y0)<0.1) return 0;
764  double cos2 = (x-x0 + dydx*(y-y0))*(x-x0 + dydx*(y-y0))/dist2(x-x0,y-y0)/dist2(1,dydx);
765  double docasqr = dist2(x-x0,y-y0)*(1-cos2);
766  // std::cout<<"docasqr: "<<docasqr<<" cos2:"<<cos2<<std::endl;
767  return docasqr;
768 }
769 
770 double zcl::FmmTrackerValidation::vfit(std::vector<TVector3> & Positions,
771  std::vector<double> & Ts,
772  std::vector<double> & Terr,
773  double & v_err,
774  double & Tmean,
775  double & Trms,
776  double & dTsqr)
777 {
778  int m = Ts.size();
779  double sumL = 0;
780  double sumT = 0;
781  double sumLT = 0;
782  double sumLL = 0;
783 
784  std::vector<double> Vs;
785 
786  for (int i = 0; i!=m; ++i ) {
787  sumL += Positions[i].Mag();
788  sumT += Ts[i];
789  sumLT += Positions[i].Mag() * Ts[i];
790  sumLL += Positions[i].Mag() * Positions[i].Mag();
791  }
792 
793  Tmean = sumT/double(m);
794  double v = (sumL*sumL-m*sumLL)/(sumT*sumL-m*sumLT);
795  if (v > 30) v = 30;
796  if (v <-30) v = -30;
797 
798  Vs.emplace_back(v);
799 
800  for (int t = 0; t!=ntrials; ++t) {
801  sumL = 0;
802  sumT = 0;
803  sumLT = 0;
804  sumLL = 0;
805  for (int i = 0; i!=m; ++i ) {
806  double xi = genR->Uniform(Positions[i].X()-X_err,Positions[i].X()+X_err);
807  double yi = genR->Uniform(Positions[i].Y()-Y_err,Positions[i].Y()+Y_err);
808  double zi = genR->Uniform(Positions[i].Z()-Z_err,Positions[i].Z()+Z_err);
809  double ti = genR->Gaus(Ts[i], Terr[i]);
810  double li = sqrt(xi*xi+yi*yi+zi*zi);
811 
812  sumL += li;
813  sumT += ti;
814  sumLT += li*ti;
815  sumLL += li*li;
816  }
817  v = (sumL*sumL-m*sumLL)/(sumT*sumL-m*sumLT);
818  if (v > 30) v = 30;
819  if (v <-30) v = -30;
820  Vs.emplace_back(v);
821  }
822 
823  //Calculate the mean of v, and standard error:
824  double v_mean = 0;
825  v_err = 0;
826  for (unsigned j = 0; j!=Vs.size(); ++j)
827  v_mean += Vs[j];
828  v_mean /= double(Vs.size());
829 
830  for (unsigned j = 0; j!=Vs.size(); ++j)
831  v_err += (Vs[j]-v_mean)*(Vs[j]-v_mean);
832 
833  v_err = sqrt(v_err)/double(Vs.size());
834 
835  if (v_mean < 0.003 && v_mean > 0)
836  v_mean = 0.003;
837 
838  else if (v_mean > -0.003 && v_mean < 0)
839  v_mean = -0.003;
840 
841  double t0 = Tmean-sumL/(m+0.)/v_mean;
842 
843  Trms = 0;
844  dTsqr = 0;
845  for (int i = 0; i!=m; ++i ) {
846  Trms += (Ts[i]-Tmean)*(Ts[i]-Tmean);
847  dTsqr+= (Ts[i]-t0-Positions[i].Mag()/v_mean)* (Ts[i]-t0-Positions[i].Mag()/v_mean);
848  }
849  Trms = sqrt(Trms/(m+0.));
850  dTsqr = sqrt(dTsqr/(m+0.));
851 
852  return v_mean;
853 }
854 
855 void zcl::FmmTrackerValidation::endJob(){
856 }
857 DEFINE_ART_MODULE(zcl::FmmTrackerValidation)
rb::CellHit::TNS
float TNS() const
Definition: CellHit.h:46
std::max
T max(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:416
BackTracker.h
back track the reconstruction to the simulation
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:211
calib::Calibrator::GetTimeRes
double GetTimeRes(rb::CellHit const &cellhit)
Definition: Calibrator_service.cc:716
zcl::FmmTrackerValidation::analyze
void analyze(const art::Event &evt)
Definition: FmmTrackerValidation_module.cc:222
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
xx
Double_t xx
Definition: macro.C:12
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:230
boost::math::isfinite
bool isfinite(const stan::math::var &v)
Definition: boost_isfinite.hpp:21
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
cheat::BackTracker::ParticleNavigator
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
geo::CellGeo::LocalToWorld
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
zcl::FmmTrackerValidation::FmmTrackerValidation
FmmTrackerValidation(fhicl::ParameterSet const &pset)
Definition: FmmTrackerValidation_module.cc:107
rb::CellHit::View
geo::View_t View() const
Definition: CellHit.h:41
geo::PlaneGeo::Cell
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
p
const char * p
Definition: xmltok.h:285
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cheat::BackTracker::ParticleToFLSHit
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Definition: BackTracker_service.cc:1514
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:238
sim::ParticleNavigator::const_iterator
list_type::const_iterator const_iterator
Definition: ParticleNavigator.h:86
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:229
calib::Calibrator::MakeRecoHit
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
Definition: Calibrator_service.cc:575
zcl::FmmTrackerValidation::LinFit
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
Definition: FmmTrackerValidation_module.cc:744
track
Definition: event.h:19
site_stats_from_log.get
get
Definition: site_stats_from_log.py:192
make_training.momentum
float momentum
Definition: make_training.py:78
rb::Track
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
DEFINE_ART_MODULE
DEFINE_ART_MODULE(TestTMapFile)
tf
Timing fit.
Definition: TimingFit_module.cc:23
compare_h5_caf.idx
int idx
Definition: compare_h5_caf.py:69
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
submit_syst.y
y
Definition: submit_syst.py:249
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:209
Y
Float_t Y
Definition: plot.C:38
rb::RecoHit
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
m
#define m
Definition: MakeGENIEReweightTable_module.cc:420
Z
Float_t Z
Definition: plot.C:38
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
dt
Definition: Cand.cxx:23
calib
CDPStorage service.
Definition: AbsCalib_module.cc:25
submit_syst.x
x
Definition: submit_syst.py:249
dy
double dy[NP][NC]
Definition: checkMeasured_Target.C:19
dx
double dx[NP][NC]
Definition: checkMeasured_Target.C:18
local
#define local
Definition: gzguts.h:107
a
const double a
Definition: gUpMuFluxGen.cxx:164
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:233
evt
int evt
Definition: FillTrainingTrees.h:22
rawdata::RawDigit::IsMC
bool IsMC() const
Definition: RawDigit.h:108
Geo.h
Collect Geo headers and supply basic geometry functions.
zcl::FmmTrackerValidation::vfit
double vfit(std::vector< TVector3 > &Positions, std::vector< double > &Ts, std::vector< double > &Terr, double &v_err, double &Tmean, double &Trms, double &dTsqr)
Definition: FmmTrackerValidation_module.cc:770
calib::j
const double j
Definition: BetheBloch.cxx:29
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Event.h
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
h1
TH1F * h1
Definition: nslices_potspill.C:4
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:220
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
rb::RecoHit::GeV
float GeV() const
Definition: RecoHit.cxx:69
zcl::FmmTrackerValidation::dist2
double dist2(double dx, double dy)
Definition: FmmTrackerValidation_module.cc:741
rawdata::RawDigit::ADC
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
geom
void geom(int which=0)
Definition: geom.C:163
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:231
Mag
float Mag() const
Definition: SRLorentzVectorAddons.h:10
zcl::FmmTrackerValidation::doca2
double doca2(double x, double y, double x0, double y0, double dydx)
Definition: FmmTrackerValidation_module.cc:761
b
const hit & b
Definition: hits.cxx:21
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:222
r
TRandom3 r(0)
registry_explorer.v
v
Definition: registry_explorer.py:22
std::min
T min(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:406
zcl::FastMonopoleTriggers::is_trigger_by_epoch1_fmmtrigger
bool is_trigger_by_epoch1_fmmtrigger(art::Handle< std::vector< rb::Cluster > > slices, int &_penetrated, double &_nSatHits, double &_meanADC)
Definition: FastMonopoleTriggers.cxx:3
geo::CellGeo
Encapsulate the cell geometry.
Definition: CellGeo.h:25
X
Float_t X
Definition: plot.C:38
W
#define W(x)
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:221
Geometry.h
Encapsulate the geometry of one entire detector (near, far, ndos)
string
enum BeamMode string
Generated on Mon Mar 1 2021 14:14:25 for NOvA by   doxygen 1.8.11