FMMTracker_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FMMTracker
3 // Module Type: producer
4 // File: FMMTracker_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 
36 #include "TVector3.h"
37 #include "TTree.h"
38 
39 #include <algorithm>
40 #include <utility>
41 #include <cmath>
42 #include <iostream>
43 
44 namespace zcl {
45  class FMMTracker;
46 }
47 
49 public:
50  explicit FMMTracker(fhicl::ParameterSet const & pset);
51  virtual ~FMMTracker();
52 
53  void produce(art::Event & evt);
54 
56  const art::PtrVector<rb::CellHit> & hitlist,
57  int view,
58  rb::Track & track,
59  std::vector<TVector3> & poslist);
60 
61  double recoT(TVector3 position,
62  std::vector<TVector3>& poslist,
63  std::vector<double>& PEs,
64  std::vector<double>& tlist);
65 
66  double calcW(double z, double w1, double w2, double z1, double z2);
67 
68  private:
69  std::string fClusterInput; ///< Input folder from cluster reco
70  double _dSqrMax; //Determine whether this hit will be included in the track
71  int _PEcut; //whether this hit is qualified to place a weight in the timing of start or stop point
72  unsigned _minHits;
73  //float _deltaT;
74  //float _deltaZ;
75 
76 };
77 
79  fClusterInput (p.get< std::string > ("ClusterInput") ),
80  _dSqrMax (200),
81  _PEcut (500),
82  _minHits (3)//, //At least 3 hits on each view
83  // _deltaT (1500), //ns
84  // _deltaZ (30) //cm
85 {
86  // Call appropriate Produces<>() functions here.
87  produces< std::vector<rb::Track> >();
88 }
89 
91 {
92  // Clean up dynamic memory and other resources here.
93 }
94 
96 {
99 
100  //Now, start reco tracks!
102  evt.getByLabel(fClusterInput, slices);
103  std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
104 
105  //Instead of splitting the hits into 2 views, just directly take the entire 3D slice:
106  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
107  const rb::Cluster& slice = (*slices)[sliceIdx];
108 
109  const art::PtrVector<rb::CellHit>& Xhitlist = slice.XCells();
110  //During the stage of the MMSlicer, the hits should already be sorted by plane number
111  const art::PtrVector<rb::CellHit>& Yhitlist = slice.YCells();
112 
113  int px_min = 999;
114  int px_max = -999;
115  double tx_min = 9e9;
116  double tx_max = -9e9;
117 
118  for (unsigned i=0; i!= Xhitlist.size(); ++i) {
119  const art::Ptr<rb::CellHit> chit = Xhitlist[i];
120  if (tx_min > chit->TNS()) tx_min = chit->TNS();
121  if (tx_max < chit->TNS()) tx_max = chit->TNS();
122  if (px_min > chit->Plane()) px_min = chit->Plane();
123  if (px_max < chit->Plane()) px_max = chit->Plane();
124  }
125 
126  int py_min = 999;
127  int py_max = -999;
128  double ty_min = 9e9;
129  double ty_max = -9e9;
130 
131  for (unsigned i=0; i!= Yhitlist.size(); ++i) {
132  const art::Ptr<rb::CellHit> chit = Yhitlist[i];
133  if (ty_min > chit->TNS()) ty_min = chit->TNS();
134  if (ty_max < chit->TNS()) ty_max = chit->TNS();
135  if (py_min > chit->Plane()) py_min = chit->Plane();
136  if (py_max < chit->Plane()) py_max = chit->Plane();
137  }
138 
139  //Matching test:
140  if (std::min(tx_max, ty_max) - std::max(tx_min,ty_min) < 0) continue;
141  if (std::min(px_max, py_max) - std::max(px_min,py_min) < 3) continue;
142 
143  // rb::SortByPlane(Yhitlist);
144  rb::Track xtrack, ytrack;
145  // TVector3 pX1, pX2, pY1, pY2;
146  std::vector<TVector3> Xp_list,Yp_list; //X(Y),Z
147  std::vector<TVector3> trjList; //X,Y,Z
148  std::vector<double> t_list;
149 
150  //Either X or Y track is not reconstructed, just skip to the next slice
151  // geom: hitlist: view: track: positions list:
152  if (!Tracking2D(geom,Xhitlist,0,xtrack,Xp_list) )
153  continue;
154 
155  if (!Tracking2D(geom,Yhitlist,1,ytrack,Yp_list) )
156  continue;
157 
158  //Merging: and setting W position for each reco hit, to get the calibrated time and energy
159  rb::Track FMMtrack3D;
160  //Reconstruct the 3D info of the start and end points
161 
162  TVector3 lowZ_P, highZ_P;
163  // double lowZ_T, highZ_T;
164 
165  unsigned nxhits_RC = Xp_list.size();
166  unsigned nyhits_RC = Yp_list.size();
167 
168  for (unsigned i=0; i!= nxhits_RC; ++i) {
169  const art::Ptr<rb::CellHit> chit = xtrack.Cell(i);
170 
171  rb::RecoHit rhit(calib->MakeRecoHit(*chit,
172  calcW(Xp_list[i].Y(),
173  Yp_list[0].X(),
174  Yp_list[nyhits_RC-1].X(),
175  Yp_list[0].Y(),
176  Yp_list[nyhits_RC-1].Y() ) ) );
177  if (!rhit.IsCalibrated() )
178  continue;
179 
180  trjList.emplace_back(TVector3(rhit.X(),rhit.Y(),rhit.Z()));
181  t_list.emplace_back(rhit.T());
182  FMMtrack3D.AppendTrajectoryPoint(TVector3(rhit.X(),rhit.Y(),rhit.Z()) );
183 
184  // std::cout<<"Xview: "<<"X: "<<rhit.X()<<" Y:"<<rhit.Y()<<" Z: "<<rhit.Z()<<std::endl;
185 
186  FMMtrack3D.Add(chit);
187  }
188 
189  //X view finished, predetermine the lowZ and highZ:
190  /*
191  lowZ_P = trjList[0];
192  lowZ_T = t_list[0];
193  highZ_P = trjList[nxhits_RC-1];
194  highZ_T = t_list[nxhits_RC-1];
195  */
196  // std::cout<<"Y hits on track: "<<Yp_list.size()<<std::endl;
197  for (unsigned i=0; i!= nyhits_RC; ++i) {
198  const art::Ptr<rb::CellHit> chit = ytrack.Cell(i);
199  rb::RecoHit rhit(calib->MakeRecoHit(*chit,
200  calcW(Yp_list[i].Y(),
201  Xp_list[0].X(),
202  Xp_list[nxhits_RC-1].X(),
203  Xp_list[0].Y(),
204  Xp_list[nxhits_RC-1].Y() ) ) );
205 
206  if (!rhit.IsCalibrated() )
207  continue;
208  // std::cout<<"Yview: "<<"X: "<<rhit.X()<<" Y:"<<rhit.Y()<<" Z: "<<rhit.Z()<<std::endl;
209 
210  trjList.emplace_back(TVector3(rhit.X(),rhit.Y(),rhit.Z()));
211  t_list.emplace_back(rhit.T());
212  FMMtrack3D.AppendTrajectoryPoint(TVector3(rhit.X(),rhit.Y(),rhit.Z()) );
213  FMMtrack3D.Add(chit);
214  }
215  /*
216  if (lowZ_P.Z() > trjList[nxhits_RC].Z()) {
217  lowZ_P = trjList[nxhits_RC];
218  lowZ_T = t_list[nxhits_RC];
219  }
220 
221  if (highZ_P.Z() < trjList[nxhits_RC+nyhits_RC-1].Z()) {
222  highZ_P = trjList[nxhits_RC+nyhits_RC-1];
223  highZ_T = t_list[nxhits_RC+nyhits_RC-1];
224  }
225 
226  //Now, determine the start and end time (simultaneously you know which is the start point):
227  lowZ_T = recoT(lowZ_P,trjList,pe_list,t_list);
228  highZ_T = recoT(highZ_P,trjList,pe_list,t_list);
229 
230 
231  if (lowZ_T < highZ_T) {
232  FMMtrack3D.SetStart(lowZ_P);
233  FMMtrack3D.SetDir((highZ_P-lowZ_P).Unit());
234  }
235  else {
236  FMMtrack3D.SetStart(highZ_P);
237  FMMtrack3D.SetDir((lowZ_P-highZ_P).Unit());
238  }
239  */
240  trackcol->emplace_back(FMMtrack3D);
241  // std::cout<<"reconstructed 3D tracks: "<<trackcol->size()<<std::endl;
242  }
243  evt.put(std::move(trackcol));
244 } //end production function
245 
246 
247 //Predetermine the first and last hit in this 2D track
249  const art::PtrVector<rb::CellHit> & hitlist,
250  int view, //For index convenience
251  rb::Track & track,
252  std::vector<TVector3> & plist //X(Y), Z, T
253  )
254 {
255  double xyz[3];
256  std::vector<double> x,y,w,t;
257 
258  if (hitlist.size()<_minHits) return false;
259 
260  //First time loop: store xyw
261  for (unsigned i=0; i!= hitlist.size(); ++i) {
262  const art::Ptr<rb::CellHit> chit = hitlist[i];
263 
264 
265  w.emplace_back(chit->ADC() );
266  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
267  x.emplace_back(xyz[2]);
268  y.emplace_back(xyz[view]);
269  t.emplace_back(chit->TNS());
270  }
271  //Doing Linfit to determine the start and end position:
272  //But actually I just use this to drop non-in-line hits in the cluster
273  double x1,y1,x2,y2,dSq;
274  dSq=geo::LinFitMinDperp(x,y,w,&x1,&x2,&y1,&y2);
275 
276  //Second time loop: drop hits not on the line
277  for (unsigned i=0; i!= hitlist.size(); ++i) {
278  dSq = geo::DsqrToLine(x[i],y[i],x1,x2,y1,y2);
279  if (dSq<_dSqrMax) {
280  track.Add(hitlist[i]);
281  plist.emplace_back( TVector3(y[i], x[i], t[i]) );
282  }
283  }
284 
285  if (track.NCell()<_minHits) return false;
286  else return true;
287 
288 }
289 
291  std::vector<TVector3>& poslist,
292  std::vector<double>& PEs,
293  std::vector<double>& tlist)
294 {
295  double T = 0;
296  double wt = 0;
297  for (unsigned i=0; i!= PEs.size(); ++i) {
298  if (PEs[i] < _PEcut) continue;
299  double weight = exp(-(poslist[i]-position).Mag2());
300  T = weight*tlist[i];
301  wt += weight;
302  }
303  T = T/wt;
304  return T;
305 }
306 
307 double zcl::FMMTracker::calcW(double z, double w1, double w2, double z1, double z2) {
308  double w = (w2*(z-z1)+w1*(z2-z))/(z2-z1);
309  if (w > 761) w = 761;
310  else if (w < -761) w = -761;
311 
312  return w;
313 }
314 
315 
316 
float TNS() const
Definition: CellHit.h:46
T max(const caf::Proxy< T > &a, T b)
bool Tracking2D(art::ServiceHandle< geo::Geometry > const &geom, const art::PtrVector< rb::CellHit > &hitlist, int view, rb::Track &track, std::vector< TVector3 > &poslist)
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const art::PtrVector< rb::CellHit > & XCells() const
Get all cells from the x-view.
Definition: Cluster.h:124
Float_t y1[n_points_granero]
Definition: compare.C:5
const Var weight
unsigned short Plane() const
Definition: CellHit.h:39
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
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
std::string fClusterInput
Input folder from cluster reco.
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
CDPStorage service.
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
const art::PtrVector< rb::CellHit > & YCells() const
Get all cells from the x-view.
Definition: Cluster.h:126
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
int evt
Collect Geo headers and supply basic geometry functions.
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float Mag2() const
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double recoT(TVector3 position, std::vector< TVector3 > &poslist, std::vector< double > &PEs, std::vector< double > &tlist)
void geom(int which=0)
Definition: geom.C:163
FMMTracker(fhicl::ParameterSet const &pset)
std::initializer_list< hep_hpc::hdf5::PropertyList > plist
double calcW(double z, double w1, double w2, double z1, double z2)
double T
Definition: Xdiff_gwt.C:5
void produce(art::Event &evt)
T min(const caf::Proxy< T > &a, T b)
Float_t X
Definition: plot.C:38
Float_t w
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)