All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ICARUSHitFinder_module.cc
Go to the documentation of this file.
1 #ifndef ICARUSHitFINDER_H
2 #define ICARUSHitFINDER_H
3 
4 ////////////////////////////////////////////////////////////////////
5 //HIT FINDER THAT RUNS ON RAW SIGNALS INSTEAD OF DECONVOLUTED ONES.
6 //DEVELOPED INITIALLY FOR ICARUS-T600 OLD ELECTRONICS AT LNGS
7 //filippo.varanini@pd.infn.it .
8 ////////////////////////////////////////////////////////////////////
9 
10 // C/C++ standard libraries
11 #include <string>
12 #include <vector>
13 #include <fstream>
14 #include <set>
15 #include <cassert>
16 
17 //Framework
18 #include "fhiclcpp/ParameterSet.h"
19 #include "messagefacility/MessageLogger/MessageLogger.h"
20 #include "art/Framework/Core/ModuleMacros.h"
21 #include "art/Framework/Core/EDProducer.h"
22 #include "art/Framework/Principal/Event.h"
23 #include "art/Framework/Principal/Handle.h"
24 #include "canvas/Persistency/Common/Ptr.h"
25 #include "art/Framework/Services/Registry/ServiceHandle.h"
26 #include "canvas/Utilities/InputTag.h"
27 #include "canvas/Persistency/Common/FindOneP.h"
28 #include "art_root_io/TFileService.h"
29 #include "art/Utilities/ToolMacros.h"
30 #include "art/Utilities/make_tool.h"
31 
32 
33 //LArSoft
34 #include "larcore/Geometry/Geometry.h"
35 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
36 #include "larevt/Filters/ChannelFilter.h"
37 #include "lardataobj/RawData/RawDigit.h"
38 #include "lardataobj/RawData/raw.h"
39 #include "lardataobj/RecoBase/Hit.h"
40 #include "lardata/ArtDataHelper/HitCreator.h"
41 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
42 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
43 
44 //LArSoft From FFT
45 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
46 #include "larcore/Geometry/Geometry.h"
47 #include "larcorealg/Geometry/CryostatGeo.h"
48 #include "larcorealg/Geometry/TPCGeo.h"
49 #include "larcorealg/Geometry/PlaneGeo.h"
50 
51 #include "larreco/RecoAlg/GausFitCache.h" // hit::GausFitCache
52 #include "larreco/HitFinder/HitFinderTools/ICandidateHitFinder.h"
53 //#include "icaruscode/HitFinder/PeakFitterICARUS.h"
54 
55 //ROOT from CalData
56 #include "TComplex.h"
57 #include "TFile.h"
58 #include "TH2D.h"
59 #include "TF1.h"
60 
61 //ROOT From Gauss
62 #include "TH1D.h"
63 #include "TDecompSVD.h"
64 #include "TMath.h"
65 
66 namespace hit {
67  /// Customized function cache for ICARUS hit shape.
68  class ICARUShitFitCache: public hit::GausFitCache {
69 
70  public:
71  /// Constructor (see base class constructor).
72  ICARUShitFitCache(std::string const& new_name="ICARUShitFitCache")
73  : hit::GausFitCache(new_name)
74  {}
75 
76  /// ICARUS hit shape.
77  static Double_t fitf(Double_t const* x, Double_t const* par);
78 
79  protected:
80 
81  /// Creates and returns the function with specified number of peaks.
82  virtual TF1* CreateFunction(size_t nFunc) const
83  {
84  unsigned int const nPeaks = nFunc;
85  std::string const func_name = FunctionName(nFunc);
86  auto* pF = new TF1(func_name.c_str(), fitf, 0.0, 1.0, 1 + nFunc * 5);
87  pF->SetParName(0, "NPeaks");
88  pF->FixParameter(0, (double) nPeaks);
89  return pF;
90  } // CreateFunction()
91 
92  }; // ICARUShitFitCache
93 
94  /// Customized function cache for ICARUS long hit shape.
95  class ICARUSlongHitFitCache: public hit::GausFitCache {
96 
97  public:
98  /// Constructor (see base class constructor).
99  ICARUSlongHitFitCache(std::string const& new_name="ICARUSlongHitFitCache")
100  : hit::GausFitCache(new_name)
101  {}
102 
103  /// ICARUS hit shape.
104  static Double_t fitlong(Double_t const* x, Double_t const* par);
105 
106  protected:
107 
108  /// Creates and returns the function with specified number of peaks.
109  virtual TF1* CreateFunction(size_t nFunc) const
110  {
111  unsigned int const nPeaks = nFunc;
112  std::string const func_name = FunctionName(nFunc);
113  auto* pF = new TF1(func_name.c_str(), fitlong, 0.0, 1.0, 1 + nFunc * 7);
114  pF->SetParName(0, "NPeaks");
115  pF->FixParameter(0, (double) nPeaks);
116  return pF;
117  } // CreateFunction()
118 
119  }; // ICARUSlongHitFitCache
120 
121 
122 
123  class ICARUSHitFinder : public art::EDProducer {
124 
125  public:
126 
127  explicit ICARUSHitFinder(fhicl::ParameterSet const& pset);
128 
129  void produce(art::Event& evt);
130  void beginJob();
131  void endJob();
132  void reconfigure(fhicl::ParameterSet const& p);
133 
134  void expandHit(reco_tool::ICandidateHitFinder::HitCandidate& h, std::vector<float> holder, std::vector<reco_tool::ICandidateHitFinder::HitCandidate> how );
135  void computeBestLocalMean(std::vector<reco_tool::ICandidateHitFinder::HitCandidate> h, std::vector<float> holder, reco_tool::ICandidateHitFinder::MergeHitCandidateVec how, float& localmean);
136 
137  using ICARUSPeakFitParams_t = struct ICARUSPeakFitParams
138  {
139  float peakCenter;
140  float peakCenterError;
141  float peakSigma;
142  float peakSigmaError;
143  float peakAmplitude;
144  float peakAmplitudeError;
145  float peakTauLeft;
146  float peakTauLeftError;
147  float peakTauRight;
148  float peakTauRightError;
149  float peakFitWidth;
150  float peakFitWidthError;
151  float peakSlope;
152  float peakSlopeError;
153  float peakBaseline;
154  float peakBaselineError;
155  };
156  using ICARUSPeakParamsVec = std::vector<ICARUSPeakFitParams_t>;
157  void findMultiPeakParameters(const std::vector<float>&,
158  const reco_tool::ICandidateHitFinder::HitCandidateVec&,
159  ICARUSPeakParamsVec&,
160  double&,
161  int&, int) const;
162  void findLongPeakParameters(const std::vector<float>&,
163  const reco_tool::ICandidateHitFinder::HitCandidateVec&,
164  ICARUSPeakParamsVec&,
165  double&,
166  int&, int) const;
167  double ComputeChiSquare(TF1 func, TH1 *histo) const;
168  double ComputeNullChiSquare(std::vector<float>) const;
169 
170 
171  void setWire(int i) {
172  iWire=i;
173  // std::cout << " setting iwire " << iWire << std::endl;
174  } ;
175  private:
176  std::vector<double> localmeans;
177  unsigned int fDataSize; //SIZE OF RAW DATA ON ONE WIRE.
178  art::InputTag fDigitModuleLabel; //MODULE THAT MADE DIGITS.
179  std::string fSpillName; //NOMINAL SPILL IS AN EMPTY STRING.
180 
181  //FFT COPIED VARIABLES.
182  std::string fCalDataModuleLabel;
183  std::string fHitLabelName;
184 
185  int fThetaAngle;
186  bool fUncompressWithPed; //OPTION TO UNCOMPRESS WITH PEDESTAL.
187 
188  std::unique_ptr<reco_tool::ICandidateHitFinder> fHitFinderTool; ///< For finding candidate hits
189  // PeakFitterICARUS* fPeakFitterTool; ///< Perform fit to candidate peaks
190 
191  size_t fMaxMultiHit; ///<maximum hits for multi fit
192 double fChi2NDF; ///maximum Chisquared / NDF allowed for a hit to be saved
193  std::vector<int> fLongMaxHitsVec; ///<Maximum number hits on a really long pulse train
194  std::vector<int> fLongPulseWidthVec; ///<Sets width of hits used to describe long pulses
195 
196  //histograms
197  TH1F* fFirstChi2;
198  TH1F* fNullChi2;
199 
200  TH1F* fChi2;
201  TH1F* fHeightC;
202  TH1F* fWidthC;
203  TH1F* fNoiseC;
204  TH1F* fAreaC;
205  TH1F* fnhwC;
206  TH1F* fIntegralC;
207  TH1F* fAreaInt;
208  TH1F* fBaselineC;
209 
210  TH1F* fWire2771;
211 
212  TH1F* fHeightI2;
213  TH1F* fWidthI2;
214  TH1F* fNoiseI2;
215  TH1F* fHeightI1;
216  TH1F* fWidthI1;
217  TH1F* fNoiseI1;
218 
219  // Member variables from the fhicl file
220  double fMinWidth; ///< minimum initial width for ICARUS fit
221  double fMaxWidthMult; ///< multiplier for max width for ICARUS fit
222  int fFittingRange; ///< semi-width of interval where to fit hit
223  int fIntegratingRange; ///< semi-width of interval where to integrate fitting function
224 
225 
226  int iWire;
227 
228  mutable ICARUShitFitCache fFitCache; ///< Cached functions for multi-peak fits.
229  mutable ICARUSlongHitFitCache fLongFitCache; ///< Cached functions for long hits.
230 
231  const geo::GeometryCore* fGeometry = lar::providerFrom<geo::Geometry>();
232 
233  }; // class ICARUSHitFinder
234 
235 
236  //-------------------------------------------------
237  ICARUSHitFinder::ICARUSHitFinder(fhicl::ParameterSet const& pset) : EDProducer{pset}
238  {
239  this->reconfigure(pset);
240 
241  //LET HITCOLLECTIONCREATOR DECLARE THAT WE ARE GOING TO PRODUCE
242  //HITS AND ASSOCIATIONS TO RAW DIGITS BUT NOT ASSOCIATIONS TO WIRES
243  //(WITH NO PARTICULAR PRODUCT LABEL).
244  recob::HitCollectionCreator::declare_products(producesCollector(),
245  /*instance_name*/"");
246  }
247 
248  //-------------------------------------------------
249  void ICARUSHitFinder::reconfigure(fhicl::ParameterSet const& p)
250  {
251  fUncompressWithPed = p.get< bool >("UncompressWithPed", true);
252  fDigitModuleLabel = p.get< art::InputTag >("DigitModuleLabel", "daq");
253  fCalDataModuleLabel = p.get< std::string >("CalDataModuleLabel");
254  fMaxMultiHit = p.get< int >("MaxMultiHit");
255  fChi2NDF = p.get< double >("Chi2NDF");
256  fLongPulseWidthVec = p.get< std::vector<int>>("LongPulseWidth", std::vector<int>() = {16,16,16});
257  fLongMaxHitsVec = p.get< std::vector<int>>("LongMaxHits", std::vector<int>() = {25,25,25});
258  fThetaAngle=p.get< double >("ThetaAngle");
259  fMinWidth=p.get< double >("MinWidth");
260  fMaxWidthMult=p.get< double >("MaxWidthMult");
261  fFittingRange=p.get< int >("FittingRange");
262  fIntegratingRange=p.get< int >("IntegratingRange");
263 
264 
265  fHitFinderTool = art::make_tool<reco_tool::ICandidateHitFinder>(p.get<fhicl::ParameterSet>("CandidateHits"));
266  // fPeakFitterTool = art::make_tool<reco_tool::PeakFitterICARUS>(p.get<fhicl::ParameterSet>("PeakFitter"));
267 
268  mf::LogInfo("ICARUSHitFinder_module") << "fDigitModuleLabel: " << fDigitModuleLabel << std::endl;
269  }
270 
271  //-------------------------------------------------
272  void ICARUSHitFinder::beginJob()
273  {
274  // get access to the TFile service
275  art::ServiceHandle<art::TFileService> tfs;
276 
277 
278  // ======================================
279  // === Hit Information for Histograms ===
280  fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 1000);
281  fNullChi2 = tfs->make<TH1F>("fNullChi2", "#chi^{2}", 100, 0, 10);
282 
283  fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 1000);
284  fHeightC = tfs->make<TH1F>("fHeightC", "height(ADC#)", 100, 0, 100);
285  fWidthC = tfs->make<TH1F>("fWidthC", "width(samples)", 200, 0, 200);
286  fNoiseC = tfs->make<TH1F>("fNoiseC", "Noise Area(ADC#)", 100, 0, 100);
287  fAreaC = tfs->make<TH1F>("fAreaC", "", 150, 0, 1500);
288  fAreaC->GetXaxis()->SetTitle("Area(ADC#*ticks)");
289  fAreaC->GetYaxis()->SetTitle("events");
290 
291  fWire2771 = tfs->make<TH1F>("fWire2771", "", 4096, -0.5, 4095.5);
292  fWire2771->GetXaxis()->SetTitle("tick");
293  fWire2771->GetYaxis()->SetTitle("ADC#");
294 
295  fIntegralC = tfs->make<TH1F>("fIntegralC", "", 500, 0, 3000);
296  fIntegralC->GetXaxis()->SetTitle("Area(ADC#*ticks)");
297  fIntegralC->GetYaxis()->SetTitle("events");
298  fBaselineC = tfs->make<TH1F>("fBaselineC", "", 500, -10., 10.);
299  fBaselineC->GetXaxis()->SetTitle("Area(ADC#*ticks)");
300  fBaselineC->GetYaxis()->SetTitle("events");
301  fAreaInt = tfs->make<TH1F>("fAreaInt", "Area/Int", 100,0.,2.);
302 
303  fnhwC = tfs->make<TH1F>("fnhwC", "", 10, 0, 10);
304  fnhwC->GetXaxis()->SetTitle("n hits per wire");
305  fnhwC->GetYaxis()->SetTitle("wires");
306  fHeightI2 = tfs->make<TH1F>("fHeightI2", "height(ADC#)", 100, 0, 100);
307  fWidthI2 = tfs->make<TH1F>("fWidthI2", "width(samples)", 100, 0, 100);
308  fNoiseI2 = tfs->make<TH1F>("fNoiseI2", "Noise Area(ADC#)", 100, 0, 100);
309  fHeightI1 = tfs->make<TH1F>("fHeightI1", "height(ADC#)", 100, 0, 100);
310  fWidthI1 = tfs->make<TH1F>("fWidthI1", "width(samples)", 100, 0, 100);
311  fNoiseI1 = tfs->make<TH1F>("fNoiseI1", "Noise Area(ADC#)", 100, 0, 100);
312  // std::cout << " ICARUSHitfinder begin " << std::endl;
313  }
314 
315  void ICARUSHitFinder::endJob()
316  {
317  // std::cout << " ICARUSHitFinder endjob " << std::endl;
318 
319  }
320 
321  //-------------------------------------------------
322  void ICARUSHitFinder::produce(art::Event& evt)
323  {
324 
325  //0
326  //return;
327 
328  std::ofstream output("areaFit.out");
329 
330  // std::cout << " ICARUSHitFinder produce " << std::endl;
331 
332  //GET THE GEOMETRY.
333  art::ServiceHandle<geo::Geometry> geom;
334 
335  // ###############################################
336  // ### Making a ptr vector to put on the event ###
337  // ###############################################
338  // this contains the hit collection
339  // and its associations to wires and raw digits
340 
341  // Handle the filtered hits collection...
342  recob::HitCollectionCreator hcol(evt);
343 
344  // if (fAllHitsInstanceName != "") filteredHitCol = &hcol;
345 
346  // ##########################################
347  // ### Reading in the Wire List object(s) ###
348  // ##########################################
349  art::Handle< std::vector<recob::Wire> > wireVecHandle;
350  evt.getByLabel(fCalDataModuleLabel,wireVecHandle);
351 
352  // #################################################################
353  // ### Reading in the RawDigit associated with these wires, too ###
354  // #################################################################
355  art::FindOneP<raw::RawDigit> RawDigits
356  (wireVecHandle, evt, fCalDataModuleLabel);
357 
358  // Channel Number
359  raw::ChannelID_t channel = raw::InvalidChannelID;
360 
361 
362  std::vector<float> holder; //HOLDS SIGNAL DATA.
363  std::vector<short> rawadc; //UNCOMPRESSED ADC VALUES.
364 
365  std::vector<float> startTimes; //STORES TIME OF WINDOW START.
366  std::vector<float> maxTimes; //STORES TIME OF LOCAL MAXIMUM.
367  std::vector<float> endTimes; //STORES TIME OF WINDOW END.
368  std::vector<float> peakHeight; //STORES ADC COUNT AT THE MAXIMUM.
369  std::vector<float> hitrms; //STORES CHARGE WEIGHTED RMS OF TIME ACROSS THE HIT.
370  std::vector<double> charge; //STORES THE TOTAL CHARGE ASSOCIATED WITH THE HIT.
371 
372  // geo::SigType_t sigType = geo::kInduction;
373  std::stringstream numConv;
374 
375  unsigned int hwC(0),lwC(5728); //lowest and highest wire with physical deposition in Collection
376  unsigned int hwI2(0),lwI2(5728); //same for Induction
377  unsigned int hwI1(0),lwI1(2112);
378 
379  unsigned int nw1hitC(0),nw1hitI2(0),nw1hitI1(0); //number of wires (between lowest and highest) with a single hit (rough definition of hitfinding efficiency)
380  int nhWire[5600];
381  for(int jw=0;jw<5600;jw++)
382  nhWire[jw]=0;
383  float wCharge[5600];
384  for(int jw=0;jw<5600;jw++)
385  wCharge[jw]=0;
386  float wInt[5600];
387  for(int jw=0;jw<5600;jw++)
388  wInt[jw]=0;
389 
390 
391  unsigned int nhitsC(0),nhitsI1(0),nhitsI2(0); //total number of reconstructed hits in a view
392 
393  unsigned int nnhitsC(0),nnhitsI1(0),nnhitsI2(0); //total number of reconstructed hits in a view
394 
395 
396  unsigned int minWireC,maxWireC;
397  if(fThetaAngle==45) {minWireC=2539; maxWireC=3142;}
398  if(fThetaAngle==0) {minWireC=2535; maxWireC=4486;}
399  if(fThetaAngle==20) {minWireC=2539; maxWireC=4000;}
400  if(fThetaAngle==40) {minWireC=2539; maxWireC=3190;}
401  if(fThetaAngle==60) {minWireC=2539; maxWireC=2905;}
402  if(fThetaAngle==70) {minWireC=2539; maxWireC=2805;}
403  if(fThetaAngle==80) {minWireC=2539; maxWireC=2740;}
404 
405  //GET THE LIST OF BAD CHANNELS.
406  lariov::ChannelStatusProvider const& channelStatus
407  = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
408 
409  lariov::ChannelStatusProvider::ChannelSet_t const BadChannels
410  = channelStatus.BadChannels();
411 
412  unsigned int minWireI2=2539; //empirical
413  unsigned int maxWireI2=4700;
414  unsigned int minDrift=850;
415  unsigned int maxDrift=1500;
416 
417  //### Looping over the wires ###
418  //##############################
419  for(size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++)
420  {
421  // ####################################
422  // ### Getting this particular wire ###
423  // ####################################
424  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
425  art::Ptr<raw::RawDigit> rawdigits = RawDigits.at(wireIter);
426 
427  // --- Setting Channel Number and Signal type ---
428  channel = wire->Channel();
429 
430  std::vector<float> signal(wire->Signal());
431 
432 
433  // get the WireID for this hit
434  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
435  // for now, just take the first option returned from ChannelToWire
436  geo::WireID wid = wids[0];
437  // We need to know the plane to look up parameters
438  geo::PlaneID::PlaneID_t plane = wid.Plane;
439  size_t cryostat=wid.Cryostat;
440  size_t tpc=wid.TPC;
441  size_t iWire=wid.Wire;
442 
443 
444  holder.clear();
445  localmeans.clear();
446 
447  //GET THE REFERENCE TO THE CURRENT raw::RawDigit.
448  channel = rawdigits->Channel();
449  fDataSize = rawdigits->Samples();
450 
451  std::vector<geo::WireID> widVec = geom->ChannelToWire(channel);
452  size_t iwire = widVec[0].Wire;
453  // size_t plane = widVec[0].Plane;
454 
455 
456  rawadc.resize(fDataSize);
457  holder.resize(fDataSize);
458 
459  //UNCOMPRESS THE DATA.
460  if (fUncompressWithPed) {
461  int pedestal = (int)rawdigits->GetPedestal();
462  raw::Uncompress(rawdigits->ADCs(), rawadc, pedestal, rawdigits->Compression());
463  }
464  else{
465  raw::Uncompress(rawdigits->ADCs(), rawadc, rawdigits->Compression());
466  }
467 
468  mf::LogDebug("ICARUSHitFinder") << " pedestal " <<rawdigits->GetPedestal() << std::endl;
469 
470  for(unsigned int bin = 0; bin < fDataSize; ++bin){
471  //holder[bin]=(rawadc[bin]-rawdigits->GetPedestal());
472  holder[bin]=signal[bin];
473 
474  //std::cout << " bin " << bin << " rawadc " << rawadc[bin]-rawdigits->GetPedestal() << " signal " << signal[bin] << std::endl;
475  if(plane == 0) holder[bin]=-holder[bin];
476  //if(plane == 1) holder[bin]=-holder[bin];
477  // if(plane==2&&iwire==2600) std::cout << " wire 2600 bin " << bin << " signal " << holder[bin] << std::endl;
478  }
479 
480  if(plane==0&&iwire<lwI1) lwI1=iwire;
481  if(plane==0&&iwire>hwI1) hwI1=iwire;
482  if(plane==1&&iwire<lwI2) lwI2=iwire;
483  if(plane==1&&iwire>hwI2) hwI2=iwire;
484  if(plane==2&&iwire<lwC) lwC=iwire;
485  if(plane==2&&iwire>hwC) hwC=iwire;
486 
487 
488 
489  // sigType = geom->SignalType(channel);
490 
491  peakHeight.clear();
492  endTimes.clear();
493  startTimes.clear();
494  maxTimes.clear();
495  charge.clear();
496  hitrms.clear();
497 
498  bool channelSwitch = false;
499 
500  for(auto it = BadChannels.begin(); it != BadChannels.end(); it++)
501  {
502  if(channel==*it)
503  {
504  channelSwitch = true;
505  break;
506  }
507  }
508 
509  if(channelSwitch==false)
510  { //1
511 
512  // Hit finding parameters
513  double chargeErr(0); //CHI2/NDF and error on charge.
514  // double hrms(0);
515  //double totSig;
516 
517  double chi2null=ComputeNullChiSquare(holder);
518 // std::cout << " wire " << iWire << " chi2null " << chi2null << std::endl;
519  fNullChi2->Fill(chi2null);
520 
521 
522 
523  reco_tool::ICandidateHitFinder::HitCandidateVec hitCandidateVec;
524 
525  reco_tool::ICandidateHitFinder::MergeHitCandidateVec mergedCandidateHitVec;
526 
527  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
528 
529  std::vector<float> tempVec = holder;
530  recob::Wire::RegionsOfInterest_t::datarange_t rangeData(size_t(0),std::move(tempVec));
531 
532  fHitFinderTool->findHitCandidates(rangeData, 0,channel,0,hitCandidateVec);
533  //int jc=0;
534  for(auto& hitCand : hitCandidateVec) {
535  expandHit(hitCand,holder,hitCandidateVec);
536 
537  }
538 
539 
540 
541  fHitFinderTool->MergeHitCandidates(rangeData, hitCandidateVec, mergedCandidateHitVec);
542 
543  //numHits = hits.size();
544  int nghC=0;
545  int nghI2=0;
546  int nghI1=0;
547 
548  // if(plane==0)
549  // std::cout << " plane " << plane << " Wire " << iwire << " numhits " << hitCandidateVec.size() <<" mergedhits " << mergedCandidateHitVec.size() << std::endl;
550 
551 
552  //FIT ONLY COLLECTION HITS
553  if(plane==2) {
554  //std::cout << " mergedcands size " << mergedCandidateHitVec.size() << std::endl;
555 
556  for(auto& mergedCands : mergedCandidateHitVec)
557  {
558 
559 
560  int startT= mergedCands.front().startTick-fFittingRange;
561  int endT = mergedCands.back().stopTick+fFittingRange;
562  //std::cout << " fitting range " << fFittingRange << std::endl;
563 
564  float mean;
565  computeBestLocalMean(mergedCands,holder,mergedCandidateHitVec,mean);
566  localmeans.push_back(mean);
567  mf::LogDebug("ICARUSHitFinder") << " adding localmean " << mean << std::endl;
568  // ### Putting in a protection in case things went wrong ###
569  // ### In the end, this primarily catches the case where ###
570  // ### a fake pulse is at the start of the ROI ###
571  if (endT - startT < 5) continue;
572  fWidthC->Fill(endT-startT);
573  // #######################################################
574  // ### Clearing the parameter vector for the new pulse ###
575  // #######################################################
576 
577  // === Setting the number of Gaussians to try ===
578  int nGausForFit = mergedCands.size();
579 
580  // ##################################################
581  // ### Calling the function for fitting ICARUS ###
582  // ##################################################
583  double chi2PerNDF(0.);
584  double chi2Long(0.);
585 
586  int NDF(1);
587  const int npk=mergedCands.size();
588  ICARUSPeakParamsVec peakParamsVec(npk);
589  peakParamsVec.clear();
590  int islong=0;
591  if (mergedCands.size() <= fMaxMultiHit)
592  {
593 
594  // std::cout << "setting iwire " << iwire << std::endl;
595  // fPeakFitterTool->setWire(iwire);
596  // std::cout << " fitting iwire " << iwire << std::endl;
597  // std::cout << " cryostat " << cryostat << " tpc " << tpc << " plane " << plane << " wire " << iwire << std::endl;
598  findMultiPeakParameters(signal, mergedCands, peakParamsVec, chi2PerNDF, NDF, iwire);
599 
600  if (!(chi2PerNDF < std::numeric_limits<double>::infinity()))
601  {
602  chi2PerNDF = 200.;
603  NDF = 2;
604  }
605  fFirstChi2->Fill(chi2PerNDF);
606  // std::cout << " wire " << iwire << " first chi2NDF " << chi2PerNDF << std::endl;
607  }
608 
609  // std::cout << " before longpulse chi2 " << chi2PerNDF << " threshold " << fChi2NDF << std::endl;
610  if (chi2PerNDF < 10.)
611  fChi2->Fill(chi2PerNDF);
612  // if(chi2PerNDF<0.1) // change from 10 to reduce output
613  // std::cout << " wire " << iwire << " SMALL chi2NDF " << chi2PerNDF << " thr chi2NDF " << fChi2NDF << std::endl;
614  ICARUSPeakParamsVec peakParamsLong(npk);
615  peakParamsLong.clear();
616  if (chi2PerNDF > fChi2NDF)
617  {
618  islong=1;
619  findLongPeakParameters(signal, mergedCands, peakParamsLong, chi2Long, NDF, iwire);
620  // if(chi2Long<0.3) std::cout << " small chi2long " << chi2Long << std::endl;
621  if(chi2Long<chi2PerNDF&&chi2Long>0.1) {
622  fChi2->Fill(chi2Long);
623  peakParamsVec=peakParamsLong;
624  }
625  else { fChi2->Fill(chi2PerNDF);
626  // if(chi2PerNDF<0.1) // change from 10 to reduce output
627  // std::cout << " wire " << iwire << " SMALLSMALL chi2NDF " << chi2PerNDF << " thr chi2NDF " << fChi2NDF << std::endl;
628  }
629  }
630 
631  // unsigned int jhit=0;
632 
633  //std::cout << " before peak loop" << std::endl;
634  // for(const auto& peakParams : peakParamsVec)
635 for(unsigned int jhit=0;jhit<mergedCands.size(); jhit++)
636  {
637  //float fitCharge=chargeFunc(peakMean, peakAmp, peakWidth, fAreaNormsVec[plane],startT,endT);
638  //float fitChargeErr = std::sqrt(TMath::Pi()) * (peakAmpErr*peakWidthErr + peakWidthErr*peakAmpErr);
639  unsigned int startInt=mergedCands[jhit].startTick-fIntegratingRange;
640  unsigned int endInt=mergedCands[jhit].stopTick+fIntegratingRange;
641 
642  if(jhit>=1&&startInt<mergedCands[jhit-1].stopTick) startInt=mergedCands[jhit-1].stopTick;
643  if(jhit<mergedCands.size()-1&&endInt>mergedCands[jhit+1].startTick) endInt=mergedCands[jhit+1].startTick;
644 
645  float fitCharge=0;
646  float peakAmp, peakMean, peakLeft, peakRight, peakBaseline;
647  float peakSlope=0, peakFitWidth=0;
648  float peakMeanErr, peakAmpErr;
649  if(!islong) {
650  // TF1 Func("ICARUSfunc",fitf,start,end,1+5*mergedCands.size());
651  TF1& Func = *(fFitCache.Get(mergedCands.size()));
652  assert(&Func);
653  Func.SetParameter(0, mergedCands.size());
654 
655  float intBaseline=0;
656 
657  ICARUSPeakFitParams_t peakParams=peakParamsVec[jhit];
658 
659  // Extract values for this FITTED peak
660  peakAmp = peakParams.peakAmplitude;
661  peakMean = peakParams.peakCenter;
662  //float peakWidth = peakParams.peakSigma;
663  peakLeft = peakParams.peakTauLeft;
664  peakRight = peakParams.peakTauRight;
665  peakBaseline = peakParams.peakBaseline;
666 
667 
668  // Place one bit of protection here
669  if (std::isnan(peakAmp))
670  {
671  // std::cout << "**** hit peak amplitude is a nan! Channel: " << channel << ", start tick: " << startT << std::endl;
672  continue;
673  }
674 
675  // Extract errors
676  peakAmpErr = peakParams.peakAmplitudeError;
677  peakMeanErr = peakParams.peakCenterError;
678  // float peakWidthErr = peakParams.peakSigmaError;
679  Func.SetParameter(1+5*jhit,peakBaseline);
680  Func.SetParameter(2+5*jhit,peakAmp);
681  Func.SetParameter(3+5*jhit,peakMean);
682  Func.SetParameter(4+5*jhit,peakRight);
683  Func.SetParameter(5+5*jhit,peakLeft);
684 
685  intBaseline+=(endInt-startInt)*peakBaseline;
686 
687  try
688  {
689  fitCharge=Func.Integral(startInt,endInt)-(endInt-startInt)*localmeans[jhit];
690 
691  }
692  catch(...) {
693  mf::LogWarning("ICARUSHitFinder") << "Icarus numerical 32 failed";
694  fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.)-(endInt-startInt)*localmeans[jhit];
695  }
696  }
697  else {
698  // TF1 FuncLong("ICARUSfuncLong",fitlong,start,end,1+7*mergedCands.size());
699  TF1& FuncLong = *(fLongFitCache.Get(mergedCands.size()));
700  assert(&FuncLong);
701  FuncLong.SetParameter(0, mergedCands.size());
702 
703 
704  float intBaseline=0;
705 
706  ICARUSPeakFitParams_t peakParams=peakParamsVec[jhit];
707 
708  // Extract values for this FITTED peak
709  peakAmp = peakParams.peakAmplitude;
710  peakMean = peakParams.peakCenter;
711  //float peakWidth = peakParams.peakSigma;
712  peakLeft = peakParams.peakTauLeft;
713  peakRight = peakParams.peakTauRight;
714  peakBaseline = peakParams.peakBaseline;
715  intBaseline+=(endInt-startInt)*peakBaseline;
716  peakAmpErr = peakParams.peakAmplitudeError;
717  peakMeanErr = peakParams.peakCenterError;
718 
719  FuncLong.SetParameter(1+7*jhit,peakBaseline);
720  FuncLong.SetParameter(2+7*jhit,peakAmp);
721  FuncLong.SetParameter(3+7*jhit,peakMean);
722  FuncLong.SetParameter(4+7*jhit,peakRight);
723  FuncLong.SetParameter(5+7*jhit,peakLeft);
724  FuncLong.SetParameter(6+7*jhit,peakFitWidth);
725  FuncLong.SetParameter(7+7*jhit,peakSlope);
726 
727 
728  try
729  { fitCharge=FuncLong.Integral(startInt,endInt)-(endInt-startInt)*localmeans[jhit];
730  //fitCharge=FuncLong.Integral(start-35,end+35);
731  }
732  catch(...)
733  {mf::LogWarning("ICARUSHitFinder") << "Icarus numerical integration failed";
734  fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.)-(endInt-startInt)*localmeans[jhit];
735  }
736  }
737  if(isnan(fitCharge)&&!islong) fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.);
738  if(isnan(fitCharge)&&islong) fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.);
739  //Func.Integral(start,end);
740 
741  // float totSig20=std::accumulate(holder.begin() + (int) start-35, holder.begin() + (int) end+35, 0.);
742 
743  float totSig=std::accumulate(holder.begin()+ (int) startInt, holder.begin()+ (int) endInt, 0.)-(endInt-startInt)*localmeans[jhit];
744  fBaselineC->Fill(localmeans[jhit]);
745  // float fitChargeErr=0;
746  // if(plane==2&&iWire==2842)
747 
748 
749 
750 //std::cout << " before hit creator " << std::endl;
751  recob::HitCreator hit(
752  *wire, //RAW DIGIT REFERENCE.
753  wid, //WIRE ID.
754  startInt, //START TICK.
755  endInt, //END TICK.
756  (peakLeft+peakRight)/2., //RMS.
757  peakMean, //PEAK_TIME.
758  peakMeanErr, //SIGMA_PEAK_TIME.
759  peakAmp, //PEAK_AMPLITUDE.
760  peakAmpErr, //SIGMA_PEAK_AMPLITUDE.
761  fitCharge, //HIT_INTEGRAL.
762  chargeErr, //HIT_SIGMA_INTEGRAL.
763  totSig, //SUMMED CHARGE.
764  nGausForFit, //MULTIPLICITY.
765  jhit, //LOCAL_INDEX.
766  chi2PerNDF, //WIRE ID.
767  NDF //DEGREES OF FREEDOM.
768  );
769 
770  mf::LogDebug("ICARUSHitFinder") << " fitcharge " << fitCharge << " totSig " << totSig << std::endl;
771  //std::cout << " summedADC " << hit.summedADC() << " integral " << hit.Integral() << std::endl;
772 
773  // filteredHitVec.push_back(hit.copy());
774  hcol.emplace_back(hit.move(), wire, rawdigits);
775 
776  bool wireWindowC=iWire>=minWireC&&iWire<=maxWireC;
777  bool outWireWindowC=iWire<=minWireC||iWire>=maxWireC;
778 
779  if(wireWindowC) {
780  nghC++;
781  if(nghC==1) nw1hitC++;
782  nhWire[iWire]++;
783  fHeightC->Fill(peakAmp);
784  // fWidthC->Fill(peakWidth);
785  // fAreaC->Fill(totSig);
786  wCharge[iWire]+=fitCharge;
787  wInt[iWire]+=totSig;
788  }
789 
790 
791  nhitsC++;
792 
793 
794  if(tpc==0&&cryostat==0&&outWireWindowC) {
795  nnhitsC++; fNoiseC->Fill(peakAmp); }
796 
797 
798 
799  if(cryostat==0&&tpc==0)
800  if(wCharge[iwire]>0)
801  fAreaC->Fill(wCharge[iwire]);
802  if(cryostat==0&&tpc==0)
803  if(wInt[iwire]>0)
804  fIntegralC->Fill(wInt[iwire]);
805  if(cryostat==0&&tpc==0)
806  if(wCharge[iwire]>0)
807  output << iwire << " " <<wInt[iwire] << std::endl;
808  if(wInt[iwire]>0&&cryostat==0&&tpc==0)
809  fAreaInt->Fill(wCharge[iwire]/wInt[iwire]);
810 
811  } // loop on peakparams vector
812 
813  } // loop on merged hits
814 
815  // std::cout << " after coll hit " << std::endl;
816  } //COLLECTION
817 
818  if(plane==0||plane==1) {
819  for(auto& mergedCands : mergedCandidateHitVec)
820  {
821  for(size_t jh=0;jh<mergedCands.size();jh++) {
822  // std::cout << " plane " << plane << " wire " << iwire << " hit " <<mergedCands[jh].hitCenter << std::endl;
823  //FOR INDUCTION HITS STORE RAW INFORMATION
824  recob::HitCreator hit(
825  *wire, //RAW DIGIT REFERENCE.
826  wid, //WIRE ID.
827  mergedCands[jh].startTick, //START TICK.
828  mergedCands[jh].stopTick, //END TICK.
829  mergedCands[jh].hitSigma, //RMS.
830  mergedCands[jh].hitCenter, //PEAK_TIME.
831  0, //SIGMA_PEAK_TIME.
832  mergedCands[jh].hitHeight, //PEAK_AMPLITUDE.
833  0, //SIGMA_PEAK_AMPLITUDE.
834  std::accumulate(holder.begin() + (int) mergedCands[jh].startTick, holder.begin() + (int) mergedCands[jh].stopTick, 0.), //HIT_INTEGRAL.
835  0, //HIT_SIGMA_INTEGRAL.
836  std::accumulate(holder.begin() + (int) mergedCands[jh].startTick, holder.begin() + (int) mergedCands[jh].stopTick, 0.), //SUMMED CHARGE.
837  1, //MULTIPLICITY.
838  jh, //LOCAL_INDEX.
839  0, //WIRE ID.
840  int(mergedCands[jh].stopTick-mergedCands[jh].startTick+1) //DEGREES OF FREEDOM.
841  );
842  hcol.emplace_back(hit.move(), wire, rawdigits);
843 
844  bool driftWindow=(mergedCands[jh].hitCenter)>=minDrift&&(mergedCands[jh].hitCenter)<=maxDrift;
845  bool wireWindowI2=iWire>=minWireI2&&iWire<=maxWireI2;
846  bool outDriftWindow=(mergedCands[jh].hitCenter)<=minDrift||(mergedCands[jh].hitCenter)>=maxDrift;
847  bool outWireWindowI2=iWire<=minWireI2||iWire>=maxWireI2;
848 
849 
850  if(plane==0&&driftWindow) {
851  // std::cout << " wire " << hits[i].iWire << " ngh " << ngh << std::endl;
852  nghI1++;
853  if(nghI1==1) nw1hitI1++;
854 
855  // std::cout << " wire " << hits[i].iWire << " nw1h " << nw1 << std::endl;
856  fHeightI1->Fill(mergedCands[jh].hitHeight);
857  fWidthI1->Fill(mergedCands[jh].hitSigma);
858  }
859  if(plane==1&&wireWindowI2) {
860  // std::cout << " wire " << hits[i].iWire << " ngh " << ngh << std::endl;
861  nghI2++;
862  if(nghI2==1) nw1hitI2++;
863  // std::cout << " filling height histo wire" << hits[i].iWire << " drift " << hits[i].hitCenter << " amplitude " << amplitude << std::endl;
864  fHeightI2->Fill(mergedCands[jh].hitHeight);
865  fWidthI2->Fill(mergedCands[jh].hitSigma);
866  }
867 
868 
869  if(plane==0) nhitsI1++;
870  if(plane==1) nhitsI2++;
871 
872  if(plane==0&&tpc==0&&cryostat==0&&outDriftWindow) {nnhitsI1++; fNoiseI1->Fill(mergedCands[jh].hitHeight); }
873  if(plane==1&&tpc==0&&cryostat==0&&outWireWindowI2) { nnhitsI2++; fNoiseI2->Fill(mergedCands[jh].hitHeight); }
874  }} // merged loop
875  } //INDUCTION
876 
877 
878 
879  } //end channel condition
880 
881  } //end loop on channels
882 
883 // std::cout << " nhitsI1 " << nhitsI1 <<" nhitsI2 " << nhitsI2 <<" nhitsC " << nhitsC << std::endl;
884 
885 
886  for(unsigned int jw=minWireC;jw<maxWireC;jw++)
887  fnhwC->Fill(nhWire[jw]);
888 
889 
890  hcol.put_into(evt);
891  //std::cout << " end ICARUSHitfinder " << std::endl;
892 
893 
894  } //end produce
895 
896 void ICARUSHitFinder::expandHit(reco_tool::ICandidateHitFinder::HitCandidate& h, std::vector<float> holder, std::vector<reco_tool::ICandidateHitFinder::HitCandidate> how)
897  {
898  // Given a hit or hit candidate <hit> expand its limits to the closest minima
899  int nsamp=50;
900  int cut=1;
901  int upordown;
902  int found=0;
903 
904  unsigned int first=h.startTick;
905  unsigned int last =h.stopTick;
906 
907  std::vector<reco_tool::ICandidateHitFinder::HitCandidate>::iterator hiter;
908  std::vector<reco_tool::ICandidateHitFinder::HitCandidate> hlist;
909 
910  // fill list of existing hits on this wire
911  for(unsigned int j=0;j<how.size();j++)
912  {
913  reco_tool::ICandidateHitFinder::HitCandidate h2=how[j];
914  if(h2.hitCenter!=h.hitCenter)
915  hlist.push_back(h2);
916  }
917 
918 
919  // look for first sample
920  while(!found)
921  {
922  if(first==0)
923  break;
924 
925  for(hiter=hlist.begin();hiter!=hlist.end();hiter++)
926  {
927  reco_tool::ICandidateHitFinder::HitCandidate h2=*hiter;
928  if(first==h2.stopTick)
929  found=1;
930  }
931 
932  if(found==1) break;
933 
934  upordown=0;
935  for(int l=0;l<nsamp/2;l++)
936  {
937  if(first-nsamp/2+l<4095) {
938  if(holder[first-nsamp/2+l+1]-holder[first-nsamp/2+l]>0) upordown++;
939  else if(holder[first-nsamp/2+l+1]-holder[first-nsamp/2+l]<0) upordown--;
940  } }
941  //std::cout << " checking " << first << " upordown " << upordown << std::endl;
942  if(upordown>cut)
943  first--;
944  else
945  found=1;
946  }
947 
948  // look for last sample
949  found=0;
950  while(!found)
951  {
952  if(last==4095)
953  {
954  found=1;
955  break;
956  }
957 
958  for(hiter=hlist.begin();hiter!=hlist.end();hiter++)
959  {
960  reco_tool::ICandidateHitFinder::HitCandidate h2=*hiter;
961  if(last==h2.startTick)
962  found=1;
963  }
964 
965  if(found==1) break;
966 
967  upordown=0;
968  for(int l=0;l<nsamp/2;l++)
969  {
970  if(last+nsamp/2-l>0) {
971  if(holder[last+nsamp/2-l]-holder[last+nsamp/2-l-1]>0) upordown++;
972  else if(holder[last+nsamp/2-l]-holder[last+nsamp/2-l-1]<0) upordown--;
973  } }
974 
975  if(upordown<-cut)
976  last++;
977  else
978  found=1;
979  }
980 
981  h.startTick=first;
982  h.stopTick=last;
983  }
984  void ICARUSHitFinder::computeBestLocalMean(std::vector<reco_tool::ICandidateHitFinder::HitCandidate> h, std::vector<float> holder, reco_tool::ICandidateHitFinder::MergeHitCandidateVec how, float& localmean)
985  {
986  const int bigw=130; //size of the window where to look for the minimum localmean value
987  const int meanw=70; //size of the window where the mean is calculated
988  const int outofbounds=9999;
989 
990  float samples1[bigw]; //list to contain samples bellow the startTick
991  float samples2[bigw]; //list to contain samples above the stopTick
992 
993  reco_tool::ICandidateHitFinder::MergeHitCandidateVec hlist;
994 
995  float min1;
996  float min2;
997  int startTick=h.front().startTick;
998  int stopTick=h.back().stopTick;
999  float count1=0;
1000  float count2=0;
1001  unsigned int shift1=0;
1002  unsigned int shift2=0;
1003  int foundborder1=0;
1004  int foundborder2=0;
1005 
1006  // fill list of existing hits on this wire
1007  for(unsigned int j=0;j<how.size();j++)
1008  {
1009  std::vector<reco_tool::ICandidateHitFinder::HitCandidate> h2=how[j];
1010  if(h2.front().startTick!=h.front().startTick)
1011  hlist.push_back(h2);
1012  }
1013 
1014  // fill the arrays of samples to be examined
1015  for(unsigned int i=0;i<bigw;i++)
1016  {
1017  // remove samples from other hits in the wire
1018  for(unsigned int j=0;j<hlist.size();j++)
1019  {
1020  std::vector<reco_tool::ICandidateHitFinder::HitCandidate> h2=hlist[j];
1021  if(startTick-i-shift1 >= h2.front().startTick && startTick-i-shift1 <= h2.back().stopTick)
1022  shift1+=h2.back().stopTick-h2.front().startTick+1;
1023  else if(stopTick+i+shift2 >= h2.front().startTick && stopTick+i+shift2 <= h2.back().stopTick)
1024  shift2+=h2.back().stopTick-h2.front().startTick+1;
1025  }
1026 
1027  // fill the lists
1028  //if(startTick-i-shift1>=0)
1029  samples1[i]=holder[h.front().startTick-i-shift1];
1030  //else
1031  //samples1[i]=outofbounds;
1032 
1033  if(stopTick+i+shift2<=4095)
1034  samples2[i]=holder[h.back().stopTick+i+shift2];
1035  else
1036  samples2[i]=outofbounds;
1037  }
1038 
1039  // initialize counters
1040  for(int j=0;j<meanw;j++)
1041  {
1042  if(samples1[j]==outofbounds) foundborder1=1;
1043  if(samples2[j]==outofbounds) foundborder2=1;
1044 
1045  if(!foundborder1) count1+=samples1[j];
1046  if(!foundborder2) count2+=samples2[j];
1047  }
1048  min1=count1;
1049  min2=count2;
1050 
1051  //look for the local minima from the filled lists
1052  for(int j=0;j<bigw-meanw;j++)
1053  {
1054  if(count1<min1) min1=count1;
1055  if(samples1[j+meanw]==outofbounds) foundborder1=1;
1056  if(!foundborder1) count1+=samples1[j+meanw]-samples1[j];
1057 
1058  if(count2<min2) min2=count2;
1059  if(samples2[j+meanw]==outofbounds) foundborder2=1;
1060  if(!foundborder2) count2+=samples2[j+meanw]-samples2[j];
1061 
1062  if(foundborder1 && foundborder2) break;
1063  }
1064  if(count1<min1) min1=count1;
1065  if(count2<min2) min2=count2;
1066 
1067  // std::cout << " localmean count1 " << count1 << " count2 " << count2 << std::endl;
1068  //for the moment take the highest value
1069  if((foundborder1 && !foundborder2) || (shift1 && !shift2))
1070  localmean=((float) min2)/meanw;
1071  else if((!foundborder1 && foundborder2) || (!shift1 && shift2))
1072  localmean=((float) min1)/meanw;
1073  else
1074  localmean=(min1>min2)?((float) min1)/meanw :((float) min2)/meanw;
1075 
1076  //for the moment take the average value
1077  // h->localmean=(min1+min2)/(2.*meanw);
1078  }
1079  void ICARUSHitFinder::findMultiPeakParameters(const std::vector<float>& roiSignalVec,
1080  const reco_tool::ICandidateHitFinder::HitCandidateVec& hitCandidateVec,
1081  ICARUSPeakParamsVec& peakParamsVec,
1082  double& chi2PerNDF,
1083  int& NDF, int iWire) const
1084  {
1085 
1086  ICARUSPeakParamsVec peakParamsVec0;
1087  TH1F* fHistogram=new TH1F("","",roiSignalVec.size(),0.,roiSignalVec.size());;
1088 
1089  std::string wireName = "PeakFitterHitSignal_" + std::to_string(iWire);
1090  fHistogram->SetName(wireName.c_str());
1091 
1092  if (hitCandidateVec.empty()) return;
1093 
1094  // in case of a fit failure, set the chi-square to infinity
1095  chi2PerNDF = std::numeric_limits<double>::infinity();
1096 
1097  int startTime = hitCandidateVec.front().startTick-fFittingRange;
1098  int endTime = hitCandidateVec.back().stopTick+fFittingRange;
1099  if(startTime<0) startTime=0;
1100  if(endTime>4095) endTime=4095;
1101  int roiSize = endTime - startTime;
1102 
1103  //std::cout << " roisize " << roiSize << std::endl;
1104 
1105  // Check to see if we need a bigger hianchstogram for fitting
1106  if (roiSize > fHistogram->GetNbinsX())
1107  {
1108  std::string histName = "PeakFitterHitSignal_" + std::to_string(iWire);
1109  fHistogram = new TH1F(histName.c_str(),"",roiSize,0.,roiSize);
1110  //fHistogram->Sumw2();
1111  }
1112 
1113  fHistogram->Reset();
1114  for(int idx = 0; idx < roiSize; idx++)
1115  fHistogram->SetBinContent(idx+1,roiSignalVec.at(startTime+idx));
1116  for(int idx = 0; idx < roiSize; idx++)
1117  fHistogram->SetBinError(idx+1,2.4);
1118  // for(int idx = 0; idx < roiSize; idx++)
1119  // std::cout << " bin " << idx << " Error " << fHistogram->GetBinError(idx+1) << std::endl;
1120 
1121 
1122  // Now define the complete function to fit
1123  // TF1 Func("ICARUSfunc",fitf,0,roiSize,1+5*hitCandidateVec.size());
1124  TF1& Func = *(fFitCache.Get(hitCandidateVec.size()));
1125  assert(&Func);
1126 
1127  // ### Setting the parameters for the ICARUS Fit ###
1128  Func.FixParameter(0, hitCandidateVec.size());
1129 
1130  int parIdx{0};
1131  for(auto const& candidateHit : hitCandidateVec)
1132  {
1133  double const peakMean = candidateHit.hitCenter - float(startTime);
1134  double const peakWidth = candidateHit.hitSigma;
1135  // std::cout << " peakWidth " << peakWidth << std::endl;
1136  // std::cout << " hitcenter " << candidateHit.hitCenter << " starttime " << startTime <<std::endl;
1137 
1138 
1139  double const amplitude = candidateHit.hitHeight;
1140  // double meanLowLim = std::max(peakMean - fPeakRange * peakWidth, 0.);
1141  // double meanHiLim = std::min(peakMean + fPeakRange * peakWidth, double(roiSize));
1142 
1143  Func.SetParameter(1+parIdx,0);
1144  Func.SetParameter(2+parIdx, amplitude);
1145  Func.SetParameter(3+parIdx, peakMean);
1146  Func.SetParameter(4+parIdx,peakWidth);
1147  Func.SetParameter(5+parIdx,peakWidth);
1148 
1149  Func.SetParLimits(1+parIdx, -5, 5);
1150  Func.SetParLimits(2+parIdx, 0.1 * amplitude, 10. * amplitude);
1151  Func.SetParLimits(3+parIdx, peakMean-peakWidth,peakMean+peakWidth);
1152  Func.SetParLimits(4+parIdx, std::max(fMinWidth, 0.01 * peakWidth), fMaxWidthMult * peakWidth);
1153  Func.SetParLimits(5+parIdx, std::max(fMinWidth, 0.01 * peakWidth), fMaxWidthMult * peakWidth);
1154 
1155  parIdx += 5;
1156 
1157  }
1158 
1159  int fitResult(-1);
1160  // if(hitCandidateVec.size()>2) return;
1161  try
1162  { fitResult = fHistogram->Fit(&Func,"QNWB","", 0., roiSize);
1163  }
1164  catch(...)
1165  {mf::LogWarning("GausHitFinder") << "Fitter failed finding a hit";}
1166 
1167 
1168  // if(fitResult==0)
1169  // std::cout << " icarus fit converges " << iWire << std::endl;
1170  if(fitResult!=0)
1171 // std::cout << " icarus fit cannot converge " << iWire << std::endl;
1172  // ##################################################
1173  // ### Getting the fitted parameters from the fit ###
1174  // ##################################################
1175  NDF = roiSize-5*hitCandidateVec.size();
1176  chi2PerNDF = (Func.GetChisquare() / NDF);
1177 
1178  double chi2mio=ComputeChiSquare(Func,fHistogram);
1179 // std::cout << " chi2mio " << chi2mio << std::endl;
1180  chi2PerNDF=chi2mio;
1181  // for(int idx = 0; idx < roiSize; idx++)
1182  // std::cout << " bin " << idx << " Error " << fHistogram->GetBinError(idx+1) << std::endl;
1183 
1184 // std::cout << " chi2 " << Func.GetChisquare() << std::endl;
1185 // std::cout << " ndf " << NDF << std::endl;
1186 
1187  // std::cout << " chi2ndf " << chi2PerNDF<< std::endl;
1188  parIdx = 0;
1189  peakParamsVec0.clear();
1190 
1191  for(size_t idx = 0; idx < hitCandidateVec.size(); idx++)
1192  {
1193  ICARUSPeakFitParams_t peakParams;
1194 
1195  peakParams.peakAmplitude = Func.GetParameter(2+parIdx);
1196  peakParams.peakAmplitudeError = Func.GetParError(2+parIdx);
1197  peakParams.peakCenter = Func.GetParameter(3+parIdx) + float(startTime);
1198  peakParams.peakCenterError = Func.GetParError(3+parIdx);
1199  //std::cout << " rising time " << Func.GetParameter(3) << " falling time " <<Func.GetParameter(4) << std::endl;
1200  peakParams.peakTauRight = Func.GetParameter(4+parIdx);
1201  peakParams.peakTauRightError = Func.GetParError(4+parIdx);
1202  peakParams.peakTauLeft = Func.GetParameter(5+parIdx);
1203  peakParams.peakTauLeftError = Func.GetParError(5+parIdx);
1204  peakParams.peakBaseline = Func.GetParameter(1+parIdx);
1205  peakParams.peakBaselineError = Func.GetParError(1+parIdx);
1206  peakParams.peakFitWidth =0;
1207  peakParams.peakFitWidthError = 0;
1208  peakParams.peakSlope = 0;
1209  peakParams.peakSlopeError = 0;
1210  peakParamsVec.emplace_back(peakParams);
1211  parIdx += 5;
1212  //std::cout << " first center " << Func.GetParameter(2) + float(startTime) << std::endl;
1213  // std::cout << " second center " << Func.GetParameter(7) + float(startTime) << std::endl;
1214  }
1215 
1216  //Gaus.Delete();
1217  //Func.Delete();
1218  bool writeWaveform=false;
1219  if(writeWaveform) {
1220  TFile *f = new TFile("fitICARUS.root","UPDATE");
1221  fHistogram->Write();
1222  f->Close();
1223  f->Delete();
1224  }
1225 
1226  fHistogram->Delete();
1227  return;
1228  }
1229 
1230  void ICARUSHitFinder::findLongPeakParameters(const std::vector<float>& roiSignalVec,
1231  const reco_tool::ICandidateHitFinder::HitCandidateVec& hitCandidateVec,
1232  ICARUSPeakParamsVec& peakParamsVec,
1233  double& chi2PerNDF,
1234  int& NDF, int iWire) const
1235  {
1236  TH1F* fHistogram=new TH1F("","",roiSignalVec.size(),0.,roiSignalVec.size());;
1237  std::string wireName = "PeakFitterHitSignal_" + std::to_string(iWire);
1238  fHistogram->SetName(wireName.c_str());
1239 
1240  if (hitCandidateVec.empty()) return;
1241 
1242  // in case of a fit failure, set the chi-square to infinity
1243  chi2PerNDF = std::numeric_limits<double>::infinity();
1244 
1245  int startTime = hitCandidateVec.front().startTick-fFittingRange;
1246  int endTime = hitCandidateVec.back().stopTick+fFittingRange;
1247  if(startTime<0) startTime=0;
1248  if(endTime>4095) endTime=4095;
1249 
1250  int roiSize = endTime - startTime;
1251 
1252  //std::cout << " roisize " << roiSize << std::endl;
1253 
1254  // Check to see if we need a bigger histogram for fitting
1255  if (roiSize > fHistogram->GetNbinsX())
1256  {
1257  std::string histName = "PeakFitterHitSignal_" + std::to_string(iWire);
1258  fHistogram = new TH1F(histName.c_str(),"",roiSize,0.,roiSize);
1259  fHistogram->Sumw2();
1260  }
1261 
1262  fHistogram->Reset();
1263  for(int idx = 0; idx < roiSize; idx++)
1264  fHistogram->SetBinContent(idx+1,roiSignalVec.at(startTime+idx));
1265 
1266  // Build the string to describe the fit formula
1267  std::string equation = "gaus(0)";
1268 
1269 
1270 
1271  // Now define the complete function to fit
1272  // TF1 Func("ICARUSfunc",fitlong,0,roiSize,1+7*hitCandidateVec.size());
1273  TF1& Func = *(fLongFitCache.Get(hitCandidateVec.size()));
1274  assert(&Func);
1275 
1276  // ### Setting the parameters for the ICARUS Fit ###
1277  Func.FixParameter(0,hitCandidateVec.size());
1278 
1279  int parIdx { 0 };
1280  for(auto const& candidateHit : hitCandidateVec)
1281  {
1282  double const peakMean = candidateHit.hitCenter - float(startTime);
1283  double const peakWidth = candidateHit.hitSigma;
1284  double const amplitude = candidateHit.hitHeight;
1285 
1286  Func.SetParameter(1+parIdx,0);
1287  Func.SetParameter(2+parIdx, amplitude);
1288  Func.SetParameter(3+parIdx, peakMean);
1289  Func.SetParameter(4+parIdx,peakWidth);
1290  Func.SetParameter(5+parIdx,peakWidth);
1291  Func.SetParameter(6+parIdx,2*peakWidth);
1292  Func.SetParameter(7+parIdx,0);
1293 
1294 
1295  Func.SetParLimits(1+parIdx, -5, 5);
1296  Func.SetParLimits(2+parIdx, 0.1 * amplitude, 10. * amplitude);
1297  Func.SetParLimits(3+parIdx, peakMean-peakWidth,peakMean+peakWidth);
1298  Func.SetParLimits(4+parIdx, std::max(fMinWidth, 0.01 * peakWidth), fMaxWidthMult * peakWidth);
1299  Func.SetParLimits(5+parIdx, std::max(fMinWidth, 0.01 * peakWidth), 4 * peakWidth);
1300  Func.SetParLimits(6+parIdx, 0,4*peakWidth);
1301  Func.SetParLimits(7+parIdx, -1,1);
1302 
1303  parIdx += 7;
1304 
1305  }
1306  int fitResult { -1 };
1307  try
1308  { fitResult = fHistogram->Fit(&Func,"QNWB","", 0., roiSize);
1309  }
1310  catch(...)
1311  {mf::LogWarning("GausHitFinder") << "Fitter failed finding a hit";}
1312 
1313  if(fitResult < -1)
1314  std::cout << " long fit cannot converge " << iWire << std::endl;
1315  // ##################################################
1316  // ### Getting the fitted parameters from the fit ###
1317  // ##################################################
1318  NDF = roiSize-7*hitCandidateVec.size();
1319  chi2PerNDF = (Func.GetChisquare() / NDF);
1320 
1321  parIdx = 0;
1322  peakParamsVec.clear();
1323  for(size_t idx = 0; idx < hitCandidateVec.size(); idx++)
1324  {
1325  ICARUSPeakFitParams_t peakParams;
1326 
1327  peakParams.peakAmplitude = Func.GetParameter(2+parIdx);
1328  peakParams.peakAmplitudeError = Func.GetParError(2+parIdx);
1329  peakParams.peakCenter = Func.GetParameter(3+parIdx) + float(startTime);
1330  peakParams.peakCenterError = Func.GetParError(3+parIdx);
1331 
1332  peakParams.peakTauRight = Func.GetParameter(4+parIdx);
1333  peakParams.peakTauRightError = Func.GetParError(4+parIdx);
1334  peakParams.peakTauLeft = Func.GetParameter(5+parIdx);
1335  peakParams.peakTauLeftError = Func.GetParError(5+parIdx);
1336  peakParams.peakFitWidth = Func.GetParameter(6+parIdx);
1337  peakParams.peakFitWidthError = Func.GetParError(6+parIdx);
1338  peakParams.peakSlope = Func.GetParameter(7+parIdx);
1339  peakParams.peakSlopeError = Func.GetParError(7+parIdx);
1340  peakParams.peakBaseline = Func.GetParameter(1+parIdx);
1341  peakParams.peakBaselineError = Func.GetParError(1+parIdx);
1342  // std::cout << " before adding peakparams size " << peakParamsVec.size() << std::endl;
1343  peakParamsVec.emplace_back(peakParams);
1344  // std::cout << " after adding peakparams size " << peakParamsVec.size() << std::endl;
1345 
1346  parIdx += 7;
1347 
1348  }
1349 
1350  //Func.Delete();
1351  bool writeWaveform=false;
1352  if(writeWaveform) {
1353  TFile *f = new TFile("fitICARUS.root","UPDATE");
1354  fHistogram->Write();
1355  f->Close();
1356  f->Delete();
1357  }
1358 
1359  fHistogram->Delete();
1360  return;
1361  }
1362 
1363 
1364 Double_t ICARUShitFitCache::fitf(Double_t const* x, Double_t const* par)
1365  {
1366  int const npeaks=(int)(par[0]);
1367  Double_t fitval=0;
1368  for(int jp=0;jp<npeaks;jp++)
1369  fitval += par[5*jp+1]+par[5*jp+2]*TMath::Exp(-(x[0]-par[5*jp+3])/par[5*jp+4])/(1+TMath::Exp(-(x[0]-par[5*jp+3])/par[5*jp+5]));
1370  return fitval;
1371  }
1372 Double_t ICARUSlongHitFitCache::fitlong(Double_t const* x, Double_t const* par)
1373  {
1374  auto const nPeaks = static_cast<std::size_t>(par[0]);
1375  Double_t fitval = 0.0;
1376  for(std::size_t jp = 0; jp < nPeaks; ++jp) {
1377  Double_t const* parj = par + (7 * jp);
1378  int const smax = std::floor(parj[6]);
1379  if (smax == 0) continue;
1380  double const neg_dxj = -(x[0] - parj[3]);
1381  fitval += (smax + parj[7] * (smax*(smax-1)/2))
1382  * (
1383  parj[1]+parj[2]*std::exp(neg_dxj/parj[4])
1384  / (1.0 + std::exp(neg_dxj/parj[5]))
1385  ) / (parj[6]);
1386  } // for
1387  return fitval;
1388  }
1389 
1390  double ICARUSHitFinder::ComputeChiSquare(TF1 func, TH1 *histo) const
1391  {
1392  double chi=0;
1393  int nb=histo->GetNbinsX();
1394  double wb=histo->GetBinWidth(0);
1395 
1396  int jp;
1397  for( jp=1;jp<nb;jp++) {
1398  if(histo->GetBinContent(jp)==0) break;
1399  double xb=jp*wb;
1400  double fv=(func)(xb);
1401  double hv=histo->GetBinContent(jp);
1402  double dv=hv-fv;
1403  double cv=dv/(histo->GetBinError(jp));
1404  chi+=cv*cv;
1405  //std::cout << " chi " << chi << std::endl;
1406 
1407  }
1408 // std::cout << " chi2mio" << chi << std::endl;
1409  //std::cout << " ndf " << ndf << std::endl;
1410  return chi/(jp-5);
1411  }
1412  double ICARUSHitFinder::ComputeNullChiSquare(std::vector<float> holder) const
1413  {
1414  double chi=0;
1415  int nb=33;
1416 
1417  int jp;
1418  for( jp=1;jp<nb;jp++) {
1419 
1420  double hv=holder[jp];
1421  double dv=hv;
1422  double cv=dv/2.4;
1423  chi+=cv*cv;
1424  // std::cout << " chi " << chi << std::endl;
1425 
1426  }
1427  // std::cout << " chi2null" << chi << std::endl;
1428  //std::cout << " ndf " << ndf << std::endl;
1429  return chi/(jp);
1430  }
1431 
1432 
1433  DEFINE_ART_MODULE(ICARUSHitFinder)
1434 
1435 } // end namespace hit
1436 
1437 
1438 #endif //ICARUSHitFinder_H
hit::ICARUSHitFinder::fGeometry
const geo::GeometryCore * fGeometry
Definition: ICARUSHitFinder_module.cc:231
hit::ICARUSHitFinder::fMaxMultiHit
size_t fMaxMultiHit
maximum hits for multi fit
Definition: ICARUSHitFinder_module.cc:191
hit::ICARUSlongHitFitCache::CreateFunction
virtual TF1 * CreateFunction(size_t nFunc) const
Creates and returns the function with specified number of peaks.
Definition: ICARUSHitFinder_module.cc:109
ServiceUtil.h
Utilities related to art service access.
lariov::ChannelStatusProvider::ChannelSet_t
std::set< raw::ChannelID_t > ChannelSet_t
Type of set of channel IDs.
Definition: ChannelStatusProvider.h:56
hit::ICARUSHitFinder::fLongPulseWidthVec
std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
Definition: ICARUSHitFinder_module.cc:194
CryostatGeo.h
Encapsulate the construction of a single cyostat.
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
Hit.h
Declaration of signal hit object.
reco_tool::ICandidateHitFinder::HitCandidate
Definition: ICandidateHitFinder.h:25
p
pdgs p
Definition: selectors.fcl:22
geo::PlaneID::PlaneID_t
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
hit::ICARUSHitFinder::ICARUSPeakFitParams_t
struct ICARUSPeakFitParams{float peakCenter ICARUSPeakFitParams_t
Definition: ICARUSHitFinder_module.cc:139
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
RawDigit.h
Definition of basic raw digits.
hit::ICARUSHitFinder::produce
void produce(art::Event &evt)
Definition: ICARUSHitFinder_module.cc:322
hit::ICARUSHitFinder::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: ICARUSHitFinder_module.cc:249
reconfigure
for pfile in ack l reconfigure(.*) override"` do echo "checking $
Definition: no_reconfigure_override.sh:3
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
hit
process_name hit
Definition: cheaterreco.fcl:51
hit::ICARUSHitFinder::fMinWidth
double fMinWidth
minimum initial width for ICARUS fit
Definition: ICARUSHitFinder_module.cc:220
hit::ICARUShitFitCache::ICARUShitFitCache
ICARUShitFitCache(std::string const &new_name="ICARUShitFitCache")
Constructor (see base class constructor).
Definition: ICARUSHitFinder_module.cc:72
lariov::ChannelStatusProvider::BadChannels
virtual ChannelSet_t BadChannels() const =0
Returns a copy of set of bad channel IDs for the current run.
recob::HitAndAssociationsWriterBase::declare_products
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
hit::ICARUSHitFinder::computeBestLocalMean
void computeBestLocalMean(std::vector< reco_tool::ICandidateHitFinder::HitCandidate > h, std::vector< float > holder, reco_tool::ICandidateHitFinder::MergeHitCandidateVec how, float &localmean)
Definition: ICARUSHitFinder_module.cc:984
hit::ICARUSHitFinder::fMaxWidthMult
double fMaxWidthMult
multiplier for max width for ICARUS fit
Definition: ICARUSHitFinder_module.cc:221
icarus::ns::util::bin
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
h
while getopts h
Definition: getStatistics.sh:13
hit::ICARUSHitFinder::ComputeNullChiSquare
double ComputeNullChiSquare(std::vector< float >) const
Definition: ICARUSHitFinder_module.cc:1412
hit::GausFitCache
A set of TF1 linear sum of base functions (Gaussians)
Definition: GausFitCache.h:46
hit::ICARUSHitFinder::fLongFitCache
ICARUSlongHitFitCache fLongFitCache
Cached functions for long hits.
Definition: ICARUSHitFinder_module.cc:229
reco_tool::ICandidateHitFinder::HitCandidate::stopTick
size_t stopTick
Definition: ICandidateHitFinder.h:27
recob::HitCreator
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:83
reco_tool::ICandidateHitFinder::HitCandidate::hitCenter
float hitCenter
Definition: ICandidateHitFinder.h:32
HitCreator.h
Helper functions to create a hit.
raw::InvalidChannelID
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
recob::HitCollectionCreator
A class handling a collection of hits and its associations.
Definition: HitCreator.h:508
hit::ICARUSHitFinder::fFittingRange
int fFittingRange
semi-width of interval where to fit hit
Definition: ICARUSHitFinder_module.cc:222
hit::ICARUSHitFinder::fHitFinderTool
std::unique_ptr< reco_tool::ICandidateHitFinder > fHitFinderTool
For finding candidate hits.
Definition: ICARUSHitFinder_module.cc:188
raw.h
Collect all the RawData header files together.
hit::ICARUSlongHitFitCache::ICARUSlongHitFitCache
ICARUSlongHitFitCache(std::string const &new_name="ICARUSlongHitFitCache")
Constructor (see base class constructor).
Definition: ICARUSHitFinder_module.cc:99
hit::ICARUSHitFinder::localmeans
std::vector< double > localmeans
Definition: ICARUSHitFinder_module.cc:174
reco_tool::ICandidateHitFinder::HitCandidate::startTick
size_t startTick
Definition: ICandidateHitFinder.h:26
recob::HitCollectionCreator::emplace_back
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
lariov::ChannelStatusProvider
Class providing information about the quality of channels.
Definition: ChannelStatusProvider.h:49
hit::GausFitCache::Get
virtual TF1 * Get(size_t nFunc)
Returns a function sum of nFunc base functions.
Definition: GausFitCache.cxx:43
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
hit::ICARUSHitFinder::ICARUSHitFinder
ICARUSHitFinder(fhicl::ParameterSet const &pset)
Definition: ICARUSHitFinder_module.cc:237
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
hit::ICARUSlongHitFitCache::fitlong
static Double_t fitlong(Double_t const *x, Double_t const *par)
ICARUS hit shape.
Definition: ICARUSHitFinder_module.cc:1372
geo_types.h
Definition of data types for geometry description.
recob::HitCollectionCreator::put_into
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:652
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
hit::ICARUSHitFinder::fFitCache
ICARUShitFitCache fFitCache
Cached functions for multi-peak fits.
Definition: ICARUSHitFinder_module.cc:228
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
hit::ICARUShitFitCache::CreateFunction
virtual TF1 * CreateFunction(size_t nFunc) const
Creates and returns the function with specified number of peaks.
Definition: ICARUSHitFinder_module.cc:82
PlaneGeo.h
Encapsulate the construction of a single detector plane.
hit::ICARUSHitFinder::ICARUSPeakParamsVec
std::vector< ICARUSPeakFitParams_t > ICARUSPeakParamsVec
Definition: ICARUSHitFinder_module.cc:156
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
hit::GausFitCache::FunctionName
virtual std::string FunctionName(size_t nFunc) const
Returns a name for the function with nFunc base functions.
Definition: GausFitCache.cxx:80
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
hit::ICARUSHitFinder::expandHit
void expandHit(reco_tool::ICandidateHitFinder::HitCandidate &h, std::vector< float > holder, std::vector< reco_tool::ICandidateHitFinder::HitCandidate > how)
Definition: ICARUSHitFinder_module.cc:896
ChannelStatusProvider.h
Interface for experiment-specific channel quality info provider.
ICandidateHitFinder.h
This provides an interface for tools which are tasked with finding candidate hits on input waveforms...
GausFitCache.h
Provide caches for TF1 functions to be used with ROOT fitters.
hit::ICARUSHitFinder::ComputeChiSquare
double ComputeChiSquare(TF1 func, TH1 *histo) const
Definition: ICARUSHitFinder_module.cc:1390
hit::ICARUSHitFinder::fLongMaxHitsVec
std::vector< int > fLongMaxHitsVec
maximum Chisquared / NDF allowed for a hit to be saved
Definition: ICARUSHitFinder_module.cc:193
hit::ICARUSlongHitFitCache
Customized function cache for ICARUS long hit shape.
Definition: ICARUSHitFinder_module.cc:95
hit::ICARUShitFitCache
Customized function cache for ICARUS hit shape.
Definition: ICARUSHitFinder_module.cc:68
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
hit::ICARUSHitFinder::findLongPeakParameters
void findLongPeakParameters(const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
Definition: ICARUSHitFinder_module.cc:1230
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
ChannelStatusService.h
Interface for experiment-specific service for channel quality info.
reco_tool::ICandidateHitFinder::MergeHitCandidateVec
std::vector< HitCandidateVec > MergeHitCandidateVec
Definition: ICandidateHitFinder.h:38
raw::Uncompress
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:776
hit::ICARUSHitFinder::fIntegratingRange
int fIntegratingRange
semi-width of interval where to integrate fitting function
Definition: ICARUSHitFinder_module.cc:223
recob::HitCreator::move
recob::Hit && move()
Prepares the constructed hit to be moved away.
Definition: HitCreator.h:343
hit::ICARUSHitFinder::findMultiPeakParameters
void findMultiPeakParameters(const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
Definition: ICARUSHitFinder_module.cc:1079
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
TPCGeo.h
Encapsulate the construction of a single detector plane.
hit::ICARUShitFitCache::fitf
static Double_t fitf(Double_t const *x, Double_t const *par)
ICARUS hit shape.
Definition: ICARUSHitFinder_module.cc:1364
reco_tool::ICandidateHitFinder::HitCandidateVec
std::vector< HitCandidate > HitCandidateVec
Definition: ICandidateHitFinder.h:37
Generated by   doxygen 1.8.5