Inheritance diagram for deconvgaushf::DeconvGausHFDUNE35t:

Public Member Functions
	DeconvGausHFDUNE35t (fhicl::ParameterSet const &pset)

virtual	~DeconvGausHFDUNE35t ()

void	produce (art::Event &evt)

void	beginJob ()

void	endJob ()

void	reconfigure (fhicl::ParameterSet const &p)

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
int	fDeconvKernSize

int	fPostsample
	number of postsample bins More...

std::string	fDigitModuleLabel
	constants More...

std::string	fSpillName

std::string	fCalDataModuleLabel

double	fMinSigInd
	Induction signal height threshold. More...

double	fMinSigCol
	Collection signal height threshold. More...

double	fIndWidth
	Initial width for induction fit. More...

double	fColWidth
	Initial width for collection fit. More...

double	fIndMinWidth
	Minimum induction hit width. More...

double	fColMinWidth
	Minimum collection hit width. More...

int	fMaxMultiHit
	maximum hits for multi fit More...

int	fAreaMethod
	Type of area calculation. More...

std::vector< double >	fAreaNorms
	factors for converting area to same units as peak height More...

double	fChi2NDF

double	StartIime
	maximum Chisquared / NDF allowed for a hit to be saved More...

double	StartTimeError

double	EndTime

double	EndTimeError

int	NumOfHits

double	RMS

double	MeanPosition

double	MeanPosError

double	Amp

double	AmpError

double	Charge

double	ChargeError

double	FitGoodnes

int	FitNDF

int	collectioncount

int	inductioncount

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 86 of file DeconvGausHFDUNE35t_module.cc.

Constructor & Destructor Documentation

deconvgaushf::DeconvGausHFDUNE35t::DeconvGausHFDUNE35t ( fhicl::ParameterSet const & pset )

explicit

Definition at line 149 of file DeconvGausHFDUNE35t_module.cc.

                                                                         : EDProducer{pset}
   {
     fSpillName="";
     this->reconfigure(pset);
 
     // let HitCollectionCreator declare that we are going to produce
     // hits and associations with wires and raw digits
     recob::HitCollectionCreator::declare_products
       (producesCollector(), fSpillName, false /* doWireAssns */, true /* doRawDigitAssns */);
     
   }

deconvgaushf::DeconvGausHFDUNE35t::~DeconvGausHFDUNE35t ( )

virtual

Definition at line 162 of file DeconvGausHFDUNE35t_module.cc.

163 {

164 }

Member Function Documentation

void deconvgaushf::DeconvGausHFDUNE35t::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 198 of file DeconvGausHFDUNE35t_module.cc.

   {  
     collectioncount = 0;
     inductioncount = 0;
 
   }

void deconvgaushf::DeconvGausHFDUNE35t::endJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 206 of file DeconvGausHFDUNE35t_module.cc.

207 {

208 }

void deconvgaushf::DeconvGausHFDUNE35t::produce ( art::Event & evt )

virtual

Todo:: need to have a disambiguation algorithm somewhere in here

Implements art::EDProducer.

