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