#include "AliCaloFitResults.h"
#include "AliLog.h"
#include "TMath.h"
+#include <stdexcept>
#include <iostream>
#include "TF1.h"
#include "TGraph.h"
#define BAD 4 //CRAP PTH
+
+ClassImp( AliCaloRawAnalyzerLMS )
-AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer(),
+AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit", "LMS"),
fkEulerSquared(7.389056098930650227),
- fSig(0),
- fTf1(0)
+ fTf1(0),
+ fTau(2.35),
+ fFixTau(kTRUE)
{
+
+ fAlgo = Algo::kLMS;
//comment
for(int i=0; i < MAXSAMPLES; i++)
{
fXaxis[i] = i;
}
- fTf1 = new TF1( "myformula", "[0]*((x - [1])/2.35)^2*exp(-2*(x -[1])/2.35)", 0, 30 );
- // fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^[3]*exp(-[3]*(x -[1])/[2])", 0, 30 );
-
+ fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])", 0, 30 );
+ if (fFixTau)
+ {
+ fTf1->FixParameter(2, fTau);
+ }
+ else
+ {
+ fTf1->ReleaseParameter(2); // allow par. to vary
+ fTf1->SetParameter(2, fTau);
+ }
}
if( index >= 0)
{
Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
- int first;
- int last;
Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
-
- if ( maxf > fAmpCut )
+ 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 || ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
+ {
+ return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
+ }
+ else if ( maxf >= fAmpCut )
{
- SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex - bunchvector.at(index).GetStartBin(), &first, &last);
- int nsamples = last - first;
+ int first = 0;
+ int last = 0;
+ SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
+ int nsamples = last - first + 1;
- if( ( nsamples ) > fNsampleCut )
+ if( ( nsamples ) >= fNsampleCut )
{
-
- TGraph *graph = new TGraph( last - first, fXaxis, &fReversed[first] );
+ 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(2, 2);
- // fTf1->SetParameter(2, 3);
-
- // return AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1);
-
- Short_t tmpStatus = graph->Fit(fTf1, "Q0RW");
-
- // return AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1);
-
+ 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);
+ }
+ else {
+ fTf1->ReleaseParameter(2); // allow par. to vary
+ fTf1->SetParameter(2, fTau);
+ }
+
+ 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, Ret::kNoFit, maxf, timebinOffset,
+ timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
+ }
+
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,
- fTf1->GetParameter(0)/fkEulerSquared,
- fTf1->GetParameter(1) + maxampindex,
- fTf1->GetChisquare(),
- fTf1->GetNDF());
-
-
+ return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
+ fTf1->GetParameter(0)/fkEulerSquared,
+ tmax,
+ timebinOffset,
+ fTf1->GetChisquare(),
+ fTf1->GetNDF(),
+ Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
+
// delete graph;
}
else
{
- return AliCaloFitResults( maxamp, ped, -1, maxf, -1, -1, -1 );
+ Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
+ Int_t ndf = last - first - 1; // nsamples - 2
+ return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
+ timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
}
- }
- else
- {
- return AliCaloFitResults( maxamp , ped, -1, maxf, -1, -1, -1 );
- }
+ } // ampcut
}
- else
- {
- return AliCaloFitResults( -99, -99 );
- }
-
- return AliCaloFitResults( -99, -99 );
+ return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
}
// cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
cout << endl << endl;
}
+