FiberCalibration_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: FiberCalibration
3 // Module Type: producer
4 // File: FiberCalibration_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 "Geometry/LiveGeometry.h"
13 #include "RecoBase/CellHit.h"
14 #include "RecoBase/Cluster.h"
15 #include "RecoBase/Track.h"
16 #include "Utilities/func/MathUtil.h"
19 #include "NovaDAQConventions/DAQConventions.h"
20 
21 //framework includes
29 #include "Utilities/AssociationUtil.h"
30 
31 //ROOT includes
32 #include "TTree.h"
33 #include "TH2F.h"
34 #include "TH1F.h"
35 
36 namespace calib {
37 
39  public:
40  explicit FiberCalibration(fhicl::ParameterSet const & p);
41  virtual ~FiberCalibration();
42 
43  void produce(art::Event & e) override;
44 
45  void reconfigure(const fhicl::ParameterSet& pset);
46  void beginJob() override;
47  void beginRun(art::Run& run) override;
48  void clearNTuple();
49 
50 
51  private:
52 
53  std::string fPCHitLabel; ///< Where to find calibration hits
54  std::string fQualXYLabel; ///< Instance label for hits with xz quality
55  std::string fQualZLabel; ///< Instance label for hits with z quality
56  std::string fQualAvgLabel; ///< Instance label for hits with avg quality
57  std::string fTrackLabel; ///< Where to find Tracks
58 
59  //variables to select track quality
60  double fFitFrac; ///<fraction of cells on track kept in final fiber speed fit
61  int fFitCells; ///<number of cells used in final fiber speed fit for a track
62  double fAsy; ///<maximum asymmetry between the x and y cells used in the final fit
63  int fTrackPlanes; ///<number of planes track passes through
64  double fFiducialX; ///<maximum distance in cm muon endpoints can be from detector edge in X
65  double fFiducialY; ///<maximum distance in cm muon endpoints can be from detector edge in Y
66  double fFiducialZ; ///<maximum distance in cm muon endpoints can be from detector edge in Z
67  int fStartPlaneDiff; ///<maximum difference between minimum plane in each view associated with the track
68  int fEndPlaneDiff; ///<maximum difference between maximum plane in each view associated with the track
69  double fMinZ; /// hack for minimum z detector boundary for situation were front diblocks are missing
70  double fMaxZ; /// hack for maximum z in detector, which could be incorrectly calculated depending on missing diblocks
71 
72  double fTimeRes; ///<number of sigma timing resolution used to keep/reject hits in speed fit
73  bool fUseLiveGeo; ///<use the live geometry to find beginning and ending diblock
74 
75  //ntuple variables
76  //event statistics
77  int run;
78  int subrun;
79  int event;
80  int cosmictr;
81  //track start coordinate
82  float sX;
83  float sY;
84  float sZ;
85  //track end coordinate
86  float eX;
87  float eY;
88  float eZ;
89  //pathlength
90  float plength;
91  //min, max plane
92  int minPl;
93  int maxPl;
94  int minPlX;
95  int minPlY;
96  int maxPlX;
97  int maxPlY;
98  //min, max cell
99  int minCX;
100  int minCY;
101  int maxCX;
102  int maxCY;
103  //number of cells
104  int ncellsT;
105  int ncellsX;
106  int ncellsY;
107  //number of cells with good timing and calibrated
111  //number of cells used in final fiber speed fit
115  //fiber speed calculated from fit
116  float fiberspeed;
117  // 1/beta calculated from fit of muon path using the fiber speed
118  float beta;
119 
120  //tree to store ntuple statistics
121  TTree* fOutTree;
122 
124 
125 
126 
127 
128 
129 
130 
131  };
132 }
133 
134 namespace calib {
135 
136  //define speed of light in cm/ns
137  //static const double sCmPerNsec = 29.9792458;
138 
140  // Initialize member data here.
141  {
142  this->reconfigure(pset);
143 
144  produces< std::vector<caldp::TCTrack> >();
145  produces< art::Assns<caldp::TCTrack, rb::Track> >();
146  }
147 
148  //---------------------------------------------------------------------------
150  {
151  fPCHitLabel = pset.get< std::string >("PCHitLabel");
152  fTrackLabel = pset.get< std::string >("TrackLabel");
153  fQualXYLabel = pset.get< std::string >("QualXYLabel");
154  fQualZLabel = pset.get< std::string >("QualZLabel");
155  fQualAvgLabel = pset.get< std::string >("QualAvgLabel");
156  fFitFrac = 0;
157  fFitCells = 0;
158  fAsy = 0;
159  fTrackPlanes = 0;
160  fFiducialX = 0;
161  fFiducialY = 0;
162  fFiducialZ = 0;
163  fStartPlaneDiff = 0;
164  fEndPlaneDiff = 0;
165  fTimeRes = 0;
166  fMinZ = 0;
167  fMaxZ = 0;
168  fUseLiveGeo = false;
169  fPSetND = pset.get<fhicl::ParameterSet>("nd");
170  fPSetFD = pset.get<fhicl::ParameterSet>("fd");
171  fPSetTB = pset.get<fhicl::ParameterSet>("tb");
172  }
173 
174  //---------------------------------------------------------------------------
176  {
177  // Clean up dynamic memory and other resources here.
178  }
179 
180  //---------------------------------------------------------------------------
182  {
183  TH1::SetDefaultSumw2();
184  TH2::SetDefaultSumw2();
185 
187 
188  // define the tree
189  fOutTree = tfs->make<TTree>("FiberCalib","Fiber Calibaration");
190 
191  //define ntuple branches
192  fOutTree->Branch("run", &run, "run/I");
193  fOutTree->Branch("subrun", &subrun, "subrun/I");
194  fOutTree->Branch("event", &event, "event/I");
195  fOutTree->Branch("cosmictr", &cosmictr, "cosmictr/I");
196  fOutTree->Branch("sX", &sX, "sX/F");
197  fOutTree->Branch("sY", &sY, "sY/F");
198  fOutTree->Branch("sZ", &sZ, "sZ/F");
199  fOutTree->Branch("eX", &eX, "eX/F");
200  fOutTree->Branch("eY", &eY, "eY/F");
201  fOutTree->Branch("eZ", &eZ, "eZ/F");
202  fOutTree->Branch("plength", &plength, "plength/F");
203  fOutTree->Branch("minPl", &minPl, "minPl/I");
204  fOutTree->Branch("maxPl", &maxPl, "maxPl/I");
205  fOutTree->Branch("minPlX", &minPlX, "minPlX/I");
206  fOutTree->Branch("maxPlX", &maxPlX, "maxPlX/I");
207  fOutTree->Branch("minPlY", &minPlY, "minPlY/I");
208  fOutTree->Branch("maxPlY", &maxPlY, "maxPlY/I");
209  //fOutTree->Branch("minCX", &minCX, "minCX/I");
210  //fOutTree->Branch("minCY", &minCY, "minCY/I");
211  //fOutTree->Branch("maxCX", &maxCX, "maxCX/I");
212  //fOutTree->Branch("maxCY", &maxCY, "maxCY/I");
213  //fOutTree->Branch("ncellsT", &ncellsT, "ncellsT/I");
214  //fOutTree->Branch("ncellsX", &ncellsX, "ncellsX/I");
215  //fOutTree->Branch("ncellsY", &ncellsY, "ncellsY/I");
216  fOutTree->Branch("ncellscalT", &ncellscalT, "ncellscalT/I");
217  fOutTree->Branch("ncellscalX", &ncellscalX, "ncellscalX/I");
218  fOutTree->Branch("ncellscalY", &ncellscalY, "ncellscalY/I");
219  fOutTree->Branch("ncellsfitT", &ncellsfitT, "ncellsfitT/I");
220  fOutTree->Branch("ncellsfitX", &ncellsfitX, "ncellsfitX/I");
221  fOutTree->Branch("ncellsfitY", &ncellsfitY, "ncellsfitY/I");
222  fOutTree->Branch("fiberspeed", &fiberspeed, "fiberspeed/F");
223  fOutTree->Branch("beta", &beta, "beta/F");
224 
225  }
226 
227  //----------------------------------------------------------------------
229  {
231 
232  fhicl::ParameterSet pset;
233  switch(geom->DetId()){
235  pset = fPSetND;
236  break;
238  pset = fPSetFD;
239  break;
241  pset = fPSetTB;
242  break;
243  default:
244  assert(0 && "Unknown detector");
245  }
246 
247  fFitFrac = pset.get< double >("FitFrac");
248  fFitCells = pset.get< int >("FitCells");
249  fAsy = pset.get< double >("Asy");
250  fTrackPlanes = pset.get< int >("TrackPlanes");
251  fFiducialX = pset.get< double >("FiducialX");
252  fFiducialY = pset.get< double >("FiducialY");
253  fFiducialZ = pset.get< double >("FiducialZ");
254  fStartPlaneDiff = pset.get< int >("StartPlaneDiff");
255  fEndPlaneDiff = pset.get< int >("EndPlaneDiff");
256  fTimeRes = pset.get< double >("TimeRes");
257  fMinZ = pset.get< double >("MinZ");
258  fMaxZ = pset.get< double >("MaxZ");
259  fUseLiveGeo = pset.get< bool >("UseLiveGeo");
260  }
261 
262  //---------------------------------------------------------------------------
264  {
265 
266  std::unique_ptr< std::vector<caldp::TCTrack> > tctracks(new std::vector<caldp::TCTrack>);
267 
268  std::unique_ptr<art::Assns<caldp::TCTrack, rb::Track> > tcAssns(new art::Assns<caldp::TCTrack, rb::Track>);
269 
272  novadaq::cnv::DetId detId = geom->DetId();
273 
274  //get the cosmic tracks associated with each event
276  evt.getByLabel(fTrackLabel,trackcol);
277 
278  //get associationgs between cosmic tracks and pchits
282 
283  std::vector< art::Ptr<rb::Track> > tracks;
284  art::fill_ptr_vector(tracks, trackcol);
285 
286  //store event info in ntuple
287  run = evt.run();
288  subrun = evt.subRun();
289  event = evt.id().event();
290 
291  //if we are in data use the live geometry to find front and back of the volume.
292  if((geom->DetId() == novadaq::cnv::kFARDET) && (fUseLiveGeo)){
293  fMinZ = livegeom->InstrumentedDetFront();
294  fMaxZ = livegeom->InstrumentedDetBack();
295  }
296 
297  //now loop over the tracks
298  for (unsigned int i=0; i<tracks.size(); ++i){
299 
300  this->clearNTuple();
301 
302  const art::Ptr<rb::Track> tr = tracks.at(i);
303  //if ((tr->NXCell()==0) || (tr->NYCell()==0)) continue;
304  //get the pc hits associated with the cosmic track
305  std::vector<art::Ptr<caldp::PCHit> > pcxy = fmpcxy.at(i);
306  std::vector<art::Ptr<caldp::PCHit> > pcz = fmpcz.at(i);
307  std::vector<art::Ptr<caldp::PCHit> > pcavg = fmpcavg.at(i);
308  //empty pchit container to combine hits from the lists in
309  std::vector<art::Ptr<caldp::PCHit> > pch;
310  //if we failed to get pc hits skip slice
311  if (fmpcxy.isValid()){
312  for (unsigned int j=0; j<pcxy.size(); ++j) pch.push_back(pcxy[j]);
313  }
314  if (fmpcz.isValid()){
315  for (unsigned int j=0; j<pcz.size(); ++j) pch.push_back(pcz[j]);
316  }
317  if (fmpcavg.isValid()){
318  for (unsigned int j=0; j<pcavg.size(); ++j) pch.push_back(pcavg[j]);
319  }
320  //must be at least two pchits with this track
321  if (pch.size() < 2) continue;
322 
323 
324  cosmictr = i; //store track number
325  //passed initial round of cuts, now get track information
326  TVector3 start(tr->Start());
327  TVector3 end(tr->Stop());
328  sX = start.X();
329  sY = start.Y();
330  sZ = start.Z();
331  eX = end.X();
332  eY = end.Y();
333  eZ = end.Z();
334  plength = (end - start).Mag();
335  /*minPl = tr->MinPlane();
336  maxPl = tr->MaxPlane();
337  minCX = tr->MinCell(geo::kX);
338  minCY = tr->MinCell(geo::kY);
339  maxCX = tr->MaxCell(geo::kX);
340  maxCY = tr->MaxCell(geo::kY);
341  ncellsT = tr->NCell();
342  ncellsX = tr->NXCell();
343  ncellsY = tr->NYCell();*/
344 
345  //track quality cuts
346  //muon starpoint must be near either the ceiling or one of the walls of the detector
347  //cut assumes muons cannot enter from bottom of detector
348  //this selects muons that are not stoppers and pass all the way through the detector
349  if ((sX < geom->DetHalfWidth() - fFiducialX) &&
350  (sX > -(geom->DetHalfWidth() - fFiducialX)) &&
351  (sY < geom->DetHalfHeight() - fFiducialY) &&
352  (sZ < fMaxZ - fFiducialZ) &&
353  (sZ > std::max(0.0, fMinZ) + fFiducialZ) ) continue;
354  if ((eX < geom->DetHalfWidth() - fFiducialX) &&
355  (eX > -(geom->DetHalfWidth() - fFiducialX)) &&
356  (eY > -(geom->DetHalfHeight() - fFiducialY)) &&
357  (eZ < fMaxZ - fFiducialZ) &&
358  (eZ > std::max(0.0, fMinZ) + fFiducialZ) ) continue;
359 
360  caldp::TCTrack tc;
361 
362  tc.SetStart(start);
363  tc.SetStop(end);
365 
366  //now loop over cell hits to get the necessary ones for fiber calibration
367  for(unsigned int j=0; j<pch.size(); ++j){
368  //if the hit had a bad timing fit, throw it out
369  //if (!(pch[j]->GoodTime())) continue;
370 
371  //if the distance to readout is less then 0, don't use
372  if (pch[j]->ReadDist() < 0.0) continue;
373 
374  tc.Add(pch[j]);
375 
376  }//end loop over cells
377 
378  ncellscalT = tc.NPCHit();
379  ncellscalX = tc.NXPCHit();
380  ncellscalY = tc.NYPCHit();
381 
382  minPlX = tc.MinPlane(geo::kX);
383  minPlY = tc.MinPlane(geo::kY);
384  maxPlX = tc.MaxPlane(geo::kX);
385  maxPlY = tc.MaxPlane(geo::kY);
386  minPl = tc.MinPlane(geo::kXorY);
387  maxPl = tc.MaxPlane(geo::kXorY);
388 
389  //track must have traveled through a minimum number of planes
390  if ((maxPl-minPl) < fTrackPlanes) continue;
391 
392  //Require the minimum and maximum planes in each view to agree
393  //Cut in place to reduce case of poor reconstruction or view matching
394  //In regions with lots of dead channels this cut would need to be reconsidered
395  if (abs(minPlX-minPlY) > fStartPlaneDiff ) continue;
396  if (abs(maxPlX-maxPlY) > fEndPlaneDiff ) continue;
397 
398  //perform the fit for fiber speed for the track, do NOT drop hits outside specified residual window
399  //NOTE: When DCMs are in sync this can be configured to drop outlying hits and calculate the fiber speed. This is not desired when calculating dcm offsets
400  std::vector<double> residuals;
401  double ficept = 0.0;
402  double chisqr = 0.0;
403  fiberspeed = tc.CalculateFiberVelocity(fTimeRes,&residuals, &ficept, &chisqr, detId, false,true);
404 
405  ncellsfitT = residuals.size();
406 
407  //more quality cuts on the fit
408  //make sure cells used in fit are above threshold
409  if (ncellscalT < fFitCells) continue;
410  //make sure cells used in the fit represent a certain fraction of total hits
411  if (((double)ncellsfitT/(double)ncellscalT) < fFitFrac) continue;
412 
413  ncellsfitX = tc.NXPCHit();
414  ncellsfitY = tc.NYPCHit();
415 
416  //cut to check difference between number of x and y cells used in the fit
417  double asy = ((double)(ncellsfitY-ncellsfitX))/((double)(ncellsfitY+ncellsfitX));
418  if (fabs(asy) > fAsy) continue;
419 
420  double dummy=0.0;
421  residuals.clear();
422  beta = tc.CalculateMuonVelocity(&residuals, &dummy, &chisqr, detId, tc.FiberVelocity());
423 
424  tctracks->push_back(tc);
425 
426  util::CreateAssn(*this, evt, *tctracks, tr, *tcAssns);
427 
428  fOutTree->Fill();
429 
430 
431  }//end loop over slices
432 
433  evt.put(std::move(tctracks));
434  evt.put(std::move(tcAssns));
435  }
436 
437  //---------------------------------------------------------------------------
439 
440  //track start coordinate
441  sX = -999.0;
442  sY = -999.0;
443  sZ = -999.0;
444  //track end coordinate
445  eX = -999.0;
446  eY = -999.0;
447  eZ = -999.0;
448  //pathlength variables
449  plength = -5.0;
450  //min, max plane
451  minPl = -5;
452  maxPl = -5;
453  //min, max cell
454  minCX = -5;
455  minCY = -5;
456  maxCX = -5;
457  maxCY = -5;
458  //number of cells
459  ncellsT = 0;
460  ncellsX = 0;
461  ncellsY = 0;
462  //number of cells with good timing and calibrated
463  ncellscalT = 0;
464  ncellscalX = 0;
465  ncellscalY = 0;
466  //number of cells used in final fiber speed fit
467  ncellsfitT = 0;
468  ncellsfitX = 0;
469  ncellsfitY = 0;
470  //fit parameters
471  fiberspeed = -99.0;
472  beta = -5.0;
473 
474 
475  }
476 
477 } //end namespace
478 
479 //-----------------------------------------------------------------------------
480 
T max(const caf::Proxy< T > &a, T b)
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
SubRunNumber_t subRun() const
Definition: Event.h:72
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.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::string fTrackLabel
Where to find Tracks.
virtual void SetStop(TVector3 stop)
set track stop point
Definition: TCTrack.h:94
X or Y views.
Definition: PlaneGeo.h:30
int fStartPlaneDiff
maximum difference between minimum plane in each view associated with the track
const char * p
Definition: xmltok.h:285
double fFiducialX
maximum distance in cm muon endpoints can be from detector edge in X
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::string fPCHitLabel
Where to find calibration hits.
FiberCalibration(fhicl::ParameterSet const &p)
void abs(TH1 *hist)
virtual void SetStart(TVector3 start)
set track start point
Definition: TCTrack.h:92
int fEndPlaneDiff
maximum difference between maximum plane in each view associated with the track
double fFiducialZ
maximum distance in cm muon endpoints can be from detector edge in Z
bool fUseLiveGeo
use the live geometry to find beginning and ending diblock
int fFitCells
number of cells used in final fiber speed fit for a track
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
double FiberVelocity() const
return the velocity of the optical fiber, measured for this track
Definition: TCTrack.h:99
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
Definition: Run.h:31
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
std::string fQualXYLabel
Instance label for hits with xz quality.
Far Detector at Ash River, MN.
CDPStorage service.
int fTrackPlanes
number of planes track passes through
virtual void Add(const art::Ptr< caldp::PCHit > &pchit)
add a pchit to the track
Definition: TCTrack.cxx:25
double fAsy
maximum asymmetry between the x and y cells used in the final fit
double fFitFrac
fraction of cells on track kept in final fiber speed fit
double InstrumentedDetFront()
get instrumented detector front of downstream part
virtual void SetTotalLength(double length)
set track total length
Definition: TCTrack.h:96
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
int evt
Near Detector in the NuMI cavern.
unsigned int NPCHit(geo::View_t view) const
return the number of hits in a view
Definition: TCTrack.cxx:43
void beginRun(art::Run &run) override
double InstrumentedDetBack()
get instrumented detector back of downstream part
const double j
Definition: BetheBloch.cxx:29
double DetHalfHeight() const
Definition: run.py:1
unsigned int NXPCHit() const
return number of hits in x view
Definition: TCTrack.h:37
void reconfigure(const fhicl::ParameterSet &pset)
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
double DetHalfWidth() const
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
std::string fQualZLabel
Instance label for hits with z quality.
assert(nhit_max >=nhit_nbins)
double fTimeRes
hack for maximum z in detector, which could be incorrectly calculated depending on missing diblocks ...
double CalculateMuonVelocity(std::vector< double > *resids, double *icept, double *chisqr, novadaq::cnv::DetId detectorID, double fiberVel=9999.0) const
Definition: TCTrack.cxx:371
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
std::string fQualAvgLabel
Instance label for hits with avg quality.
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
double fFiducialY
maximum distance in cm muon endpoints can be from detector edge in Y
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
double fMaxZ
hack for minimum z detector boundary for situation were front diblocks are missing ...
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
void produce(art::Event &e) override