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