#include <PulseShapeFitOOTPileupCorrection.h>

Public Member Functions
void	phase1Apply (const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2, const HcalTimeSlew *hcalTimeSlew_delay) const

	PulseShapeFitOOTPileupCorrection ()

void	resetPulseShapeTemplate (const HcalPulseShapes::Shape &ps, unsigned nSamples)

void	setChi2Term (bool isHPD)

void	setPulseShapeTemplate (const HcalPulseShapes::Shape &ps, bool isHPD, unsigned nSamples)

void	setPUParams (bool iPedestalConstraint, bool iTimeConstraint, bool iAddPulseJitter, bool iApplyTimeSlew, double iTS4Min, const std::vector< double > &iTS4Max, double iPulseJitter, double iTimeMean, double iTimeSigHPD, double iTimeSigSiPM, double iPedMean, double iTMin, double iTMax, const std::vector< double > &its4Chi2, HcalTimeSlew::BiasSetting slewFlavor, int iFitTimes)

	~PulseShapeFitOOTPileupCorrection ()

Public Attributes
const HcalPulseShapes::Shape *	currentPulseShape_ =nullptr

Private Member Functions
void	fit (int iFit, float &timevalfit, float &chargevalfit, float &pedvalfit, float &chi2, bool &fitStatus, double &iTSMax, const double &iTSTOTen, double *iEnArr, unsigned(&iBX)[3]) const

double	getSiPMDarkCurrent (double darkCurrent, double fcByPE, double lambda) const

int	pulseShapeFit (const double energyArr, const double pedenArr, const double chargeArr, const double pedArr, const double gainArr, const double tsTOTen, std::vector< float > &fitParsVec, const double ADCnoise, unsigned int soi, const HcalTimeSlew *hcalTimeSlew_delay) const

Private Attributes
bool	addPulseJitter_

bool	applyTimeSlew_

double	chargeThreshold_

int	cntsetPulseShape

std::unique_ptr< ROOT::Math::Functor >	dpfunctor_

int	fitTimes_

PSFitter::HybridMinimizer *	hybridfitter

std::array< double, HcalConst::maxSamples >	iniTimesArr

bool	isCurrentChannelHPD_

double	noise_

double	noiseHPD_

double	noiseSiPM_

bool	pedestalConstraint_

double	pedMean_

double	pedSig_

double	pedSigHPD_

double	pedSigSiPM_

std::unique_ptr< FitterFuncs::PulseShapeFunctor >	psfPtr_

double	pulseJitter_

HcalTimeSlew::BiasSetting	slewFlavor_

std::unique_ptr< ROOT::Math::Functor >	spfunctor_

bool	timeConstraint_

double	timeMean_

double	timeSig_

double	timeSigHPD_

double	timeSigSiPM_

std::unique_ptr< ROOT::Math::Functor >	tpfunctor_

double	ts4Chi2_

double	ts4Max_

double	ts4Min_

int	TSMax_

int	TSMin_

bool	unConstrainedFit_

std::vector< double >	vts4Chi2_

std::vector< double >	vts4Max_

Detailed Description

Definition at line 23 of file PulseShapeFitOOTPileupCorrection.h.

Constructor & Destructor Documentation

PulseShapeFitOOTPileupCorrection::PulseShapeFitOOTPileupCorrection ( )

Definition at line 8 of file PulseShapeFitOOTPileupCorrection.cc.

References hybridfitter, iniTimesArr, and PSFitter::HybridMinimizer::kMigrad.

                                                                    : cntsetPulseShape(0),
                                        psfPtr_(nullptr), spfunctor_(nullptr), dpfunctor_(nullptr), tpfunctor_(nullptr),
                                        TSMin_(0), TSMax_(0), vts4Chi2_(0), pedestalConstraint_(false),
                                        timeConstraint_(false), addPulseJitter_(false), applyTimeSlew_(false),
                                        ts4Min_(0), vts4Max_(0), pulseJitter_(0), timeMean_(0), timeSig_(0), pedMean_(0), pedSig_(0),
                                        noise_(0) {
    hybridfitter = new PSFitter::HybridMinimizer(PSFitter::HybridMinimizer::kMigrad);
    iniTimesArr = { {-100,-75,-50,-25,0,25,50,75,100,125} };
 }

