HcalSimpleRecAlgoImpl Namespace Reference

Functions
template<class Digi , class RecHit >
RecHit	reco (const Digi &digi, const HcalCoder &coder, const HcalCalibrations &calibs, const int ifirst, const int n, const bool slewCorrect, const bool pulseCorrect, const HcalPulseContainmentCorrection *corr, const HcalTimeSlew::BiasSetting slewFlavor, const int runnum, const bool useLeak, const AbsOOTPileupCorrection *pileupCorrection, const BunchXParameter *bxInfo, const unsigned lenInfo, const int puCorrMethod, const PulseShapeFitOOTPileupCorrection *psFitOOTpuCorr, HcalDeterministicFit *hltOOTpuCorr, PedestalSub *hltPedSub)
template<class Digi , class RecHit >
RecHit	recoHBHE (const Digi &digi, const HcalCoder &coder, const HcalCalibrations &calibs, const int ifirst, const int n, const bool slewCorrect, const bool pulseCorrect, const HcalPulseContainmentCorrection *corr, const HcalTimeSlew::BiasSetting slewFlavor, const int runnum, const bool useLeak, const AbsOOTPileupCorrection *pileupCorrection, const BunchXParameter *bxInfo, const unsigned lenInfo, const int puCorrMethod, const PulseShapeFitOOTPileupCorrection *psFitOOTpuCorr, HcalDeterministicFit *hltOOTpuCorr, PedestalSub *hltPedSub)
template<class Digi >
float	recoHFTime (const Digi &digi, const int maxI, const double amp_fC, const bool slewCorrect, double maxA, float t0, float t2)
template<class Digi >
void	removePileup (const Digi &digi, const HcalCoder &coder, const HcalCalibrations &calibs, const int ifirst, const int n, const bool pulseCorrect, const HcalPulseContainmentCorrection *corr, const AbsOOTPileupCorrection *pileupCorrection, const BunchXParameter *bxInfo, const unsigned lenInfo, double *p_maxA, double *p_ampl, double *p_uncorr_ampl, double *p_fc_ampl, int *p_nRead, int *p_maxI, bool *leakCorrApplied, float *p_t0, float *p_t2)

Function Documentation

template<class Digi , class RecHit >

RecHit HcalSimpleRecAlgoImpl::reco ( const Digi &  digi, const HcalCoder &  coder, const HcalCalibrations &  calibs, const int  ifirst, const int  n, const bool  slewCorrect, const bool  pulseCorrect, const HcalPulseContainmentCorrection *  corr, const HcalTimeSlew::BiasSetting  slewFlavor, const int  runnum, const bool  useLeak, const AbsOOTPileupCorrection *  pileupCorrection, const BunchXParameter *  bxInfo, const unsigned  lenInfo, const int  puCorrMethod, const PulseShapeFitOOTPileupCorrection *  psFitOOTpuCorr, HcalDeterministicFit *  hltOOTpuCorr, PedestalSub *  hltPedSub  )

inline

Definition at line 313 of file HcalSimpleRecAlgo.cc.

References HcalCoder::adc2fC(), HcalDeterministicFit::apply(), PulseShapeFitOOTPileupCorrection::apply(), fwrapper::cs, HcalTimeSlew::delay(), eCorr(), HcalBarrel, HcalDetId::ieta(), HcalDetId::iphi(), leakCorr(), bookConverter::max, findQualityFiles::maxI, removePileup(), setAuxEnergy(), setRawEnergy(), CaloSamples::size(), HcalDetId::subdet(), FrontierCondition_GT_autoExpress_cfi::t0, reco::t2, HcalCalibrations::timecorr(), and timeshift_ns_hbheho().

  {
    double fc_ampl =0, ampl =0, uncorr_ampl =0, m3_ampl =0, maxA = -1.e300;
    int nRead = 0, maxI = -1;
    bool leakCorrApplied = false;
    float t0 =0, t2 =0;
    float time = -9999;

// Disable method 1 inside the removePileup function this way!
// Some code in removePileup does NOT do pileup correction & to make sure maximum share of code
    const AbsOOTPileupCorrection * inputAbsOOTpuCorr = ( puCorrMethod == 1 ? pileupCorrection: 0 );

    removePileup(digi, coder, calibs, ifirst, n,
        pulseCorrect, corr, inputAbsOOTpuCorr,
        bxInfo, lenInfo, &maxA, &ampl,
        &uncorr_ampl, &fc_ampl, &nRead, &maxI,
        &leakCorrApplied, &t0, &t2);
      
    if (maxI > 0 && maxI < (nRead - 1))
    {
      // Handle negative excursions by moving "zero":
      float minA=t0;
      if (maxA<minA) minA=maxA;
      if (t2<minA)   minA=t2;
      if (minA<0) { maxA-=minA; t0-=minA; t2-=minA; } // positivizes all samples
      
      float wpksamp = (t0 + maxA + t2);
      if (wpksamp!=0) wpksamp=(maxA + 2.0*t2) / wpksamp; 
      time = (maxI - digi.presamples())*25.0 + timeshift_ns_hbheho(wpksamp);
      
      if (slewCorrect) time-=HcalTimeSlew::delay(std::max(1.0,fc_ampl),slewFlavor);
      
      time=time-calibs.timecorr(); // time calibration
    }

    // Note that uncorr_ampl is always set from outside of method 2!
    if( puCorrMethod == 2 ){
      std::vector<double> correctedOutput;

      CaloSamples cs;
      coder.adc2fC(digi,cs);
      std::vector<int> capidvec;
      for(int ip=0; ip<cs.size(); ip++){
        const int capid = digi[ip].capid();
        capidvec.push_back(capid);
      }
      psFitOOTpuCorr->apply(cs, capidvec, calibs, correctedOutput);
      if( correctedOutput.back() == 0 && correctedOutput.size() >1 ){
    time = correctedOutput[1]; ampl = correctedOutput[0];
      }
    }
    
    // S. Brandt - Feb 19th : Adding Section for HLT
    // Run "Method 3" all the time.
    {
      std::vector<double> hltCorrOutput;

      CaloSamples cs;
      coder.adc2fC(digi,cs);
      std::vector<int> capidvec;
      for(int ip=0; ip<cs.size(); ip++){
        const int capid = digi[ip].capid();
        capidvec.push_back(capid);
      }
      hltOOTpuCorr->apply(cs, capidvec, calibs, digi, hltCorrOutput);
      if( hltCorrOutput.size() > 1 ){
        m3_ampl = hltCorrOutput[0];
        if (puCorrMethod == 3) {
      time = hltCorrOutput[1]; ampl = hltCorrOutput[0];
        }
      }
    }

    // Temporary hack to apply energy-dependent corrections to some HB- cells
    if (runnum > 0) {
      const HcalDetId& cell = digi.id();
      if (cell.subdet() == HcalBarrel) {
        const int ieta = cell.ieta();
        const int iphi = cell.iphi();
        uncorr_ampl *= eCorr(ieta, iphi, uncorr_ampl, runnum);
        ampl *= eCorr(ieta, iphi, ampl, runnum);
        m3_ampl *= eCorr(ieta, iphi, m3_ampl, runnum);
      }
    }

    // Correction for a leak to pre-sample
    if(useLeak && !leakCorrApplied) {
      uncorr_ampl *= leakCorr(uncorr_ampl);
      if (puCorrMethod < 2)
        ampl *= leakCorr(ampl); 
    }

    RecHit rh(digi.id(),ampl,time);
    setRawEnergy(rh, static_cast<float>(uncorr_ampl));
    setAuxEnergy(rh, static_cast<float>(m3_ampl));
    return rh;
  }

template<class Digi , class RecHit >

RecHit HcalSimpleRecAlgoImpl::recoHBHE ( const Digi &  digi, const HcalCoder &  coder, const HcalCalibrations &  calibs, const int  ifirst, const int  n, const bool  slewCorrect, const bool  pulseCorrect, const HcalPulseContainmentCorrection *  corr, const HcalTimeSlew::BiasSetting  slewFlavor, const int  runnum, const bool  useLeak, const AbsOOTPileupCorrection *  pileupCorrection, const BunchXParameter *  bxInfo, const unsigned  lenInfo, const int  puCorrMethod, const PulseShapeFitOOTPileupCorrection *  psFitOOTpuCorr, HcalDeterministicFit *  hltOOTpuCorr, PedestalSub *  hltPedSub  )

inline

Definition at line 419 of file HcalSimpleRecAlgo.cc.

References HcalCoder::adc2fC(), HcalDeterministicFit::apply(), PulseShapeFitOOTPileupCorrection::apply(), fwrapper::cs, HcalTimeSlew::delay(), eCorr(), HcalCaloFlagLabels::HBHEPulseFitBit, HcalBarrel, HcalDetId::ieta(), HcalDetId::iphi(), leakCorr(), bookConverter::max, findQualityFiles::maxI, removePileup(), setAuxEnergy(), CaloRecHit::setFlagField(), setRawEnergy(), CaloSamples::size(), HcalDetId::subdet(), FrontierCondition_GT_autoExpress_cfi::t0, reco::t2, HcalCalibrations::timecorr(), and timeshift_ns_hbheho().

  {
    double fc_ampl =0, ampl =0, uncorr_ampl =0, m3_ampl =0, maxA = -1.e300;
    int nRead = 0, maxI = -1;
    bool leakCorrApplied = false;
    float t0 =0, t2 =0;
    float time = -9999;
    bool useTriple = false;
    
    // Disable method 1 inside the removePileup function this way!
    // Some code in removePileup does NOT do pileup correction & to make sure maximum share of code
    const AbsOOTPileupCorrection * inputAbsOOTpuCorr = ( puCorrMethod == 1 ? pileupCorrection: 0 );
    
    removePileup(digi, coder, calibs, ifirst, n,
         pulseCorrect, corr, inputAbsOOTpuCorr,
         bxInfo, lenInfo, &maxA, &ampl,
         &uncorr_ampl, &fc_ampl, &nRead, &maxI,
         &leakCorrApplied, &t0, &t2);
    
    if (maxI > 0 && maxI < (nRead - 1))
      {
    // Handle negative excursions by moving "zero":
    float minA=t0;
    if (maxA<minA) minA=maxA;
    if (t2<minA)   minA=t2;
    if (minA<0) { maxA-=minA; t0-=minA; t2-=minA; } // positivizes all samples
    
    float wpksamp = (t0 + maxA + t2);
    if (wpksamp!=0) wpksamp=(maxA + 2.0*t2) / wpksamp; 
    time = (maxI - digi.presamples())*25.0 + timeshift_ns_hbheho(wpksamp);
    
    if (slewCorrect) time-=HcalTimeSlew::delay(std::max(1.0,fc_ampl),slewFlavor);
    
    time=time-calibs.timecorr(); // time calibration
      }
    
    // Note that uncorr_ampl is always set from outside of method 2!
    if( puCorrMethod == 2 ){
      std::vector<double> correctedOutput;

      CaloSamples cs;
      coder.adc2fC(digi,cs);
      std::vector<int> capidvec;
      for(int ip=0; ip<cs.size(); ip++){
    const int capid = digi[ip].capid();
    capidvec.push_back(capid);
      }
      psFitOOTpuCorr->apply(cs, capidvec, calibs, correctedOutput);
      if( correctedOutput.back() == 0 && correctedOutput.size() >1 ){
    time = correctedOutput[1]; ampl = correctedOutput[0]; 
    useTriple = correctedOutput[4];
      }
    }
    
    // S. Brandt - Feb 19th : Adding Section for HLT
    // Run "Method 3" all the time.
    {
      std::vector<double> hltCorrOutput;

      CaloSamples cs;
      coder.adc2fC(digi,cs);
      std::vector<int> capidvec;
      for(int ip=0; ip<cs.size(); ip++){
    const int capid = digi[ip].capid();
    capidvec.push_back(capid);
      }
      hltOOTpuCorr->apply(cs, capidvec, calibs, digi, hltCorrOutput);
      if (hltCorrOutput.size() > 1) {
        m3_ampl = hltCorrOutput[0];
        if (puCorrMethod == 3) {
      time = hltCorrOutput[1]; ampl = hltCorrOutput[0];
        }
      }
    }

    // Temporary hack to apply energy-dependent corrections to some HB- cells
    if (runnum > 0) {
      const HcalDetId& cell = digi.id();
      if (cell.subdet() == HcalBarrel) {
    const int ieta = cell.ieta();
    const int iphi = cell.iphi();
    uncorr_ampl *= eCorr(ieta, iphi, uncorr_ampl, runnum);
        ampl *= eCorr(ieta, iphi, ampl, runnum);
        m3_ampl *= eCorr(ieta, iphi, m3_ampl, runnum);
      }
    }

    // Correction for a leak to pre-sample
    if(useLeak && !leakCorrApplied) {
      uncorr_ampl *= leakCorr(uncorr_ampl); 
      if (puCorrMethod < 2)
        ampl *= leakCorr(ampl); 
    }

    HBHERecHit rh(digi.id(),ampl,time);
    if(useTriple)
      {
    rh.setFlagField(1, HcalCaloFlagLabels::HBHEPulseFitBit);
      }
    setRawEnergy(rh, static_cast<float>(uncorr_ampl));
    setAuxEnergy(rh, static_cast<float>(m3_ampl));
    return rh;
  }

template<class Digi >

float HcalSimpleRecAlgoImpl::recoHFTime ( const Digi &  digi, const int  maxI, const double  amp_fC, const bool  slewCorrect, double  maxA, float  t0, float  t2  )

inline

Definition at line 129 of file HcalSimpleRecAlgo.cc.

References alignCSCRings::e, create_public_lumi_plots::exp, f, min(), reco::t2, and timeshift_ns_hf().

Referenced by HcalSimpleRecAlgo::reconstruct(), HcalSimpleRecAlgo::reconstructHFUpgrade(), and HcalSimpleRecAlgo::reconstructQIE10().

  {
    // Handle negative excursions by moving "zero":
    float zerocorr=std::min(t0,t2);
    if (zerocorr<0.f) {
      t0   -= zerocorr;
      t2   -= zerocorr;
      maxA -= zerocorr;
    }
    
    // pair the peak with the larger of the two neighboring time samples
    float wpksamp=0.f;
    if (t0>t2) {
      wpksamp = t0+maxA;
      if (wpksamp != 0.f) wpksamp = maxA/wpksamp;
    } else {
      wpksamp = maxA+t2;
      if (wpksamp != 0.f) wpksamp = 1.+(t2/wpksamp);
    }

    float time = (maxI - digi.presamples())*25.0 + timeshift_ns_hf(wpksamp);

    if (slewCorrect && amp_fC > 0.0) {
      // -5.12327 - put in calibs.timecorr()
      double tslew=exp(0.337681-5.94689e-4*amp_fC)+exp(2.44628-1.34888e-2*amp_fC);
      time -= (float)tslew;
    }

    return time;
  }

template<class Digi >

void HcalSimpleRecAlgoImpl::removePileup ( const Digi &  digi, const HcalCoder &  coder, const HcalCalibrations &  calibs, const int  ifirst, const int  n, const bool  pulseCorrect, const HcalPulseContainmentCorrection *  corr, const AbsOOTPileupCorrection *  pileupCorrection, const BunchXParameter *  bxInfo, const unsigned  lenInfo, double *  p_maxA, double *  p_ampl, double *  p_uncorr_ampl, double *  p_fc_ampl, int *  p_nRead, int *  p_maxI, bool *  leakCorrApplied, float *  p_t0, float *  p_t2  )

inline

Definition at line 163 of file HcalSimpleRecAlgo.cc.

References HcalCoder::adc2fC(), AbsOOTPileupCorrection::apply(), fwrapper::cs, HcalPulseContainmentCorrection::getCorrection(), i, AbsOOTPileupCorrection::inputIsEnergy(), LogDebug, bookConverter::max, findQualityFiles::maxI, CaloSamples::MAXSAMPLES, min(), gen::n, HcalCalibrations::pedestal(), HcalCalibrations::respcorrgain(), and CaloSamples::size().

Referenced by reco(), recoHBHE(), HcalSimpleRecAlgo::reconstruct(), HcalSimpleRecAlgo::reconstructHFUpgrade(), and HcalSimpleRecAlgo::reconstructQIE10().

  {
    CaloSamples cs;
    coder.adc2fC(digi,cs);
    const int nRead = cs.size();
    const int iStop = std::min(nRead, n + ifirst);

    // Signal energy will be calculated both with
    // and without OOT pileup corrections. Try to
    // arrange the calculations so that we do not
    // repeat them.
    double uncorrectedEnergy[CaloSamples::MAXSAMPLES] {}, buf[CaloSamples::MAXSAMPLES] {};
    double* correctedEnergy = 0;
    double fc_ampl = 0.0, corr_fc_ampl = 0.0;
    bool pulseShapeCorrApplied = false, readjustTiming = false;
    *leakCorrApplied = false;

    if (pileupCorrection)
    {
        correctedEnergy = &buf[0];

        double correctionInput[CaloSamples::MAXSAMPLES] {};
        double gains[CaloSamples::MAXSAMPLES] {};

        for (int i=0; i<nRead; ++i)
        {
            const int capid = digi[i].capid();
            correctionInput[i] = cs[i] - calibs.pedestal(capid);
            gains[i] = calibs.respcorrgain(capid);
        }

        for (int i=ifirst; i<iStop; ++i)
            fc_ampl += correctionInput[i];

        const bool useGain = pileupCorrection->inputIsEnergy();
        for (int i=0; i<nRead; ++i)
        {
            uncorrectedEnergy[i] = correctionInput[i]*gains[i];
            if (useGain)
                correctionInput[i] = uncorrectedEnergy[i];
        }

        pileupCorrection->apply(digi.id(), correctionInput, nRead,
                                bxInfo, lenInfo, ifirst, n,
                                correctedEnergy, CaloSamples::MAXSAMPLES,
                                &pulseShapeCorrApplied, leakCorrApplied,
                                &readjustTiming);
        if (useGain)
        {
            // Gain factors have been already applied.
            // Divide by them for accumulating corr_fc_ampl.
            for (int i=ifirst; i<iStop; ++i)
                if (gains[i])
                    corr_fc_ampl += correctedEnergy[i]/gains[i];
        }
        else
        {
            for (int i=ifirst; i<iStop; ++i)
                corr_fc_ampl += correctedEnergy[i];
            for (int i=0; i<nRead; ++i)
                correctedEnergy[i] *= gains[i];
        }
    }
    else
    {
        correctedEnergy = &uncorrectedEnergy[0];

        // In this situation, we do not need to process all time slices
        const int istart = std::max(ifirst - 1, 0);
        const int iend = std::min(n + ifirst + 1, nRead);
        for (int i=istart; i<iend; ++i)
        {
            const int capid = digi[i].capid();
            float ta = cs[i] - calibs.pedestal(capid);
            if (i >= ifirst && i < iStop)
                fc_ampl += ta;
            ta *= calibs.respcorrgain(capid);
            uncorrectedEnergy[i] = ta;
        }
        corr_fc_ampl = fc_ampl;
    }

    // Uncorrected and corrected energies
    double ampl = 0.0, corr_ampl = 0.0;
    for (int i=ifirst; i<iStop; ++i)
    {
        ampl += uncorrectedEnergy[i];
        corr_ampl += correctedEnergy[i];
    }

    // Apply phase-based amplitude correction:
    if (corr && pulseCorrect)
    {
        ampl *= corr->getCorrection(fc_ampl);
        if (pileupCorrection)
        {
            if (!pulseShapeCorrApplied)
                corr_ampl *= corr->getCorrection(corr_fc_ampl);
        }
        else
            corr_ampl = ampl;
    }

    // Which energies we want to use for timing?
    const double *etime = readjustTiming ? &correctedEnergy[0] : &uncorrectedEnergy[0];
    int maxI = -1; double maxA = -1.e300;
    for (int i=ifirst; i<iStop; ++i)
        if (etime[i] > maxA)
        {
            maxA = etime[i];
            maxI = i;
        }

    // Fill out the output
    *p_maxA = maxA;
    *p_ampl = corr_ampl;
    *p_uncorr_ampl = ampl;
    *p_fc_ampl = readjustTiming ? corr_fc_ampl : fc_ampl;
    *p_nRead = nRead;
    *p_maxI = maxI;

    if (maxI <= 0 || maxI >= (nRead-1))
    {
      LogDebug("HCAL Pulse") << "HcalSimpleRecAlgoImpl::removePileup :" 
                           << " Invalid max amplitude position, " 
                           << " max Amplitude: " << maxI
                           << " first: " << ifirst
                           << " last: " << ifirst + n
                           << std::endl;
      *p_t0 = 0.f;
      *p_t2 = 0.f;
    }
    else
    {
      *p_t0 = etime[maxI - 1];
      *p_t2 = etime[maxI + 1];
    }
  }