From 947c8e28ae00737bbcf18a4d27247c97bcd06715 Mon Sep 17 00:00:00 2001 From: dsilverm Date: Tue, 11 May 2010 13:11:30 +0000 Subject: [PATCH] chi2, mean and rms checks for input samples for fitting --- EMCAL/AliCaloFitResults.h | 2 +- EMCAL/AliCaloRawAnalyzer.cxx | 126 ++++++++++++++++++++++++++++++++--- EMCAL/AliCaloRawAnalyzer.h | 8 +++ 3 files changed, 124 insertions(+), 12 deletions(-) diff --git a/EMCAL/AliCaloFitResults.h b/EMCAL/AliCaloFitResults.h index 771d15025ea..8e5faaa469a 100644 --- a/EMCAL/AliCaloFitResults.h +++ b/EMCAL/AliCaloFitResults.h @@ -76,7 +76,7 @@ class AliCaloFitResults Float_t GetTof() const { return fTime; }; Float_t GetTime() const { return fTime; }; Int_t GetMaxTimebin() const { return fMaxTimebin; }; - Float_t GetChisSquare() const { return fChi2Sig;}; + Float_t GetChi2() const { return fChi2Sig;}; Int_t GetNdf() const { return fNdfSig; }; AliCaloFitSubarray GetFitSubarray() const { return fFitSubarray; }; diff --git a/EMCAL/AliCaloRawAnalyzer.cxx b/EMCAL/AliCaloRawAnalyzer.cxx index d989f395cf5..6ce22c5bec7 100644 --- a/EMCAL/AliCaloRawAnalyzer.cxx +++ b/EMCAL/AliCaloRawAnalyzer.cxx @@ -101,25 +101,42 @@ AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, con //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump) int tmpfirst = maxindex; int tmplast = maxindex; - Double_t valFirst = fData[maxindex]; - Double_t valLast = fData[maxindex]; - - while( (tmpfirst >= 0) && (fData[tmpfirst] >= fFitArrayCut) && - (fData[tmpfirst]= 0) && (fData[tmpfirst] >= fFitArrayCut) && (!firstJump) ) { - valFirst = fData[tmpfirst]; + // jump check: + if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex + if (fData[tmpfirst] >= prevFirst) { + firstJump = true; + } + } + prevFirst = fData[tmpfirst]; tmpfirst -- ; } - while( (tmplast < length) && (fData[tmplast] >= fFitArrayCut) && - (fData[tmplast]= fFitArrayCut) && (!lastJump) ) { - valLast = fData[tmplast]; + // jump check: + if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex + if (fData[tmplast] >= prevLast) { + lastJump = true; + } + } + prevLast = fData[tmplast]; tmplast ++; } - *first = tmpfirst +1; - *last = tmplast -1; + // we keep one pre- or post- sample if we can (as in online) + // check though if we ended up on a 'jump', or out of bounds: if so, back up + if (firstJump || tmpfirst<0) tmpfirst ++; + if (lastJump || tmplast>=length) tmplast --; + + *first = tmpfirst; + *last = tmplast; return; } @@ -296,6 +313,93 @@ AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const } +Double_t +AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time, + const Int_t first, const Int_t last, + const Double_t adcErr, + const Double_t tau) +{ + // Input: + // amp - max amplitude; + // time - time of max amplitude; + // first, last - sample array indices to be used + // adcErr - nominal error of amplitude measurement (one value for all channels) + // if adcErr<0 that mean adcErr=1. + // tau - filter time response (in timebin units) + // Output: + // chi2 - chi2 + + if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. + // or, first is negative, the indices are not valid + return AliCaloFitResults::kDummy; + } + + int nsamples = last - first + 1; + // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i \n", first, last, nsamples); + + Int_t x = 0; + Double_t chi2 = 0; + Double_t dy = 0.0, xx = 0.0, f=0.0; + + for (Int_t i=0; i 0) { + f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ; + } + dy = fReversed[x] - f; + chi2 += dy*dy; + // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); + } + + if (adcErr>0.0) { // weight chi2 + chi2 /= (adcErr*adcErr); + } + return chi2; +} + + +void +AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last, + Double_t & mean, Double_t & rms) +{ + // Input: + // first, last - sample array indices to be used + // Output: + // mean and RMS of samples + // + // To possibly be used to differentiate good signals from bad before fitting + // + mean = AliCaloFitResults::kDummy; + rms = AliCaloFitResults::kDummy; + + if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. + // or, first is negative, the indices are not valid + return; + } + + int nsamples = last - first + 1; + // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples); + + int x = 0; + Double_t sampleSum = 0; // sum of samples + Double_t squaredSampleSum = 0; // sum of samples squared + + for (Int_t i=0; i &/*bunchvector*/, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/) { // method to do the selection of what should possibly be fitted diff --git a/EMCAL/AliCaloRawAnalyzer.h b/EMCAL/AliCaloRawAnalyzer.h index 2994bb66cf6..a9df4e6fbb3 100644 --- a/EMCAL/AliCaloRawAnalyzer.h +++ b/EMCAL/AliCaloRawAnalyzer.h @@ -67,6 +67,14 @@ class AliCaloRawAnalyzer : public TObject const char * GetAlgoName() const { return fName; }; const char * GetAlgoAbbr() const { return fNameShort; }; + Double_t CalculateChi2(const Double_t amp, const Double_t time, + const Int_t first, const Int_t last, + const Double_t adcErr=1, + const Double_t tau=2.35); + + void CalculateMeanAndRMS(const Int_t first, const Int_t last, + Double_t & mean, Double_t & rms); + protected: short Max( const AliCaloBunchInfo *const bunch, int *const maxindex) const; UShort_t Max(const UShort_t *data, const int length ) const; -- 2.43.0