|
|
|
#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 double *pedestals, const double *pedestalRMSes, const double *gainRatios) |
| virtual EcalUncalibratedRecHit | makeRecHit (const C &dataFrame, 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_ |
| 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 491 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 187 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 156 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++) {
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 57 of file EcalUncalibRecHitRatioMethodAlgo.h.
Referenced by EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerRatio::run().
{ return calculatedRechit_; };
| void EcalUncalibRecHitRatioMethodAlgo< C >::init | ( | const C & | dataFrame, |
| const double * | pedestals, | ||
| const double * | pedestalRMSes, | ||
| const double * | gainRatios | ||
| ) |
Definition at line 75 of file EcalUncalibRecHitRatioMethodAlgo.h.
References alignCSCRings::e.
Referenced by EcalUncalibRecHitWorkerGlobal::run().
{
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) {
pedestal_ += double (dataFrame.sample(0).adc());
num_++;
}
if (num_!=0 &&
dataFrame.sample(1).gainId() == 1 &&
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();
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 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 542 of file EcalUncalibRecHitRatioMethodAlgo.h.
References init.
Referenced by EcalUncalibRecHitWorkerRatio::run().
{
init( dataFrame, 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 62 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< double > EcalUncalibRecHitRatioMethodAlgo< C >::amplitudes_ [protected] |
Definition at line 61 of file EcalUncalibRecHitRatioMethodAlgo.h.
double EcalUncalibRecHitRatioMethodAlgo< C >::ampMaxError_ [protected] |
Definition at line 69 of file EcalUncalibRecHitRatioMethodAlgo.h.
CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::calculatedRechit_ [protected] |
Definition at line 71 of file EcalUncalibRecHitRatioMethodAlgo.h.
Referenced by EcalUncalibRecHitRatioMethodAlgo< EBDataFrame >::getCalculatedRecHit().
int EcalUncalibRecHitRatioMethodAlgo< C >::num_ [protected] |
Definition at line 68 of file EcalUncalibRecHitRatioMethodAlgo.h.
double EcalUncalibRecHitRatioMethodAlgo< C >::pedestal_ [protected] |
Definition at line 67 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< Ratio > EcalUncalibRecHitRatioMethodAlgo< C >::ratios_ [protected] |
Definition at line 63 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< Tmax > EcalUncalibRecHitRatioMethodAlgo< C >::times_ [protected] |
Definition at line 64 of file EcalUncalibRecHitRatioMethodAlgo.h.
std::vector< Tmax > EcalUncalibRecHitRatioMethodAlgo< C >::timesAB_ [protected] |
Definition at line 65 of file EcalUncalibRecHitRatioMethodAlgo.h.
1.7.3