FmmTrackerAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FmmTrackerAna
3 // Module Type: analyzer
4 // File: FmmTrackerAna_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 ////////////////////////////////////////////////////////////////////////
14 
15 #include "MCCheater/BackTracker.h"
16 #include "Calibrator/Calibrator.h"
17 #include "Simulation/Particle.h"
21 #include "Utilities/func/MathUtil.h"
22 #include "fhiclcpp/ParameterSet.h"
24 
25 #include "Geometry/Geometry.h"
27 #include "GeometryObjects/Geo.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 
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 FmmTrackerAna;
50 }
51 
53 public:
54  explicit FmmTrackerAna(fhicl::ParameterSet const & pset);
55  virtual ~FmmTrackerAna();
56 
57  void analyze(const art::Event & evt);
58  void beginJob();
59  void endJob();
60 
61 private:
63 
64  const int ntrials = 100; //Number of MC trilas in determining the fitting error
65  const double X_err = 1.8;
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 * RC_Tree;
71  TRandom *genR = new TRandom3();
73  int nxhits, nyhits;
74  int nhits_sat;
75  double MIP_tot;
76  double ADC_tot;
77  double GeV_tot;
78 
80 
82 
83  //Notice the start and end point in the fmmtrack might be flipped:
84  double Xi, Yi, Zi, Ti;
85  double Xf, Yf, Zf, Tf;
86 
87  double t_min, t_max;
89 
90  double lxsqr, lysqr;
91 
92  double vrc_sat; // cm/ns
93  double tSat_mean;
94  double tSat_rms;
95  double vSat_err;
96  double vrc_nonsat; // cm/ns
97  double tnonSat_mean;
98  double tnonSat_rms;
99  double vnonSat_err;
100  double Chisqr_TS;
101  double Chisqr_TN;
102 
103  double adc_W;
104 
105  double dist2(double dx, double dy);
106  double doca2(double x, double y, double x0, double y0, double dydx);
107 
108  //Linear Regression function
109  double vfit(std::vector<TVector3> & Positions,
110  std::vector<double> & Ts,
111  std::vector<double> & Terr,
112  double & v_err,
113  double & Tmean,
114  double & Trms,
115  double & dTsqr);
117 
118 };
119 
121  EDAnalyzer(p),
122  fTrackerInput (p.get< std::string > ("TrackerInput") )
123 {
124 }
125 
127 {
128  // Clean up dynamic memory and other resources here.
129 }
130 
132 {
134 
135  //Combining cheater info and the RC info to make the data structure better
136  //One RC track info per entry, if no tracks reconstructed in this event, fill in 1 entry anyway tagging the
137  //reconstruction failed
138  RC_Tree = tfs->make<TTree>("RC_Tree","Information from RC track info");
139  RC_Tree -> Branch("event", &evtID, "event/I");
140  RC_Tree -> Branch("subRun", &subRunID, "subRun/I");
141  RC_Tree -> Branch("Run", &RunID, "Run/I");
142  RC_Tree -> Branch("nxhits", &nxhits, "nxhits/I");
143  RC_Tree -> Branch("nyhits", &nyhits, "nyhits/I");
144  RC_Tree -> Branch("nhits_sat", &nhits_sat, "nhits_sat/I");
145  RC_Tree -> Branch("MIP_tot", &MIP_tot, "MIP_tot/D");
146  RC_Tree -> Branch("ADC_tot", &ADC_tot, "ADC_tot/D");
147  RC_Tree -> Branch("GeV_tot", &GeV_tot, "GeV_tot/D");
148 
149  //Missing Planes and the max gap:
150  RC_Tree -> Branch("max_gap", &max_gap, "max_gap/I");
151  RC_Tree -> Branch("miss_plane", &miss_plane, "miss_plane/I");
152  RC_Tree -> Branch("p_cross", &p_cross, "p_cross/I");
153  RC_Tree -> Branch("p_overlap", &p_overlap, "p_overlap/I");
154 
155  RC_Tree -> Branch("Xi", &Xi, "Xi/D");
156  RC_Tree -> Branch("Xf", &Xf, "Xf/D");
157  RC_Tree -> Branch("Yi", &Yi, "Yi/D");
158  RC_Tree -> Branch("Yf", &Yf, "Yf/D");
159  RC_Tree -> Branch("Zi", &Zi, "Zi/D");
160  RC_Tree -> Branch("Zf", &Zf, "Zf/D");
161  RC_Tree -> Branch("Ti", &Ti, "Ti/D");
162  RC_Tree -> Branch("Tf", &Tf, "Tf/D");
163 
164  RC_Tree -> Branch("tns_mean", &tns_mean, "tns_mean/D");
165  RC_Tree -> Branch("tns_rms", &tns_rms, "tns_rms/D");
166  RC_Tree -> Branch("tns_min", &tns_min, "tns_min/D");
167  RC_Tree -> Branch("tns_max", &tns_max, "tns_max/D");
168 
169  //This is double checking the velocity check:
170  RC_Tree -> Branch("t_min", &t_min, "t_min/D");
171  RC_Tree -> Branch("t_max", &t_max, "t_max/D");
172  RC_Tree -> Branch("x_tmin", &x_tmin, "x_tmin/D");
173  RC_Tree -> Branch("x_tmax", &x_tmax, "x_tmax/D");
174  RC_Tree -> Branch("y_tmin", &y_tmin, "y_tmin/D");
175  RC_Tree -> Branch("y_tmax", &y_tmax, "y_tmax/D");
176  RC_Tree -> Branch("z_tmin", &z_tmin, "z_tmin/D");
177  RC_Tree -> Branch("z_tmax", &z_tmax, "z_tmax/D");
178 
179  //fitting parameters:
180  RC_Tree -> Branch("lxsqr", &lxsqr, "lxsqr/D");
181  RC_Tree -> Branch("lysqr", &lysqr, "lysqr/D");
182 
183  RC_Tree -> Branch("vrc_sat", &vrc_sat, "vrc_sat/D");
184  RC_Tree -> Branch("tSat_mean", &tSat_mean,"tSat_mean/D");
185  RC_Tree -> Branch("tSat_rms", &tSat_rms, "tSat_rms/D");
186  RC_Tree -> Branch("vSat_err", &vSat_err, "vsat_err/D");
187  RC_Tree -> Branch("vrc_nonsat", &vrc_nonsat, "vrc_nonsat/D");
188  RC_Tree -> Branch("tnonSat_mean", &tnonSat_mean,"tnonSat_mean/D");
189  RC_Tree -> Branch("tnonSat_rms", &tnonSat_rms, "tnonSat_rms/D");
190  RC_Tree -> Branch("vnonSat_err", &vnonSat_err, "vnonSat_err/D");
191  RC_Tree -> Branch("Chisqr_TS", &Chisqr_TS, "Chisqr_TS/D");
192  RC_Tree -> Branch("Chisqr_TN", &Chisqr_TN, "Chisqr_TN/D");
193 
194  RC_Tree -> Branch("adc_W", &adc_W, "adc_W/D");
195 }
196 
198 {
199  evtID = evt.id().event();
200  subRunID = evt.subRun();
201  RunID = evt.run();
202 
205 
206  //=====================================================================================
208  evt.getByLabel(fTrackerInput,trks);
209 
210  if (trks->size() < 1) {
211  std::cout<<"There is no track reconstructed in this event!"<<std::endl;
212  return;
213  }
214 
215  //take the entire 3D track:
216  for (unsigned trkIdx = 0; trkIdx != trks->size(); ++trkIdx ) {
217  const rb::Track& track = (*trks)[trkIdx];
218 
219  if (track.NCell()<6) continue;
220 
221  tns_mean = 0; //ADC Weighed average time
222  double ADC_sum = 0;
223 
224  //With the first round of loop, can determine the following info: t_mean
225 
226  for (unsigned i = 0; i != track.NCell(); ++i) {
227  const art::Ptr<rb::CellHit> chit = track.Cell(i);
228  double tns = chit->TNS();
229  double adc = chit->ADC();
230  tns_mean += tns*adc;
231  ADC_sum += adc;
232  }
233 
234  tns_mean /= ADC_sum;
235  //=============================================================
236 
237  //Just determine the RMS of the time:
238  tns_rms = 0;
239  for (unsigned i = 0; i != track.NCell(); ++i) {
240  const art::Ptr<rb::CellHit> chit = track.Cell(i);
241  double tns = chit->TNS();
242  double adc = chit->ADC();
243  tns_rms = (tns-tns_mean)*adc*(tns-tns_mean)*adc;
244  }
245  tns_rms = sqrt(tns_rms/ADC_sum+250000/12.);
246 
247  //==============================================================
248  //Now, loop through the hits again to throw away the hits outside the 4 sigma range,
249  //Make sure they get updated!
250  tns_min = 99999999;
251  tns_max = -99999999;
252  nxhits = 0;
253  nyhits = 0;
254  nhits_sat = 0;
255  ADC_tot = 0;
256  std::vector<int> planes;
257 
258  int P_x0=9999;
259  int P_x1=0;
260  int P_y0=9999;
261  int P_y1=0;
262 
263  for (unsigned i = 0; i != track.NCell(); ++i) {
264  const art::Ptr<rb::CellHit> chit = track.Cell(i);
265  double tns = chit->TNS();
266  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
267  if (tns > tns_max) tns_max = tns;
268  if (tns < tns_min) tns_min = tns;
269  planes.emplace_back(chit->Plane());
270  ADC_tot += chit->ADC();
271 
272  if (chit->View()==geo::kX) {
273  ++nxhits;
274  if (chit->Plane() < P_x0) P_x0 = chit->Plane();
275  if (chit->Plane() > P_x1) P_x1 = chit->Plane();
276  }
277  else {
278  ++nyhits;
279  if (chit->Plane() < P_y0) P_y0 = chit->Plane();
280  if (chit->Plane() > P_y1) P_y1 = chit->Plane();
281  }
282  }
283 
284  if (planes.size()<3) continue;
285 
286  if (P_x1-P_x0<3) continue;
287  if (P_y1-P_y0<3) continue;
288 
289  //Now calculate the missing planes and the max gap:
290  std::sort (planes.begin(), planes.end());
291  max_gap = 0;
292  miss_plane = 0;
293  p_cross = planes.back()-planes.front()+1;
294  p_overlap = 1;
295  bool start_flag = false;
296  for (unsigned pidx =0; pidx+1<planes.size(); ++pidx) {
297  int dP = planes[pidx+1] - planes[pidx] - 1;
298  //Calculate the number of planes in the overlapped region:
299  if (!start_flag) {
300  if ( (planes[pidx+1] - planes[pidx]) % 2 ) start_flag = true;
301  }
302  else {
303  if ( (planes[pidx+1] - planes[pidx]) % 2 ) ++p_overlap;
304  }
305  if (dP > 0 ) {
306  miss_plane += dP;
307  if (dP > max_gap) max_gap = dP;
308  }
309  }
310  //========================================================
311  //Loop through the hits again: this time, determine dxdz and dydz
312 
313  std::vector<double> x0s_x, x1s_x, y0s_y, y1s_y;
314  std::vector<double> x0s_z, x1s_z, y0s_z, y1s_z;
315  for (unsigned i = 0; i != track.NCell(); ++i) {
316  const art::Ptr<rb::CellHit> chit = track.Cell(i);
317  double tns = chit->TNS();
318  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
319 
320  double xyz[3];
321  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
322 
323  if (chit->Plane() == P_x0 ) {
324  x0s_x.emplace_back(xyz[0]);
325  x0s_z.emplace_back(xyz[2]);
326  }
327  if (chit->Plane() == P_x1 ) {
328  x1s_x.emplace_back(xyz[0]);
329  x1s_z.emplace_back(xyz[2]);
330  }
331 
332  if (chit->Plane() == P_y0 ) {
333  y0s_y.emplace_back(xyz[1]);
334  y0s_z.emplace_back(xyz[2]);
335  }
336  if (chit->Plane() == P_y1 ) {
337  y1s_y.emplace_back(xyz[1]);
338  y1s_z.emplace_back(xyz[2]);
339  }
340  }
341 
342  double z_x0=0;
343  double x_x0=0;
344  double z_y0=0;
345  double y_y0=0;
346  double z_x1=0;
347  double z_y1=0;
348 
349  double dx = 0;
350  double dy = 0;
351  double dxdz = 0;
352  double dydz = 0;
353 
354  for (unsigned h0 = 0; h0 != x0s_x.size(); ++h0) {
355  for (unsigned h1 = 0; h1 != x1s_x.size(); ++h1) {
356  if (fabs(dx) < fabs(x1s_x[h1]-x0s_x[h0])) {
357  dx = x1s_x[h1]-x0s_x[h0];
358  dxdz = dx/(x1s_z[h1]-x0s_z[h0]);
359  z_x0 = x0s_z[h0];
360  x_x0 = x0s_x[h0];
361  z_x1 = x1s_z[h1];
362  }
363  }
364  }
365 
366  for (unsigned h0 = 0; h0 != y0s_y.size(); ++h0) {
367  for (unsigned h1 = 0; h1 != y1s_y.size(); ++h1) {
368  if (fabs(dy) < fabs(y1s_y[h1]-y0s_y[h0])) {
369  dy = y1s_y[h1]-y0s_y[h0];
370  dydz = dy/(y1s_z[h1]-y0s_z[h0]);
371 
372  z_y0 = y0s_z[h0];
373  y_y0 = y0s_y[h0];
374  z_y1 = y1s_z[h1];
375  }
376  }
377  }
378 
379  //Now you can determine the two ends:
380  double zi = std::min(z_x0,z_y0);
381  double zf = std::max(z_x1,z_y1);
382 
383  double xi = dxdz*(zi-z_x0)+x_x0;
384  double xf = dxdz*(zf-z_x0)+x_x0;
385  double yi = dydz*(zi-z_y0)+y_y0;
386  double yf = dydz*(zf-z_y0)+y_y0;
387 
388  TVector3 Pi(xi,yi,zi); // Low Z end
389  TVector3 Pf(xf,yf,zf); // High Z end
390 
391  //Now, we know the direction of the track, can provide the W position,
392  //Loop through the hits attached to the track again!
393  MIP_tot = 0;
394  GeV_tot = 0;
395  double t_min = 9999999;
396  double t_max = -9999999;
397 
398  std::vector<double> ts_sat,terrs_sat;
399  std::vector<double> ts_unsat,terrs_unsat;
400  std::vector<TVector3> pos_sat, pos_unsat;
401 
402  lxsqr = 0;
403  lysqr = 0;
404 
405  //In the following loop, determine the ADC center from the geometry center:
406  TVector3 adc_Center(0,0,0);
407  TVector3 geo_Center(0,0,0);
408  ADC_tot = 0;
409 
410  for (unsigned i = 0; i != track.NCell(); ++i) {
411  const art::Ptr<rb::CellHit> chit = track.Cell(i);
412  double tns = chit->TNS();
413  if (tns > tns_mean + 4*tns_rms || tns < tns_mean - 4*tns_rms) continue;
414 
415  //Interpolate the W position:
416  double W;
417  double xyz[3];
418  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
419 
420  if (chit->View()==geo::kX) {
421  W = dydz*(xyz[2]-z_y0)+y_y0;
422  lxsqr += doca2(xyz[2],xyz[0],z_x0,x_x0, dxdz);
423  // std::cout<<"lxsqr: "<<lxsqr<<std::endl;
424  }
425  else {
426  W = dxdz*(xyz[2]-z_x0)+x_x0;
427  lysqr += doca2(xyz[2],xyz[1],z_y0,y_y0, dydz);
428  // std::cout<<"lysqr: "<<lysqr<<std::endl;
429  }
430 
431  rb::RecoHit rhit(calib->MakeRecoHit(*chit,W));
432 
433  if (!rhit.IsCalibrated()) continue;
434 
435  MIP_tot += rhit.MIP();
436  GeV_tot += rhit.GeV();
437 
438  TVector3 pos(rhit.X(), rhit.Y(), rhit.Z());
439 
440  adc_Center += (chit->ADC())*pos;
441  geo_Center += pos;
442  ADC_tot += chit->ADC();
443 
444  double st = calib->GetTimeRes(*chit);
445 
446  if (chit->ADC(3) == 4095) {
447  pos_sat.emplace_back(pos-Pi);
448  ts_sat.emplace_back(rhit.T()-tns_min);
449  terrs_sat.emplace_back(st);
450  }
451 
452  else {
453  pos_unsat.emplace_back(pos-Pi);
454  ts_unsat.emplace_back(rhit.T()-tns_min);
455  terrs_unsat.emplace_back(st);
456  }
457 
458  if (rhit.T() < t_min ) {
459  x_tmin = rhit.X();
460  y_tmin = rhit.Y();
461  z_tmin = rhit.Z();
462  t_min = rhit.T();
463  }
464 
465  if (rhit.T() > t_max ) {
466  x_tmax = rhit.X();
467  y_tmax = rhit.Y();
468  z_tmax = rhit.Z();
469  t_max = rhit.T();
470  }
471  }
472 
473  adc_Center = adc_Center*(1/double(ADC_tot));
474  geo_Center = geo_Center*(1/double(ts_sat.size()+ts_unsat.size()));
475 
476  adc_W = (adc_Center-geo_Center).Mag();
477 
478  lxsqr /= double(nxhits);
479  lysqr /= double(nyhits);
480  nhits_sat = pos_sat.size();
481 
482 
483  // std::cout<<"final lxsqr: "<<lxsqr<<" final lysqr: "<<lysqr<<std::endl;
484 
485  if (nhits_sat > 1) {
486  vrc_sat = vfit(pos_sat, ts_sat, terrs_sat, vSat_err, tSat_mean, tSat_rms,Chisqr_TS);
487  }
488 
489  else {
490  //Make them apparently invalid
491  vrc_sat = 60;
492  vSat_err = -1;
493  tSat_rms = -1;
494  tSat_mean = -1e6;
495  }
496 
497  if (pos_unsat.size() > 1) {
498  vrc_nonsat = vfit(pos_unsat, ts_unsat, terrs_unsat, vnonSat_err, tnonSat_mean, tnonSat_rms,Chisqr_TN);
499  }
500 
501  else {
502  //Make them apparently invalid
503  vrc_nonsat = 60;
504  vnonSat_err = -1;
505  tnonSat_rms = -1;
506  }
507 
508  //Just set the lowZ end to be the starting point:
509  Xi = xi;
510  Yi = yi;
511  Zi = zi;
512  Ti = t_min;
513  Xf = xf;
514  Yf = yf;
515  Zf = zf;
516  Tf = t_max;
517 
518  RC_Tree->Fill();
519  }
520 
521 } //end analyzing
522 
523 double zcl::FmmTrackerAna::dist2(double dx, double dy) {
524  return dx*dx+dy*dy;
525 }
526 
527 double zcl::FmmTrackerAna::doca2(double x, double y, double x0, double y0, double dydx)
528 {
529  if (fabs(x-x0) < 0.1 && fabs(y-y0)<0.1) return 0;
530  double cos2 = (x-x0 + dydx*(y-y0))*(x-x0 + dydx*(y-y0))/dist2(x-x0,y-y0)/dist2(1,dydx);
531  if (cos2==1) return (y-y0)*(y-y0);
532  double docasqr = dist2(x-x0,y-y0)*(1-cos2);
533  /*
534  std::cout<<"x: "<<x<<" y: "<<y<<" x0: "<<x0<<" y0: "<<y0
535  <<" dydx: "<<dydx<<" cos2: "<<cos2<<" docasqr: "<<docasqr<<std::endl;
536  */
537  return docasqr;
538 }
539 
540 double zcl::FmmTrackerAna::vfit(std::vector<TVector3> & Positions,
541  std::vector<double> & Ts,
542  std::vector<double> & Terr,
543  double & v_err,
544  double & Tmean,
545  double & Trms,
546  double & dTsqr)
547 {
548  int m = Ts.size();
549  double sumL = 0;
550  double sumT = 0;
551  double sumLT = 0;
552  double sumLL = 0;
553 
554  std::vector<double> Vs;
555 
556  for (int i = 0; i!=m; ++i ) {
557  sumL += Positions[i].Mag();
558  sumT += Ts[i];
559  sumLT += Positions[i].Mag() * Ts[i];
560  sumLL += Positions[i].Mag() * Positions[i].Mag();
561  }
562 
563  Tmean = sumT/double(m);
564  double v = (sumL*sumL-m*sumLL)/(sumT*sumL-m*sumLT);
565 
566  Vs.emplace_back(v);
567 
568  for (int t = 0; t!=ntrials; ++t) {
569  sumL = 0;
570  sumT = 0;
571  sumLT = 0;
572  sumLL = 0;
573  for (int i = 0; i!=m; ++i ) {
574  double xi = genR->Uniform(Positions[i].X()-X_err,Positions[i].X()+X_err);
575  double yi = genR->Uniform(Positions[i].Y()-Y_err,Positions[i].Y()+Y_err);
576  double zi = genR->Uniform(Positions[i].Z()-Z_err,Positions[i].Z()+Z_err);
577  double ti = genR->Gaus(Ts[i], Terr[i]);
578  double li = sqrt(xi*xi+yi*yi+zi*zi);
579 
580  sumL += li;
581  sumT += ti;
582  sumLT += li*ti;
583  sumLL += li*li;
584  }
585  v = (sumL*sumL-m*sumLL)/(sumT*sumL-m*sumLT);
586 
587  Vs.emplace_back(v);
588  }
589 
590  //Calculate the mean of v, and standard error:
591 
592  double v_mean = 0;
593  v_err = 0;
594 
595  for (unsigned j = 0; j!=Vs.size(); ++j)
596  v_mean += Vs[j];
597  v_mean /= double(Vs.size());
598 
599  for (unsigned j = 0; j!=Vs.size(); ++j)
600  v_err += (Vs[j]-v_mean)*(Vs[j]-v_mean);
601 
602  v_err = sqrt(v_err)/double(Vs.size());
603 
604  if (v_mean < 0.003 && v_mean > 0)
605  v_mean = 0.003;
606 
607  else if (v_mean > -0.003 && v_mean < 0)
608  v_mean = -0.003;
609 
610  if (v_mean > 30) v_mean = 30;
611  if (v_mean <-30) v_mean = -30;
612 
613  double t0 = Tmean-sumL/(m+0.)/v_mean;
614 
615  Trms = 0;
616  dTsqr = 0;
617 
618  for (int i = 0; i!=m; ++i ) {
619  Trms += (Ts[i]-Tmean)*(Ts[i]-Tmean);
620  dTsqr+= (Ts[i]-t0-Positions[i].Mag()/v_mean)* (Ts[i]-t0-Positions[i].Mag()/v_mean);
621  }
622 
623  Trms = sqrt(Trms/(m+0.));
624  dTsqr = sqrt(dTsqr/(m+0.));
625  return v_mean;
626 }
627 
629  auto nodes = rh_->GetBNEVBList();
630  for (auto node : nodes){
631  std::cout<<node.appname<<" ";
632  }
633 }
634 
float TNS() const
Definition: CellHit.h:46
T max(const caf::Proxy< T > &a, T b)
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
double GetTimeRes(rb::CellHit const &cellhit)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned short Plane() const
Definition: CellHit.h:39
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
Vertical planes which measure X.
Definition: PlaneGeo.h:28
void analyze(const art::Event &evt)
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
Definition: event.h:19
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
double dist2(double dx, double dy)
Float_t Y
Definition: plot.C:38
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
Float_t Z
Definition: plot.C:38
unsigned short Cell() const
Definition: CellHit.h:40
CDPStorage service.
double dy[NP][NC]
double dx[NP][NC]
int evt
std::vector< BNEVB > GetBNEVBList()
Definition: RunHistory.h:398
Collect Geo headers and supply basic geometry functions.
const double j
Definition: BetheBloch.cxx:29
dictionary nodes
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
art::ServiceHandle< nova::dbi::RunHistoryService > rh_
OStream cout
Definition: OStream.cxx:6
FmmTrackerAna(fhicl::ParameterSet const &pset)
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double vfit(std::vector< TVector3 > &Positions, std::vector< double > &Ts, std::vector< double > &Terr, double &v_err, double &Tmean, double &Trms, double &dTsqr)
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
double doca2(double x, double y, double x0, double y0, double dydx)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
float Mag() const
T min(const caf::Proxy< T > &a, T b)
RunNumber_t run() const
Definition: Event.h:77
Float_t X
Definition: plot.C:38
#define W(x)
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)