Definition at line 211 of file DeconvGausHFDUNE35t_module.cc.

   {      
     // get the geometry
     art::ServiceHandle<geo::Geometry> geom;
 
     // get the FFT service to have access to the FFT size
     art::ServiceHandle<util::LArFFT> fFFT;
     int transformSize = fFFT->FFTSize();
 
     // Get signal shaping service.
     art::ServiceHandle<util::SignalShapingServiceDUNE> sss;
 
     // Make file for induction and collection plane histograms.
     art::ServiceHandle<art::TFileService> tfs;
     art::TFileDirectory dirc = tfs->mkdir("DeconvolutionKernels", "DeconvolutionKernels");
     
  
     // don't make a collection of Wires
     //std::unique_ptr<std::vector<recob::Wire> > wirecol(new std::vector<recob::Wire>);
     
 
     //Gaussian hit finder initializations etc.
     TH1::AddDirectory(kFALSE);
     
     
     // ---------------------------
     // --- Seed Fit Functions  --- (This function used to "seed" the mean peak position)
     // ---------------------------
     TF1 *hit    = new TF1("hit","gaus",0,3200); // FIXME use DetectorProperties
     
     // ###############################################
     // ### Making a ptr vector to put on the event ###
     // ###############################################
     // this contains the hit collection
     // and its associations to raw digits (not to wires)
     recob::HitCollectionCreator hcol(
       evt, fSpillName, false /* doWireAssns */, true /* doRawDigitAssns */
       );
     
     
     // ##########################################
     // ### Reading in the Wire List object(s) ###
     // ##########################################
     //auto wireVecHandle = evt.getHandle< std::vector<recob::Wire> >(fCalDataModuleLabel);
     //art::ServiceHandle<geo::Geometry> geom;
     
     // #########################################################
     // ### List of useful variables used throughout the code ###
     // #########################################################  
     double threshold              = 0.;               // minimum signal size for id'ing a hit
     double fitWidth               = 0.;               //hit fit width initial value
     double minWidth               = 0.;               //minimum hit width
     //channel              = 0;                // channel number
     geo::SigType_t sigType;                           // Signal Type (Collection or Induction)
     geo::View_t view;                                // view
     std::vector<int> startTimes;                    // stores time of 1st local minimum
     std::vector<int> maxTimes;                      // stores time of local maximum    
     std::vector<int> endTimes;                      // stores time of 2nd local minimum
     int time             = 0;                         // current time bin
     int minTimeHolder    = 0;                         // current start time
     
     std::string eqn        = "gaus(0)";      // string for equation for gaus fit
     //std::string seed        = "gaus(0)";      // string for equation seed gaus fit
     
     
     bool maxFound        = false;            // Flag for whether a peak > threshold has been found
     std::stringstream numConv;
           
 
     // Read in the digit List object(s). 
     art::InputTag itag1(fDigitModuleLabel, fSpillName);
     art::Handle< std::vector<raw::RawDigit> > digitVecHandle;
     if (fSpillName.size()>0) digitVecHandle = evt.getHandle< std::vector<raw::RawDigit> >(itag1);
     else digitVecHandle = evt.getHandle< std::vector<raw::RawDigit> >(fDigitModuleLabel);
 
     if (!digitVecHandle->size())  return;
     mf::LogInfo("DeconvGausHFDUNE35t") << "DeconvGausHFDUNE35t:: digitVecHandle size is " << digitVecHandle->size();
 
     // Use the handle to get a particular (0th) element of collection.
     art::Ptr<raw::RawDigit> digitVec0(digitVecHandle, 0);
         
     unsigned int dataSize = digitVec0->Samples(); //size of raw data vectors
     
     if (digitVec0->Compression() != raw::kZeroSuppression) {
       throw art::Exception(art::errors::UnimplementedFeature)
         << "CalGausHFDUNE only supports zero-suppressed raw digit input!";
     } // if
 
 
 
     raw::ChannelID_t channel = raw::InvalidChannelID; // channel number
     unsigned int bin(0);     // time bin loop variable
     
     filter::ChannelFilter *chanFilt = new filter::ChannelFilter();  
 
     std::vector<float> holder;                // holds signal data
     std::vector<short> rawadc(transformSize);  // vector holding uncompressed adc values
     std::vector<TComplex> freqHolder(transformSize+1); // temporary frequency data
     
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
 
     // loop over all wires    
     for(size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter){ // ++ move
       holder.clear();
       
       // get the reference to the current raw::RawDigit
       art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
       channel = digitVec->Channel();
 
       // skip bad channels
       if(!chanFilt->BadChannel(channel)) {
         holder.resize(transformSize);
         
         // uncompress the data
         raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->Compression());
         
         // loop over all adc values and subtract the pedestal
         for(bin = 0; bin < dataSize; ++bin) 
           holder[bin]=(rawadc[bin]-digitVec->GetPedestal());
 
         // Do deconvolution.
         sss->Deconvolute(clockData, channel, holder);
 
         //} // end if not a bad channel 
       
       holder.resize(dataSize,1e-5);
 
 
       //This restores the DC component to signal removed by the deconvolution.
       if(fPostsample) {
         double average=0.0;
         for(bin=0; bin < (unsigned int)fPostsample; ++bin) 
           average+=holder[holder.size()-1-bin]/(double)fPostsample;
         for(bin = 0; bin < holder.size(); ++bin) holder[bin]-=average;
       }  
 
       
       //wirecol->push_back(recob::Wire(holder,digitVec));
       
       //if(wirecol->size() == 0)
       //mf::LogWarning("DeconvGausHFDUNE35t") << "No wires made for this event.";
       
       // Don't save wire collection in event ROOT file
       //if(fSpillName.size()>0)
       //evt.put(std::move(wirecol), fSpillName);
       //else evt.put(std::move(wirecol));
       
       //Beginning of Gaussian hit finder loop over signal blocks
       
       
       //##############################
       //### Looping over the wires ###
       //############################## 
       
           //for(size_t wireIter = 0; wireIter < wirecol->size(); wireIter++){
           
           //art::Ptr<recob::Wire> wire(wirecol, wireIter);
           
           //std::vector<recob::Wire>::iterator wire;
           //for(wire=wirecol->begin(); wire != wirecol->end(); wire++){
           // --- Setting Channel Number and Wire Number as well as signal type ---
           //channel = wire->RawDigit()->Channel();
           sigType = geom->SignalType(channel);
           view = geom->View(channel);
           
           // -----------------------------------------------------------
           // -- Clearing variables at the start of looping over wires --
           // -----------------------------------------------------------
           startTimes.clear();
           maxTimes.clear();
           endTimes.clear();
           //std::vector<float> signal(wire->Signal());
           std::vector<float> signal(holder);
           std::vector<float>::iterator timeIter;            // iterator for time bins
           time          = 0;
           minTimeHolder = 0;
           maxFound      = false;
           // -- Setting the appropriate signal widths and thresholds --
           // --    for either the collection or induction plane      --
           if(sigType == geo::kInduction){
             threshold     = fMinSigInd;
             fitWidth      = fIndWidth;
             minWidth      = fIndMinWidth;
           }//<-- End if Induction Plane
           else if(sigType == geo::kCollection){
             threshold = fMinSigCol;
             fitWidth  = fColWidth;
             minWidth  = fColMinWidth;
           }//<-- End if Collection Plane
           
           
           // ##################################
           // ### Looping over Signal Vector ###
           // ##################################
           for(timeIter = signal.begin();timeIter+2<signal.end();timeIter++){
             // ##########################################################
             // ###                LOOK FOR A MINIMUM                  ###
             // ### Testing if the point timeIter+1 is at a minimum by ###
             // ###  checking if timeIter and timeIter+1 and if it is  ###
             // ###         then we add this to the endTimes           ###
             // ##########################################################
             if(*timeIter > *(timeIter+1) && *(timeIter+1) < *(timeIter+2)){
               //--- Note: We only keep the a minimum if we've already ---
               //---          found a point above threshold            ---
               if(maxFound){
                 endTimes.push_back(time+1);
                 maxFound = false;
                 //keep these in case new hit starts right away
                 minTimeHolder = time+2; 
               }
               else minTimeHolder = time+1; 
               
             }//<---End Checking if this is a minimum
             
             
             // ######################################################## 
             // ### Testing if the point timeIter+1 is a maximum and ###
             // ###  if it and is above threshold then we add it to  ###
             // ###                  the startTime                   ###
             // ########################################################
             //if not a minimum, test if we are at a local maximum 
             //if so, and the max value is above threshold, add it and proceed.
             else if(*timeIter < *(timeIter+1) && *(timeIter+1) > *(timeIter+2) && *(timeIter+1) > threshold){ 
               maxFound = true;
               maxTimes.push_back(time+1);
               startTimes.push_back(minTimeHolder);          
             }
             
             time++;
             
           }//end loop over signal vec 
           
           // ###########################################################    
           // ### If there was no minimum found before the end, but a ### 
           // ###  maximum was found then add an end point to the end ###
           // ###########################################################
           while( maxTimes.size() > endTimes.size() ){ 
             endTimes.push_back(signal.size()-1);
             
           }//<---End maxTimes.size > endTimes.size
           
           // ####################################################
           // ### If no startTime hit was found skip this wire ###
           // ####################################################
           if( startTimes.size() == 0 ){
             continue;
           }  
           
           /////////////////////////////////////////////////////////////////////////////
           /////////////////////////////////////////////////////////////////////////////
           
           // ##########################################################
           // ### PERFORM THE FITTING, ADDING HITS TO THE HIT VECTOR ###
           // ##########################################################
           
           //All code below does the fitting, adding of hits
           //to the hit vector and when all wires are complete 
           //saving them 
           
           double totSig(0);                                             //stores the total hit signal
           double startT(0);                                             //stores the start time
           double endT(0);                                               //stores the end time
           int numHits(0);                                               //number of consecutive hits being fitted
           int size(0);                                                  //size of data vector for fit
           int hitIndex(0);                                              //index of current hit in sequence
           double amplitude(0);                                          //fit parameters
           double minPeakHeight(0);                                      //lowest peak height in multi-hit fit
           
           
           StartIime = 0;                                                        // stores the start time of the hit
           StartTimeError = 0;                                           // stores the error assoc. with the start time of the hit
           EndTime = 0;                                                  // stores the end time of the hit
           EndTimeError = 0;                                             // stores the error assoc. with the end time of the hit
           RMS = 0.;                                                     // stores the width of the hit
           MeanPosition = 0;                                             // stores the peak time position of the hit
           MeanPosError = 0;                                             // stores the error assoc. with thte peak time of the hit
           Charge = 0;                                                   // stores the total charge assoc. with the hit
           ChargeError = 0;                                                      // stores the error on the charge
           Amp = 0;                                                      // stores the amplitude of the hit
           AmpError = 0;                                                 // stores the error assoc. with the amplitude
           NumOfHits = 0;                                                // stores the multiplicity of the hit
           FitGoodnes = 0;                                                       // stores the Chi2/NDF of the hit
           FitNDF = -1;                                                          // stores the NDF of the hit
           
           
           
           //stores gaussian paramters first index is the hit number
           //the second refers to height, position, and width respectively
           std::vector<double>  hitSig;  
           
           // ########################################
           // ### Looping over the hitIndex Vector ###
           // ########################################
           while( hitIndex < (signed)startTimes.size() ) {
             // --- Zeroing Start Times and End Times ---
             startT = 0;
             endT   = 0;
             // Note: numHits is hardcoded to one here (Need to fix!)
             numHits=1;
             
             minPeakHeight = signal[maxTimes[hitIndex]];
             
             // consider adding pulse to group of consecutive hits if:
             // 1 less than max consecutive hits
             // 2 we are not at the last point in the signal vector
             // 3 the height of the dip between the two is greater than threshold/2
             // 4 and there is no gap between them
             while(numHits < fMaxMultiHit && numHits+hitIndex < (signed)endTimes.size() && 
                   signal[endTimes[hitIndex+numHits-1]] >threshold/2.0 &&  
                   startTimes[hitIndex+numHits] - endTimes[hitIndex+numHits-1] < 2){
               if(signal[maxTimes[hitIndex+numHits]] < minPeakHeight){ 
                 minPeakHeight=signal[maxTimes[hitIndex+numHits]];
                 
               }//<---Only move the mean peakHeight in this hit
               
               numHits++;//<---Interate the multihit function
               
             }//<---End multihit while loop
             
             
             
             // ###########################################################      
             // ### Finding the first point from the beginning (startT) ###
             // ###     which is greater then 1/2 the smallest peak     ###
             // ###########################################################
             startT=startTimes[hitIndex];
             while(signal[(int)startT] < minPeakHeight/2.0)  {startT++;}
             
             // #########################################################        
             // ### Finding the first point from the end (endT) which ###
             // ###      is greater then 1/2 the smallest peak        ###
             // #########################################################
             endT=endTimes[hitIndex+numHits-1];
             while(signal[(int)endT] <minPeakHeight/2.0) {endT--;}
             
             
             // ###############################################################
             // ### Putting the current considered hit into a 1-D Histogram ###
             // ###############################################################
             //--- Size of hit = endT - startT ---
             size = (int)(endT-startT);
             
             // ###########################################################################
             // ###    Bug Fix: For ADC counts occuring at the end of the ticks range   ###
             // ### the hitfinder incorrectly assigns the size of the hit as a negative ###
             // ###      number...so we fix this to be 0 so that this hit is skipped    ###
             // ###########################################################################
             if(size < 0){size = 0;}
             
             
             // --- TH1D HitSignal ---
             TH1D hitSignal("hitSignal","",size,startT,endT);
             
             for(int i = (int)startT; i < (int)endT; i++){
               // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
               // ~~~ Filling the pulse signals into TH1D histograms ~~~
               // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
               hitSignal.Fill(i,signal[i]);
               
               //std::cout<<"i = "<<i<<" , signal[i] = "<<signal[i]<<std::endl;
             }//<---End for looping over startT and endT
             
             
             // ##################################################################################################
             // ### Trying to utilize the old approach of builing my TF1 and then implementing this new method ###
             // ##################################################################################################
             
             // ########################################################
             // ### Building TFormula for seedHit and basic Gaussian ###
             // ########################################################
             eqn = "gaus(0)";
             for(int i = 3; i < numHits*3; i+=3){
               eqn.append("+gaus(");
               numConv.str("");
               numConv << i;
               eqn.append(numConv.str());
               eqn.append(")");
             }
             
             
             
             // --- TF1 function for GausHit  ---
             TF1 Gaus("Gaus",eqn.c_str(),0,size);
             
             // ##############################################################################
             // ### For multi-pulse hits we loop over all the hits and fit N Gaus function ###
             // ##############################################################################
             if(numHits > 1) {
               // --- Loop over numHits ---
               for(int i = 0; i < numHits; ++i){
                 
                 // ###################################################################
                 // ### Set the parameters of the Gaus hit based on the seedHit fit ###
                 // ###################################################################
                 //amplitude = 0.5*(threshold+signal[maxTimes[hitIndex+i]]);
                 
                 amplitude = signal[maxTimes[hitIndex+i]];
                 Gaus.SetParameter(3*i,amplitude);
                 Gaus.SetParameter(1+3*i, maxTimes[hitIndex+i]);
                 Gaus.SetParameter(2+3*i, fitWidth);
                 Gaus.SetParLimits(3*i, 0.0, 3.0*amplitude);
                 Gaus.SetParLimits(1+3*i, startT , endT);
                 Gaus.SetParLimits(2+3*i, 0.0, 10.0*fitWidth);
               }//end loop over hits
             }//end if numHits > 1
             
             // ######################################################################
             // ### For single pulse hits we perform a fit using a single Gaussian ###
             // ######################################################################
             else{
               
               Gaus.SetParameters(signal[maxTimes[hitIndex]],maxTimes[hitIndex],fitWidth);
               Gaus.SetParLimits(0,0.0,1.5*signal[maxTimes[hitIndex]]);
               Gaus.SetParLimits(1, startT , endT);
               Gaus.SetParLimits(2,0.0,10.0*fitWidth);
             }
             
             // ####################################################
             // ### PERFORMING THE TOTAL GAUSSIAN FIT OF THE HIT ###
             // ####################################################
             hitSignal.Sumw2();
             
             hitSignal.Fit(&Gaus,"QNRW","", startT, endT);
             //hitSignal.Fit(&Gaus,"QR0LLi","", startT, endT);
             
             for(int hitNumber = 0; hitNumber < numHits; ++hitNumber){
               totSig = 0;       
               
               
               // ###########################################################
               // ### Record this hit if the amplitude is > threshold/2.0 ###
               // ###             and the width is > minWidth             ###
               // ###########################################################
               if(Gaus.GetParameter(3*hitNumber) > threshold/2.0 && 
                  Gaus.GetParameter(3*hitNumber+2) > minWidth){
                 StartIime                       = Gaus.GetParameter(3*hitNumber+1) - Gaus.GetParameter(3*hitNumber+2); // position - width
                 
                 EndTime                 = Gaus.GetParameter(3*hitNumber+1) + Gaus.GetParameter(3*hitNumber+2); // position + width
                 
                 MeanPosition            = Gaus.GetParameter(3*hitNumber+1);
                 RMS                     = Gaus.GetParameter(3*hitNumber+2);
                 
                 //StartTimeError                        = TMath::Sqrt( (Gaus.GetParError(3*hitNumber+1)*Gaus.GetParError(3*hitNumber+1)) + 
                 //                                             (Gaus.GetParError(3*hitNumber+2)*Gaus.GetParError(3*hitNumber+2)));
                 
                 //EndTimeError                  = TMath::Sqrt( (Gaus.GetParError(3*hitNumber+1)*Gaus.GetParError(3*hitNumber+1)) + 
                 //                                             (Gaus.GetParError(3*hitNumber+2)*Gaus.GetParError(3*hitNumber+2)));
                 
                 
                 //MeanPosError                  = Gaus.GetParError(3*hitNumber+1);
                 
                 
                 
                 
                 hitSig.resize(size);
                 for(int sigPos = 0; sigPos<size; sigPos++){ //<---Loop over the size (endT - startT)
                   
                   hitSig[sigPos] = Gaus.GetParameter(3*hitNumber)*TMath::Gaus(sigPos+startT,Gaus.GetParameter(3*hitNumber+1), Gaus.GetParameter(3*hitNumber+2));
                   totSig+=hitSig[(int)sigPos];
                   
                 }//<---End Signal postion loop
                 
                 // --- Getting the total charge using the area method ---
                 if(fAreaMethod) {
                   totSig = std::sqrt(2*TMath::Pi())*Gaus.GetParameter(3*hitNumber)*Gaus.GetParameter(3*hitNumber+2)/fAreaNorms[(size_t)sigType];
                   
                 }//<---End Area Method
                 Charge                          = totSig;
                 // ---------------------------------------------------------------------------------
                 // --- chargeErr=TMath::Sqrt(TMath::Pi())*(amplitudeErr*width+widthErr*amplitude); ---
                 // -----------------------------------------------------------------------------------
                 ChargeError                     = TMath::Sqrt(TMath::Pi())*(Gaus.GetParError(3*hitNumber+0)*Gaus.GetParameter(3*hitNumber+2)+
                                                                             Gaus.GetParError(3*hitNumber+2)*Gaus.GetParameter(3*hitNumber+0));   //estimate from area of Gaussian
             Amp                         = Gaus.GetParameter(3*hitNumber);
             AmpError                    = Gaus.GetParError(3*hitNumber);
             NumOfHits                   = numHits;
             
             // #######################################################
             // ### Using seeded values to get a better estimate of ###
             // ###        the errors associated with the fit       ###
             // #######################################################
             // -----------------------------------------
             // --- Determing the goodness of the fit ---
             // -----------------------------------------
             hit->SetParameters(Amp,MeanPosition, Gaus.GetParameter(3*hitNumber+2));
             hit->FixParameter(0,Amp);
             //hit->SetParLimits(0,Amp/2, Amp*2);
             //hit->FixParameter(1,MeanPosition);
             hit->SetParLimits(1,MeanPosition - 3, MeanPosition + 3);
             //hit->FixParameter(2,Gaus.GetParameter(3*hitNumber+2));//<---Should I be setting the fitWidth initally
             // #############################
             // ### Perform Hit on Signal ###
             // #############################
             
             
             float TempStartTime = MeanPosition - 4;
             float TempEndTime   = MeanPosition  + 4;
             
             //hitSignal.Fit(hit,"QNRLLi","", TempStartTime, TempEndTime);
 
             char outputfilename[100];
             if(size>0){
               if(sigType == geo::kInduction && inductioncount < 100){
                 sprintf(outputfilename,"/dune/data/users/jti3/fits/induction_deconvgaus_fit_%d.eps",inductioncount);
                 TCanvas *c1 = new TCanvas("c1","plot");
                 c1->SetFillColor(0);
                 gStyle->SetOptFit(1);
                 hitSignal.Draw();
                 hitSignal.Fit(hit,"QRLLi","", TempStartTime, TempEndTime);
                 //hitSignal.Draw();
                 c1->Print(outputfilename);
                 c1->Close();
                 inductioncount++;
                 
               }
               else if(sigType == geo::kCollection && collectioncount < 100){
                 sprintf(outputfilename,"/dune/data/users/jti3/fits/collection_deconvgaus_fit_%d.eps",collectioncount);
                 TCanvas *c1 = new TCanvas("c1","plot");
                 c1->SetFillColor(0);
                 gStyle->SetOptFit(1);
                 hitSignal.Draw();
                 hitSignal.Fit(hit,"QRLLi","", TempStartTime, TempEndTime);
                 //hitSignal.Draw();
                 c1->Print(outputfilename);
                 c1->Close();
                 collectioncount++;
               }
               
               else
                 hitSignal.Fit(hit,"QNRLLi","", TempStartTime, TempEndTime);
             }
             else hitSignal.Fit(hit,"QNRLLi","", TempStartTime, TempEndTime);
                     
             
             FitGoodnes                  = hit->GetChisquare() / hit->GetNDF();
             FitNDF                      = hit->GetNDF();
             
             StartTimeError                      = TMath::Sqrt( (hit->GetParError(1)*hit->GetParError(1)) + 
                                                                (hit->GetParError(2)*hit->GetParError(2)));
             
             EndTimeError                        = TMath::Sqrt( (hit->GetParError(1)*hit->GetParError(1)) + 
                                                                (hit->GetParError(2)*hit->GetParError(2)));
             
             
             MeanPosError                        = hit->GetParError(1);
             
             // ####################################################################################
             // ### Skipping hits with a chi2/NDF larger than the value defined in the .fcl file ###
             // ####################################################################################
             if(FitGoodnes > fChi2NDF){continue;}
             
             // get the WireID for this hit
             std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
             ///\todo need to have a disambiguation algorithm somewhere in here
             // for now, just take the first option returned from ChannelToWire
             geo::WireID wid = wids[0];
             
             recob::Hit hit(
               channel,                 // channel
               (raw::TDCtick_t) startT, // start_tick
               (raw::TDCtick_t) endT,   // end_tick
               MeanPosition,            // peak_time
               MeanPosError,            // sigma_peak_time
               RMS,                     // rms,
               Amp,                     // peak_amplitude,
               AmpError,                // sigma_peak_amplitude
               std::accumulate          // summedADC
                 (signal.begin() + TempStartTime, signal.begin() + TempEndTime, 0.),
               Charge,                  // hit_integral
               ChargeError,             // hit_sigma_integral
               NumOfHits,               // multiplicity
               hitNumber,               // local_index
               FitGoodnes,              // goodness_of_fit
               FitNDF,                  // dof
               view,                    // view
               sigType,                 // signal_type
               wid                      // wireID
               );
             hcol.emplace_back(std::move(hit), digitVec);
             
               }//<---End looking at hits over threshold
               
             }//<---End Looping over numHits
             
             
             hitIndex+=numHits;
           }//<---End while hitIndex < startTimes.size()
           
                   
           //}//<---End looping over wireIter
                  
           
           //} // end loop over nonzero blocks
         
         
       } // end if not a bad channel 
       
       
     } // end loop over wires
     
     // move the hit collection and the associations into the event
     hcol.put_into(evt);
 
     delete chanFilt;      
     
     
     return;
   }

