From 507751cef0c49a3df98b4e4eaf48e5812c4b8c5b Mon Sep 17 00:00:00 2001 From: dsilverm Date: Fri, 19 Mar 2010 00:34:32 +0000 Subject: [PATCH] more standardizing of result return codes + always return tmax info + loose parameter limits and consistency checks + ADC smearing for no fits + try-catch for possible minuit exceptions --- EMCAL/AliCaloFitResults.cxx | 42 ++++++++---- EMCAL/AliCaloFitResults.h | 29 ++++++-- EMCAL/AliCaloRawAnalyzer.cxx | 40 +++++++++-- EMCAL/AliCaloRawAnalyzer.h | 1 + EMCAL/AliCaloRawAnalyzerCrude.cxx | 6 +- EMCAL/AliCaloRawAnalyzerFastFit.cxx | 14 ++-- EMCAL/AliCaloRawAnalyzerLMS.cxx | 34 +++++++--- EMCAL/AliCaloRawAnalyzerNN.cxx | 10 ++- EMCAL/AliCaloRawAnalyzerPeakFinder.cxx | 9 ++- EMCAL/AliEMCALRawUtils.cxx | 92 +++++++++++++++----------- EMCAL/AliEMCALRawUtils.h | 2 +- 11 files changed, 187 insertions(+), 92 deletions(-) diff --git a/EMCAL/AliCaloFitResults.cxx b/EMCAL/AliCaloFitResults.cxx index 2d0348b90ac..37b149e9f6b 100644 --- a/EMCAL/AliCaloFitResults.cxx +++ b/EMCAL/AliCaloFitResults.cxx @@ -27,16 +27,18 @@ // Applied. fStatus holds information on wether or not // The signal was fitted sucessfully. fStatus might have a different meaning If other // procedures than A different meaning Fitting is applied + AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, const Short_t fitstatus, const Float_t amp, - const Float_t t0, const Float_t chi, + const Float_t time, const Int_t maxTimebin, const Float_t chi, const Int_t ndf, Int_t minSig, const AliCaloFitSubarray fitSubarray) : fMaxSig(maxSig), fPed(ped), fStatus(fitstatus), fAmpSig(amp), - fT0(t0), + fTime(time), + fMaxTimebin(maxTimebin), fChi2Sig(chi), fNdfSig(ndf), fMinSig(minSig), @@ -47,34 +49,50 @@ AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, const Short_t fitstatus, const Float_t amp, - const Float_t t0, const Float_t chi, + const Float_t time, const Int_t maxTimebin, const Float_t chi, const Int_t ndf, Int_t minSig ) : fMaxSig(maxSig), fPed(ped), fStatus(fitstatus), fAmpSig(amp), - fT0(t0), + fTime(time), + fMaxTimebin(maxTimebin), fChi2Sig(chi), fNdfSig(ndf), fMinSig(minSig), fFitSubarray(kDummy) { - } +AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, + const Short_t fitstatus, const Float_t amp, + const Int_t maxTimebin) : + fMaxSig(maxSig), + fPed(ped), + fStatus(fitstatus), + fAmpSig(amp), + fTime(maxTimebin), + fMaxTimebin(maxTimebin), + fChi2Sig( kNoFit ), + fNdfSig( kNoFit ), + fMinSig( kNoFit ), + fFitSubarray( kNoFit ) +{ +} AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Int_t minSig) : fMaxSig(maxSig), - fPed(kNoFit), - fStatus( kNoFit ), - fAmpSig( kNoFit ), - fT0(kNoFit), - fChi2Sig( kNoFit ), - fNdfSig( kNoFit), + fPed( kInvalid ), + fStatus( kInvalid ), + fAmpSig( kInvalid ), + fTime( kInvalid ), + fMaxTimebin( kInvalid ), + fChi2Sig( kInvalid ), + fNdfSig( kInvalid), fMinSig (minSig), - fFitSubarray(kNoFit) + fFitSubarray(kInvalid) { } diff --git a/EMCAL/AliCaloFitResults.h b/EMCAL/AliCaloFitResults.h index bd141bff409..3a6cfced71a 100644 --- a/EMCAL/AliCaloFitResults.h +++ b/EMCAL/AliCaloFitResults.h @@ -30,13 +30,19 @@ class AliCaloFitResults { public: - enum kReturnCode {kDummy=-1, kCrude=-9, kNoFit=-99, kInvalid=-9999};// possible return values + enum kReturnCode {kFitPar=1, kDummy=-1, kCrude=-9, kNoFit=-99, kInvalid=-9999};// possible return values + // kFitPar: method fit or parametrization was used + // kDummy: just a filler parameter, if e.g. chi2 not available + // kCrude: maximum was used + // kNoFit: maximum was used, exception handling for fit invoked + // kInvalid: could not even look for maximum explicit AliCaloFitResults( const Int_t maxSig, const Float_t ped, const Short_t fitStatus, const Float_t amp, - const Float_t t0, + const Float_t time, + const Int_t maxTimebin, const Float_t chi, const Int_t ndf, const Int_t minSig, @@ -46,13 +52,21 @@ class AliCaloFitResults const Float_t ped, const Short_t fitStatus, const Float_t amp, - const Float_t t0, + const Float_t time, + const Int_t maxTimebin, const Float_t chi, const Int_t ndf, const Int_t minSig = kDummy); + // shorter interface when no fit is done + explicit AliCaloFitResults( const Int_t maxSig, + const Float_t ped, + const Short_t fitStatus, + const Float_t amp, + const Int_t maxTimebin); + + // minimum interface explicit AliCaloFitResults( const Int_t maxSig, const Int_t minSig ); - //AliCaloFitResults( const Int_t maxSig, const Int_t minSig ); virtual ~AliCaloFitResults(); @@ -61,7 +75,9 @@ class AliCaloFitResults Int_t GetMinSig() const { return fMinSig;}; Int_t GetStatus() const { return fStatus;}; Float_t GetAmp() const { return fAmpSig; }; - Float_t GetTof() const { return fT0; }; + Float_t GetTof() const { return fTime; }; + Float_t GetTime() const { return fTime; }; + Int_t GetMaxTimebin() const { return fMaxTimebin; }; Float_t GetChisSquare() const { return fChi2Sig;}; Int_t GetNdf() const { return fNdfSig; }; AliCaloFitSubarray GetFitSubarray() const { return fFitSubarray; }; @@ -72,7 +88,8 @@ class AliCaloFitResults Float_t fPed; //Pedestal Int_t fStatus; //Sucess or failure of fitting pocedure Float_t fAmpSig; //Amplitude in entities of ADC counts - Float_t fT0; //Start time of signal in entities of sample intervals + Float_t fTime; //peak/max time of signal in entities of sample intervals + Int_t fMaxTimebin; // timebin with maximum ADC value Float_t fChi2Sig; //Chi Square of fit Int_t fNdfSig; //Number of degrees of freedom of fit Int_t fMinSig; //Pedestal diff --git a/EMCAL/AliCaloRawAnalyzer.cxx b/EMCAL/AliCaloRawAnalyzer.cxx index 4dbd2c6e190..38b317b7d25 100644 --- a/EMCAL/AliCaloRawAnalyzer.cxx +++ b/EMCAL/AliCaloRawAnalyzer.cxx @@ -192,6 +192,33 @@ AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxind } +bool +AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const +{ + // a bunch is considered invalid if the maximum is in the first or last time-bin + short tmpmax = -1; + int tmpindex = -1; + const UShort_t *sig = bunch->GetData(); + + for(int i=0; i < bunch->GetLength(); i++ ) + { + if( sig[i] > tmpmax ) + { + tmpmax = sig[i]; + tmpindex = i; + } + } + + bool bunchOK = true; + if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) ) + { + bunchOK = false; + } + + return bunchOK; +} + + int AliCaloRawAnalyzer::SelectBunch( const vector &bunchvector,short *const maxampbin, short *const maxamplitude ) const { @@ -216,10 +243,13 @@ AliCaloRawAnalyzer::SelectBunch( const vector &bunchvector,sho } } - - // *maxampbin = indx; - // *maxamplitude = max; - + if (bunchindex >= 0) { + bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) ); + if (! bunchOK) { + bunchindex = -1; + } + } + return bunchindex; } @@ -272,7 +302,7 @@ AliCaloFitResults AliCaloRawAnalyzer::Evaluate( const vector &/*bunchvector*/, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/) { // method to do the selection of what should possibly be fitted // not implemented for base class - return AliCaloFitResults( 0, 0, 0, 0, 0, 0, 0 ); + return AliCaloFitResults( 0, 0 ); } diff --git a/EMCAL/AliCaloRawAnalyzer.h b/EMCAL/AliCaloRawAnalyzer.h index 2575d3b159b..7fab35e66de 100644 --- a/EMCAL/AliCaloRawAnalyzer.h +++ b/EMCAL/AliCaloRawAnalyzer.h @@ -82,6 +82,7 @@ class AliCaloRawAnalyzer : public TObject protected: short Max( const AliCaloBunchInfo *const bunch, int *const maxindex) const; UShort_t Max(const UShort_t *data, const int length ) const; + bool CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch) const; bool IsInTimeRange( const int maxindex ) const; Float_t ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t altrocfg1, const UInt_t altrocfg2, double *outarray ) const; int SelectBunch( const vector &bunchvector, short *const maxampbin, short *const maxamplitude ) const; diff --git a/EMCAL/AliCaloRawAnalyzerCrude.cxx b/EMCAL/AliCaloRawAnalyzerCrude.cxx index 2cf79c86148..191a60d71c0 100644 --- a/EMCAL/AliCaloRawAnalyzerCrude.cxx +++ b/EMCAL/AliCaloRawAnalyzerCrude.cxx @@ -48,11 +48,11 @@ AliCaloRawAnalyzerCrude::Evaluate(const vector &bunchvector, c // Evaluation of signal parameters if( bunchvector.size() <= 0 ) { - return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid , AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid ); + return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid); } Int_t amp = 0; - Float_t tof = -99; + Int_t tof = -99; const UShort_t *sig; double ped = EvaluatePedestal( bunchvector.at(0).GetData(), bunchvector.at(0).GetLength() ) ; @@ -72,7 +72,7 @@ AliCaloRawAnalyzerCrude::Evaluate(const vector &bunchvector, c //:EvaluatePedestal(const UShort_t * const data, const int length ) // double ped = EvaluatePedestal(sig, length) ; - return AliCaloFitResults(amp, ped, AliCaloFitResults::kNoFit, amp - ped, tof, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit ); + return AliCaloFitResults(amp, ped, AliCaloFitResults::kCrude, amp - ped, tof); } //end Crude diff --git a/EMCAL/AliCaloRawAnalyzerFastFit.cxx b/EMCAL/AliCaloRawAnalyzerFastFit.cxx index ae3a20dbda0..48e2b5e0f1e 100644 --- a/EMCAL/AliCaloRawAnalyzerFastFit.cxx +++ b/EMCAL/AliCaloRawAnalyzerFastFit.cxx @@ -69,7 +69,6 @@ AliCaloRawAnalyzerFastFit::Evaluate( const vector &bunchvector Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed ); int maxrev = maxampindex - bunchvector.at(index).GetStartBin(); short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1); - if ( maxf > fAmpCut ) { SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev , &first, &last); @@ -87,26 +86,27 @@ AliCaloRawAnalyzerFastFit::Evaluate( const vector &bunchvector Double_t eSignal = 1; // nominal 1 ADC error Double_t dAmp = maxf; Double_t eAmp = 0; - Double_t dTime = 0; + Double_t dTime0 = 0; Double_t eTime = 0; Double_t chi2 = 0; Double_t dTau = 2.35; // time-bin units AliCaloFastAltroFitv0::FastFit(fXAxis, ordered , nsamples, - eSignal, dTau, dAmp, eAmp, dTime, eTime, chi2); + eSignal, dTau, dAmp, eAmp, dTime0, eTime, chi2); - return AliCaloFitResults(maxamp, ped, AliCaloFitResults::kDummy, dAmp, dTime+timebinOffset, chi2, AliCaloFitResults::kDummy, + Double_t dTimeMax = dTime0 + timebinOffset - (maxrev - first) // abs. t0 + + dTau; // +tau, makes sum tmax + return AliCaloFitResults(maxamp, ped, AliCaloFitResults::kFitPar, dAmp, dTimeMax, timebinOffset, chi2, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); } // samplecut else { - return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit, - AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); + return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } } // ampcut else { - return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit); + return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } } // bunch index diff --git a/EMCAL/AliCaloRawAnalyzerLMS.cxx b/EMCAL/AliCaloRawAnalyzerLMS.cxx index 0be5fd0a2fa..ddeadb9eae8 100644 --- a/EMCAL/AliCaloRawAnalyzerLMS.cxx +++ b/EMCAL/AliCaloRawAnalyzerLMS.cxx @@ -29,6 +29,7 @@ #include "AliCaloFitResults.h" #include "AliLog.h" #include "TMath.h" +#include #include #include "TF1.h" #include "TGraph.h" @@ -87,8 +88,8 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector &bunchvector, c int last; Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed ); short maxrev = maxampindex - bunchvector.at(index).GetStartBin(); + // timebinOffset is timebin value at maximum (maxrev) short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1); - if ( maxf > fAmpCut ) { SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last); @@ -96,10 +97,13 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector &bunchvector, c if( ( nsamples ) >= fNsampleCut ) { - + Float_t tmax = (maxrev - first); // local tmax estimate TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] ); fTf1->SetParameter(0, maxf*fkEulerSquared ); - fTf1->SetParameter(1, 0.2); + fTf1->SetParameter(1, tmax - fTau); + // set rather loose parameter limits + fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared ); + fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4); if (fFixTau) { fTf1->FixParameter(2, fTau); @@ -109,18 +113,29 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector &bunchvector, c fTf1->SetParameter(2, fTau); } - Short_t tmpStatus = graph->Fit(fTf1, "Q0RW"); - + Short_t tmpStatus = 0; + try { + tmpStatus = graph->Fit(fTf1, "Q0RW"); + } + catch (const std::exception & e) { + AliError( Form("TGraph Fit exception %s", e.what()) ); + return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, timebinOffset); + } + if( fVerbose == true ) { AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) ); PrintFitResult( fTf1 ) ; } + // global tmax + tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0 + + fTf1->GetParameter(2); // +tau, makes sum tmax delete graph; - return AliCaloFitResults( maxamp, ped , tmpStatus, + return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, fTf1->GetParameter(0)/fkEulerSquared, - fTf1->GetParameter(1) + timebinOffset, + tmax, + timebinOffset, fTf1->GetChisquare(), fTf1->GetNDF(), AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); @@ -130,13 +145,12 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector &bunchvector, c } else { - return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit, - AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); + return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } } else { - return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit); + return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } } diff --git a/EMCAL/AliCaloRawAnalyzerNN.cxx b/EMCAL/AliCaloRawAnalyzerNN.cxx index ac4565ae2e7..5d5f9fba0d4 100644 --- a/EMCAL/AliCaloRawAnalyzerNN.cxx +++ b/EMCAL/AliCaloRawAnalyzerNN.cxx @@ -62,7 +62,7 @@ AliCaloRawAnalyzerNN::Evaluate( const vector &bunchvector, // The eveluation of Peak position and amplitude using the Neural Network if( bunchvector.size() <= 0 ) { - return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid , AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid ); + return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid); } short maxampindex; @@ -90,8 +90,7 @@ AliCaloRawAnalyzerNN::Evaluate( const vector &bunchvector, { if ( ( maxrev - first) < 2 && (last - maxrev ) < 2) { - return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit, - AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); + return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } else { @@ -105,13 +104,12 @@ AliCaloRawAnalyzerNN::Evaluate( const vector &bunchvector, double amp = (maxamp - ped)*fNeuralNet->Value( 0, fNNInput[0], fNNInput[1], fNNInput[2], fNNInput[3], fNNInput[4]); double tof = (fNeuralNet->Value( 1, fNNInput[0], fNNInput[1], fNNInput[2], fNNInput[3], fNNInput[4]) + timebinOffset ) ; - return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kDummy, amp , tof, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, + return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, amp , tof, timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); } } - return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit, - AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); + return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } diff --git a/EMCAL/AliCaloRawAnalyzerPeakFinder.cxx b/EMCAL/AliCaloRawAnalyzerPeakFinder.cxx index 78910119f59..be45052dcf1 100644 --- a/EMCAL/AliCaloRawAnalyzerPeakFinder.cxx +++ b/EMCAL/AliCaloRawAnalyzerPeakFinder.cxx @@ -133,7 +133,7 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector &bunchvec if( maxf < fAmpCut || ( maxamp - ped) > 900 ) // (maxamp - ped) > 900 = Close to saturation (use low gain then) { // cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex << endl; - return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit ); + return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } int first; @@ -221,18 +221,17 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector &bunchvec // tof = tof/fAmpA[tmpindex]; - return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kDummy, fAmpA[tmpindex], tof, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, + return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, fAmpA[tmpindex], tof, timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); } else { - return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit, - AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); + return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kCrude, maxf, timebinOffset); } } else { - return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit); + return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset); } } // cout << __FILE__ << __LINE__ << "WARNING, returning amp = -1 " << endl; diff --git a/EMCAL/AliEMCALRawUtils.cxx b/EMCAL/AliEMCALRawUtils.cxx index 09b92c752f1..4367eb79cc3 100644 --- a/EMCAL/AliEMCALRawUtils.cxx +++ b/EMCAL/AliEMCALRawUtils.cxx @@ -27,6 +27,7 @@ //*-- Author: Marco van Leeuwen (LBL) #include "AliEMCALRawUtils.h" +#include #include "TF1.h" #include "TGraph.h" @@ -376,27 +377,33 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, } // loop over bunches - if ( caloFlag < 2 ){ // ALTRO + if ( caloFlag < 2 ){ // ALTRO - Float_t time = 0; - Float_t amp = 0; + Float_t time = 0; + Float_t amp = 0; + short timeEstimate = 0; + Float_t ampEstimate = 0; + Bool_t fitDone = kFALSE; if ( fFittingAlgorithm == kFastFit || fFittingAlgorithm == kNeuralNet || fFittingAlgorithm == kLMS || fFittingAlgorithm == kPeakFinder || fFittingAlgorithm == kCrude) { // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2()); amp = fitResults.GetAmp(); - time = fitResults.GetTof(); + time = fitResults.GetTime(); + timeEstimate = fitResults.GetMaxTimebin(); + ampEstimate = fitResults.GetMaxSig(); + if (fitResults.GetStatus() == AliCaloFitResults::kFitPar) { + fitDone = kTRUE; + } } else { // for the other methods we for now use the functionality of // AliCaloRawAnalyzer as well, to select samples and prepare for fits, // if it looks like there is something to fit // parameters init. - Float_t ampEstimate = 0; + Float_t pedEstimate = 0; short maxADC = 0; - short timeEstimate = 0; - Float_t pedEstimate = 0; Int_t first = 0; Int_t last = 0; Int_t bunchIndex = 0; @@ -412,35 +419,35 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr, time = timeEstimate; // maxrev in AliCaloRawAnalyzer speak; comes with an offset w.r.t. real timebin Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1); - amp = ampEstimate; if ( nsamples > 1 ) { // possibly something to fit - FitRaw(first, last, amp, time); + FitRaw(first, last, amp, time, fitDone); time += timebinOffset; + timeEstimate += timebinOffset; } - if ( amp>0 && time>0 ) { // brief sanity check of fit results - - // check fit results: should be consistent with initial estimates - // more magic numbers, but very loose cuts, for now.. - // We have checked that amp and ampEstimate values are positive so division for assymmetry - // calculation should be OK/safe - Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); - if ( (TMath::Abs(ampAsymm) > 0.1) ) { - AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations ped %f max-ped %f time %d", - amp, time, pedEstimate, ampEstimate, timeEstimate)); - - // what should do we do then? skip this channel or assign the simple estimate? - // for now just overwrite the fit results with the simple estimate - amp = ampEstimate; - time = timeEstimate; - } // asymm check - } // amp & time check } // ampEstimate check } // method selection + + if ( fitDone ) { // brief sanity check of fit results + Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate); + Float_t timeDiff = time - timeEstimate; + if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) { + // AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations amp %f time %d", amp, time, ampEstimate, timeEstimate)); + + // for now just overwrite the fit results with the simple/initial estimate + amp = ampEstimate; + time = timeEstimate; + fitDone = kFALSE; + } + } // fitDone if (amp > fNoiseThreshold && ampRndm()); // Rndm generates a number in ]0,1] + } + Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; lowGain = in.IsLowGain(); @@ -544,13 +551,14 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain } //____________________________________________________________________________ -void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time) const +void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const { // Fits the raw signal time distribution //-------------------------------------------------- //Do the fit, different fitting algorithms available //-------------------------------------------------- int nsamples = lastTimeBin - firstTimeBin + 1; + fitDone = kFALSE; switch(fFittingAlgorithm) { case kStandard: @@ -573,18 +581,27 @@ void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here signalF->SetParameter(1, time); signalF->SetParameter(0, amp); - - gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points - - // assign fit results - amp = signalF->GetParameter(0); - time = signalF->GetParameter(1); - + // set rather loose parameter limits + signalF->SetParLimits(0, 0.5*amp, 2*amp ); + signalF->SetParLimits(1, time - 4, time + 4); + + try { + gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points + // assign fit results + amp = signalF->GetParameter(0); + time = signalF->GetParameter(1); + + // cross-check with ParabolaFit to see if the results make sense + FitParabola(gSig, amp); // amp is possibly updated + fitDone = kTRUE; + } + catch (const std::exception & e) { + AliError( Form("TGraph Fit exception %s", e.what()) ); + // stay with default amp and time in case of exception, i.e. no special action required + fitDone = kFALSE; + } delete signalF; - // cross-check with ParabolaFit to see if the results make sense - FitParabola(gSig, amp); // amp is possibly updated - //printf("Std : Amp %f, time %g\n",amp, time); delete gSig; // delete TGraph @@ -620,6 +637,7 @@ void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp amp = TMath::Exp(amplog); time = signalFLog->GetParameter(1); + fitDone = kTRUE; delete signalFLog; //printf("LogFit: Amp %f, time %g\n",amp, time); diff --git a/EMCAL/AliEMCALRawUtils.h b/EMCAL/AliEMCALRawUtils.h index 6d456764bc6..371c56c5a28 100644 --- a/EMCAL/AliEMCALRawUtils.h +++ b/EMCAL/AliEMCALRawUtils.h @@ -88,7 +88,7 @@ class AliEMCALRawUtils : public TObject { // Signal shape functions - void FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time) const ; + void FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const ; void FitParabola(const TGraph *gSig, Float_t & amp) const ; static Double_t RawResponseFunction(Double_t *x, Double_t *par); static Double_t RawResponseFunctionLog(Double_t *x, Double_t *par); -- 2.43.0