Merge pull request cms-sw#5321 from bendavid/ecalMultifit_optimized

Ecal multifit optimized
Resolved Conflicts:
	RecoLocalCalo/EcalRecAlgos/BuildFile.xml
	RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitProducer.h
	RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.cc
	RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.h

RecoLocalCalo/EcalRecAlgos/BuildFile.xml

@@ -9,6 +9,7 @@
 <use   name="CondFormats/ESObjects"/>
 <use   name="CondFormats/EcalObjects"/>
 <use   name="CondFormats/DataRecord"/>
+<use   name="eigen"/>
 <export>
   <lib   name="1"/>
 </export>

RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitMultiFitAlgo.h

@@ -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:

RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimeWeightsAlgo.h

@@ -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>
 
@@ -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;
 
@@ -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;
 
       TVectorD pulse;
       pulse.ResizeTo(nsample);

RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h

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

RecoLocalCalo/EcalRecAlgos/interface/PulseChiSqSNNLS.h

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

RecoLocalCalo/EcalRecAlgos/src/EcalUncalibRecHitMultiFitAlgo.cc

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