void deconvgaushf::DeconvGausHFDUNE35t::reconfigure ( fhicl::ParameterSet const & p )

Definition at line 167 of file DeconvGausHFDUNE35t_module.cc.

   {
     fDigitModuleLabel = p.get< std::string >("DigitModuleLabel", "daq");
     fPostsample       = p.get< int >        ("PostsampleBins");
     fDeconvKernSize   = p.get< int >        ("DeconvKernSize");
     
     fSpillName="";
     
     size_t pos = fDigitModuleLabel.find(":");
     if( pos!=std::string::npos ) {
       fSpillName = fDigitModuleLabel.substr( pos+1 );
       fDigitModuleLabel = fDigitModuleLabel.substr( 0, pos );
     }
   
     // Implementation of optional member function here.
     //fCalDataModuleLabel = p.get< std::string  >("CalDataModuleLabel");
     fMinSigInd          = p.get< double       >("MinSigInd");
     fMinSigCol          = p.get< double       >("MinSigCol"); 
     fIndWidth           = p.get< double       >("IndWidth");  
     fColWidth           = p.get< double       >("ColWidth");
     fIndMinWidth        = p.get< double       >("IndMinWidth");
     fColMinWidth        = p.get< double       >("ColMinWidth");                 
     fMaxMultiHit        = p.get< int          >("MaxMultiHit");
     fAreaMethod         = p.get< int          >("AreaMethod");
     fAreaNorms          = p.get< std::vector< double > >("AreaNorms");
     fChi2NDF          = p.get< double       >("Chi2NDF");
   
   
   }

