#include <PulseShapeFitOOTPileupCorrection.h>

Public Member Functions
void	apply (const CaloSamples &cs, const std::vector< int > &capidvec, const HcalCalibrations &calibs, std::vector< double > &correctedOutput) const

	PulseShapeFitOOTPileupCorrection ()

void	resetPulseShapeTemplate (const HcalPulseShapes::Shape &ps)

void	setPulseShapeTemplate (const HcalPulseShapes::Shape &ps)

void	setPUParams (bool iPedestalConstraint, bool iTimeConstraint, bool iAddPulseJitter, bool iUnConstrainedFit, bool iApplyTimeSlew, double iTS4Min, double iTS4Max, double iPulseJitter, double iTimeMean, double iTimeSig, double iPedMean, double iPedSig, double iNoise, double iTMin, double iTMax, double its3Chi2, double its4Chi2, double its345Chi2, double iChargeThreshold, HcalTimeSlew::BiasSetting slewFlavor, int iFitTimes)

	~PulseShapeFitOOTPileupCorrection ()

Private Member Functions
void	fit (int iFit, float &timevalfit, float &chargevalfit, float &pedvalfit, float &chi2, bool &fitStatus, double &iTSMax, const double &iTSTOTen, double *iEnArr, int(&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< double > &fitParsVec) const

Private Attributes
bool	addPulseJitter_

bool	applyTimeSlew_

double	chargeThreshold_

int	cntsetPulseShape

ROOT::Math::Functor *	dpfunctor_

int	fitTimes_

PSFitter::HybridMinimizer *	hybridfitter

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

double	noise_

bool	pedestalConstraint_

double	pedMean_

double	pedSig_

std::auto_ptr < FitterFuncs::PulseShapeFunctor >	psfPtr_

double	pulseJitter_

HcalTimeSlew::BiasSetting	slewFlavor_

ROOT::Math::Functor *	spfunctor_

bool	timeConstraint_

double	timeMean_

double	timeSig_

ROOT::Math::Functor *	tpfunctor_

double	ts345Chi2_

double	ts3Chi2_

double	ts4Chi2_

double	ts4Max_

double	ts4Min_

int	TSMax_

int	TSMin_

bool	unConstrainedFit_

Detailed Description

Definition at line 86 of file PulseShapeFitOOTPileupCorrection.h.

Constructor & Destructor Documentation

PulseShapeFitOOTPileupCorrection::PulseShapeFitOOTPileupCorrection ( )

Definition at line 199 of file PulseShapeFitOOTPileupCorrection.cc.

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

                                                                    : cntsetPulseShape(0), chargeThreshold_(6.),
                                        psfPtr_(nullptr), spfunctor_(nullptr), dpfunctor_(nullptr), tpfunctor_(nullptr),
                                        TSMin_(0), TSMax_(0), ts4Chi2_(0), ts3Chi2_(0), ts345Chi2_(0), pedestalConstraint_(0),
                                        timeConstraint_(0), addPulseJitter_(0), unConstrainedFit_(0), applyTimeSlew_(0),
                                        ts4Min_(0), ts4Max_(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 209 of file PulseShapeFitOOTPileupCorrection.cc.

References dpfunctor_, hybridfitter, spfunctor_, and tpfunctor_.

                                                                     { 
    if(hybridfitter) delete hybridfitter;
    if(spfunctor_)   delete spfunctor_;
    if(dpfunctor_)   delete dpfunctor_;
    if(tpfunctor_)   delete tpfunctor_;
 }

Member Function Documentation

void PulseShapeFitOOTPileupCorrection::apply	(	const CaloSamples &	cs,
		const std::vector< int > &	capidvec,
		const HcalCalibrations &	calibs,
		std::vector< double > &	correctedOutput
	)		const

Definition at line 272 of file PulseShapeFitOOTPileupCorrection.cc.

References RecoTauCleanerPlugins::charge, relval_parameters_module::energy, HcalConst::maxSamples, HcalCalibrations::pedestal(), psfPtr_, pulseShapeFit(), HcalCalibrations::respcorrgain(), CaloSamples::size(), and ts4Min_.

Referenced by heavyIonTools.ConfigureHeavyIons::__call__(), editorTools.UserCodeTool::__call__(), HiCoreTools.RestrictInputToAOD::__call__(), coreTools.RunOnData::__call__(), trackTools.MakeAODTrackCandidates::__call__(), runJetUncertainties.RunJetUncertainties::__call__(), metTools.AddMETCollection::__call__(), heavyIonTools.ProductionDefaults::__call__(), editorTools.ChangeSource::__call__(), HiCoreTools.RemoveMCMatching::__call__(), coreTools.RemoveMCMatching::__call__(), trackTools.MakePATTrackCandidates::__call__(), trigTools.SwitchOnTrigger::__call__(), runMETCorrectionsAndUncertainties.RunMETCorrectionsAndUncertainties::__call__(), heavyIonTools.SelectionDefaults::__call__(), HiCoreTools.RemoveAllPATObjectsBut::__call__(), heavyIonTools.DisbaleMonteCarloDeps::__call__(), HiCoreTools.RemoveSpecificPATObjects::__call__(), trigTools.SwitchOnTriggerStandAlone::__call__(), trackTools.MakeTrackCandidates::__call__(), tauTools.AddTauCollection::__call__(), trigTools.SwitchOnTriggerMatching::__call__(), HiCoreTools.RemoveCleaning::__call__(), HiCoreTools.AddCleaning::__call__(), trigTools.SwitchOnTriggerMatchingStandAlone::__call__(), jetTools.AddJetCollection::__call__(), trigTools.SwitchOnTriggerMatchEmbedding::__call__(), jetTools.SwitchJetCollection::__call__(), jetTools.UpdateJetCollection::__call__(), jetTools.AddJetID::__call__(), jetTools.SetTagInfos::__call__(), HcalSimpleRecAlgoImpl::reco(), and HcalSimpleRecAlgoImpl::recoHBHE().

 {
    psfPtr_->setDefaultcntNANinfit();
 
    const unsigned int cssize = cs.size();
 // initialize arrays to be zero
    double chargeArr[HcalConst::maxSamples]={}, pedArr[HcalConst::maxSamples]={}, gainArr[HcalConst::maxSamples]={};
    double energyArr[HcalConst::maxSamples]={}, pedenArr[HcalConst::maxSamples]={};
    double tsTOT = 0, tstrig = 0; // in fC
    double tsTOTen = 0; // in GeV
    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 = capidvec[ip];
       double charge = cs[ip];
       double ped = calibs.pedestal(capid);
       double gain = calibs.respcorrgain(capid);
 
       double energy = charge*gain;
       double peden = ped*gain;
 
       chargeArr[ip] = charge; pedArr[ip] = ped; gainArr[ip] = gain;
       energyArr[ip] = energy; pedenArr[ip] = peden;
       
       tsTOT += charge - ped;
       tsTOTen += energy - peden;
       if( ip ==4 || ip==5 ){
          tstrig += charge - ped;
       }
    }
    
    std::vector<double> fitParsVec;
    if(tstrig >= ts4Min_&& tsTOTen > 0.) { //Two sigma from 0 
      pulseShapeFit(energyArr, pedenArr, chargeArr, pedArr, gainArr, tsTOTen, fitParsVec);
    }
    else if((tstrig < ts4Min_||tsTOTen < 0.)&&(ts4Min_==0)){
      fitParsVec.clear();
      fitParsVec.push_back(0.);
      fitParsVec.push_back(0.);
      fitParsVec.push_back(0.);
      fitParsVec.push_back(999.);
      fitParsVec.push_back(false);
    }
    correctedOutput.swap(fitParsVec); correctedOutput.push_back(psfPtr_->getcntNANinfit());
 }

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

private

Definition at line 387 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 = 0;
    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];
    if(!(unConstrainedFit_ && iFit == 2)) return;
    //Add the option of the old method 2
    float timeval2fit   = results[2];
    float chargeval2fit = results[3];
    if(std::abs(timevalfit)>std::abs(timeval2fit)) {// if timeval1fit and timeval2fit are differnt, choose the one which is closer to zero
      timevalfit=timeval2fit;
      chargevalfit=chargeval2fit;
    } else if(timevalfit==timeval2fit) { // if the two times are the same, then for charge we just sum the two  
      timevalfit=(timevalfit+timeval2fit)/2;
      chargevalfit=chargevalfit+chargeval2fit;
    } else {
      timevalfit=-999.;
      chargevalfit=-999.;
    }
 }

