Raw data decoder is migrated from AliCaloRawStream to AliCaloRawStreamV3
[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 plots samples and qualify their quality according to their shape
19 // 
20 // Typical use case:
21 //     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
22 //     fitter->SetSamples(sig,sigStart,sigLength);
23 //     fitter->SetNBunches(nBunches);
24 //     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
25 //     fitter->SetCalibData(fgCalibData) ;
26 //     fitter->Eval();
27 //     Double_t amplitude = fitter.GetEnergy();
28 //     Double_t time      = fitter.GetTime();
29 //     Bool_t   isLowGain = fitter.GetCaloFlag()==0;
30
31 // Author: Dmitri Peressounko (Oct.2008)
32 // Modified: Yuri Kharlov (Jul.2009)
33
34 // --- ROOT system ---
35 #include "TList.h"
36 #include "TMath.h"
37 #include "TMinuit.h"
38 #include "TCanvas.h"
39 #include "TH1.h"
40 #include "TH2.h"
41 #include "TF1.h"
42 #include "TROOT.h"
43
44 // --- AliRoot header files ---
45 #include "AliLog.h"
46 #include "AliPHOSCalibData.h"
47 #include "AliPHOSRawFitterv2.h"
48 #include "AliPHOSPulseGenerator.h"
49
50 ClassImp(AliPHOSRawFitterv2)
51
52 //-----------------------------------------------------------------------------
53 AliPHOSRawFitterv2::AliPHOSRawFitterv2():
54   AliPHOSRawFitterv0(),
55   fNtimeSamples(25),
56   fRMScut(11.)
57 {
58   //Default constructor.
59   fLGpar[0]=0.971 ;
60   fLGpar[1]=0.0465;
61   fLGpar[2]=1.56  ;
62   fHGpar[0]=0.941 ; 
63   fHGpar[1]=0.0436;
64   fHGpar[2]=1.65  ;
65 }
66
67 //-----------------------------------------------------------------------------
68 AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
69 {
70   //Destructor.
71 }
72
73 //-----------------------------------------------------------------------------
74 AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
75   AliPHOSRawFitterv0(phosFitter), 
76   fNtimeSamples(25),
77   fRMScut(11.)
78 {
79   //Copy constructor.
80   fNtimeSamples=phosFitter.fNtimeSamples ;
81   for(Int_t i=0; i<3;i++){
82     fLGpar[i]=phosFitter.fLGpar[i] ;
83     fHGpar[i]=phosFitter.fHGpar[i] ;
84   }
85   fRMScut=phosFitter.fRMScut ;
86 }
87
88 //-----------------------------------------------------------------------------
89 AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter)
90 {
91   //Assignment operator.
92
93   fNtimeSamples=phosFitter.fNtimeSamples ;
94   for(Int_t i=0; i<3;i++){
95     fLGpar[i]=phosFitter.fLGpar[i] ;
96     fHGpar[i]=phosFitter.fHGpar[i] ;
97   }
98   fRMScut=phosFitter.fRMScut ;
99   return *this;
100 }
101
102 //-----------------------------------------------------------------------------
103 Bool_t AliPHOSRawFitterv2::Eval()
104 {
105   //Extract an energy deposited in the crystal,
106   //crystal' position (module,column,row),
107   //time and gain (high or low).
108   //First collects sample, then evaluates it and if it has
109   //reasonable shape, fits it with Gamma2 function and extracts 
110   //energy and time.
111
112   if (fNBunches > 1) {
113     fQuality = 1000;
114     return kTRUE;
115   }
116
117   const Float_t kBaseLine   = 1.0;
118   const Int_t   kPreSamples = 10;
119
120   Float_t  pedMean   = 0;
121   Float_t  pedRMS    = 0;
122   Int_t    nPed      = 0;
123   UShort_t maxSample = 0;
124   Int_t    nMax      = 0;
125
126   TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ;
127   if(!cs)
128     cs = new TCanvas("CSample","CSample") ;
129
130   TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ;
131   if(!h) h = new TH1D("hSample","",200,0.,200.) ;
132
133   Double_t pedestal = 0;
134   for (Int_t i=0; i<fLength; i++) {
135     if (i<kPreSamples) {
136       nPed++;
137       pedMean += fSignal[i];
138       pedRMS  += fSignal[i]*fSignal[i] ;
139     }
140     if(fSignal[i] > maxSample) maxSample = fSignal[i];
141     if(fSignal[i] == maxSample) nMax++;
142     h->SetBinContent(i+1,fSignal[i]) ;
143   }
144   fEnergy = (Double_t)maxSample;
145   if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
146
147   if (fPedSubtract) {
148     if (nPed > 0) {
149       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
150       if(fPedestalRMS > 0.) 
151         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
152       fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
153     }
154     else
155       return kFALSE;
156   }
157   else {
158     //take pedestals from DB
159     pedestal = (Double_t) fAmpOffset ;
160     if (fCalibData) {
161       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
162       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
163       pedestal += truePed - altroSettings ;
164     }
165     else{
166       AliWarning(Form("Can not read data from OCDB")) ;
167     }
168     fEnergy-=pedestal ;
169   }
170
171   if (fEnergy < kBaseLine) {
172     fEnergy = 0;
173     return kTRUE;
174   }
175   
176   // calculate time
177   fTime=0. ;
178   Double_t tRMS = 0. ;
179   Double_t tW = 0. ;
180   Int_t cnts=0 ;
181   Double_t a=0,b=0,c=0 ;
182   if(fCaloFlag == 0){ // Low gain
183     a=fLGpar[0] ; 
184     b=fLGpar[1] ; 
185     c=fLGpar[2] ; 
186   }
187   else if(fCaloFlag == 1){ // High gain
188     a=fHGpar[0] ; 
189     b=fHGpar[1] ; 
190     c=fHGpar[2] ; 
191   }
192
193
194   fQuality = 0. ;
195   
196   for(Int_t i=1; i<fLength && cnts<fNtimeSamples; i++){
197     if(fSignal[i] < pedestal)
198       continue ;
199     Double_t de = fSignal[i]   - pedestal ;
200     Double_t av = fSignal[i-1] - pedestal + de;
201     if(av<=0.) //this is fluctuation around pedestal, skip it
202       continue ;
203     Double_t ds = fSignal[i] - fSignal[i-1] ;
204     Double_t ti = ds/av ;       // calculate log. derivative
205     ti = a/(ti+b)-c*ti ;        // and compare with parameterization
206     ti = i - ti ; 
207     Double_t wi = TMath::Abs(ds) ;
208     fTime += ti*wi ;
209     tW    += wi;
210     tRMS  += ti*ti*wi ;
211     cnts++ ;
212   } 
213
214   if(tW>0.){
215     fTime/=tW ;
216     fQuality = tRMS/tW-fTime*fTime ;
217   }
218   else{
219     fTime=-999. ;
220     fQuality=999. ;
221   }
222
223   Bool_t isBad = 0 ;
224   for(Int_t i=1; i<fLength-1&&!isBad; i++){
225     if(fSignal[i] > fSignal[i-1]+5 && fSignal[i] > fSignal[i+1]+5) { //single jump
226       isBad=1 ;
227     }
228   }
229   if(pedestal < 10.)
230     isBad=1 ;
231   
232   if(fPedestalRMS > 0.1)
233     isBad=1 ;
234   
235   for(Int_t i=1; i<fLength-1&&!isBad; i++){
236     if(fSignal[i] < pedestal-1)
237       isBad=1 ;
238   }
239   
240   if(fEnergy>10. && !isBad ){
241     printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ;
242     if(fOverflow)printf(" Overflow \n") ;
243     if(isBad)printf("bad") ;
244     cs->cd() ;
245     h->Draw() ;
246     cs->Update() ;
247     getchar() ;
248   }
249
250   return kTRUE;
251 }