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