X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliCaloRawAnalyzerKStandard.cxx;h=0abcd7b3f3403317c8312af37fb1dfe12b274bcd;hb=32bff2dafd88f9fe4929c25cc49a2b1ef2acb4ef;hp=fd5bd0d09ff095058468e5da5093074375ebdd16;hpb=92d9f31756fbed289bb86e1ca378c8db95224dd0;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliCaloRawAnalyzerKStandard.cxx b/EMCAL/AliCaloRawAnalyzerKStandard.cxx index fd5bd0d09ff..0abcd7b3f34 100644 --- a/EMCAL/AliCaloRawAnalyzerKStandard.cxx +++ b/EMCAL/AliCaloRawAnalyzerKStandard.cxx @@ -16,7 +16,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ - +// // Extraction of amplitude and peak position // FRom CALO raw data using // least square fit for the @@ -35,53 +35,31 @@ #include "TGraph.h" #include "TRandom.h" +#include "AliEMCALRawResponse.h" -using namespace std; +using namespace std; -#define BAD 4 //CRAP PTH - ClassImp( AliCaloRawAnalyzerKStandard ) -AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzer("Chi Square ( kStandard )", "KStandard"), - fkEulerSquared(7.389056098930650227), - fTf1(0), - fTau(2.35), - fFixTau(kTRUE) +AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard") { 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 ); - if (fFixTau) - { - fTf1->FixParameter(2, fTau); - } - else - { - fTf1->ReleaseParameter(2); // allow par. to vary - fTf1->SetParameter(2, fTau); - } } AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard() { - delete fTf1; + // delete fTf1; } - AliCaloFitResults -AliCaloRawAnalyzerKStandard::Evaluate( const vector &bunchlist, const UInt_t altrocfg1, const UInt_t altrocfg2 ) +AliCaloRawAnalyzerKStandard::Evaluate( const vector &bunchlist, UInt_t altrocfg1, UInt_t altrocfg2 ) { - + //Evaluation Amplitude and TOF Float_t pedEstimate = 0; short maxADC = 0; Int_t first = 0; @@ -94,10 +72,8 @@ AliCaloRawAnalyzerKStandard::Evaluate( const vector &bunchlis 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 ) @@ -132,173 +108,24 @@ AliCaloRawAnalyzerKStandard::Evaluate( const vector &bunchlis { 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)); - - } - - - return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid ); -} - - - - -/* - return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset, - timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); - } - } // ampcut - } - return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid ); -*/ - - -/* - // Extracting signal parameters using fitting - short maxampindex; //index of maximum amplitude - short maxamp; //Maximum amplitude - int index = SelectBunch( bunchvector, &maxampindex, &maxamp ); - - if( index >= 0) - { - Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed ); - 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 || ( 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 ) - { - int first = 0; - int last = 0; - SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut); - int nsamples = last - first + 1; - - 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, 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 , Ret::kFitPar, - fTf1->GetParameter(0)/fkEulerSquared, - tmax, - timebinOffset, - fTf1->GetChisquare(), - fTf1->GetNDF(), - Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); - - // delete graph; - - } - else - { - - 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) ); - } - } // ampcut - } + (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__ ); - + AliCaloRawAnalyzerKStandard::FitRaw( Int_t firstTimeBin, Int_t lastTimeBin, + Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const +{ + // 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); @@ -307,182 +134,36 @@ void 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; iSetPoint(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 ; -} -