Possibility to set offset in SZ runs
[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       else{
219         pedestal = fAmpOffset ;
220       }
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. ;                                                                                                                   
227       Double_t wts=0 ;                                                                                                                       
228       Int_t tStart = 0 ;                                                                                                                   
229       for(Int_t i=iBin; i<fSamples->GetSize(); i++){
230         if(fSamples->At(i)>pedestal){                                                                                                             
231           Double_t de=fSamples->At(i)-pedestal ;                                                                                           
232           if(de>1.){
233             aMean+=de*i ;                                                                                                                      
234             aRMS+=de*i*i ;                                                                                                                    
235             wts+=de; 
236           }                                                                                                                         
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       }
245       if(maxBin==fSamples->GetSize()-1){//bad "rising" sample
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 ;
256       }
257       
258       if(wts>0){
259         aMean/=wts; 
260         aRMS=aRMS/wts-aMean*aMean;
261       }
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
277       
278 //Debug:=====Draw sample
279 //if(fEnergy>pedestal+10.){
280 //if(fLowGainFlag && fEnergy>2){
281 //  if(!c)
282 //    if(!fLowGainFlag && fRow==32 && fColumn==18){
283 //    TCanvas *c = new TCanvas("CSample","CSample") ;
284 //    c->cd() ;
285 //    h->Draw() ;
286 //    c->Update() ;
287 // printf("fEnergy=%f, aRMS=%f \n",fEnergy,aRMS) ;   
288 //getchar() ;
289 //}
290 //======================
291
292       //IF sample has reasonable mean and RMS, try to fit it with gamma2
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 
298         
299        fToFit->Clear("nodelete") ;
300        Double_t b,bmin,bmax ;
301        if(fLowGainFlag){
302          fSampleParamsLow->AddAt(pedestal,4) ;
303          if(fOverflow)
304            fSampleParamsLow->AddAt(double(maxAmp),5) ;
305          else
306            fSampleParamsLow->AddAt(double(1023),5) ;
307          fSampleParamsLow->AddAt(double(iBin),6) ;
308          fToFit->AddFirst((TObject*)fSampleParamsLow) ; 
309          b=fSampleParamsLow->At(2) ;
310          bmin=0.5 ;
311          bmax=10. ;
312        }
313        else{
314          fSampleParamsHigh->AddAt(pedestal,4) ;
315          if(fOverflow)
316            fSampleParamsHigh->AddAt(double(maxAmp),5) ;
317          else
318            fSampleParamsHigh->AddAt(double(1023),5);
319          fSampleParamsHigh->AddAt(double(iBin),6);
320          fToFit->AddFirst((TObject*)fSampleParamsHigh) ; 
321          b=fSampleParamsHigh->At(2) ;
322          bmin=0.05 ;
323          bmax=0.4 ;
324         }
325         fToFit->AddLast((TObject*)fSamples) ;
326         fToFit->AddLast((TObject*)fTimes) ;
327
328         gMinuit->SetObjectFit((TObject*)fToFit) ;         // To tranfer pointer to UnfoldingChiSquare
329         Int_t ierflg ;
330         gMinuit->mnparm(0, "t0",  1.*tStart, 0.01, -500., 500., ierflg) ;
331         if(ierflg != 0){
332 //        AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
333           fEnergy=0. ;
334           fTime=-999. ;
335           fQuality=999 ;
336           return kTRUE ; //will scan further
337         }
338         Double_t amp0; 
339         if(fLowGainFlag)
340           amp0=fEnergy/sampleMaxLG;
341         else
342           amp0=fEnergy/sampleMaxHG;
343
344         gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ;
345         if(ierflg != 0){
346 //        AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
347           fEnergy=0. ;
348           fTime=-999. ;
349           fQuality=999 ;
350           return kTRUE ; //will scan further
351         }
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  
362         
363         Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
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         
374         Double_t err,t0err ;
375         Double_t t0,efit ;
376         gMinuit->GetParameter(0,t0, t0err) ;    
377         gMinuit->GetParameter(1,efit, err) ;    
378
379         Double_t bfit, berr ;
380         gMinuit->GetParameter(2,bfit,berr) ;
381
382         //Calculate total energy
383         //this isparameterization of depetendence of pulse height on parameter b
384         if(fLowGainFlag)
385           efit*=99.54910 + 78.65038*bfit ;
386         else
387           efit*=80.33109+128.6433*bfit ;
388
389         if(efit<0. || efit > 10000.){                                                                          
390 //set energy to previously found max
391 //          fEnergy=0 ; //bad sample                                                    
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
403         if(fLowGainFlag){
404           fQuality/=2.+0.002*fEnergy ;
405         }
406         else{
407           fQuality/=0.75+0.0025*fEnergy ;
408         }
409
410 //Debug================
411 //        Double_t n,alpha,beta ;
412 //        Double_t en ;
413 //       if(fLowGainFlag){
414 //          n=fSampleParamsLow->At(0) ;
415 //          alpha=fSampleParamsLow->At(1) ;
416 //          beta=fSampleParamsLow->At(3) ;
417 //          en=efit/(99.54910 + 78.65038*bfit) ;
418 //        }
419 //        else{
420 //          n=fSampleParamsHigh->At(0) ;
421 //          alpha=fSampleParamsHigh->At(1) ;
422 //          beta=fSampleParamsHigh->At(3) ;
423 //          en=efit/(80.33109+128.6433*bfit) ;
424 //        }
425 //
426 ////    if( fQuality > 1 && fEnergy > 20. && !fOverflow){
427 ////    if(!fLowGainFlag && fRow==32 && fColumn==18){
428 //{
429 //    printf("Col=%d, row=%d, qual=%f, E=%f, t0=%f, b=%f\n",fColumn,fRow,fQuality,efit,t0,bfit) ;
430 //    printf("    Energy = %f \n",fEnergy) ;
431 //    TCanvas * c = new TCanvas("samp") ;
432 //    c->cd() ;
433 //    h->Draw() ;
434 //    if(fLowGainFlag){
435 //      fff->SetParameters(pedestal,en,t0,n,alpha,bfit,beta) ;
436 //    }
437 //    else{
438 //     fff->SetParameters(pedestal,en,t0,n,alpha,bfit,beta) ;
439 //    }
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.) ;
446 //    h->Draw() ;
447 //    fff->Draw("same") ;
448 //    c->Update();
449 //    getchar() ;
450 //    }
451 //====================
452
453       fEnergy=efit ;
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 ;
457       
458       return kTRUE;
459     }
460     
461     fLowGainFlag = in->IsLowGain();
462     fModule = in->GetModule()+1;
463     fRow    = in->GetRow()   +1;
464     fColumn = in->GetColumn()+1;
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();                                                                                                          
472          pedRMS += in->GetSignal()*in->GetSignal() ;
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     }
481     
482     // Fill array with samples
483     iBin--;
484     if(fPedSubtract && (in->GetTime() < nPreSamples)) {
485       pedMean += in->GetSignal();
486       pedRMS += in->GetSignal()*in->GetSignal() ;
487       nPed++;
488     }
489     fSamples->AddAt(in->GetSignal()-10,iBin);
490     fTimes->AddAt(in->GetTime(),iBin);
491  
492 //Debug==============
493 //    h->SetBinContent(in->GetTime(),in->GetSignal()) ;
494 //EndDebug==============
495     
496   } // in.Next()
497   
498   return kFALSE;
499 }
500 //_____________________________________________________________________________
501 void AliPHOSRawDecoderv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
502 {
503   // Number of parameters, Gradient, Chi squared, parameters, what to do
504
505   TList * toFit= (TList*)gMinuit->GetObjectFit() ;
506   TArrayD * params=(TArrayD*)toFit->At(0) ; 
507   TArrayI * samples = (TArrayI*)toFit->At(1) ;
508   TArrayI * times = (TArrayI*)toFit->At(2) ;
509
510   fret = 0. ;     
511   if(iflag == 2)
512     for(Int_t iparam = 0 ; iparam < 3 ; iparam++)    
513       Grad[iparam] = 0 ; // Will evaluate gradient
514   
515   Double_t t0=x[0] ;
516   Double_t en=x[1] ;
517   Double_t b=x[2] ;
518   Double_t n=params->At(0) ;
519   Double_t alpha=params->At(1) ;
520   Double_t beta=params->At(3) ;
521   Double_t ped=params->At(4) ;
522
523   Double_t overflow=params->At(5) ;
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)
526   // iBin - first non-zero sample 
527   Int_t tStep=times->At(iBin+1)-times->At(iBin) ;
528   Double_t ddt=times->At(iBin)-t0-tStep ;
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++) {
534     Double_t dt=double(times->At(i))-t0 ;
535     Double_t fsample = double(samples->At(i)) ;
536     if(fsample>=overflow)
537       continue ;
538     Double_t diff ;
539     exp1*=dexp1 ;
540     exp2*=dexp2 ;
541     if(dt<=0.){
542       diff=fsample - ped ; 
543       fret += diff*diff ;
544       continue ;
545     }
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) ;
550 //    if(fit>=overflow){
551 //      diff=fsample-overflow ;
552 //      fret += diff*diff ;
553 //      //zero gradient here
554 //    }
555 //    else{
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 ;
562       }
563 //    }
564   }
565   if(iflag == 2)
566     for(Int_t iparam = 0 ; iparam < 3 ; iparam++)    
567       Grad[iparam] *= 2. ; 
568 }
569 //-----------------------------------------------------------------------------
570 Double_t AliPHOSRawDecoderv1::Gamma2(Double_t dt,Double_t en,Double_t b,TArrayD * params){  //Function for fitting samples
571   //parameters:
572   //dt-time after start
573   //en-amplutude
574   //function parameters
575   
576   Double_t ped=params->At(4) ;
577   if(dt<0.)
578     return ped ; //pedestal
579   else
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))) ;
581 }
582