Member Data Documentation

double deconvgaushf::DeconvGausHFDUNE35t::Amp

private

Definition at line 131 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::AmpError

private

Definition at line 132 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::Charge

private

Definition at line 133 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::ChargeError

private

Definition at line 134 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::collectioncount

private

Definition at line 138 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::EndTime

private

Definition at line 125 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::EndTimeError

private

Definition at line 126 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::fAreaMethod

private

Type of area calculation.

Definition at line 117 of file DeconvGausHFDUNE35t_module.cc.

std::vector<double> deconvgaushf::DeconvGausHFDUNE35t::fAreaNorms

private

factors for converting area to same units as peak height

Definition at line 118 of file DeconvGausHFDUNE35t_module.cc.

std::string deconvgaushf::DeconvGausHFDUNE35t::fCalDataModuleLabel

private

Definition at line 109 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fChi2NDF

private

Definition at line 119 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fColMinWidth

private

Minimum collection hit width.

Definition at line 115 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fColWidth

private

Initial width for collection fit.

Definition at line 113 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::fDeconvKernSize

private

Definition at line 101 of file DeconvGausHFDUNE35t_module.cc.

std::string deconvgaushf::DeconvGausHFDUNE35t::fDigitModuleLabel

private

constants

module that made digits

Definition at line 104 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fIndMinWidth

