Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

HBHE: M2 small fixed #21821

Merged
merged 5 commits into from Jan 26, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
19 changes: 18 additions & 1 deletion HLTrigger/Configuration/python/customizeHLTforCMSSW.py
Expand Up @@ -17,6 +17,21 @@
# pset.minGoodStripCharge = cms.PSet(refToPSet_ = cms.string('HLTSiStripClusterChargeCutNone'))
# return process

def customiseFor21821(process):
for producer in producers_by_type(process, "HBHEPhase1Reconstructor"):
producer.algorithm.ts4Max = cms.vdouble(100., 20000., 30000)
del producer.algorithm.pedestalUpperLimit
del producer.algorithm.pedSigmaHPD
del producer.algorithm.pedSigmaSiPM
del producer.algorithm.noiseHPD
del producer.algorithm.noiseSiPM

for producer in producers_by_type(process, "HcalHitReconstructor"):
if hasattr(producer,"puCorrMethod"):
del producer.puCorrMethod

return process

def customiseFor21664_forMahiOn(process):
for producer in producers_by_type(process, "HBHEPhase1Reconstructor"):
producer.algorithm.useMahi = cms.bool(True)
Expand Down Expand Up @@ -66,5 +81,7 @@ def customizeHLTforCMSSW(process, menuType="GRun"):

# add call to action function in proper order: newest last!
# process = customiseFor12718(process)


process = customiseFor21821(process)

