#include <PulseShapeFitOOTPileupCorrection.h>

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

	PulseShapeFitOOTPileupCorrection ()

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

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

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_ = 0

const HcalTimeSlew *	hcalTimeSlewDelay_ = 0

double	tsDelay1GeV_ =0

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

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

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_

bool	pedestalConstraint_

double	pedMean_

double	pedSig_

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)
                                        {
    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 154 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];
 
 }

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

Definition at line 231 of file PulseShapeFitOOTPileupCorrection.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, HcalDetId::depth(), HBHEChannelInfo::fcByPE(), muonCSCDigis_cfi::gain, HBHEChannelInfo::hasTimeInfo(), HBHEChannelInfo::id(), HcalConst::maxSamples, HBHEChannelInfo::nSamples(), psfPtr_, pulseShapeFit(), HBHEChannelInfo::soi(), mathSSE::sqrt(), ts4Chi2_, ts4Max_, ts4Min_, HBHEChannelInfo::tsDFcPerADC(), HBHEChannelInfo::tsGain(), HBHEChannelInfo::tsPedestal(), HBHEChannelInfo::tsPedestalWidth(), HBHEChannelInfo::tsRawCharge(), vts4Chi2_, and vts4Max_.

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 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);
 
     // 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] + 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()) { 
 
     ts4Chi2_=vts4Chi2_[1];
     if(channelData.id().depth()==1) ts4Max_=vts4Max_[1];
     else ts4Max_=vts4Max_[2];
 
   } 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());
   } 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

private

Definition at line 84 of file PulseShapeFitOOTPileupCorrection.cc.

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

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_) {
     if(chargeArr[i]<=1.0) tmpslew[i] = tsDelay1GeV_;
     else tmpslew[i] = hcalTimeSlewDelay_->delay(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 72 of file PulseShapeFitOOTPileupCorrection.cc.

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

Referenced by setPulseShapeTemplate().

                                                                                                                 {
    ++ cntsetPulseShape;
    psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,pedestalConstraint_,timeConstraint_,addPulseJitter_,
                             pulseJitter_,timeMean_,pedMean_,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::setPulseShapeTemplate	(	const HcalPulseShapes::Shape &	ps,
		bool	isHPD,
		unsigned	nSamples,
		const HcalTimeSlew *	hcalTimeSlewDelay
	)

Definition at line 52 of file PulseShapeFitOOTPileupCorrection.cc.

References currentPulseShape_, HcalTimeSlew::delay(), hcalTimeSlewDelay_, isCurrentChannelHPD_, psfPtr_, resetPulseShapeTemplate(), slewFlavor_, timeSigHPD_, timeSigSiPM_, and tsDelay1GeV_.

                                                                                                                                                                  {
   // initialize for every different channel types (HPD vs SiPM)
 
   if (!(&ps == currentPulseShape_ && isHPD == isCurrentChannelHPD_))
     {
       resetPulseShapeTemplate(ps,nSamples);
 
       // redefine the inverttimeSig2
       if(!isHPD) psfPtr_->setinverttimeSig2(1./(timeSigSiPM_*timeSigSiPM_));
       else psfPtr_->setinverttimeSig2(1./(timeSigHPD_*timeSigHPD_));
 
       currentPulseShape_ = &ps;
       isCurrentChannelHPD_ = isHPD;
 
       hcalTimeSlewDelay_ = hcalTimeSlewDelay;
       tsDelay1GeV_= hcalTimeSlewDelay->delay(1.0, slewFlavor_);
 
     }
 }

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 22 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;
 
 }