Raw fitter based on FastFitting algorithm
authorprsnko <prsnko@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 9 Nov 2009 13:05:49 +0000 (13:05 +0000)
committerprsnko <prsnko@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 9 Nov 2009 13:05:49 +0000 (13:05 +0000)
PHOS/AliPHOSRawFitterv3.cxx [new file with mode: 0644]
PHOS/AliPHOSRawFitterv3.h [new file with mode: 0644]

diff --git a/PHOS/AliPHOSRawFitterv3.cxx b/PHOS/AliPHOSRawFitterv3.cxx
new file mode 100644 (file)
index 0000000..9507dd5
--- /dev/null
@@ -0,0 +1,389 @@
+/**************************************************************************
+ * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: $ */
+
+// This class extracts the signal parameters (energy, time, quality)
+// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
+// Class uses FastFitting algorithm to fit sample and extract time and Amplitude 
+// and evaluate sample quality = (chi^2/NDF)/some parameterization providing 
+// efficiency close to 100%
+// 
+// Typical use case:
+//     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
+//     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
+//     fitter->SetCalibData(fgCalibData) ;
+//     fitter->Eval(sig,sigStart,sigLength);
+//     Double_t amplitude = fitter.GetEnergy();
+//     Double_t time      = fitter.GetTime();
+//     Bool_t   isLowGain = fitter.GetCaloFlag()==0;
+
+// Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
+
+// --- ROOT system ---
+#include "TArrayI.h"
+#include "TList.h"
+#include "TMath.h"
+#include "TH1I.h"
+#include "TF1.h"
+#include "TCanvas.h"
+#include "TFile.h"
+#include "TROOT.h"
+
+// --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliPHOSCalibData.h"
+#include "AliPHOSRawFitterv3.h"
+#include "AliPHOSPulseGenerator.h"
+
+ClassImp(AliPHOSRawFitterv3)
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv3::AliPHOSRawFitterv3():
+  AliPHOSRawFitterv0()
+{
+  //Default constructor.
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
+{
+  //Destructor.
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
+  AliPHOSRawFitterv0(phosFitter) 
+{
+  //Copy constructor.
+}
+
+//-----------------------------------------------------------------------------
+AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
+{
+  //Assignment operator.
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliPHOSRawFitterv3::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
+{
+  //Extract an energy deposited in the crystal,
+  //crystal' position (module,column,row),
+  //time and gain (high or low).
+  //First collects sample, then evaluates it and if it has
+  //reasonable shape, fits it with Gamma2 function and extracts 
+  //energy and time.
+
+  if (fCaloFlag == 2 || fNBunches > 1) {
+    fQuality = 1000;
+    return kTRUE;
+  }
+  if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
+    fQuality=2000;
+    fEnergy=0 ;
+    return kTRUE;
+  }
+
+  Float_t pedMean = 0;
+  Float_t pedRMS  = 0;
+  Int_t   nPed    = 0;
+  const Float_t kBaseLine   = 1.0;
+  const Int_t   kPreSamples = 10;
+  
+
+  //We tryed individual taus for each channel, but
+  //this approach seems to be unstable. Much better results are obtaned with
+  //fixed decay time for all channels.
+  const Double_t tau=21. ;
+
+  TArrayD *samples = new TArrayD(sigLength); // array of sample amplitudes
+  TArrayD *times   = new TArrayD(sigLength); // array of sample time stamps
+  for (Int_t i=0; i<sigLength; i++) {
+    if (i<kPreSamples) {
+      nPed++;
+      pedMean += signal[i];
+      pedRMS  += signal[i]*signal[i] ;
+    }
+    samples->AddAt(signal[i],sigLength-i-1);
+    times  ->AddAt(i/tau ,i);
+  }
+
+  fEnergy = -111;
+  fQuality= 999. ;
+  const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
+  Double_t pedestal = 0;
+
+  if (fPedSubtract) {
+    if (nPed > 0) {
+      fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
+      if(fPedestalRMS > 0.) 
+       fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
+      fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
+    }
+    else
+      return kFALSE;
+  }
+  else {
+    //take pedestals from DB
+    pedestal = (Double_t) fAmpOffset ;
+    if (fCalibData) {
+      Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
+      Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
+      pedestal += truePed - altroSettings ;
+    }
+    else{
+      AliWarning(Form("Can not read data from OCDB")) ;
+    }
+    fEnergy-=pedestal ;
+  }
+
+  if (fEnergy < kBaseLine) fEnergy = 0;
+  //Evaluate time
+  Int_t iStart = 0;
+  while(iStart<sigLength && samples->At(iStart)-pedestal <kBaseLine) iStart++ ;
+  fTime = sigStart+iStart; 
+  
+  //calculate time and energy
+  Int_t    maxBin=0 ;
+  Int_t    maxAmp=0 ;
+  Int_t    nMax = 0 ; //number of points in plato
+  Double_t aMean =0. ;
+  Double_t aRMS  =0. ;
+  Double_t wts   =0 ;
+  Int_t    tStart=0 ;
+
+  for (Int_t i=0; i<sigLength; i++){
+    if(signal[i] > pedestal){
+      Double_t de = signal[i] - pedestal ;
+      if(de > 1.) {
+       aMean += de*i ;
+       aRMS  += de*i*i ;
+       wts   += de; 
+      }
+      if(de > 2 && tStart==0) 
+       tStart = i ;
+      if(signal[i] >  maxAmp){
+        maxAmp = signal[i]; 
+        nMax=0;
+       maxBin = i ;
+      }
+      if(signal[i] == maxAmp)
+        nMax++;
+    }
+  }
+
+  if (maxBin==sigLength-1){//bad "rising" sample
+    fEnergy =    0. ;
+    fTime   = -999. ;
+    fQuality=  999. ;
+    return kTRUE ;
+  }
+
+  fEnergy=Double_t(maxAmp)-pedestal ;
+  fOverflow =0 ;  //look for plato on the top of sample
+  if (fEnergy>500 &&  //this is not fluctuation of soft sample
+     nMax>2){ //and there is a plato
+    fOverflow = kTRUE ;
+  }
+  
+  if (wts > 0) {
+    aMean /= wts; 
+    aRMS   = aRMS/wts - aMean*aMean;
+  }
+
+  //do not take too small energies
+  if (fEnergy < kBaseLine) 
+    fEnergy = 0;
+  
+  //do not test quality of too soft samples
+  if (fEnergy < maxEtoFit){
+    fTime = tStart;
+    if (aRMS < 2.) //sigle peak
+      fQuality = 999. ;
+    else
+      fQuality =   0. ;
+    return kTRUE ;
+  }
+      
+  // if sample has reasonable mean and RMS, try to fit it with gamma2
+  // First estimate t0
+  Double_t a=0,b=0,c=0 ;
+  for(Int_t i=0; i<sigLength; i++){
+    if(samples->At(i)<pedestal)
+      continue ;
+    Double_t t= times->At(i) ;
+    Double_t f02 = TMath::Exp(-2.*t);
+    Double_t f12 = t*f02;
+    Double_t f22 = t*f12;
+    // Derivatives
+    Double_t f02d = -2.*f02;
+    Double_t f12d = f02 - 2.*f12;
+    Double_t f22d = 2.*(f12 - f22);
+    a += f02d * (samples->At(i)-pedestal) ;
+    b -= f12d * (samples->At(i)-pedestal) ;
+    c += f22d * (samples->At(i)-pedestal) ;
+  }
+  
+  //Find 2 roots
+  Double_t det = b*b - a*c;
+  if(det>=1.e-6 && det<0.0) {
+    det = 0.0; //rounding error
+  }
+  if(det<0.){ //Problem
+    fQuality = 1500. ;
+    if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
+      printf(" det=%e \n",det) ;
+      goto plot ;
+    }
+    return kTRUE ;
+  }
+
+  if(det>0.0 && a!=0.) {
+    det = TMath::Sqrt(det);
+    Double_t t1 = (-b + det) / a;
+//    Double_t t2 = (-b - det) / a; //second root is wrong one
+    Double_t amp1=0., den1=0. ;
+    for(Int_t i=0; i<sigLength; i++){
+       Double_t dt1 = times->At(i) - t1;
+       Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
+       amp1 += f01*(samples->At(i)-pedestal);
+       den1 += f01*f01;
+     }
+     if(den1>0.0) amp1 /= den1;
+     Double_t chi1=0.; // chi2=0. ;
+     for(Int_t i=0; i<sigLength; i++){
+       Double_t dt1 = times->At(i)- t1;
+       Double_t dy1 = samples->At(i)-pedestal- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
+       chi1 += dy1*dy1;
+     }
+     fEnergy=amp1*TMath::Exp(-2.) ; ; 
+     fTime=t1*tau ;
+     fQuality=chi1/sigLength ;
+  } 
+  else { 
+    Double_t t1 ;
+    if(a!=0.)
+      t1=-b/a ;
+    else
+      t1=-c/b ;
+    Double_t amp=0.,den=0.; ;
+    for(Int_t i=0; i<sigLength; i++){
+       Double_t dt = times->At(i) - t1;
+       Double_t f = dt*dt*TMath::Exp(-2.*dt);
+       amp += f*samples->At(i);
+       den += f*f;
+     }
+     if(den>0.0) amp /= den;
+     // chi2 calculation
+     fQuality=0.;
+     for(Int_t i=0; i<sigLength; i++){
+       Double_t t = times->At(i)- t1;
+       Double_t dy = samples->At(i)- amp*t*t*TMath::Exp(-2.*t) ;
+       fQuality += dy*dy;
+     }
+     fTime=t1*tau ;
+     fEnergy = amp*TMath::Exp(-2.);
+     fQuality/= sigLength ;
+  }
+
+  //Impose cut on quality
+  fQuality/=2.+0.004*fEnergy*fEnergy ;
+
+  //Draw corrupted samples
+  if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
+    plot:
+    if(fQuality >1. && !fOverflow ){ //Draw only bad samples
+      printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
+      TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
+      if(!h) h = new TH1I("h","Samples",50,0.,50.) ;
+      h->Reset() ;
+      for (Int_t i=0; i<sigLength; i++) {
+        h->SetBinContent(i,samples->At(i)) ;
+      }
+      TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
+      fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
+      fffit->SetLineColor(2) ;
+      TCanvas * c = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
+      if(!c){
+        c = new TCanvas("cSamples","cSamples",10,10,400,400) ;
+        c->SetFillColor(0) ;
+        c->SetFillStyle(0) ;
+        c->Range(0,0,1,1);
+        c->SetBorderSize(0);
+      }
+      c->cd() ;
+  
+      TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
+      spectrum_1->Draw();
+      spectrum_1->cd();
+      spectrum_1->Range(0,0,1,1);
+      spectrum_1->SetFillColor(0);
+      spectrum_1->SetFillStyle(4000);
+      spectrum_1->SetBorderSize(1);
+      spectrum_1->SetBottomMargin(0.012);
+      spectrum_1->SetTopMargin(0.03);
+      spectrum_1->SetLeftMargin(0.10);
+      spectrum_1->SetRightMargin(0.05);
+
+      char title[155] ;
+      sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
+      h->SetTitle(title) ;
+      h->Draw() ;
+      fffit->Draw("same") ;
+
+      sprintf(title,"mod%d_x%d_z%d_HG",fModule,fCellX,fCellZ) ;
+      TFile fout("samples_bad.root","update") ;
+      h->Write(title);
+      fout.Close() ;
+
+      c->cd() ;
+      TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
+      spectrum_2->SetFillColor(0) ;
+      spectrum_2->SetFillStyle(0) ;
+      spectrum_2->SetGridy() ;
+      spectrum_2->Draw();
+      spectrum_2->Range(0,0,1,1);
+      spectrum_2->SetFillColor(0);
+      spectrum_2->SetBorderSize(1);
+      spectrum_2->SetTopMargin(0.01);
+      spectrum_2->SetBottomMargin(0.25);
+      spectrum_2->SetLeftMargin(0.10);
+      spectrum_2->SetRightMargin(0.05);
+      spectrum_2->cd() ;
+
+      TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
+      if(!hd) hd = new TH1I("hd","Samples",50,0.,50.) ;
+      hd->Reset() ;
+      for (Int_t i=0; i<sigLength; i++) {
+        hd->SetBinContent(i,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
+      }
+      hd->Draw();
+      c->Update() ;
+      printf("Press <enter> to continue\n") ;
+      getchar();
+
+
+      delete fffit ;
+      delete spectrum_1 ;
+      delete spectrum_2 ;
+    }
+  }
+  
+  delete samples ;
+  delete times ;
+  return kTRUE;
+}
diff --git a/PHOS/AliPHOSRawFitterv3.h b/PHOS/AliPHOSRawFitterv3.h
new file mode 100644 (file)
index 0000000..ef463da
--- /dev/null
@@ -0,0 +1,31 @@
+#ifndef ALIPHOSRAWFITTERV3_H
+#define ALIPHOSRAWFITTERV3_H
+/* Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                          */
+
+/* $Id: $ */
+
+// This class extracts the PHOS "digits" of current event
+// (amplitude,time, position,gain) from the raw stream 
+// provided by AliRawReader. 
+// Uses FastFitting procedure to evaluate time and estimate sample quality
+
+#include "AliPHOSRawFitterv0.h"
+
+class AliPHOSRawFitterv3 : public AliPHOSRawFitterv0 {
+
+public:
+
+  AliPHOSRawFitterv3();
+  AliPHOSRawFitterv3(const AliPHOSRawFitterv3& rawFitter);
+  AliPHOSRawFitterv3& operator = (const AliPHOSRawFitterv3& rawFitter);
+  virtual ~AliPHOSRawFitterv3();
+
+  virtual Bool_t Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength);
+
+private:
+  
+  ClassDef(AliPHOSRawFitterv3,1)
+};
+
+#endif