PulseShapeFitOOTPileupCorrection::~PulseShapeFitOOTPileupCorrection ( )

Definition at line 18 of file PulseShapeFitOOTPileupCorrection.cc.

References hybridfitter.

                                                                     { 
    if(hybridfitter) delete hybridfitter;
 }

Member Function Documentation

void PulseShapeFitOOTPileupCorrection::fit	(	int	iFit,
		float &	timevalfit,
		float &	chargevalfit,
		float &	pedvalfit,
		float &	chi2,
		bool &	fitStatus,
		double &	iTSMax,
		const double &	iTSTOTen,
		double *	iEnArr,
		unsigned(&)	iBX[3]
	)		const

private

Definition at line 156 of file PulseShapeFitOOTPileupCorrection.cc.

Referenced by trackingPlots.Iteration::modules(), and pulseShapeFit().

                                                                                                                                                                                                                     {
   int n = 3;
   if(iFit == 2) n = 5; //Two   Pulse Fit 
   if(iFit == 3) n = 7; //Three Pulse Fit 
   //Step 1 Single Pulse fit
    float pedMax =  iTSMax;   //=> max timeslice
    float tMin   =  TSMin_;   //Fitting Time Min
    float tMax   =  TSMax_;   //Fitting Time Max
    //Checks to make sure fitting happens
    if(pedMax   < 1.) pedMax = 1.;
    // Set starting values andf step sizes for parameters
    double vstart[n];
    for(int i = 0; i < int((n-1)/2); i++) { 
      vstart[2*i+0] = iniTimesArr[iBX[i]]+timeMean_;
      vstart[2*i+1] = iEnArr[iBX[i]];
    }
    vstart[n-1] = pedMean_;
 
    double step[n];
    for(int i = 0; i < n; i++) step[i] = 0.1;
       
    if(iFit == 1) hybridfitter->SetFunction(*spfunctor_);
    if(iFit == 2) hybridfitter->SetFunction(*dpfunctor_);
    if(iFit == 3) hybridfitter->SetFunction(*tpfunctor_);
    hybridfitter->Clear();
    //Times and amplitudes
    for(int i = 0; i < int((n-1)/2); i++) {
      hybridfitter->SetLimitedVariable(0+i*2, varNames[2*i+0]  , vstart[0+i*2],   step[0+i*2],iniTimesArr[iBX[i]]+tMin, iniTimesArr[ iBX[i] ]+tMax);
      hybridfitter->SetLimitedVariable(1+i*2, varNames[2*i+1]  , vstart[1+i*2],   step[1+i*2],      0, 1.2*iTSTOTEn);
      //Secret Option to fix the time 
      if(timeSig_ < 0) hybridfitter->SetFixedVariable(0+i*2, varNames[2*i+0],vstart[0+i*2]);
    }
    //Pedestal
    if(vstart[n-1] > std::abs(pedMax)) vstart[n-1] = pedMax;
    hybridfitter->SetLimitedVariable(n-1, varNames[n-1], vstart[n-1], step[n-1],-pedMax,pedMax);
    //Secret Option to fix the pedestal
    if(pedSig_ < 0) hybridfitter->SetFixedVariable(n-1,varNames[n-1],vstart[n-1]);
    //a special number to label the initial condition
    chi2=-1;
    //3 fits why?!
    const double *results = nullptr;
    for(int tries=0; tries<=3;++tries){
      if( fitTimes_ != 2 || tries !=1 ){
         hybridfitter->SetMinimizerType(PSFitter::HybridMinimizer::kMigrad);
         fitStatus = hybridfitter->Minimize();
      }
      double chi2valfit = hybridfitter->MinValue();
      const double *newresults = hybridfitter->X();
      if(chi2 == -1 || chi2>chi2valfit+0.01) {
        results=newresults;
        chi2=chi2valfit;
        if( tries == 0 && fitTimes_ == 1 ) break;
        if( tries == 1 && (fitTimes_ == 2 || fitTimes_ ==3 ) ) break;
        if( tries == 2 && fitTimes_ == 4 ) break;
        if( tries == 3 && fitTimes_ == 5 ) break; 
        //Secret option to speed up the fit => perhaps we should drop this
        if(timeSig_ < 0 || pedSig_ < 0) break;
        if(tries==0){
      hybridfitter->SetMinimizerType(PSFitter::HybridMinimizer::kScan);
      fitStatus = hybridfitter->Minimize();
        } else if(tries==1){
      hybridfitter->SetStrategy(1);
        } else if(tries==2){
      hybridfitter->SetStrategy(2);
        }
      } else {
        break;
      }
    }
    assert(results);
 
    timevalfit   = results[0];
    chargevalfit = results[1];
    pedvalfit    = results[n-1];
 
 }

