PHOS/AliPHOSRawFitterv3.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 "AliPHOSRawFitterv3.h"
  49 #include "AliPHOSPulseGenerator.h"
  50
  51 ClassImp(AliPHOSRawFitterv3)
  52
  53 //-----------------------------------------------------------------------------
  54 AliPHOSRawFitterv3::AliPHOSRawFitterv3():
  55   AliPHOSRawFitterv0()
  56 {
  57   //Default constructor.
  58 }
  59
  60 //-----------------------------------------------------------------------------
  61 AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
  62 {
  63   //Destructor.
  64 }
  65
  66 //-----------------------------------------------------------------------------
  67 AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
  68   AliPHOSRawFitterv0(phosFitter)
  69 {
  70   //Copy constructor.
  71 }
  72
  73 //-----------------------------------------------------------------------------
  74 AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
  75 {
  76   //Assignment operator.
  77   return *this;
  78 }
  79
  80 //-----------------------------------------------------------------------------
  81 Bool_t AliPHOSRawFitterv3::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
  82 {
  83   //Extract an energy deposited in the crystal,
  84   //crystal' position (module,column,row),
  85   //time and gain (high or low).
  86   //First collects sample, then evaluates it and if it has
  87   //reasonable shape, fits it with Gamma2 function and extracts
  88   //energy and time.
  89
  90   if (fCaloFlag == 2 || fNBunches > 1) {
  91     fQuality = 1000;
  92     return kTRUE;
  93   }
  94   if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
  95     fQuality=2000;
  96     fEnergy=0 ;
  97     return kTRUE;
  98   }
  99
 100   Float_t pedMean = 0;
 101   Float_t pedRMS  = 0;
 102   Int_t   nPed    = 0;
 103   const Float_t kBaseLine   = 1.0;
 104   const Int_t   kPreSamples = 10;
 105
 106
 107   //We tryed individual taus for each channel, but
 108   //this approach seems to be unstable. Much better results are obtaned with
 109   //fixed decay time for all channels.
 110   const Double_t tau=22.18 ;
 111
 112   TArrayD samples(sigLength); // array of sample amplitudes
 113   TArrayD times(sigLength); // array of sample time stamps
 114   for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
 115     nPed++;
 116     pedMean += signal[i];
 117     pedRMS  += signal[i]*signal[i] ;
 118   }
 119
 120   fEnergy = -111;
 121   fQuality= 999. ;
 122   const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
 123   Double_t pedestal = 0;
 124
 125   if (fPedSubtract) {
 126     if (nPed > 0) {
 127       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
 128       if(fPedestalRMS > 0.)
 129         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
 130       pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
 131       fEnergy -= pedestal; // pedestal subtraction
 132     }
 133     else
 134       return kFALSE;
 135   }
 136   else {
 137     //take pedestals from DB
 138     pedestal = (Double_t) fAmpOffset ;
 139     if (fCalibData) {
 140       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
 141       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
 142       pedestal += truePed - altroSettings ;
 143     }
 144     else{
 145 //      AliWarning(Form("Can not read data from OCDB")) ;
 146     }
 147     fEnergy-=pedestal ;
 148   }
 149
 150   if (fEnergy < kBaseLine) fEnergy = 0;
 151   //Evaluate time
 152    fTime = sigStart-sigLength-3;
 153   for (Int_t i=0; i<sigLength; i++) {
 154     samples.AddAt(signal[i]-pedestal,sigLength-i-1);
 155     times.AddAt(i/tau ,i);
 156   }
 157
 158   //calculate time and energy
 159   Int_t    maxBin=0 ;
 160   Int_t    maxAmp=0 ;
 161   Int_t    nMax = 0 ; //number of points in plato
 162   Double_t aMean =0. ;
 163   Double_t aRMS  =0. ;
 164   Double_t wts   =0 ;
 165   for (Int_t i=0; i<sigLength; i++){
 166     if(signal[i] > pedestal){
 167       Double_t de = signal[i] - pedestal ;
 168       if(de > 1.) {
 169         aMean += de*i ;
 170         aRMS  += de*i*i ;
 171         wts   += de;
 172       }
 173       if(signal[i] >  maxAmp){
 174         maxAmp = signal[i];
 175         nMax=0;
 176         maxBin = i ;
 177       }
 178       if(signal[i] == maxAmp){
 179         nMax++;
 180       }
 181     }
 182   }
 183
 184   if (maxBin==sigLength-1){//bad "rising" sample
 185     fEnergy =    0. ;
 186     fTime   = -999. ;
 187     fQuality=  999. ;
 188     return kTRUE ;
 189   }
 190
 191   fEnergy=Double_t(maxAmp)-pedestal ;
 192   fOverflow =0 ;  //look for plato on the top of sample
 193   if (fEnergy>500 &&  //this is not fluctuation of soft sample
 194      nMax>2){ //and there is a plato
 195     fOverflow = kTRUE ;
 196   }
 197
 198   if (wts > 0) {
 199     aMean /= wts;
 200     aRMS   = aRMS/wts - aMean*aMean;
 201   }
 202
 203   //do not take too small energies
 204   if (fEnergy < kBaseLine)
 205     fEnergy = 0;
 206
 207   //do not test quality of too soft samples
 208   if (fEnergy < maxEtoFit){
 209     if (aRMS < 2.) //sigle peak
 210       fQuality = 999. ;
 211     else
 212       fQuality =   0. ;
 213     //Evaluate time of signal arriving
 214     return kTRUE ;
 215   }
 216
 217   // if sample has reasonable mean and RMS, try to fit it with gamma2
 218   //This method can not analyse overflow samples
 219   if(fOverflow){
 220     fQuality = 99. ;
 221     return kTRUE ;
 222   }
 223   // First estimate t0
 224   Double_t a=0,b=0,c=0 ;
 225   Int_t minI=0 ;
 226   if (fPedSubtract)
 227     minI=kPreSamples ;
 228   for(Int_t i=minI; i<sigLength; i++){
 229     if(samples.At(i)<=0.)
 230       continue ;
 231     Double_t t= times.At(i) ;
 232     Double_t f02 = TMath::Exp(-2.*t);
 233     Double_t f12 = t*f02;
 234     Double_t f22 = t*f12;
 235     // Derivatives
 236     Double_t f02d = -2.*f02;
 237     Double_t f12d = f02 - 2.*f12;
 238     Double_t f22d = 2.*(f12 - f22);
 239     a += f02d * samples.At(i) ;
 240     b -= f12d * samples.At(i) ;
 241     c += f22d * samples.At(i) ;
 242   }
 243
 244   //Find roots
 245   if(a==0.){
 246     if(b==0.){ //no roots
 247       fQuality = 2000. ;
 248       if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
 249         printf(" a=%f, b=%f, c=%f \n",a,b,c) ;
 250         goto plot ;
 251       }
 252       return kTRUE ;
 253     }
 254     Double_t t1=-c/b ;
 255     Double_t amp=0.,den=0.; ;
 256     for(Int_t i=minI; i<sigLength; i++){
 257       if(samples.At(i)<=0.)
 258         continue ;
 259       Double_t dt = times.At(i) - t1;
 260       Double_t f = dt*dt*TMath::Exp(-2.*dt);
 261       amp += f*samples.At(i);
 262       den += f*f;
 263     }
 264     if(den>0.0) amp /= den;
 265     // chi2 calculation
 266     fQuality=0.;
 267     for(Int_t i=minI; i<sigLength; i++){
 268       if(samples.At(i)<=0.)
 269         continue ;
 270       Double_t t = times.At(i)- t1;
 271       Double_t dy = samples.At(i)- amp*t*t*TMath::Exp(-2.*t) ;
 272       fQuality += dy*dy;
 273     }
 274     fTime+=t1*tau ;
 275     fEnergy = amp*TMath::Exp(-2.);
 276     fQuality/= sigLength ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
 277   }
 278   else{
 279     Double_t det = b*b - a*c;
 280     if(det>=1.e-6 && det<0.0) {
 281       det = 0.0; //rounding error
 282     }
 283     if(det<0.){ //Problem
 284       fQuality = 1500. ;
 285       if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
 286         printf(" det=%e \n",det) ;
 287         goto plot ;
 288       }
 289       return kTRUE ;
 290     }
 291
 292     det = TMath::Sqrt(det);
 293     Double_t t1 = (-b + det) / a;
 294 //    Double_t t2 = (-b - det) / a; //second root is wrong one
 295     Double_t amp1=0., den1=0. ;
 296     for(Int_t i=minI; i<sigLength; i++){
 297       if(samples.At(i)<=0.)
 298         continue ;
 299       Double_t dt1 = times.At(i) - t1;
 300       Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
 301       amp1 += f01*samples.At(i);
 302       den1 += f01*f01;
 303     }
 304     if(den1>0.0) amp1 /= den1;
 305     Double_t chi1=0.; // chi2=0. ;
 306     for(Int_t i=minI; i<sigLength; i++){
 307       if(samples.At(i)<=0.)
 308         continue ;
 309       Double_t dt1 = times.At(i)- t1;
 310       Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
 311       chi1 += dy1*dy1;
 312     }
 313     fEnergy=amp1*TMath::Exp(-2.) ; ;
 314     fTime+=t1*tau ;
 315     fQuality=chi1/sigLength ;
 316   }
 317
 318   //Impose cut on quality
 319   fQuality/=2.+0.004*fEnergy*fEnergy ;
 320
 321   //Draw corrupted samples
 322   if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
 323     plot:
 324     if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples
 325 //    if(!fOverflow ){ //Draw only bad samples
 326       printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
 327       TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
 328       if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
 329       h->Reset() ;
 330       for (Int_t i=0; i<sigLength; i++) {
 331         h->SetBinContent(i+1,samples.At(i)+pedestal) ;
 332       }
 333       TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
 334       fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
 335       fffit->SetLineColor(2) ;
 336       TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
 337       if(!can){
 338         can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
 339         can->SetFillColor(0) ;
 340         can->SetFillStyle(0) ;
 341         can->Range(0,0,1,1);
 342         can->SetBorderSize(0);
 343       }
 344       can->cd() ;
 345
 346       TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
 347       spectrum_1->Draw();
 348       spectrum_1->cd();
 349       spectrum_1->Range(0,0,1,1);
 350       spectrum_1->SetFillColor(0);
 351       spectrum_1->SetFillStyle(4000);
 352       spectrum_1->SetBorderSize(1);
 353       spectrum_1->SetBottomMargin(0.012);
 354       spectrum_1->SetTopMargin(0.03);
 355       spectrum_1->SetLeftMargin(0.10);
 356       spectrum_1->SetRightMargin(0.05);
 357
 358       char title[155] ;
 359       snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
 360       h->SetTitle(title) ;
 361       h->Draw() ;
 362       fffit->Draw("same") ;
 363
 364       snprintf(title,155,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
 365       TFile fout("samples_bad.root","update") ;
 366       h->Write(title);
 367       fout.Close() ;
 368
 369       can->cd() ;
 370       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
 371       spectrum_2->SetFillColor(0) ;
 372       spectrum_2->SetFillStyle(0) ;
 373       spectrum_2->SetGridy() ;
 374       spectrum_2->Draw();
 375       spectrum_2->Range(0,0,1,1);
 376       spectrum_2->SetFillColor(0);
 377       spectrum_2->SetBorderSize(1);
 378       spectrum_2->SetTopMargin(0.01);
 379       spectrum_2->SetBottomMargin(0.25);
 380       spectrum_2->SetLeftMargin(0.10);
 381       spectrum_2->SetRightMargin(0.05);
 382       spectrum_2->cd() ;
 383
 384       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
 385       if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
 386       hd->Reset() ;
 387       for (Int_t i=0; i<sigLength; i++) {
 388         hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
 389       }
 390       hd->Draw();
 391 /*
 392       can->Update() ;
 393       printf("Press <enter> to continue\n") ;
 394       getchar();
 395 */
 396
 397       delete fffit ;
 398       delete spectrum_1 ;
 399       delete spectrum_2 ;
 400     }
 401   }
 402
 403   return kTRUE;
 404 }