Merge branch 'master_patch'
[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=22.18 ;
111
112   TArrayD samples(sigLength); // array of sample amplitudes
113   TArrayD times(sigLength); // array of sample time stamps
114   for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
115     nPed++;
116     pedMean += signal[i];
117     pedRMS  += signal[i]*signal[i] ;
118   }
119
120   fEnergy = -111;
121   fQuality= 999. ;
122   const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
123   Double_t pedestal = 0;
124
125   if (fPedSubtract) {
126     if (nPed > 0) {
127       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
128       if(fPedestalRMS > 0.) 
129         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
130       pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
131       fEnergy -= pedestal; // pedestal subtraction
132     }
133     else
134       return kFALSE;
135   }
136   else {
137     //take pedestals from DB
138     pedestal = (Double_t) fAmpOffset ;
139     if (fCalibData) {
140       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
141       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
142       pedestal += truePed - altroSettings ;
143     }
144     else{
145 //      AliWarning(Form("Can not read data from OCDB")) ;
146     }
147     fEnergy-=pedestal ;
148   }
149
150   if (fEnergy < kBaseLine) fEnergy = 0;
151   //Evaluate time
152    fTime = sigStart-sigLength-3;
153   for (Int_t i=0; i<sigLength; i++) {
154     samples.AddAt(signal[i]-pedestal,sigLength-i-1);
155     times.AddAt(i/tau ,i);
156   }
157   
158   //calculate time and energy
159   Int_t    maxBin=0 ;
160   Int_t    maxAmp=0 ;
161   Int_t    nMax = 0 ; //number of points in plato
162   Double_t aMean =0. ;
163   Double_t aRMS  =0. ;
164   Double_t wts   =0 ;
165   for (Int_t i=0; i<sigLength; i++){
166     if(signal[i] > pedestal){
167       Double_t de = signal[i] - pedestal ;
168       if(de > 1.) {
169         aMean += de*i ;
170         aRMS  += de*i*i ;
171         wts   += de; 
172       }
173       if(signal[i] >  maxAmp){
174         maxAmp = signal[i]; 
175         nMax=0;
176         maxBin = i ;
177       }
178       if(signal[i] == maxAmp){
179         nMax++;
180       }
181     }
182   }
183
184   if (maxBin==sigLength-1){//bad "rising" sample
185     fEnergy =    0. ;
186     fTime   = -999. ;
187     fQuality=  999. ;
188     return kTRUE ;
189   }
190
191   fEnergy=Double_t(maxAmp)-pedestal ;
192   fOverflow =0 ;  //look for plato on the top of sample
193   if (fEnergy>500 &&  //this is not fluctuation of soft sample
194      nMax>2){ //and there is a plato
195     fOverflow = kTRUE ;
196   }
197   
198   if (wts > 0) {
199     aMean /= wts; 
200     aRMS   = aRMS/wts - aMean*aMean;
201   }
202
203   //do not take too small energies
204   if (fEnergy < kBaseLine) 
205     fEnergy = 0;
206   
207   //do not test quality of too soft samples
208   if (fEnergy < maxEtoFit){
209     if (aRMS < 2.) //sigle peak
210       fQuality = 999. ;
211     else
212       fQuality =   0. ;
213     //Evaluate time of signal arriving
214     return kTRUE ;
215   }
216       
217   // if sample has reasonable mean and RMS, try to fit it with gamma2
218   //This method can not analyse overflow samples
219   if(fOverflow){
220     fQuality = 99. ;
221     return kTRUE ;
222   }
223   // First estimate t0
224   Double_t a=0,b=0,c=0 ;
225   Int_t minI=0 ;
226   if (fPedSubtract) 
227     minI=kPreSamples ;
228   for(Int_t i=minI; i<sigLength; i++){
229     if(samples.At(i)<=0.)
230       continue ;
231     Double_t t= times.At(i) ;
232     Double_t f02 = TMath::Exp(-2.*t);
233     Double_t f12 = t*f02;
234     Double_t f22 = t*f12;
235     // Derivatives
236     Double_t f02d = -2.*f02;
237     Double_t f12d = f02 - 2.*f12;
238     Double_t f22d = 2.*(f12 - f22);
239     a += f02d * samples.At(i) ;
240     b -= f12d * samples.At(i) ;
241     c += f22d * samples.At(i) ;
242   }
243   
244   //Find roots
245   if(a==0.){
246     if(b==0.){ //no roots
247       fQuality = 2000. ;
248       if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
249         printf(" a=%f, b=%f, c=%f \n",a,b,c) ;
250         goto plot ;
251       }
252       return kTRUE ;
253     }
254     Double_t t1=-c/b ;
255     Double_t amp=0.,den=0.; ;
256     for(Int_t i=minI; i<sigLength; i++){
257       if(samples.At(i)<=0.)
258         continue ;
259       Double_t dt = times.At(i) - t1;
260       Double_t f = dt*dt*TMath::Exp(-2.*dt);
261       amp += f*samples.At(i);
262       den += f*f;
263     }
264     if(den>0.0) amp /= den;
265     // chi2 calculation
266     fQuality=0.;
267     for(Int_t i=minI; i<sigLength; i++){
268       if(samples.At(i)<=0.)
269         continue ;
270       Double_t t = times.At(i)- t1;
271       Double_t dy = samples.At(i)- amp*t*t*TMath::Exp(-2.*t) ;
272       fQuality += dy*dy;
273     }
274     fTime+=t1*tau ;
275     fEnergy = amp*TMath::Exp(-2.);
276     fQuality/= sigLength ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
277   }
278   else{
279     Double_t det = b*b - a*c;
280     if(det>=1.e-6 && det<0.0) {
281       det = 0.0; //rounding error
282     }
283     if(det<0.){ //Problem
284       fQuality = 1500. ;
285       if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
286         printf(" det=%e \n",det) ;
287         goto plot ;
288       }
289       return kTRUE ;
290     }
291
292     det = TMath::Sqrt(det);
293     Double_t t1 = (-b + det) / a;
294 //    Double_t t2 = (-b - det) / a; //second root is wrong one
295     Double_t amp1=0., den1=0. ;
296     for(Int_t i=minI; i<sigLength; i++){
297       if(samples.At(i)<=0.)
298         continue ;
299       Double_t dt1 = times.At(i) - t1;
300       Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
301       amp1 += f01*samples.At(i);
302       den1 += f01*f01;
303     }
304     if(den1>0.0) amp1 /= den1;
305     Double_t chi1=0.; // chi2=0. ;
306     for(Int_t i=minI; i<sigLength; i++){
307       if(samples.At(i)<=0.)
308         continue ;
309       Double_t dt1 = times.At(i)- t1;
310       Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
311       chi1 += dy1*dy1;
312     }
313     fEnergy=amp1*TMath::Exp(-2.) ; ; 
314     fTime+=t1*tau ;
315     fQuality=chi1/sigLength ;
316   } 
317
318   //Impose cut on quality
319   fQuality/=2.+0.004*fEnergy*fEnergy ;
320
321   //Draw corrupted samples
322   if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
323     plot:
324     if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples
325 //    if(!fOverflow ){ //Draw only bad samples
326       printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
327       TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
328       if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
329       h->Reset() ;
330       for (Int_t i=0; i<sigLength; i++) {
331         h->SetBinContent(i+1,samples.At(i)+pedestal) ;
332       }
333       TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
334       fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
335       fffit->SetLineColor(2) ;
336       TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
337       if(!can){
338         can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
339         can->SetFillColor(0) ;
340         can->SetFillStyle(0) ;
341         can->Range(0,0,1,1);
342         can->SetBorderSize(0);
343       }
344       can->cd() ;
345   
346       TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
347       spectrum_1->Draw();
348       spectrum_1->cd();
349       spectrum_1->Range(0,0,1,1);
350       spectrum_1->SetFillColor(0);
351       spectrum_1->SetFillStyle(4000);
352       spectrum_1->SetBorderSize(1);
353       spectrum_1->SetBottomMargin(0.012);
354       spectrum_1->SetTopMargin(0.03);
355       spectrum_1->SetLeftMargin(0.10);
356       spectrum_1->SetRightMargin(0.05);
357
358       char title[155] ;
359       snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
360       h->SetTitle(title) ;
361       h->Draw() ;
362       fffit->Draw("same") ;
363
364       snprintf(title,155,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
365       TFile fout("samples_bad.root","update") ;
366       h->Write(title);
367       fout.Close() ;
368
369       can->cd() ;
370       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
371       spectrum_2->SetFillColor(0) ;
372       spectrum_2->SetFillStyle(0) ;
373       spectrum_2->SetGridy() ;
374       spectrum_2->Draw();
375       spectrum_2->Range(0,0,1,1);
376       spectrum_2->SetFillColor(0);
377       spectrum_2->SetBorderSize(1);
378       spectrum_2->SetTopMargin(0.01);
379       spectrum_2->SetBottomMargin(0.25);
380       spectrum_2->SetLeftMargin(0.10);
381       spectrum_2->SetRightMargin(0.05);
382       spectrum_2->cd() ;
383
384       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
385       if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
386       hd->Reset() ;
387       for (Int_t i=0; i<sigLength; i++) {
388         hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
389       }
390       hd->Draw();
391 /* 
392       can->Update() ;
393       printf("Press <enter> to continue\n") ;
394       getchar();
395 */
396
397       delete fffit ;
398       delete spectrum_1 ;
399       delete spectrum_2 ;
400     }
401   }
402   
403   return kTRUE;
404 }