Raw fitter based on FastFitting algorithm
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv3.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 "AliPHOSRawFitterv3.h"
49 #include "AliPHOSPulseGenerator.h"
50
51 ClassImp(AliPHOSRawFitterv3)
52
53 //-----------------------------------------------------------------------------
54 AliPHOSRawFitterv3::AliPHOSRawFitterv3():
55   AliPHOSRawFitterv0()
56 {
57   //Default constructor.
58 }
59
60 //-----------------------------------------------------------------------------
61 AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
62 {
63   //Destructor.
64 }
65
66 //-----------------------------------------------------------------------------
67 AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
68   AliPHOSRawFitterv0(phosFitter) 
69 {
70   //Copy constructor.
71 }
72
73 //-----------------------------------------------------------------------------
74 AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
75 {
76   //Assignment operator.
77   return *this;
78 }
79
80 //-----------------------------------------------------------------------------
81 Bool_t AliPHOSRawFitterv3::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
82 {
83   //Extract an energy deposited in the crystal,
84   //crystal' position (module,column,row),
85   //time and gain (high or low).
86   //First collects sample, then evaluates it and if it has
87   //reasonable shape, fits it with Gamma2 function and extracts 
88   //energy and time.
89
90   if (fCaloFlag == 2 || fNBunches > 1) {
91     fQuality = 1000;
92     return kTRUE;
93   }
94   if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
95     fQuality=2000;
96     fEnergy=0 ;
97     return kTRUE;
98   }
99
100   Float_t pedMean = 0;
101   Float_t pedRMS  = 0;
102   Int_t   nPed    = 0;
103   const Float_t kBaseLine   = 1.0;
104   const Int_t   kPreSamples = 10;
105   
106
107   //We tryed individual taus for each channel, but
108   //this approach seems to be unstable. Much better results are obtaned with
109   //fixed decay time for all channels.
110   const Double_t tau=21. ;
111
112   TArrayD *samples = new TArrayD(sigLength); // array of sample amplitudes
113   TArrayD *times   = new TArrayD(sigLength); // array of sample time stamps
114   for (Int_t i=0; i<sigLength; i++) {
115     if (i<kPreSamples) {
116       nPed++;
117       pedMean += signal[i];
118       pedRMS  += signal[i]*signal[i] ;
119     }
120     samples->AddAt(signal[i],sigLength-i-1);
121     times  ->AddAt(i/tau ,i);
122   }
123
124   fEnergy = -111;
125   fQuality= 999. ;
126   const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
127   Double_t pedestal = 0;
128
129   if (fPedSubtract) {
130     if (nPed > 0) {
131       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
132       if(fPedestalRMS > 0.) 
133         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
134       fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
135     }
136     else
137       return kFALSE;
138   }
139   else {
140     //take pedestals from DB
141     pedestal = (Double_t) fAmpOffset ;
142     if (fCalibData) {
143       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
144       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
145       pedestal += truePed - altroSettings ;
146     }
147     else{
148       AliWarning(Form("Can not read data from OCDB")) ;
149     }
150     fEnergy-=pedestal ;
151   }
152
153   if (fEnergy < kBaseLine) fEnergy = 0;
154   //Evaluate time
155   Int_t iStart = 0;
156   while(iStart<sigLength && samples->At(iStart)-pedestal <kBaseLine) iStart++ ;
157   fTime = sigStart+iStart; 
158   
159   //calculate time and energy
160   Int_t    maxBin=0 ;
161   Int_t    maxAmp=0 ;
162   Int_t    nMax = 0 ; //number of points in plato
163   Double_t aMean =0. ;
164   Double_t aRMS  =0. ;
165   Double_t wts   =0 ;
166   Int_t    tStart=0 ;
167
168   for (Int_t i=0; i<sigLength; i++){
169     if(signal[i] > pedestal){
170       Double_t de = signal[i] - pedestal ;
171       if(de > 1.) {
172         aMean += de*i ;
173         aRMS  += de*i*i ;
174         wts   += de; 
175       }
176       if(de > 2 && tStart==0) 
177         tStart = i ;
178       if(signal[i] >  maxAmp){
179         maxAmp = signal[i]; 
180         nMax=0;
181         maxBin = i ;
182       }
183       if(signal[i] == maxAmp)
184         nMax++;
185     }
186   }
187
188   if (maxBin==sigLength-1){//bad "rising" sample
189     fEnergy =    0. ;
190     fTime   = -999. ;
191     fQuality=  999. ;
192     return kTRUE ;
193   }
194
195   fEnergy=Double_t(maxAmp)-pedestal ;
196   fOverflow =0 ;  //look for plato on the top of sample
197   if (fEnergy>500 &&  //this is not fluctuation of soft sample
198      nMax>2){ //and there is a plato
199     fOverflow = kTRUE ;
200   }
201   
202   if (wts > 0) {
203     aMean /= wts; 
204     aRMS   = aRMS/wts - aMean*aMean;
205   }
206
207   //do not take too small energies
208   if (fEnergy < kBaseLine) 
209     fEnergy = 0;
210   
211   //do not test quality of too soft samples
212   if (fEnergy < maxEtoFit){
213     fTime = tStart;
214     if (aRMS < 2.) //sigle peak
215       fQuality = 999. ;
216     else
217       fQuality =   0. ;
218     return kTRUE ;
219   }
220       
221   // if sample has reasonable mean and RMS, try to fit it with gamma2
222   // First estimate t0
223   Double_t a=0,b=0,c=0 ;
224   for(Int_t i=0; i<sigLength; i++){
225     if(samples->At(i)<pedestal)
226       continue ;
227     Double_t t= times->At(i) ;
228     Double_t f02 = TMath::Exp(-2.*t);
229     Double_t f12 = t*f02;
230     Double_t f22 = t*f12;
231     // Derivatives
232     Double_t f02d = -2.*f02;
233     Double_t f12d = f02 - 2.*f12;
234     Double_t f22d = 2.*(f12 - f22);
235     a += f02d * (samples->At(i)-pedestal) ;
236     b -= f12d * (samples->At(i)-pedestal) ;
237     c += f22d * (samples->At(i)-pedestal) ;
238   }
239   
240   //Find 2 roots
241   Double_t det = b*b - a*c;
242   if(det>=1.e-6 && det<0.0) {
243     det = 0.0; //rounding error
244   }
245   if(det<0.){ //Problem
246     fQuality = 1500. ;
247     if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
248       printf(" det=%e \n",det) ;
249       goto plot ;
250     }
251     return kTRUE ;
252   }
253
254   if(det>0.0 && a!=0.) {
255     det = TMath::Sqrt(det);
256     Double_t t1 = (-b + det) / a;
257 //    Double_t t2 = (-b - det) / a; //second root is wrong one
258     Double_t amp1=0., den1=0. ;
259     for(Int_t i=0; i<sigLength; i++){
260        Double_t dt1 = times->At(i) - t1;
261        Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
262        amp1 += f01*(samples->At(i)-pedestal);
263        den1 += f01*f01;
264      }
265      if(den1>0.0) amp1 /= den1;
266      Double_t chi1=0.; // chi2=0. ;
267      for(Int_t i=0; i<sigLength; i++){
268        Double_t dt1 = times->At(i)- t1;
269        Double_t dy1 = samples->At(i)-pedestal- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
270        chi1 += dy1*dy1;
271      }
272      fEnergy=amp1*TMath::Exp(-2.) ; ; 
273      fTime=t1*tau ;
274      fQuality=chi1/sigLength ;
275   } 
276   else { 
277     Double_t t1 ;
278     if(a!=0.)
279       t1=-b/a ;
280     else
281       t1=-c/b ;
282     Double_t amp=0.,den=0.; ;
283     for(Int_t i=0; i<sigLength; i++){
284        Double_t dt = times->At(i) - t1;
285        Double_t f = dt*dt*TMath::Exp(-2.*dt);
286        amp += f*samples->At(i);
287        den += f*f;
288      }
289      if(den>0.0) amp /= den;
290      // chi2 calculation
291      fQuality=0.;
292      for(Int_t i=0; i<sigLength; i++){
293        Double_t t = times->At(i)- t1;
294        Double_t dy = samples->At(i)- amp*t*t*TMath::Exp(-2.*t) ;
295        fQuality += dy*dy;
296      }
297      fTime=t1*tau ;
298      fEnergy = amp*TMath::Exp(-2.);
299      fQuality/= sigLength ;
300   }
301
302   //Impose cut on quality
303   fQuality/=2.+0.004*fEnergy*fEnergy ;
304
305   //Draw corrupted samples
306   if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
307     plot:
308     if(fQuality >1. && !fOverflow ){ //Draw only bad samples
309       printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
310       TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
311       if(!h) h = new TH1I("h","Samples",50,0.,50.) ;
312       h->Reset() ;
313       for (Int_t i=0; i<sigLength; i++) {
314         h->SetBinContent(i,samples->At(i)) ;
315       }
316       TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
317       fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
318       fffit->SetLineColor(2) ;
319       TCanvas * c = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
320       if(!c){
321         c = new TCanvas("cSamples","cSamples",10,10,400,400) ;
322         c->SetFillColor(0) ;
323         c->SetFillStyle(0) ;
324         c->Range(0,0,1,1);
325         c->SetBorderSize(0);
326       }
327       c->cd() ;
328   
329       TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
330       spectrum_1->Draw();
331       spectrum_1->cd();
332       spectrum_1->Range(0,0,1,1);
333       spectrum_1->SetFillColor(0);
334       spectrum_1->SetFillStyle(4000);
335       spectrum_1->SetBorderSize(1);
336       spectrum_1->SetBottomMargin(0.012);
337       spectrum_1->SetTopMargin(0.03);
338       spectrum_1->SetLeftMargin(0.10);
339       spectrum_1->SetRightMargin(0.05);
340
341       char title[155] ;
342       sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
343       h->SetTitle(title) ;
344       h->Draw() ;
345       fffit->Draw("same") ;
346
347       sprintf(title,"mod%d_x%d_z%d_HG",fModule,fCellX,fCellZ) ;
348       TFile fout("samples_bad.root","update") ;
349       h->Write(title);
350       fout.Close() ;
351
352       c->cd() ;
353       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
354       spectrum_2->SetFillColor(0) ;
355       spectrum_2->SetFillStyle(0) ;
356       spectrum_2->SetGridy() ;
357       spectrum_2->Draw();
358       spectrum_2->Range(0,0,1,1);
359       spectrum_2->SetFillColor(0);
360       spectrum_2->SetBorderSize(1);
361       spectrum_2->SetTopMargin(0.01);
362       spectrum_2->SetBottomMargin(0.25);
363       spectrum_2->SetLeftMargin(0.10);
364       spectrum_2->SetRightMargin(0.05);
365       spectrum_2->cd() ;
366
367       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
368       if(!hd) hd = new TH1I("hd","Samples",50,0.,50.) ;
369       hd->Reset() ;
370       for (Int_t i=0; i<sigLength; i++) {
371         hd->SetBinContent(i,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
372       }
373       hd->Draw();
374  
375       c->Update() ;
376       printf("Press <enter> to continue\n") ;
377       getchar();
378
379
380       delete fffit ;
381       delete spectrum_1 ;
382       delete spectrum_2 ;
383     }
384   }
385   
386   delete samples ;
387   delete times ;
388   return kTRUE;
389 }