From 6a265faf9bc1ec37f0b8b255250850f15c69c48e Mon Sep 17 00:00:00 2001 From: prsnko Date: Wed, 11 Aug 2010 15:26:40 +0000 Subject: [PATCH] Improved efficiency; use same signal shape for all channels --- PHOS/AliPHOSRawFitterv2.cxx | 479 +++++++++++++++++++++++++----------- PHOS/AliPHOSRawFitterv2.h | 24 +- 2 files changed, 351 insertions(+), 152 deletions(-) diff --git a/PHOS/AliPHOSRawFitterv2.cxx b/PHOS/AliPHOSRawFitterv2.cxx index 1a1a6a4f545..04a8b40995c 100644 --- a/PHOS/AliPHOSRawFitterv2.cxx +++ b/PHOS/AliPHOSRawFitterv2.cxx @@ -15,7 +15,11 @@ /* $Id: $ */ -// This class plots samples and qualify their quality according to their shape +// This class extracts the signal parameters (energy, time, quality) +// from ALTRO samples. Energy is in ADC counts, time is in time bin units. +// Class uses FastFitting algorithm to fit sample and extract time and Amplitude +// and evaluate sample quality = (chi^2/NDF)/some parameterization providing +// efficiency close to 100% // // Typical use case: // AliPHOSRawFitter *fitter=new AliPHOSRawFitter(); @@ -26,17 +30,16 @@ // Double_t time = fitter.GetTime(); // Bool_t isLowGain = fitter.GetCaloFlag()==0; -// Author: Dmitri Peressounko (Oct.2008) -// Modified: Yuri Kharlov (Jul.2009) +// Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx) // --- ROOT system --- +#include "TArrayI.h" #include "TList.h" #include "TMath.h" -#include "TMinuit.h" -#include "TCanvas.h" -#include "TH1.h" -#include "TH2.h" +#include "TH1I.h" #include "TF1.h" +#include "TCanvas.h" +#include "TFile.h" #include "TROOT.h" // --- AliRoot header files --- @@ -50,16 +53,9 @@ ClassImp(AliPHOSRawFitterv2) //----------------------------------------------------------------------------- AliPHOSRawFitterv2::AliPHOSRawFitterv2(): AliPHOSRawFitterv0(), - fNtimeSamples(25), - fRMScut(11.) + fAlpha(0.1),fBeta(0.035),fMax(0) { //Default constructor. - fLGpar[0]=0.971 ; - fLGpar[1]=0.0465; - fLGpar[2]=1.56 ; - fHGpar[0]=0.941 ; - fHGpar[1]=0.0436; - fHGpar[2]=1.65 ; } //----------------------------------------------------------------------------- @@ -70,30 +66,16 @@ AliPHOSRawFitterv2::~AliPHOSRawFitterv2() //----------------------------------------------------------------------------- AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ): - AliPHOSRawFitterv0(phosFitter), - fNtimeSamples(25), - fRMScut(11.) + AliPHOSRawFitterv0(phosFitter), + fAlpha(0.1),fBeta(0.035),fMax(0) { //Copy constructor. - fNtimeSamples=phosFitter.fNtimeSamples ; - for(Int_t i=0; i<3;i++){ - fLGpar[i]=phosFitter.fLGpar[i] ; - fHGpar[i]=phosFitter.fHGpar[i] ; - } - fRMScut=phosFitter.fRMScut ; } //----------------------------------------------------------------------------- -AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter) +AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 & /*phosFitter*/) { //Assignment operator. - - fNtimeSamples=phosFitter.fNtimeSamples ; - for(Int_t i=0; i<3;i++){ - fLGpar[i]=phosFitter.fLGpar[i] ; - fHGpar[i]=phosFitter.fHGpar[i] ; - } - fRMScut=phosFitter.fRMScut ; return *this; } @@ -107,47 +89,42 @@ Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t si //reasonable shape, fits it with Gamma2 function and extracts //energy and time. - if (fNBunches > 1) { - fQuality = 1000; - return kTRUE; - } - + const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it const Float_t kBaseLine = 1.0; const Int_t kPreSamples = 10; - Float_t pedMean = 0; - Float_t pedRMS = 0; - Int_t nPed = 0; - UShort_t maxSample = 0; - Int_t nMax = 0; - - TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ; - if(!cs) - cs = new TCanvas("CSample","CSample") ; + fOverflow = kFALSE ; + fEnergy=0 ; + if (fCaloFlag == 2 || fNBunches > 1) { + fQuality = 150; + return kTRUE; + } + if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample + fQuality=200; + fEnergy=0 ; + return kTRUE; + } - TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ; - if(!h) h = new TH1D("hSample","",200,0.,200.) ; + //Evaluate pedestals + Float_t pedMean = 0; + Float_t pedRMS = 0; + Int_t nPed = 0; + for (Int_t i=sigLength-kPreSamples; i maxSample) maxSample = signal[i]; - if(signal[i] == maxSample) nMax++; - h->SetBinContent(i+1,signal[i]) ; - } - fEnergy = (Double_t)maxSample; - if (maxSample > 900 && nMax > 2) fOverflow = kTRUE; if (fPedSubtract) { if (nPed > 0) { fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ; if(fPedestalRMS > 0.) fPedestalRMS = TMath::Sqrt(fPedestalRMS) ; - fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction + pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction } else return kFALSE; @@ -155,96 +132,320 @@ Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t si else { //take pedestals from DB pedestal = (Double_t) fAmpOffset ; - if (fCalibData) { - Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ; - Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ; - pedestal += truePed - altroSettings ; - } - else{ - AliWarning(Form("Can not read data from OCDB")) ; - } - fEnergy-=pedestal ; } - if (fEnergy < kBaseLine) { - fEnergy = 0; - return kTRUE; - } - // calculate time - fTime=0. ; - Double_t tRMS = 0. ; - Double_t tW = 0. ; - Int_t cnts=0 ; - Double_t a=0,b=0,c=0 ; - if(fCaloFlag == 0){ // Low gain - a=fLGpar[0] ; - b=fLGpar[1] ; - c=fLGpar[2] ; + //calculate rough quality of the sample and check for overflow + Int_t maxBin=0 ; + Int_t maxAmp=0 ; + Int_t minAmp= signal[0] ; + Int_t nMax = 0 ; //number of points in plato + Double_t aMean =0. ; + Double_t aRMS =0. ; + Double_t wts =0 ; + Bool_t falling = kTRUE ; //Bad monotoneusly falling sample + Bool_t rising = kTRUE ; //Bad monotoneusly riging sample + for (Int_t i=0; i pedestal){ + Double_t de = signal[i] - pedestal ; + if(de > 1.) { + aMean += de*i ; + aRMS += de*i*i ; + wts += de; + } + if(signal[i] > maxAmp){ + maxAmp = signal[i]; + nMax=0; + maxBin = i ; + } + if(signal[i] == maxAmp){ + nMax++; + } + if(signal[i] < minAmp) + minAmp=signal[i] ; + if(falling && i>0 && signal[i]0 && signal[i]>signal[i-1]) + rising=kFALSE ; + } + } + + if(rising || falling){//bad "rising" or falling sample + fEnergy = 0. ; + fTime = 0. ; //-999. ; + fQuality= 250. ; + return kTRUE ; } - else if(fCaloFlag == 1){ // High gain - a=fHGpar[0] ; - b=fHGpar[1] ; - c=fHGpar[2] ; + if(maxAmp-minAmp<3 && maxAmp>7 && sigLength>20){ //bad flat sample + fEnergy = 0. ; + fTime = 0; //-999. ; + fQuality= 260. ; + return kTRUE ; } + fEnergy=Double_t(maxAmp)-pedestal ; + if (fEnergy < kBaseLine) fEnergy = 0; + fTime = sigStart-sigLength-3; + + //do not test quality of too soft samples + if (wts > 0) { + aMean /= wts; + aRMS = aRMS/wts - aMean*aMean; + } + if (fEnergy <= maxEtoFit){ + if (aRMS < 2.) //sigle peak + fQuality = 299. ; + else + fQuality = 0. ; + //Evaluate time of signal arriving + return kTRUE ; + } - fQuality = 0. ; + //look for plato on the top of sample + if (fEnergy>500 && //this is not fluctuation of soft sample + nMax>2){ //and there is a plato + fOverflow = kTRUE ; + } - for(Int_t i=1; i0.){ - fTime/=tW ; - fQuality = tRMS/tW-fTime*fTime ; - fTime+=sigStart; - } - else{ - fTime=-999. ; - fQuality=999. ; - } - - Bool_t isBad = 0 ; - for(Int_t i=1; i signal[i-1]+5 && signal[i] > signal[i+1]+5) { //single jump - isBad=1 ; + + //do not fit High Gain samples with overflow + if(fCaloFlag==1 && fOverflow){ + fQuality = 99. ; + return kTRUE; + + } + + //----Now fit sample with reasonable shape------ + TArrayD samples(sigLength); // array of sample amplitudes + TArrayD times(sigLength); // array of sample time stamps + for (Int_t i=0; i3){ + goto plot ; + } + else{ + return kFALSE ; } } - if(pedestal < 10.) - isBad=1 ; - - if(fPedestalRMS > 0.1) - isBad=1 ; + fEnergy*=fMax ; + fTime += sigStart-sigLength-3; + + + //Impose cut on quality +// fQuality/=4. ; + fQuality/=1.+0.005*fEnergy ; + + //Draw corrupted samples + if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){ + if(fEnergy > 50. ){ + plot: + printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ; + TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ; + if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ; + h->Reset() ; + for (Int_t i=0; iSetBinContent(i+1,float(samples.At(i))) ; + } +// TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ; + TF1 * fffit = new TF1("fffit","[0]*((x-[1])*(x-[1])*exp(-[2]*(x-[1]))+(x-[1])*exp(-[3]*(x-[1])))",0.,60.) ; + fffit->SetParameters(fEnergy/fMax,fTime-(sigStart-sigLength-3),fAlpha,fBeta) ; + fffit->SetLineColor(2) ; + TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ; + if(!can){ + can = new TCanvas("cSamples","cSamples",10,10,600,600) ; + can->SetFillColor(0) ; + can->SetFillStyle(0) ; + can->Range(0,0,1,1); + can->SetBorderSize(0); + } + can->cd() ; - for(Int_t i=1; iDraw(); + spectrum_1->cd(); + spectrum_1->Range(0,0,1,1); + spectrum_1->SetFillColor(0); + spectrum_1->SetFillStyle(4000); + spectrum_1->SetBorderSize(1); + spectrum_1->SetBottomMargin(0.012); + spectrum_1->SetTopMargin(0.03); + spectrum_1->SetLeftMargin(0.10); + spectrum_1->SetRightMargin(0.05); + + char title[155] ; + sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ; + h->SetTitle(title) ; +// h->Fit(fffit,"","",0.,51.) ; + h->Draw() ; + fffit->Draw("same") ; +/* + sprintf(title,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ; + TFile fout("samples_bad.root","update") ; + h->Write(title); + fout.Close() ; +*/ + can->cd() ; + TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32); + spectrum_2->SetFillColor(0) ; + spectrum_2->SetFillStyle(0) ; + spectrum_2->SetGridy() ; + spectrum_2->Draw(); + spectrum_2->Range(0,0,1,1); + spectrum_2->SetFillColor(0); + spectrum_2->SetBorderSize(1); + spectrum_2->SetTopMargin(0.01); + spectrum_2->SetBottomMargin(0.25); + spectrum_2->SetLeftMargin(0.10); + spectrum_2->SetRightMargin(0.05); + spectrum_2->cd() ; + + TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ; + if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ; + hd->Reset() ; + for (Int_t i=0; iSetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ; + } + hd->Draw(); + + can->Update() ; + printf("Press to continue\n") ; + getchar(); + + + delete fffit ; + delete spectrum_1 ; + delete spectrum_2 ; + } } - if(fEnergy>10. && !isBad ){ - printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ; - if(fOverflow)printf(" Overflow \n") ; - if(isBad)printf("bad") ; - cs->cd() ; - h->Draw() ; - cs->Update() ; - getchar() ; - } - return kTRUE; } +//------------------------------------------------------------------ +Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){ +// makes fit + + const Int_t nMaxIter=50 ; //Maximal number of iterations + const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation + + Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation +//printf(" start fit... \n") ; + + Int_t nPoints = samples.GetSize() ; + Double_t dea=TMath::Exp(-fAlpha) ; + Double_t deb=TMath::Exp(-fBeta) ; + Double_t dt=1.,timeOld=dTime,dfOld=0. ; + for(Int_t iter=0; iter0 + if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size + dt*=0.5 ; + dTime=timeOld+dt ; + continue ; + } + if(f<0){ //f<0 => dTime is too small and we still do not know root region + dTime+=2. ; + continue ; + } + else{ //dTime is too large, we are beyond the root region + dTime-=2. ; + continue ; + } + } + dt=-f/df ; + if(TMath::Abs(dt)10.) dt=10. ; //restrict step size + if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another + timeOld=dTime ; //remember current position for the case + dfOld=df ; //of reduction of dt step size + dTime+=dt ; + + if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy. + fQuality=(yy-yf*yf/ff)/nfit ; + fEnergy=yf/ff ; //ff!=0 already tested + fTime=dTime ; + return kFALSE ; + } + + } + //failed to find a root, too many iterations + fQuality=99.; + fEnergy=0 ; + return kFALSE ; +} +//_________________________________________ +void AliPHOSRawFitterv2::FindMax(){ + //Finds maxumum of currecnt parameterization + Double_t t=2./fAlpha ; + fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ; + Double_t dt=15 ; + while(dt>0.01){ + Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ; + if(dfdt>0.) + t+=dt ; + else + t-=dt ; + Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ; + if(maxNew>fMax) + fMax=maxNew ; + else{ + dt/=2 ; + if(dfdt<0.) + t+=dt ; + else + t-=dt ; + } + } +} + diff --git a/PHOS/AliPHOSRawFitterv2.h b/PHOS/AliPHOSRawFitterv2.h index 997ed356060..d0c30921d43 100644 --- a/PHOS/AliPHOSRawFitterv2.h +++ b/PHOS/AliPHOSRawFitterv2.h @@ -5,12 +5,11 @@ /* $Id: $ */ -// This class extracts the PHOS "digits" of current event -// (amplitude,time, position,gain) from the raw stream -// provided by AliRawReader. See cxx source for use case. +// This class extracts amplitude, t0 and quality of the PHOS "samples" +// ising FastFit and two-exponent parameterization #include "AliPHOSRawFitterv0.h" -class TList; +class TArrayD ; class AliPHOSRawFitterv2 : public AliPHOSRawFitterv0 { @@ -22,19 +21,18 @@ public: virtual ~AliPHOSRawFitterv2(); virtual Bool_t Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength); + void SetRawParams(Double_t alpha, Double_t beta){fAlpha=alpha; fBeta=beta;} - void SetNTimeSamples(Short_t n=25) { fNtimeSamples=n ;} - void SetLowGainTParams (Double_t *pars){ for(Int_t i=0;i<3;i++) fLGpar[i]=pars[i] ;} - void SetHighGainTParams(Double_t *pars){ for(Int_t i=0;i<3;i++) fHGpar[i]=pars[i] ;} - void SetRMScut(Double_t cut=2.) { fRMScut = cut ;} +private: + Bool_t FindAmpT(TArrayD samples, TArrayD times) ; + void FindMax() ; private: - Short_t fNtimeSamples ; //Number of samples (after start) used to extract time - Double_t fLGpar[3] ; //parameters for shape parameterization - Double_t fHGpar[3] ; //parameters for shape parameterization - Double_t fRMScut ; //cut to estmate goodness of sample + Double_t fAlpha ; //Parameter of sample shape + Double_t fBeta ; //Parameter of sample shape + Double_t fMax ; //Maximum of parameterization - ClassDef(AliPHOSRawFitterv2,1) + ClassDef(AliPHOSRawFitterv2,2) }; #endif -- 2.43.0