* provided "as is" without express or implied warranty. *
**************************************************************************/
-
+//
// Extraction of amplitude and peak position
// FRom CALO raw data using
// least square fit for the
#include "TGraph.h"
#include "TRandom.h"
+#include "AliEMCALRawResponse.h"
-using namespace std;
+using namespace std;
-#define BAD 4 //CRAP PTH
-
ClassImp( AliCaloRawAnalyzerKStandard )
AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard")
-// fkEulerSquared(7.389056098930650227),
-// fTf1(0)
-// fTau(2.35),
-// fFixTau(kTRUE)
{
fAlgo = Algo::kStandard;
- //comment
-
- /*
- for(int i=0; i < ALTROMAXSAMPLES; i++)
- {
- fXaxis[i] = i;
- }
- */
-
- // fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])", 0, 30 );
-
- // InitFormula(fTf1);
-
- /*
- if (fFixTau)
- {
- fTf1->FixParameter(2, fTau);
- }
- else
- {
- fTf1->ReleaseParameter(2); // allow par. to vary
- fTf1->SetParameter(2, fTau);
- }
- */
}
AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
{
- delete fTf1;
+ // delete fTf1;
}
-/*
-void
-AliCaloRawAnalyzerKStandard::InitFormula( TF1* f)
-{
- f = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])", 0, 30 );
-}
-*/
-
AliCaloFitResults
AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo> &bunchlist, const UInt_t altrocfg1, const UInt_t altrocfg2 )
{
-
+ //Evaluation Amplitude and TOF
Float_t pedEstimate = 0;
short maxADC = 0;
Int_t first = 0;
Float_t chi2 = 0;
Int_t ndf = 0;
Bool_t fitDone = kFALSE;
-
-
int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
- maxADC, timeEstimate, pedEstimate, first, last, fAmpCut );
+ maxADC, timeEstimate, pedEstimate, first, last, (int)fAmpCut );
if (ampEstimate >= fAmpCut )
{
amp += (0.5 - gRandom->Rndm());
}
- //Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
- // lowGain = in.IsLowGain();
-
time = time * TIMEBINWITH;
-
- /////////////!!!!!!!!!time -= in.GetL1Phase();
-
time -= fL1Phase;
- // AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
- // AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf);
-
-
return AliCaloFitResults( -99, -99, fAlgo , amp, time,
- time, chi2, ndf, Ret::kDummy );
-
-
- // AliCaloFitSubarray(index, maxrev, first, last));
-
- }
-
-
+ (int)time, chi2, ndf, Ret::kDummy );
+ }
return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
}
-
-/*
-void
-AliCaloRawAnalyzerKStandard::PrintFitResult(const TF1 *f) const
-{
- //comment
- cout << endl;
- cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
- cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
- cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
- cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
- // cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
- cout << endl << endl;
-}
-*/
-
-
-
-//____________________________________________________________________________
void
AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
-{ // Fits the raw signal time distribution
-
- //--------------------------------------------------
- //Do the fit, different fitting algorithms available
- //--------------------------------------------------
-
- // fprintf(fp, "%s:%d:%s\n", __FILE__, __LINE__, __FUNCTION__ );
-
+{
+ // Fits the raw signal time distribution
int nsamples = lastTimeBin - firstTimeBin + 1;
fitDone = kFALSE;
-
- // switch(fFittingAlgorithm)
- // {
- // case Algo::kStandard:
- // {
- if (nsamples < 3) { return; } // nothing much to fit
- //printf("Standard fitter \n");
-
- // Create Graph to hold data we will fit
+ if (nsamples < 3) { return; }
TGraph *gSig = new TGraph( nsamples);
Int_t timebin = firstTimeBin + i;
gSig->SetPoint(i, timebin, GetReversed(timebin));
}
-
- TF1 * signalF = new TF1("signal", RawResponseFunction, 0, TIMEBINS , 5);
+
+ TF1 * signalF = new TF1("signal", AliEMCALRawResponse::RawResponseFunction, 0, TIMEBINS , 5);
+
signalF->SetParameters(10.,5., TAU ,ORDER,0.); //set all defaults once, just to be safe
signalF->SetParNames("amp","t0","tau","N","ped");
- signalF->FixParameter(2,TAU); // tau in units of time bin
- signalF->FixParameter(3,ORDER); // order
- signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here
+ signalF->FixParameter(2,TAU);
+ signalF->FixParameter(3,ORDER);
+ signalF->FixParameter(4, 0);
signalF->SetParameter(1, time);
signalF->SetParameter(0, amp);
- // 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);
chi2 = signalF->GetChisquare();
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;
+ 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;
-
- //printf("Std : Amp %f, time %g\n",amp, time);
delete gSig; // delete TGraph
-
- // break;
- // }//kStandard Fitter
- //----------------------------
-
- /*
- case Algo::kLogFit:
- {
- if (nsamples < 3) { return; } // nothing much to fit
- //printf("LogFit \n");
-
- // Create Graph to hold data we will fit
- TGraph *gSigLog = new TGraph( nsamples);
- for (int i=0; i<nsamples; i++) {
- Int_t timebin = firstTimeBin + i;
- gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) );
- }
-
- TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0, TIMEBINS , 5);
- signalFLog->SetParameters(2.3, 5.,TAU,ORDER,0.); //set all defaults once, just to be safe
- signalFLog->SetParNames("amplog","t0","tau","N","ped");
- signalFLog->FixParameter(2,TAU); // tau in units of time bin
- signalFLog->FixParameter(3, ORDER); // order
- signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here
- signalFLog->SetParameter(1, time);
- if (amp>=1) {
- signalFLog->SetParameter(0, TMath::Log(amp));
- }
-
- gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points
-
- // assign fit results
- 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);
- delete gSigLog;
- break;
- } //kLogFit
- //----------------------------
- //----------------------------
- }//switch fitting algorithms
- */
- return;
-}
-
-
-//__________________________________________________________________
-void
-AliCaloRawAnalyzerKStandard::FitParabola(const TGraph *gSig, Float_t & amp) const
-{
- //BEG YS alternative methods to calculate the amplitude
- Double_t * ymx = gSig->GetX() ;
- Double_t * ymy = gSig->GetY() ;
- const Int_t kN = 3 ;
- Double_t ymMaxX[kN] = {0., 0., 0.} ;
- Double_t ymMaxY[kN] = {0., 0., 0.} ;
- Double_t ymax = 0. ;
- // find the maximum amplitude
- Int_t ymiMax = 0 ;
- for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
- if (ymy[ymi] > ymMaxY[0] ) {
- ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
- ymMaxX[0] = ymx[ymi] ;
- ymiMax = ymi ;
- }
- }
- // find the maximum by fitting a parabola through the max and the two adjacent samples
- if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
- ymMaxY[1] = ymy[ymiMax+1] ;
- ymMaxY[2] = ymy[ymiMax-1] ;
- ymMaxX[1] = ymx[ymiMax+1] ;
- ymMaxX[2] = ymx[ymiMax-1] ;
- if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
- //fit a parabola through the 3 points y= a+bx+x*x*x
- Double_t sy = 0 ;
- Double_t sx = 0 ;
- Double_t sx2 = 0 ;
- Double_t sx3 = 0 ;
- Double_t sx4 = 0 ;
- Double_t sxy = 0 ;
- Double_t sx2y = 0 ;
- for (Int_t i = 0; i < kN ; i++) {
- sy += ymMaxY[i] ;
- sx += ymMaxX[i] ;
- sx2 += ymMaxX[i]*ymMaxX[i] ;
- sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
- sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
- sxy += ymMaxX[i]*ymMaxY[i] ;
- sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
- }
- Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
- Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
- Double_t c = cN / cD ;
- Double_t b = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
- Double_t a = (sy-b*sx-c*sx2)/kN ;
- Double_t xmax = -b/(2*c) ;
- ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
- amp = ymax;
- }
- }
-
- Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
- if (diff > 0.1)
- amp = ymMaxY[0] ;
- //printf("Yves : Amp %f, time %g\n",amp, time);
- //END YS
return;
}
-
-//__________________________________________________________________
-Double_t
-AliCaloRawAnalyzerKStandard::RawResponseFunction(Double_t *x, Double_t *par)
-{
- // Matches version used in 2007 beam test
- //
- // Shape of the electronics raw reponse:
- // It is a semi-gaussian, 2nd order Gamma function of the general form
- //
- // xx = (t - t0 + tau) / tau [xx is just a convenient help variable]
- // F = A * (xx**N * exp( N * ( 1 - xx) ) for xx >= 0
- // F = 0 for xx < 0
- //
- // parameters:
- // A: par[0] // Amplitude = peak value
- // t0: par[1]
- // tau: par[2]
- // N: par[3]
- // ped: par[4]
- //
- Double_t signal = 0.;
- Double_t tau = par[2];
- Double_t n = par[3];
- Double_t ped = par[4];
- Double_t xx = ( x[0] - par[1] + tau ) / tau ;
-
- if (xx <= 0)
- signal = ped ;
- else {
- signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
- }
- return signal ;
-}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff