#include <HcalHFStatusBitFromDigis.h>

Public Member Functions
	HcalHFStatusBitFromDigis ()

	HcalHFStatusBitFromDigis (int recoFirstSample, int recoSamplesToAdd, const edm::ParameterSet &HFDigiTimeParams, const edm::ParameterSet &HFTimeInWindowParams)

void	hfSetFlagFromDigi (HFRecHit &hf, const HFDataFrame &digi, const HcalCoder &coder, const HcalCalibrations &calib)

void	resetTimeSamples (int firstSample, int samplesToAdd)

	~HcalHFStatusBitFromDigis ()

Private Attributes
double	coef0_

double	coef1_

double	coef2_

int	expectedPeak_

int	firstSample_

double	HFlongwindowEthresh_

std::vector< double >	HFlongwindowMaxTime_

std::vector< double >	HFlongwindowMinTime_

double	HFshortwindowEthresh_

std::vector< double >	HFshortwindowMaxTime_

std::vector< double >	HFshortwindowMinTime_

double	minthreshold_

int	recoFirstSample_

int	recoSamplesToAdd_

int	samplesToAdd_

Detailed Description

This class sets status bit in the status words for the revised CaloRecHit objets according to informatino from the digi associated to the hit.

Date:: 2010/06/04 14:13:51

Revision:: 1.11

Author: J. Temple – University of Maryland and E. Yazgan

Definition at line 20 of file HcalHFStatusBitFromDigis.h.

Constructor & Destructor Documentation

HcalHFStatusBitFromDigis::HcalHFStatusBitFromDigis ( )

Full featured constructor for HB/HE and HO (HPD-based detectors)

Definition at line 8 of file HcalHFStatusBitFromDigis.cc.

References coef0_, coef1_, coef2_, expectedPeak_, firstSample_, HFlongwindowEthresh_, HFlongwindowMaxTime_, HFlongwindowMinTime_, HFshortwindowEthresh_, HFshortwindowMaxTime_, HFshortwindowMinTime_, minthreshold_, recoFirstSample_, recoSamplesToAdd_, and samplesToAdd_.

 {
   // use simple values in default constructor
   minthreshold_=40; // minimum energy threshold (in GeV)
   recoFirstSample_=0;
   recoSamplesToAdd_=10;
   firstSample_=3;
   samplesToAdd_=4;
   expectedPeak_=4;
   // Based on Igor V's algorithm:
   //TS4/(TS3+TS4+TS5+TS6) > 0.93 - exp(-0.38275-0.012667*E)
   coef0_= 0.93;
   coef1_ = -0.38275;
   coef2_ = -0.012667;
   // Minimum energy of 10 GeV required
   HFlongwindowEthresh_=40;
   HFlongwindowMinTime_.clear();
   HFlongwindowMinTime_.push_back(-10);
   HFlongwindowMaxTime_.clear();
   HFlongwindowMaxTime_.push_back(8);
   HFshortwindowEthresh_=40;
   HFshortwindowMinTime_.clear();
   HFshortwindowMinTime_.push_back(-10);
   HFshortwindowMaxTime_.clear();
   HFshortwindowMaxTime_.push_back(8);
 }

HcalHFStatusBitFromDigis::HcalHFStatusBitFromDigis	(	int	recoFirstSample,
		int	recoSamplesToAdd,
		const edm::ParameterSet &	HFDigiTimeParams,
		const edm::ParameterSet &	HFTimeInWindowParams
	)

Definition at line 35 of file HcalHFStatusBitFromDigis.cc.

