Skip to content

Commit

Permalink
Merge pull request #6432 from lihux25/forpre3_Hcal25nsConstrained
Browse files Browse the repository at this point in the history 
Hcal 25ns constrained Method 2 for pre3
  • Loading branch information
@cmsbuild
cmsbuild committed Nov 19, 2014
2 parents 04b5d0d + b88853a commit e5b18cb
Show file tree
Hide file tree
Showing 11 changed files with 569 additions and 474 deletions.
9 changes: 6 additions & 3 deletions RecoLocalCalo/HcalRecAlgos/interface/HcalSimpleRecAlgo.h
View file
Expand Up @@ -52,7 +52,7 @@ class HcalSimpleRecAlgo {
void setRecoParams(bool correctForTimeslew, bool correctForPulse, bool setLeakCorrection, int pileupCleaningID, float phaseNS);

// ugly hack related to HB- e-dependent corrections
void setForData(int runnum);
void setForData (int runnum);

// usage of leak correction
void setLeakCorrection();
Expand All @@ -77,8 +77,11 @@ class HcalSimpleRecAlgo {
HcalCalibRecHit reconstruct(const HcalCalibDataFrame& digi, int first, int toadd, const HcalCoder& coder, const HcalCalibrations& calibs) const;

void setpuCorrMethod(int method){ puCorrMethod_ = method; if( puCorrMethod_ ==2 ) psFitOOTpuCorr_ = std::auto_ptr<PulseShapeFitOOTPileupCorrection>(new PulseShapeFitOOTPileupCorrection()); }
void setpufitChrgThr(double chrgThrInput){ if( puCorrMethod_ ==2 ) psFitOOTpuCorr_->setChargeThreshold(chrgThrInput); }

void setpuCorrParams(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, int iFitTimes);

private:
bool correctForTimeslew_;
bool correctForPulse_;
Expand Down
127 changes: 89 additions & 38 deletions RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFitOOTPileupCorrection.h
View file
Expand Up @@ -5,6 +5,7 @@

#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
#include "CalibFormats/HcalObjects/interface/HcalCoder.h"
#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"

Expand All @@ -17,38 +18,64 @@

#include "RecoLocalCalo/HcalRecAlgos/src/HybridMinimizer.h"

namespace HcalConst{

constexpr int maxSamples = 10;
constexpr int maxPSshapeBin = 256;
constexpr int nsPerBX = 25;
constexpr float iniTimeShift = 98.5f;
constexpr double invertnsPerBx = 0.04;

}

namespace FitterFuncs{

class PulseShapeFunctor {
public:
PulseShapeFunctor(const HcalPulseShapes::Shape& pulse);
~PulseShapeFunctor();

double EvalSinglePulse(const std::vector<double>& pars);
double EvalDoublePulse(const std::vector<double>& pars);

void setDefaultcntNANinfit(){ cntNANinfit =0; }
int getcntNANinfit(){ return cntNANinfit; }

void setpsFitx(double *x ){ for(int i=0; i<10; ++i) psFit_x[i] = x[i]; }
void setpsFity(double *y ){ for(int i=0; i<10; ++i) psFit_y[i] = y[i]; }
void setpsFiterry(double *erry ){ for(int i=0; i<10; ++i) psFit_erry[i] = erry[i]; }

double singlePulseShapeFunc( const double *x );
double doublePulseShapeFunc( const double *x );
private:
std::array<float,256> pulse_hist;

int cntNANinfit;

std::vector<float> acc25nsVec, diff25nsItvlVec;
std::vector<float> accVarLenIdxZEROVec, diffVarItvlIdxZEROVec;
std::vector<float> accVarLenIdxMinusOneVec, diffVarItvlIdxMinusOneVec;

std::array<float,10> funcHPDShape(const std::vector<double>& pars);
std::array<float,10> func_DoublePulse_HPDShape(const std::vector<double>& pars);

double psFit_x[10], psFit_y[10], psFit_erry[10];
PulseShapeFunctor(const HcalPulseShapes::Shape& pulse,bool iPedestalConstraint, bool iTimeConstraint,bool iAddPulseJitter,bool iAddTimeSlew,
double iPulseJitter,double iTimeMean,double iTimeSig,double iPedMean,double iPedSig,
double iNoise);
~PulseShapeFunctor();

double EvalPulse(const std::vector<double>& pars);

void setDefaultcntNANinfit(){ cntNANinfit =0; }
int getcntNANinfit(){ return cntNANinfit; }

void setpsFitx(double *x ){ for(int i=0; i<HcalConst::maxSamples; ++i) psFit_x[i] = x[i]; }
void setpsFity(double *y ){ for(int i=0; i<HcalConst::maxSamples; ++i) psFit_y[i] = y[i]; }
void setpsFiterry (double *erry ){ for(int i=0; i<HcalConst::maxSamples; ++i) psFit_erry [i] = erry [i]; }
void setpsFiterry2(double *erry2 ){ for(int i=0; i<HcalConst::maxSamples; ++i) psFit_erry2 [i] = erry2[i]; }
void setpsFitslew (double *slew ){ for(int i=0; i<HcalConst::maxSamples; ++i) {psFit_slew [i] = slew [i]; } }
double sigma(double ifC);
double singlePulseShapeFunc( const double *x );
double doublePulseShapeFunc( const double *x );
double triplePulseShapeFunc( const double *x );

private:
std::array<float,HcalConst::maxPSshapeBin> pulse_hist;

int cntNANinfit;
std::vector<float> acc25nsVec, diff25nsItvlVec;
std::vector<float> accVarLenIdxZEROVec, diffVarItvlIdxZEROVec;
std::vector<float> accVarLenIdxMinusOneVec, diffVarItvlIdxMinusOneVec;
void funcHPDShape(std::array<float,HcalConst::maxSamples> & ntmpbin, const double &pulseTime, const double &pulseHeight,const double &slew);
double psFit_x[HcalConst::maxSamples], psFit_y[HcalConst::maxSamples], psFit_erry[HcalConst::maxSamples], psFit_erry2[HcalConst::maxSamples], psFit_slew[HcalConst::maxSamples];

bool pedestalConstraint_;
bool timeConstraint_;
bool addPulseJitter_;
bool unConstrainedFit_;
double pulseJitter_;
double timeMean_;
double timeSig_;
double pedMean_;
double pedSig_;
double noise_;
double timeShift_;

double inverttimeSig_, inverttimeSig2_;
double invertpedSig_, invertpedSig2_;
};

}
Expand All @@ -60,25 +87,49 @@ class PulseShapeFitOOTPileupCorrection
~PulseShapeFitOOTPileupCorrection();

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

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);

void setPulseShapeTemplate (const HcalPulseShapes::Shape& ps);
void resetPulseShapeTemplate(const HcalPulseShapes::Shape& ps);
void setChargeThreshold(double chargeThrInput){ chargeThreshold_ = chargeThrInput; }

private:

int pulseShapeFit(const double * energyArr, const double * pedenArr, const double *chargeArr, const double *pedArr, const double tsTOTen, std::vector<double> &fitParsVec) 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;
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;

PSFitter::HybridMinimizer * hybridfitter;

int cntsetPulseShape;

std::array<double,10> iniTimesArr;

std::array<double,HcalConst::maxSamples> iniTimesArr;
double chargeThreshold_;
int fitTimes_;

std::auto_ptr<FitterFuncs::PulseShapeFunctor> psfPtr_;

ROOT::Math::Functor *spfunctor_;
ROOT::Math::Functor *dpfunctor_;
ROOT::Math::Functor *tpfunctor_;
int TSMin_;
int TSMax_;
double ts4Chi2_;
double ts3Chi2_;
double ts345Chi2_;
bool pedestalConstraint_;
bool timeConstraint_;
bool addPulseJitter_;
bool unConstrainedFit_;
bool applyTimeSlew_;
double ts4Min_;
double ts4Max_;
double pulseJitter_;
double timeMean_;
double timeSig_;
double pedMean_;
double pedSig_;
double noise_;
HcalTimeSlew::BiasSetting slewFlavor_;
};

