+//------------------------------------------------------------------
+Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){
+// makes fit
+
+ const Int_t nMaxIter=50 ; //Maximal number of iterations
+ const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation
+
+ Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation
+//printf(" start fit... \n") ;
+
+ Int_t nPoints = samples.GetSize() ;
+ Double_t dea=TMath::Exp(-fAlpha) ;
+ Double_t deb=TMath::Exp(-fBeta) ;
+ Double_t dt=1.,timeOld=dTime,dfOld=0. ;
+ for(Int_t iter=0; iter<nMaxIter; iter++){
+ Double_t yy=0.;
+ Double_t yf=0. ;
+ Double_t ydf=0. ;
+ Double_t yddf=0. ;
+ Double_t ff=0. ;
+ Double_t fdf=0. ;
+ Double_t dfdf=0. ;
+ Double_t fddf=0. ;
+ Int_t nfit=0 ;
+ Double_t aexp=TMath::Exp(-fAlpha*(times.At(0)-1.-dTime)) ;
+ Double_t bexp=TMath::Exp(-fBeta*(times.At(0)-1.-dTime)) ;
+ for(Int_t i=0; i<nPoints; i++){
+ Double_t t= times.At(i)-dTime ;
+ aexp*=dea ;
+ bexp*=deb ;
+ if(t<0.) continue ;
+ Double_t y=samples.At(i) ;
+ if(y<=fAmpThreshold)
+ continue ;
+ nfit++ ;
+ Double_t at=fAlpha*t ;
+ Double_t bt = fBeta*t ;
+ Double_t phi=t*(t*aexp+bexp) ;
+ Double_t dphi=t*aexp*(2.-at)+(1.-bt)*bexp ;
+ Double_t ddphi=aexp*(2.-4.*at+at*at)+bexp*fBeta*(bt-2.) ;
+ yy+=y*y ;
+ yf+=y*phi ;
+ ydf+=y*dphi ;
+ yddf+=y*ddphi ;
+ ff+=phi*phi ;
+ fdf+=phi*dphi ;
+ dfdf+=dphi*dphi ;
+ fddf+=phi*ddphi ;
+ }
+
+ if(ff==0.||nfit==0. ){
+ fQuality=199 ;
+ return kFALSE ;
+ }
+ Double_t f=ydf*ff-yf*fdf ; //d(chi2)/dt
+ Double_t df=yf*(dfdf+fddf)-yddf*ff-ydf*fdf;
+ if(df<=0.){ //we are too far from the root. In the wicinity of root df>0
+ if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size
+ dt*=0.5 ;
+ dTime=timeOld+dt ;
+ continue ;
+ }
+ if(f<0){ //f<0 => dTime is too small and we still do not know root region
+ dTime+=2. ;
+ continue ;
+ }
+ else{ //dTime is too large, we are beyond the root region
+ dTime-=2. ;
+ continue ;
+ }
+ }
+ dt=-f/df ;
+ if(TMath::Abs(dt)<epsdt){
+ fQuality=(yy-yf*yf/ff)/nfit ;
+ fEnergy=yf/ff ; //ff!=0 already tested
+ fTime=dTime ;
+ return kTRUE ;
+ }
+ //In some cases time steps are huge (derivative ~0)
+ if(dt>10.) dt=10. ; //restrict step size
+ if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another
+ timeOld=dTime ; //remember current position for the case
+ dfOld=df ; //of reduction of dt step size
+ dTime+=dt ;
+
+ if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy.
+ fQuality=(yy-yf*yf/ff)/nfit ;
+ fEnergy=yf/ff ; //ff!=0 already tested
+ fTime=dTime ;
+ return kFALSE ;
+ }
+
+ }
+ //failed to find a root, too many iterations
+ fQuality=99.;
+ fEnergy=0 ;
+ return kFALSE ;
+}
+//_________________________________________
+void AliPHOSRawFitterv2::FindMax(){
+ //Finds maxumum of currecnt parameterization
+ Double_t t=2./fAlpha ;
+ fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
+ Double_t dt=15 ;
+ while(dt>0.01){
+ Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ;
+ if(dfdt>0.)
+ t+=dt ;
+ else
+ t-=dt ;
+ Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
+ if(maxNew>fMax)
+ fMax=maxNew ;
+ else{
+ dt/=2 ;
+ if(dfdt<0.)
+ t+=dt ;
+ else
+ t-=dt ;
+ }
+ }
+}
+