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.

Sign up for GitHub

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

Jump to bottom

Ecal multifit optimized #5321

Merged
merged 8 commits into from Sep 23, 2014
Merged

Ecal multifit optimized #5321

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
155fd7f
first implementation with Eigen
bendavid Sep 7, 2014
ad40c7e
adapt to inclusion of eigen as an external
bendavid Sep 12, 2014
4942e05
optimizations for symmetric matrices and avoiding the use of std set
bendavid Sep 15, 2014
4bd417c
further optimizations and fixes, plus consistent uncertainty calculat…
bendavid Sep 15, 2014
ce58d14
bring setting of OOT energy back in line with latest changes and sile…
bendavid Sep 15, 2014
a095900
make uncalib rechit producer multithreaded
bendavid Sep 15, 2014
3ea6128
remove no longer needed code
bendavid Sep 16, 2014
a57e67a
cleanup debug printouts
bendavid Sep 16, 2014
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
2 changes: 2 additions & 0 deletions RecoLocalCalo/EcalRecAlgos/BuildFile.xml
View file
Open in desktop
Expand Up @@ -12,6 +12,8 @@
<use name="root"/>
<use name="rootminuit"/>

<use name="eigen"/>

<export>
<lib name="1"/>
</export>
2 changes: 1 addition & 1 deletion RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitMultiFitAlgo.h
View file
Open in desktop
Expand Up @@ -25,7 +25,7 @@ class EcalUncalibRecHitMultiFitAlgo

EcalUncalibRecHitMultiFitAlgo();
~EcalUncalibRecHitMultiFitAlgo() { };
EcalUncalibratedRecHit makeRecHit(const EcalDataFrame& dataFrame, const EcalPedestals::Item * aped, const EcalMGPAGainRatio * aGain, const TMatrixDSym &noisecor, const TVectorD &fullpulse, const TMatrixDSym &fullpulsecov, const std::set<int> &activeBX);
EcalUncalibratedRecHit makeRecHit(const EcalDataFrame& dataFrame, const EcalPedestals::Item * aped, const EcalMGPAGainRatio * aGain, const SampleMatrix &noisecor, const FullSampleVector &fullpulse, const FullSampleMatrix &fullpulsecov, const BXVector &activeBX);
void disableErrorCalculation() { _computeErrors = false; }

private:
Expand Down
7 changes: 5 additions & 2 deletions RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimeWeightsAlgo.h
View file
Open in desktop
Expand Up @@ -17,6 +17,9 @@
#include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
#include "CondFormats/EcalObjects/interface/EcalWeightSet.h"

#include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"


#include "TVectorD.h"
#include <vector>

Expand All @@ -28,7 +31,7 @@ template<class C> class EcalUncalibRecHitTimeWeightsAlgo
virtual ~EcalUncalibRecHitTimeWeightsAlgo<C>() { };

