#include <ESRecHitSimAlgo.h>

Public Member Functions
EcalRecHit	reconstruct (const ESDataFrame &digi) const

void	setAngleCorrectionFactors (const ESAngleCorrectionFactors *ang)

void	setChannelStatus (const ESChannelStatus *status)

void	setESGain (float value)

void	setIntercalibConstants (const ESIntercalibConstants *mips)

void	setMIPGeV (float value)

void	setPedestals (const ESPedestals *peds)

void	setRatioCuts (const ESRecHitRatioCuts *ratioCuts)

void	setW0 (float value)

void	setW1 (float value)

void	setW2 (float value)

Private Member Functions
EcalRecHit::ESFlags	evalAmplitude (float *result, const ESDataFrame &digi, float ped) const

double *	oldEvalAmplitude (const ESDataFrame &digi, const double &ped, const double &w0, const double &w1, const double &w2) const

EcalRecHit	oldreconstruct (const ESDataFrame &digi) const

Private Attributes
const ESAngleCorrectionFactors *	ang_

const ESChannelStatus *	channelStatus_

int	gain_

float	MIPGeV_

const ESIntercalibConstants *	mips_

const ESPedestals *	peds_

const ESRecHitRatioCuts *	ratioCuts_

float	w0_

float	w1_

float	w2_

Detailed Description

Definition at line 12 of file ESRecHitSimAlgo.h.

Member Function Documentation

◆ evalAmplitude()

EcalRecHit::ESFlags ESRecHitSimAlgo::evalAmplitude	(	float *	result,
		const ESDataFrame &	digi,
		float	ped
	)		const

private

Definition at line 12 of file ESRecHitSimAlgo.cc.

                                                                                                          {
   float energy = 0;
   float adc[3]{};
   float pw[3]{};
   pw[0] = w0_;
   pw[1] = w1_;
   pw[2] = w2_;
  
   for (int i = 0; i < digi.size(); i++) {
     energy += pw[i] * (digi.sample(i).adc() - ped);
     LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : Digi " << i << " ADC counts " << digi.sample(i).adc() << " Ped "
                                 << ped;
     //std::cout<<i<<" "<<digi.sample(i).adc()<<" "<<ped<<" "<<pw[i]<<std::endl;
     adc[i] = digi.sample(i).adc() - ped;
   }
  
   EcalRecHit::ESFlags status = EcalRecHit::kESGood;
   if (adc[0] > 20.f)
     status = EcalRecHit::kESTS13Sigmas;  // 14;
   if (adc[1] <= 0 || adc[2] <= 0)
     status = EcalRecHit::kESTS3Negative;  // 10;
   if (adc[0] > adc[1] && adc[0] > adc[2])
     status = EcalRecHit::kESTS1Largest;  // 8;
   if (adc[2] > adc[1] && adc[2] > adc[0])
     status = EcalRecHit::kESTS3Largest;  // 9;
   auto r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99.f;
   auto r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99.f;
   if (r12 > ratioCuts_->getR12High())
     status = EcalRecHit::kESBadRatioFor12;  // 5;
   if (r23 > ratioCuts_->getR23High())
     status = EcalRecHit::kESBadRatioFor23Upper;  // 6;
   if (r23 < ratioCuts_->getR23Low())
     status = EcalRecHit::kESBadRatioFor23Lower;  // 7
  
   auto A1 = adc[1];
   auto A2 = adc[2];
  
   // t0 from analytical formula:
   constexpr float n = 1.798;
   constexpr float w = 0.07291;
   constexpr float DeltaT = 25.;
   auto aaa = (A2 > 0 && A1 > 0) ? std::log(A2 / A1) / n : 20.f;  // if A1=0, t0=20
   constexpr float bbb = w / n * DeltaT;
   auto ccc = std::exp(aaa + bbb);
  
   auto t0 = (2.f - ccc) / (1.f - ccc) * DeltaT - 5.f;
  
   // A from analytical formula:
   constexpr float t1 = 20.;
 #if defined(__clang__) || defined(__INTEL_COMPILER)
   const float A_1 = 1. / (std::pow(w / n * (t1), n) * std::exp(n - w * (t1)));
 #else
   constexpr float A_1 = 1. / (std::pow(w / n * (t1), n) * std::exp(n - w * (t1)));
 #endif
   auto AA1 = A1 * A_1;
  
   if (adc[1] > 2800.f && adc[2] > 2800.f)
     status = EcalRecHit::kESSaturated;
   else if (adc[1] > 2800.f)
     status = EcalRecHit::kESTS2Saturated;
   else if (adc[2] > 2800.f)
     status = EcalRecHit::kESTS3Saturated;
  
   results[0] = energy;  // energy with weight method
   results[1] = t0;      // timing
   results[2] = AA1;     // energy with analytic method
  
   return status;
 }

