#include <EcalUncalibRecHitRatioMethodAlgo.h>
Classes | |
struct | CalculatedRecHit |
struct | Ratio |
struct | Tmax |
Public Member Functions | |
void | computeAmplitude (std::vector< double > &litudeFitParameters) |
void | computeTime (std::vector< double > &timeFitParameters, std::pair< double, double > &timeFitLimits, std::vector< double > &litudeFitParameters) |
bool | fixMGPAslew (const C &dataFrame) |
CalculatedRecHit | getCalculatedRecHit () |
void | init (const C &dataFrame, const EcalSampleMask &sampleMask, const double *pedestals, const double *pedestalRMSes, const double *gainRatios) |
virtual EcalUncalibratedRecHit | makeRecHit (const C &dataFrame, const EcalSampleMask &sampleMask, const double *pedestals, const double *pedestalRMSes, const double *gainRatios, std::vector< double > &timeFitParameters, std::vector< double > &litudeFitParameters, std::pair< double, double > &timeFitLimits) |
virtual | ~EcalUncalibRecHitRatioMethodAlgo () |
Protected Attributes | |
std::vector< double > | amplitudeErrors_ |
std::vector< double > | amplitudes_ |
double | ampMaxError_ |
CalculatedRecHit | calculatedRechit_ |
int | num_ |
double | pedestal_ |
std::vector< Ratio > | ratios_ |
EcalSampleMask | sampleMask_ |
DetId | theDetId_ |
std::vector< Tmax > | times_ |
std::vector< Tmax > | timesAB_ |
Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse using a ratio method
Definition at line 20 of file EcalUncalibRecHitRatioMethodAlgo.h.
virtual EcalUncalibRecHitRatioMethodAlgo< C >::~EcalUncalibRecHitRatioMethodAlgo | ( | ) | [inline, virtual] |
Definition at line 43 of file EcalUncalibRecHitRatioMethodAlgo.h.
{ };
void EcalUncalibRecHitRatioMethodAlgo< C >::computeAmplitude | ( | std::vector< double > & | amplitudeFitParameters | ) |
Definition at line 518 of file EcalUncalibRecHitRatioMethodAlgo.h.
References beta, alignCSCRings::e, funct::exp(), f, i, and funct::log().
Referenced by EcalUncalibRecHitWorkerGlobal::run().
{ // // // CALCULATE AMPLITUDE // // // double alpha = amplitudeFitParameters[0]; double beta = amplitudeFitParameters[1]; // calculate pedestal, again double pedestalLimit = calculatedRechit_.timeMax - (alpha * beta) - 1.0; double sumA = 0; double sumF = 0; double sumAF = 0; double sumFF = 0; double sum1 = 0; for (unsigned int i = 0; i < amplitudes_.size(); i++) { double err2 = amplitudeErrors_[i]*amplitudeErrors_[i]; double f = 0; double termOne = 1 + (i - calculatedRechit_.timeMax) / (alpha * beta); if (termOne > 1.e-5) f = exp(alpha * log(termOne)) * exp(-(i - calculatedRechit_.timeMax) / beta); // apply range of interesting samples if ( (i < pedestalLimit) || (f > 0.6 && i <= calculatedRechit_.timeMax) || (f > 0.4 && i >= calculatedRechit_.timeMax)) { sum1 += 1/err2; sumA += amplitudes_[i]/err2; sumF += f/err2; sumAF += f*amplitudes_[i]/err2; sumFF += f*f/err2; } } calculatedRechit_.amplitudeMax = 0; if(sum1 > 0){ double denom = sumFF*sum1 - sumF*sumF; if(fabs(denom)>1.0e-20){ calculatedRechit_.amplitudeMax = (sumAF*sum1 - sumA*sumF)/denom; } } }
void EcalUncalibRecHitRatioMethodAlgo< C >::computeTime | ( | std::vector< double > & | timeFitParameters, |
std::pair< double, double > & | timeFitLimits, | ||
std::vector< double > & | amplitudeFitParameters | ||
) |
Definition at line 214 of file EcalUncalibRecHitRatioMethodAlgo.h.
References alpha, beta, alignCSCRings::e, funct::exp(), f, i, getHLTprescales::index, EcalUncalibRecHitRatioMethodAlgo< C >::Tmax::index, j, gen::k, funct::log(), evf::evtn::offset(), Gflash::Rmax, Gflash::Rmin, mathSSE::sqrt(), testRegression::stat, launcher::step, tmax, and relativeConstraints::value.
Referenced by EcalUncalibRecHitWorkerGlobal::run().
{ // // // RATIO METHOD FOR TIME STARTS HERE // // // double ampMaxAlphaBeta = 0; double tMaxAlphaBeta = 5; double tMaxErrorAlphaBeta = 999; double tMaxRatio = 5; double tMaxErrorRatio = 999; double sumAA = 0; double sumA = 0; double sum1 = 0; double sum0 = 0; double sumAf = 0; double sumff = 0; double NullChi2 = 0; // null hypothesis = no pulse, pedestal only for(unsigned int i = 0; i < amplitudes_.size(); i++){ double err2 = amplitudeErrors_[i]*amplitudeErrors_[i]; sum0 += 1; sum1 += 1/err2; sumA += amplitudes_[i]/err2; sumAA += amplitudes_[i]*amplitudes_[i]/err2; } if(sum0>0){ NullChi2 = (sumAA - sumA*sumA/sum1)/sum0; }else{ // not enough samples to reconstruct the pulse return; } // Make all possible Ratio's based on any pair of samples i and j // (j>i) with positive amplitudes_ // // Ratio[k] = Amp[i]/Amp[j] // where Amp[i] is pedestal subtracted ADC value in a time sample [i] // double alphabeta = amplitudeFitParameters[0]*amplitudeFitParameters[1]; double alpha = amplitudeFitParameters[0]; double beta = amplitudeFitParameters[1]; for(unsigned int i = 0; i < amplitudes_.size()-1; i++){ for(unsigned int j = i+1; j < amplitudes_.size(); j++){ if(amplitudes_[i]>1 && amplitudes_[j]>1){ // ratio double Rtmp = amplitudes_[i]/amplitudes_[j]; // error^2 due to stat fluctuations of time samples // (uncorrelated for both samples) double err1 = Rtmp*Rtmp*( (amplitudeErrors_[i]*amplitudeErrors_[i]/(amplitudes_[i]*amplitudes_[i])) + (amplitudeErrors_[j]*amplitudeErrors_[j]/(amplitudes_[j]*amplitudes_[j])) ); // error due to fluctuations of pedestal (common to both samples) double stat; if(num_>0) stat = num_; // num presampeles used to compute pedestal else stat = 1; // pedestal from db double err2 = amplitudeErrors_[j]*(amplitudes_[i]-amplitudes_[j])/(amplitudes_[j]*amplitudes_[j])/sqrt(stat); //error due to integer round-down. It is relevant to low //amplitudes_ in gainID=1 and negligible otherwise. double err3 = 0.289/amplitudes_[j]; double totalError = sqrt(err1 + err2*err2 +err3*err3); // don't include useless ratios if(totalError < 1.0 && Rtmp>0.001 && Rtmp<exp(double(j-i)/beta)-0.001 ){ Ratio currentRatio = { i, (j-i), Rtmp, totalError }; ratios_.push_back(currentRatio); } } } } // No useful ratios, return zero amplitude and no time measurement if(!ratios_.size() >0) return; // make a vector of Tmax measurements that correspond to each ratio // and based on Alpha-Beta parameterization of the pulse shape for(unsigned int i = 0; i < ratios_.size(); i++){ double stepOverBeta = double(ratios_[i].step)/beta; double offset = double(ratios_[i].index) + alphabeta; double Rmin = ratios_[i].value - ratios_[i].error; if(Rmin<0.001) Rmin=0.001; double Rmax = ratios_[i].value + ratios_[i].error; double RLimit = exp(stepOverBeta)-0.001; if( Rmax > RLimit ) Rmax = RLimit; double time1 = offset - ratios_[i].step/(exp((stepOverBeta-log(Rmin))/alpha)-1.0); double time2 = offset - ratios_[i].step/(exp((stepOverBeta-log(Rmax))/alpha)-1.0); // this is the time measurement based on the ratios[i] double tmax = 0.5 * (time1 + time2); double tmaxerr = 0.5 * sqrt( (time1 - time2)*(time1 - time2) ); // calculate chi2 sumAf = 0; sumff = 0; for(unsigned int it = 0; it < amplitudes_.size(); it++){ double err2 = amplitudeErrors_[it]*amplitudeErrors_[it]; double offset = (double(it) - tmax)/alphabeta; double term1 = 1.0 + offset; if(term1>1e-6){ double f = exp( alpha*(log(1.0+offset) - offset) ); sumAf += amplitudes_[it]*f/err2; sumff += f*f/err2; } } double chi2 = sumAA; double amp = 0; if( sumff > 0 ){ chi2 = sumAA - sumAf*sumAf/sumff; amp = sumAf/sumff; } chi2 /= sum0; // choose reasonable measurements. One might argue what is // reasonable and what is not. if(chi2 > 0 && tmaxerr > 0 && tmax > 0){ Tmax currentTmax={ ratios_[i].index, ratios_[i].step, tmax, tmaxerr, amp, chi2 }; timesAB_.push_back(currentTmax); } } // no reasonable time measurements! if( !(timesAB_.size()> 0)) return; // find minimum chi2 double chi2min = 1.0e+9; //double timeMinimum = 5; //double errorMinimum = 999; for(unsigned int i = 0; i < timesAB_.size(); i++){ if( timesAB_[i].chi2 <= chi2min ){ chi2min = timesAB_[i].chi2; //timeMinimum = timesAB_[i].value; //errorMinimum = timesAB_[i].error; } } // make a weighted average of tmax measurements with "small" chi2 // (within 1 sigma of statistical uncertainty :-) double chi2Limit = chi2min + 1.0; double time_max = 0; double time_wgt = 0; for(unsigned int i = 0; i < timesAB_.size(); i++){ if( timesAB_[i].chi2 < chi2Limit ){ double inverseSigmaSquared = 1.0/(timesAB_[i].error*timesAB_[i].error); time_wgt += inverseSigmaSquared; time_max += timesAB_[i].value*inverseSigmaSquared; } } tMaxAlphaBeta = time_max/time_wgt; tMaxErrorAlphaBeta = 1.0/sqrt(time_wgt); // find amplitude and chi2 sumAf = 0; sumff = 0; for(unsigned int i = 0; i < amplitudes_.size(); i++){ double err2 = amplitudeErrors_[i]*amplitudeErrors_[i]; double offset = (double(i) - tMaxAlphaBeta)/alphabeta; double term1 = 1.0 + offset; if(term1>1e-6){ double f = exp( alpha*(log(1.0+offset) - offset) ); sumAf += amplitudes_[i]*f/err2; sumff += f*f/err2; } } if( sumff > 0 ){ ampMaxAlphaBeta = sumAf/sumff; double chi2AlphaBeta = (sumAA - sumAf*sumAf/sumff)/sum0; if(chi2AlphaBeta > NullChi2){ // null hypothesis is better return; } }else{ // no visible pulse here return; } // if we got to this point, we have a reconstructied Tmax // using RatioAlphaBeta Method. To summarize: // // tMaxAlphaBeta - Tmax value // tMaxErrorAlphaBeta - error on Tmax, but I would not trust it // ampMaxAlphaBeta - amplitude of the pulse // ampMaxError_ - uncertainty of the time sample with max amplitude // // Do Ratio's Method with "large" pulses if( ampMaxAlphaBeta/ampMaxError_ > 5.0 ){ // make a vector of Tmax measurements that correspond to each // ratio. Each measurement have it's value and the error double time_max = 0; double time_wgt = 0; for (unsigned int i = 0; i < ratios_.size(); i++) { if(ratios_[i].step == 1 && ratios_[i].value >= timeFitLimits.first && ratios_[i].value <= timeFitLimits.second ){ double time_max_i = ratios_[i].index; // calculate polynomial for Tmax double u = timeFitParameters[timeFitParameters.size() - 1]; for (int k = timeFitParameters.size() - 2; k >= 0; k--) { u = u * ratios_[i].value + timeFitParameters[k]; } // calculate derivative for Tmax error double du = (timeFitParameters.size() - 1) * timeFitParameters[timeFitParameters.size() - 1]; for (int k = timeFitParameters.size() - 2; k >= 1; k--) { du = du * ratios_[i].value + k * timeFitParameters[k]; } // running sums for weighted average double errorsquared = ratios_[i].error * ratios_[i].error * du * du; if (errorsquared > 0) { time_max += (time_max_i - u) / errorsquared; time_wgt += 1.0 / errorsquared; Tmax currentTmax = { ratios_[i].index, 1, (time_max_i - u), sqrt(errorsquared),0,1 }; times_.push_back(currentTmax); } } } // calculate weighted average of all Tmax measurements if (time_wgt > 0) { tMaxRatio = time_max/time_wgt; tMaxErrorRatio = 1.0/sqrt(time_wgt); // combine RatioAlphaBeta and Ratio Methods if( ampMaxAlphaBeta/ampMaxError_ > 10.0 ){ // use pure Ratio Method calculatedRechit_.timeMax = tMaxRatio; calculatedRechit_.timeError = tMaxErrorRatio; }else{ // combine two methods calculatedRechit_.timeMax = ( tMaxAlphaBeta*(10.0-ampMaxAlphaBeta/ampMaxError_) + tMaxRatio*(ampMaxAlphaBeta/ampMaxError_ - 5.0) )/5.0; calculatedRechit_.timeError = ( tMaxErrorAlphaBeta*(10.0-ampMaxAlphaBeta/ampMaxError_) + tMaxErrorRatio*(ampMaxAlphaBeta/ampMaxError_ - 5.0) )/5.0; } }else{ // use RatioAlphaBeta Method calculatedRechit_.timeMax = tMaxAlphaBeta; calculatedRechit_.timeError = tMaxErrorAlphaBeta; } }else{ // use RatioAlphaBeta Method calculatedRechit_.timeMax = tMaxAlphaBeta; calculatedRechit_.timeError = tMaxErrorAlphaBeta; } }
bool EcalUncalibRecHitRatioMethodAlgo< C >::fixMGPAslew | ( | const C & | dataFrame | ) |
Definition at line 179 of file EcalUncalibRecHitRatioMethodAlgo.h.
References alignCSCRings::e, and query::result.
Referenced by EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerRatio::run().
{ // This fuction finds sample(s) preceeding gain switching and // inflates errors on this sample, therefore, making this sample // invisible for Ratio Method. Only gain switching DOWN is // considered Only gainID=1,2,3 are considered. In case of the // saturation (gainID=0), we keep "distorted" sample because it is // the only chance to make time measurement; the qualilty of it will // be bad anyway. bool result = false; int GainIdPrev; int GainIdNext; for (int iSample = 1; iSample < C::MAXSAMPLES; iSample++) { // only use samples which are desired if (!sampleMask_.useSample(iSample, theDetId_) ) continue; GainIdPrev = dataFrame.sample(iSample-1).gainId(); GainIdNext = dataFrame.sample(iSample).gainId(); if( GainIdPrev>=1 && GainIdPrev<=3 && GainIdNext>=1 && GainIdNext<=3 && GainIdPrev<GainIdNext ){ amplitudes_[iSample-1]=1e-9; amplitudeErrors_[iSample-1]=1e+9; result = true; } } return result; }
CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::getCalculatedRecHit | ( | ) | [inline] |
Definition at line 58 of file EcalUncalibRecHitRatioMethodAlgo.h.
Referenced by EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerRatio::run().
{ return calculatedRechit_; };
void EcalUncalibRecHitRatioMethodAlgo< C >::init | ( | const C & | dataFrame, |
const EcalSampleMask & | sampleMask, | ||
const double * | pedestals, | ||
const double * | pedestalRMSes, | ||
const double * | gainRatios | ||
) |
Definition at line 79 of file EcalUncalibRecHitRatioMethodAlgo.h.
References alignCSCRings::e, and compare_using_db::sample.
Referenced by EcalUncalibRecHitWorkerGlobal::run().
{ sampleMask_ = sampleMask; theDetId_ = DetId(dataFrame.id().rawId()); calculatedRechit_.timeMax = 5; calculatedRechit_.amplitudeMax = 0; calculatedRechit_.timeError = -999; amplitudes_.clear(); amplitudes_.reserve(C::MAXSAMPLES); amplitudeErrors_.clear(); amplitudeErrors_.reserve(C::MAXSAMPLES); ratios_.clear(); ratios_.reserve(C::MAXSAMPLES*(C::MAXSAMPLES-1)/2); times_.clear(); times_.reserve(C::MAXSAMPLES*(C::MAXSAMPLES-1)/2); timesAB_.clear(); timesAB_.reserve(C::MAXSAMPLES*(C::MAXSAMPLES-1)/2); // to obtain gain 12 pedestal: // -> if it's in gain 12, use first sample // --> average it with second sample if in gain 12 and 3-sigma-noise compatible (better LF noise cancellation) // -> else use pedestal from database pedestal_ = 0; num_ = 0; if (dataFrame.sample(0).gainId() == 1 && sampleMask_.useSample(0, theDetId_ ) ) { pedestal_ += double (dataFrame.sample(0).adc()); num_++; } if (num_!=0 && dataFrame.sample(1).gainId() == 1 && sampleMask_.useSample(1, theDetId_) && fabs(dataFrame.sample(1).adc()-dataFrame.sample(0).adc())<3*pedestalRMSes[0]) { pedestal_ += double (dataFrame.sample(1).adc()); num_++; } if (num_ != 0) pedestal_ /= num_; else pedestal_ = pedestals[0]; // fill vector of amplitudes, pedestal subtracted and vector // of amplitude uncertainties Also, find the uncertainty of a // sample with max amplitude. We will use it later. ampMaxError_ = 0; double ampMaxValue = -1000; // ped-subtracted and gain-renormalized samples. It is VERY // IMPORTANT to have samples one clock apart which means to // have vector size equal to MAXSAMPLES double sample; double sampleError; int GainId; for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) { GainId = dataFrame.sample(iSample).gainId(); // only use normally samples which are desired; if sample not to be used // inflate error so won't generate ratio considered for the measurement if (!sampleMask_.useSample(iSample, theDetId_ ) ) { sample = 1e-9; sampleError = 1e+9; } else if (GainId == 1) { sample = double (dataFrame.sample(iSample).adc() - pedestal_); sampleError = pedestalRMSes[0]; } else if (GainId == 2 || GainId == 3){ sample = (double (dataFrame.sample(iSample).adc() - pedestals[GainId - 1])) *gainRatios[GainId - 1]; sampleError = pedestalRMSes[GainId-1]*gainRatios[GainId-1]; } else { sample = 1e-9; // GainId=0 case falls here, from saturation sampleError = 1e+9; // inflate error so won't generate ratio considered for the measurement } if(sampleError>0){ amplitudes_.push_back(sample); amplitudeErrors_.push_back(sampleError); if(ampMaxValue < sample){ ampMaxValue = sample; ampMaxError_ = sampleError; } }else{ // inflate error for useless samples amplitudes_.push_back(sample); amplitudeErrors_.push_back(1e+9); } } }
EcalUncalibratedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::makeRecHit | ( | const C & | dataFrame, |
const EcalSampleMask & | sampleMask, | ||
const double * | pedestals, | ||
const double * | pedestalRMSes, | ||
const double * | gainRatios, | ||
std::vector< double > & | timeFitParameters, | ||
std::vector< double > & | amplitudeFitParameters, | ||
std::pair< double, double > & | timeFitLimits | ||
) | [virtual] |
Definition at line 569 of file EcalUncalibRecHitRatioMethodAlgo.h.
References init.
Referenced by EcalUncalibRecHitWorkerRatio::run().
{ init( dataFrame, sampleMask, pedestals, pedestalRMSes, gainRatios ); computeTime( timeFitParameters, timeFitLimits, amplitudeFitParameters ); computeAmplitude( amplitudeFitParameters ); // 1st parameters is id // // 2nd parameters is amplitude. It is calculated by this method. // // 3rd parameter is pedestal. It is not calculated. This method // relies on input parameters for pedestals and gain ratio. Return // zero. // // 4th parameter is jitter which is a bad choice to call Tmax. It is // calculated by this method (in 25 nsec clock units) // // GF subtract 5 so that jitter==0 corresponds to synchronous hit // // // 5th parameter is chi2. It is possible to calculate chi2 for // Tmax. It is possible to calculate chi2 for Amax. However, these // values are not very useful without TmaxErr and AmaxErr. This // method can return one value for chi2 but there are 4 different // parameters that have useful information about the quality of Amax // ans Tmax. For now we can return TmaxErr. The quality of Tmax and // Amax can be judged from the magnitude of TmaxErr return EcalUncalibratedRecHit(dataFrame.id(), calculatedRechit_.amplitudeMax, pedestal_, calculatedRechit_.timeMax - 5, calculatedRechit_.timeError); }
std::vector< double > EcalUncalibRecHitRatioMethodAlgo< C >::amplitudeErrors_ [protected] |
Definition at line 66 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< double > EcalUncalibRecHitRatioMethodAlgo< C >::amplitudes_ [protected] |
Definition at line 65 of file EcalUncalibRecHitRatioMethodAlgo.h.
double EcalUncalibRecHitRatioMethodAlgo< C >::ampMaxError_ [protected] |
Definition at line 73 of file EcalUncalibRecHitRatioMethodAlgo.h.
CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::calculatedRechit_ [protected] |
Definition at line 75 of file EcalUncalibRecHitRatioMethodAlgo.h.
Referenced by EcalUncalibRecHitRatioMethodAlgo< EBDataFrame >::getCalculatedRecHit().
int EcalUncalibRecHitRatioMethodAlgo< C >::num_ [protected] |
Definition at line 72 of file EcalUncalibRecHitRatioMethodAlgo.h.
double EcalUncalibRecHitRatioMethodAlgo< C >::pedestal_ [protected] |
Definition at line 71 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< Ratio > EcalUncalibRecHitRatioMethodAlgo< C >::ratios_ [protected] |
Definition at line 67 of file EcalUncalibRecHitRatioMethodAlgo.h.
EcalSampleMask EcalUncalibRecHitRatioMethodAlgo< C >::sampleMask_ [protected] |
Definition at line 63 of file EcalUncalibRecHitRatioMethodAlgo.h.
DetId EcalUncalibRecHitRatioMethodAlgo< C >::theDetId_ [protected] |
Definition at line 64 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< Tmax > EcalUncalibRecHitRatioMethodAlgo< C >::times_ [protected] |
Definition at line 68 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< Tmax > EcalUncalibRecHitRatioMethodAlgo< C >::timesAB_ [protected] |
Definition at line 69 of file EcalUncalibRecHitRatioMethodAlgo.h.