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