private

Minimum induction hit width.

Definition at line 114 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fIndWidth

private

Initial width for induction fit.

Definition at line 112 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::FitGoodnes

private

Definition at line 135 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::FitNDF

private

Definition at line 136 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::fMaxMultiHit

private

maximum hits for multi fit

Definition at line 116 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fMinSigCol

private

Collection signal height threshold.

Definition at line 111 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::fMinSigInd

private

Induction signal height threshold.

Definition at line 110 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::fPostsample

private

number of postsample bins

Definition at line 103 of file DeconvGausHFDUNE35t_module.cc.

std::string deconvgaushf::DeconvGausHFDUNE35t::fSpillName

private

nominal spill is an empty string it is set by the DigitModuleLabel ex.: "daq:preSpill" for prespill data

Definition at line 106 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::inductioncount

private

Definition at line 139 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::MeanPosError

private

Definition at line 130 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::MeanPosition

private

Definition at line 129 of file DeconvGausHFDUNE35t_module.cc.

int deconvgaushf::DeconvGausHFDUNE35t::NumOfHits

private

Definition at line 127 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::RMS

private

Definition at line 128 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::StartIime

private

maximum Chisquared / NDF allowed for a hit to be saved

Definition at line 123 of file DeconvGausHFDUNE35t_module.cc.

double deconvgaushf::DeconvGausHFDUNE35t::StartTimeError

private

Definition at line 124 of file DeconvGausHFDUNE35t_module.cc.

The documentation for this class was generated from the following file:

dunedataprep/dunedataprep/CalData/DeconvGausHFDUNE35t_module.cc

Public Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation