1 /**************************************************************************
2 * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 // This class extracts the signal parameters (energy, time, quality)
19 // from ALTRO samples. Energy is in ADC counts, time is in time bin units.
20 // Class uses FastFitting algorithm to fit sample and extract time and Amplitude
21 // and evaluate sample quality = (chi^2/NDF)/some parameterization providing
22 // efficiency close to 100%
25 // AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
26 // fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
27 // fitter->SetCalibData(fgCalibData) ;
28 // fitter->Eval(sig,sigStart,sigLength);
29 // Double_t amplitude = fitter.GetEnergy();
30 // Double_t time = fitter.GetTime();
31 // Bool_t isLowGain = fitter.GetCaloFlag()==0;
33 // Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
35 // --- ROOT system ---
45 // --- AliRoot header files ---
47 #include "AliPHOSCalibData.h"
48 #include "AliPHOSRawFitterv2.h"
49 #include "AliPHOSPulseGenerator.h"
51 ClassImp(AliPHOSRawFitterv2)
53 //-----------------------------------------------------------------------------
54 AliPHOSRawFitterv2::AliPHOSRawFitterv2():
56 fAlpha(0.1),fBeta(0.035),fMax(0)
58 //Default constructor.
61 //-----------------------------------------------------------------------------
62 AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
67 //-----------------------------------------------------------------------------
68 AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
69 AliPHOSRawFitterv0(phosFitter),
70 fAlpha(0.1),fBeta(0.035),fMax(0)
75 //-----------------------------------------------------------------------------
76 AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 & /*phosFitter*/)
78 //Assignment operator.
82 //-----------------------------------------------------------------------------
83 Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
85 //Extract an energy deposited in the crystal,
86 //crystal' position (module,column,row),
87 //time and gain (high or low).
88 //First collects sample, then evaluates it and if it has
89 //reasonable shape, fits it with Gamma2 function and extracts
92 const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
93 const Float_t kBaseLine = 1.0;
94 const Int_t kPreSamples = 10;
98 if (fCaloFlag == 2 || fNBunches > 1) {
102 if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
112 for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
114 pedMean += signal[i];
115 pedRMS += signal[i]*signal[i] ;
120 Double_t pedestal = 0;
124 fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
125 if(fPedestalRMS > 0.)
126 fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
127 pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
133 //take pedestals from DB
134 pedestal = (Double_t) fAmpOffset ;
138 //calculate rough quality of the sample and check for overflow
140 Int_t minAmp= signal[0] ;
141 Int_t nMax = 0 ; //number of points in plato
145 Bool_t falling = kTRUE ; //Bad monotoneusly falling sample
146 Bool_t rising = kTRUE ; //Bad monotoneusly riging sample
147 for (Int_t i=0; i<sigLength; i++){
148 if(signal[i] > pedestal){
149 Double_t de = signal[i] - pedestal ;
155 if(signal[i] > maxAmp){
159 if(signal[i] == maxAmp){
162 if(signal[i] < minAmp)
164 if(falling && i>0 && signal[i]<signal[i-1])
166 if(rising && i>0 && signal[i]>signal[i-1])
171 if(rising || falling){//bad "rising" or falling sample
173 fTime = 0. ; //-999. ;
177 if(maxAmp-minAmp<3 && maxAmp>7 && sigLength>20){ //bad flat sample
184 fEnergy=Double_t(maxAmp)-pedestal ;
185 if (fEnergy < kBaseLine) fEnergy = 0;
186 fTime = sigStart-sigLength-3;
188 //do not test quality of too soft samples
191 aRMS = aRMS/wts - aMean*aMean;
193 if (fEnergy <= maxEtoFit){
194 if (aRMS < 2.) //sigle peak
198 //Evaluate time of signal arriving
202 //look for plato on the top of sample
203 if (fEnergy>500 && //this is not fluctuation of soft sample
204 nMax>2){ //and there is a plato
209 //do not fit High Gain samples with overflow
210 if(fCaloFlag==1 && fOverflow){
216 //----Now fit sample with reasonable shape------
217 TArrayD samples(sigLength); // array of sample amplitudes
218 TArrayD times(sigLength); // array of sample time stamps
219 for (Int_t i=0; i<sigLength; i++) {
220 samples.AddAt(signal[i]-pedestal,sigLength-i-1);
221 times.AddAt(double(i),i);
226 if(!FindAmpT(samples,times)){
227 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
235 fTime += sigStart-sigLength-3;
238 //Impose cut on quality
240 fQuality/=1.+0.005*fEnergy ;
242 //Draw corrupted samples
243 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
246 printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
247 TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
248 if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
250 for (Int_t i=0; i<sigLength; i++) {
251 h->SetBinContent(i+1,float(samples.At(i))) ;
253 // TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
254 TF1 * fffit = new TF1("fffit","[0]*((x-[1])*(x-[1])*exp(-[2]*(x-[1]))+(x-[1])*exp(-[3]*(x-[1])))",0.,60.) ;
255 fffit->SetParameters(fEnergy/fMax,fTime-(sigStart-sigLength-3),fAlpha,fBeta) ;
256 fffit->SetLineColor(2) ;
257 TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
259 can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
260 can->SetFillColor(0) ;
261 can->SetFillStyle(0) ;
263 can->SetBorderSize(0);
267 TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
270 spectrum_1->Range(0,0,1,1);
271 spectrum_1->SetFillColor(0);
272 spectrum_1->SetFillStyle(4000);
273 spectrum_1->SetBorderSize(1);
274 spectrum_1->SetBottomMargin(0.012);
275 spectrum_1->SetTopMargin(0.03);
276 spectrum_1->SetLeftMargin(0.10);
277 spectrum_1->SetRightMargin(0.05);
280 snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
282 // h->Fit(fffit,"","",0.,51.) ;
284 fffit->Draw("same") ;
286 sprintf(title,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
287 TFile fout("samples_bad.root","update") ;
292 TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
293 spectrum_2->SetFillColor(0) ;
294 spectrum_2->SetFillStyle(0) ;
295 spectrum_2->SetGridy() ;
297 spectrum_2->Range(0,0,1,1);
298 spectrum_2->SetFillColor(0);
299 spectrum_2->SetBorderSize(1);
300 spectrum_2->SetTopMargin(0.01);
301 spectrum_2->SetBottomMargin(0.25);
302 spectrum_2->SetLeftMargin(0.10);
303 spectrum_2->SetRightMargin(0.05);
306 TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
307 if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
309 for (Int_t i=0; i<sigLength; i++) {
310 hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
315 printf("Press <enter> to continue\n") ;
327 //------------------------------------------------------------------
328 Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){
331 const Int_t nMaxIter=50 ; //Maximal number of iterations
332 const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation
334 Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation
335 //printf(" start fit... \n") ;
337 Int_t nPoints = samples.GetSize() ;
338 Double_t dea=TMath::Exp(-fAlpha) ;
339 Double_t deb=TMath::Exp(-fBeta) ;
340 Double_t dt=1.,timeOld=dTime,dfOld=0. ;
341 for(Int_t iter=0; iter<nMaxIter; iter++){
351 Double_t aexp=TMath::Exp(-fAlpha*(times.At(0)-1.-dTime)) ;
352 Double_t bexp=TMath::Exp(-fBeta*(times.At(0)-1.-dTime)) ;
353 for(Int_t i=0; i<nPoints; i++){
354 Double_t t= times.At(i)-dTime ;
358 Double_t y=samples.At(i) ;
362 Double_t at=fAlpha*t ;
363 Double_t bt = fBeta*t ;
364 Double_t phi=t*(t*aexp+bexp) ;
365 Double_t dphi=t*aexp*(2.-at)+(1.-bt)*bexp ;
366 Double_t ddphi=aexp*(2.-4.*at+at*at)+bexp*fBeta*(bt-2.) ;
377 if(ff<1.e-09||nfit==0 ){
381 Double_t f=ydf*ff-yf*fdf ; //d(chi2)/dt
382 Double_t df=yf*(dfdf+fddf)-yddf*ff-ydf*fdf;
383 if(df<=0.){ //we are too far from the root. In the wicinity of root df>0
384 if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size
389 if(f<0){ //f<0 => dTime is too small and we still do not know root region
393 else{ //dTime is too large, we are beyond the root region
399 if(TMath::Abs(dt)<epsdt){
400 fQuality=(yy-yf*yf/ff)/nfit ;
401 fEnergy=yf/ff ; //ff!=0 already tested
405 //In some cases time steps are huge (derivative ~0)
406 if(dt>10.) dt=10. ; //restrict step size
407 if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another
408 timeOld=dTime ; //remember current position for the case
409 dfOld=df ; //of reduction of dt step size
412 if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy.
413 fQuality=(yy-yf*yf/ff)/nfit ;
414 fEnergy=yf/ff ; //ff!=0 already tested
420 //failed to find a root, too many iterations
425 //_________________________________________
426 void AliPHOSRawFitterv2::FindMax(){
427 //Finds maxumum of currecnt parameterization
428 Double_t t=2./fAlpha ;
429 fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
432 Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ;
437 Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;