EcalUncalibRecHitTimeWeightsAlgo.h

Go to the documentation of this file.

#ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimeWeightsAlgo_HH
#define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitTimeWeightsAlgo_HH

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "CondFormats/EcalObjects/interface/EcalPedestals.h"
#include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
#include "CondFormats/EcalObjects/interface/EcalWeightSet.h"

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


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

template<class C> class EcalUncalibRecHitTimeWeightsAlgo 
{
 public:

  EcalUncalibRecHitTimeWeightsAlgo<C>() { };
  virtual ~EcalUncalibRecHitTimeWeightsAlgo<C>() { };

  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;
  
    double maxamplitude = -std::numeric_limits<double>::max();
  
    double pulsenorm = 0.;
    int iGainSwitch = 0;

    ROOT::Math::SVector<double,nsample> pedSubSamples;
    for(unsigned int iSample = 0; iSample < nsample; iSample++) {
    
      const EcalMGPASample &sample = dataFrame.sample(iSample);
    
      double amplitude = 0.;
      int gainId = sample.gainId();
    
      double pedestal = 0.;
      double gainratio = 1.;
        
      if (gainId==0 || gainId==3) {
        pedestal = aped->mean_x1;
        gainratio = aGain->gain6Over1()*aGain->gain12Over6();
        iGainSwitch = 1;
      }
      else if (gainId==1) {
        pedestal = aped->mean_x12;
        gainratio = 1.;
        iGainSwitch = 0;
      }
      else if (gainId==2) {
        pedestal = aped->mean_x6;
        gainratio = aGain->gain12Over6();
        iGainSwitch = 1;
      }

      amplitude = ((double)(sample.adc()) - pedestal) * gainratio;
    
      if (gainId == 0) {
        //saturation
        amplitude = (4095. - pedestal) * gainratio;
      }
        
      pedSubSamples(iSample) = amplitude;

      if (amplitude>maxamplitude) {
        maxamplitude = amplitude;
      }    
      pulsenorm += fullpulse(iSample);
    }

    std::vector<double>::const_iterator amplit;
    for(amplit=amplitudes.begin(); amplit<amplitudes.end(); ++amplit) {
      int ipulse = std::distance(amplitudes.begin(),amplit);
      int bx = ipulse - 5;
      int firstsamplet = std::max(0,bx + 3);
      int offset = 7-3-bx;

      TVectorD pulse;
      pulse.ResizeTo(nsample);
      for (unsigned int isample = firstsamplet; isample<nsample; ++isample) {
        pulse(isample) = fullpulse(isample+offset);
        pedSubSamples(isample) = std::max(0., pedSubSamples(isample) - amplitudes[ipulse]*pulse(isample)/pulsenorm);
      }
    }

    // Compute parameters
    double amplitude_(-1.), jitter_(-1.);
    ROOT::Math::SVector <double,3> param = (*(weights[iGainSwitch])) * pedSubSamples;
    amplitude_ = param(EcalUncalibRecHitRecAbsAlgo<C>::iAmplitude);
    if (amplitude_) jitter_ = -param(EcalUncalibRecHitRecAbsAlgo<C>::iTime) / amplitude_;
    else jitter_ = 0.;
  
    return jitter_;
  }

};
#endif