References coef0_, coef1_, coef2_, expectedPeak_, firstSample_, edm::ParameterSet::getParameter(), HFlongwindowEthresh_, HFlongwindowMaxTime_, HFlongwindowMinTime_, HFshortwindowEthresh_, HFshortwindowMaxTime_, HFshortwindowMinTime_, minthreshold_, recoFirstSample_, recoSamplesToAdd_, and samplesToAdd_.

 {
   recoFirstSample_    = recoFirstSample;
   recoSamplesToAdd_   = recoSamplesToAdd;
 
   // Specify parameters used in forming the HFDigiTime flag
   firstSample_        = HFDigiTimeParams.getParameter<int>("HFdigiflagFirstSample");
   samplesToAdd_       = HFDigiTimeParams.getParameter<int>("HFdigiflagSamplesToAdd");
   expectedPeak_       = HFDigiTimeParams.getParameter<int>("HFdigiflagExpectedPeak");
   minthreshold_       = HFDigiTimeParams.getParameter<double>("HFdigiflagMinEthreshold");
   coef0_              = HFDigiTimeParams.getParameter<double>("HFdigiflagCoef0");
   coef1_              = HFDigiTimeParams.getParameter<double>("HFdigiflagCoef1");
   coef2_              = HFDigiTimeParams.getParameter<double>("HFdigiflagCoef2");
 
   // Specify parameters used in forming HFInTimeWindow flag
   HFlongwindowMinTime_    = HFTimeInWindowParams.getParameter<std::vector<double> >("hflongMinWindowTime");
   HFlongwindowMaxTime_    = HFTimeInWindowParams.getParameter<std::vector<double> >("hflongMaxWindowTime");
   HFlongwindowEthresh_    = HFTimeInWindowParams.getParameter<double>("hflongEthresh");
   HFshortwindowMinTime_    = HFTimeInWindowParams.getParameter<std::vector<double> >("hfshortMinWindowTime");
   HFshortwindowMaxTime_    = HFTimeInWindowParams.getParameter<std::vector<double> >("hfshortMaxWindowTime");
   HFshortwindowEthresh_    = HFTimeInWindowParams.getParameter<double>("hfshortEthresh");
 }

HcalHFStatusBitFromDigis::~HcalHFStatusBitFromDigis ( )

Definition at line 61 of file HcalHFStatusBitFromDigis.cc.

61 {}

Member Function Documentation

void HcalHFStatusBitFromDigis::hfSetFlagFromDigi	(	HFRecHit &	hf,
		const HFDataFrame &	digi,
		const HcalCoder &	coder,
		const HcalCalibrations &	calib
	)

Definition at line 63 of file HcalHFStatusBitFromDigis.cc.

References HcalCoder::adc2fC(), HcalQIESample::capid(), coef0_, coef1_, coef2_, prof2calltree::cutoff, HcalDetId::depth(), CaloRecHit::energy(), funct::exp(), expectedPeak_, firstSample_, HcalCaloFlagLabels::Fraction2TS, HcalCaloFlagLabels::HFDigiTime, HcalCaloFlagLabels::HFInTimeWindow, HFlongwindowEthresh_, HFlongwindowMaxTime_, HFlongwindowMinTime_, HFshortwindowEthresh_, HFshortwindowMaxTime_, HFshortwindowMinTime_, i, HFRecHit::id(), minthreshold_, VarParsing::mult, HcalCalibrations::pedestal(), recoFirstSample_, recoSamplesToAdd_, HcalCalibrations::respcorrgain(), HFDataFrame::sample(), samplesToAdd_, CaloRecHit::setFlagField(), HFDataFrame::size(), CaloRecHit::time(), and relativeConstraints::value.

