Possibility to set offset in SZ runs
[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()) {
77ea1c6f 26// dc.SubtractPedestals(kTRUE);
27// while ( dc.NextDigit() ) {
28// Int_t module = dc.GetModule();
29// Int_t column = dc.GetColumn();
30// Int_t row = dc.GetRow();
31// Double_t amplitude = dc.GetEnergy();
32// Double_t time = dc.GetTime();
33// Bool_t IsLowGain = dc.IsLowGain();
34// ..........
35// }
36// }
37
38// Author: Dmitri Peressounko
39
40// --- ROOT system ---
11f2ec15 41#include "TList.h"
77ea1c6f 42#include "TMath.h"
43#include "TMinuit.h"
44
11f2ec15 45#include "TCanvas.h"
46#include "TH1.h"
47#include "TH2.h"
48#include "TF1.h"
49#include "TROOT.h"
50
77ea1c6f 51// --- AliRoot header files ---
11f2ec15 52//#include "AliLog.h"
77ea1c6f 53#include "AliPHOSRawDecoderv1.h"
54#include "AliPHOSPulseGenerator.h"
55
56
57ClassImp(AliPHOSRawDecoderv1)
58
59//-----------------------------------------------------------------------------
11f2ec15 60 AliPHOSRawDecoderv1::AliPHOSRawDecoderv1():AliPHOSRawDecoder(),
61fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
77ea1c6f 62{
63 //Default constructor.
64}
65
66//-----------------------------------------------------------------------------
67AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(AliRawReader* rawReader, AliAltroMapping **mapping):
11f2ec15 68 AliPHOSRawDecoder(rawReader,mapping),
69 fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
77ea1c6f 70{
71 //Construct a decoder object.
72 //Is is user responsibility to provide next raw event
73 //using AliRawReader::NextEvent().
74
75 if(!gMinuit)
76 gMinuit = new TMinuit(100);
70d93620 77 fSampleParamsHigh =new TArrayD(7) ;
78 fSampleParamsHigh->AddAt(2.174,0) ;
79 fSampleParamsHigh->AddAt(0.106,1) ;
80 fSampleParamsHigh->AddAt(0.173,2) ;
81 fSampleParamsHigh->AddAt(0.06106,3) ;
82 //last two parameters are pedestal and overflow
83 fSampleParamsLow=new TArrayD(7) ;
84 fSampleParamsLow->AddAt(2.456,0) ;
85 fSampleParamsLow->AddAt(0.137,1) ;
86 fSampleParamsLow->AddAt(2.276,2) ;
87 fSampleParamsLow->AddAt(0.08246,3) ;
11f2ec15 88 fToFit = new TList() ;
77ea1c6f 89}
90
91//-----------------------------------------------------------------------------
92AliPHOSRawDecoderv1::~AliPHOSRawDecoderv1()
93{
94 //Destructor.
11f2ec15 95 if(fSampleParamsLow){
96 delete fSampleParamsLow ;
97 fSampleParamsLow=0 ;
98 }
99 if(fSampleParamsHigh){
100 delete fSampleParamsHigh ;
101 fSampleParamsHigh=0;
102 }
103 if(fToFit){
104 delete fToFit ;
105 fToFit=0 ;
106 }
77ea1c6f 107}
108
109//-----------------------------------------------------------------------------
110AliPHOSRawDecoderv1::AliPHOSRawDecoderv1(const AliPHOSRawDecoderv1 &phosDecoder ):
11f2ec15 111 AliPHOSRawDecoder(phosDecoder),
112 fSampleParamsLow(0x0),fSampleParamsHigh(0x0),fToFit(0x0)
77ea1c6f 113{
114 //Copy constructor.
11f2ec15 115 fToFit = new TList() ;
116 fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ;
117 fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ;
77ea1c6f 118}
119
120//-----------------------------------------------------------------------------
11f2ec15 121AliPHOSRawDecoderv1& AliPHOSRawDecoderv1::operator = (const AliPHOSRawDecoderv1 &phosDecoder)
77ea1c6f 122{
123 //Assignment operator.
124
11f2ec15 125 // if(this != &phosDecoder) {
77ea1c6f 126 // }
11f2ec15 127 fToFit = new TList() ;
128 if(fSampleParamsLow){
129 fSampleParamsLow = phosDecoder.fSampleParamsLow ;
130 fSampleParamsHigh= phosDecoder.fSampleParamsHigh ;
131 }
132 else{
133 fSampleParamsLow =new TArrayD(*(phosDecoder.fSampleParamsLow)) ;
134 fSampleParamsHigh=new TArrayD(*(phosDecoder.fSampleParamsHigh)) ;
135 }
77ea1c6f 136 return *this;
137}
138
139//-----------------------------------------------------------------------------
140Bool_t AliPHOSRawDecoderv1::NextDigit()
141{
142 //Extract an energy deposited in the crystal,
143 //crystal' position (module,column,row),
144 //time and gain (high or low).
145 //First collects sample, then evaluates it and if it has
146 //reasonable shape, fits it with Gamma2 function and extracts
147 //energy and time.
11f2ec15 148
70d93620 149//Debug=====================
39569d36 150// TCanvas * c = 0; //(TCanvas*)gROOT->FindObjectAny("CSample") ;
70d93620 151// if(!c)
152// c = new TCanvas("CSample","CSample") ;
153//
39569d36 154// TH1D * h = 0 ; //(TH1D*)gROOT->FindObjectAny("hSample") ;
70d93620 155// if(!h)
156// h=new TH1D("hSample","",200,0.,200.) ;
157//
39569d36 158// TF1 * fff = 0 ; //(TF1*)gROOT->FindObjectAny("fff") ;
70d93620 159// if(!fff)
160// fff = new TF1("fff","[0]+[1]*((abs(x-[2]))^[3]*exp(-(x-[2])*[4])+[5]*(x-[2])*(x-[2])*exp(-(x-[2])*[6]))",0.,1000.) ;
161//End debug===========
77ea1c6f 162
163 AliCaloRawStream* in = fCaloStream;
164
70d93620 165 Int_t iBin = fSamples->GetSize() ;
77ea1c6f 166 Int_t tLength = 0;
167 fEnergy = -111;
168 Float_t pedMean = 0;
3876e91d 169 Float_t pedRMS = 0;
77ea1c6f 170 Int_t nPed = 0;
171 Float_t baseLine = 1.0;
172 const Float_t nPreSamples = 10;
70d93620 173 fQuality= 999. ;
174 const Float_t sampleMaxHG=102.332 ; //maximal height of HG sample with given parameterization
175 const Float_t sampleMaxLG=277.196 ; //maximal height of HG sample with given parameterization
39569d36 176 const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
1defa3f3 177
77ea1c6f 178 while ( in->Next() ) {
1defa3f3 179
77ea1c6f 180 if(!tLength) {
181 tLength = in->GetTimeLength();
182 if(tLength!=fSamples->GetSize()) {
183 delete fSamples ;
70d93620 184 delete fTimes ;
77ea1c6f 185 fSamples = new TArrayI(tLength);
70d93620 186 fTimes = new TArrayI(tLength);
187 iBin= fSamples->GetSize() ;
77ea1c6f 188 }
189 else{
70d93620 190 fSamples->Reset() ;
77ea1c6f 191 }
192 }
193
194 // Fit the full sample
1defa3f3 195 if((in->IsNewHWAddress() && iBin != fSamples->GetSize()) //new HW address
196 ||(iBin<=0)) { //or new signal in same address
39569d36 197
198 //First remember new sample
199 fNewLowGainFlag = in->IsLowGain();
200 fNewModule = in->GetModule()+1;
201 fNewRow = in->GetRow() +1;
202 fNewColumn = in->GetColumn()+1;
203 fNewAmp = in->GetSignal() ;
1defa3f3 204 fNewTime=in->GetTime() ;
205
3876e91d 206 //now handle already collected
77ea1c6f 207 Double_t pedestal =0. ;
3876e91d 208 fPedestalRMS=0. ;
77ea1c6f 209 if(fPedSubtract){
3876e91d 210 if (nPed > 0){
77ea1c6f 211 pedestal = (Double_t)(pedMean/nPed);
3876e91d 212 fPedestalRMS=pedRMS/nPed-pedestal*pedestal ;
213 if(fPedestalRMS>0.) fPedestalRMS=TMath::Sqrt(fPedestalRMS) ;
214 }
77ea1c6f 215 else
216 return kFALSE;
217 }
60144e5f 218 else{
219 pedestal = fAmpOffset ;
220 }
70d93620 221
222 //calculate time and energy
223 Int_t maxBin=0 ;
224 Int_t maxAmp=0 ;
225 Double_t aMean=0. ;
226 Double_t aRMS=0. ;
39569d36 227 Double_t wts=0 ;
70d93620 228 Int_t tStart = 0 ;
70d93620 229 for(Int_t i=iBin; i<fSamples->GetSize(); i++){
60144e5f 230 if(fSamples->At(i)>pedestal){
70d93620 231 Double_t de=fSamples->At(i)-pedestal ;
39569d36 232 if(de>1.){
233 aMean+=de*i ;
234 aRMS+=de*i*i ;
235 wts+=de;
236 }
70d93620 237 if(de>2 && tStart==0)
238 tStart=i ;
239 if(maxAmp<fSamples->At(i)){
240 maxBin=i ;
241 maxAmp=fSamples->At(i) ;
242 }
243 }
244 }
39569d36 245 if(maxBin==fSamples->GetSize()-1){//bad "rising" sample
70d93620 246 fEnergy=0. ;
247 fTime=-999.;
248 fQuality= 999. ;
249 return kTRUE ;
250 }
251 fEnergy=Double_t(maxAmp)-pedestal ;
252 fOverflow =0 ; //look for plato on the top of sample
253 if(fEnergy>500 && //this is not fluctuation of soft sample
254 maxBin<fSamples->GetSize()-1 && fSamples->At(maxBin+1)==maxAmp){ //and there is a plato
255 fOverflow = kTRUE ;
77ea1c6f 256 }
70d93620 257
39569d36 258 if(wts>0){
259 aMean/=wts;
260 aRMS=aRMS/wts-aMean*aMean;
77ea1c6f 261 }
39569d36 262
263 //do not take too small energies
264 if(fEnergy < baseLine)
265 fEnergy = 0;
266
267 //do not test quality of too soft samples
268 if(fEnergy<maxEtoFit){
269 fTime=fTimes->At(tStart);
270 if(aRMS<2.) //sigle peak
271 fQuality=999. ;
272 else
273 fQuality= 0. ;
274 return kTRUE ;
275 }
276
77ea1c6f 277
70d93620 278//Debug:=====Draw sample
279//if(fEnergy>pedestal+10.){
39569d36 280//if(fLowGainFlag && fEnergy>2){
281// if(!c)
113139da 282// if(!fLowGainFlag && fRow==32 && fColumn==18){
283// TCanvas *c = new TCanvas("CSample","CSample") ;
70d93620 284// c->cd() ;
285// h->Draw() ;
286// c->Update() ;
39569d36 287// printf("fEnergy=%f, aRMS=%f \n",fEnergy,aRMS) ;
288//getchar() ;
70d93620 289//}
290//======================
291
77ea1c6f 292 //IF sample has reasonable mean and RMS, try to fit it with gamma2
77ea1c6f 293
294 gMinuit->mncler(); // Reset Minuit's list of paramters
295 gMinuit->SetPrintLevel(-1) ; // No Printout
296 gMinuit->SetFCN(AliPHOSRawDecoderv1::UnfoldingChiSquare) ;
297 // To set the address of the minimization function
11f2ec15 298
60144e5f 299 fToFit->Clear("nodelete") ;
113139da 300 Double_t b,bmin,bmax ;
11f2ec15 301 if(fLowGainFlag){
302 fSampleParamsLow->AddAt(pedestal,4) ;
70d93620 303 if(fOverflow)
304 fSampleParamsLow->AddAt(double(maxAmp),5) ;
305 else
306 fSampleParamsLow->AddAt(double(1023),5) ;
307 fSampleParamsLow->AddAt(double(iBin),6) ;
11f2ec15 308 fToFit->AddFirst((TObject*)fSampleParamsLow) ;
113139da 309 b=fSampleParamsLow->At(2) ;
310 bmin=0.5 ;
311 bmax=10. ;
312 }
313 else{
11f2ec15 314 fSampleParamsHigh->AddAt(pedestal,4) ;
70d93620 315 if(fOverflow)
316 fSampleParamsHigh->AddAt(double(maxAmp),5) ;
317 else
318 fSampleParamsHigh->AddAt(double(1023),5);
319 fSampleParamsHigh->AddAt(double(iBin),6);
11f2ec15 320 fToFit->AddFirst((TObject*)fSampleParamsHigh) ;
113139da 321 b=fSampleParamsHigh->At(2) ;
322 bmin=0.05 ;
323 bmax=0.4 ;
11f2ec15 324 }
325 fToFit->AddLast((TObject*)fSamples) ;
70d93620 326 fToFit->AddLast((TObject*)fTimes) ;
11f2ec15 327
328 gMinuit->SetObjectFit((TObject*)fToFit) ; // To tranfer pointer to UnfoldingChiSquare
77ea1c6f 329 Int_t ierflg ;
1defa3f3 330 gMinuit->mnparm(0, "t0", 1.*tStart, 0.01, -500., 500., ierflg) ;
77ea1c6f 331 if(ierflg != 0){
11f2ec15 332// AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
77ea1c6f 333 fEnergy=0. ;
334 fTime=-999. ;
70d93620 335 fQuality=999 ;
77ea1c6f 336 return kTRUE ; //will scan further
337 }
70d93620 338 Double_t amp0;
339 if(fLowGainFlag)
340 amp0=fEnergy/sampleMaxLG;
341 else
342 amp0=fEnergy/sampleMaxHG;
11f2ec15 343
113139da 344 gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ;
77ea1c6f 345 if(ierflg != 0){
11f2ec15 346// AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
77ea1c6f 347 fEnergy=0. ;
348 fTime=-999. ;
70d93620 349 fQuality=999 ;
77ea1c6f 350 return kTRUE ; //will scan further
351 }
113139da 352
353 gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ;
354 if(ierflg != 0){
355// AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
356 fEnergy=0. ;
357 fTime=-999. ;
358 fQuality=999 ;
359 return kTRUE ; //will scan further
360 }
361
77ea1c6f 362
11f2ec15 363 Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
77ea1c6f 364 // depends on it.
365 Double_t p1 = 1.0 ;
366 Double_t p2 = 0.0 ;
367 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
368 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
369 // gMinuit->SetMaxIterations(100);
370 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
371
372 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
373
11f2ec15 374 Double_t err,t0err ;
375 Double_t t0,efit ;
376 gMinuit->GetParameter(0,t0, t0err) ;
377 gMinuit->GetParameter(1,efit, err) ;
11f2ec15 378
113139da 379 Double_t bfit, berr ;
380 gMinuit->GetParameter(2,bfit,berr) ;
381
70d93620 382 //Calculate total energy
113139da 383 //this isparameterization of depetendence of pulse height on parameter b
70d93620 384 if(fLowGainFlag)
113139da 385 efit*=99.54910 + 78.65038*bfit ;
70d93620 386 else
113139da 387 efit*=80.33109+128.6433*bfit ;
70d93620 388
389 if(efit<0. || efit > 10000.){
e6fe1709 390//set energy to previously found max
b20d9db0 391// fEnergy=0 ; //bad sample
70d93620 392 fTime=-999.;
393 fQuality=999 ;
394 return kTRUE;
395 }
396
397 //evaluate fit quality
398 Double_t fmin,fedm,errdef ;
399 Int_t npari,nparx,istat;
400 gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
401 fQuality=fmin/(fSamples->GetSize()-iBin) ;
402 //compare quality with some parameterization
113139da 403 if(fLowGainFlag){
404 fQuality/=2.+0.002*fEnergy ;
405 }
406 else{
407 fQuality/=0.75+0.0025*fEnergy ;
408 }
70d93620 409
410//Debug================
113139da 411// Double_t n,alpha,beta ;
e6fe1709 412// Double_t en ;
413// if(fLowGainFlag){
70d93620 414// n=fSampleParamsLow->At(0) ;
415// alpha=fSampleParamsLow->At(1) ;
70d93620 416// beta=fSampleParamsLow->At(3) ;
e6fe1709 417// en=efit/(99.54910 + 78.65038*bfit) ;
70d93620 418// }
419// else{
420// n=fSampleParamsHigh->At(0) ;
421// alpha=fSampleParamsHigh->At(1) ;
70d93620 422// beta=fSampleParamsHigh->At(3) ;
e6fe1709 423// en=efit/(80.33109+128.6433*bfit) ;
70d93620 424// }
425//
e6fe1709 426//// if( fQuality > 1 && fEnergy > 20. && !fOverflow){
427//// if(!fLowGainFlag && fRow==32 && fColumn==18){
428//{
113139da 429// printf("Col=%d, row=%d, qual=%f, E=%f, t0=%f, b=%f\n",fColumn,fRow,fQuality,efit,t0,bfit) ;
e6fe1709 430// printf(" Energy = %f \n",fEnergy) ;
113139da 431// TCanvas * c = new TCanvas("samp") ;
11f2ec15 432// c->cd() ;
433// h->Draw() ;
434// if(fLowGainFlag){
e6fe1709 435// fff->SetParameters(pedestal,en,t0,n,alpha,bfit,beta) ;
11f2ec15 436// }
437// else{
e6fe1709 438// fff->SetParameters(pedestal,en,t0,n,alpha,bfit,beta) ;
113139da 439// }
e6fe1709 440////// for(Int_t i=1;i<=h->GetNbinsX(); i++){
441//// Double_t x=h->GetBinCenter(i) ;
442//// h->SetBinContent(i,h->GetBinContent(i)-fff->Eval(x)) ;
443//// }
444//// h->SetMinimum(-15.) ;
445//// h->SetMaximum(15.) ;
113139da 446// h->Draw() ;
11f2ec15 447// fff->Draw("same") ;
448// c->Update();
70d93620 449// getchar() ;
450// }
451//====================
11f2ec15 452
39569d36 453 fEnergy=efit ;
1defa3f3 454 fTime=t0-4.024*bfit ; //-10.402*bfit+4.669*bfit*bfit ; //Correction for 70 samples
455// fTime=t0+2.8*bfit ; //-10.402*bfit+4.669*bfit*bfit ; //Correction for 50 samples
456// fQuality = bfit ;
77ea1c6f 457
458 return kTRUE;
459 }
460
461 fLowGainFlag = in->IsLowGain();
77ea1c6f 462 fModule = in->GetModule()+1;
463 fRow = in->GetRow() +1;
464 fColumn = in->GetColumn()+1;
39569d36 465
466 //add previouly taken if coincides
467 if(fLowGainFlag==fNewLowGainFlag && fModule==fNewModule &&
468 fRow==fNewRow && fColumn==fNewColumn){
469 iBin--;
470 if(fPedSubtract && fNewTime < nPreSamples) {
471 pedMean += in->GetSignal();
3876e91d 472 pedRMS += in->GetSignal()*in->GetSignal() ;
39569d36 473 nPed++;
474 }
475 fSamples->AddAt(fNewAmp,iBin);
476 fTimes->AddAt(fNewTime,iBin);
477
478 //Mark that we already take it
479 fNewModule=-1 ;
480 }
77ea1c6f 481
77ea1c6f 482 // Fill array with samples
70d93620 483 iBin--;
484 if(fPedSubtract && (in->GetTime() < nPreSamples)) {
77ea1c6f 485 pedMean += in->GetSignal();
3876e91d 486 pedRMS += in->GetSignal()*in->GetSignal() ;
77ea1c6f 487 nPed++;
488 }
60144e5f 489 fSamples->AddAt(in->GetSignal()-10,iBin);
70d93620 490 fTimes->AddAt(in->GetTime(),iBin);
491
492//Debug==============
e6fe1709 493// h->SetBinContent(in->GetTime(),in->GetSignal()) ;
494//EndDebug==============
77ea1c6f 495
496 } // in.Next()
497
498 return kFALSE;
499}
500//_____________________________________________________________________________
501void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
502{
77ea1c6f 503 // Number of parameters, Gradient, Chi squared, parameters, what to do
504
11f2ec15 505 TList * toFit= (TList*)gMinuit->GetObjectFit() ;
506 TArrayD * params=(TArrayD*)toFit->At(0) ;
507 TArrayI * samples = (TArrayI*)toFit->At(1) ;
70d93620 508 TArrayI * times = (TArrayI*)toFit->At(2) ;
77ea1c6f 509
510 fret = 0. ;
511 if(iflag == 2)
113139da 512 for(Int_t iparam = 0 ; iparam < 3 ; iparam++)
77ea1c6f 513 Grad[iparam] = 0 ; // Will evaluate gradient
514
11f2ec15 515 Double_t t0=x[0] ;
516 Double_t en=x[1] ;
113139da 517 Double_t b=x[2] ;
70d93620 518 Double_t n=params->At(0) ;
11f2ec15 519 Double_t alpha=params->At(1) ;
11f2ec15 520 Double_t beta=params->At(3) ;
113139da 521 Double_t ped=params->At(4) ;
522
523 Double_t overflow=params->At(5) ;
70d93620 524 Int_t iBin = (Int_t) params->At(6) ;
525 Int_t nSamples=TMath::Min(iBin+70,samples->GetSize()) ; //Here we set number of points to fit (70)
e6fe1709 526 // iBin - first non-zero sample
113139da 527 Int_t tStep=times->At(iBin+1)-times->At(iBin) ;
e6fe1709 528 Double_t ddt=times->At(iBin)-t0-tStep ;
113139da 529 Double_t exp1=TMath::Exp(-alpha*ddt) ;
530 Double_t exp2=TMath::Exp(-beta*ddt) ;
531 Double_t dexp1=TMath::Exp(-alpha*tStep) ;
532 Double_t dexp2=TMath::Exp(-beta*tStep) ;
533 for(Int_t i = iBin; i<nSamples ; i++) {
113139da 534 Double_t dt=double(times->At(i))-t0 ;
e6fe1709 535 Double_t fsample = double(samples->At(i)) ;
1defa3f3 536 if(fsample>=overflow)
537 continue ;
113139da 538 Double_t diff ;
e6fe1709 539 exp1*=dexp1 ;
540 exp2*=dexp2 ;
1defa3f3 541 if(dt<=0.){
113139da 542 diff=fsample - ped ;
543 fret += diff*diff ;
77ea1c6f 544 continue ;
113139da 545 }
113139da 546 Double_t dtn=TMath::Power(dt,n) ;
547 Double_t dtnE=dtn*exp1 ;
548 Double_t dt2E=dt*dt*exp2 ;
549 Double_t fit=ped+en*(dtnE + b*dt2E) ;
1defa3f3 550// if(fit>=overflow){
551// diff=fsample-overflow ;
552// fret += diff*diff ;
553// //zero gradient here
554// }
555// else{
113139da 556 diff = fsample - fit ;
557 fret += diff*diff ;
558 if(iflag == 2){ // calculate gradient
559 Grad[0] += en*diff*(dtnE*(n/dt-alpha)+b*dt2E*(2./dt-beta)) ; //derivative over t0
560 Grad[1] -= diff*(dtnE+b*dt2E) ;
561 Grad[2] -= en*diff*dt2E ;
77ea1c6f 562 }
1defa3f3 563// }
77ea1c6f 564 }
113139da 565 if(iflag == 2)
566 for(Int_t iparam = 0 ; iparam < 3 ; iparam++)
567 Grad[iparam] *= 2. ;
77ea1c6f 568}
569//-----------------------------------------------------------------------------
113139da 570Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t en,Double_t b,TArrayD * params){ //Function for fitting samples
77ea1c6f 571 //parameters:
11f2ec15 572 //dt-time after start
77ea1c6f 573 //en-amplutude
11f2ec15 574 //function parameters
575
576 Double_t ped=params->At(4) ;
77ea1c6f 577 if(dt<0.)
11f2ec15 578 return ped ; //pedestal
77ea1c6f 579 else
113139da 580 return ped+en*(TMath::Power(dt,params->At(0))*TMath::Exp(-dt*params->At(1))+b*dt*dt*TMath::Exp(-dt*params->At(3))) ;
77ea1c6f 581}
582