int PulseShapeFitOOTPileupCorrection::pulseShapeFit	(	const double *	energyArr,
		const double *	pedenArr,
		const double *	chargeArr,
		const double *	pedArr,
		const double *	gainArr,
		const double	tsTOTen,
		std::vector< double > &	fitParsVec
	)		const

private

Definition at line 319 of file PulseShapeFitOOTPileupCorrection.cc.

References funct::abs(), applyTimeSlew_, rpcdqm::BX, chargeThreshold_, beam_dqm_sourceclient-live_cfg::chi2, HcalTimeSlew::delay(), fit(), i, bookConverter::max, HcalConst::maxSamples, noise_, psfPtr_, slewFlavor_, mathSSE::sqrt(), ts4Chi2_, ts4Max_, and unConstrainedFit_.

Referenced by apply().

                                                                                                                                                                                                                                        {
    double tsMAX=0;
    int    nAboveThreshold = 0;
    double tmpx[HcalConst::maxSamples], tmpy[HcalConst::maxSamples], tmperry[HcalConst::maxSamples],tmperry2[HcalConst::maxSamples],tmpslew[HcalConst::maxSamples];
    double tstrig = 0; // in fC
    for(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(std::max(1.0,chargeArr[i]),slewFlavor_); 
       //Add Greg's channel discretization
       double sigmaBin =  psfPtr_->sigma(chargeArr[i]);
       tmperry2[i]=noise_*noise_+ sigmaBin*sigmaBin; //Greg's Granularity
       //Propagate it through
       tmperry2[i]*=(gainArr[i]*gainArr[i]); //Convert from fC to GeV
       tmperry [i]=sqrt(tmperry2[i]);
       //Add the Uncosntrained Double Pulse Switch
       if((chargeArr[i])>chargeThreshold_) nAboveThreshold++;
       if(std::abs(energyArr[i])>tsMAX) tsMAX=std::abs(tmpy[i]);
       if( i ==4 || i ==5 ){
          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;
 
    int BX[3] = {4,5,3};
    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[2] > 3.*tmpy[3]) BX[2] = 2;
 // Only do three-pulse fit when tstrig < ts4Max_, otherwise one-pulse fit is used (above)
    if(chi2 > ts4Chi2_ && !unConstrainedFit_ && tstrig < ts4Max_)   { //fails chi2 cut goes straight to 3 Pulse fit
      fit(3,timevalfit,chargevalfit,pedvalfit,chi2,fitStatus,tsMAX,tsTOTen,tmpy,BX);
      useTriple=true;
    }
    if(unConstrainedFit_ && nAboveThreshold > 5) { //For the old method 2 do double pulse fit if values above a threshold
      fit(2,timevalfit,chargevalfit,pedvalfit,chi2,fitStatus,tsMAX,tsTOTen,tmpy,BX); 
    }
    /*
    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 )

Definition at line 263 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().

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

void PulseShapeFitOOTPileupCorrection::setPulseShapeTemplate ( const HcalPulseShapes::Shape & ps )

Definition at line 252 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().

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

void PulseShapeFitOOTPileupCorrection::setPUParams	(	bool	iPedestalConstraint,
		bool	iTimeConstraint,
		bool	iAddPulseJitter,
		bool	iUnConstrainedFit,
		bool	iApplyTimeSlew,
		double	iTS4Min,
		double	iTS4Max,
		double	iPulseJitter,
		double	iTimeMean,
		double	iTimeSig,
		double	iPedMean,
		double	iPedSig,
		double	iNoise,
		double	iTMin,
		double	iTMax,
		double	its3Chi2,
		double	its4Chi2,
		double	its345Chi2,
		double	iChargeThreshold,
		HcalTimeSlew::BiasSetting	slewFlavor,
		int	iFitTimes
	)

Definition at line 216 of file PulseShapeFitOOTPileupCorrection.cc.

References addPulseJitter_, applyTimeSlew_, chargeThreshold_, fitTimes_, noise_, pedestalConstraint_, pedMean_, pedSig_, pulseJitter_, slewFlavor_, timeConstraint_, timeMean_, timeSig_, ts345Chi2_, ts3Chi2_, ts4Chi2_, ts4Max_, ts4Min_, TSMax_, TSMin_, and unConstrainedFit_.

                                                                                                       { 
 
   TSMin_ = iTMin;
   TSMax_ = iTMax;
   ts3Chi2_   = its3Chi2;
   ts4Chi2_   = its4Chi2;
   ts345Chi2_ = its345Chi2;
   pedestalConstraint_ = iPedestalConstraint;
   timeConstraint_     = iTimeConstraint;
   addPulseJitter_     = iAddPulseJitter;
   unConstrainedFit_   = iUnConstrainedFit;
   applyTimeSlew_      = iApplyTimeSlew;
   ts4Min_             = iTS4Min;
   ts4Max_             = iTS4Max;
   pulseJitter_        = iPulseJitter*iPulseJitter;
   timeMean_           = iTimeMean;
   timeSig_            = iTimeSig;
   pedMean_            = iPedMean;
   pedSig_             = iPedSig;
   noise_              = iNoise;
   slewFlavor_         = slewFlavor;
   chargeThreshold_    = iChargeThreshold;
   fitTimes_           = iFitTimes;
   if(unConstrainedFit_) { //Turn off all Constraints
     //pedestalConstraint_ = false; => Leaving this as tunable
     //timeConstraint_     = false;
     TSMin_ = -100.;
     TSMax_ =   75.;
   }
 }