double PulseShapeFitOOTPileupCorrection::getSiPMDarkCurrent	(	double	darkCurrent,
		double	fcByPE,
		double	lambda
	)		const

private

Definition at line 23 of file PulseShapeFitOOTPileupCorrection.cc.

References RPCpg::mu, funct::pow(), and mathSSE::sqrt().

Referenced by phase1Apply().

                                                                                                                   {
   double mu = darkCurrent * 25 / fcByPE;
   return sqrt(mu/pow(1-lambda,3)) * fcByPE;
 }

void PulseShapeFitOOTPileupCorrection::phase1Apply	(	const HBHEChannelInfo &	channelData,
		float &	reconstructedEnergy,
		float &	reconstructedTime,
		bool &	useTriple,
		float &	chi2,
		const HcalTimeSlew *	hcalTimeSlew_delay
	)		const

Definition at line 233 of file PulseShapeFitOOTPileupCorrection.cc.

Referenced by SimpleHBHEPhase1Algo::reconstruct().

 {
 
   psfPtr_->setDefaultcntNANinfit();
 
   const unsigned cssize = channelData.nSamples();
   const unsigned int soi = channelData.soi();
 
   // initialize arrays to be zero
   double chargeArr[HcalConst::maxSamples]={}, pedArr[HcalConst::maxSamples]={}, gainArr[HcalConst::maxSamples]={};
   double energyArr[HcalConst::maxSamples]={}, pedenArr[HcalConst::maxSamples]={};
   double noiseADCArr[HcalConst::maxSamples]={};
   double noiseDCArr[HcalConst::maxSamples]={};
   double noiseArrSq[HcalConst::maxSamples]={};
   double noisePHArr[HcalConst::maxSamples]={};
   double tsTOT = 0, tstrig = 0; // in fC
   double tsTOTen = 0; // in GeV
 
 
   // go over the time slices
   for(unsigned int ip=0; ip<cssize; ++ip){
     if( ip >= (unsigned) HcalConst::maxSamples ) continue; // Too many samples than what we wanna fit (10 is enough...) -> skip them
 
     //      const int capid = channelData.capid(); // not needed
     double charge = channelData.tsRawCharge(ip);
     double ped = channelData.tsPedestal(ip);
     double gain = channelData.tsGain(ip);
 
     double energy = charge*gain;
     double peden = ped*gain;
 
     chargeArr[ip] = charge; pedArr[ip] = ped; gainArr[ip] = gain;
     energyArr[ip] = energy; pedenArr[ip] = peden;
 
     // quantization noise from the ADC (QIE8 or QIE10/11)
     noiseADCArr[ip] = (1./sqrt(12))*channelData.tsDFcPerADC(ip);
 
     // dark current noise relevant for siPM (only if effective pedestal not used)
     noiseDCArr[ip] = 0;
     if(channelData.hasTimeInfo() && !channelData.hasEffectivePedestals() && (charge-ped)>channelData.tsPedestalWidth(ip)) {
       noiseDCArr[ip] = getSiPMDarkCurrent(channelData.darkCurrent(),channelData.fcByPE(),channelData.lambda());
     }
 
     // Photo statistics uncertainties
     //      sigmaFC/FC = 1/sqrt(Ne);
     // Note2. (from kPedro): the output number of photoelectrons after smearing is treated very differently for SiPMs: *each* pe is assigned a different time based on a random generation from the Y11 pulse plus the SimHit time. In HPDs, the overall pulse is shaped all at once using just the SimHit time.
 
     noisePHArr[ip] = 0;
     if((charge-ped)>channelData.tsPedestalWidth(ip)) {
       noisePHArr[ip] = sqrt((charge-ped)*channelData.fcByPE());
     }
 
     // sum all in quadrature
     noiseArrSq[ip]= noiseADCArr[ip]*noiseADCArr[ip] + noiseDCArr[ip]*noiseDCArr[ip] + channelData.tsPedestalWidth(ip)*channelData.tsPedestalWidth(ip) +  noisePHArr[ip]*noisePHArr[ip];
 
     tsTOT += charge - ped;
     tsTOTen += energy - peden;
     if( ip == soi || ip == soi+1 ){
       tstrig += charge - ped;
     }
   }
 
   double averagePedSig2GeV=0.25*(channelData.tsPedestalWidth(0)*channelData.tsPedestalWidth(0)*channelData.tsGain(0)*channelData.tsGain(0) +
                  channelData.tsPedestalWidth(1)*channelData.tsPedestalWidth(1)*channelData.tsGain(1)*channelData.tsGain(1) +
                  channelData.tsPedestalWidth(2)*channelData.tsPedestalWidth(2)*channelData.tsGain(2)*channelData.tsGain(2) +
                  channelData.tsPedestalWidth(3)*channelData.tsPedestalWidth(3)*channelData.tsGain(3)*channelData.tsGain(3));
 
   // redefine the invertpedSig2
   psfPtr_->setinvertpedSig2(1./(averagePedSig2GeV));
 
   if(channelData.hasTimeInfo()) { 
     ts4Max_=vts4Max_[1]; ts4Chi2_=vts4Chi2_[1]; 
   } else {
     ts4Max_=vts4Max_[0]; ts4Chi2_=vts4Chi2_[0];
   }
 
   std::vector<float> fitParsVec;
   if(tstrig >= ts4Min_ && tsTOTen > 0.) { //Two sigma from 0
     pulseShapeFit(energyArr, pedenArr, chargeArr, pedArr, gainArr, tsTOTen, fitParsVec, noiseArrSq, channelData.soi(), hcalTimeSlew_delay);
   } else {
     fitParsVec.clear();
     fitParsVec.push_back(0.); //charge
     fitParsVec.push_back(-9999); // time
     fitParsVec.push_back(0.); // ped
     fitParsVec.push_back(-9999); // chi2
     fitParsVec.push_back(false); // triple
    }
 
 
   reconstructedEnergy = fitParsVec[0];
   reconstructedTime = fitParsVec[1];
   chi2 = fitParsVec[3];
   useTriple = fitParsVec[4];
 
 }

