PHOS/AliPHOSRawFitterv2.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 extracts the signal parameters (energy, time, quality)
  19 // from ALTRO samples. Energy is in ADC counts, time is in time bin units.
  20 // Class uses FastFitting algorithm to fit sample and extract time and Amplitude
  21 // and evaluate sample quality = (chi^2/NDF)/some parameterization providing
  22 // efficiency close to 100%
  23 //
  24 // Typical use case:
  25 //     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
  26 //     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
  27 //     fitter->SetCalibData(fgCalibData) ;
  28 //     fitter->Eval(sig,sigStart,sigLength);
  29 //     Double_t amplitude = fitter.GetEnergy();
  30 //     Double_t time      = fitter.GetTime();
  31 //     Bool_t   isLowGain = fitter.GetCaloFlag()==0;
  32
  33 // Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
  34
  35 // --- ROOT system ---
  36 #include "TArrayI.h"
  37 #include "TList.h"
  38 #include "TMath.h"
  39 #include "TH1I.h"
  40 #include "TF1.h"
  41 #include "TCanvas.h"
  42 #include "TFile.h"
  43 #include "TROOT.h"
  44
  45 // --- AliRoot header files ---
  46 #include "AliLog.h"
  47 #include "AliPHOSCalibData.h"
  48 #include "AliPHOSRawFitterv2.h"
  49 #include "AliPHOSPulseGenerator.h"
  50
  51 ClassImp(AliPHOSRawFitterv2)
  52
  53 //-----------------------------------------------------------------------------
  54 AliPHOSRawFitterv2::AliPHOSRawFitterv2():
  55   AliPHOSRawFitterv0(),
  56   fAlpha(0.1),fBeta(0.035),fMax(0)
  57 {
  58   //Default constructor.
  59 }
  60
  61 //-----------------------------------------------------------------------------
  62 AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
  63 {
  64   //Destructor.
  65 }
  66
  67 //-----------------------------------------------------------------------------
  68 AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
  69   AliPHOSRawFitterv0(phosFitter),
  70   fAlpha(0.1),fBeta(0.035),fMax(0)
  71 {
  72   //Copy constructor.
  73 }
  74
  75 //-----------------------------------------------------------------------------
  76 AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 & /*phosFitter*/)
  77 {
  78   //Assignment operator.
  79   return *this;
  80 }
  81
  82 //-----------------------------------------------------------------------------
  83 Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
  84 {
  85   //Extract an energy deposited in the crystal,
  86   //crystal' position (module,column,row),
  87   //time and gain (high or low).
  88   //First collects sample, then evaluates it and if it has
  89   //reasonable shape, fits it with Gamma2 function and extracts
  90   //energy and time.
  91
  92   const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
  93   const Float_t kBaseLine   = 1.0;
  94   const Int_t   kPreSamples = 10;
  95
  96   fOverflow = kFALSE ;
  97   fEnergy=0 ;
  98   if (fCaloFlag == 2 || fNBunches > 1) {
  99     fQuality = 150;
 100     return kTRUE;
 101   }
 102   if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
 103     fQuality=200;
 104     fEnergy=0 ;
 105     return kTRUE;
 106   }
 107
 108   //Evaluate pedestals
 109   Float_t pedMean = 0;
 110   Float_t pedRMS  = 0;
 111   Int_t   nPed    = 0;
 112   for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
 113     nPed++;
 114     pedMean += signal[i];
 115     pedRMS  += signal[i]*signal[i] ;
 116   }
 117
 118   fEnergy = -111;
 119   fQuality= 999. ;
 120   Double_t pedestal = 0;
 121
 122   if (fPedSubtract) {
 123     if (nPed > 0) {
 124       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
 125       if(fPedestalRMS > 0.)
 126         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
 127       pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
 128     }
 129     else
 130       return kFALSE;
 131   }
 132   else {
 133     //take pedestals from DB
 134     pedestal = (Double_t) fAmpOffset ;
 135   }
 136
 137
 138   //calculate rough quality of the sample and check for overflow
 139   Int_t    maxAmp=0 ;
 140   Int_t    minAmp= signal[0] ;
 141   Int_t    nMax = 0 ; //number of points in plato
 142   Double_t aMean =0. ;
 143   Double_t aRMS  =0. ;
 144   Double_t wts   =0 ;
 145   Bool_t falling = kTRUE ; //Bad monotoneusly falling sample
 146   Bool_t rising = kTRUE ; //Bad monotoneusly riging sample
 147   for (Int_t i=0; i<sigLength; i++){
 148     if(signal[i] > pedestal){
 149       Double_t de = signal[i] - pedestal ;
 150       if(de > 1.) {
 151         aMean += de*i ;
 152         aRMS  += de*i*i ;
 153         wts   += de;
 154       }
 155       if(signal[i] >  maxAmp){
 156         maxAmp = signal[i];
 157         nMax=0;
 158       }
 159       if(signal[i] == maxAmp){
 160         nMax++;
 161       }
 162       if(signal[i] <  minAmp)
 163         minAmp=signal[i] ;
 164       if(falling && i>0 && signal[i]<signal[i-1])
 165         falling=kFALSE ;
 166       if(rising && i>0 && signal[i]>signal[i-1])
 167         rising=kFALSE ;
 168     }
 169   }
 170
 171   if(rising || falling){//bad "rising" or falling  sample
 172     fEnergy =    0. ;
 173     fTime   = 0. ; //-999. ;
 174     fQuality=  250. ;
 175     return kTRUE ;
 176   }
 177   if(maxAmp-minAmp<3 && maxAmp>7 && sigLength>20){ //bad flat sample
 178     fEnergy =    0. ;
 179     fTime   = 0; //-999. ;
 180     fQuality=  260. ;
 181     return kTRUE ;
 182   }
 183
 184   fEnergy=Double_t(maxAmp)-pedestal ;
 185   if (fEnergy < kBaseLine) fEnergy = 0;
 186   fTime = sigStart-sigLength-3;
 187
 188   //do not test quality of too soft samples
 189   if (wts > 0) {
 190     aMean /= wts;
 191     aRMS   = aRMS/wts - aMean*aMean;
 192   }
 193   if (fEnergy <= maxEtoFit){
 194     if (aRMS < 2.) //sigle peak
 195       fQuality = 299. ;
 196     else
 197       fQuality =   0. ;
 198     //Evaluate time of signal arriving
 199     return kTRUE ;
 200   }
 201
 202   //look for plato on the top of sample
 203   if (fEnergy>500 &&  //this is not fluctuation of soft sample
 204      nMax>2){ //and there is a plato
 205     fOverflow = kTRUE ;
 206   }
 207
 208
 209   //do not fit High Gain samples with overflow
 210   if(fCaloFlag==1 && fOverflow){
 211     fQuality = 99. ;
 212     return kTRUE;
 213
 214   }
 215
 216   //----Now fit sample with reasonable shape------
 217   TArrayD samples(sigLength); // array of sample amplitudes
 218   TArrayD times(sigLength); // array of sample time stamps
 219   for (Int_t i=0; i<sigLength; i++) {
 220     samples.AddAt(signal[i]-pedestal,sigLength-i-1);
 221     times.AddAt(double(i),i);
 222   }
 223
 224   if(fMax==0)
 225     FindMax() ;
 226   if(!FindAmpT(samples,times)){
 227     if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
 228       goto plot ;
 229     }
 230     else{
 231       return kFALSE ;
 232     }
 233   }
 234   fEnergy*=fMax ;
 235   fTime += sigStart-sigLength-3;
 236
 237
 238   //Impose cut on quality
 239 //  fQuality/=4. ;
 240   fQuality/=1.+0.005*fEnergy ;
 241
 242   //Draw corrupted samples
 243   if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
 244     if(fEnergy > 50. ){
 245     plot:
 246       printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
 247       TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
 248       if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
 249       h->Reset() ;
 250       for (Int_t i=0; i<sigLength; i++) {
 251         h->SetBinContent(i+1,float(samples.At(i))) ;
 252       }
 253 //      TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
 254       TF1 * fffit = new TF1("fffit","[0]*((x-[1])*(x-[1])*exp(-[2]*(x-[1]))+(x-[1])*exp(-[3]*(x-[1])))",0.,60.) ;
 255       fffit->SetParameters(fEnergy/fMax,fTime-(sigStart-sigLength-3),fAlpha,fBeta) ;
 256       fffit->SetLineColor(2) ;
 257       TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
 258       if(!can){
 259         can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
 260         can->SetFillColor(0) ;
 261         can->SetFillStyle(0) ;
 262         can->Range(0,0,1,1);
 263         can->SetBorderSize(0);
 264       }
 265       can->cd() ;
 266
 267       TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
 268       spectrum_1->Draw();
 269       spectrum_1->cd();
 270       spectrum_1->Range(0,0,1,1);
 271       spectrum_1->SetFillColor(0);
 272       spectrum_1->SetFillStyle(4000);
 273       spectrum_1->SetBorderSize(1);
 274       spectrum_1->SetBottomMargin(0.012);
 275       spectrum_1->SetTopMargin(0.03);
 276       spectrum_1->SetLeftMargin(0.10);
 277       spectrum_1->SetRightMargin(0.05);
 278
 279       char title[155] ;
 280       snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
 281       h->SetTitle(title) ;
 282 //      h->Fit(fffit,"","",0.,51.) ;
 283       h->Draw() ;
 284       fffit->Draw("same") ;
 285 /*
 286       sprintf(title,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
 287       TFile fout("samples_bad.root","update") ;
 288       h->Write(title);
 289       fout.Close() ;
 290 */
 291       can->cd() ;
 292       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
 293       spectrum_2->SetFillColor(0) ;
 294       spectrum_2->SetFillStyle(0) ;
 295       spectrum_2->SetGridy() ;
 296       spectrum_2->Draw();
 297       spectrum_2->Range(0,0,1,1);
 298       spectrum_2->SetFillColor(0);
 299       spectrum_2->SetBorderSize(1);
 300       spectrum_2->SetTopMargin(0.01);
 301       spectrum_2->SetBottomMargin(0.25);
 302       spectrum_2->SetLeftMargin(0.10);
 303       spectrum_2->SetRightMargin(0.05);
 304       spectrum_2->cd() ;
 305
 306       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
 307       if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
 308       hd->Reset() ;
 309       for (Int_t i=0; i<sigLength; i++) {
 310         hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
 311       }
 312       hd->Draw();
 313
 314       can->Update() ;
 315       printf("Press <enter> to continue\n") ;
 316       getchar();
 317
 318
 319       delete fffit ;
 320       delete spectrum_1 ;
 321       delete spectrum_2 ;
 322     }
 323   }
 324
 325   return kTRUE;
 326 }
 327 //------------------------------------------------------------------
 328 Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){
 329 // makes fit
 330
 331   const Int_t nMaxIter=50 ;   //Maximal number of iterations
 332   const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation
 333
 334   Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation
 335 //printf(" start fit... \n") ;
 336
 337   Int_t nPoints = samples.GetSize() ;
 338   Double_t dea=TMath::Exp(-fAlpha) ;
 339   Double_t deb=TMath::Exp(-fBeta) ;
 340   Double_t dt=1.,timeOld=dTime,dfOld=0. ;
 341   for(Int_t iter=0; iter<nMaxIter; iter++){
 342     Double_t yy=0.;
 343     Double_t yf=0. ;
 344     Double_t ydf=0. ;
 345     Double_t yddf=0. ;
 346     Double_t ff=0. ;
 347     Double_t fdf=0. ;
 348     Double_t dfdf=0. ;
 349     Double_t fddf=0. ;
 350     Int_t nfit=0 ;
 351     Double_t aexp=TMath::Exp(-fAlpha*(times.At(0)-1.-dTime)) ;
 352     Double_t bexp=TMath::Exp(-fBeta*(times.At(0)-1.-dTime)) ;
 353     for(Int_t i=0; i<nPoints; i++){
 354       Double_t t= times.At(i)-dTime ;
 355       aexp*=dea ;
 356       bexp*=deb ;
 357       if(t<0.) continue ;
 358       Double_t y=samples.At(i) ;
 359       if(y<=fAmpThreshold)
 360         continue ;
 361       nfit++ ;
 362       Double_t at=fAlpha*t ;
 363       Double_t bt = fBeta*t ;
 364       Double_t phi=t*(t*aexp+bexp) ;
 365       Double_t dphi=t*aexp*(2.-at)+(1.-bt)*bexp ;
 366       Double_t ddphi=aexp*(2.-4.*at+at*at)+bexp*fBeta*(bt-2.) ;
 367       yy+=y*y ;
 368       yf+=y*phi ;
 369       ydf+=y*dphi ;
 370       yddf+=y*ddphi ;
 371       ff+=phi*phi ;
 372       fdf+=phi*dphi ;
 373       dfdf+=dphi*dphi ;
 374       fddf+=phi*ddphi ;
 375     }
 376
 377     if(ff<1.e-09||nfit==0 ){
 378       fQuality=199 ;
 379       return kFALSE ;
 380     }
 381     Double_t f=ydf*ff-yf*fdf ;     //d(chi2)/dt
 382     Double_t df=yf*(dfdf+fddf)-yddf*ff-ydf*fdf;
 383     if(df<=0.){ //we are too far from the root. In the wicinity of root df>0
 384       if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size
 385         dt*=0.5 ;
 386         dTime=timeOld+dt ;
 387         continue ;
 388       }
 389       if(f<0){ //f<0 => dTime is too small and we still do not know root region
 390         dTime+=2. ;
 391         continue ;
 392       }
 393       else{ //dTime is too large, we are beyond the root region
 394         dTime-=2. ;
 395         continue ;
 396       }
 397     }
 398     dt=-f/df ;
 399     if(TMath::Abs(dt)<epsdt){
 400       fQuality=(yy-yf*yf/ff)/nfit ;
 401       fEnergy=yf/ff ;  //ff!=0 already tested
 402       fTime=dTime ;
 403       return kTRUE ;
 404     }
 405     //In some cases time steps are huge (derivative ~0)
 406     if(dt>10.) dt=10. ;   //restrict step size
 407     if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another
 408     timeOld=dTime ;  //remember current position for the case
 409     dfOld=df ;       //of reduction of dt step size
 410     dTime+=dt ;
 411
 412     if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy.
 413       fQuality=(yy-yf*yf/ff)/nfit ;
 414       fEnergy=yf/ff ;  //ff!=0 already tested
 415       fTime=dTime ;
 416       return kFALSE ;
 417     }
 418
 419   }
 420   //failed to find a root, too many iterations
 421   fQuality=99.;
 422   fEnergy=0 ;
 423   return kFALSE ;
 424 }
 425 //_________________________________________
 426 void AliPHOSRawFitterv2::FindMax(){
 427   //Finds maxumum of currecnt parameterization
 428   Double_t t=2./fAlpha ;
 429   fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
 430   Double_t dt=15 ;
 431   while(dt>0.01){
 432      Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ;
 433      if(dfdt>0.)
 434         t+=dt ;
 435      else
 436        t-=dt ;
 437      Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
 438      if(maxNew>fMax)
 439         fMax=maxNew ;
 440      else{
 441        dt/=2 ;
 442        if(dfdt<0.)
 443          t+=dt ;
 444        else
 445         t-=dt ;
 446      }
 447   }
 448 }
 449