#endif // PulseShapeFitOOTPileupCorrection_h
46 changes: 32 additions & 14 deletions RecoLocalCalo/HcalRecAlgos/src/HcalSimpleRecAlgo.cc
View file
Expand Up @@ -12,7 +12,7 @@

constexpr double MaximumFractionalError = 0.002; // 0.2% error allowed from this source
constexpr int HPDShapev3DataNum = 105;
constexpr int HPDShapev3MCNum = 125;
constexpr int HPDShapev3MCNum = 105;

HcalSimpleRecAlgo::HcalSimpleRecAlgo(bool correctForTimeslew, bool correctForPulse, float phaseNS) :
correctForTimeslew_(correctForTimeslew),
Expand All @@ -21,7 +21,7 @@ HcalSimpleRecAlgo::HcalSimpleRecAlgo(bool correctForTimeslew, bool correctForPul
{

pulseCorr_ = std::auto_ptr<HcalPulseContainmentManager>(
new HcalPulseContainmentManager(MaximumFractionalError)
new HcalPulseContainmentManager(MaximumFractionalError)
);
}

Expand Down Expand Up @@ -55,7 +55,21 @@ void HcalSimpleRecAlgo::setRecoParams(bool correctForTimeslew, bool correctForPu
pileupCleaningID_=pileupCleaningID;
}

void HcalSimpleRecAlgo::setForData (int runnum) {
void HcalSimpleRecAlgo::setpuCorrParams(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, int iFitTimes) {
if( iPedestalConstraint ) assert ( iPedSig );
if( iTimeConstraint ) assert( iTimeSig );
psFitOOTpuCorr_->setPUParams(iPedestalConstraint,iTimeConstraint,iAddPulseJitter,iUnConstrainedFit,iApplyTimeSlew,
iTS4Min, iTS4Max, iPulseJitter,iTimeMean,iTimeSig,iPedMean,iPedSig,iNoise,iTMin,iTMax,its3Chi2,its4Chi2,its345Chi2,
iChargeThreshold,HcalTimeSlew::Medium, iFitTimes);
// int shapeNum = HPDShapev3MCNum;
// psFitOOTpuCorr_->setPulseShapeTemplate(theHcalPulseShapes_.getShape(shapeNum));
}

void HcalSimpleRecAlgo::setForData (int runnum) {
runnum_ = runnum;
if( puCorrMethod_ ==2 ){
int shapeNum = HPDShapev3MCNum;
Expand Down Expand Up @@ -303,37 +317,42 @@ namespace HcalSimpleRecAlgoImpl {
const HcalTimeSlew::BiasSetting slewFlavor,
const int runnum, const bool useLeak,
const AbsOOTPileupCorrection* pileupCorrection,
const BunchXParameter* bxInfo, const unsigned lenInfo, const int puCorrMethod, const PulseShapeFitOOTPileupCorrection * psFitOOTpuCorr)
const BunchXParameter* bxInfo, const unsigned lenInfo, const int puCorrMethod, const PulseShapeFitOOTPileupCorrection * psFitOOTpuCorr)// const on end
{
double fc_ampl =0, ampl =0, uncorr_ampl =0, maxA = -1.e300;
int nRead = 0, maxI = -1;
bool leakCorrApplied = false;
float t0 =0, t2 =0;
float time = -9999;

removePileup(digi, coder, calibs, ifirst, n,
pulseCorrect, corr, pileupCorrection,
bxInfo, lenInfo, &maxA, &ampl,
&uncorr_ampl, &fc_ampl, &nRead, &maxI,
&leakCorrApplied, &t0, &t2);
// Disable method 1 inside the removePileup function this way!
// Some code in removePileup does NOT do pileup correction & to make sure maximum share of code
const AbsOOTPileupCorrection * inputAbsOOTpuCorr = ( puCorrMethod == 1 ? pileupCorrection: 0 );

float time = -9999;
removePileup(digi, coder, calibs, ifirst, n,
pulseCorrect, corr, inputAbsOOTpuCorr,
bxInfo, lenInfo, &maxA, &ampl,
&uncorr_ampl, &fc_ampl, &nRead, &maxI,
&leakCorrApplied, &t0, &t2);

if (maxI > 0 && maxI < (nRead - 1))
{
// Handle negative excursions by moving "zero":
float minA=t0;
if (maxA<minA) minA=maxA;
if (t2<minA) minA=t2;
if (minA<0) { maxA-=minA; t0-=minA; t2-=minA; } // positivizes all samples

float wpksamp = (t0 + maxA + t2);
if (wpksamp!=0) wpksamp=(maxA + 2.0*t2) / wpksamp;
time = (maxI - digi.presamples())*25.0 + timeshift_ns_hbheho(wpksamp);

if (slewCorrect) time-=HcalTimeSlew::delay(std::max(1.0,fc_ampl),slewFlavor);

time=time-calibs.timecorr(); // time calibration
}

// Note that uncorr_ampl is always set from outside of method 2!
if( puCorrMethod == 2 ){
std::vector<double> correctedOutput;

Expand All @@ -344,7 +363,6 @@ namespace HcalSimpleRecAlgoImpl {
const int capid = digi[ip].capid();
capidvec.push_back(capid);
}

psFitOOTpuCorr->apply(cs, capidvec, calibs, correctedOutput);
if( correctedOutput.back() == 0 && correctedOutput.size() >1 ){
time = correctedOutput[1]; ampl = correctedOutput[0];
Expand Down

0 comments on commit e5b18cb

Please sign in to comment.