int PulseShapeFitOOTPileupCorrection::pulseShapeFit	(	const double *	energyArr,
		const double *	pedenArr,
		const double *	chargeArr,
		const double *	pedArr,
		const double *	gainArr,
		const double	tsTOTen,
		std::vector< float > &	fitParsVec,
		const double *	ADCnoise,
		unsigned int	soi,
		const HcalTimeSlew *	hcalTimeSlew_delay
	)		const

private

Definition at line 90 of file PulseShapeFitOOTPileupCorrection.cc.

References funct::abs(), applyTimeSlew_, rpcdqm::BX, vertices_cff::chi2, HcalTimeSlew::delay(), fit(), mps_fire::i, createfilelist::int, SiStripPI::max, HcalConst::maxSamples, HcalConst::nsPerBX, psfPtr_, HcalConst::shiftTS, slewFlavor_, mathSSE::sqrt(), ts4Chi2_, and ts4Max_.

Referenced by phase1Apply().

                                                                                                                                                                                                                                                                                                                            {
    double tsMAX=0;
    double tmpx[HcalConst::maxSamples], tmpy[HcalConst::maxSamples], tmperry[HcalConst::maxSamples],tmperry2[HcalConst::maxSamples],tmpslew[HcalConst::maxSamples];
    double tstrig = 0; // in fC
    for(unsigned int i=0;i<HcalConst::maxSamples;++i){
       tmpx[i]=i;
       tmpy[i]=energyArr[i]-pedenArr[i];
       //Add Time Slew !!! does this need to be pedestal subtracted
       tmpslew[i] = 0;
       if(applyTimeSlew_) tmpslew[i] = hcalTimeSlew_delay->delay(std::max(1.0,chargeArr[i]),slewFlavor_); 
       // add the noise components
       tmperry2[i]=noiseArrSq[i];
 
       //Propagate it through
       tmperry2[i]*=(gainArr[i]*gainArr[i]); //Convert from fC to GeV
       tmperry [i]=sqrt(tmperry2[i]);
 
       if(std::abs(energyArr[i])>tsMAX) tsMAX=std::abs(tmpy[i]);
       if( i ==soi || i ==(soi+1) ){
          tstrig += chargeArr[i] - pedArr[i];
       }
    }
    psfPtr_->setpsFitx    (tmpx);
    psfPtr_->setpsFity    (tmpy);
    psfPtr_->setpsFiterry (tmperry);
    psfPtr_->setpsFiterry2(tmperry2);
    psfPtr_->setpsFitslew (tmpslew);
    
    //Fit 1 single pulse
    float timevalfit  = 0;
    float chargevalfit= 0;
    float pedvalfit   = 0;
    float chi2        = 999; //cannot be zero
    bool  fitStatus   = false;
    bool useTriple = false;
 
    unsigned BX[3] = {soi,soi+1,soi-1};
    if(ts4Chi2_ != 0) fit(1,timevalfit,chargevalfit,pedvalfit,chi2,fitStatus,tsMAX,tsTOTen,tmpy,BX);
 // Based on the pulse shape ( 2. likely gives the same performance )
    if(tmpy[soi-2] > 3.*tmpy[soi-1]) BX[2] = soi-2;
 // Only do three-pulse fit when tstrig < ts4Max_, otherwise one-pulse fit is used (above)
    if(chi2 > ts4Chi2_ && tstrig < ts4Max_)   { //fails chi2 cut goes straight to 3 Pulse fit
      fit(3,timevalfit,chargevalfit,pedvalfit,chi2,fitStatus,tsMAX,tsTOTen,tmpy,BX);
      useTriple=true;
    }
 
    timevalfit -= (int(soi)-HcalConst::shiftTS)*HcalConst::nsPerBX;
 
    /*
    if(chi2 > ts345Chi2_)   { //fails do two pulse chi2 for TS5 
      BX[1] = 5;
      fit(3,timevalfit,chargevalfit,pedvalfit,chi2,fitStatus,tsMAX,tsTOTen,BX);
    }
    */
    //Fix back the timeslew
    //if(applyTimeSlew_) timevalfit+=HcalTimeSlew::delay(std::max(1.0,chargeArr[4]),slewFlavor_);
    int outfitStatus = (fitStatus ? 1: 0 );
    fitParsVec.clear();
    fitParsVec.push_back(chargevalfit);
    fitParsVec.push_back(timevalfit);
    fitParsVec.push_back(pedvalfit);
    fitParsVec.push_back(chi2);
    fitParsVec.push_back(useTriple);
    return outfitStatus;
 }

