TCTrack.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////
2 /// \brief A TCHit contains pchits and information related to performing a timing calibration
3 /// \author edniner@indiana.edu
4 /// \date Feb. 2013
5 /////////////////////////////////////////////////////////////////////////////
7 
9 #include "Utilities/func/MathUtil.h"
10 
11 #include <cfloat>
12 #include <climits>
13 
14 namespace caldp
15 {
17  fTotalLength(0.0),
18  fFiberVelocity(0.0)
19  { }
20 
21 
23 
24  //--------------------------------------------------------------------------
26  {
27  //if the nit view was never we want this to fail
28  assert(pchit->View() != geo::kXorY);
29 
30  if (pchit->View() == geo::kX) fXPCHits.push_back(pchit);
31  if (pchit->View() == geo::kY) fYPCHits.push_back(pchit);
32  }
33 
34  //--------------------------------------------------------------------------
35  void TCTrack::Add(const std::vector<art::Ptr<caldp::PCHit> >& pchits)
36  {
37  for (unsigned int n=0; n<pchits.size(); ++n){
38  Add(pchits[n]);
39  }
40  }
41 
42  //--------------------------------------------------------------------------
43  unsigned int TCTrack::NPCHit(geo::View_t view) const
44  {
45  switch(view){
46  case geo::kX: return NXPCHit();
47  case geo::kY: return NYPCHit();
48  case geo::kXorY: return NPCHit();
49  default: assert(0 && "Unknown view");
50  }
51  }
52 
53  //--------------------------------------------------------------------------
55  {
56  switch(view){
57  case geo::kX: return XPCHit(viewIdx);
58  case geo::kY: return YPCHit(viewIdx);
59  case geo::kXorY: return PCHit(viewIdx);
60  default: assert(0 && "Unknown view");
61  }
62  }
63 
64  //--------------------------------------------------------------------------
65  art::Ptr<caldp::PCHit> TCTrack::XPCHit(unsigned int xIdx) const
66  {
67  assert(xIdx < fXPCHits.size());
68  return fXPCHits[xIdx];
69  }
70 
71  //--------------------------------------------------------------------------
72  art::Ptr<caldp::PCHit> TCTrack::YPCHit(unsigned int yIdx) const
73  {
74  assert(yIdx < fYPCHits.size());
75  return fYPCHits[yIdx];
76  }
77 
78  //--------------------------------------------------------------------------
79  art::Ptr<caldp::PCHit> TCTrack::PCHit(unsigned int globalIdx) const
80  {
81  if (globalIdx < fXPCHits.size()) return XPCHit(globalIdx);
82  return YPCHit(globalIdx - NXPCHit());
83  }
84 
85  //--------------------------------------------------------------------------
86  std::vector<art::Ptr<caldp::PCHit> > TCTrack::AllPCHits() const
87  {
88  std::vector<art::Ptr<caldp::PCHit> > all;
89  for(size_t c=0; c<fXPCHits.size(); ++c) all.push_back(fXPCHits[c]);
90  for(size_t c=0; c<fYPCHits.size(); ++c) all.push_back(fYPCHits[c]);
91  return all;
92  }
93 
94  //--------------------------------------------------------------------------
96  {
97  if(pch->View() == geo::kX){
98  for(unsigned int i=0; i<NXPCHit(); ++i){
99  if(fXPCHits[i] == pch) fXPCHits.erase(fXPCHits.begin()+i);
100  break;
101  }
102  }
103  if(pch->View() == geo::kY){
104  for(unsigned int i=0; i<NYPCHit(); ++i){
105  if(fYPCHits[i] == pch) fYPCHits.erase(fYPCHits.begin()+i);
106  break;
107  }
108  }
109  }
110 
111  //--------------------------------------------------------------------------
112  void TCTrack::RemovePCHit(unsigned int globalIdx)
113  {
114  assert(globalIdx < NPCHit());
115  if(globalIdx < NXPCHit()) fXPCHits.erase(fXPCHits.begin()+globalIdx);
116  else fYPCHits.erase(fYPCHits.begin()+globalIdx-NXPCHit());
117  }
118 
119  //--------------------------------------------------------------------------
121  {
122  int ret = INT_MAX;
123  for(unsigned int i=0; i<NPCHit(view); ++i){
124  if(PCHit(view, i)->Plane() < ret) ret = PCHit(view, i)->Plane();
125  }
126  return ret;
127  }
128 
129  //--------------------------------------------------------------------------
131  {
132  int ret = INT_MAX;
133  for(unsigned int i=0; i<NPCHit(view); ++i){
134  if(PCHit(view, i)->Cell() < ret) ret = PCHit(view, i)->Cell();
135  }
136  return ret;
137  }
138 
139  //--------------------------------------------------------------------------
141  {
142  int ret = -INT_MAX;
143  for(unsigned int i=0; i<NPCHit(view); ++i){
144  if(PCHit(view, i)->Plane() > ret) ret = PCHit(view, i)->Plane();
145  }
146  return ret;
147  }
148 
149  //--------------------------------------------------------------------------
151  {
152  int ret = -INT_MAX;
153  for(unsigned int i=0; i<NPCHit(view); ++i){
154  if(PCHit(view, i)->Cell() > ret) ret = PCHit(view, i)->Cell();
155  }
156  return ret;
157  }
158 
159  //--------------------------------------------------------------------------
161  {
162  double ret = DBL_MAX;
163  for(unsigned int i=0; i<NPCHit(view); ++i){
164  if(PCHit(view, i)->TNS() < ret) ret = PCHit(view, i)->TNS();
165  }
166  return ret;
167  }
168 
169  //--------------------------------------------------------------------------
171  {
172  double ret = -DBL_MAX;
173  for(unsigned int i=0; i<NPCHit(view); ++i){
174  if(PCHit(view, i)->TNS() > ret) ret = PCHit(view, i)->TNS();
175  }
176  return ret;
177  }
178 
179  //--------------------------------------------------------------------------
181  {
182  double ret = DBL_MAX;
183  for(unsigned int i=0; i<NPCHit(view); ++i){
184  if(PCHit(view, i)->FlightLen() < ret) ret = PCHit(view, i)->FlightLen();
185  }
186  return ret;
187  }
188 
189  //--------------------------------------------------------------------------
191  {
192  double ret = -DBL_MAX;
193  for(unsigned int i=0; i<NPCHit(view); ++i){
194  if(PCHit(view, i)->FlightLen() > ret) ret = PCHit(view, i)->FlightLen();
195  }
196  return ret;
197  }
198 
199  //--------------------------------------------------------------------------
200  double TCTrack::TNSUncertainty(float pe, bool goodTime, novadaq::cnv::DetId detectorID) const
201  {
202 
203  //This is a fit for the timing resolution as a function of PE.
204  //Current values of the fit come from S14-04-14 data-like MC. This fit can be found in
205  //slide 4 of docdb-11252.
206  //TODO: The constants in this function should be moved into the DB for ease of retrival
207 
208  // uncertainty derived separately for single and multipoint readout. If timing is good, use multip, otherwise
209  // stay with the single point value. These values are for the far detector, not checking detector type yet.
210 
211  double stdev = 144.0; //simple assumption of 500ns/sqrt(12)
212  if (detectorID == novadaq::cnv::kFARDET) {
213  if (goodTime){
214  stdev = 231594.0/(2267.16+pow(pe,2.01011)) + 9.55689; //2014-10-13 fit, multipoint readout from data, docdb-12122, with 460 ns rise time
215  }
216  else {
217  stdev = 848357.0/(5734.75+pow(pe,2.20060)) + 142.221; //2014-06-13 fit, from data
218  }
219  }
220  else if (detectorID == novadaq::cnv::kNEARDET) {
221  if (goodTime){
222  stdev = 2229.53/(77.0043+pow(pe,1.22743)) + 4.89355; //2014-09-15 fit from data with 140 ns rise time, 4.5 us fall time
223  }
224  else {
225  stdev = 6136470.0/(185518.0+pow(pe,2.92360)) + 31.6071; //2014-09-15 singe point fit of data
226  }
227  }
228 
229  return stdev;
230  }
231 
232  //--------------------------------------------------------------------------
233  double TCTrack::ReadoutDistCorrection(float x, novadaq::cnv::DetId detectorID) const
234  {
235 
236  //This is a fit for the timewalk as a function of distance of the hit to the readout electronics.
237  //This effect comes from the fiber being looped, causing early and late arrival of light which can distort the timing fit
238  //The timewalk has been parameterized as an 8 degree polynomial which is ok since extrapolation is never needed.
239  //Current values of this fit come from MC with multipoint readout looking at the difference between expected and reconstructed photon times.
240  //This fit was made simultaneously with a fit of fiberspeed as a function of distance to readout. More details can be found on slide 3 of docdb-11252.
241  //TODO: The constants in this function should be moved into the DB for ease of retrival
242  //TODO: Incorporate timewalk for NearDet, where the effect will be different. For now apply FD values for both detectors.
243 
244  double tw = 0.0; //default to a value of 0 for NDOS
245  if (detectorID == novadaq::cnv::DetId::kFARDET) {
246  tw = -9.63368
247  + x*(0.0185403
248  + x*(-(1.29139e-4)
249  + x*((3.98478e-7)
250  + x*(-(7.66614e-10)
251  + x*((8.9566e-13)
252  + x*(-(6.2088e-16)
253  + x*((2.37032e-19)
254  + x*(-3.81279e-23))))))));
255  }
256  else if (detectorID == novadaq::cnv::DetId::kNEARDET) {
257  /*tw = -9.63368
258  + x*(0.0185403
259  + x*(-(1.29139e-4)
260  + x*((3.98478e-7)
261  + x*(-(7.66614e-10)
262  + x*((8.9566e-13)
263  + x*(-(6.2088e-16)
264  + x*((2.37032e-19)
265  + x*(-3.81279e-23))))))));*/
266  }
267  return tw;
268 
269  }
270 
271  //--------------------------------------------------------------------------
272  double TCTrack::CalcFiberVelocity(float d, novadaq::cnv::DetId detectorID) const
273  {
274 
275  //The effective fiber speed has some dependence on the location of light within the cell. Near the
276  //top of the cell the speed will be slower since light cannot travel both ways in scintillator before being absorbed
277  //in the fiber. This fit is done for now from multipoint MC and is documented in docdb-10863 and the final fit result is on slide 3 of docdb-11252
278  //TODO: This fit is specific to the FD, will have to be redone for ND. For now assume a flat fiber speed for ND
279  double fspeed = 15.7446; //cm/ns
280  if (detectorID == novadaq::cnv::DetId::kFARDET) {
281  fspeed = 15.7446/(((15.7446*1.63517)/d)+1.0);
282  }
283  else if (detectorID == novadaq::cnv::DetId::kNEARDET) {
284  //fspeed = 15.7446/(((15.7446*1.63517)/d)+1.0);
285  }
286 
287  return fspeed;
288 
289  }
290 
291  //--------------------------------------------------------------------------
292  double TCTrack::CalculateFiberVelocity(double timeResid,
293  std::vector<double> *resids,
294  double *icept,
295  double *chisqr,
296  novadaq::cnv::DetId detectorID,
297  bool removeHits,
298  bool overwrite)
299  {
300  //define the speed of light in cm/ns
301  const double sCmPerNsec = 29.9792458;
302 
303  //Containers to store pieces needed for fiber speed fit.
304  //Times are the raw TNS for a cell - TOF correction
305  //Where TOF is the time-of-flight from the start of the track to a hit
306  std::vector<double> time;
307  std::vector<double> rDist;
308  std::vector<double> sterror;
309  std::vector<double> weight;
310 
311  for(unsigned int i=0; i<NPCHit(); ++i){
313  double tns;
314  double pe = pch->PE();
315  double x = pch->ReadDist();
316  double tw = ReadoutDistCorrection(x, detectorID);
317 
318  tns = pch->TNS() + tw;
319  time.push_back(tns - pch->FlightLen()/sCmPerNsec);
320  rDist.push_back(pch->ReadDist());
321  double stdev = TNSUncertainty(pe, pch->GoodTime(), detectorID);
322  sterror.push_back(stdev);
323  weight.push_back(1.0/(stdev*stdev));
324  }
325  double slope, intercept, x1, y1;
326  double fspeed = -99.0;
327  //do iterative fit for the speed light in the optical fiber
328  while (true){
329  //don't let fit go below two points
330  if (rDist.size() < 2) break;
331  resids->clear();
332  *chisqr = util::LinFit(rDist, time, weight, slope, intercept);
333  fspeed = slope;
334  x1 = rDist[0];
335  y1 = slope*x1 + intercept;
336  double maxresid = 0;
337  unsigned int maxpos = 0;
338  for(unsigned int j=0; j<rDist.size(); ++j){
339  double recresid = time[j] - (fspeed*(rDist[j]-x1)+y1);
340  resids->push_back(recresid);
341  if (fabs(recresid/sterror[j]) > fabs(maxresid)){
342  maxresid = fabs(recresid/sterror[j]);
343  maxpos = j;
344  }
345  }
346  //if fit only has two points, dont remove any
347  if (rDist.size() == 2) break;
348  //if there was a hit with a residual above the maximum
349  //remove that hit and recalculate
350  if (fabs(maxresid) > timeResid){
351  weight.erase(weight.begin()+maxpos);
352  rDist.erase(rDist.begin()+maxpos);
353  time.erase(time.begin()+maxpos);
354  sterror.erase(sterror.begin()+maxpos);
355  resids->erase(resids->begin()+maxpos);
356  if (removeHits == true){
357  fDropPCHits.push_back(PCHit(maxpos));
358  RemovePCHit(maxpos);
359  }
360  continue;
361  }
362  //stop iterating if no timing residuals were too large
363  break;
364  }
365  *icept = intercept;
366  if (overwrite == true) fFiberVelocity = 1.0/fspeed;
367  return 1.0/fspeed;
368  }
369 
370  //--------------------------------------------------------------------------
371  double TCTrack::CalculateMuonVelocity(std::vector<double> *resids,
372  double *icept,
373  double *chisqr,
374  novadaq::cnv::DetId detectorID,
375  double fiberVel) const
376  {
377  //define the speed of light in cm/ns
378  const double sCmPerNsec = 29.9792458;
379 
380 
381 
382  //Containers to store pieces needed for fiber speed fit.
383  //Times are the raw TNS for a cell - Fiber speed correction
384  //Where this correction is the time for light to travel up the optical fiber
385  std::vector<double> time;
386  std::vector<double> rDist;
387  std::vector<double> weight;
388 
389  for(unsigned int i=0; i<NPCHit(); ++i){
391  double tns;
392  double pe = pch->PE();
393  double x = pch->ReadDist();
394  double tw = ReadoutDistCorrection(x, detectorID);
395 
396  //if a specific fibr speed is not specified, use parameterization as a function of distance to readout
397  if (fiberVel == 9999.0) fiberVel = CalcFiberVelocity(x, detectorID);
398 
399  tns = pch->TNS() + tw;
400  time.push_back(tns - pch->ReadDist()/fiberVel);
401  rDist.push_back(pch->FlightLen());
402  double stdev = TNSUncertainty(pe, pch->GoodTime(), detectorID);
403 
404  weight.push_back(1.0/(stdev*stdev));
405  }
406 
407  double slope, intercept, x1, y1;
408  double mspeed = -5.0;
409 
410  //do the fit one time
411  resids->clear();
412  //make sure there are at least two points before doing fit;
413  if (rDist.size() < 2) return -5.0;
414  *chisqr = util::LinFit(rDist, time, weight, slope, intercept);
415  mspeed = slope;
416  x1 = rDist[0];
417  y1 = slope*x1 + intercept;
418  for(unsigned int j=0; j<rDist.size(); ++j){
419  double recresid = time[j] - (mspeed*(rDist[j]-x1)+y1);
420  resids->push_back(recresid);
421  }
422 
423  *icept = intercept;
424  //return muon velocity in terms of 1/Beta
425  return mspeed*sCmPerNsec;
426 
427  }
428 
429 
430 
431 } // end namespace calib
double MinTNS(geo::View_t view=geo::kXorY) const
return min pchit time for a specified view
Definition: TCTrack.cxx:160
int MaxCell(geo::View_t view) const
return max cell number in the specified view
Definition: TCTrack.cxx:150
double CalculateFiberVelocity(double timeResid, std::vector< double > *resids, double *icept, double *chisqr, novadaq::cnv::DetId detectorID, bool removeHits=false, bool overwrite=false)
Definition: TCTrack.cxx:292
float TNS() const
Return uncorrected hit time.
Definition: PCHit.h:52
float FlightLen() const
Return path from start of track to current cell.
Definition: PCHit.h:56
def stdev(lst)
Definition: HandyFuncs.py:286
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int NPCHit() const
return total number of hits
Definition: TCTrack.h:41
Float_t y1[n_points_granero]
Definition: compare.C:5
const Var weight
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Float_t x1[n_points_granero]
Definition: compare.C:5
bool GoodTime() const
Return quality of timing fit for cell.
Definition: PCHit.h:54
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
art::Ptr< caldp::PCHit > YPCHit(unsigned int yIdx) const
get the ith y view pchit
Definition: TCTrack.cxx:72
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
double MaxTNS(geo::View_t view=geo::kXorY) const
return max pchit time for a specified view
Definition: TCTrack.cxx:170
int MaxPlane(geo::View_t view) const
return max plane number in the specified view
Definition: TCTrack.cxx:140
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Far Detector at Ash River, MN.
virtual void Add(const art::Ptr< caldp::PCHit > &pchit)
add a pchit to the track
Definition: TCTrack.cxx:25
float PE() const
Return PE value.
Definition: PCHit.h:38
std::vector< art::Ptr< caldp::PCHit > > fDropPCHits
contains collection of pchits dropped by fiber fit
Definition: TCTrack.h:153
void RemovePCHit(const art::Ptr< caldp::PCHit > pch)
remove a pchit from the current cluster
Definition: TCTrack.cxx:95
virtual ~TCTrack()
Definition: TCTrack.cxx:22
double MinPath(geo::View_t view=geo::kXorY) const
return shortest path from start of track to a hit in a view
Definition: TCTrack.cxx:180
Near Detector in the NuMI cavern.
int MinCell(geo::View_t view) const
return min cell number in the specified view
Definition: TCTrack.cxx:130
Float_t d
Definition: plot.C:236
double MaxPath(geo::View_t view=geo::kXorY) const
return longest path from start of track to a hit in a view
Definition: TCTrack.cxx:190
const double j
Definition: BetheBloch.cxx:29
geo::View_t View() const
Return view.
Definition: PCHit.h:36
unsigned int NXPCHit() const
return number of hits in x view
Definition: TCTrack.h:37
Histograms used by attenuation calibration.
art::Ptr< caldp::PCHit > PCHit(geo::View_t view, unsigned int viewIdx) const
get the ith pchit from a specified view
Definition: TCTrack.cxx:54
art::Ptr< caldp::PCHit > XPCHit(unsigned int xIdx) const
get the ith x view pchit
Definition: TCTrack.cxx:65
float ReadDist() const
Return distance to readout, with pigtail.
Definition: PCHit.h:58
static const double sCmPerNsec
double CalcFiberVelocity(float d, novadaq::cnv::DetId detectorID) const
Calculate fiber speed as a function of distance from readout.
Definition: TCTrack.cxx:272
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
assert(nhit_max >=nhit_nbins)
double ReadoutDistCorrection(float x, novadaq::cnv::DetId detectorID) const
Correct for timewalk effect as a function of distance to readout.
Definition: TCTrack.cxx:233
std::vector< art::Ptr< caldp::PCHit > > fYPCHits
Definition: TCTrack.h:150
double fFiberVelocity
optical fiber velocity for track
Definition: TCTrack.h:163
double CalculateMuonVelocity(std::vector< double > *resids, double *icept, double *chisqr, novadaq::cnv::DetId detectorID, double fiberVel=9999.0) const
Definition: TCTrack.cxx:371
Float_t e
Definition: plot.C:35
std::vector< art::Ptr< caldp::PCHit > > AllPCHits() const
return the collection of all pchits
Definition: TCTrack.cxx:86
Definition: fwd.h:28
int MinPlane(geo::View_t view) const
return min plane number in the specified view
Definition: TCTrack.cxx:120
unsigned int NYPCHit() const
return number of hits in y view
Definition: TCTrack.h:39
std::vector< art::Ptr< caldp::PCHit > > fXPCHits
contain collections of pchits
Definition: TCTrack.h:149
double TNSUncertainty(float pe, bool goodTime, novadaq::cnv::DetId detectorID) const
Calculate the uncertainty in the time based on pulse height.
Definition: TCTrack.cxx:200