/// Compute time
double time(const C& dataFrame, const std::vector<double> &amplitudes, const EcalPedestals::Item * aped, const EcalMGPAGainRatio * aGain, const TVectorD &fullpulse, const EcalWeightSet::EcalWeightMatrix** weights) {
double time(const C& dataFrame, const std::vector<double> &amplitudes, const EcalPedestals::Item * aped, const EcalMGPAGainRatio * aGain, const FullSampleVector &fullpulse, const EcalWeightSet::EcalWeightMatrix** weights) {

const unsigned int nsample = EcalDataFrame::MAXSAMPLES;

Expand Down Expand Up @@ -84,7 +87,7 @@ template<class C> class EcalUncalibRecHitTimeWeightsAlgo
int ipulse = std::distance(amplitudes.begin(),amplit);
int bx = ipulse - 5;
int firstsamplet = std::max(0,bx + 3);
int offset = -3-bx;
int offset = 7-3-bx;
Copy link
Contributor

Choose a reason for hiding this comment

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

This offset appears in exactly the same form three times. Can you make it a simple function of bx or at least define "7-3" once somewhere (could you elaborate for me where is the transition -3 -> 7-3 coming from? having to do with time slots I guess but still not clear why the result).


TVectorD pulse;
pulse.ResizeTo(nsample);
Expand Down
11 changes: 11 additions & 0 deletions RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h
View file
Open in desktop
@@ -0,0 +1,11 @@
#include <Eigen/Dense>

typedef Eigen::Matrix<double,10,1> SampleVector;
typedef Eigen::Matrix<double,19,1> FullSampleVector;
typedef Eigen::Matrix<double,Eigen::Dynamic,1,0,10,1> PulseVector;
typedef Eigen::Matrix<char,Eigen::Dynamic,1,0,10,1> BXVector;
typedef Eigen::Matrix<double,10,10> SampleMatrix;
typedef Eigen::Matrix<double,19,19> FullSampleMatrix;
typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,0,10,10> PulseMatrix;
typedef Eigen::Matrix<double,10,Eigen::Dynamic,0,10,10> SamplePulseMatrix;
typedef Eigen::LLT<SampleMatrix> SampleDecompLLT;
68 changes: 26 additions & 42 deletions RecoLocalCalo/EcalRecAlgos/interface/PulseChiSqSNNLS.h
View file
Open in desktop
@@ -1,69 +1,53 @@
#ifndef PulseChiSqSNNLS_h
#define PulseChiSqSNNLS_h

#include "TMatrixD.h"
#include "TVectorD.h"
#include "TMatrixDSym.h"
#include "Math/Minimizer.h"
#include "Math/IFunction.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/TDecompCholFast.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"

#include <set>
#include <array>

class PulseChiSqSNNLS {
public:

typedef BXVector::Index Index;

PulseChiSqSNNLS();
~PulseChiSqSNNLS();


bool DoFit(const std::vector<double> &samples, const TMatrixDSym &samplecor, double pederr, const std::set<int> &bxs, const TVectorD &fullpulse, const TMatrixDSym &fullpulsecov);
bool DoFit(const SampleVector &samples, const SampleMatrix &samplecor, double pederr, const BXVector &bxs, const FullSampleVector &fullpulse, const FullSampleMatrix &fullpulsecov);

const TMatrixD &pulsemat() const { return _pulsemat; }
const TMatrixDSym &invcov() const { return _invcov; }
const SamplePulseMatrix &pulsemat() const { return _pulsemat; }
const SampleMatrix &invcov() const { return _invcov; }

const TVectorD &X() const { return _ampvecmin; }
const TVectorD &Errors() const { return _errvec; }
const PulseVector &X() const { return _ampvecmin; }
const PulseVector &Errors() const { return _errvec; }
const BXVector &BXs() const { return _bxsmin; }

double ChiSq() const { return _chisq; }
void disableErrorCalculation() { _computeErrors = false; }

protected:

bool Minimize(const TMatrixDSym &samplecor, double pederr, const std::set<int> &bxs, const TMatrixDSym &fullpulsecov);
bool Minimize(const SampleMatrix &samplecor, double pederr, const FullSampleMatrix &fullpulsecov);
bool NNLS();
bool updateCov(const TMatrixDSym &samplecor, double pederr, const std::set<int> &bxs, const TMatrixDSym &fullpulsecov);
bool updateCov(const SampleMatrix &samplecor, double pederr, const FullSampleMatrix &fullpulsecov);
double ComputeChiSq();
double ComputeApproxUncertainty(unsigned int ipulse);


TVectorD _sampvec;
TMatrixDSym _invcov;
TVectorD _workvec;
TMatrixD _pulsemat;
TVectorD _ampvec;
TVectorD _ampvecmin;
TVectorD _errvec;
TMatrixD _workmat;
TMatrixD _aTamat;
TVectorD _wvec;
TVectorD _aTbvec;
TVectorD _aTbcorvec;
TMatrixDSym _aPmat;
TVectorD _sPvec;
TDecompCholFast _decompP;
std::array<double,10*10> _pulsematstorage;
std::array<double,10> _ampvecstorage;
std::array<double,10> _ampvecminstorage;
std::array<double,10> _errvecstorage;
std::array<double,10*10> _workmatstorage;
std::array<double,10*10> _aTamatstorage;
std::array<double,10> _wvecstorage;
std::array<double,10> _aTbvecstorage;
std::array<double,10> _aTbcorvecstorage;
std::array<double,10*10> _aPstorage;
std::array<double,10> _sPstorage;
std::array<double,10*10> _decompPstorage;
std::set<unsigned int> _idxsP;
std::set<unsigned int> _idxsFixed;
SampleVector _sampvec;
SampleMatrix _invcov;
SamplePulseMatrix _pulsemat;
PulseVector _ampvec;
PulseVector _errvec;
PulseVector _ampvecmin;

SampleDecompLLT _covdecomp;

BXVector _bxs;
BXVector _bxsmin;
unsigned int _nP;

double _chisq;
bool _computeErrors;
Expand Down
13 changes: 0 additions & 13 deletions RecoLocalCalo/EcalRecAlgos/interface/TDecompCholFast.h
View file
Open in desktop

This file was deleted.

31 changes: 17 additions & 14 deletions RecoLocalCalo/EcalRecAlgos/src/EcalUncalibRecHitMultiFitAlgo.cc
View file
Open in desktop
Expand Up @@ -10,7 +10,7 @@ EcalUncalibRecHitMultiFitAlgo::EcalUncalibRecHitMultiFitAlgo() :
}

/// compute rechits
EcalUncalibratedRecHit EcalUncalibRecHitMultiFitAlgo::makeRecHit(const EcalDataFrame& dataFrame, const EcalPedestals::Item * aped, const EcalMGPAGainRatio * aGain, const TMatrixDSym &noisecor, const TVectorD &fullpulse, const TMatrixDSym &fullpulsecov, const std::set<int> &activeBX) {
EcalUncalibratedRecHit EcalUncalibRecHitMultiFitAlgo::makeRecHit(const EcalDataFrame& dataFrame, const EcalPedestals::Item * aped, const EcalMGPAGainRatio * aGain, const SampleMatrix &noisecor, const FullSampleVector &fullpulse, const FullSampleMatrix &fullpulsecov, const BXVector &activeBX) {

uint32_t flags = 0;

Expand All @@ -21,7 +21,7 @@ EcalUncalibratedRecHit EcalUncalibRecHitMultiFitAlgo::makeRecHit(const EcalDataF
double pedval = 0.;
double pedrms = 0.;

std::vector<double> amplitudes(nsample);
SampleVector amplitudes;
for(unsigned int iSample = 0; iSample < nsample; iSample++) {

const EcalMGPASample &sample = dataFrame.sample(iSample);
Expand Down Expand Up @@ -66,10 +66,6 @@ EcalUncalibratedRecHit EcalUncalibRecHitMultiFitAlgo::makeRecHit(const EcalDataF
}

}


std::vector<double> fitvals;
std::vector<double> fiterrs;

if(!_computeErrors) _pulsefunc.disableErrorCalculation();
bool status = _pulsefunc.DoFit(amplitudes,noisecor,pedrms,activeBX,fullpulse,fullpulsecov);
Expand All @@ -79,25 +75,32 @@ EcalUncalibratedRecHit EcalUncalibRecHitMultiFitAlgo::makeRecHit(const EcalDataF
edm::LogWarning("EcalUncalibRecHitMultiFitAlgo::makeRecHit") << "Failed Fit" << std::endl;
}

unsigned int ipulseintime = std::distance(activeBX.begin(),activeBX.find(0));
unsigned int ipulseintime = 0;
for (unsigned int ipulse=0; ipulse<_pulsefunc.BXs().rows(); ++ipulse) {
if (_pulsefunc.BXs().coeff(ipulse)==0) {
ipulseintime = ipulse;
break;
}
}

double amplitude = status ? _pulsefunc.X()[ipulseintime] : 0.;
double amperr = status ? _pulsefunc.Errors()[ipulseintime] : 0.;

double jitter = 0.;

//printf("status = %i\n",int(status));
//printf("amplitude = %5f +- %5f, chisq = %5f\n",amplitude,amperr,chisq);

EcalUncalibratedRecHit rh( dataFrame.id(), amplitude , pedval, jitter, chisq, flags );
rh.setAmplitudeError(amperr);
for (std::set<int>::const_iterator bxit = activeBX.begin(); bxit!=activeBX.end(); ++bxit) {
int ipulse = std::distance(activeBX.begin(),bxit);
if(*bxit==0) {
rh.setOutOfTimeAmplitude(*bxit+5,0.);
} else {
rh.setOutOfTimeAmplitude(*bxit+5, status ? _pulsefunc.X()[ipulse] : 0.);

for (unsigned int ipulse=0; ipulse<_pulsefunc.BXs().rows(); ++ipulse) {
int bx = _pulsefunc.BXs().coeff(ipulse);
if (bx!=0) {
rh.setOutOfTimeAmplitude(bx+5, status ? _pulsefunc.X().coeff(ipulse) : 0.);
}
}

return rh;
}