]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PHOS/AliPHOSRawDecoderv1.cxx
Ignoring smell subdet
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawDecoderv1.cxx
CommitLineData
77ea1c6f 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 decodes the stream of ALTRO samples to extract
19// the PHOS "digits" of current event. Uses fitting procedure
20// to separate reasonable samples
21//
22// Typical use case:
23// AliRawReader* rf = new AliRawReaderDate("2006run2211.raw");
24// AliPHOSRawDecoder dc(rf);
25// while (rf->NextEvent()) {
26// dc.SetOldRCUFormat(kTRUE);
27// dc.SubtractPedestals(kTRUE);
28// while ( dc.NextDigit() ) {
29// Int_t module = dc.GetModule();
30// Int_t column = dc.GetColumn();
31// Int_t row = dc.GetRow();
32// Double_t amplitude = dc.GetEnergy();
33// Double_t time = dc.GetTime();
34// Bool_t IsLowGain = dc.IsLowGain();
35// ..........
36// }
37// }
38
39// Author: Dmitri Peressounko
40
41// --- ROOT system ---
42#include "TMath.h"
43#include "TMinuit.h"
44
45// --- AliRoot header files ---
46#include "AliPHOSRawDecoderv1.h"
47#include "AliPHOSPulseGenerator.h"
48
49
50ClassImp(AliPHOSRawDecoderv1)
51
52//-----------------------------------------------------------------------------
53 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1():AliPHOSRawDecoder()
54{
55 //Default constructor.
56}
57
58//-----------------------------------------------------------------------------
59AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(AliRawReader* rawReader, AliAltroMapping **mapping):
60 AliPHOSRawDecoder(rawReader,mapping)
61{
62 //Construct a decoder object.
63 //Is is user responsibility to provide next raw event
64 //using AliRawReader::NextEvent().
65
66 if(!gMinuit)
67 gMinuit = new TMinuit(100);
68
69}
70
71//-----------------------------------------------------------------------------
72AliPHOSRawDecoderv1::~AliPHOSRawDecoderv1()
73{
74 //Destructor.
75}
76
77//-----------------------------------------------------------------------------
78AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(const AliPHOSRawDecoderv1 &phosDecoder ):
79 AliPHOSRawDecoder(phosDecoder)
80{
81 //Copy constructor.
82}
83
84//-----------------------------------------------------------------------------
85AliPHOSRawDecoderv1& AliPHOSRawDecoderv1::operator = (const AliPHOSRawDecoderv1 &phosDecode)
86{
87 //Assignment operator.
88
89 // if(this != &phosDecode) {
90 // }
91
92 return *this;
93}
94
95//-----------------------------------------------------------------------------
96Bool_t AliPHOSRawDecoderv1::NextDigit()
97{
98 //Extract an energy deposited in the crystal,
99 //crystal' position (module,column,row),
100 //time and gain (high or low).
101 //First collects sample, then evaluates it and if it has
102 //reasonable shape, fits it with Gamma2 function and extracts
103 //energy and time.
104
105 AliCaloRawStream* in = fCaloStream;
106
107 Int_t iBin = 0;
108 Int_t tLength = 0;
109 fEnergy = -111;
110 Float_t pedMean = 0;
111 Int_t nPed = 0;
112 Float_t baseLine = 1.0;
113 const Float_t nPreSamples = 10;
114
115 while ( in->Next() ) {
116
117 if(!tLength) {
118 tLength = in->GetTimeLength();
119 if(tLength!=fSamples->GetSize()) {
120 delete fSamples ;
121 fSamples = new TArrayI(tLength);
122 }
123 else{
124 for(Int_t i=0; i<fSamples->GetSize(); i++){
125 fSamples->AddAt(0,i) ;
126 }
127 }
128 }
129
130 // Fit the full sample
131 if(in->IsNewHWAddress() && iBin>0) {
132
133 Double_t pedestal =0. ;
134 if(fPedSubtract){
135 if (nPed > 0)
136 pedestal = (Double_t)(pedMean/nPed);
137 else
138 return kFALSE;
139 }
140
141 //calculate energy
142 //first estimate if this sample looks like gamma2 function
143 Double_t aMean=0. ;
144 Double_t aRMS=0. ;
145 Int_t tStart = 0 ;
146 Int_t cnts=0 ;
147 for(Int_t i=0; i<fSamples->GetSize(); i++){
148 if(fSamples->At(i)>0){
149 Double_t de=fSamples->At(i)-pedestal ;
150 aMean+=de ;
151 aRMS+=de*de ;
152 cnts++;
153 if(de>2 && tStart==0)
154 tStart=i ;
155 if(fSamples->At(i)>fEnergy)
156 fEnergy=fSamples->At(i) ;
157 }
158 }
159 if(cnts>0){
160 aMean/=cnts;
161 aRMS=aRMS/cnts-aMean*aMean;
162 }
163
164 if(fPedSubtract){
165 fEnergy-=pedestal ;
166 }
167
168 //IF sample has reasonable mean and RMS, try to fit it with gamma2
169 if(fEnergy>2.&& cnts >20 && aMean>0. && aRMS>2.){ //more or less reasonable sample
170
171 gMinuit->mncler(); // Reset Minuit's list of paramters
172 gMinuit->SetPrintLevel(-1) ; // No Printout
173 gMinuit->SetFCN(AliPHOSRawDecoderv1::UnfoldingChiSquare) ;
174 // To set the address of the minimization function
175
176 gMinuit->SetObjectFit((TObject*)fSamples) ; // To tranfer pointer to UnfoldingChiSquare
177 Int_t ierflg ;
178 gMinuit->mnparm(0, "p", pedestal, 0.1, 0, 0, ierflg) ;
179 if(ierflg != 0){
180 // Warning("NextDigit", "Uunable to set initial value for fit procedure : ped = %f\n", pedestal ) ;
181 fEnergy=0. ;
182 fTime=-999. ;
183 return kTRUE ; //will scan further
184 }
185 gMinuit->mnparm(1, "t0", 1.*tStart, 0.1, 0, 0, ierflg) ;
186 if(ierflg != 0){
187 // Warning("NextDigit", "Unable to set initial value for fit procedure : t0\n" ) ;
188 fEnergy=0. ;
189 fTime=-999. ;
190 return kTRUE ; //will scan further
191 }
192 gMinuit->mnparm(2, "Energy", fEnergy*0.018 , 0.001*fEnergy, 0, 0, ierflg) ;
193 if(ierflg != 0){
194 // Warning("NextDigit", "Unable to set initial value for fit procedure : E=%f\n", fEnergy*0.018) ;
195 fEnergy=0. ;
196 fTime=-999. ;
197 return kTRUE ; //will scan further
198 }
199 gMinuit->mnparm(3, "Slope", 0.09 , 0.001, 0.001, 1., ierflg) ;
200 if(ierflg != 0){
201 // Warning("NextDigit", "Unable to set initial value for fit procedure : Slope\n") ;
202 fEnergy=0. ;
203 fTime=-999. ;
204 return kTRUE ; //will scan further
205 }
206
207 Double_t p0 = 1. ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
208 // depends on it.
209 Double_t p1 = 1.0 ;
210 Double_t p2 = 0.0 ;
211 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
212 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
213 // gMinuit->SetMaxIterations(100);
214 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
215
216 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
217
218 Double_t err ;
219 Double_t p,t0,efit,slope ;
220 gMinuit->GetParameter(0,p, err) ;
221 gMinuit->GetParameter(1,t0, err) ;
222 gMinuit->GetParameter(2,efit, err) ;
223 gMinuit->GetParameter(3,slope, err) ;
224
225 efit=efit*4.*TMath::Exp(-2.)/slope/slope ; //slope can not be zero - par limits
226
227 Double_t fmin,fedm,errdef ;
228 Int_t npari,nparx,istat;
229 gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
230
231 if(fmin < cnts*(efit/10.) ){ //Chi^2 of a good sample
232 if(efit<1050.&& efit>0.){
233 fEnergy=efit ;
234 fTime=t0 ;
235 }
236 }
237 else{
238 fEnergy=0 ; //bad sample
239 fTime=-999.;
240 }
241 if(fLowGainFlag)
242 fEnergy *= fPulseGenerator->GetRawFormatHighLowGainFactor(); // *16
243
244 if (fEnergy < baseLine) fEnergy = 0;
245 }
246 else{ //bad sample
247 fEnergy=0. ;
248 fTime=-999. ;
249 }
250
251 return kTRUE;
252 }
253
254 fLowGainFlag = in->IsLowGain();
255 fTime = fPulseGenerator->GetRawFormatTimeTrigger() * in->GetTime();
256 fModule = in->GetModule()+1;
257 fRow = in->GetRow() +1;
258 fColumn = in->GetColumn()+1;
259
260
261 // Fill array with samples
262 iBin++;
263 if(tLength-iBin < nPreSamples) {
264 pedMean += in->GetSignal();
265 nPed++;
266 }
267 fSamples->AddAt(in->GetSignal(),tLength-iBin);
268
269 } // in.Next()
270
271 return kFALSE;
272}
273//_____________________________________________________________________________
274void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
275{
276 // Calculates the Chi square for the samples minimization
277 // Number of parameters, Gradient, Chi squared, parameters, what to do
278
279 TArrayI * samples = (TArrayI*)gMinuit->GetObjectFit() ;
280
281 fret = 0. ;
282 if(iflag == 2)
283 for(Int_t iparam = 0 ; iparam < 4 ; iparam++)
284 Grad[iparam] = 0 ; // Will evaluate gradient
285
286 Int_t nSamples=samples->GetSize() ;
287 Double_t p =x[0] ;
288 Double_t t0=x[1] ;
289 Double_t en=x[2] ;
290 Double_t a =x[3] ;
291
292 for(Int_t i = 0 ; i < nSamples ; i++) {
293 if(samples->At(i)==0)
294 continue ;
295 Double_t dt=i*1.-t0 ;
296 Double_t diff=samples->At(i)-Gamma2(dt,p,en,a) ;
297 Double_t w=1. ; //TMath::Ceil(TMath::Abs(samples->At(i)-ped)/10.) ;
298 // if(w==0)w=1. ;
299 diff/=w ;
300 if(iflag == 2){ // calculate gradient
301 Grad[0] += -2.*diff/w ; //derivative over pedestal
302 if(dt>=0.){
303 Grad[1] += -2.*en*diff*(a*dt-2.*dt)*dt*TMath::Exp(-a*dt)/w ; //derivative over t0
304 Grad[2] += -2.*diff*dt*dt*TMath::Exp(-a*dt)/w ;
305 Grad[3] += 2.*en*diff*dt*dt*dt*TMath::Exp(-a*dt)/w ;
306 }
307 }
308 fret += diff*diff ;
309 }
310 /*
311 if(nSamples){
312 fret/=nSamples ;
313 if(iflag == 2){
314 for(Int_t iparam = 0 ; iparam < 4 ; iparam++)
315 Grad[iparam] /= nSamples ;
316 }
317 }
318 */
319
320
321}
322//-----------------------------------------------------------------------------
323Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t p,Double_t en,Double_t a){
324 //Function for fitting samples
325 //parameters:
326 //p-pedestal
327 //en-amplutude
328 //a-decay time
329
330 if(dt<0.)
331 return p ; //pedestal
332 else
333 return p+en*dt*dt*TMath::Exp(-dt*a) ;
334}
335