Updates for the ROOT trunk
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv2.cxx
1 /**************************************************************************
2  * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16 /* $Id: $ */
17
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%
23 // 
24 // Typical use case:
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;
32
33 // Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
34
35 // --- ROOT system ---
36 #include "TArrayI.h"
37 #include "TList.h"
38 #include "TMath.h"
39 #include "TH1I.h"
40 #include "TF1.h"
41 #include "TCanvas.h"
42 #include "TFile.h"
43 #include "TROOT.h"
44
45 // --- AliRoot header files ---
46 #include "AliLog.h"
47 #include "AliPHOSCalibData.h"
48 #include "AliPHOSRawFitterv2.h"
49 #include "AliPHOSPulseGenerator.h"
50
51 ClassImp(AliPHOSRawFitterv2)
52
53 //-----------------------------------------------------------------------------
54 AliPHOSRawFitterv2::AliPHOSRawFitterv2():
55   AliPHOSRawFitterv0(),
56   fAlpha(0.1),fBeta(0.035),fMax(0) 
57 {
58   //Default constructor.
59 }
60
61 //-----------------------------------------------------------------------------
62 AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
63 {
64   //Destructor.
65 }
66
67 //-----------------------------------------------------------------------------
68 AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
69   AliPHOSRawFitterv0(phosFitter),
70   fAlpha(0.1),fBeta(0.035),fMax(0)
71 {
72   //Copy constructor.
73 }
74
75 //-----------------------------------------------------------------------------
76 AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 & /*phosFitter*/)
77 {
78   //Assignment operator.
79   return *this;
80 }
81
82 //-----------------------------------------------------------------------------
83 Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
84 {
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 
90   //energy and time.
91
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;
95
96   fOverflow = kFALSE ;  
97   fEnergy=0 ;
98   if (fCaloFlag == 2 || fNBunches > 1) {
99     fQuality = 150;
100     return kTRUE;
101   }
102   if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
103     fQuality=200;
104     fEnergy=0 ;
105     return kTRUE;
106   }
107
108   //Evaluate pedestals 
109   Float_t pedMean = 0;
110   Float_t pedRMS  = 0;
111   Int_t   nPed    = 0;
112   for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
113     nPed++;
114     pedMean += signal[i];
115     pedRMS  += signal[i]*signal[i] ;
116   }
117
118   fEnergy = -111;
119   fQuality= 999. ;
120   Double_t pedestal = 0;
121
122   if (fPedSubtract) {
123     if (nPed > 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
128     }
129     else
130       return kFALSE;
131   }
132   else {
133     //take pedestals from DB
134     pedestal = (Double_t) fAmpOffset ;
135   }
136
137   
138   //calculate rough quality of the sample and check for overflow
139   Int_t    maxBin=0 ;
140   Int_t    maxAmp=0 ;
141   Int_t    minAmp= signal[0] ;
142   Int_t    nMax = 0 ; //number of points in plato
143   Double_t aMean =0. ;
144   Double_t aRMS  =0. ;
145   Double_t wts   =0 ;
146   Bool_t falling = kTRUE ; //Bad monotoneusly falling sample
147   Bool_t rising = kTRUE ; //Bad monotoneusly riging sample
148   for (Int_t i=0; i<sigLength; i++){
149     if(signal[i] > pedestal){
150       Double_t de = signal[i] - pedestal ;
151       if(de > 1.) {
152         aMean += de*i ;
153         aRMS  += de*i*i ;
154         wts   += de; 
155       }
156       if(signal[i] >  maxAmp){
157         maxAmp = signal[i]; 
158         nMax=0;
159         maxBin = i ;
160       }
161       if(signal[i] == maxAmp){
162         nMax++;
163       }
164       if(signal[i] <  minAmp)
165         minAmp=signal[i] ;
166       if(falling && i>0 && signal[i]<signal[i-1])
167         falling=kFALSE ;
168       if(rising && i>0 && signal[i]>signal[i-1])
169         rising=kFALSE ;
170     }
171   }
172
173   if(rising || falling){//bad "rising" or falling  sample
174     fEnergy =    0. ;
175     fTime   = 0. ; //-999. ;
176     fQuality=  250. ;
177     return kTRUE ;
178   }
179   if(maxAmp-minAmp<3 && maxAmp>7 && sigLength>20){ //bad flat sample
180     fEnergy =    0. ;
181     fTime   = 0; //-999. ;
182     fQuality=  260. ;
183     return kTRUE ;
184   }
185
186   fEnergy=Double_t(maxAmp)-pedestal ;
187   if (fEnergy < kBaseLine) fEnergy = 0;
188   fTime = sigStart-sigLength-3;
189
190   //do not test quality of too soft samples
191   if (wts > 0) {
192     aMean /= wts; 
193     aRMS   = aRMS/wts - aMean*aMean;
194   }
195   if (fEnergy <= maxEtoFit){
196     if (aRMS < 2.) //sigle peak
197       fQuality = 299. ;
198     else
199       fQuality =   0. ;
200     //Evaluate time of signal arriving
201     return kTRUE ;
202   }
203
204   //look for plato on the top of sample
205   if (fEnergy>500 &&  //this is not fluctuation of soft sample
206      nMax>2){ //and there is a plato
207     fOverflow = kTRUE ;
208   }
209   
210
211   //do not fit High Gain samples with overflow
212   if(fCaloFlag==1 && fOverflow){
213     fQuality = 99. ;
214     return kTRUE;
215
216   }
217
218   //----Now fit sample with reasonable shape------
219   TArrayD samples(sigLength); // array of sample amplitudes
220   TArrayD times(sigLength); // array of sample time stamps
221   for (Int_t i=0; i<sigLength; i++) {
222     samples.AddAt(signal[i]-pedestal,sigLength-i-1);
223     times.AddAt(double(i),i);
224   }
225       
226   if(fMax==0)
227     FindMax() ;
228   if(!FindAmpT(samples,times)){
229     if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
230       goto plot ;
231     }
232     else{
233       return kFALSE ;
234     }
235   }
236   fEnergy*=fMax ;
237   fTime += sigStart-sigLength-3;
238
239
240   //Impose cut on quality
241 //  fQuality/=4. ;
242   fQuality/=1.+0.005*fEnergy ;
243
244   //Draw corrupted samples
245   if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
246     if(fEnergy > 50. ){
247     plot:
248       printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
249       TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
250       if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
251       h->Reset() ;
252       for (Int_t i=0; i<sigLength; i++) {
253         h->SetBinContent(i+1,float(samples.At(i))) ;
254       }
255 //      TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
256       TF1 * fffit = new TF1("fffit","[0]*((x-[1])*(x-[1])*exp(-[2]*(x-[1]))+(x-[1])*exp(-[3]*(x-[1])))",0.,60.) ;
257       fffit->SetParameters(fEnergy/fMax,fTime-(sigStart-sigLength-3),fAlpha,fBeta) ;
258       fffit->SetLineColor(2) ;
259       TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
260       if(!can){
261         can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
262         can->SetFillColor(0) ;
263         can->SetFillStyle(0) ;
264         can->Range(0,0,1,1);
265         can->SetBorderSize(0);
266       }
267       can->cd() ;
268   
269       TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
270       spectrum_1->Draw();
271       spectrum_1->cd();
272       spectrum_1->Range(0,0,1,1);
273       spectrum_1->SetFillColor(0);
274       spectrum_1->SetFillStyle(4000);
275       spectrum_1->SetBorderSize(1);
276       spectrum_1->SetBottomMargin(0.012);
277       spectrum_1->SetTopMargin(0.03);
278       spectrum_1->SetLeftMargin(0.10);
279       spectrum_1->SetRightMargin(0.05);
280
281       char title[155] ;
282       snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
283       h->SetTitle(title) ;
284 //      h->Fit(fffit,"","",0.,51.) ;
285       h->Draw() ;
286       fffit->Draw("same") ;
287 /*
288       sprintf(title,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
289       TFile fout("samples_bad.root","update") ;
290       h->Write(title);
291       fout.Close() ;
292 */
293       can->cd() ;
294       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
295       spectrum_2->SetFillColor(0) ;
296       spectrum_2->SetFillStyle(0) ;
297       spectrum_2->SetGridy() ;
298       spectrum_2->Draw();
299       spectrum_2->Range(0,0,1,1);
300       spectrum_2->SetFillColor(0);
301       spectrum_2->SetBorderSize(1);
302       spectrum_2->SetTopMargin(0.01);
303       spectrum_2->SetBottomMargin(0.25);
304       spectrum_2->SetLeftMargin(0.10);
305       spectrum_2->SetRightMargin(0.05);
306       spectrum_2->cd() ;
307
308       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
309       if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
310       hd->Reset() ;
311       for (Int_t i=0; i<sigLength; i++) {
312         hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
313       }
314       hd->Draw();
315  
316       can->Update() ;
317       printf("Press <enter> to continue\n") ;
318       getchar();
319
320
321       delete fffit ;
322       delete spectrum_1 ;
323       delete spectrum_2 ;
324     }
325   }
326   
327   return kTRUE;
328 }
329 //------------------------------------------------------------------
330 Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){
331 // makes fit
332
333   const Int_t nMaxIter=50 ;   //Maximal number of iterations
334   const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation
335
336   Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation
337 //printf(" start fit... \n") ;
338
339   Int_t nPoints = samples.GetSize() ;
340   Double_t dea=TMath::Exp(-fAlpha) ;
341   Double_t deb=TMath::Exp(-fBeta) ;
342   Double_t dt=1.,timeOld=dTime,dfOld=0. ; 
343   for(Int_t iter=0; iter<nMaxIter; iter++){
344     Double_t yy=0.;
345     Double_t yf=0. ;
346     Double_t ydf=0. ;
347     Double_t yddf=0. ;
348     Double_t ff=0. ;
349     Double_t fdf=0. ;
350     Double_t dfdf=0. ;
351     Double_t fddf=0. ;
352     Int_t nfit=0 ;
353     Double_t aexp=TMath::Exp(-fAlpha*(times.At(0)-1.-dTime)) ;
354     Double_t bexp=TMath::Exp(-fBeta*(times.At(0)-1.-dTime)) ;
355     for(Int_t i=0; i<nPoints; i++){
356       Double_t t= times.At(i)-dTime ;
357       aexp*=dea ;
358       bexp*=deb ;
359       if(t<0.) continue ;
360       Double_t y=samples.At(i) ;
361       if(y<=fAmpThreshold)
362         continue ;
363       nfit++ ;
364       Double_t at=fAlpha*t ;
365       Double_t bt = fBeta*t ;
366       Double_t phi=t*(t*aexp+bexp) ;
367       Double_t dphi=t*aexp*(2.-at)+(1.-bt)*bexp ;
368       Double_t ddphi=aexp*(2.-4.*at+at*at)+bexp*fBeta*(bt-2.) ;
369       yy+=y*y ;
370       yf+=y*phi ;
371       ydf+=y*dphi ;
372       yddf+=y*ddphi ;
373       ff+=phi*phi ;
374       fdf+=phi*dphi ;
375       dfdf+=dphi*dphi ;
376       fddf+=phi*ddphi ;
377     }
378
379     if(ff<1.e-09||nfit==0 ){
380       fQuality=199 ;
381       return kFALSE ;
382     }
383     Double_t f=ydf*ff-yf*fdf ;     //d(chi2)/dt
384     Double_t df=yf*(dfdf+fddf)-yddf*ff-ydf*fdf;
385     if(df<=0.){ //we are too far from the root. In the wicinity of root df>0
386       if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size
387         dt*=0.5 ;
388         dTime=timeOld+dt ;  
389         continue ;
390       }
391       if(f<0){ //f<0 => dTime is too small and we still do not know root region
392         dTime+=2. ;
393         continue ;
394       }
395       else{ //dTime is too large, we are beyond the root region
396         dTime-=2. ;
397         continue ;
398       }
399     }
400     dt=-f/df ; 
401     if(TMath::Abs(dt)<epsdt){
402       fQuality=(yy-yf*yf/ff)/nfit ;
403       fEnergy=yf/ff ;  //ff!=0 already tested
404       fTime=dTime ;
405       return kTRUE ;
406     }
407     //In some cases time steps are huge (derivative ~0)
408     if(dt>10.) dt=10. ;   //restrict step size
409     if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another
410     timeOld=dTime ;  //remember current position for the case
411     dfOld=df ;       //of reduction of dt step size
412     dTime+=dt ;
413
414     if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy.
415       fQuality=(yy-yf*yf/ff)/nfit ;
416       fEnergy=yf/ff ;  //ff!=0 already tested
417       fTime=dTime ;
418       return kFALSE ;
419     }
420
421   }
422   //failed to find a root, too many iterations
423   fQuality=99.;
424   fEnergy=0 ; 
425   return kFALSE ;
426 }
427 //_________________________________________
428 void AliPHOSRawFitterv2::FindMax(){
429   //Finds maxumum of currecnt parameterization
430   Double_t t=2./fAlpha ;
431   fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
432   Double_t dt=15 ;
433   while(dt>0.01){
434      Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ;
435      if(dfdt>0.)
436         t+=dt ;
437      else
438        t-=dt ;
439      Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
440      if(maxNew>fMax)
441         fMax=maxNew ;
442      else{
443        dt/=2 ;
444        if(dfdt<0.)
445          t+=dt ;
446        else
447         t-=dt ;
448      }
449   }   
450 }
451