return process
3 changes: 2 additions & 1 deletion RecoLocalCalo/HcalRecAlgos/interface/HcalDeterministicFit.h
Expand Up @@ -19,7 +19,8 @@ class HcalDeterministicFit {
~HcalDeterministicFit();

enum FType {shapeLandau, shape205, shape206, shape207};
void init(HcalTimeSlew::ParaSource tsParam, HcalTimeSlew::BiasSetting bias, bool iApplyTimeSlew, PedestalSub pedSubFxn, double respCorr);

void init(HcalTimeSlew::ParaSource tsParam, HcalTimeSlew::BiasSetting bias, bool iApplyTimeSlew, double respCorr);

void phase1Apply(const HBHEChannelInfo& channelData,
float& reconstructedEnergy,
Expand Down
12 changes: 0 additions & 12 deletions RecoLocalCalo/HcalRecAlgos/interface/HcalSimpleRecAlgo.h
Expand Up @@ -18,11 +18,6 @@
#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseContainmentManager.h"
#include "CondFormats/HcalObjects/interface/AbsOOTPileupCorrection.h"

#include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFitOOTPileupCorrection.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalDeterministicFit.h"

#include "RecoLocalCalo/HcalRecAlgos/interface/PedestalSub.h"

#include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"

/** \class HcalSimpleRecAlgo
Expand Down Expand Up @@ -88,13 +83,6 @@ class HcalSimpleRecAlgo {
int puCorrMethod_;

const HcalTimeSlew* hcalTimeSlew_delay_;

std::unique_ptr<PulseShapeFitOOTPileupCorrection> psFitOOTpuCorr_;

std::unique_ptr<PedestalSub> pedSubFxn_;

// S.Brandt Feb19 : Add a pointer to the HLT algo
std::unique_ptr<HcalDeterministicFit> hltOOTpuCorr_;
};

#endif
7 changes: 0 additions & 7 deletions RecoLocalCalo/HcalRecAlgos/interface/PedestalSub.h
Expand Up @@ -11,18 +11,11 @@ class PedestalSub
PedestalSub();
~PedestalSub();

void init(int runCond, float threshold, float quantile);

void calculate(const std::vector<double> & inputCharge, const std::vector<double> & inputPedestal, const std::vector<double> & inputNoise, std::vector<double> & corrCharge, int soi, int nSample) const;

double getCorrection(const std::vector<double> & inputCharge, const std::vector<double> & inputPedestal, const std::vector<double> & inputNoise, int soi, int nSample) const;


private:
float fThreshold;
float fQuantile;
float fCondition;

};

#endif
Expand Up @@ -42,7 +42,6 @@ class PulseShapeFitOOTPileupCorrection
const std::vector<double> & its4Chi2, HcalTimeSlew::BiasSetting slewFlavor, int iFitTimes);

const HcalPulseShapes::Shape* currentPulseShape_=nullptr;
void setChi2Term( bool isHPD );

void setPulseShapeTemplate (const HcalPulseShapes::Shape& ps, bool isHPD, unsigned nSamples);
void resetPulseShapeTemplate(const HcalPulseShapes::Shape& ps, unsigned nSamples);
Expand Down Expand Up @@ -85,11 +84,6 @@ class PulseShapeFitOOTPileupCorrection
double timeSigSiPM_;
double pedMean_;
double pedSig_;
double pedSigHPD_;
double pedSigSiPM_;
double noise_;
double noiseHPD_;
double noiseSiPM_;
HcalTimeSlew::BiasSetting slewFlavor_;

bool isCurrentChannelHPD_;
Expand Down
13 changes: 6 additions & 7 deletions RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h
Expand Up @@ -18,9 +18,9 @@ namespace FitterFuncs{

class PulseShapeFunctor {
public:
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, unsigned int nSamplesToFit);
PulseShapeFunctor(const HcalPulseShapes::Shape& pulse,bool iPedestalConstraint, bool iTimeConstraint,bool iAddPulseJitter,
double iPulseJitter,double iTimeMean,double iPedMean,
unsigned int nSamplesToFit);
~PulseShapeFunctor();

double EvalPulse(const double *pars, unsigned int nPar);
Expand All @@ -35,6 +35,7 @@ namespace FitterFuncs{
void setpsFitslew (double *slew ){ for(int i=0; i<HcalConst::maxSamples; ++i) {psFit_slew [i] = slew [i]; } }
double getSiPMDarkCurrent(double darkCurrent, double fcByPE, double lambda);
void setinvertpedSig2(double x) { invertpedSig2_ = x; }
void setinverttimeSig2(double x) { inverttimeSig2_ = x; }

double singlePulseShapeFunc( const double *x );
double doublePulseShapeFunc( const double *x );
Expand Down Expand Up @@ -63,12 +64,10 @@ namespace FitterFuncs{
double timeMean_;
double timeSig_;
double pedMean_;
double pedSig_;
double noise_;
double timeShift_;

double inverttimeSig_, inverttimeSig2_;
double invertpedSig_, invertpedSig2_;
double inverttimeSig2_;
double invertpedSig2_;
std::array<double,HcalConst::maxSamples> pulse_shape_;
std::array<double,HcalConst::maxSamples> pulse_shape_sum_;

Expand Down
5 changes: 2 additions & 3 deletions RecoLocalCalo/HcalRecAlgos/src/HcalDeterministicFit.cc
Expand Up @@ -17,11 +17,10 @@ HcalDeterministicFit::HcalDeterministicFit() {
HcalDeterministicFit::~HcalDeterministicFit() {
}

void HcalDeterministicFit::init(HcalTimeSlew::ParaSource tsParam, HcalTimeSlew::BiasSetting bias, bool iApplyTimeSlew, PedestalSub pedSubFxn, double respCorr) {
void HcalDeterministicFit::init(HcalTimeSlew::ParaSource tsParam, HcalTimeSlew::BiasSetting bias, bool iApplyTimeSlew, double respCorr) {
fTimeSlew_=tsParam;
fTimeSlewBias_=bias;
applyTimeSlew_=iApplyTimeSlew;
fPedestalSubFxn_=pedSubFxn;
frespCorr_=respCorr;
}

Expand Down Expand Up @@ -130,7 +129,7 @@ void HcalDeterministicFit::phase1Apply(const HBHEChannelInfo& channelData,
}

float ch3,ch4,ch5, i3,n3,nn3, i4,n4,i5,n5;
ch3=0.f,ch4=0.f,ch5=0.f,i3=0.f,n3=0.f,nn3=0.f,i4=0.f,n4=0.f,i5=0.f,n5=0.f;
ch4=0.f,i3=0.f,n3=0.f,nn3=0.f,i4=0.f,n4=0.f,i5=0.f,n5=0.f;

FType fType;
if(channelData.hasTimeInfo() && channelData.recoShape()==205) fType = shape205;
Expand Down
9 changes: 3 additions & 6 deletions RecoLocalCalo/HcalRecAlgos/src/HcalSimpleRecAlgo.cc
Expand Up @@ -24,8 +24,6 @@ HcalSimpleRecAlgo::HcalSimpleRecAlgo(bool correctForTimeslew, bool correctForPul
{
hcalTimeSlew_delay_ = nullptr;
pulseCorr_ = std::make_unique<HcalPulseContainmentManager>(MaximumFractionalError);
pedSubFxn_ = std::make_unique<PedestalSub>();
hltOOTpuCorr_ = std::make_unique<HcalDeterministicFit>();
}


Expand Down Expand Up @@ -278,8 +276,7 @@ namespace HcalSimpleRecAlgoImpl {
const int runnum, const bool useLeak,
const AbsOOTPileupCorrection* pileupCorrection,
const BunchXParameter* bxInfo, const unsigned lenInfo,
const int puCorrMethod, const PulseShapeFitOOTPileupCorrection * psFitOOTpuCorr,
HcalDeterministicFit * hltOOTpuCorr, PedestalSub * hltPedSub, /* whatever don't know what to do with the pointer...*/
const int puCorrMethod,
const HcalTimeSlew* hcalTimeSlew_delay_)
{
double fc_ampl =0, ampl =0, uncorr_ampl =0, m3_ampl =0, maxA = -1.e300;
Expand Down Expand Up @@ -336,7 +333,7 @@ HORecHit HcalSimpleRecAlgo::reconstruct(const HODataFrame& digi, int first, int
HcalTimeSlew::Slow,
runnum_, false, hoPileupCorr_.get(),
bunchCrossingInfo_, lenBunchCrossingInfo_,
puCorrMethod_, psFitOOTpuCorr_.get(),/*hlt*/hltOOTpuCorr_.get(),pedSubFxn_.get(),
puCorrMethod_,
hcalTimeSlew_delay_);
}

Expand All @@ -348,7 +345,7 @@ HcalCalibRecHit HcalSimpleRecAlgo::reconstruct(const HcalCalibDataFrame& digi, i
HcalTimeSlew::Fast,
runnum_, false, nullptr,
bunchCrossingInfo_, lenBunchCrossingInfo_,
puCorrMethod_, psFitOOTpuCorr_.get(),/*hlt*/hltOOTpuCorr_.get(),pedSubFxn_.get(),
puCorrMethod_,
hcalTimeSlew_delay_);
}

Expand Down
14 changes: 8 additions & 6 deletions RecoLocalCalo/HcalRecAlgos/src/HybridMinimizer.cc
Expand Up @@ -190,7 +190,6 @@ bool HybridMinimizer::SetVariable(unsigned int ivar, const std::string & name, d
std::string txtmsg("Wrong index used for the variable " + name);
MN_INFO_MSG2("HybridMinimizer::SetVariable",txtmsg);
MN_INFO_VAL2("HybridMinimizer::SetVariable",minuit2Index);
ivar = minuit2Index;
return false;
}
fState.RemoveLimits(ivar);
Expand Down Expand Up @@ -320,17 +319,19 @@ bool HybridMinimizer::Minimize() {
int strategyLevel = Strategy();
fMinuitFCN->SetErrorDef(ErrorDef() );

/*
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All the commented out parts of code in this HybridMinimizer.cc should be better removed, or alternatively some additional line of comment should be added to explain why they must remain there

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this file there are large chunks of code commented out, as they are apparently not needed: they should be better removed. Nonetheless, I will give my "+1" to this PR just in case it can still be allowed (as "bug fix") in 10_0_X: in that case HybridMinimizer.cc will be adjusted in a forthcoming pull request. If this PR does not enter 10_0_X, or (even better) if it can be adjusted quickly now, I will remove the "+1" waiting for it

if (PrintLevel() >=1) {
// print the real number of maxfcn used (defined in ModularFuncitonMinimizer)
int maxfcn_used = maxfcn;
if (maxfcn_used == 0) {
int nvar = fState.VariableParameters();
maxfcn_used = 200 + 100*nvar + 5*nvar*nvar;
int nvar = fState.VariableParameters();
maxfcn_used = 200 + 100*nvar + 5*nvar*nvar;
}
// std::cout << "HybridMinimizer: Minimize with max-calls " << maxfcn_used
// << " convergence for edm < " << tol << " strategy "
// << strategyLevel << std::endl;
}
*/

// internal minuit messages
MnPrint::SetLevel(PrintLevel() );
Expand Down Expand Up @@ -701,17 +702,18 @@ bool HybridMinimizer::GetMinosError(unsigned int i, double & errLow, double & er
// cut off too small tolerance (they are not needed)
tol = std::max(tol, 0.01);

/*
if (PrintLevel() >=1) {
// print the real number of maxfcn used (defined in MnMinos)
int maxfcn_used = maxfcn;
if (maxfcn_used == 0) {
int nvar = fState.VariableParameters();
maxfcn_used = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
int nvar = fState.VariableParameters();
maxfcn_used = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
}
// std::cout << "HybridMinimizer::GetMinosError for parameter " << i << " " << par_name
// << " using max-calls " << maxfcn_used << ", tolerance " << tol << std::endl;
}

*/

if (runLower) low = minos.Loval(i,maxfcn,tol);
if (runUpper) up = minos.Upval(i,maxfcn,tol);
Expand Down
5 changes: 3 additions & 2 deletions RecoLocalCalo/HcalRecAlgos/src/MahiFit.cc
Expand Up @@ -498,10 +498,11 @@ void MahiFit::resetPulseShapeTemplate(const HcalPulseShapes::Shape& ps) {

// only the pulse shape itself from PulseShapeFunctor is used for Mahi
// the uncertainty terms calculated inside PulseShapeFunctor are used for Method 2 only
psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,false,false,false,false,
1,0,2.5,0,0.00065,1,10));
psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,false,false,false,
1,0,0,10));
pfunctor_ = std::unique_ptr<ROOT::Math::Functor>( new ROOT::Math::Functor(psfPtr_.get(),&FitterFuncs::PulseShapeFunctor::singlePulseShapeFunc, 3) );


}