void PulseShapeFitOOTPileupCorrection::resetPulseShapeTemplate	(	const HcalPulseShapes::Shape &	ps,
		unsigned	nSamples
	)

Definition at line 78 of file PulseShapeFitOOTPileupCorrection.cc.

References addPulseJitter_, applyTimeSlew_, cntsetPulseShape, FitterFuncs::PulseShapeFunctor::doublePulseShapeFunc(), dpfunctor_, noise_, pedestalConstraint_, pedMean_, pedSig_, psfPtr_, pulseJitter_, FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc(), spfunctor_, timeConstraint_, timeMean_, timeSig_, tpfunctor_, and FitterFuncs::PulseShapeFunctor::triplePulseShapeFunc().

Referenced by setPulseShapeTemplate().

                                                                                                                 {
    ++ cntsetPulseShape;
    psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,pedestalConstraint_,timeConstraint_,addPulseJitter_,applyTimeSlew_,
                             pulseJitter_,timeMean_,timeSig_,pedMean_,pedSig_,noise_,nSamples));
    spfunctor_    = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3) );
    dpfunctor_    = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::doublePulseShapeFunc, 5) );
    tpfunctor_    = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::triplePulseShapeFunc, 7) );
 
 }

void PulseShapeFitOOTPileupCorrection::setChi2Term ( bool isHPD )

