ZdcSimpleRecAlgo.cc

Go to the documentation of this file.

#include "RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
#include <algorithm> // for "max"
#include <iostream>
#include <cmath>



constexpr double MaximumFractionalError = 0.0005; // 0.05% error allowed from this source

ZdcSimpleRecAlgo::ZdcSimpleRecAlgo(bool correctForTimeslew, bool correctForPulse, float phaseNS, int recoMethod, int lowGainOffset, double lowGainFrac) : 
  recoMethod_(recoMethod),
  correctForTimeslew_(correctForTimeslew),
  correctForPulse_(correctForPulse),
  phaseNS_(phaseNS),
  lowGainOffset_(lowGainOffset),
  lowGainFrac_(lowGainFrac)
  {
}

ZdcSimpleRecAlgo::ZdcSimpleRecAlgo(int recoMethod) : 
  recoMethod_(recoMethod),
  correctForTimeslew_(false) {
}

void ZdcSimpleRecAlgo::initPulseCorr(int toadd, const HcalTimeSlew* hcalTimeSlew_delay) {
  if (correctForPulse_) {    
    pulseCorr_ = std::make_unique<HcalPulseContainmentCorrection>(toadd,phaseNS_,MaximumFractionalError, hcalTimeSlew_delay);
  }
}
//static float timeshift_ns_zdc(float wpksamp);

namespace ZdcSimpleRecAlgoImpl {
  template<class Digi, class RecHit>
  inline RecHit reco1(const Digi& digi, const HcalCoder& coder, const HcalCalibrations& calibs, 
              const std::vector<unsigned int>& myNoiseTS, const std::vector<unsigned int>& mySignalTS, int lowGainOffset, double lowGainFrac, bool slewCorrect, const HcalPulseContainmentCorrection* corr, HcalTimeSlew::BiasSetting slewFlavor) {
    CaloSamples tool;
    coder.adc2fC(digi,tool);
    int ifirst = mySignalTS[0];
    int n = mySignalTS.size();
    double ampl=0; int maxI = -1; double maxA = -1e10; double ta=0;
    double fc_ampl=0;
    double lowGEnergy=0; double lowGfc_ampl=0; double TempLGAmp=0;
// TS increment for regular energy to lowGainEnergy      
// Signal in higher TS (effective "low Gain") has a fraction of the whole signal
// This constant for fC --> GeV is dervied from 2010 PbPb analysis of single neutrons
// assumed similar fraction for EM and HAD sections
// this variable converts from current assumed TestBeam values for fC--> GeV
// to the lowGain TS region fraction value (based on 1N Had, assume EM same response)
// regular energy    
    for (int i=ifirst; i<tool.size() && i<n+ifirst; i++) {
      int capid=digi[i].capid();
      ta = (tool[i]-calibs.pedestal(capid)); // pedestal subtraction
      fc_ampl+=ta; 
      ta*= calibs.respcorrgain(capid) ; // fC --> GeV
      ampl+=ta;
      if(ta>maxA){
    maxA=ta;
    maxI=i;
      }
    }
// calculate low Gain Energy (in 2010 PbPb, signal TS 4,5,6, lowGain TS: 6,7,8) 
    int topLowGain=10;
    if((n+ifirst+lowGainOffset)<=10){
      topLowGain=n+ifirst+lowGainOffset;
    } else {
      topLowGain=10;
    }
    for (int iLG=(ifirst+lowGainOffset); iLG<tool.size() && iLG<topLowGain; iLG++) {
      int capid=digi[iLG].capid();
      TempLGAmp = (tool[iLG]-calibs.pedestal(capid)); // pedestal subtraction
      lowGfc_ampl+=TempLGAmp; 
      TempLGAmp*= calibs.respcorrgain(capid) ; // fC --> GeV
      TempLGAmp*= lowGainFrac ; // TS (signalRegion) --> TS (lowGainRegion)
      lowGEnergy+=TempLGAmp;
    }    
    double time=-9999;
    // Time based on regular energy (lowGainEnergy signal assumed to happen at same Time)
    if(maxI==0 || maxI==(tool.size()-1)) {      
      LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco1 :" 
                           << " Invalid max amplitude position, " 
                           << " max Amplitude: "<< maxI
                           << " first: "<<ifirst
                           << " last: "<<(tool.size()-1)
                           << std::endl;
    } else {
      int capid=digi[maxI-1].capid();
      double Energy0 = ((tool[maxI-1])*calibs.respcorrgain(capid) );
// if any of the energies used in the weight are negative, make them 0 instead
// these are actually QIE values, not energy
      if(Energy0<0){Energy0=0.;}
      capid=digi[maxI].capid();
      double Energy1 = ((tool[maxI])*calibs.respcorrgain(capid) ) ;
      if(Energy1<0){Energy1=0.;}
      capid=digi[maxI+1].capid();
      double Energy2 = ((tool[maxI+1])*calibs.respcorrgain(capid) );
      if(Energy2<0){Energy2=0.;}
//
      double TSWeightEnergy = ((maxI-1)*Energy0 + maxI*Energy1 + (maxI+1)*Energy2);
      double EnergySum=Energy0+Energy1+Energy2;
      double AvgTSPos=0.;
      if (EnergySum!=0) AvgTSPos=TSWeightEnergy/ EnergySum; 
// If time is zero, set it to the "nonsensical" -99
// Time should be between 75ns and 175ns (Timeslices 3-7)
      if(AvgTSPos==0){
         time=-99;
      } else {
         time = (AvgTSPos*25.0);
      }
      if (corr!=nullptr) {
        // Apply phase-based amplitude correction:
           ampl *= corr->getCorrection(fc_ampl);     
      }
    }
    return RecHit(digi.id(),ampl,time,lowGEnergy);    
  }
}

namespace ZdcSimpleRecAlgoImpl {
  template<class Digi, class RecHit>
  inline RecHit reco2(const Digi& digi, const HcalCoder& coder, const HcalCalibrations& calibs, 
             const std::vector<unsigned int>& myNoiseTS, const std::vector<unsigned int>& mySignalTS, int lowGainOffset, double lowGainFrac, bool slewCorrect, const HcalPulseContainmentCorrection* corr, HcalTimeSlew::BiasSetting slewFlavor) {
    CaloSamples tool;
    coder.adc2fC(digi,tool);
    // Reads noiseTS and signalTS from database
    int ifirst = mySignalTS[0];
//    int n = mySignalTS.size();
    double ampl=0; int maxI = -1; double maxA = -1e10; double ta=0;
    double fc_ampl=0;
    double lowGEnergy=0; double lowGfc_ampl=0; double TempLGAmp=0;
//  TS increment for regular energy to lowGainEnergy      
// Signal in higher TS (effective "low Gain") has a fraction of the whole signal
// This constant for fC --> GeV is dervied from 2010 PbPb analysis of single neutrons
// assumed similar fraction for EM and HAD sections
// this variable converts from current assumed TestBeam values for fC--> GeV
// to the lowGain TS region fraction value (based on 1N Had, assume EM same response)
    double Allnoise = 0; 
    int noiseslices = 0;
    int CurrentTS = 0;
    double noise = 0;
// regular energy (both use same noise)    
    for(unsigned int iv = 0; iv<myNoiseTS.size(); ++iv)
    {
      CurrentTS = myNoiseTS[iv];
      Allnoise += tool[CurrentTS];
      noiseslices++;
    }
    if(noiseslices != 0) {
      noise = (Allnoise)/double(noiseslices);
    } else {
      noise = 0;
    }
    for(unsigned int ivs = 0; ivs<mySignalTS.size(); ++ivs)
    {
      CurrentTS = mySignalTS[ivs];
      int capid=digi[CurrentTS].capid();
//       if(noise<0){
//       // flag hit as having negative noise, and don't subtract anything, because
//       // it will falsely increase the energy
//          noisefactor=0.;
//       } 
      ta = tool[CurrentTS]-noise;
      fc_ampl+=ta; 
      ta*= calibs.respcorrgain(capid) ; // fC --> GeV
      ampl+=ta;
      if(ta>maxA){
         maxA=ta;
         maxI=CurrentTS;
      }
    }   
// calculate low Gain Energy (in 2010 PbPb, signal TS 4,5,6, lowGain TS: 6,7,8)    
    for(unsigned int iLGvs = 0; iLGvs<mySignalTS.size(); ++iLGvs)
    {
      CurrentTS = mySignalTS[iLGvs]+lowGainOffset;
      int capid=digi[CurrentTS].capid();
      TempLGAmp = tool[CurrentTS]-noise;
      lowGfc_ampl+=TempLGAmp; 
      TempLGAmp*= calibs.respcorrgain(capid) ; // fC --> GeV
      TempLGAmp*= lowGainFrac ; // TS (signalRegion) --> TS (lowGainRegion)
      lowGEnergy+=TempLGAmp;
    }      
//    if(ta<0){
//      // flag hits that have negative energy
//    }

    double time=-9999;
    // Time based on regular energy (lowGainEnergy signal assumed to happen at same Time)
    if(maxI==0 || maxI==(tool.size()-1)) {      
      LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco2 :" 
                           << " Invalid max amplitude position, " 
                           << " max Amplitude: "<< maxI
                           << " first: "<<ifirst
                           << " last: "<<(tool.size()-1)
                           << std::endl;
    } else {
      int capid=digi[maxI-1].capid();
      double Energy0 = ((tool[maxI-1])*calibs.respcorrgain(capid) );
// if any of the energies used in the weight are negative, make them 0 instead
// these are actually QIE values, not energy
      if(Energy0<0){Energy0=0.;}
      capid=digi[maxI].capid();
      double Energy1 = ((tool[maxI])*calibs.respcorrgain(capid) ) ;
      if(Energy1<0){Energy1=0.;}
      capid=digi[maxI+1].capid();
      double Energy2 = ((tool[maxI+1])*calibs.respcorrgain(capid) );
      if(Energy2<0){Energy2=0.;}
//
      double TSWeightEnergy = ((maxI-1)*Energy0 + maxI*Energy1 + (maxI+1)*Energy2);
      double EnergySum=Energy0+Energy1+Energy2;
      double AvgTSPos=0.;
      if (EnergySum!=0) AvgTSPos=TSWeightEnergy/ EnergySum; 
// If time is zero, set it to the "nonsensical" -99
// Time should be between 75ns and 175ns (Timeslices 3-7)
      if(AvgTSPos==0){
         time=-99;
      } else {
         time = (AvgTSPos*25.0);
      }
      if (corr!=nullptr) {
    // Apply phase-based amplitude correction:
           ampl *= corr->getCorrection(fc_ampl);     
      }
    }
    return RecHit(digi.id(),ampl,time,lowGEnergy);    
  }
}

ZDCRecHit ZdcSimpleRecAlgo::reconstruct(const ZDCDataFrame& digi, const std::vector<unsigned int>& myNoiseTS, const std::vector<unsigned int>& mySignalTS, const HcalCoder& coder, const HcalCalibrations& calibs) const {

  if(recoMethod_ == 1)
   return ZdcSimpleRecAlgoImpl::reco1<ZDCDataFrame,ZDCRecHit>(digi,coder,calibs,
                                  myNoiseTS,mySignalTS,lowGainOffset_,lowGainFrac_,false,
                                  nullptr,
                                  HcalTimeSlew::Fast);
  if(recoMethod_ == 2)
   return ZdcSimpleRecAlgoImpl::reco2<ZDCDataFrame,ZDCRecHit>(digi,coder,calibs,
                                  myNoiseTS,mySignalTS,lowGainOffset_,lowGainFrac_,false,
                                  nullptr,HcalTimeSlew::Fast);

     edm::LogError("ZDCSimpleRecAlgoImpl::reconstruct, recoMethod was not declared");
     throw cms::Exception("ZDCSimpleRecoAlgos::exiting process");
   
}