TimingCalibration_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: TimingCalibration
3 // Module Type: producer
4 // File: TimingCalibration_module.cc
5 //
6 // Generated at Mon Feb 11 09:01:32 2013 by Evan Niner using artmod
7 // from art v1_02_06.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //NOvASoft includes
11 #include "Geometry/Geometry.h"
12 #include "RecoBase/CellHit.h"
13 #include "RecoBase/Cluster.h"
14 #include "RecoBase/Track.h"
15 #include "Utilities/func/MathUtil.h"
19 #include "NovaDAQConventions/DAQConventions.h"
20 #include "DAQChannelMap/DAQChannelMap.h"
21 #include "Calibrator/Calibrator.h"
22 #include "MCCheater/BackTracker.h"
23 
24 //framework includes
32 #include "Utilities/AssociationUtil.h"
35 
36 
37 //ROOT includes
38 #include "TTree.h"
39 #include "TH2F.h"
40 #include "TH1F.h"
41 
42 namespace calib {
43 
45  public:
46  explicit TimingCalibration(fhicl::ParameterSet const & p);
47  virtual ~TimingCalibration();
48 
49  void produce(art::Event & e) override;
50 
51  void reconfigure(const fhicl::ParameterSet& pset);
52  void beginJob() override;
53  void clearNTuple();
54  void beginRun(art::Run& run) override;
55 
56 
57  private:
58 
59  std::string fTCTrackLabel; ///< Where to find calibration tracks
60  std::string fTrackLabel; ///< Where to find reco Tracks
61 
62  //variables to select track quality
63  int fMinDCMHits; ///<minimum number of hits on a dcm to keep results
64  int fMinDCM; ///<minimum number of good dcms on a track
65  bool fDebug; ///<Add exta debugging info to ntuple
66 
67  float fFiberSpeed; ///<fixed fiberspeed
68  bool fRemoveOutlier; ///<remove outlier hits in DCMs based on timing resolution
69  int fRemoveThreshold; ///<number of sigma away from mean to remove hits
70  bool fCorrectTimewalk; ///<adjust hit times based on timewalk effect for distance to readout
71  bool fAdjustSpeed; ///<adjust fiberspeed based on distance to readout instead of assuming flat speed
72 
73 
74  //map relating unique dcm to collection of pchits
75  std::map<int, std::vector<art::Ptr<caldp::PCHit> > > fDCMHitMap;
76 
77  //ntuple variables
78  //event statistics
79  int run;
80  int subrun;
81  int event;
82  int cosmictr;
83  //track start coordinate
84  double sX;
85  double sY;
86  double sZ;
87  //track end coordinate
88  double eX;
89  double eY;
90  double eZ;
91 
92  double plength;
93 
94  //min, max plane
95  int minPl;
96  int maxPl;
97  //min, max cell
98  int minPlX;
99  int minPlY;
100  int maxPlX;
101  int maxPlY;
102  //number of cells used in final fiber speed fit
106 
107  uint32_t evtTime;
108 
109  //vectors filled for each track, one entry per dcm
110  std::vector<float> dcmnum; ///<dcm index
111  std::vector<float> dcmhits; ///<number of hits in dcm
112  std::vector<float> dcmsimtime; ///<simultaneous hit time in dcm
113 
114  //vectors filled for each track, one entry per hit, sorted by dcm
115  std::vector<float> cellnum; ///<cell number
116  std::vector<float> planenum; ///<plane number
117  std::vector<float> dcmsimresid; ///<difference between simultaneous hit time and average
118  std::vector<float> dcmtime; ///<uncorrected times of hits
119  std::vector<float> dcmstime; ///<times with distance to readout, time-of-flight, timewalk corrections
120  std::vector<float> dcmrd; ///<distance to readout
121  std::vector<float> dcmpe; ///<number of photoelectrons
122  std::vector<float> dcmpl; ///<path length along track
123 
124 
125  //tree to store ntuple statistics
126  TTree* fOutTree;
127 
129 
130 
131 
132 
133 
134 
135  };
136 }
137 
138 namespace calib {
139 
140  //define speed of light in cm/ns
141  static const double sCmPerNsec = 29.9792458;
142 
144  // :
145  // Initialize member data here.
146  {
147  this->reconfigure(pset);
148 
149  produces< std::vector<caldp::DCMStat> >();
150  produces< art::Assns<caldp::DCMStat, rb::Track> >();
151  }
152 
153 //---------------------------------------------------------------------------
155  {
156  fTCTrackLabel = pset.get< std::string >("TCTrackLabel");
157  fTrackLabel = pset.get< std::string >("TrackLabel");
158  fMinDCMHits = 0;
159  fMinDCM = 0;
160  fDebug = pset.get< bool >("Debug");
161  fFiberSpeed = pset.get< float >("FiberSpeed");
162  fRemoveOutlier = pset.get< bool >("RemoveOutlier");
163  fRemoveThreshold = pset.get< int >("RemoveThreshold");
164  fCorrectTimewalk = pset.get< bool >("CorrectTimewalk");
165  fAdjustSpeed = pset.get< bool >("AdjustSpeed");
166  fPSetND = pset.get< fhicl::ParameterSet >("nd");
167  fPSetFD = pset.get< fhicl::ParameterSet >("fd");
168  fPSetTB = pset.get< fhicl::ParameterSet >("tb");
169 
170  }
171 
172  //---------------------------------------------------------------------------
174  {
175  }
176 
177  //---------------------------------------------------------------------------
179  {
180 
182 
183  // define the tree
184  fOutTree = tfs->make<TTree>("TimingCalib","Timing Calibaration");
185 
186  //define ntuple branches
187  fOutTree->Branch("run", &run, "run/I");
188  fOutTree->Branch("subrun", &subrun, "subrun/I");
189  fOutTree->Branch("event", &event, "event/I");
190  fOutTree->Branch("cosmictr", &cosmictr, "cosmictr/I");
191  fOutTree->Branch("sX", &sX, "sX/D");
192  fOutTree->Branch("sY", &sY, "sY/D");
193  fOutTree->Branch("sZ", &sZ, "sZ/D");
194  fOutTree->Branch("eX", &eX, "eX/D");
195  fOutTree->Branch("eY", &eY, "eY/D");
196  fOutTree->Branch("eZ", &eZ, "eZ/D");
197  fOutTree->Branch("plength", &plength, "plength/D");
198  fOutTree->Branch("minPl", &minPl, "minPl/I");
199  fOutTree->Branch("maxPl", &maxPl, "maxPl/I");
200  fOutTree->Branch("minPlX", &minPlX, "minPlX/I");
201  fOutTree->Branch("minPlY", &minPlY, "minPlY/I");
202  fOutTree->Branch("maxPlX", &maxPlX, "maxPlX/I");
203  fOutTree->Branch("maxPlY", &maxPlY, "maxPlY/I");
204  fOutTree->Branch("ncellsfitT", &ncellsfitT, "ncellsfitT/I");
205  fOutTree->Branch("ncellsfitX", &ncellsfitX, "ncellsfitX/I");
206  fOutTree->Branch("ncellsfitY", &ncellsfitY, "ncellsfitY/I");
207  fOutTree->Branch("evtTime", &evtTime, "evtTime/i");
208 
209  fOutTree->Branch("dcmnum", &dcmnum);
210  fOutTree->Branch("dcmhits", &dcmhits);
211  fOutTree->Branch("dcmsimtime", &dcmsimtime);
212  if (fDebug){
213  fOutTree->Branch("planenum", &planenum);
214  fOutTree->Branch("cellnum", &cellnum);
215  fOutTree->Branch("dcmsimresid", &dcmsimresid);
216  fOutTree->Branch("dcmtime", &dcmtime);
217  fOutTree->Branch("dcmstime", &dcmstime);
218  fOutTree->Branch("dcmrd", &dcmrd);
219  fOutTree->Branch("dcmpe", &dcmpe);
220  fOutTree->Branch("dcmpl",&dcmpl);
221  }
222 
223 
224  }
225 
226  //----------------------------------------------------------------------
228  {
230 
231  fhicl::ParameterSet pset;
232  switch(geom->DetId()){
234  pset = fPSetND;
235  break;
237  pset = fPSetFD;
238  break;
240  pset = fPSetTB;
241  break;
242  default:
243  assert(0 && "Unknown detector");
244  }
245 
246  fMinDCMHits = pset.get< int >("MinDCMHits");
247  fMinDCM = pset.get< int >("MinDCM");
248 
249  }
250 
251  //---------------------------------------------------------------------------
253  {
256 
257  std::unique_ptr< std::vector<caldp::DCMStat> > dcmstats(new std::vector<caldp::DCMStat>);
258 
259  std::unique_ptr<art::Assns<caldp::DCMStat, rb::Track> > dcmAssns(new art::Assns<caldp::DCMStat, rb::Track>);
260 
263  novadaq::cnv::DetId detId = geom->DetId();
264 
265  //get the cosmic tracks associated with each event
267  evt.getByLabel(fTrackLabel,trackcol);
268 
269  //get associationgs between slices and tracks
270  art::FindManyP<caldp::TCTrack> fmtc(trackcol, evt, fTCTrackLabel);
271 
272  std::vector< art::Ptr<rb::Track> > tracks;
273  art::fill_ptr_vector(tracks, trackcol);
274 
275  //store event info in ntuple
276  run = evt.run();
277  subrun = evt.subRun();
278  event = evt.id().event();
279 
280  evtTime = evt.time().timeHigh();
281 
282  //now loop over the tracks
283  for (unsigned int i=0; i<tracks.size(); ++i){
284 
285  this->clearNTuple();
286 
287  const art::Ptr<rb::Track> tr = tracks.at(i);
288 
289  //if ((tr->NXCell()==0) || (tr->NYCell()==0)) continue;
290  //get the tctrack associated with the cosmic track
291  std::vector<art::Ptr<caldp::TCTrack> > tct = fmtc.at(i);
292  //if we failed to get tctracks skip s
293  if (!(fmtc.isValid())) continue;
294  //there should only be one tctrack associated with a reco track
295  if (tct.size() != 1) continue;
296 
297  //make selection cuts on tracks
298  fDCMHitMap.clear();
299 
300 
301  cosmictr = i; //store track number
302  //passed initial round of cuts, now get track information
303  TVector3 start(tr->Start());
304  TVector3 end(tr->Stop());
305  sX = start.X();
306  sY = start.Y();
307  sZ = start.Z();
308  eX = end.X();
309  eY = end.Y();
310  eZ = end.Z();
311  plength = (end - start).Mag();
312  minPl = tct[0]->MinPlane(geo::kXorY);
313  maxPl = tct[0]->MaxPlane(geo::kXorY);
314 
315  minPlX = tct[0]->MinPlane(geo::kX);
316  minPlY = tct[0]->MinPlane(geo::kY);
317  maxPlX = tct[0]->MaxPlane(geo::kX);
318  maxPlY = tct[0]->MaxPlane(geo::kY);
319 
320  //make other selection cuts if desired
321 
322 
323  //put hits into map connected to DCM
324  for(unsigned int j=0; j<tct[0]->NPCHit(); ++j){
325  art::Ptr<caldp::PCHit> phit = tct[0]->PCHit(j);
326  int key = cmap->getNumberOfDCM(daqchannelmap::AFULLAFULL_DIBLOCK)*(phit->Diblock()-1) + phit->DCM();
327  fDCMHitMap[key].push_back(phit);
328  }
329 
330  //for each dcm calculate timing offsets
331  caldp::DCMStat dcms; //create object to hold dcm info for the track
332  for(std::map<int, std::vector<art::Ptr<caldp::PCHit> > >::iterator itr = fDCMHitMap.begin(); itr!=fDCMHitMap.end(); ++itr){
333  std::vector<art::Ptr<caldp::PCHit> > dcmcells = itr->second;
334  if ((int)dcmcells.size() < fMinDCMHits) continue;
335 
336  caldp::DCMTime dcmTime;
337 
338  std::vector<float> cellnumtemp;
339  std::vector<float> planenumtemp;
340  std::vector<float> dcmtimetemp;
341  std::vector<float> dcmstimetemp;
342  std::vector<float> dcmrdtemp;
343  std::vector<float> dcmpetemp;
344  std::vector<float> dcmpltemp;
345  std::vector<double> weights;
346  std::vector<double> sigma;
347  for(unsigned int k=0; k<dcmcells.size(); ++k){
348  double tns = dcmcells[k]->TNS();
349  double rd = dcmcells[k]->ReadDist();
350  double pl = dcmcells[k]->FlightLen();
351  if (fCorrectTimewalk) tns += tct[0]->ReadoutDistCorrection(dcmcells[k]->ReadDist(), detId);
352  //reverse any timing calibration that has been baked in
353  tns -= cal->GetTimingOffset(dcmcells[k]->Plane(), dcmcells[k]->Cell(), !bt->HaveTruthInfo());
354  float tmpFiberSpeed = fFiberSpeed;
355  if (fAdjustSpeed) tmpFiberSpeed = tct[0]->CalcFiberVelocity(rd, detId);
356  dcmstimetemp.push_back(tns - pl/sCmPerNsec - rd/tmpFiberSpeed);
357  dcmtimetemp.push_back(tns);
358  dcmrdtemp.push_back(rd);
359  dcmpetemp.push_back(dcmcells[k]->PE());
360  dcmpltemp.push_back(pl);
361  cellnumtemp.push_back(dcmcells[k]->Cell());
362  planenumtemp.push_back(dcmcells[k]->Plane());
363  double stdev = tct[0]->TNSUncertainty(dcmcells[k]->PE(), dcmcells[k]->GoodTime(), detId);
364  sigma.push_back(stdev);
365  double w = 1.0/(stdev*stdev);
366  weights.push_back(w);
367  }
368 
369  while (true){
370  double st = 0.0;
371  double weight = 0.0;
372  for(unsigned int k=0; k<dcmstimetemp.size(); ++k){
373  st += dcmstimetemp[k]*weights[k];
374  weight += weights[k];
375  }
376  st /= weight;
377  double maxresid = 0;
378  unsigned int maxpos = 0;
379  for(unsigned int k=0; k<dcmstimetemp.size(); ++k){
380  float recresid = st - dcmstimetemp[k];
381  if (fabs(recresid/sigma[k]) > fabs(maxresid)){
382  maxresid = fabs(recresid/sigma[k]);
383  maxpos = k;
384  }
385  }
386  if (dcmcells.size() <= 2) break;
387  if (maxresid > 2.0){
388  if (fRemoveOutlier){
389  dcmcells.erase(dcmcells.begin()+maxpos);
390  cellnumtemp.erase(cellnumtemp.begin()+maxpos);
391  planenumtemp.erase(planenumtemp.begin()+maxpos);
392  dcmtimetemp.erase(dcmtimetemp.begin()+maxpos);
393  dcmstimetemp.erase(dcmstimetemp.begin()+maxpos);
394  dcmrdtemp.erase(dcmrdtemp.begin()+maxpos);
395  dcmpetemp.erase(dcmpetemp.begin()+maxpos);
396  dcmpltemp.erase(dcmpltemp.begin()+maxpos);
397  weights.erase(weights.begin()+maxpos);
398  sigma.erase(sigma.begin()+maxpos);
399  continue;
400  }
401  }
402  if ((int)dcmcells.size() < fMinDCMHits) break;
403  dcmnum.push_back(itr->first);
404  dcmhits.push_back(dcmstimetemp.size());
405  dcmsimtime.push_back(st);
406  dcmTime.dcmTime = st;
407  dcmTime.dcmNum = itr->first;
408  dcmTime.dcmHits = dcmcells.size();
409  double sum = 0.0;
410  for(unsigned int k=0; k<dcmstimetemp.size(); ++k){
411  double recresid = st - dcmstimetemp[k];
412  dcmsimresid.push_back(recresid);
413  sum += weights[k]*(dcmstimetemp[k] - st)*(dcmstimetemp[k]-st);
414  cellnum.push_back(cellnumtemp[k]);
415  planenum.push_back(planenumtemp[k]);
416  dcmtime.push_back(dcmtimetemp[k]);
417  dcmstime.push_back(dcmstimetemp[k]);
418  dcmrd.push_back(dcmrdtemp[k]);
419  dcmpe.push_back(dcmpetemp[k]);
420  dcmpl.push_back(dcmpltemp[k]);
421  }
422  dcmTime.sigma = sqrt(sum/weight);
423  dcms.dcmList.push_back(dcmTime);
424  break;
425  }
426 
427  }
428 
429  if ((int)dcms.dcmList.size() < fMinDCM) continue;
430 
431  dcmstats->push_back(dcms);
432 
433  util::CreateAssn(*this, evt, *dcmstats, tr, *dcmAssns);
434 
435 
436 
437  ncellsfitT = tct[0]->NPCHit();
438  ncellsfitX = tct[0]->NXPCHit();
439  ncellsfitY = tct[0]->NYPCHit();
440 
441 
442  fOutTree->Fill();
443 
444 
445  }//end loop over tracks
446 
447  evt.put(std::move(dcmstats));
448  evt.put(std::move(dcmAssns));
449 
450  }
451 
452  //---------------------------------------------------------------------------
454 
455  //track start coordinate
456  sX = -999.0;
457  sY = -999.0;
458  sZ = -999.0;
459  //track end coordinate
460  eX = -999.0;
461  eY = -999.0;
462  eZ = -999.0;
463  plength = -5.0;
464  //min, max plane
465  minPl = -5;
466  maxPl = -5;
467  //min, max cell
468  minPlX = -5;
469  minPlY = -5;
470  maxPlX = -5;
471  maxPlY = -5;
472  //number of cells used in final fiber speed fit
473  ncellsfitT = 0;
474  ncellsfitX = 0;
475  ncellsfitY = 0;
476 
477 
478  dcmnum.clear();
479  dcmhits.clear();
480  dcmsimtime.clear();
481  planenum.clear();
482  cellnum.clear();
483  dcmsimresid.clear();
484  dcmtime.clear();
485  dcmstime.clear();
486  dcmrd.clear();
487  dcmpe.clear();
488  dcmpl.clear();
489 
490  }
491 
492 } //end namespace
493 
494 //-----------------------------------------------------------------------------
495 
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
def stdev(lst)
Definition: HandyFuncs.py:286
int fMinDCMHits
minimum number of hits on a dcm to keep results
std::vector< float > dcmrd
distance to readout
virtual unsigned int getNumberOfDCM(DiBlock_TYPE dbt) const =0
How many DCMs does this diblock have?
std::vector< float > planenum
plane number
float sigma
uncertainty in average hit time
Definition: DCMStat.h:21
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int fMinDCM
minimum number of good dcms on a track
const Var weight
X or Y views.
Definition: PlaneGeo.h:30
bool fDebug
Add exta debugging info to ntuple.
TimingCalibration(fhicl::ParameterSet const &p)
std::vector< float > dcmtime
uncorrected times of hits
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
std::vector< float > cellnum
cell number
Vertical planes which measure X.
Definition: PlaneGeo.h:28
bool fCorrectTimewalk
adjust hit times based on timewalk effect for distance to readout
void reconfigure(const fhicl::ParameterSet &pset)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void beginRun(art::Run &run) override
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::map< int, std::vector< art::Ptr< caldp::PCHit > > > fDCMHitMap
Definition: Run.h:31
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
std::string fTCTrackLabel
Where to find calibration tracks.
std::string fTrackLabel
Where to find reco Tracks.
int dcmHits
number of hits in dcm
Definition: DCMStat.h:19
Far Detector at Ash River, MN.
CDPStorage service.
std::vector< float > dcmsimtime
simultaneous hit time in dcm
static DAQChannelMap * getInstance(int detID)
int fRemoveThreshold
number of sigma away from mean to remove hits
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
float dcmTime
weighted average time of hits in dcm, correcting for time-of-flight and path lenght ...
Definition: DCMStat.h:20
std::vector< float > dcmstime
times with distance to readout, time-of-flight, timewalk corrections
std::vector< float > dcmnum
dcm index
Var weights
const double j
Definition: BetheBloch.cxx:29
std::vector< float > dcmpl
path length along track
Identifier for diblocks using a 32/32 configuration.
int Diblock() const
Return diblock value.
Definition: PCHit.h:28
double sigma(TH1F *hist, double percentile)
Definition: run.py:1
EventNumber_t event() const
Definition: EventID.h:116
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
static const double sCmPerNsec
bool fRemoveOutlier
remove outlier hits in DCMs based on timing resolution
double GetTimingOffset(unsigned int const &plane, unsigned int const &cell, bool const &isData)
Get the timing offset for a given plane, cell. Useful downstream to check calibration.
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
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
float Mag() const
std::vector< float > dcmsimresid
difference between simultaneous hit time and average
assert(nhit_max >=nhit_nbins)
std::vector< float > dcmhits
number of hits in dcm
bool fAdjustSpeed
adjust fiberspeed based on distance to readout instead of assuming flat speed
Timestamp time() const
Definition: Event.h:61
cmap::CMap class source code
Definition: CMap.cxx:17
std::vector< DCMTime > dcmList
Definition: DCMStat.h:31
int DCM() const
Return dcm value.
Definition: PCHit.h:30
float fFiberSpeed
fixed fiberspeed
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
void produce(art::Event &e) override
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
std::vector< float > dcmpe
number of photoelectrons
Float_t w
Definition: plot.C:20
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
int dcmNum
Number of the dcm.
Definition: DCMStat.h:18