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