Definition at line 28 of file PulseShapeFitOOTPileupCorrection.cc.

References timeSig_, timeSigHPD_, and timeSigSiPM_.

Referenced by setPulseShapeTemplate().

                                                                {
 
   if(isHPD) timeSig_       = timeSigHPD_;
   else timeSig_            = timeSigSiPM_;
 
 }

void PulseShapeFitOOTPileupCorrection::setPulseShapeTemplate	(	const HcalPulseShapes::Shape &	ps,
		bool	isHPD,
		unsigned	nSamples
	)

Definition at line 66 of file PulseShapeFitOOTPileupCorrection.cc.

References currentPulseShape_, isCurrentChannelHPD_, resetPulseShapeTemplate(), and setChi2Term().

                                                                                                                           {
   // initialize for every different channel types (HPD vs SiPM)
 
   if (!(&ps == currentPulseShape_ && isHPD == isCurrentChannelHPD_))
     {
       setChi2Term(isHPD);
       resetPulseShapeTemplate(ps,nSamples);
       currentPulseShape_ = &ps;
       isCurrentChannelHPD_ = isHPD;
     }
 }

void PulseShapeFitOOTPileupCorrection::setPUParams	(	bool	iPedestalConstraint,
		bool	iTimeConstraint,
		bool	iAddPulseJitter,
		bool	iApplyTimeSlew,
		double	iTS4Min,
		const std::vector< double > &	iTS4Max,
		double	iPulseJitter,
		double	iTimeMean,
		double	iTimeSigHPD,
		double	iTimeSigSiPM,
		double	iPedMean,
		double	iTMin,
		double	iTMax,
		const std::vector< double > &	its4Chi2,
		HcalTimeSlew::BiasSetting	slewFlavor,
		int	iFitTimes
	)

Definition at line 36 of file PulseShapeFitOOTPileupCorrection.cc.

References addPulseJitter_, applyTimeSlew_, fitTimes_, pedestalConstraint_, pedMean_, pulseJitter_, slewFlavor_, timeConstraint_, timeMean_, timeSigHPD_, timeSigSiPM_, ts4Min_, TSMax_, TSMin_, vts4Chi2_, and vts4Max_.

                                                                               {
 
   TSMin_ = iTMin;
   TSMax_ = iTMax;
   //  ts4Chi2_   = its4Chi2;
   vts4Chi2_   = its4Chi2;
   pedestalConstraint_ = iPedestalConstraint;
   timeConstraint_     = iTimeConstraint;
   addPulseJitter_     = iAddPulseJitter;
   applyTimeSlew_      = iApplyTimeSlew;
   ts4Min_             = iTS4Min;
   //  ts4Max_            = iTS4Max;
   vts4Max_            = iTS4Max;
   pulseJitter_        = iPulseJitter*iPulseJitter;
   timeMean_           = iTimeMean;
   //  timeSig_            = iTimeSig;
   timeSigHPD_         = iTimeSigHPD;
   timeSigSiPM_        = iTimeSigSiPM;
   pedMean_            = iPedMean;
   slewFlavor_         = slewFlavor;
   fitTimes_           = iFitTimes;
 
 }