Referenced by HcalHitReconstructor::produce().

 {
   // The following 3 values are computed using the default reconstruction window (for Shuichi's algorithm)
   double maxInWindow=-10; // maximum value found in reco window
   int maxCapid=-1;
   int maxTS=-1;  // time slice where maximum is found
 
   // The following 3 values are computed only in the window (firstSample_, firstSample_ + samplesToAdD_), which may not be the same as the default reco window  (for Igor's algorithm)
   double totalCharge=0;
   double peakCharge=0;
   double RecomputedEnergy=0;
 
   CaloSamples tool;
   coder.adc2fC(digi,tool);
 
   // Compute quantities needed for HFDigiTime, Fraction2TS FLAGS
   for (int i=0;i<digi.size();++i)
     {
       int capid=digi.sample(i).capid();
       //double value = digi.sample(i).nominal_fC()-calib.pedestal(capid);
       double value = tool[i]-calib.pedestal(capid);
       // Find largest value within reconstruction window
       if (i>=recoFirstSample_ && i <recoFirstSample_+recoSamplesToAdd_)
     {
       // Find largest overall pulse within the full digi, or just the allowed window?
       if (value>maxInWindow) 
         {
           maxCapid=capid;
           maxInWindow=value;  
           maxTS=i;
         }
     }
 
       // Sum all charge within flagging window, find charge in expected peak time slice
       if (i >=firstSample_ && i < firstSample_+samplesToAdd_)
     {
       totalCharge+=value;
       RecomputedEnergy+=value*calib.respcorrgain(capid);
       if (i==expectedPeak_) peakCharge=value;
     }
     } // for (int i=0;i<digi.size();++i)
   
   // FLAG:  HcalCaloLabel::Fraction2TS
   // Shuichi's Algorithm:  Compare size of peak in digi to charge in TS immediately before peak
   int TSfrac_counter=1; 
   // get pedestals for each capid -- add 4 to each capid, and then check mod 4.
   // (This takes care of the case where max capid =0 , and capid-1 would then be negative)
   if (maxTS>0 &&
       tool[maxTS]!=calib.pedestal(maxCapid))
     TSfrac_counter=int(50*((tool[maxTS-1]-calib.pedestal((maxCapid+3)%4))/(tool[maxTS]-calib.pedestal((maxCapid+4)%4)))+1); // 6-bit counter to hold peak ratio info
   hf.setFlagField(TSfrac_counter, HcalCaloFlagLabels::Fraction2TS,6);
 
   // FLAG:  HcalCaloLabels::HFDigiTime
   // Igor's algorithm:  compare charge in peak to total charge in window
   if (RecomputedEnergy>=minthreshold_)  // don't set noise flags for cells below a given threshold
     {
       // Calculate allowed minimum value of (TS4/TS3+4+5+6):
       double cutoff=coef0_-exp(coef1_+coef2_*RecomputedEnergy);
       
       if (peakCharge/totalCharge<cutoff)
     hf.setFlagField(1,HcalCaloFlagLabels::HFDigiTime);
     }
 
   // FLAG:  HcalCaloLabels:: HFInTimeWindow
   // Timing algorithm
   if (hf.id().depth()==1)
     {
       if (hf.energy()>=HFlongwindowEthresh_)
     {
       float mult=1./hf.energy();
       float enPow=1.;
       float mintime=0;
       float maxtime=0;
       for (unsigned int i=0;i<HFlongwindowMinTime_.size();++i)
         {
           mintime+=HFlongwindowMinTime_[i]*enPow;
           maxtime+=HFlongwindowMaxTime_[i]*enPow;
           enPow*=mult;
         }
       if (hf.time()<mintime || hf.time()>maxtime)
         hf.setFlagField(1,HcalCaloFlagLabels::HFInTimeWindow);
     }
     }
   else if (hf.id().depth()==2)
     {
       if (hf.energy()>=HFshortwindowEthresh_)
     {
       float mult=1./hf.energy();
       float enPow=1.;
       float mintime=0;
       float maxtime=0;
       for (unsigned int i=0;i<HFshortwindowMinTime_.size();++i)
         {
           mintime+=HFshortwindowMinTime_[i]*enPow;
           maxtime+=HFshortwindowMaxTime_[i]*enPow;
           enPow*=mult;
         }
       if (hf.time()<mintime || hf.time()>maxtime)
         hf.setFlagField(1,HcalCaloFlagLabels::HFInTimeWindow);
     }
     }
 
   return;
 }