HcalDeterministicFit.cc

Go to the documentation of this file.

#include <iostream>
#include <cmath>
#include <climits>
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalDeterministicFit.h"

constexpr float HcalDeterministicFit::invGpar[3];
constexpr float HcalDeterministicFit::negThresh[2];
constexpr int HcalDeterministicFit::HcalRegion[2];
constexpr float HcalDeterministicFit::rCorr[2];
constexpr float HcalDeterministicFit::rCorrSiPM[2];

using namespace std;

HcalDeterministicFit::HcalDeterministicFit() {
}

HcalDeterministicFit::~HcalDeterministicFit() { 
}

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

constexpr float HcalDeterministicFit::landauFrac[];
// Landau function integrated in 1 ns intervals
//Landau pulse shape from https://indico.cern.ch/event/345283/contribution/3/material/slides/0.pdf
//Landau turn on by default at left edge of time slice 
// normalized to 1 on [0,10000]
void HcalDeterministicFit::getLandauFrac(float tStart, float tEnd, float &sum) const{

  if (std::abs(tStart-tEnd-tsWidth)<0.1f) {
    sum=0.f;
    return;
  }
  sum= landauFrac[int(ceil(tStart+tsWidth))];
  return;
}

constexpr float HcalDeterministicFit::siPM205Frac[];
void HcalDeterministicFit::get205Frac(float tStart, float tEnd, float &sum) const{

  if (std::abs(tStart-tEnd-tsWidth)<0.1f) {
    sum=0.f;
    return;
  }
  sum= siPM205Frac[int(ceil(tStart+tsWidth))];
  return;
}

constexpr float HcalDeterministicFit::siPM206Frac[];
void HcalDeterministicFit::get206Frac(float tStart, float tEnd, float &sum) const{

  if (std::abs(tStart-tEnd-tsWidth)<0.1f) {
    sum=0.f;
    return;
  }
  sum= siPM206Frac[int(ceil(tStart+tsWidth))];
  return;
}


constexpr float HcalDeterministicFit::siPM207Frac[];
void HcalDeterministicFit::get207Frac(float tStart, float tEnd, float &sum) const{

  if (std::abs(tStart-tEnd-tsWidth)<0.1f) {
    sum=0.f;
    return;
  }
  sum= siPM207Frac[int(ceil(tStart+tsWidth))];
  return;
}

void HcalDeterministicFit::getFrac(float tStart, float tEnd, float &sum, FType fType) const{ 
  switch(fType){
    case shape205:    get205Frac(tStart,tEnd,sum);    break; 
    case shape206:    get206Frac(tStart,tEnd,sum);    break; 
    case shape207:    get207Frac(tStart,tEnd,sum);    break; 
    case shapeLandau: getLandauFrac(tStart,tEnd,sum); break; 
  }
}

void HcalDeterministicFit::phase1Apply(const HBHEChannelInfo& channelData,
                       float& reconstructedEnergy,
                       float& reconstructedTime,
                       const HcalTimeSlew* hcalTimeSlew_delay) const
{


  unsigned int soi=channelData.soi();

  std::vector<double> corrCharge;
  std::vector<double> inputCharge;
  std::vector<double> inputPedestal;
  std::vector<double> inputNoise;
  double gainCorr = 0;
  double respCorr = 0;

  for(unsigned int ip=0; ip<channelData.nSamples(); ip++){

    double charge = channelData.tsRawCharge(ip);
    double ped = channelData.tsPedestal(ip); 
    double noise = channelData.tsPedestalWidth(ip);
    double gain = channelData.tsGain(ip);

    gainCorr = gain;
    inputCharge.push_back(charge);
    inputPedestal.push_back(ped);
    inputNoise.push_back(noise);

  }

  fPedestalSubFxn_.calculate(inputCharge, inputPedestal, inputNoise, corrCharge, soi, channelData.nSamples());

  if      (fTimeSlew_==0) respCorr=1.0;
  else if (fTimeSlew_==1) channelData.hasTimeInfo()?respCorr=rCorrSiPM[0]:respCorr=rCorr[0];
  else if (fTimeSlew_==2) channelData.hasTimeInfo()?respCorr=rCorrSiPM[1]:respCorr=rCorr[1];
  else if (fTimeSlew_==3) respCorr=frespCorr_;

  float tsShift3,tsShift4,tsShift5;
  tsShift3=0.f,tsShift4=0.f,tsShift5=0.f;

  if(applyTimeSlew_) {
    tsShift3=hcalTimeSlew_delay->delay(inputCharge[soi-1], fTimeSlew_, fTimeSlewBias_, !channelData.hasTimeInfo());
    tsShift4=hcalTimeSlew_delay->delay(inputCharge[soi],   fTimeSlew_, fTimeSlewBias_, !channelData.hasTimeInfo());
    tsShift5=hcalTimeSlew_delay->delay(inputCharge[soi+1], fTimeSlew_, fTimeSlewBias_, !channelData.hasTimeInfo());
  }

  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;

  FType fType;
  if(channelData.hasTimeInfo() && channelData.recoShape()==205) fType = shape205;
  else if(channelData.hasTimeInfo() && channelData.recoShape()==206) fType = shape206;
  else if(channelData.hasTimeInfo() && channelData.recoShape()==207) fType = shape207;
  else fType = shapeLandau;

  getFrac(-tsShift3,-tsShift3+tsWidth,i3,fType);
  getFrac(-tsShift3+tsWidth,-tsShift3+tsWidth*2,n3,fType);
  getFrac(-tsShift3+tsWidth*2,-tsShift3+tsWidth*3,nn3,fType);

  getFrac(-tsShift4,-tsShift4+tsWidth,i4,fType);
  getFrac(-tsShift4+tsWidth,-tsShift4+tsWidth*2,n4,fType);

  getFrac(-tsShift5,-tsShift5+tsWidth,i5,fType);
  getFrac(-tsShift5+tsWidth,-tsShift5+tsWidth*2,n5,fType);
  
  if (i3 != 0 && i4 != 0 && i5 != 0) {

    ch3=corrCharge[soi-1]/i3;
    ch4=(i3*corrCharge[soi]-n3*corrCharge[soi-1])/(i3*i4);
    ch5=(n3*n4*corrCharge[soi-1]-i4*nn3*corrCharge[soi-1]-i3*n4*corrCharge[soi]+i3*i4*corrCharge[soi+1])/(i3*i4*i5);

    if (ch3<negThresh[0]) {
      ch3=negThresh[0];
      ch4=corrCharge[soi]/i4;
      ch5=(i4*corrCharge[soi+1]-n4*corrCharge[soi])/(i4*i5);
    }
    if (ch5<negThresh[0] && ch4>negThresh[1]) {
      double ratio = (corrCharge[soi]-ch3*i3)/(corrCharge[soi+1]-negThresh[0]*i5);
      if (ratio < 5 && ratio > 0.5) {
        double invG = invGpar[0]+invGpar[1]*std::sqrt(2*std::log(invGpar[2]/ratio));
        float iG=0.f;
    getFrac(-invG,-invG+tsWidth,iG,fType);
    if (iG != 0 ) {
      ch4=(corrCharge[soi]-ch3*n3)/(iG);
      tsShift4=invG;
    }
      }
    }
  }

  if (ch4<1) {
    ch4=0.f;
  }

  reconstructedEnergy=ch4*gainCorr*respCorr;
  reconstructedTime=tsShift4;
}