EMCAL
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv0.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 // A coarse algorithm takes the energy as the maximum
21 // sample, time is the first sample index, and the quality is the number
22 // of bunches in the signal.
23 // 
24 //     AliPHOSRawFitterv0 *fitterv0=new AliPHOSRawFitterv0();
25 //     fitterv0->SetChannelGeo(module,cellX,cellZ,caloFlag);
26 //     fitterv0->SetCalibData(fgCalibData) ;
27 //     fitterv0->Eval(sig,sigStart,sigLength);
28 //     Double_t amplitude = fitterv0.GetEnergy();
29 //     Double_t time      = fitterv0.GetTime();
30 //     Bool_t   isLowGain = fitterv0.GetCaloFlag()==0;
31
32 // Author: Yuri Kharlov
33
34 // --- ROOT system ---
35 #include "TArrayI.h"
36 #include "TMath.h"
37 #include "TObject.h"
38
39 // --- AliRoot header files ---
40 #include "AliPHOSRawFitterv0.h"
41 #include "AliPHOSCalibData.h"
42 #include "AliLog.h"
43
44 ClassImp(AliPHOSRawFitterv0)
45
46 //-----------------------------------------------------------------------------
47 AliPHOSRawFitterv0::AliPHOSRawFitterv0():
48   TObject(),
49   fModule(0),
50   fCellX(0),
51   fCellZ(0),
52   fCaloFlag(0),
53   fNBunches(0),
54   fPedSubtract(kFALSE),
55   fEnergy(-111),
56   fTime(-111),
57   fQuality(0.),
58   fPedestalRMS(0.),
59   fAmpOffset(0),
60   fAmpThreshold(0),
61   fOverflow(kFALSE),
62   fCalibData(0)
63 {
64   //Default constructor
65 }
66
67 //-----------------------------------------------------------------------------
68 AliPHOSRawFitterv0::~AliPHOSRawFitterv0()
69 {
70   //Destructor
71 }
72
73 //-----------------------------------------------------------------------------
74 AliPHOSRawFitterv0::AliPHOSRawFitterv0(const AliPHOSRawFitterv0 &phosFitter ):
75   TObject(),
76   fModule      (phosFitter.fModule),
77   fCellX       (phosFitter.fCellX),
78   fCellZ       (phosFitter.fCellZ),
79   fCaloFlag    (phosFitter.fCaloFlag),
80   fNBunches    (phosFitter.fNBunches),
81   fPedSubtract (phosFitter.fPedSubtract),
82   fEnergy      (phosFitter.fEnergy),
83   fTime        (phosFitter.fTime),
84   fQuality     (phosFitter.fQuality),
85   fPedestalRMS (phosFitter.fPedestalRMS),
86   fAmpOffset   (phosFitter.fAmpOffset),
87   fAmpThreshold(phosFitter.fAmpThreshold),
88   fOverflow    (phosFitter.fOverflow),
89   fCalibData   (phosFitter.fCalibData)
90 {
91   //Copy constructor
92 }
93
94 //-----------------------------------------------------------------------------
95 AliPHOSRawFitterv0& AliPHOSRawFitterv0::operator = (const AliPHOSRawFitterv0 &phosFitter)
96 {
97   //Assignment operator.
98
99   if(this != &phosFitter) {
100     fModule       = phosFitter.fModule;
101     fCellX        = phosFitter.fCellX;
102     fCellZ        = phosFitter.fCellZ;
103     fCaloFlag     = phosFitter.fCaloFlag;
104     fNBunches     = phosFitter.fNBunches;
105     fPedSubtract  = phosFitter.fPedSubtract;
106     fEnergy       = phosFitter.fEnergy;
107     fTime         = phosFitter.fTime;
108     fQuality      = phosFitter.fQuality;
109     fPedestalRMS  = phosFitter.fPedestalRMS;
110     fAmpOffset    = phosFitter.fAmpOffset;
111     fAmpThreshold = phosFitter.fAmpThreshold;
112     fOverflow     = phosFitter.fOverflow;
113     fCalibData    = phosFitter.fCalibData;
114   }
115
116   return *this;
117 }
118
119 //-----------------------------------------------------------------------------
120
121 void AliPHOSRawFitterv0::SetChannelGeo(const Int_t module, const Int_t cellX,
122                                      const Int_t cellZ,  const Int_t caloFlag)
123 {
124   // Set geometry address of the channel
125   // (for a case if fitting parameters are different for different channels)
126   
127   fModule   = module;
128   fCellX    = cellX;
129   fCellZ    = cellZ;
130   fCaloFlag = caloFlag;
131 }
132 //-----------------------------------------------------------------------------
133
134 Bool_t AliPHOSRawFitterv0::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
135 {
136   // Calculate signal parameters (energy, time, quality) from array of samples
137   // Energy is a maximum sample minus pedestal 9
138   // Time is the first time bin
139   // Signal overflows is there are at least 3 samples of the same amplitude above 900
140
141   fOverflow= kFALSE ;
142   fEnergy  = 0;
143   if (fNBunches > 1) {
144     fQuality = 1000;
145     return kTRUE;
146   }
147   
148   const Float_t kBaseLine   = 1.0;
149   const Int_t   kPreSamples = 10;
150
151   Float_t  pedMean   = 0;
152   Float_t  pedRMS    = 0;
153   Int_t    nPed      = 0;
154   UShort_t maxSample = 0;
155   Int_t    nMax      = 0;
156
157   for (Int_t i=0; i<sigLength; i++) {
158     if (i>sigLength-kPreSamples) { //inverse signal time order
159       nPed++;
160       pedMean += signal[i];
161       pedRMS  += signal[i]*signal[i] ;
162     }
163     if(signal[i] >  maxSample){ maxSample = signal[i]; nMax=0;}
164     if(signal[i] == maxSample) nMax++;
165
166   }
167   
168   
169   fEnergy = (Double_t)maxSample;
170   if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
171
172   Double_t pedestal = 0 ;
173   if (fPedSubtract) {
174     if (nPed > 0) {
175       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
176       if(fPedestalRMS > 0.) 
177         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
178       pedestal = (Double_t)(pedMean/nPed);
179     }
180     else
181       return kFALSE;
182   }
183   else {
184     //take pedestals from DB
185     pedestal = (Double_t) fAmpOffset ;
186     if (fCalibData) {
187       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
188       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
189       pedestal += truePed - altroSettings ;
190     }
191     else{
192       AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
193     }
194   }
195   fEnergy-=pedestal ;
196   if (fEnergy < kBaseLine) fEnergy = 0;
197
198   //Evaluate time
199   fTime = sigStart-sigLength-3; 
200   const Int_t nLine= 6 ;        //Parameters of fitting
201   const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
202   const Float_t kAmp=0.35 ;     //Result slightly depends on them, so no getters
203   // Avoid too low peak:
204   if(fEnergy < eMinTOF){
205      return kTRUE;
206   }
207
208   // Find index posK (kLevel is a level of "timestamp" point Tk):
209   Int_t posK =sigLength-1 ; //last point before crossing k-level
210   Double_t levelK = pedestal + kAmp*fEnergy;
211   while(signal[posK] <= levelK && posK>=0){
212      posK-- ;
213   }
214   posK++ ;
215
216   if(posK == 0 || posK==sigLength-1){
217     return kTRUE; 
218   }
219
220   // Find crosing point by solving linear equation (least squares method)
221   Int_t np = 0;
222   Int_t iup=posK-1 ;
223   Int_t idn=posK ;
224   Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
225   Double_t x,y ;
226
227   while(np<nLine){
228     //point above crossing point
229     if(iup>=0){
230       x = sigLength-iup-1;
231       y = signal[iup];
232       sx += x;
233       sy += y;
234       sxx += (x*x);
235       sxy += (x*y);
236       np++ ;
237       iup-- ;
238     }
239     //Point below crossing point
240     if(idn<sigLength){
241       if(signal[idn]<pedestal){
242         idn=sigLength-1 ; //do not scan further
243         idn++ ;
244         continue ;
245       }
246       x = sigLength-idn-1;
247       y = signal[idn];
248       sx += x;
249       sy += y;
250       sxx += (x*x);
251       sxy += (x*y);
252       np++;
253       idn++ ;
254     }
255     if(idn>=sigLength && iup<0){
256       break ; //can not fit futher
257     }
258   }
259
260   Double_t det = np*sxx - sx*sx;
261   if(det == 0){
262     return kTRUE;
263   }
264   if(np == 0){
265     return kFALSE;
266   }
267   Double_t c1 = (np*sxy - sx*sy)/det;  //slope
268   Double_t c0 = (sy-c1*sx)/np; //offset
269   if(c1 == 0){
270     return kTRUE;
271   }
272
273   // Find where the line cross kLevel:
274   fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
275   return kTRUE;
276
277 }