SMMCluster_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SMMCluster
3 // Module Type: producer
4 // File: SMMCluster_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 "fhiclcpp/ParameterSet.h"
23 
24 #include "Geometry/Geometry.h"
27 
28 #include "NovaDAQConventions/DAQConventions.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/RecoHit.h"
31 #include "RecoBase/Cluster.h"
32 
33 #include "TH3.h"
34 #include "TVector3.h"
35 #include "TTree.h"
36 
37 #include <algorithm>
38 #include <utility>
39 #include <cmath>
40 #include <iostream>
41 
42 namespace zcl {
43  class SMMCluster;
44 }
45 
47 public:
48  explicit SMMCluster(fhicl::ParameterSet const & pset);
49  virtual ~SMMCluster();
50 
51  void produce(art::Event & evt);
52  void beginJob();
53  void endJob();
54 
56  art::PtrVector<rb::CellHit> & hitlist,
57  TVector3 & startC, //C, Z, T, sigmaT... I am just using this as a container
58  TVector3 & stopC);
59 
60  void LinFit(std::vector<TVector3> const & hits,
61  TVector3 const h0,
62  TVector3 & vinfo);
63 
64  double calcW(double z, double w1, double w2, double z1, double z2);
65 
66  private:
67  std::string fClusterInput; ///< Input folder from cluster reco
68  double _sigmaT;
69  double _deltaZ;
70 
71  TTree* RC_Tree;
73  double beta_MC;
74  double Edep_MC;
75  double xi_MC, xf_MC;
76  double yi_MC, yf_MC;
77  double zi_MC, zf_MC;
78  double ti_MC, tf_MC;
80  double ti_RC;
81  double tf_RC;
86  int noise_Xhits, noise_Yhits, noise_XPE, noise_YPE; //This is actually cheater's info
87  double Vx_RC, Vy_RC, Vz_xRC, Vz_yRC;
90 
91 };
92 
94  fClusterInput (p.get< std::string > ("ClusterInput") ),
95  _sigmaT (4687.5),
96  _deltaZ (50)
97 {
98  // Call appropriate Produces<>() functions here.
99  produces< std::vector<rb::Cluster> >();
100 }
101 
103 {
104  // Clean up dynamic memory and other resources here.
105 }
106 
108 {
110 
111  RC_Tree = tfs->make<TTree>("SMMRC_Tree","Information from cluster RC info");
112  RC_Tree -> Branch("triggered", &triggered, "triggered/O");
113  RC_Tree -> Branch("reconstructed", &reconstructed, "reconstructed/O");
114  RC_Tree -> Branch("event", &evtID, "event/I");
115  RC_Tree -> Branch("subRun", &subRunID, "subRun/I");
116  RC_Tree -> Branch("Run", &RunID, "Run/I");
117  RC_Tree -> Branch("ti", &ti_RC, "ti/D");
118  RC_Tree -> Branch("tf", &tf_RC, "tf/D");
119  RC_Tree -> Branch("xi", &xi_RC, "xi/D");
120  RC_Tree -> Branch("yi", &yi_RC, "yi/D");
121  RC_Tree -> Branch("zi", &zi_RC, "zi/D");
122  RC_Tree -> Branch("xf", &xf_RC, "xf/D");
123  RC_Tree -> Branch("yf", &yf_RC, "yf/D");
124  RC_Tree -> Branch("zf", &zf_RC, "zf/D");
125  RC_Tree -> Branch("sigmaxi", &sigmaXi, "sigmaxi/D");
126  RC_Tree -> Branch("sigmayi", &sigmaYi, "sigmayi/D");
127  RC_Tree -> Branch("sigmazi", &sigmaZi, "sigmazi/D");
128  RC_Tree -> Branch("sigmaxf", &sigmaXf, "sigmaxf/D");
129  RC_Tree -> Branch("sigmayf", &sigmaYf, "sigmayf/D");
130  RC_Tree -> Branch("sigmazf", &sigmaZf, "sigmazf/D");
131  RC_Tree -> Branch("missH_X", &missH_X, "missHX/I"); //Missing hits: overall discountinuity in plane
132  RC_Tree -> Branch("missH_Y", &missH_Y, "missHY/I");
133  RC_Tree -> Branch("nHits_X", &nxhits_RC, "nXHits/I");
134  RC_Tree -> Branch("nHits_Y", &nyhits_RC, "nYHits/I");
135  RC_Tree -> Branch("Vx", &Vx_RC, "Vx/D");
136  RC_Tree -> Branch("Vy", &Vy_RC, "Vy/D");
137  RC_Tree -> Branch("Vz_x", &Vz_xRC, "Vz_x/D");
138  RC_Tree -> Branch("Vz_y", &Vz_yRC, "Vz_y/D");
139  RC_Tree -> Branch("sigmaT2_X", &sigmaT2_X, "sigmaT2_X/D");
140  RC_Tree -> Branch("sigmaT2_Y", &sigmaT2_Y, "sigmaT2_Y/D");
141 
142  //Cheater info part
143  RC_Tree -> Branch("betaMC", &beta_MC, "beta/D");
144  RC_Tree -> Branch("edepMC", &Edep_MC, "edep/D");
145  RC_Tree -> Branch("tiMC", &ti_MC, "ti/D");
146  RC_Tree -> Branch("xiMC", &xi_MC, "xi/D");
147  RC_Tree -> Branch("yiMC", &yi_MC, "yi/D");
148  RC_Tree -> Branch("ziMC", &zi_MC, "zi/D");
149  RC_Tree -> Branch("tfMC", &tf_MC, "tf/D");
150  RC_Tree -> Branch("xfMC", &xf_MC, "xf/D");
151  RC_Tree -> Branch("yfMC", &yf_MC, "yf/D");
152  RC_Tree -> Branch("zfMC", &zf_MC, "zf/D");
153  RC_Tree -> Branch("xhits", &nxhits_MC, "xhits/I");
154  RC_Tree -> Branch("yhits", &nyhits_MC, "yhits/I");
155  RC_Tree -> Branch("noise_Xhits", &noise_Xhits, "noise_Xhits/I");
156  RC_Tree -> Branch("noise_Yhits", &noise_Yhits, "noise_Yhits/I");
157  RC_Tree -> Branch("noise_XPE", &noise_XPE, "noise_XPE/D");
158  RC_Tree -> Branch("noise_YPE", &noise_YPE, "noise_YPE/D");
159 
160 }
161 
163 {
164  evtID = evt.id().event();
165  subRunID = evt.subRun();
166  RunID = evt.run();
167 
171 
172  //Get MC info first:
173  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
174  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
175  nxhits_MC = 0;
176  nyhits_MC = 0;
177 
178  const sim::Particle *p = (*i).second;
179  if (p->PdgCode()!=42 ) continue;
180  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
181  unsigned ntruehits_Tot = flshits.size();
182  if ( ntruehits_Tot==0 ) continue;
183  double momentum2 = (p->P())*(p->P());
184  double mass2 = (p->Mass())*(p->Mass());
185  beta_MC = sqrt(momentum2 / (momentum2+mass2) );
186 
187  ti_MC = flshits[0].GetEntryT();
188  xi_MC = flshits[0].GetEntryX();
189  yi_MC = flshits[0].GetEntryY();
190  zi_MC = flshits[0].GetEntryZ();
191  tf_MC = flshits[0].GetExitT();
192  xf_MC = flshits[0].GetExitX();
193  yf_MC = flshits[0].GetExitY();
194  zf_MC = flshits[0].GetExitZ();
195  Edep_MC = flshits[0].GetEdep();
196 
197  unsigned int planeID_i = flshits[0].GetPlaneID();
198  unsigned int cellID_i = flshits[0].GetCellID();
199  unsigned int planeID_f = planeID_i;
200  unsigned int cellID_f = cellID_i;
201 
202  if (planeID_i%2) ++nxhits_MC;
203  else ++nyhits_MC;
204 
205  for (size_t h =1; h!= flshits.size(); ++h) {
206  double ti = flshits[h].GetEntryT();
207  double tf = flshits[h].GetExitT();
208  if (ti < ti_MC) {
209  ti_MC = ti;
210  xi_MC = flshits[h].GetEntryX();
211  yi_MC = flshits[h].GetEntryY();
212  zi_MC = flshits[h].GetEntryZ();
213  planeID_i = flshits[h].GetPlaneID();
214  cellID_i = flshits[h].GetCellID();
215  }
216  if (tf > tf_MC) {
217  tf_MC = tf;
218  xf_MC = flshits[h].GetExitX();
219  yf_MC = flshits[h].GetExitY();
220  zf_MC = flshits[h].GetExitZ();
221  planeID_f = flshits[h].GetPlaneID();
222  cellID_f = flshits[h].GetCellID();
223  }
224 
225  if (flshits[h].GetPlaneID() % 2) ++nxhits_MC;
226  else ++nyhits_MC;
227 
228  Edep_MC += flshits[h].GetEdep();
229  }
230 
231  //By now, you have got the local coordinate, you have to translate it to world:
232  const geo::CellGeo* startCell = geom->Plane(planeID_i)->Cell(cellID_i);
233  double local[3], world[3];
234  local[0] = xi_MC;
235  local[1] = yi_MC;
236  local[2] = zi_MC;
237  startCell->LocalToWorld(local,world);
238  xi_MC = world[0];
239  yi_MC = world[1];
240  zi_MC = world[2];
241 
242  const geo::CellGeo* stopCell = geom->Plane(planeID_f)->Cell(cellID_f);
243  local[0] = xf_MC;
244  local[1] = yf_MC;
245  local[2] = zf_MC;
246  stopCell->LocalToWorld(local,world);
247  xf_MC = world[0];
248  yf_MC = world[1];
249  zf_MC = world[2];
250 
251  //1 monopole per event, so you can break out from the loop:
252  break;
253  }
254 
255  //MC info ready, fill the tree:
256 
257 
258  //=============================MC/RC Separation Line============================================
259 
260  //initialization of triggering and reconstruction status:
261  triggered = false;
262  reconstructed = false;
263 
265  evt.getByLabel(fClusterInput, slices);
266 
267  //Devide the clusters from different views:
268  std::vector<rb::Cluster> Xclusters;
269  std::vector<rb::Cluster> Yclusters;
270 
271  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
272  const rb::Cluster& slice = (*slices)[sliceIdx];
273  if ((slice.Cell(0)) -> View()==geo::kX) Xclusters.emplace_back(slice);
274  else Yclusters.emplace_back(slice);
275  }
276 
277  if (Xclusters.empty())
278  RC_Tree -> Fill();
279 
280  else triggered = true;
281 
282  /*
283  std::cout<<"X sp clusters: "<<Xclusters.size()
284  <<" Y sp clusters: "<<Yclusters.size()<<std::endl;
285  */
286  std::unique_ptr<std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);
287 
288  //Your mission: merge 2D tracks and meanwhile record the information from both cheater and RC
289  for (unsigned xsid = 0 ; xsid!= Xclusters.size(); ++xsid) {
290  TVector3 X_start, X_stop;
291  art::PtrVector<rb::CellHit> xhitlist = Xclusters[xsid].AllCells();
292  missH_X = MissC(geom, xhitlist, X_start, X_stop);
293  nxhits_RC = xhitlist.size();
294 
295  for (unsigned ysid = 0; ysid!= Yclusters.size(); ++ysid) {
296  TVector3 Y_start, Y_stop;
297  art::PtrVector<rb::CellHit> yhitlist = Yclusters[ysid].AllCells();
298  missH_Y = MissC(geom, yhitlist, Y_start, Y_stop);
299  nyhits_RC = yhitlist.size();
300 
301  //Check whether to merge these 2:
302  if (X_start.Z() > Y_start.Z()+_sigmaT || X_start.Z() < Y_start.Z()-_sigmaT) continue;
303  if (X_stop.Z() > Y_stop.Z() +_sigmaT || X_stop.Z() < Y_stop.Z() -_sigmaT) continue;
304  if (X_start.Y() > Y_start.Y()+_deltaZ || X_start.Y() < Y_start.Y()-_deltaZ) continue;
305  if (X_stop.Y() > Y_stop.Y() +_deltaZ || X_stop.Y() < Y_stop.Y() -_deltaZ) continue;
306 
307  //if all these requirements have been reached, merge!
308  reconstructed = true;
309  rb::Cluster SMMcluster3D;
310  std::vector<TVector3> xhits;
311 
312  noise_Xhits = 0;
313  noise_XPE = 0;
314  for (int i= 0; i!=nxhits_RC; ++i) {
315  const art::Ptr<rb::CellHit> chit = xhitlist[i];
316 
317  double xyz[3];
318  int p = chit->Plane();
319  int c = chit->Cell();
320  geom->Plane(p)->Cell(c)->GetCenter(xyz);
321 
322  rb::RecoHit rhit(calib->MakeRecoHit( *chit,
323  calcW(xyz[2],
324  Y_start.X(),
325  Y_stop.X(),
326  Y_start.Y(),
327  Y_stop.Y()) ) );
328 
329  if (!rhit.IsCalibrated()) continue;
330 
331  TVector3 hit_st(xyz[0],xyz[2],rhit.T());
332  xhits.emplace_back(hit_st);
333  SMMcluster3D.Add(chit);
334  if (bt->IsNoise(chit)) {
335  ++noise_Xhits;
336  noise_XPE += rhit.PECorr();
337  }
338  }
339 
340  std::vector<TVector3> yhits;
341  noise_Yhits = 0;
342  noise_YPE = 0;
343  for (int i= 0; i!=nyhits_RC; ++i) {
344  const art::Ptr<rb::CellHit> chit = yhitlist[i];
345  double xyz[3];
346  int p = chit->Plane();
347  int c = chit->Cell();
348  geom->Plane(p)->Cell(c)->GetCenter(xyz);
349 
350  rb::RecoHit rhit(calib->MakeRecoHit(*chit,
351  calcW(xyz[2],
352  X_start.X(),
353  X_stop.X(),
354  X_start.Y(),
355  X_stop.Y()) ) );
356 
357  if (!rhit.IsCalibrated()) continue;
358 
359  TVector3 hit_st(xyz[1],xyz[2],rhit.T());
360  yhits.emplace_back(hit_st);
361  SMMcluster3D.Add(chit);
362  if (bt->IsNoise(chit)) {
363  ++noise_Yhits;
364  noise_YPE += rhit.PECorr();
365  }
366  }
367 
368  clustercol->emplace_back(SMMcluster3D);
369 
370  //Now, do the vx ,vy,vz fit:
371  TVector3 X_vinfo, Y_vinfo;
372  LinFit(xhits,X_start,X_vinfo);
373  LinFit(yhits,Y_start,Y_vinfo);
374 
375  Vx_RC = X_vinfo.X();
376  Vz_xRC = X_vinfo.Y();
377  sigmaT2_X = X_vinfo.Z();
378  Vy_RC = Y_vinfo.X();
379  Vz_yRC = Y_vinfo.Y();
380  sigmaT2_Y = Y_vinfo.Z();
381 
382  //By now, the 3D cluster has been reconstructed, calculate the info and fill the RC tree:
383  //The remaining infomation to be calculated:
384  //sigmaX,Y,Z,T, (XYZ)_if
385  sigmaXi = 1.8; //1.8cm is the half size of a cell width
386  sigmaYi = 1.8;
387 
388  double sigmaT2 = (sigmaT2_X + sigmaT2_Y)/(nxhits_RC+nyhits_RC+0.);
389 
390  //If the track is forward going:
391  if (Vz_xRC+Vz_yRC > 0) {
392  //First, start point:
393  if (X_start.Y() < Y_start.Y()) {
394  //This means, the first hit is in X view:
395  xi_RC = X_start.X();
396  zi_RC = X_start.Y();
397  sigmaZi = Y_start.Y()-zi_RC-6.6;
398  yi_RC = Y_start.X() - (Y_start.Y() - zi_RC)*Vy_RC/Vz_yRC;
399  ti_RC = X_start.Z();
400  sigmaYi = sqrt(sigmaYi*sigmaYi + Vy_RC*Vy_RC*sigmaT2);
401  }
402  else {
403  //This means, the first hit is in Y view:
404  yi_RC = Y_start.X();
405  zi_RC = Y_start.Y();
406  sigmaZi = X_start.Y() - zi_RC - 6.6;
407  xi_RC = X_start.X() - (X_start.Y() - zi_RC)*Vx_RC/Vz_xRC;
408  ti_RC = Y_start.Z();
409  sigmaXi = sqrt(sigmaXi*sigmaXi + Vx_RC*Vx_RC*sigmaT2);
410  }
411 
412  //Now, stop point:
413  if (X_stop.Y() > Y_stop.Y()) {
414  //This means, the last hit is in X view:
415  xf_RC = X_stop.X();
416  zf_RC = X_stop.Y();
417  sigmaZf = zf_RC - Y_stop.Y() - 6.6;
418  yf_RC = Y_stop.X() + (zf_RC-Y_stop.Y() )*Vy_RC/Vz_yRC;
419  ti_RC = X_stop.Z();
420  sigmaYi = sqrt(sigmaYi*sigmaYi+ Vy_RC*Vy_RC*sigmaT2);
421  }
422  else {
423  //This means, the last hit is in Y view:
424  yf_RC = Y_stop.X();
425  zf_RC = Y_stop.Y();
426  sigmaZf = zf_RC-Y_stop.Y()-6.6;
427  xf_RC = X_stop.X() + (X_stop.Y() - zi_RC)*Vx_RC/Vz_xRC;
428  tf_RC = Y_stop.Z();
429  sigmaXi = sqrt(sigmaXi*sigmaXi+ Vx_RC*Vx_RC*sigmaT2);
430  }
431  }
432 
433  //Or, this track is going backwards:
434  else {
435  if (X_start.Y() > Y_start.Y()) {
436  //This means, the first hit is in X view: Xz > Yz
437  xi_RC = X_start.X();
438  zi_RC = X_start.Y();
439  sigmaZi = zi_RC-Y_start.Y()-6.6;
440  yi_RC = Y_start.X() + (sigmaZi+6.6)*Vy_RC/Vz_yRC;
441  ti_RC = X_start.Z();
442  sigmaYi = sqrt(sigmaYi*sigmaYi + Vy_RC*Vy_RC*sigmaT2);
443  }
444  else {
445  //This means, the first hit is in Y view: Yz > Xz
446  yi_RC = Y_start.X();
447  zi_RC = Y_start.Y();
448  sigmaZi = zi_RC - X_start.Y() - 6.6;
449  xi_RC = X_start.X() + (sigmaZi+6.6)*Vx_RC/Vz_xRC;
450  ti_RC = Y_start.Z();
451  sigmaXi = sqrt(sigmaXi*sigmaXi + Vx_RC*Vx_RC*sigmaT2);
452  }
453 
454  //Now, stop point:
455  if (X_stop.Y() < Y_stop.Y()) {
456  //This means, the last hit is in X view: Xz < Yz
457  xf_RC = X_stop.X();
458  zf_RC = X_stop.Y();
459  sigmaZf = Y_stop.Y() - zf_RC - 6.6;
460  yf_RC = Y_stop.X() - (sigmaZf+6.6)*Vy_RC/Vz_yRC;
461  ti_RC = X_stop.Z();
462  sigmaYi = sqrt(sigmaYi*sigmaYi+ Vy_RC*Vy_RC*sigmaT2);
463  }
464  else {
465  //This means, the last hit is in Y view: Xz > Yz
466  yf_RC = Y_stop.X();
467  zf_RC = Y_stop.Y();
468  sigmaZf = Y_stop.Y()-zf_RC-6.6;
469  xf_RC = X_stop.X() - (sigmaZf+6.6)*Vx_RC/Vz_xRC;
470  tf_RC = Y_stop.Z();
471  sigmaXi = sqrt(sigmaXi*sigmaXi+ Vx_RC*Vx_RC*sigmaT2);
472  }
473  }
474 
475  sigmaZi = sqrt(sigmaZi*sigmaZi+7.9524);
476  sigmaZf = sqrt(sigmaZf*sigmaZf+7.9524);
477 
478  RC_Tree -> Fill();
479 
480  }
481  }
482 
483  //If triggered but not reconstructed, fill the mc info anyway:
484  if (triggered && !reconstructed) RC_Tree -> Fill();
485 
486  //std::cout<<"beta: "<< beta_MC<<std::endl;
487  evt.put(std::move(clustercol));
488 } //end production function
489 
491 
492 }
493 
494 //Predetermine the first and last hit in this 2D cluster
496  art::PtrVector<rb::CellHit> & hitlist,
497  TVector3 & startC, //x(y)/cm, z/cm, tns
498  TVector3 & stopC
499  )
500 {
501  int missHits = 0;
502  int nhits = hitlist.size();
503  rb::SortByPlane(hitlist);
504 
505  const art::Ptr<rb::CellHit> chit0 = hitlist[0];
506  const art::Ptr<rb::CellHit> chitf = hitlist[nhits-1];
507 
508  float t0 = chit0->TNS();
509  int c0 = chit0->Cell();
510  int p0 = chit0->Plane();
511  float tf = chitf->TNS();
512  int cf = chitf->Cell();
513  int pf = chitf->Plane();
514 
515  int p_pre = p0;
516 
517  if (c0>cf) //largest c corresponds to c0 in that plane
518  {
519  for (int i=0; i!= nhits; ++i) {
520  const art::Ptr<rb::CellHit> chiti = hitlist[i];
521  int pi = chiti->Plane();
522  int ci = chiti->Cell();
523  if (pi==p0 && ci>c0) c0 = ci;
524  if (pi==pf && ci<cf) c0 = cf;
525  if (pi>p_pre+2) {
526  missHits+=(pi-p_pre)/2;
527  }
528  p_pre=pi;
529  }
530  }
531 
532  else
533  {
534  for (int i=0; i!= nhits; ++i) {
535  const art::Ptr<rb::CellHit> chiti = hitlist[i];
536  int pi = chiti->Plane();
537  int ci = chiti->Cell();
538  if (pi==p0 && ci<c0) c0 = ci;
539  if (pi==pf && ci>cf) c0 = cf;
540  if (pi>p_pre+2) missHits+=(pi-p_pre)/2;
541  p_pre=pi;
542  }
543  }
544 
545  double xyz0[3],xyzf[3];
546  geom->Plane(p0)->Cell(c0)->GetCenter(xyz0);
547  geom->Plane(pf)->Cell(cf)->GetCenter(xyzf);
548 
549  double x0, xf;
550  if (chit0->View()==geo::kX){
551  x0 = xyz0[0];
552  xf = xyzf[0];
553  }
554 
555  else {
556  x0 = xyz0[1];
557  xf = xyzf[1];
558  }
559 
560  if (t0<tf) {
561  startC.SetXYZ(x0,xyz0[2],t0);
562  stopC.SetXYZ(xf,xyzf[2],tf);
563  }
564  else {
565  startC.SetXYZ(xf,xyzf[2],tf);
566  stopC.SetXYZ(x0,xyz0[2],t0);
567  }
568 
569  return missHits;
570 }
571 
572 void zcl::SMMCluster::LinFit(std::vector<TVector3> const & hits,
573  TVector3 const h0, //First hit
574  TVector3 & vinfo) //vx, vy, sigmaT2
575 {
576  float xt = 0;
577  float x0 = h0.X();
578  float yt = 0;
579  float y0 = h0.Y();
580  float tt = 0;
581  float t0 = h0.Z();
582 
583  for (unsigned i = 0; i!=hits.size(); ++i) {
584  xt += (hits[i].X()-x0)*(hits[i].Z()-t0);
585  yt += (hits[i].Y()-y0)*(hits[i].Z()-t0);
586  tt += (hits[i].Z()-t0)*(hits[i].Z()-t0);
587  }
588 
589  double vx = xt/tt;
590  double vy = yt/tt;
591  double sigmaT2=0;
592 
593  for (unsigned i = 0; i!=hits.size(); ++i) {
594  float A = hits[i].Z()-t0-(hits[i].X()-x0)/vx;
595  float B = hits[i].Z()-t0-(hits[i].Y()-y0)/vy;
596  sigmaT2 += A*A+B*B;
597  }
598 
599  vinfo.SetXYZ(vx,vy,sigmaT2);
600 }
601 
602 double zcl::SMMCluster::calcW(double z, double w1, double w2, double z1, double z2) {
603  double w = (w2*(z-z1)+w1*(z2-z))/(z2-z1);
604  if (w > 761) w = 761;
605  else if (w < -761) w = -761;
606 
607  return w;
608 }
609 
610 
float TNS() const
Definition: CellHit.h:46
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
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...
double Mass() const
Definition: MCParticle.h:238
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
SMMCluster(fhicl::ParameterSet const &pset)
A collection of associated CellHits.
Definition: Cluster.h:47
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Timing fit.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int MissC(art::ServiceHandle< geo::Geometry > const &geom, art::PtrVector< rb::CellHit > &hitlist, TVector3 &startC, TVector3 &stopC)
int TrackId() const
Definition: MCParticle.h:209
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
Definition: type_traits.h:56
CDPStorage service.
void hits()
Definition: readHits.C:15
#define local
Definition: gzguts.h:107
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double P(const int i=0) const
Definition: MCParticle.h:233
int evt
Definition: View.py:1
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
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
void produce(art::Event &evt)
static const double A
Definition: Units.h:82
double calcW(double z, double w1, double w2, double z1, double z2)
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fClusterInput
Input folder from cluster reco.
void geom(int which=0)
Definition: geom.C:163
RunNumber_t run() const
Definition: Event.h:77
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Float_t w
Definition: plot.C:20
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
void LinFit(std::vector< TVector3 > const &hits, TVector3 const h0, TVector3 &vinfo)
enum BeamMode string