RobustHitFinder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: RobustHitFinder
3 // Module Type: producer
4 // File: RobustHitFinder_module.cc
5 //
6 // Generated at Thu Aug 18 04:21:10 2016 by Matthew Thiesse using artmod
7 // from cetpkgsupport v1_10_02.
8 ////////////////////////////////////////////////////////////////////////
9 
10 /*************************************************
11 
12  October 2016
13 
14  m.thiesse@sheffield.ac.uk
15 
16  **************************************************/
17 
20 #include "canvas/Persistency/Common/FindOneP.h"
21 #include "canvas/Persistency/Common/FindManyP.h"
27 #include "art_root_io/TFileService.h"
29 #include "fhiclcpp/ParameterSet.h"
31 
32 #include "nurandom/RandomUtils/NuRandomService.h"
42 
50 
57 
58 #include <math.h>
59 #include <memory>
60 #include <algorithm>
61 #include <fstream>
62 #include <iostream>
63 #include <vector>
64 #include <map>
65 
66 #include "TMath.h"
67 #include "TVector3.h"
68 #include "TTree.h"
69 #include "TF1.h"
70 #include "TH1.h"
71 
72 #include "CLHEP/Random/RandomEngine.h"
73 
74 //for the counter position map
76 
77 #include "RobustHitFinderSupport.h"
78 #include "HitLineFitAlg.h"
79 #include "RMSHitFinderAlg.h"
80 
81 namespace dune {
82  class RobustHitFinder;
83 
84  // Struct to contain combined IDE information
85  struct MCHit {
87  int pdg;
88  std::vector<sim::IDE> idevec;
89  std::vector<unsigned short> tdcvec;
90  std::vector<double> chargevec;
91  std::vector<const simb::MCParticle *> particlevec;
92  };
93 }
94 
96 public:
97  explicit RobustHitFinder(fhicl::ParameterSet const & p);
98  // The destructor generated by the compiler is fine for classes
99  // without bare pointers or other resource use.
100 
101  // Plugins should not be copied or assigned.
102  RobustHitFinder(RobustHitFinder const &) = delete;
103  RobustHitFinder(RobustHitFinder &&) = delete;
104  RobustHitFinder & operator = (RobustHitFinder const &) = delete;
105  RobustHitFinder & operator = (RobustHitFinder &&) = delete;
106 
107  // Required functions.
108  void produce(art::Event & e) override;
109 
110  // Selected optional functions.
111  void beginJob() override;
112  void reconfigure(fhicl::ParameterSet const & p) ;
113 
114 private:
115  void SetTreeVariables(const dune::ChannelInformation & chan, const dune::HitInformation & hit, const dune::HitLineFitAlg::HitLineFitResults & fitresult);
116  void SetTreeVariablesMCTruth(const ChanMap_t & chanMap, art::Event & e);
117  void MakeupMissedHits(ChanMap_t & chanMap, HitVec_t & hitVec);
118  void FillHitInformation(dune::ChannelInformation & chan, dune::HitVec_t & hitVec, bool assumedHit);
119  bool ValidTrigger(std::vector<unsigned int> evtTriggers, unsigned int & c1arg, unsigned int & c2arg, unsigned int & trignumarg);
120  float TimeToDriftDist(float thistime, unsigned int thistpc);
121  float TimeToDisplacement(float thistime);
122  float TimeToX(float thistime, unsigned int thistpc);
123  float hitGeomDist(TVector3 hitloc, TVector3 trigloc1, TVector3 trigloc2);
124  void Reset();
125 
128 
129  TTree * fAllTree;
130  //int run;
131  //int subrun;
132  //int event;
133  //unsigned int trignum;
134  //unsigned int c1;
135  //unsigned int c2;
136 
137  TTree * fGoodTree;
138  //int run;
139  //int subrun;
140  //int event;
141  //unsigned int trignum;
142  //unsigned int c1;
143  //unsigned int c2;
144 
145  TTree * fTree;
146  int run;
147  int subrun;
148  int event;
149  double t0;
150  unsigned int c1;
151  unsigned int c2;
152  unsigned int trignum;
153  float c1x;
154  float c1y;
155  float c1z;
156  float c2x;
157  float c2y;
158  float c2z;
159  float distancecut;
160  std::vector<int> allchannels;
161  std::vector<int> channel;
162  std::vector<int> wire;
163  std::vector<int> tpc;
164  std::vector<int> signalsize;
165  std::vector<std::vector<float> > signal;
166  std::vector<std::vector<float> > signalFilter;
167  std::vector<float> baseline;
168  std::vector<float> rms;
169  std::vector<float> baselineFilter;
170  std::vector<float> rmsFilter;
171  std::vector<float> pedmean;
172  std::vector<float> pedrms;
173  std::vector<float> integral;
174  std::vector<float> integralFilter;
175  std::vector<float> sumADC;
176  std::vector<float> sigmaintegral;
177  std::vector<float> sigmaintegralFilter;
178  std::vector<float> amplitude;
179  std::vector<float> amplitudeFilter;
180  std::vector<float> peaktick;
181  std::vector<float> peaktickFilter;
182  std::vector<float> peaktime;
183  std::vector<float> peaktimeFilter;
184  std::vector<int> begintick;
185  std::vector<int> endtick;
186  std::vector<int> width;
187  std::vector<float> hitx;
188  std::vector<float> hity;
189  std::vector<float> hitz;
190  std::vector<float> hiterrxlo;
191  std::vector<float> hiterrxhi;
192  std::vector<float> hiterrylo;
193  std::vector<float> hiterryhi;
194  std::vector<float> hiterrzlo;
195  std::vector<float> hiterrzhi;
196  std::vector<float> perpdist;
197  std::vector<float> hitt;
198  std::vector<float> driftdist;
199  std::vector<bool> countercut;
200  float fitconstant;
202  float fitlinear;
206  float fitchi2;
208  float fitndf;
209  float fitmle;
211  std::vector<bool> fitrealhit;
212  std::vector<float> segmentlength;
213  std::vector<bool> assumedhit;
214  std::vector<int> numGoodHitsChan;
215  std::vector<bool> ismctruth;
224 
225  bool fMakeTree;
229  std::vector<float> fEfield;
235  float fMinPedMean;
236  float fMaxPedMean;
237  float fMinPedRms;
238  float fMaxPedRms;
244  float fMCScale;
245 
250  float fC1Vert;
251  float fC1Horiz;
252  float fC2Vert;
253  float fC2Horiz;
254  //UInt_t fSeedValue; // unused
255  std::map<unsigned int, std::pair<TVector3, std::vector<TVector3> > > fCounterPositionMap;
256 
259 
264  const lariov::DetPedestalProvider& fPedestalRetrievalAlg = *(lar::providerFrom<lariov::DetPedestalService>());
265 
266  //throwaway variables
267  float temppedmean;
268  float temppedrms;
269 
270  CLHEP::HepRandomEngine& fEngine;
271 };
272 
273 
275  : EDProducer{p},
276  fFitAlg(p.get<fhicl::ParameterSet>("HitLineFitAlg")),
277  fHitFinderAlg(p.get<fhicl::ParameterSet>("RMSHitFinderAlg")),
278  fClks(lar::providerFrom<detinfo::DetectorClocksService>()),
279  fDetProp(lar::providerFrom<detinfo::DetectorPropertiesService>()),
280  fEngine(art::ServiceHandle<rndm::NuRandomService>{}->createEngine(*this,"HepJamesRandom","Seed"))
281 {
282  this->reconfigure(p);
283 
285 }
286 
288 {
289  Reset();
290 
292 
293  fFitAlg.SetSeed(fEngine.getSeed());
294 
295  // get recob::Wires
297  if (!e.getByLabel(fWireModuleLabel,wireHandle))
298  {
299  mf::LogError("RobustHitFinder") << "No recob::Wires found when expected";
300  hcol.put_into(e);
301  return;
302  }
303 
304  // get associated raw digits
305  art::FindOneP<raw::RawDigit> rawdigits(wireHandle,e,fWireModuleLabel);
306 
307  // Get all T0 products
309  if (!e.getByLabel(fCounterT0ModuleLabel,t0Handle))
310  {
311  mf::LogError("RobustHitFinder") << "No T0s found in this event";
312  hcol.put_into(e);
313  return;
314  }
315 
316  // get associated external triggers
317  art::FindManyP<raw::ExternalTrigger> triggers(t0Handle,e,fCounterT0ModuleLabel);
318 
319  run = e.run();
320  subrun = e.subRun();
321  event = e.event();
322 
324  for (size_t iwire = 0; iwire < wireHandle->size(); ++iwire)
325  {
326  art::Ptr<recob::Wire> wire(wireHandle,iwire);
327  if (wire->View() != geo::kZ) continue;
328  if (fCSP->IsBad(wire->Channel())) continue;
329  int temptpc = fGeom->ChannelToWire(wire->Channel())[0].TPC;
330  if (temptpc==0) nwiresTPC0++;
331  else if (temptpc==1) nwiresTPC1++;
332  else if (temptpc==2) nwiresTPC2++;
333  else if (temptpc==3) nwiresTPC3++;
334  else if (temptpc==4) nwiresTPC4++;
335  else if (temptpc==5) nwiresTPC5++;
336  else if (temptpc==6) nwiresTPC6++;
337  else if (temptpc==7) nwiresTPC7++;
338  else mf::LogVerbatim("RobustHitFinder") << "Unknown TPC";
339  allchannels.push_back(wire->Channel());
340  }
341 
342  for (size_t i_t0 = 0; i_t0 < t0Handle->size(); i_t0++)
343  {
344  art::Ptr<anab::T0> pt0(t0Handle,i_t0);
345  t0 = pt0->Time();
346 
347  unsigned int tick0 = static_cast<unsigned int>(fClks->TPCG4Time2Tick(t0));
348  mf::LogVerbatim("RobustHitFinder") << "Found t0 at TPC tick = " << tick0;
349 
350  std::vector<art::Ptr<raw::ExternalTrigger> > trigvec = triggers.at(i_t0);
351 
352  std::vector<unsigned int> evtTriggers;
353  for (auto const &trig : trigvec) evtTriggers.push_back(trig->GetTrigID());
354 
355  if (!ValidTrigger(evtTriggers,c1,c2,trignum))
356  {
357  continue;
358  }
359 
361  {
362  c1x = fCounterPositionMap[c1].first.X();
363  c1y = fCounterPositionMap[c1].first.Y();
364  c1z = fCounterPositionMap[c1].first.Z();
365  c2x = fCounterPositionMap[c2].first.X();
366  c2y = fCounterPositionMap[c2].first.Y();
367  c2z = fCounterPositionMap[c2].first.Z();
368  }
369  else
370  {
371  try
372  {
373  const geo::AuxDetGeo & ad = fGeom->AuxDet(c1);
374  double center[3] = {0.,0.,0.5*ad.Length()};
375  ad.GetCenter(center);
376  c1x = center[0];
377  c1y = center[1];
378  c1z = center[2];
379  }
380  catch (...)
381  {
382  std::cout << "AuxDet " << c1 << " not found. Aborting..." << std::endl;
383  continue;
384  }
385  try
386  {
387  const geo::AuxDetGeo & ad = fGeom->AuxDet(c2);
388  double center[3] = {0.,0.,0.5*ad.Length()};
389  ad.GetCenter(center);
390  c2x = center[0];
391  c2y = center[1];
392  c2z = center[2];
393  }
394  catch (...)
395  {
396  std::cout << "AuxDet " << c2 << " not found. Aborting..." << std::endl;
397  continue;
398  }
399  }
400 
401  int triggercode = 100*c1+c2; //gives a 4 digit number, ABCD, where AB is c1 ID and CD is c2 ID
402  hEventsByTrigger->Fill(triggercode);
403  fAllTree->Fill();
404 
406  if (trignum == 111)
407  {
408  fHorizRangeMin = -10;
409  fHorizRangeMax = 170;
410  fVertRangeMin = -50;
411  fVertRangeMax = 250;
412  fC1Vert = c1x;
413  fC1Horiz = c1z;
414  fC2Vert = c2x;
415  fC2Horiz = c2z;
416  distancecut *= fabs(TMath::Sin(TMath::ATan2(c1z-c2z,c1x-c2x)));
417  }
418  if (trignum == 112 || trignum == 113)
419  {
420  fHorizRangeMin = -50;
421  fHorizRangeMax = 250;
422  fVertRangeMin = -10;
423  fVertRangeMax = 170;
424  fC1Vert = c1z;
425  fC1Horiz = c1x;
426  fC2Vert = c2z;
427  fC2Horiz = c2x;
428  distancecut *= fabs(TMath::Cos(TMath::ATan2(c1z-c2z,c1x-c2x)));
429  }
430  double counterslope = ((fC1Vert-fC2Vert)/(fC1Horiz-fC2Horiz));
432  fFitAlg.SetParameter(1,1,counterslope-0.15,counterslope+0.15);
433  fFitAlg.SetParameter(2,0,-0.0002,0.0002);
435 
436  dune::ChanMap_t chanMap;
437 
438  for (size_t iwire = 0; iwire < wireHandle->size(); ++iwire)
439  {
440  art::Ptr<recob::Wire> pwire(wireHandle,iwire);
441  if (fCSP->IsBad(pwire->Channel())) continue;
442  art::Ptr<raw::RawDigit> prawdigit = rawdigits.at(iwire);
443 
444  if (pwire->View() != geo::kZ) continue;
445 
447 
448  chan.artWire = pwire;
449  chan.artRawDigit = prawdigit;
450 
451  chan.signalVec.clear();
452  chan.signalFilterVec.clear();
453 
454  chan.channelID = pwire->Channel();
455  chan.wireID = fGeom->ChannelToWire(chan.channelID)[0].Wire;
456  chan.tpcNum = fGeom->ChannelToWire(chan.channelID)[0].TPC;
457 
458  chan.signalVec = pwire->Signal();
459  chan.signalSize = chan.signalVec.size();
460 
461  double wirexyz[3];
462  fGeom->Wire(*(fGeom->ChannelToWire(chan.channelID).begin())).GetCenter(wirexyz);
463  chan.chanz = static_cast<float>(wirexyz[2]);
464 
466  tick0 + ((fSearchPostTicks>0) ? fSearchPostTicks : chan.signalSize));
470 
473 
474  if (temppedmean > fMaxPedMean || temppedmean < fMinPedMean || temppedrms > fMaxPedRms || temppedrms < fMinPedRms) continue;
475 
476  fHitFinderAlg.FindHits(chan);
477 
478  chanMap.emplace(std::make_pair(chan.channelID,chan));
479  }
480 
481  if (chanMap.size() == 0)
482  {
483  continue;
484  }
485 
486  dune::HitVec_t hitVec;
487  for (auto & chanmapitr : chanMap)
488  {
489  FillHitInformation(chanmapitr.second,hitVec,false);
490  }
491  if (hitVec.size() == 0)
492  {
493  continue;
494  }
495 
496  std::vector<dune::HitLineFitAlg::HitLineFitData> fitdata;
497  for (auto & hit : hitVec)
498  {
500  if (trignum == 111)
501  {
502  hlfd.hitHoriz = hit.hitz;
503  hlfd.hitVert = hit.hitx;
504  hlfd.hitHorizErrLo = hit.hiterrzlo;
505  hlfd.hitHorizErrHi = hit.hiterrzhi;
506  hlfd.hitVertErrLo = hit.hiterrxlo;
507  hlfd.hitVertErrHi = hit.hiterrxhi;
508  }
509  else if (trignum == 112 || trignum == 113)
510  {
511  hlfd.hitHoriz = hit.hitx;
512  hlfd.hitVert = hit.hitz;
513  hlfd.hitHorizErrLo = hit.hiterrxlo;
514  hlfd.hitHorizErrHi = hit.hiterrxhi;
515  hlfd.hitVertErrLo = hit.hiterrzlo;
516  hlfd.hitVertErrHi = hit.hiterrzhi;
517  }
518  hlfd.hitREAL = false;
519  fitdata.push_back(hlfd);
520  }
521 
523  int retval = fFitAlg.FitLine(fitdata,fitresult);
524 
525  for (size_t i_h = 0; i_h < fitdata.size(); ++i_h)
526  {
527  hitVec[i_h].fitrealhit = fitdata[i_h].hitREAL;
528  }
529 
530  if (fMakeupMissedHits && retval == 1) MakeupMissedHits(chanMap,hitVec);
531 
532  if (retval == 1 || retval == 0)
533  {
534  if (fitresult.fitsuccess)
535  {
536  hTracksByTrigger->Fill(triggercode);
537  fGoodTree->Fill();
538  }
539 
540  for (auto & hit : hitVec)
541  {
542  dune::ChannelInformation & chan = chanMap[hit.channelID];
543 
544  SetTreeVariables(chan,hit,fitresult);
545 
546  if (hit.fitrealhit && fitresult.fitsuccess)
547  {
548  hcol.emplace_back(hit.artHit,chan.artWire,chan.artRawDigit);
549  }
550  }
551  }
552  if (fSimulation) SetTreeVariablesMCTruth(chanMap,e);
553  if (fMakeTree) fTree->Fill();
554  }
555  hcol.put_into(e);
556 }
557 
559 {
560  DAQToOffline::MakeCounterPositionMap("","counterInformation.txt",fCounterPositionMap,0,0,0);
561 
562  if (fMakeTree)
563  {
564  hEventsByTrigger = fTfs->make<TH1I>("hEventsByTrigger","Num events by trigger",3716,0,3716);
565  hTracksByTrigger = fTfs->make<TH1I>("hTracksByTrigger","Num reconstructed tracks by trigger",3716,0,3716);
566 
567  fAllTree = fTfs->make<TTree>("AllTree","AllTree");
568  fAllTree->Branch("run",&run,"run/I");
569  fAllTree->Branch("subrun",&subrun,"subrun/I");
570  fAllTree->Branch("event",&event,"event/I");
571  fAllTree->Branch("c1",&c1,"c1/i");
572  fAllTree->Branch("c2",&c2,"c2/i");
573  fAllTree->Branch("trignum",&trignum,"trignum/i");
574 
575  fGoodTree = fTfs->make<TTree>("GoodTree","GoodTree");
576  fGoodTree->Branch("run",&run,"run/I");
577  fGoodTree->Branch("subrun",&subrun,"subrun/I");
578  fGoodTree->Branch("event",&event,"event/I");
579  fGoodTree->Branch("c1",&c1,"c1/i");
580  fGoodTree->Branch("c2",&c2,"c2/i");
581  fGoodTree->Branch("trignum",&trignum,"trignum/i");
582 
583  fTree = fTfs->make<TTree>("RobustHitFinder","RobustHitFinder");
584  fTree->Branch("run",&run,"run/I");
585  fTree->Branch("subrun",&subrun,"subrun/I");
586  fTree->Branch("event",&event,"event/I");
587  fTree->Branch("t0",&t0,"t0/D");
588  fTree->Branch("c1",&c1,"c1/i");
589  fTree->Branch("c2",&c2,"c2/i");
590  fTree->Branch("trignum",&trignum,"trignum/i");
591  fTree->Branch("c1x",&c1x,"c1x/F");
592  fTree->Branch("c1y",&c1y,"c1y/F");
593  fTree->Branch("c1z",&c1z,"c1z/F");
594  fTree->Branch("c2x",&c2x,"c2x/F");
595  fTree->Branch("c2y",&c2y,"c2y/F");
596  fTree->Branch("c2z",&c2z,"c2z/F");
597  fTree->Branch("distancecut",&distancecut,"distancecut/F");
598  fTree->Branch("allchannels",&allchannels);
599  fTree->Branch("channel",&channel);
600  fTree->Branch("wire",&wire);
601  fTree->Branch("tpc",&tpc);
602  fTree->Branch("signalsize",&signalsize);
603  //fTree->Branch("signal",&signal);
604  //fTree->Branch("signalFilter",&signalFilter);
605  fTree->Branch("baseline",&baseline);
606  fTree->Branch("rms",&rms);
607  fTree->Branch("baselineFilter",&baselineFilter);
608  fTree->Branch("rmsFilter",&rmsFilter);
609  fTree->Branch("pedmean",&pedmean);
610  fTree->Branch("pedrms",&pedrms);
611  fTree->Branch("integral",&integral);
612  fTree->Branch("integralFilter",&integralFilter);
613  fTree->Branch("sumADC",&sumADC);
614  fTree->Branch("sigmaintegral",&sigmaintegral);
615  fTree->Branch("sigmaintegralFilter",&sigmaintegralFilter);
616  fTree->Branch("amplitude",&amplitude);
617  fTree->Branch("amplitudeFilter",&amplitudeFilter);
618  fTree->Branch("peaktick",&peaktick);
619  fTree->Branch("peaktickFilter",&peaktickFilter);
620  fTree->Branch("peaktime",&peaktime);
621  fTree->Branch("peaktimeFilter",&peaktimeFilter);
622  fTree->Branch("begintick",&begintick);
623  fTree->Branch("endtick",&endtick);
624  fTree->Branch("width",&width);
625  fTree->Branch("hitx",&hitx);
626  fTree->Branch("hity",&hity);
627  fTree->Branch("hitz",&hitz);
628  fTree->Branch("hiterrxlo",&hiterrxlo);
629  fTree->Branch("hiterrxhi",&hiterrxhi);
630  fTree->Branch("hiterrylo",&hiterrylo);
631  fTree->Branch("hiterryhi",&hiterryhi);
632  fTree->Branch("hiterrzlo",&hiterrzlo);
633  fTree->Branch("hiterrzhi",&hiterrzhi);
634  fTree->Branch("perpdist",&perpdist);
635  fTree->Branch("hitt",&hitt);
636  fTree->Branch("driftdist",&driftdist);
637  fTree->Branch("countercut",&countercut);
638  fTree->Branch("fitconstant",&fitconstant,"fitconstant/F");
639  fTree->Branch("fitconstanterr",&fitconstanterr,"fitconstanterr/F");
640  fTree->Branch("fitlinear",&fitlinear,"fitlinear/F");
641  fTree->Branch("fitlinearerr",&fitlinearerr,"fitlinearerr/F");
642  fTree->Branch("fitquadratic",&fitquadratic,"fitquadratic/F");
643  fTree->Branch("fitquadraticerr",&fitquadraticerr,"fitquadraticerr/F");
644  fTree->Branch("fitchi2",&fitchi2,"fitchi2/F");
645  fTree->Branch("fitsumsqrresidual",&fitsumsqrresidual,"fitsumsqrresidual/F");
646  fTree->Branch("fitndf",&fitndf,"fitndf/F");
647  fTree->Branch("fitmle",&fitmle,"fitmle/F");
648  fTree->Branch("fitsuccess",&fitsuccess,"fitsuccess/O");
649  fTree->Branch("fitrealhit",&fitrealhit);
650  fTree->Branch("assumedhit",&assumedhit);
651  fTree->Branch("segmentlength",&segmentlength);
652  fTree->Branch("numGoodHitsChan",&numGoodHitsChan);
653  fTree->Branch("ismctruth",&ismctruth);
654  fTree->Branch("nwiresTPC0",&nwiresTPC0,"nwiresTPC0/I");
655  fTree->Branch("nwiresTPC1",&nwiresTPC1,"nwiresTPC1/I");
656  fTree->Branch("nwiresTPC2",&nwiresTPC2,"nwiresTPC2/I");
657  fTree->Branch("nwiresTPC3",&nwiresTPC3,"nwiresTPC3/I");
658  fTree->Branch("nwiresTPC4",&nwiresTPC4,"nwiresTPC4/I");
659  fTree->Branch("nwiresTPC5",&nwiresTPC5,"nwiresTPC5/I");
660  fTree->Branch("nwiresTPC6",&nwiresTPC6,"nwiresTPC6/I");
661  fTree->Branch("nwiresTPC7",&nwiresTPC7,"nwiresTPC7/I");
662  }
663 }
664 
666 {
667  fMakeTree = p.get<bool>("MakeTree");
668  fHitGeomDistanceCut = p.get<float>("HitGeomDistanceCut");
669  fCounterT0ModuleLabel = p.get<std::string>("CounterT0ModuleLabel");
670  fWireModuleLabel = p.get<std::string>("WireModuleLabel");
671  fEfield = p.get<std::vector<float> >("Efield");
672  fSearchPreTicks = p.get<int>("SearchPreTicks",-1);
673  fSearchPostTicks = p.get<int>("SearchPostTicks",-1);
674  fMinPedMean = p.get<float>("MinPedMean");
675  fMaxPedMean = p.get<float>("MaxPedMean");
676  fMinPedRms = p.get<float>("MinPedRms");
677  fMaxPedRms = p.get<float>("MaxPedRms");
678  fMakeupMissedHits = p.get<bool>("MakeupMissedHits",true);
679  fMissedBufferTicksLow = p.get<int>("MissedBufferTicksLow",0);
680  fMissedBufferTicksHigh = p.get<int>("MissedBufferTicksHigh",0);
681  fUseMeasuredCounterPositions = p.get<bool>("UseMeasuredCounterPositions");
682  fPreHitTicks = p.get<int>("PreHitTicks",100);
683  fPostHitTicks = p.get<int>("PostHitTicks",200);
684  fConstantHitWidth = p.get<bool>("ConstantHitWidth",true);
685  fSimulation = p.get<bool>("Simulation");
686  fMCScale = p.get<float>("PreviousMCScale",1.0);
687 }
688 
690 {
691  for (int i_tpc = 0; i_tpc<8; ++i_tpc)
692  {
693  float chanbeginz = 99999;
694  float chanendz = -99999;
695  int chanbeginid = 99999;
696  int chanendid = -99999;
697  for (auto & chanitr : chanMap)
698  {
699  if (chanitr.second.tpcNum != i_tpc) continue;
700  chanitr.second.nGoodHits = 0;
701  for (auto & hit : hitVec)
702  {
703  if (hit.channelID == chanitr.first && hit.fitrealhit)
704  {
705  chanitr.second.goodHitStartTick = hit.hitBeginTick;
706  chanitr.second.goodHitEndTick = hit.hitEndTick;
707  chanitr.second.nGoodHits++;
708  }
709  }
710  if (chanitr.second.nGoodHits > 0 && chanitr.second.chanz < chanbeginz)
711  {
712  chanbeginz = chanitr.second.chanz;
713  chanbeginid = chanitr.first;
714  }
715  if (chanitr.second.nGoodHits > 0 && chanitr.second.chanz > chanendz)
716  {
717  chanendz = chanitr.second.chanz;
718  chanendid = chanitr.first;
719  }
720  }
721 
722  if (chanendz < chanbeginz)
723  {
724  mf::LogError("RobustHitFinder") << "Problem: chanbeginz=" << chanbeginz << " is greater than chanendz=" << chanendz;
725  continue;
726  }
727 
728  if (chanbeginid == 99999 || chanendid == -99999)
729  {
730  mf::LogError("RobustHitFinder") << "Track not found. No assumed hits added." << std::endl;
731  continue;
732  }
733 
734  if (chanMap[chanbeginid].tpcNum % 2 != chanMap[chanendid].tpcNum % 2)
735  {
736  mf::LogError("RobustHitFinder") << "Track crosses APA. Don't know how to deal with this yet.";
737  continue;
738  }
739 
740  std::vector<std::pair<float,int> > trackchans;
741  for (auto const & chanitr : chanMap)
742  {
743  if (chanitr.second.tpcNum != i_tpc) continue;
744  if (chanitr.second.chanz >= chanbeginz && chanitr.second.chanz <= chanendz && chanitr.second.tpcNum % 2 == chanMap[chanbeginid].tpcNum % 2)
745  {
746  trackchans.push_back(std::make_pair(chanitr.second.chanz,chanitr.first));
747  }
748  }
749  std::sort(trackchans.begin(),trackchans.end());
750 
751  for (size_t i_tc = 0; i_tc < trackchans.size(); ++i_tc)
752  {
753  dune::ChannelInformation ch = chanMap[trackchans[i_tc].second];
754 
755  if (ch.nGoodHits == 0 && ch.chanz >= chanbeginz && ch.chanz <= chanendz)
756  {
757  int sub = 0;
758  dune::ChannelInformation chnearlow = ch;
759  while (chnearlow.nGoodHits != 1 && chnearlow.chanz >= chanbeginz)
760  {
761  ++sub;
762  chnearlow = chanMap[trackchans[i_tc-sub].second];
763  }
764 
765  int add = 0;
766  dune::ChannelInformation chnearhigh = ch;
767  while (chnearhigh.nGoodHits != 1 && chnearhigh.chanz <= chanendz)
768  {
769  ++add;
770  chnearhigh = chanMap[trackchans[i_tc+add].second];
771  }
772 
773  if (ch.tpcNum != chnearlow.tpcNum || ch.tpcNum != chnearhigh.tpcNum) continue;
774 
775  int cls = chnearlow.goodHitStartTick;
776  int chs = chnearhigh.goodHitStartTick;
777  int cle = chnearlow.goodHitEndTick;
778  int che = chnearhigh.goodHitEndTick;
779  double clz = chnearlow.chanz;
780  double chz = chnearhigh.chanz;
781 
782  int chstart = cls + (((chs - cls) / (chz - clz)) * (ch.chanz - clz));
783  int chend = cle + (((che - cle) / (chz - clz)) * (ch.chanz - clz));
784 
785  if (!fConstantHitWidth)
786  {
787  chstart -= fMissedBufferTicksLow;
788  chend += fMissedBufferTicksHigh;
789  }
790 
791  if (chz-clz < 0) continue; // something really weird must've happened
792 
793  ch.pulse_ends.push_back(std::make_pair(chstart,chend));
794 
795  FillHitInformation(ch,hitVec,true);
796  }
797  }
798  }
799 }
800 
802 {
803  fitsuccess = fitresult.fitsuccess;
804  TF1 * model = new TF1("model","pol2",fHorizRangeMin,fHorizRangeMax);
805  if (fitsuccess)
806  {
807  fitconstant = fitresult.bestVal.at(0);
808  fitconstanterr = fitresult.bestValError.at(0);
809  fitlinear = fitresult.bestVal.at(1);
810  fitlinearerr = fitresult.bestValError.at(1);
811  fitquadratic = fitresult.bestVal.at(2);
812  fitquadraticerr = fitresult.bestValError.at(2);
813  fitchi2 = fitresult.chi2;
814  fitsumsqrresidual = fitresult.sum2resid;
815  fitmle = fitresult.mle;
816  fitndf = fitresult.ndf;
817  model->SetParameters(fitconstant,fitlinear,fitquadratic);
818  }
819 
820  pedmean.push_back(temppedmean);
821  pedrms.push_back(temppedrms);
822  channel.push_back(chan.channelID);
823  wire.push_back(chan.wireID);
824  tpc.push_back(chan.tpcNum);
825  signalsize.push_back(chan.signalSize);
826  //signal = chan.signalVec;
827  //signalFilter = chan.signalFilterVec;
828  baseline.push_back(chan.baseline);
829  rms.push_back(chan.rms);
830  baselineFilter.push_back(chan.baselineFilter);
831  rmsFilter.push_back(chan.rmsFilter);
832  integral.push_back(hit.hitIntegral);
833  integralFilter.push_back(hit.hitIntegralFilter);
834  sumADC.push_back(hit.hitSumADC);
835  sigmaintegral.push_back(hit.hitSigmaIntegral);
837  amplitude.push_back(hit.hitAmplitude);
838  amplitudeFilter.push_back(hit.hitAmplitudeFilter);
839  peaktick.push_back(hit.hitPeakTick);
840  peaktickFilter.push_back(hit.hitPeakTickFilter);
841  peaktime.push_back(hit.hitPeakTime);
842  peaktimeFilter.push_back(hit.hitPeakTimeFilter);
843  begintick.push_back(hit.hitBeginTick);
844  endtick.push_back(hit.hitEndTick);
845  width.push_back(hit.hitWidth);
846  hitx.push_back(hit.hitx);
847  hity.push_back(hit.hity);
848  hitz.push_back(hit.hitz);
849  hiterrxlo.push_back(hit.hiterrxlo);
850  hiterrxhi.push_back(hit.hiterrxhi);
851  hiterrylo.push_back(hit.hiterrylo);
852  hiterryhi.push_back(hit.hiterryhi);
853  hiterrzlo.push_back(hit.hiterrzlo);
854  hiterrzhi.push_back(hit.hiterrzhi);
855  perpdist.push_back(hit.perpdist);
856  hitt.push_back(hit.hitt);
857  driftdist.push_back(hit.driftdist);
858  countercut.push_back(hit.countercut);
859  fitrealhit.push_back(hit.fitrealhit);
860  assumedhit.push_back(hit.assumedhit);
861  numGoodHitsChan.push_back(chan.nGoodHits);
862  ismctruth.push_back(false);
863 
864  float tempsegmentlength = 0.449;
865  if (fitsuccess && hit.fitrealhit)
866  {
867  if (trignum == 111)
868  {
869  double thetayz = TMath::ATan2(model->Eval(hit.hitz+1)-model->Eval(hit.hitz-1),2);
870  double tan2thetayz = TMath::Power(TMath::Tan(thetayz),2);
871  double y2z2 = ((c1y-c2y)*(c1y-c2y))/((c1z-c2z)*(c1z-c2z));
872  double projL = sqrt(1+tan2thetayz+y2z2);
873  tempsegmentlength *= static_cast<float>(projL);
874  }
875  else if (trignum == 112 || trignum == 113)
876  {
877  double thetayx = TMath::ATan2(2,model->Eval(hit.hitx+1)-model->Eval(hit.hitx-1));
878  double tan2thetayx = TMath::Power(TMath::Tan(thetayx),2);
879  double y2x2 = ((c1y-c2y)*(c1y-c2y))/((c1x-c2x)*(c1x-c2x));
880  double projL = sqrt(1+tan2thetayx*(1+y2x2));
881  tempsegmentlength *= static_cast<float>(projL);
882  }
883  }
884  segmentlength.push_back(tempsegmentlength);
885 }
886 
888 {
889  //TF1 * gaus = new TF1("gaus","([0]/([2]*sqrt(2*3.1415926)))*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))+[3]+x*[4]",0,32000);
890  //gaus->SetNpx(64000);
891 
892  for (size_t i_hit = 0; i_hit < chan.pulse_ends.size(); i_hit++)
893  {
895  hit.channelID = chan.channelID;
896  int begin_index,end_index;
897  if (fConstantHitWidth)
898  {
899  std::vector<float>::iterator bi = chan.signalVec.begin()+chan.pulse_ends[i_hit].first;
900  std::vector<float>::iterator ei = chan.signalVec.begin()+chan.pulse_ends[i_hit].second;
901  float peak = std::distance(chan.signalVec.begin(),std::max_element(bi,ei));
902  begin_index = (int)(peak-fPreHitTicks);
903  end_index = (int)(peak+fPostHitTicks);
904  }
905  else
906  {
907  begin_index = chan.pulse_ends[i_hit].first;
908  end_index = chan.pulse_ends[i_hit].second;
909  }
910  hit.hitBeginTick = begin_index;
911  hit.hitEndTick = end_index;
912  std::vector<float>::iterator beginitr = chan.signalVec.begin()+begin_index;
913  std::vector<float>::iterator enditr = chan.signalVec.begin()+end_index;
914  std::vector<float> pulse(beginitr,enditr);
915  std::vector<float>::iterator fbeginitr = chan.signalFilterVec.begin()+begin_index;
916  std::vector<float>::iterator fenditr = chan.signalFilterVec.begin()+end_index;
917  std::vector<float> pulseFilter(fbeginitr,fenditr);
918  hit.hitAmplitude = *std::max_element(beginitr,enditr)-chan.baseline;
919  hit.hitAmplitudeFilter = *std::max_element(fbeginitr,fenditr)-chan.baselineFilter;
920  hit.hitWidth = end_index-begin_index;
921  hit.hitIntegral = std::accumulate(beginitr,enditr,0)-(hit.hitWidth*chan.baseline);
922  hit.hitSumADC = hit.hitIntegral;
923  hit.hitIntegralFilter = std::accumulate(fbeginitr,fenditr,0)-(hit.hitWidth*chan.baselineFilter);
924  hit.hitSigmaIntegral = TMath::Sqrt(pulse.size())*TMath::RMS(pulse.size(),pulse.data());
925  hit.hitSigmaIntegralFilter = TMath::Sqrt(pulseFilter.size())*TMath::RMS(pulseFilter.size(),pulseFilter.data());
926  hit.hitPeakTick = std::distance(chan.signalVec.begin(),std::max_element(beginitr,enditr));
927  hit.hitPeakTickFilter = std::distance(chan.signalFilterVec.begin(),std::max_element(fbeginitr,fenditr));
930 
931  /*
932  gaus->SetParameter(1,hit.hitPeakTick);
933  gaus->SetParLimits(1,hit.hitPeakTick-10,hit.hitPeakTick+10);
934  gaus->SetParameter(2,5);
935  gaus->SetParLimits(2,1,end_index-begin_index);
936  TGraph * gr = new TGraph();
937  Int_t tick = 0;
938  for (auto adc : chan.signalVec)
939  {
940  gr->SetPoint(tick,(Double_t)tick,(Double_t)adc);
941  ++tick;
942  }
943  gr->Fit(gaus,"BQ0");
944 
945  hit.hitAmplitude = gaus->GetParameter(0)/(gaus->GetParameter(2)*sqrt(2*3.1415926));
946  hit.hitIntegral = gaus->GetParameter(0);
947  hit.hitSigmaIntegral = gaus->GetParError(0);
948  hit.hitSumADC = std::accumulate(beginitr,enditr,0)-((end_index-begin_index)*chan.baseline);
949  hit.hitWidth = gaus->GetParameter(2);
950  hit.hitPeakTick = gaus->GetParameter(1);
951  hit.hitPeakTime = fClks->TPCTick2TrigTime(hit.hitPeakTick);
952 
953  delete gr;
954  */
955 
956  hit.hitt = hit.hitPeakTime - t0/1000;
957 
958  hit.hitx = TimeToX(hit.hitt,chan.tpcNum);
959  hit.hity = 2.;
960  hit.hitz = chan.chanz;
961  hit.hiterrxlo = TMath::Sqrt(fabs(hit.hitx-TimeToDisplacement(fClks->TPCTick2TrigTime(hit.hitBeginTick)-t0/1000)));
962  hit.hiterrxhi = TMath::Sqrt(fabs(hit.hitx-TimeToDisplacement(fClks->TPCTick2TrigTime(hit.hitEndTick)-t0/1000)));
963  hit.hiterrylo = 0.;
964  hit.hiterryhi = 0.;
965  hit.hiterrzlo = 0.5/sqrt(12.0);
966  hit.hiterrzhi = 0.5/sqrt(12.0);
967  hit.driftdist = TimeToDriftDist(hit.hitt,chan.tpcNum);
968  hit.perpdist = hitGeomDist(TVector3(hit.hitx,2.,hit.hitz),TVector3(c1x,2.,c1z),TVector3(c2x,2.,c2z));
969  hit.countercut = (hit.perpdist < distancecut) ? true : false;
970  hit.fitrealhit = false;
971  hit.assumedhit = assumedHit;
972 
973  float tempsegmentlength = 0.449;
974  if (trignum == 111)
975  {
976  double thetayz = TMath::ATan2((hit.hitx+1)-(hit.hitx-1),2);
977  double tan2thetayz = TMath::Power(TMath::Tan(thetayz),2);
978  double y2z2 = ((c1y-c2y)*(c1y-c2y))/((c1z-c2z)*(c1z-c2z));
979  double projL = sqrt(1+tan2thetayz+y2z2);
980  tempsegmentlength *= static_cast<float>(projL);
981  }
982  else if (trignum == 112 || trignum == 113)
983  {
984  double thetayx = TMath::ATan2(2,(hit.hitz+1)-(hit.hitz-1));
985  double tan2thetayx = TMath::Power(TMath::Tan(thetayx),2);
986  double y2x2 = ((c1y-c2y)*(c1y-c2y))/((c1x-c2x)*(c1x-c2x));
987  double projL = sqrt(1+tan2thetayx*(1+y2x2));
988  tempsegmentlength *= static_cast<float>(projL);
989  }
990  hit.pitchgamma = tempsegmentlength;
991 
992 
993  recob::HitCreator temphit(*(chan.artWire),
994  fGeom->ChannelToWire(chan.channelID)[0],
995  hit.hitBeginTick,
996  hit.hitEndTick,
997  hit.hitWidth/sqrt(12.0),
998  hit.hitPeakTick,
999  hit.hitWidth/sqrt(12.0),
1000  hit.hitAmplitude,
1001  hit.hitAmplitude/sqrt(12.0),
1002  hit.hitIntegral,
1003  hit.hitSigmaIntegral,
1004  hit.hitSumADC,
1005  1,
1006  -1,
1007  0,
1008  (int)(hit.hitWidth-2));
1009  hit.artHit = temphit.move();
1010  if (hit.hitx > -400 && hit.countercut) hitVec.push_back(hit);
1011  }
1012 }
1013 
1015 {
1016  std::vector<MCHit> MCHitVec;
1020  if (!e.getByLabel(fWireModuleLabel,wireHandle)) return;
1022 
1023  for (auto sc : bt_serv->SimChannels())
1024  {
1025  bool channelexists = false;
1026  for (size_t i_wire = 0; i_wire < wireHandle->size(); ++i_wire)
1027  {
1028  art::Ptr<recob::Wire> pwire(wireHandle,i_wire);
1029  if (pwire->Channel() == sc->Channel())
1030  {
1031  channelexists = true;
1032  break;
1033  }
1034  }
1035  if (!channelexists) continue;
1036 
1037  // also remove channels from absentChannels list if it has a Sim Hit
1038  if (fCSP->IsBad(sc->Channel())) continue;
1039  if (fGeom->SignalType(sc->Channel()) != geo::kCollection) continue;
1040 
1041  // Initialize a MCHit object for this channel
1042  dune::MCHit mchit;
1043  auto const& tdcidemap = sc->TDCIDEMap();
1044  // Loop over the track id and energies deposited on this channel
1045  for (auto const& tdcIt : tdcidemap)
1046  {
1047  // Get the vector of IDEs
1048  auto const& ideVec = tdcIt.second;
1049 
1050  // Get the TDC (time) of the deposit
1051  unsigned short tdc = tdcIt.first;
1052 
1053  // Loop over IDEs
1054  for (auto const& ideIt : ideVec)
1055  {
1056  // We want only primary muons which cross the detector.
1057  // Don't include deltas and other secondary particles because we want to
1058  // see how well the reconstruction discounts these other tracks
1059  if (abs(ideIt.trackID) != 1) continue;
1060  const simb::MCParticle * part = pi_serv->TrackIdToParticle_P(abs(ideIt.trackID));
1061  if (part->Mother() != 0) continue;
1062 
1063  // If the track is a primary muon, then collect this IDE for this channel
1064  mchit.idevec.push_back(ideIt);
1065 
1066  // Sanity check: make sure channel numbers match
1067  if (mchit.idevec.size() > 1 && mchit.channel != sc->Channel())
1068  throw cet::exception("RobustHitFinder") << "Channel IDs don't match!";
1069 
1070  // Check : if the IDE vector contains different particle codes, the something is weird
1071  mchit.channel = sc->Channel();
1072  if (mchit.idevec.size() > 1 && mchit.pdg != part->PdgCode())
1073  throw cet::exception("RobustHitFinder") << "Pdg from previous IDE (" << mchit.pdg << ") doesn't equal this IDE Pdg (" << part->PdgCode() << ")";
1074 
1075  // Fill the MCHit object, appending energy deposits in the case where more than one IDE exists on the channel
1076  // Usually, it's either one or two IDEs per channel with signal. Very rarely more than 2
1077  mchit.pdg = part->PdgCode();
1078  mchit.particlevec.push_back(part);
1079  mchit.tdcvec.push_back(tdc);
1080  // Here, adjust the signal size for the MCScale value used in the DataOverlay step
1081  // Otherwise, the charges wouldn't make any sense
1082  mchit.chargevec.push_back(ideIt.numElectrons*fMCScale);
1083  }
1084  }
1085  if (mchit.idevec.size()==0) continue;
1086  MCHitVec.push_back(mchit);
1087  }
1088 
1089  for (auto mchit : MCHitVec)
1090  {
1091  // Collect information about MCHit
1092  int mcchannel = mchit.channel;
1093  //int numIDEs = mchit.idevec.size();
1094  double meanSimTime = TMath::Mean(mchit.tdcvec.begin(),mchit.tdcvec.end(),mchit.chargevec.begin()); // charge weighted mean
1095  double sigmaSimTime = TMath::RMS(mchit.tdcvec.begin(),mchit.tdcvec.end(),mchit.chargevec.begin()); // charge weighted RMS
1096  double totalSimCharge = std::accumulate(mchit.chargevec.begin(),mchit.chargevec.end(),0);
1097  totalSimCharge *= 1.602e-4; // convert # electrons to fC
1098  totalSimCharge *= 14.0*7.615; // convert fC to mV
1099  totalSimCharge *= 2.808; // convert to ADC
1100 
1101  pedmean.push_back(temppedmean);
1102  pedrms.push_back(temppedrms);
1103  channel.push_back(mcchannel); ///
1104  wire.push_back(chanMap.at(mcchannel).wireID);
1105  tpc.push_back(chanMap.at(mcchannel).tpcNum); ///
1106  //signalsize.push_back(chanMap.at(mcchannel).signalSize);
1107  //signal = chanMap.at(mcchannel).signalVec;
1108  //signalFilter = chanMap.at(mcchannel).signalFilterVec;
1109  baseline.push_back(chanMap.at(mcchannel).baseline); ///
1110  rms.push_back(chanMap.at(mcchannel).rms); ///
1111  baselineFilter.push_back(chanMap.at(mcchannel).baselineFilter);
1112  rmsFilter.push_back(chanMap.at(mcchannel).rmsFilter);
1113  integral.push_back(totalSimCharge); ///
1114  integralFilter.push_back(0);
1115  sumADC.push_back(0);
1116  sigmaintegral.push_back(0);
1117  sigmaintegralFilter.push_back(0);
1118  amplitude.push_back(0);
1119  amplitudeFilter.push_back(0);
1120  peaktick.push_back(0);
1121  peaktickFilter.push_back(0);
1122  peaktime.push_back(0);
1123  peaktimeFilter.push_back(0);
1124  begintick.push_back(0);
1125  endtick.push_back(0);
1126  width.push_back(sigmaSimTime); ///
1127  hitx.push_back(TimeToX(meanSimTime,chanMap.at(mcchannel).tpcNum));
1128  hity.push_back(0);
1129  hitz.push_back(chanMap.at(mcchannel).chanz);
1130  hiterrxlo.push_back(0);
1131  hiterrxhi.push_back(0);
1132  hiterrylo.push_back(0);
1133  hiterryhi.push_back(0);
1134  hiterrzlo.push_back(0);
1135  hiterrzhi.push_back(0);
1136  perpdist.push_back(0);
1137  hitt.push_back(meanSimTime); ///
1138  driftdist.push_back(0);
1139  countercut.push_back(0);
1140  fitrealhit.push_back(0);
1141  assumedhit.push_back(0);
1142  numGoodHitsChan.push_back(chanMap.at(mcchannel).nGoodHits);
1143  ismctruth.push_back(true);
1144 
1145  float tempsegmentlength = 0.449;
1146  if (trignum == 111)
1147  {
1148  double thetayz = TMath::ATan2(c1x-c2x,c1z-c2z);
1149  double tan2thetayz = TMath::Power(TMath::Tan(thetayz),2);
1150  double y2z2 = ((c1y-c2y)*(c1y-c2y))/((c1z-c2z)*(c1z-c2z));
1151  double projL = sqrt(1+tan2thetayz+y2z2);
1152  tempsegmentlength *= static_cast<float>(projL);
1153  }
1154  else if (trignum == 112 || trignum == 113)
1155  {
1156  double thetayx = TMath::ATan2(c1x-c2x,c1z-c2z);
1157  double tan2thetayx = TMath::Power(TMath::Tan(thetayx),2);
1158  double y2x2 = ((c1y-c2y)*(c1y-c2y))/((c1x-c2x)*(c1x-c2x));
1159  double projL = sqrt(1+tan2thetayx*(1+y2x2));
1160  tempsegmentlength *= static_cast<float>(projL);
1161  }
1162  segmentlength.push_back(tempsegmentlength); ///
1163  }
1164 }
1165 
1166 bool dune::RobustHitFinder::ValidTrigger(std::vector<unsigned int> evtTriggers, unsigned int & c1arg, unsigned int & c2arg, unsigned int & trignumarg)
1167 {
1168  c1arg=999; c2arg=999;
1169  int contains_111 = 0, contains_112 = 0, contains_113 = 0;
1170  int contains_Ntrigs = 0, contains_NU = 0, contains_NL = 0, contains_SU = 0, contains_SL = 0;
1171  int contains_EL = 0, contains_WU = 0, contains_TEL = 0;
1172  for (size_t i_c = 0; i_c < evtTriggers.size(); i_c++)
1173  {
1174  unsigned int trigID = evtTriggers[i_c];
1175  // for c2: trigID is an unsigned int and always >= 0
1176  //if (trigID >= 0 && trigID <= 5) contains_SL++;
1177  if (trigID <= 5) contains_SL++;
1178  if (trigID >= 6 && trigID <= 15) contains_EL++;
1179  if (trigID >= 16 && trigID <= 21) contains_NL++;
1180  if (trigID >= 22 && trigID <= 27) contains_NU++;
1181  if (trigID >= 28 && trigID <= 37) contains_WU++;
1182  if (trigID >= 38 && trigID <= 43) contains_SU++;
1183  if (trigID >= 44 && trigID <= 92) contains_TEL++;
1184  if (trigID == 111) contains_111++;
1185  if (trigID == 112) contains_112++;
1186  if (trigID == 113) contains_113++;
1187  contains_Ntrigs++;
1188  }
1189  if (contains_111 + contains_112 + contains_113 != 1) return false; // too many/few coincidences!
1190  if (contains_TEL &&
1191  (contains_NU || contains_NL || contains_SU || contains_SL || contains_EL || contains_WU)) return false; // track probably doesn't go through detector
1192  if (contains_Ntrigs != 3) return false; // too much/little going on!
1193  if (contains_111 && (contains_NU || contains_NL || contains_SU || contains_SL)) return false; // 111 should not have NU/NL/SU/SL
1194  if (contains_112 && (contains_EL || contains_WU || contains_SU || contains_NL)) return false; // 112 should not have EL/WU/SU/NL
1195  if (contains_113 && (contains_EL || contains_WU || contains_NU || contains_SL)) return false; // 113 should not have EL/WU/NU/SL
1196  if (contains_111 && (!contains_EL || !contains_WU)) return false; // incomplete trigger
1197  if (contains_112 && (!contains_NU || !contains_SL)) return false; // incomplete trigger
1198  if (contains_113 && (!contains_SU || !contains_NL)) return false; // incomplete trigger
1199 
1200  std::vector<unsigned int> counterIDs;
1201  trignumarg = 0;
1202  for (size_t i_c = 0; i_c < evtTriggers.size(); i_c++)
1203  {
1204  unsigned int trigID = evtTriggers[i_c];
1205  if (trigID >= 44 && trigID <= 100) continue;
1206  if (trigID >= 111 && trigID <= 113)
1207  {
1208  trignumarg = trigID;
1209  continue;
1210  }
1211  counterIDs.push_back(trigID);
1212  }
1213  if (counterIDs.size() != 2) return false;
1214  if (trignumarg == 0) return false;
1215 
1216  if (trignumarg == 112 || trignumarg == 113)
1217  {
1218  if (fCounterPositionMap[counterIDs[0]].first.X() > fCounterPositionMap[counterIDs[1]].first.X())
1219  {
1220  c1arg = counterIDs[0];
1221  c2arg = counterIDs[1];
1222  }
1223  else
1224  {
1225  c1arg = counterIDs[1];
1226  c2arg = counterIDs[0];
1227  }
1228  }
1229  else if (trignumarg == 111)
1230  {
1231  if (fCounterPositionMap[counterIDs[0]].first.Z() > fCounterPositionMap[counterIDs[1]].first.Z())
1232  {
1233  c1arg = counterIDs[0];
1234  c2arg = counterIDs[1];
1235  }
1236  else
1237  {
1238  c1arg = counterIDs[1];
1239  c2arg = counterIDs[0];
1240  }
1241  }
1242  if (c1arg == c2arg) return false;
1243  if (c1arg == 999 || c2arg == 999) return false;
1244 
1245  return true;
1246 }
1247 
1248 float dune::RobustHitFinder::TimeToDriftDist(float thistime, unsigned int thistpc)
1249 {
1250  float vd = fDetProp->DriftVelocity(fEfield[0]);
1251  float v1 = fDetProp->DriftVelocity(fEfield[1]);
1252  float v2 = fDetProp->DriftVelocity(fEfield[2]);
1253  float v3 = fDetProp->DriftVelocity(fEfield[3]);
1254 
1255  float ld = -1, l1 = -1, l2 = -1, l3 = -1;
1256 
1257  try
1258  {
1259  if (thistpc == 1 || thistpc == 3 || thistpc == 5 || thistpc == 7)
1260  {
1261  ld = fGeom->TPC(thistpc).MaxX()-0.511;
1262  l1 = 0.511-fGeom->TPC(thistpc).PlaneLocation(0)[0];
1263  l2 = fGeom->TPC(thistpc).PlaneLocation(0)[0]-fGeom->TPC(thistpc).PlaneLocation(1)[0];
1264  l3 = fGeom->TPC(thistpc).PlaneLocation(1)[0]-fGeom->TPC(thistpc).PlaneLocation(2)[0];
1265  }
1266  else if (thistpc == 0 || thistpc == 2 || thistpc == 4 || thistpc == 6)
1267  {
1268  ld = -8.490-fGeom->TPC(thistpc).MinX();
1269  l1 = fGeom->TPC(thistpc).PlaneLocation(0)[0]+8.490;
1270  l2 = fGeom->TPC(thistpc).PlaneLocation(1)[0]-fGeom->TPC(thistpc).PlaneLocation(0)[0];
1271  l3 = fGeom->TPC(thistpc).PlaneLocation(2)[0]-fGeom->TPC(thistpc).PlaneLocation(1)[0];
1272  }
1273  }
1274  catch (cet::exception &e)
1275  {
1276  mf::LogError("RobustHitFinder") << e;
1277  return -99999;
1278  }
1279 
1280  if (ld < 0 || l1 < 0 || l2 < 0 || l3 < 0) return -99998;
1281 
1282  float t3max = l3/v3;
1283  float t2max = t3max+(l2/v2);
1284  float t1max = t2max+(l1/v1);
1285  float tdmax = t1max+(ld/vd);
1286 
1287  if (0 <= thistime && thistime < t3max) return v3*thistime;
1288  else if (t3max <= thistime && thistime < t2max) return l3+v2*(thistime-t3max);
1289  else if (t2max <= thistime && thistime < t1max) return l3+l2+v1*(thistime-t2max);
1290  else if (t1max <= thistime && thistime < tdmax) return l3+l2+l1+vd*(thistime-t1max);
1291 
1292  return -99997;
1293 }
1294 
1296 {
1297  return thistime*fDetProp->DriftVelocity(fEfield[0]);
1298 }
1299 
1300 float dune::RobustHitFinder::TimeToX(float thistime, unsigned int thistpc)
1301 {
1302  float driftdistance = TimeToDriftDist(thistime,thistpc);
1303  if (driftdistance < -89999) return driftdistance;
1304  try
1305  {
1306  if (thistpc == 1 || thistpc == 3 || thistpc == 5 || thistpc == 7)
1307  {
1308  return (fGeom->TPC(thistpc).PlaneLocation(2)[0]+driftdistance);
1309  }
1310  else if (thistpc == 0 || thistpc == 2 || thistpc == 4 || thistpc == 6)
1311  {
1312  return (fGeom->TPC(thistpc).PlaneLocation(2)[0]-driftdistance);
1313  }
1314  }
1315  catch (cet::exception &e)
1316  {
1317  mf::LogError("RobustHitFinder") << e;
1318  }
1319  return -99995;
1320 }
1321 
1322 float dune::RobustHitFinder::hitGeomDist(TVector3 hitloc, TVector3 trigloc1, TVector3 trigloc2)
1323 {
1324  return (((hitloc-trigloc1).Cross(hitloc-trigloc2)).Mag()/(trigloc2-trigloc1).Mag());
1325 }
1326 
1328 {
1329  run = -99999;
1330  subrun = -99999;
1331  event = -99999;
1332  t0 = -99999;
1333  c1 = 99999;
1334  c2 = 99999;
1335  trignum = 99999;
1336  c1x = -99999;
1337  c1y = -99999;
1338  c1z = -99999;
1339  c2x = -99999;
1340  c2y = -99999;
1341  c2z = -99999;
1342  channel.clear();
1343  tpc.clear();
1344  wire.clear();
1345  distancecut = -99999;
1346  signalsize.clear();
1347  signal.clear();
1348  signalFilter.clear();
1349  baseline.clear();
1350  rms.clear();
1351  baselineFilter.clear();
1352  rmsFilter.clear();
1353  pedmean.clear();
1354  pedrms.clear();
1355  integral.clear();
1356  integralFilter.clear();
1357  sumADC.clear();
1358  sigmaintegral.clear();
1359  sigmaintegralFilter.clear();
1360  amplitude.clear();
1361  amplitudeFilter.clear();
1362  peaktick.clear();
1363  peaktickFilter.clear();
1364  peaktime.clear();
1365  peaktimeFilter.clear();
1366  begintick.clear();
1367  endtick.clear();
1368  width.clear();
1369  hitx.clear();
1370  hity.clear();
1371  hitz.clear();
1372  hiterrxlo.clear();
1373  hiterrxhi.clear();
1374  hiterrylo.clear();
1375  hiterryhi.clear();
1376  hiterrzlo.clear();
1377  hiterrzhi.clear();
1378  perpdist.clear();
1379  hitt.clear();
1380  driftdist.clear();
1381  countercut.clear();
1382  fitconstant = -99999;
1383  fitconstanterr = -99999;
1384  fitlinear = -99999;
1385  fitlinearerr = -99999;
1386  fitquadratic = -99999;
1387  fitquadraticerr = -99999;
1388  fitchi2 = -99999;
1389  fitsumsqrresidual = -99999;
1390  fitndf = -99999;
1391  fitmle = -99999;
1392  fitsuccess = false;
1393  fitrealhit.clear();
1394  assumedhit.clear();
1395  segmentlength.clear();
1396  nwiresTPC0 = 0;
1397  nwiresTPC1 = 0;
1398  nwiresTPC2 = 0;
1399  nwiresTPC3 = 0;
1400  nwiresTPC4 = 0;
1401  nwiresTPC5 = 0;
1402  nwiresTPC6 = 0;
1403  nwiresTPC7 = 0;
1404  numGoodHitsChan.clear();
1405 }
1406 
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
Pedestal provider class for DUNE.
code to link reconstructed objects back to the MC truth information
virtual float PedRms(raw::ChannelID_t ch) const =0
intermediate_table::iterator iterator
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
RobustHitFinder(fhicl::ParameterSet const &p)
std::vector< float > hiterrzhi
art::Ptr< raw::RawDigit > artRawDigit
std::vector< float > hiterrxhi
int PdgCode() const
Definition: MCParticle.h:211
void reconfigure(fhicl::ParameterSet const &p)
EventNumber_t event() const
Definition: DataViewImpl.cc:96
void SetTreeVariablesMCTruth(const ChanMap_t &chanMap, art::Event &e)
CLHEP::HepRandomEngine & fEngine
void FilterWaveform(std::vector< float > wf, std::vector< float > &fwf)
std::vector< std::vector< float > > signal
std::vector< float > hiterrxlo
std::vector< float > peaktime
std::vector< float > pedmean
std::vector< int > numGoodHitsChan
std::vector< float > hiterrylo
const double & Time() const
Definition: T0.h:44
std::string string
Definition: nybbler.cc:12
std::vector< bool > ismctruth
art::Ptr< recob::Wire > artWire
float TimeToDisplacement(float thistime)
const simb::MCParticle * TrackIdToParticle_P(int id) const
std::vector< float > driftdist
int Mother() const
Definition: MCParticle.h:212
void RobustRMSBase(std::vector< float > wf, float &bl, float &r)
std::vector< float > integralFilter
Declaration of signal hit object.
std::vector< float > signalFilterVec
std::vector< float > perpdist
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
Coord add(Coord c1, Coord c2)
Definition: restypedef.cpp:23
int FitLine(std::vector< HitLineFitData > &data, HitLineFitResults &bestfit)
std::vector< float > peaktick
std::vector< float > amplitudeFilter
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition of basic raw digits.
std::vector< std::pair< int, int > > pulse_ends
Planes which measure Z direction.
Definition: geo_types.h:128
raw::ChannelID_t channel
std::vector< float > integral
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
std::vector< float > rmsFilter
bool ValidTrigger(std::vector< unsigned int > evtTriggers, unsigned int &c1arg, unsigned int &c2arg, unsigned int &trignumarg)
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Particle class.
std::vector< float > baseline
std::vector< float > hiterrzlo
virtual double TPCG4Time2Tick(double g4time) const =0
Converts simulation time into a TPC electronics time tick.
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248
float TimeToDriftDist(float thistime, unsigned int thistpc)
void SetSeed(UInt_t seed)
Definition: HitLineFitAlg.h:58
std::vector< sim::IDE > idevec
std::map< int, float > bestValError
Definition: HitLineFitAlg.h:37
detinfo::DetectorClocks const * fClks
geo::View_t View() const
Returns the view the channel belongs to.
Definition: Wire.h:230
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:83
T abs(T value)
std::vector< bool > countercut
Helper functions to create a hit.
std::vector< float > peaktickFilter
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
const double e
A class handling a collection of hits and its associations.
Definition: HitCreator.h:508
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:69
void produce(art::Event &e) override
void beginJob()
Definition: Breakpoints.cc:14
raw::ChannelID_t Channel() const
Returns the ID of the channel (or InvalidChannelID)
Definition: Wire.h:231
void FindHits(dune::ChannelInformation &chan)
float TimeToX(float thistime, unsigned int thistpc)
std::vector< float > pedrms
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::vector< std::vector< float > > signalFilter
void SetParameter(int i, double startValue, double minValue, double maxValue)
void SetTreeVariables(const dune::ChannelInformation &chan, const dune::HitInformation &hit, const dune::HitLineFitAlg::HitLineFitResults &fitresult)
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:89
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:290
Class providing information about the quality of channels.
std::vector< float > fEfield
detinfo::DetectorProperties const * fDetProp
std::vector< float > baselineFilter
RunNumber_t run() const
Definition: DataViewImpl.cc:82
std::vector< float > sigmaintegralFilter
std::vector< unsigned short > tdcvec
Definition of data types for geometry description.
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:652
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
std::vector< float > peaktimeFilter
Detector simulation of raw signals on wires.
Conversion of times between different formats and references.
std::vector< float > Signal() const
Return a zero-padded full length vector filled with RoI signal.
Definition: Wire.cxx:47
std::vector< const simb::MCParticle * > particlevec
Encapsulate the geometry of an auxiliary detector.
ProducesCollector & producesCollector() noexcept
p
Definition: test.py:223
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
std::vector< HitInformation > HitVec_t
art::ServiceHandle< geo::Geometry > fGeom
std::vector< bool > assumedhit
double Length() const
Definition: AuxDetGeo.h:102
dune::RMSHitFinderAlg fHitFinderAlg
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
std::vector< float > amplitude
std::vector< float > hiterryhi
def center(depos, point)
Definition: depos.py:117
void FillHitInformation(dune::ChannelInformation &chan, dune::HitVec_t &hitVec, bool assumedHit)
object containing MC truth information necessary for making RawDigits and doing back tracking ...
Interface for experiment-specific channel quality info provider.
art::ServiceHandle< art::TFileService > fTfs
virtual float PedMean(raw::ChannelID_t ch) const =0
Retrieve pedestal information.
void SetHorizVertRanges(float hmin, float hmax, float vmin, float vmax)
const std::vector< art::Ptr< sim::SimChannel > > & SimChannels() const
std::map< int, ChannelInformation > ChanMap_t
Declaration of basic channel signal object.
void MakeCounterPositionMap(std::string CounterDir, std::string CounterFile, std::map< unsigned int, std::pair< TVector3, std::vector< TVector3 > > > &CounterPositionMap, double fExtendCountersX=0, double fExtendCountersY=0, double fExtendCountersZ=0)
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
virtual double TPCTick2TrigTime(double tick) const =0
Converts a TPC time (in ticks) into a trigger time [µs].
void MakeupMissedHits(ChanMap_t &chanMap, HitVec_t &hitVec)
std::vector< float > sumADC
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:69
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:27
Interface for experiment-specific service for channel quality info.
float hitGeomDist(TVector3 hitloc, TVector3 trigloc1, TVector3 trigloc2)
std::vector< float > sigmaintegral
void SetSearchTicks(int s, int e)
Tools and modules for checking out the basics of the Monte Carlo.
std::vector< double > chargevec
std::map< unsigned int, std::pair< TVector3, std::vector< TVector3 > > > fCounterPositionMap
const double * PlaneLocation(unsigned int p) const
Definition: TPCGeo.cxx:382
std::vector< bool > fitrealhit
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< float > signalVec
std::vector< float > segmentlength
QTextStream & endl(QTextStream &s)
Event finding and building.
void GetCenter(double *xyz, double localz=0.0) const
Return the center position of an AuxDet.
Definition: AuxDetGeo.cxx:63
std::vector< int > allchannels
Signal from collection planes.
Definition: geo_types.h:142
const lariov::DetPedestalProvider & fPedestalRetrievalAlg