References ESSample::adc(), gpuClustering::adc, HCALHighEnergyHPDFilter_cfi::energy, JetChargeProducer_cfi::exp, f, ESRecHitRatioCuts::getR12High(), ESRecHitRatioCuts::getR23High(), mps_fire::i, EcalRecHit::kESBadRatioFor12, EcalRecHit::kESBadRatioFor23Lower, EcalRecHit::kESBadRatioFor23Upper, EcalRecHit::kESGood, EcalRecHit::kESSaturated, EcalRecHit::kESTS13Sigmas, EcalRecHit::kESTS1Largest, EcalRecHit::kESTS2Saturated, EcalRecHit::kESTS3Largest, EcalRecHit::kESTS3Negative, EcalRecHit::kESTS3Saturated, dqm-mbProfile::log, LogDebug, dqmiodumpmetadata::n, funct::pow(), ratioCuts_, bookConverter::results, ESDataFrame::sample(), ESDataFrame::size(), mps_update::status, FrontierCondition_GT_autoExpress_cfi::t0, RandomServiceHelper::t1, w, w0_, w1_, and w2_.

Referenced by reconstruct().

◆ oldEvalAmplitude()

double * ESRecHitSimAlgo::oldEvalAmplitude	(	const ESDataFrame &	digi,
		const double &	ped,
		const double &	w0,
		const double &	w1,
		const double &	w2
	)		const

private

Definition at line 156 of file ESRecHitSimAlgo.cc.

                                                                                                             {
   double* results = new double[4]{};
   float energy = 0;
   double adc[3]{};
   float pw[3]{};
   pw[0] = w0;
   pw[1] = w1;
   pw[2] = w2;
  
   for (int i = 0; i < digi.size(); i++) {
     energy += pw[i] * (digi.sample(i).adc() - ped);
     LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : Digi " << i << " ADC counts " << digi.sample(i).adc() << " Ped "
                                 << ped;
     //std::cout<<i<<" "<<digi.sample(i).adc()<<" "<<ped<<" "<<pw[i]<<std::endl;
     adc[i] = digi.sample(i).adc() - ped;
   }
  
   double status = 0;
   if (adc[0] > 20)
     status = 14;
   if (adc[1] <= 0 || adc[2] <= 0)
     status = 10;
   if (adc[0] > adc[1] && adc[0] > adc[2])
     status = 8;
   if (adc[2] > adc[1] && adc[2] > adc[0])
     status = 9;
   double r12 = (adc[1] != 0) ? adc[0] / adc[1] : 99;
   double r23 = (adc[2] != 0) ? adc[1] / adc[2] : 99;
   if (r12 > ratioCuts_->getR12High())
     status = 5;
   if (r23 > ratioCuts_->getR23High())
     status = 6;
   if (r23 < ratioCuts_->getR23Low())
     status = 7;
  
   double A1 = adc[1];
   double A2 = adc[2];
  
   // t0 from analytical formula:
   double n = 1.798;
   double w = 0.07291;
   double DeltaT = 25.;
   double aaa = (A2 > 0 && A1 > 0) ? log(A2 / A1) / n : 20.;  // if A1=0, t0=20
   double bbb = w / n * DeltaT;
   double ccc = exp(aaa + bbb);
  
   double t0 = (2. - ccc) / (1. - ccc) * DeltaT - 5;
  
   // A from analytical formula:
   double t1 = 20.;
   double A_1 = pow(w / n * (t1), n) * exp(n - w * (t1));
   double AA1 = (A_1 != 0.) ? A1 / A_1 : 0.;
  
   if (adc[1] > 2800 && adc[2] > 2800)
     status = 11;
   else if (adc[1] > 2800)
     status = 12;
   else if (adc[2] > 2800)
     status = 13;
  
   results[0] = energy;  // energy with weight method
   results[1] = t0;      // timing
   results[2] = status;  // hit status
   results[3] = AA1;     // energy with analytic method
  
   return results;
 }

References ESSample::adc(), gpuClustering::adc, HCALHighEnergyHPDFilter_cfi::energy, JetChargeProducer_cfi::exp, ESRecHitRatioCuts::getR12High(), ESRecHitRatioCuts::getR23High(), mps_fire::i, dqm-mbProfile::log, LogDebug, dqmiodumpmetadata::n, funct::pow(), ratioCuts_, bookConverter::results, ESDataFrame::sample(), ESDataFrame::size(), mps_update::status, FrontierCondition_GT_autoExpress_cfi::t0, RandomServiceHelper::t1, w, and w2.

Referenced by oldreconstruct().

◆ oldreconstruct()

EcalRecHit ESRecHitSimAlgo::oldreconstruct ( const ESDataFrame & digi ) const

private

Definition at line 225 of file ESRecHitSimAlgo.cc.

                                                                         {
   ESPedestals::const_iterator it_ped = peds_->find(digi.id());
  
   ESIntercalibConstantMap::const_iterator it_mip = mips_->getMap().find(digi.id());
   ESAngleCorrectionFactors::const_iterator it_ang = ang_->getMap().find(digi.id());
  
   ESChannelStatusMap::const_iterator it_status = channelStatus_->getMap().find(digi.id());
  
   double* results;
  
   results = oldEvalAmplitude(digi, it_ped->getMean(), w0_, w1_, w2_);
  
   double energy = results[0];
   double t0 = results[1];
   int status = (int)results[2];
   double otenergy = results[3] * 1000000.;  // set out-of-time energy to keV
   delete[] results;
  
   double mipCalib = (fabs(cos(*it_ang)) != 0.) ? (*it_mip) / fabs(cos(*it_ang)) : 0.;
   energy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
   otenergy *= (mipCalib != 0.) ? MIPGeV_ / mipCalib : 0.;
  
   LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : reconstructed energy " << energy;
  
   EcalRecHit rechit(digi.id(), energy, t0);
   // edm: this is just a placeholder for alternative energy reconstruction,
   // so put it in the same float, with different name
   // rechit.setOutOfTimeEnergy(otenergy);
   rechit.setEnergyError(otenergy);
  
   if (it_status->getStatusCode() == 1) {
     rechit.setFlag(EcalRecHit::kESDead);
   } else {
     if (status == 0)
       rechit.setFlag(EcalRecHit::kESGood);
     else if (status == 5)
       rechit.setFlag(EcalRecHit::kESBadRatioFor12);
     else if (status == 6)
       rechit.setFlag(EcalRecHit::kESBadRatioFor23Upper);
     else if (status == 7)
       rechit.setFlag(EcalRecHit::kESBadRatioFor23Lower);
     else if (status == 8)
       rechit.setFlag(EcalRecHit::kESTS1Largest);
     else if (status == 9)
       rechit.setFlag(EcalRecHit::kESTS3Largest);
     else if (status == 10)
       rechit.setFlag(EcalRecHit::kESTS3Negative);
     else if (status == 11)
       rechit.setFlag(EcalRecHit::kESSaturated);
     else if (status == 12)
       rechit.setFlag(EcalRecHit::kESTS2Saturated);
     else if (status == 13)
       rechit.setFlag(EcalRecHit::kESTS3Saturated);
     else if (status == 14)
       rechit.setFlag(EcalRecHit::kESTS13Sigmas);
   }
  
   return rechit;
 }

References ang_, channelStatus_, funct::cos(), HCALHighEnergyHPDFilter_cfi::energy, ESCondObjectContainer< T >::find(), ESCondObjectContainer< T >::getMap(), ESDataFrame::id(), createfilelist::int, EcalRecHit::kESBadRatioFor12, EcalRecHit::kESBadRatioFor23Lower, EcalRecHit::kESBadRatioFor23Upper, EcalRecHit::kESDead, EcalRecHit::kESGood, EcalRecHit::kESSaturated, EcalRecHit::kESTS13Sigmas, EcalRecHit::kESTS1Largest, EcalRecHit::kESTS2Saturated, EcalRecHit::kESTS3Largest, EcalRecHit::kESTS3Negative, EcalRecHit::kESTS3Saturated, LogDebug, MIPGeV_, mips_, oldEvalAmplitude(), peds_, bookConverter::results, EcalRecHit::setEnergyError(), mps_update::status, FrontierCondition_GT_autoExpress_cfi::t0, w0_, w1_, and w2_.

◆ reconstruct()

EcalRecHit ESRecHitSimAlgo::reconstruct ( const ESDataFrame & digi ) const

Definition at line 82 of file ESRecHitSimAlgo.cc.

                                                                      {
   auto ind = digi.id().hashedIndex();
  
   auto const& ped = peds_->preshower(ind);
   auto const& mip = mips_->getMap().preshower(ind);
   auto const& ang = ang_->getMap().preshower(ind);
   auto const& statusCh = channelStatus_->getMap().preshower(ind);
  
   float results[3]{};
  
   auto status = evalAmplitude(results, digi, ped.getMean());
  
   auto energy = results[0];
   auto t0 = results[1];
   auto otenergy = results[2] * 1000000.f;  // set out-of-time energy to keV
  
   auto mipCalib = (mip != 0.f) ? MIPGeV_ * std::abs(vdt::fast_cosf(ang)) / (mip) : 0.f;
   energy *= mipCalib;
   otenergy *= mipCalib;
  
   LogDebug("ESRecHitSimAlgo") << "ESRecHitSimAlgo : reconstructed energy " << energy;
  
   EcalRecHit rechit(digi.id(), energy, t0);
   // edm: this is just a placeholder for alternative energy reconstruction,
   // so put it in the same float, with different name
   // rechit.setOutOfTimeEnergy(otenergy);
   rechit.setEnergyError(otenergy);
  
   rechit.setFlag(statusCh.getStatusCode() == 1 ? EcalRecHit::kESDead : status);
  
   return rechit;
 }