DeconvGausHFDUNE35t_module.cc
Go to the documentation of this file.
1 #ifndef DECONVGAUSHFDUNE35T_H
2 #define DECONVGAUSHFDUNE35T_H
3 
4 ////////////////////////////////////////////////////////////////////////
5 //
6 // DeconvGausHFDUNE35t class
7 //
8 // jti3@fnal.gov
9 //
10 // 06-21-13 Branched from CalWire_module and added GausHitFinder_module code
11 //
12 // brebel@fnal.gov
13 //
14 // 11-3-09 Pulled all FFT code out and put into Utilitiess/LArFFT
15 //
16 ////////////////////////////////////////////////////////////////////////
17 
18 #include <string>
19 #include <vector>
20 #include <algorithm> // std::accumulate
21 #include <utility> // std::move
22 
29 #include "art_root_io/TFileService.h"
30 #include "art_root_io/TFileDirectory.h"
31 #include "fhiclcpp/ParameterSet.h"
33 #include "cetlib_except/exception.h"
34 
36 #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
40 #include "lardataobj/RawData/raw.h"
45 
46 // ROOT Includes
47 #include "TH1D.h"
48 #include "TDecompSVD.h"
49 #include "TMath.h"
50 #include "TComplex.h"
51 #include "TFile.h"
52 #include "TH2D.h"
53 #include "TF1.h"
54 #include "TTree.h"
55 #include "TCanvas.h"
56 #include "TStyle.h"
57 
58 /* unused function
59 namespace{
60  // Fill histogram from vector (set underflow/overflow bins to zero).
61 
62  void vector_to_hist(const std::vector<double>& v, TH1D* h)
63  {
64  if (!h) throw cet::exception("vector_to_hist") << "no histogram specified\n";
65  int nvec = v.size();
66  int nbins = h->GetNbinsX();
67  int nfill = std::min(nvec, nbins);
68  h->SetBinContent(0, 0.);
69  int zerobin = h->GetXaxis()->FindBin(0.);
70  for(int i=1; i<=nfill; ++i) {
71  if(i >= zerobin)
72  h->SetBinContent(i, v[i - zerobin]);
73  else
74  h->SetBinContent(i, v[i - zerobin + nvec]);
75  }
76  for(int i=nfill+1; i<=nbins+1; ++i)
77  h->SetBinContent(i, 0.);
78  }
79 
80 }
81 */
82 
83 ///creation of calibrated signals on wires
84 namespace deconvgaushf {
85 
87 
88  public:
89 
90  // create calibrated signals on wires. this class runs
91  // an fft to remove the electronics shaping.
92  explicit DeconvGausHFDUNE35t(fhicl::ParameterSet const& pset);
93  virtual ~DeconvGausHFDUNE35t();
94 
95  void produce(art::Event& evt);
96  void beginJob();
97  void endJob();
98  void reconfigure(fhicl::ParameterSet const& p);
99 
100  private:
101  int fDeconvKernSize; //< Length of truncated deconvolution kernel
102  //int fDataSize; ///< size of raw data on one wire // unused
103  int fPostsample; ///< number of postsample bins
104  std::string fDigitModuleLabel; ///< module that made digits
105  ///< constants
106  std::string fSpillName; ///< nominal spill is an empty string
107  ///< it is set by the DigitModuleLabel
108  ///< ex.: "daq:preSpill" for prespill data
110  double fMinSigInd; ///<Induction signal height threshold
111  double fMinSigCol; ///<Collection signal height threshold
112  double fIndWidth; ///<Initial width for induction fit
113  double fColWidth; ///<Initial width for collection fit
114  double fIndMinWidth; ///<Minimum induction hit width
115  double fColMinWidth; ///<Minimum collection hit width
116  int fMaxMultiHit; ///<maximum hits for multi fit
117  int fAreaMethod; ///<Type of area calculation
118  std::vector<double> fAreaNorms; ///<factors for converting area to same units as peak height
119  double fChi2NDF; ///maximum Chisquared / NDF allowed for a hit to be saved
120 
121  //double WireNumber[100000];
122  //double TotalSignal[100000];
123  double StartIime;
125  double EndTime;
126  double EndTimeError;
128  double RMS;
129  double MeanPosition;
130  double MeanPosError;
131  double Amp;
132  double AmpError;
133  double Charge;
134  double ChargeError;
135  double FitGoodnes;
136  int FitNDF;
137 
140 
141 
142  protected:
143 
144  }; // class DeconvGausHFDUNE35t
145 
147 
148  //-------------------------------------------------
150  {
151  fSpillName="";
152  this->reconfigure(pset);
153 
154  // let HitCollectionCreator declare that we are going to produce
155  // hits and associations with wires and raw digits
157  (producesCollector(), fSpillName, false /* doWireAssns */, true /* doRawDigitAssns */);
158 
159  }
160 
161  //-------------------------------------------------
163  {
164  }
165 
166  //////////////////////////////////////////////////////
168  {
169  fDigitModuleLabel = p.get< std::string >("DigitModuleLabel", "daq");
170  fPostsample = p.get< int > ("PostsampleBins");
171  fDeconvKernSize = p.get< int > ("DeconvKernSize");
172 
173  fSpillName="";
174 
175  size_t pos = fDigitModuleLabel.find(":");
176  if( pos!=std::string::npos ) {
177  fSpillName = fDigitModuleLabel.substr( pos+1 );
178  fDigitModuleLabel = fDigitModuleLabel.substr( 0, pos );
179  }
180 
181  // Implementation of optional member function here.
182  //fCalDataModuleLabel = p.get< std::string >("CalDataModuleLabel");
183  fMinSigInd = p.get< double >("MinSigInd");
184  fMinSigCol = p.get< double >("MinSigCol");
185  fIndWidth = p.get< double >("IndWidth");
186  fColWidth = p.get< double >("ColWidth");
187  fIndMinWidth = p.get< double >("IndMinWidth");
188  fColMinWidth = p.get< double >("ColMinWidth");
189  fMaxMultiHit = p.get< int >("MaxMultiHit");
190  fAreaMethod = p.get< int >("AreaMethod");
191  fAreaNorms = p.get< std::vector< double > >("AreaNorms");
192  fChi2NDF = p.get< double >("Chi2NDF");
193 
194 
195  }
196 
197  //-------------------------------------------------
199  {
200  collectioncount = 0;
201  inductioncount = 0;
202 
203  }
204 
205  //////////////////////////////////////////////////////
207  {
208  }
209 
210  //////////////////////////////////////////////////////
212  {
213  // get the geometry
215 
216  // get the FFT service to have access to the FFT size
218  int transformSize = fFFT->FFTSize();
219 
220  // Get signal shaping service.
222 
223  // Make file for induction and collection plane histograms.
225  art::TFileDirectory dirc = tfs->mkdir("DeconvolutionKernels", "DeconvolutionKernels");
226 
227 
228  // don't make a collection of Wires
229  //std::unique_ptr<std::vector<recob::Wire> > wirecol(new std::vector<recob::Wire>);
230 
231 
232  //Gaussian hit finder initializations etc.
233  TH1::AddDirectory(kFALSE);
234 
235 
236  // ---------------------------
237  // --- Seed Fit Functions --- (This function used to "seed" the mean peak position)
238  // ---------------------------
239  TF1 *hit = new TF1("hit","gaus",0,3200); // FIXME use DetectorProperties
240 
241  // ###############################################
242  // ### Making a ptr vector to put on the event ###
243  // ###############################################
244  // this contains the hit collection
245  // and its associations to raw digits (not to wires)
247  evt, fSpillName, false /* doWireAssns */, true /* doRawDigitAssns */
248  );
249 
250 
251  // ##########################################
252  // ### Reading in the Wire List object(s) ###
253  // ##########################################
254  //auto wireVecHandle = evt.getHandle< std::vector<recob::Wire> >(fCalDataModuleLabel);
255  //art::ServiceHandle<geo::Geometry> geom;
256 
257  // #########################################################
258  // ### List of useful variables used throughout the code ###
259  // #########################################################
260  double threshold = 0.; // minimum signal size for id'ing a hit
261  double fitWidth = 0.; //hit fit width initial value
262  double minWidth = 0.; //minimum hit width
263  //channel = 0; // channel number
264  geo::SigType_t sigType; // Signal Type (Collection or Induction)
265  geo::View_t view; // view
266  std::vector<int> startTimes; // stores time of 1st local minimum
267  std::vector<int> maxTimes; // stores time of local maximum
268  std::vector<int> endTimes; // stores time of 2nd local minimum
269  int time = 0; // current time bin
270  int minTimeHolder = 0; // current start time
271 
272  std::string eqn = "gaus(0)"; // string for equation for gaus fit
273  //std::string seed = "gaus(0)"; // string for equation seed gaus fit
274 
275 
276  bool maxFound = false; // Flag for whether a peak > threshold has been found
277  std::stringstream numConv;
278 
279 
280  // Read in the digit List object(s).
283  if (fSpillName.size()>0) digitVecHandle = evt.getHandle< std::vector<raw::RawDigit> >(itag1);
284  else digitVecHandle = evt.getHandle< std::vector<raw::RawDigit> >(fDigitModuleLabel);
285 
286  if (!digitVecHandle->size()) return;
287  mf::LogInfo("DeconvGausHFDUNE35t") << "DeconvGausHFDUNE35t:: digitVecHandle size is " << digitVecHandle->size();
288 
289  // Use the handle to get a particular (0th) element of collection.
290  art::Ptr<raw::RawDigit> digitVec0(digitVecHandle, 0);
291 
292  unsigned int dataSize = digitVec0->Samples(); //size of raw data vectors
293 
294  if (digitVec0->Compression() != raw::kZeroSuppression) {
296  << "CalGausHFDUNE only supports zero-suppressed raw digit input!";
297  } // if
298 
299 
300 
301  raw::ChannelID_t channel = raw::InvalidChannelID; // channel number
302  unsigned int bin(0); // time bin loop variable
303 
304  filter::ChannelFilter *chanFilt = new filter::ChannelFilter();
305 
306  std::vector<float> holder; // holds signal data
307  std::vector<short> rawadc(transformSize); // vector holding uncompressed adc values
308  std::vector<TComplex> freqHolder(transformSize+1); // temporary frequency data
309 
310  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
311 
312  // loop over all wires
313  for(size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter){ // ++ move
314  holder.clear();
315 
316  // get the reference to the current raw::RawDigit
317  art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
318  channel = digitVec->Channel();
319 
320  // skip bad channels
321  if(!chanFilt->BadChannel(channel)) {
322  holder.resize(transformSize);
323 
324  // uncompress the data
325  raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->Compression());
326 
327  // loop over all adc values and subtract the pedestal
328  for(bin = 0; bin < dataSize; ++bin)
329  holder[bin]=(rawadc[bin]-digitVec->GetPedestal());
330 
331  // Do deconvolution.
332  sss->Deconvolute(clockData, channel, holder);
333 
334  //} // end if not a bad channel
335 
336  holder.resize(dataSize,1e-5);
337 
338 
339  //This restores the DC component to signal removed by the deconvolution.
340  if(fPostsample) {
341  double average=0.0;
342  for(bin=0; bin < (unsigned int)fPostsample; ++bin)
343  average+=holder[holder.size()-1-bin]/(double)fPostsample;
344  for(bin = 0; bin < holder.size(); ++bin) holder[bin]-=average;
345  }
346 
347 
348  //wirecol->push_back(recob::Wire(holder,digitVec));
349 
350  //if(wirecol->size() == 0)
351  //mf::LogWarning("DeconvGausHFDUNE35t") << "No wires made for this event.";
352 
353  // Don't save wire collection in event ROOT file
354  //if(fSpillName.size()>0)
355  //evt.put(std::move(wirecol), fSpillName);
356  //else evt.put(std::move(wirecol));
357 
358  //Beginning of Gaussian hit finder loop over signal blocks
359 
360 
361  //##############################
362  //### Looping over the wires ###
363  //##############################
364 
365  //for(size_t wireIter = 0; wireIter < wirecol->size(); wireIter++){
366 
367  //art::Ptr<recob::Wire> wire(wirecol, wireIter);
368 
369  //std::vector<recob::Wire>::iterator wire;
370  //for(wire=wirecol->begin(); wire != wirecol->end(); wire++){
371  // --- Setting Channel Number and Wire Number as well as signal type ---
372  //channel = wire->RawDigit()->Channel();
373  sigType = geom->SignalType(channel);
374  view = geom->View(channel);
375 
376  // -----------------------------------------------------------
377  // -- Clearing variables at the start of looping over wires --
378  // -----------------------------------------------------------
379  startTimes.clear();
380  maxTimes.clear();
381  endTimes.clear();
382  //std::vector<float> signal(wire->Signal());
383  std::vector<float> signal(holder);
384  std::vector<float>::iterator timeIter; // iterator for time bins
385  time = 0;
386  minTimeHolder = 0;
387  maxFound = false;
388  // -- Setting the appropriate signal widths and thresholds --
389  // -- for either the collection or induction plane --
390  if(sigType == geo::kInduction){
391  threshold = fMinSigInd;
392  fitWidth = fIndWidth;
393  minWidth = fIndMinWidth;
394  }//<-- End if Induction Plane
395  else if(sigType == geo::kCollection){
396  threshold = fMinSigCol;
397  fitWidth = fColWidth;
398  minWidth = fColMinWidth;
399  }//<-- End if Collection Plane
400 
401 
402  // ##################################
403  // ### Looping over Signal Vector ###
404  // ##################################
405  for(timeIter = signal.begin();timeIter+2<signal.end();timeIter++){
406  // ##########################################################
407  // ### LOOK FOR A MINIMUM ###
408  // ### Testing if the point timeIter+1 is at a minimum by ###
409  // ### checking if timeIter and timeIter+1 and if it is ###
410  // ### then we add this to the endTimes ###
411  // ##########################################################
412  if(*timeIter > *(timeIter+1) && *(timeIter+1) < *(timeIter+2)){
413  //--- Note: We only keep the a minimum if we've already ---
414  //--- found a point above threshold ---
415  if(maxFound){
416  endTimes.push_back(time+1);
417  maxFound = false;
418  //keep these in case new hit starts right away
419  minTimeHolder = time+2;
420  }
421  else minTimeHolder = time+1;
422 
423  }//<---End Checking if this is a minimum
424 
425 
426  // ########################################################
427  // ### Testing if the point timeIter+1 is a maximum and ###
428  // ### if it and is above threshold then we add it to ###
429  // ### the startTime ###
430  // ########################################################
431  //if not a minimum, test if we are at a local maximum
432  //if so, and the max value is above threshold, add it and proceed.
433  else if(*timeIter < *(timeIter+1) && *(timeIter+1) > *(timeIter+2) && *(timeIter+1) > threshold){
434  maxFound = true;
435  maxTimes.push_back(time+1);
436  startTimes.push_back(minTimeHolder);
437  }
438 
439  time++;
440 
441  }//end loop over signal vec
442 
443  // ###########################################################
444  // ### If there was no minimum found before the end, but a ###
445  // ### maximum was found then add an end point to the end ###
446  // ###########################################################
447  while( maxTimes.size() > endTimes.size() ){
448  endTimes.push_back(signal.size()-1);
449 
450  }//<---End maxTimes.size > endTimes.size
451 
452  // ####################################################
453  // ### If no startTime hit was found skip this wire ###
454  // ####################################################
455  if( startTimes.size() == 0 ){
456  continue;
457  }
458 
459  /////////////////////////////////////////////////////////////////////////////
460  /////////////////////////////////////////////////////////////////////////////
461 
462  // ##########################################################
463  // ### PERFORM THE FITTING, ADDING HITS TO THE HIT VECTOR ###
464  // ##########################################################
465 
466  //All code below does the fitting, adding of hits
467  //to the hit vector and when all wires are complete
468  //saving them
469 
470  double totSig(0); //stores the total hit signal
471  double startT(0); //stores the start time
472  double endT(0); //stores the end time
473  int numHits(0); //number of consecutive hits being fitted
474  int size(0); //size of data vector for fit
475  int hitIndex(0); //index of current hit in sequence
476  double amplitude(0); //fit parameters
477  double minPeakHeight(0); //lowest peak height in multi-hit fit
478 
479 
480  StartIime = 0; // stores the start time of the hit
481  StartTimeError = 0; // stores the error assoc. with the start time of the hit
482  EndTime = 0; // stores the end time of the hit
483  EndTimeError = 0; // stores the error assoc. with the end time of the hit
484  RMS = 0.; // stores the width of the hit
485  MeanPosition = 0; // stores the peak time position of the hit
486  MeanPosError = 0; // stores the error assoc. with thte peak time of the hit
487  Charge = 0; // stores the total charge assoc. with the hit
488  ChargeError = 0; // stores the error on the charge
489  Amp = 0; // stores the amplitude of the hit
490  AmpError = 0; // stores the error assoc. with the amplitude
491  NumOfHits = 0; // stores the multiplicity of the hit
492  FitGoodnes = 0; // stores the Chi2/NDF of the hit
493  FitNDF = -1; // stores the NDF of the hit
494 
495 
496 
497  //stores gaussian paramters first index is the hit number
498  //the second refers to height, position, and width respectively
499  std::vector<double> hitSig;
500 
501  // ########################################
502  // ### Looping over the hitIndex Vector ###
503  // ########################################
504  while( hitIndex < (signed)startTimes.size() ) {
505  // --- Zeroing Start Times and End Times ---
506  startT = 0;
507  endT = 0;
508  // Note: numHits is hardcoded to one here (Need to fix!)
509  numHits=1;
510 
511  minPeakHeight = signal[maxTimes[hitIndex]];
512 
513  // consider adding pulse to group of consecutive hits if:
514  // 1 less than max consecutive hits
515  // 2 we are not at the last point in the signal vector
516  // 3 the height of the dip between the two is greater than threshold/2
517  // 4 and there is no gap between them
518  while(numHits < fMaxMultiHit && numHits+hitIndex < (signed)endTimes.size() &&
519  signal[endTimes[hitIndex+numHits-1]] >threshold/2.0 &&
520  startTimes[hitIndex+numHits] - endTimes[hitIndex+numHits-1] < 2){
521  if(signal[maxTimes[hitIndex+numHits]] < minPeakHeight){
522  minPeakHeight=signal[maxTimes[hitIndex+numHits]];
523 
524  }//<---Only move the mean peakHeight in this hit
525 
526  numHits++;//<---Interate the multihit function
527 
528  }//<---End multihit while loop
529 
530 
531 
532  // ###########################################################
533  // ### Finding the first point from the beginning (startT) ###
534  // ### which is greater then 1/2 the smallest peak ###
535  // ###########################################################
536  startT=startTimes[hitIndex];
537  while(signal[(int)startT] < minPeakHeight/2.0) {startT++;}
538 
539  // #########################################################
540  // ### Finding the first point from the end (endT) which ###
541  // ### is greater then 1/2 the smallest peak ###
542  // #########################################################
543  endT=endTimes[hitIndex+numHits-1];
544  while(signal[(int)endT] <minPeakHeight/2.0) {endT--;}
545 
546 
547  // ###############################################################
548  // ### Putting the current considered hit into a 1-D Histogram ###
549  // ###############################################################
550  //--- Size of hit = endT - startT ---
551  size = (int)(endT-startT);
552 
553  // ###########################################################################
554  // ### Bug Fix: For ADC counts occuring at the end of the ticks range ###
555  // ### the hitfinder incorrectly assigns the size of the hit as a negative ###
556  // ### number...so we fix this to be 0 so that this hit is skipped ###
557  // ###########################################################################
558  if(size < 0){size = 0;}
559 
560 
561  // --- TH1D HitSignal ---
562  TH1D hitSignal("hitSignal","",size,startT,endT);
563 
564  for(int i = (int)startT; i < (int)endT; i++){
565  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
566  // ~~~ Filling the pulse signals into TH1D histograms ~~~
567  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
568  hitSignal.Fill(i,signal[i]);
569 
570  //std::cout<<"i = "<<i<<" , signal[i] = "<<signal[i]<<std::endl;
571  }//<---End for looping over startT and endT
572 
573 
574  // ##################################################################################################
575  // ### Trying to utilize the old approach of builing my TF1 and then implementing this new method ###
576  // ##################################################################################################
577 
578  // ########################################################
579  // ### Building TFormula for seedHit and basic Gaussian ###
580  // ########################################################
581  eqn = "gaus(0)";
582  for(int i = 3; i < numHits*3; i+=3){
583  eqn.append("+gaus(");
584  numConv.str("");
585  numConv << i;
586  eqn.append(numConv.str());
587  eqn.append(")");
588  }
589 
590 
591 
592  // --- TF1 function for GausHit ---
593  TF1 Gaus("Gaus",eqn.c_str(),0,size);
594 
595  // ##############################################################################
596  // ### For multi-pulse hits we loop over all the hits and fit N Gaus function ###
597  // ##############################################################################
598  if(numHits > 1) {
599  // --- Loop over numHits ---
600  for(int i = 0; i < numHits; ++i){
601 
602  // ###################################################################
603  // ### Set the parameters of the Gaus hit based on the seedHit fit ###
604  // ###################################################################
605  //amplitude = 0.5*(threshold+signal[maxTimes[hitIndex+i]]);
606 
607  amplitude = signal[maxTimes[hitIndex+i]];
608  Gaus.SetParameter(3*i,amplitude);
609  Gaus.SetParameter(1+3*i, maxTimes[hitIndex+i]);
610  Gaus.SetParameter(2+3*i, fitWidth);
611  Gaus.SetParLimits(3*i, 0.0, 3.0*amplitude);
612  Gaus.SetParLimits(1+3*i, startT , endT);
613  Gaus.SetParLimits(2+3*i, 0.0, 10.0*fitWidth);
614  }//end loop over hits
615  }//end if numHits > 1
616 
617  // ######################################################################
618  // ### For single pulse hits we perform a fit using a single Gaussian ###
619  // ######################################################################
620  else{
621 
622  Gaus.SetParameters(signal[maxTimes[hitIndex]],maxTimes[hitIndex],fitWidth);
623  Gaus.SetParLimits(0,0.0,1.5*signal[maxTimes[hitIndex]]);
624  Gaus.SetParLimits(1, startT , endT);
625  Gaus.SetParLimits(2,0.0,10.0*fitWidth);
626  }
627 
628  // ####################################################
629  // ### PERFORMING THE TOTAL GAUSSIAN FIT OF THE HIT ###
630  // ####################################################
631  hitSignal.Sumw2();
632 
633  hitSignal.Fit(&Gaus,"QNRW","", startT, endT);
634  //hitSignal.Fit(&Gaus,"QR0LLi","", startT, endT);
635 
636  for(int hitNumber = 0; hitNumber < numHits; ++hitNumber){
637  totSig = 0;
638 
639 
640  // ###########################################################
641  // ### Record this hit if the amplitude is > threshold/2.0 ###
642  // ### and the width is > minWidth ###
643  // ###########################################################
644  if(Gaus.GetParameter(3*hitNumber) > threshold/2.0 &&
645  Gaus.GetParameter(3*hitNumber+2) > minWidth){
646  StartIime = Gaus.GetParameter(3*hitNumber+1) - Gaus.GetParameter(3*hitNumber+2); // position - width
647 
648  EndTime = Gaus.GetParameter(3*hitNumber+1) + Gaus.GetParameter(3*hitNumber+2); // position + width
649 
650  MeanPosition = Gaus.GetParameter(3*hitNumber+1);
651  RMS = Gaus.GetParameter(3*hitNumber+2);
652 
653  //StartTimeError = TMath::Sqrt( (Gaus.GetParError(3*hitNumber+1)*Gaus.GetParError(3*hitNumber+1)) +
654  // (Gaus.GetParError(3*hitNumber+2)*Gaus.GetParError(3*hitNumber+2)));
655 
656  //EndTimeError = TMath::Sqrt( (Gaus.GetParError(3*hitNumber+1)*Gaus.GetParError(3*hitNumber+1)) +
657  // (Gaus.GetParError(3*hitNumber+2)*Gaus.GetParError(3*hitNumber+2)));
658 
659 
660  //MeanPosError = Gaus.GetParError(3*hitNumber+1);
661 
662 
663 
664 
665  hitSig.resize(size);
666  for(int sigPos = 0; sigPos<size; sigPos++){ //<---Loop over the size (endT - startT)
667 
668  hitSig[sigPos] = Gaus.GetParameter(3*hitNumber)*TMath::Gaus(sigPos+startT,Gaus.GetParameter(3*hitNumber+1), Gaus.GetParameter(3*hitNumber+2));
669  totSig+=hitSig[(int)sigPos];
670 
671  }//<---End Signal postion loop
672 
673  // --- Getting the total charge using the area method ---
674  if(fAreaMethod) {
675  totSig = std::sqrt(2*TMath::Pi())*Gaus.GetParameter(3*hitNumber)*Gaus.GetParameter(3*hitNumber+2)/fAreaNorms[(size_t)sigType];
676 
677  }//<---End Area Method
678  Charge = totSig;
679  // ---------------------------------------------------------------------------------
680  // --- chargeErr=TMath::Sqrt(TMath::Pi())*(amplitudeErr*width+widthErr*amplitude); ---
681  // -----------------------------------------------------------------------------------
682  ChargeError = TMath::Sqrt(TMath::Pi())*(Gaus.GetParError(3*hitNumber+0)*Gaus.GetParameter(3*hitNumber+2)+
683  Gaus.GetParError(3*hitNumber+2)*Gaus.GetParameter(3*hitNumber+0)); //estimate from area of Gaussian
684  Amp = Gaus.GetParameter(3*hitNumber);
685  AmpError = Gaus.GetParError(3*hitNumber);
686  NumOfHits = numHits;
687 
688  // #######################################################
689  // ### Using seeded values to get a better estimate of ###
690  // ### the errors associated with the fit ###
691  // #######################################################
692  // -----------------------------------------
693  // --- Determing the goodness of the fit ---
694  // -----------------------------------------
695  hit->SetParameters(Amp,MeanPosition, Gaus.GetParameter(3*hitNumber+2));
696  hit->FixParameter(0,Amp);
697  //hit->SetParLimits(0,Amp/2, Amp*2);
698  //hit->FixParameter(1,MeanPosition);
699  hit->SetParLimits(1,MeanPosition - 3, MeanPosition + 3);
700  //hit->FixParameter(2,Gaus.GetParameter(3*hitNumber+2));//<---Should I be setting the fitWidth initally
701  // #############################
702  // ### Perform Hit on Signal ###
703  // #############################
704 
705 
706  float TempStartTime = MeanPosition - 4;
707  float TempEndTime = MeanPosition + 4;
708 
709  //hitSignal.Fit(hit,"QNRLLi","", TempStartTime, TempEndTime);
710 
711  char outputfilename[100];
712  if(size>0){
713  if(sigType == geo::kInduction && inductioncount < 100){
714  sprintf(outputfilename,"/dune/data/users/jti3/fits/induction_deconvgaus_fit_%d.eps",inductioncount);
715  TCanvas *c1 = new TCanvas("c1","plot");
716  c1->SetFillColor(0);
717  gStyle->SetOptFit(1);
718  hitSignal.Draw();
719  hitSignal.Fit(hit,"QRLLi","", TempStartTime, TempEndTime);
720  //hitSignal.Draw();
721  c1->Print(outputfilename);
722  c1->Close();
723  inductioncount++;
724 
725  }
726  else if(sigType == geo::kCollection && collectioncount < 100){
727  sprintf(outputfilename,"/dune/data/users/jti3/fits/collection_deconvgaus_fit_%d.eps",collectioncount);
728  TCanvas *c1 = new TCanvas("c1","plot");
729  c1->SetFillColor(0);
730  gStyle->SetOptFit(1);
731  hitSignal.Draw();
732  hitSignal.Fit(hit,"QRLLi","", TempStartTime, TempEndTime);
733  //hitSignal.Draw();
734  c1->Print(outputfilename);
735  c1->Close();
736  collectioncount++;
737  }
738 
739  else
740  hitSignal.Fit(hit,"QNRLLi","", TempStartTime, TempEndTime);
741  }
742  else hitSignal.Fit(hit,"QNRLLi","", TempStartTime, TempEndTime);
743 
744 
745  FitGoodnes = hit->GetChisquare() / hit->GetNDF();
746  FitNDF = hit->GetNDF();
747 
748  StartTimeError = TMath::Sqrt( (hit->GetParError(1)*hit->GetParError(1)) +
749  (hit->GetParError(2)*hit->GetParError(2)));
750 
751  EndTimeError = TMath::Sqrt( (hit->GetParError(1)*hit->GetParError(1)) +
752  (hit->GetParError(2)*hit->GetParError(2)));
753 
754 
755  MeanPosError = hit->GetParError(1);
756 
757  // ####################################################################################
758  // ### Skipping hits with a chi2/NDF larger than the value defined in the .fcl file ###
759  // ####################################################################################
760  if(FitGoodnes > fChi2NDF){continue;}
761 
762  // get the WireID for this hit
763  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
764  ///\todo need to have a disambiguation algorithm somewhere in here
765  // for now, just take the first option returned from ChannelToWire
766  geo::WireID wid = wids[0];
767 
768  recob::Hit hit(
769  channel, // channel
770  (raw::TDCtick_t) startT, // start_tick
771  (raw::TDCtick_t) endT, // end_tick
772  MeanPosition, // peak_time
773  MeanPosError, // sigma_peak_time
774  RMS, // rms,
775  Amp, // peak_amplitude,
776  AmpError, // sigma_peak_amplitude
777  std::accumulate // summedADC
778  (signal.begin() + TempStartTime, signal.begin() + TempEndTime, 0.),
779  Charge, // hit_integral
780  ChargeError, // hit_sigma_integral
781  NumOfHits, // multiplicity
782  hitNumber, // local_index
783  FitGoodnes, // goodness_of_fit
784  FitNDF, // dof
785  view, // view
786  sigType, // signal_type
787  wid // wireID
788  );
789  hcol.emplace_back(std::move(hit), digitVec);
790 
791  }//<---End looking at hits over threshold
792 
793  }//<---End Looping over numHits
794 
795 
796  hitIndex+=numHits;
797  }//<---End while hitIndex < startTimes.size()
798 
799 
800  //}//<---End looping over wireIter
801 
802 
803  //} // end loop over nonzero blocks
804 
805 
806  } // end if not a bad channel
807 
808 
809  } // end loop over wires
810 
811  // move the hit collection and the associations into the event
812  hcol.put_into(evt);
813 
814  delete chanFilt;
815 
816 
817  return;
818  }
819 
820 } // end namespace deconvgaushf
821 
822 
823 #endif // DECONVGAUSHFDUNE35T_H
float GetPedestal() const
Definition: RawDigit.h:214
void Deconvolute(detinfo::DetectorClocksData const &clockData, Channel channel, std::vector< T > &func) const
intermediate_table::iterator iterator
const ADCvector_t & ADCs() const
Reference to the compressed ADC count vector.
Definition: RawDigit.h:210
ULong64_t Samples() const
Number of samples in the uncompressed ADC data.
Definition: RawDigit.h:213
double StartIime
maximum Chisquared / NDF allowed for a hit to be saved
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
bool BadChannel(uint32_t channel) const
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
int fMaxMultiHit
maximum hits for multi fit
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:212
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
uint8_t channel
Definition: CRTFragment.hh:201
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
Service to provide microboone-specific signal shaping for simulation (convolution) and reconstruction...
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
art framework interface to geometry description
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
void reconfigure(fhicl::ParameterSet const &p)
Zero Suppression algorithm.
Definition: RawTypes.h:11
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
Helper functions to create a hit.
double fColWidth
Initial width for collection fit.
double fColMinWidth
Minimum collection hit width.
int FFTSize() const
Definition: LArFFT.h:69
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
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:67
Signal from induction planes.
Definition: geo_types.h:145
enum geo::_plane_sigtype SigType_t
def move(depos, offset)
Definition: depos.py:107
T get(std::string const &key) const
Definition: ParameterSet.h:271
double fIndMinWidth
Minimum induction hit width.
int fPostsample
number of postsample bins
p
Definition: test.py:223
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
creation of calibrated signals on wires
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:652
Detector simulation of raw signals on wires.
ProducesCollector & producesCollector() noexcept
double fMinSigInd
Induction signal height threshold.
raw::Compress_t Compression() const
Compression algorithm used to store the ADC counts.
Definition: RawDigit.h:216
double fIndWidth
Initial width for induction fit.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
Declaration of signal hit object.
DeconvGausHFDUNE35t(fhicl::ParameterSet const &pset)
QTextStream & bin(QTextStream &s)
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
TCEvent evt
Definition: DataStructs.cxx:7
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:776
int fAreaMethod
Type of area calculation.
std::vector< double > fAreaNorms
factors for converting area to same units as peak height
double fMinSigCol
Collection signal height threshold.
Signal from collection planes.
Definition: geo_types.h:146