void MahiFit::nnlsUnconstrainParameter(Index idxp) const {
Expand Down
7 changes: 1 addition & 6 deletions RecoLocalCalo/HcalRecAlgos/src/PedestalSub.cc
Expand Up @@ -5,17 +5,12 @@

using namespace std;

PedestalSub::PedestalSub() : fThreshold(2.7),fQuantile(0.0),fCondition(0){
PedestalSub::PedestalSub() {
}

PedestalSub::~PedestalSub() {
}

void PedestalSub::init(int runCond=0, float threshold=0.0, float quantile=0.0) {
fThreshold=threshold;
fQuantile=quantile;
fCondition=runCond;
}

void PedestalSub::calculate(const std::vector<double> & inputCharge, const std::vector<double> & inputPedestal, const std::vector<double> & inputNoise, std::vector<double> & corrCharge, int soi, int nTS) const {

Expand Down
27 changes: 14 additions & 13 deletions RecoLocalCalo/HcalRecAlgos/src/PulseShapeFitOOTPileupCorrection.cc
Expand Up @@ -9,8 +9,8 @@ PulseShapeFitOOTPileupCorrection::PulseShapeFitOOTPileupCorrection() : cntsetPul
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) {
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} };
}
Expand All @@ -25,13 +25,6 @@ double PulseShapeFitOOTPileupCorrection::getSiPMDarkCurrent(double darkCurrent,
return sqrt(mu/pow(1-lambda,3)) * fcByPE;
}

void PulseShapeFitOOTPileupCorrection::setChi2Term( bool isHPD ) {

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

}


void PulseShapeFitOOTPileupCorrection::setPUParams(bool iPedestalConstraint, bool iTimeConstraint,bool iAddPulseJitter,
bool iApplyTimeSlew,double iTS4Min, const std::vector<double> & iTS4Max,
Expand Down Expand Up @@ -68,17 +61,21 @@ void PulseShapeFitOOTPileupCorrection::setPulseShapeTemplate(const HcalPulseShap

if (!(&ps == currentPulseShape_ && isHPD == isCurrentChannelHPD_))
{
setChi2Term(isHPD);
resetPulseShapeTemplate(ps,nSamples);

// redefine the inverttimeSig2
if(!isHPD) psfPtr_->setinverttimeSig2(1./(timeSigSiPM_*timeSigSiPM_));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Compute 1./(timeSigSiPM_*timeSigSiPM_) and 1./(timeSigHPD_*timeSigHPD_) only once (inside PulseShapeFitOOTPileupCorrection::setPUParams, I would say).

Then, timeSigSiPM_ and timeSigHPD_ do not need to be data member: you could use (e.g.) invertTimeSig2SiPM_ and invertimeSig2HPD_ instead

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@perrotta
the setPulseShapeTemplate is already called once for rechit type (once for HPD and once for siPM)

else psfPtr_->setinverttimeSig2(1./(timeSigHPD_*timeSigHPD_));

currentPulseShape_ = &ps;
isCurrentChannelHPD_ = isHPD;
}
}

void PulseShapeFitOOTPileupCorrection::resetPulseShapeTemplate(const HcalPulseShapes::Shape& ps, unsigned nSamples) {
++ cntsetPulseShape;
psfPtr_.reset(new FitterFuncs::PulseShapeFunctor(ps,pedestalConstraint_,timeConstraint_,addPulseJitter_,applyTimeSlew_,
pulseJitter_,timeMean_,timeSig_,pedMean_,pedSig_,noise_,nSamples));
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) );
Expand Down Expand Up @@ -306,7 +303,11 @@ void PulseShapeFitOOTPileupCorrection::phase1Apply(const HBHEChannelInfo& channe
psfPtr_->setinvertpedSig2(1./(averagePedSig2GeV));

if(channelData.hasTimeInfo()) {
ts4Max_=vts4Max_[1]; ts4Chi2_=vts4Chi2_[1];

ts4Chi2_=vts4Chi2_[1];
if(channelData.id().depth()==1) ts4Max_=vts4Max_[1];
else ts4Max_=vts4Max_[2];

} else {
ts4Max_=vts4Max_[0]; ts4Chi2_=vts4Chi2_[0];
}
Expand Down
16 changes: 5 additions & 11 deletions RecoLocalCalo/HcalRecAlgos/src/PulseShapeFunctor.cc
Expand Up @@ -8,9 +8,9 @@ namespace FitterFuncs{

//Decalare the Pulse object take it in from Hcal and set some options
PulseShapeFunctor::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,unsigned nSamplesToFit) :
bool iPedestalConstraint, bool iTimeConstraint,bool iAddPulseJitter,
double iPulseJitter,double iTimeMean,double iPedMean,
unsigned nSamplesToFit) :
cntNANinfit(0),
acc25nsVec(HcalConst::maxPSshapeBin), diff25nsItvlVec(HcalConst::maxPSshapeBin),
accVarLenIdxZEROVec(HcalConst::nsPerBX), diffVarItvlIdxZEROVec(HcalConst::nsPerBX),
Expand Down Expand Up @@ -48,19 +48,13 @@ namespace FitterFuncs{
timeConstraint_ = iTimeConstraint;
addPulseJitter_ = iAddPulseJitter;
pulseJitter_ = iPulseJitter*iPulseJitter;

// for M2
timeMean_ = iTimeMean;
timeSig_ = iTimeSig;
pedMean_ = iPedMean;
pedSig_ = iPedSig;
noise_ = iNoise;
timeShift_ = 100.;
timeShift_ += 12.5;//we are trying to get BX

inverttimeSig_ = 1./timeSig_;
inverttimeSig2_ = inverttimeSig_*inverttimeSig_;
invertpedSig_ = 1./pedSig_;
invertpedSig2_ = invertpedSig_*invertpedSig_;

nSamplesToFit_ = nSamplesToFit;

}
Expand Down