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 plots samples and qualify their quality according to their shape
  19 //
  20 // Typical use case:
  21 //     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
  22 //     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
  23 //     fitter->SetCalibData(fgCalibData) ;
  24 //     fitter->Eval(sig,sigStart,sigLength);
  25 //     Double_t amplitude = fitter.GetEnergy();
  26 //     Double_t time      = fitter.GetTime();
  27 //     Bool_t   isLowGain = fitter.GetCaloFlag()==0;
  28
  29 // Author: Dmitri Peressounko (Oct.2008)
  30 // Modified: Yuri Kharlov (Jul.2009)
  31
  32 // --- ROOT system ---
  33 #include "TList.h"
  34 #include "TMath.h"
  35 #include "TMinuit.h"
  36 #include "TCanvas.h"
  37 #include "TH1.h"
  38 #include "TH2.h"
  39 #include "TF1.h"
  40 #include "TROOT.h"
  41
  42 // --- AliRoot header files ---
  43 #include "AliLog.h"
  44 #include "AliPHOSCalibData.h"
  45 #include "AliPHOSRawFitterv2.h"
  46 #include "AliPHOSPulseGenerator.h"
  47
  48 ClassImp(AliPHOSRawFitterv2)
  49
  50 //-----------------------------------------------------------------------------
  51 AliPHOSRawFitterv2::AliPHOSRawFitterv2():
  52   AliPHOSRawFitterv0(),
  53   fNtimeSamples(25),
  54   fRMScut(11.)
  55 {
  56   //Default constructor.
  57   fLGpar[0]=0.971 ;
  58   fLGpar[1]=0.0465;
  59   fLGpar[2]=1.56  ;
  60   fHGpar[0]=0.941 ;
  61   fHGpar[1]=0.0436;
  62   fHGpar[2]=1.65  ;
  63 }
  64
  65 //-----------------------------------------------------------------------------
  66 AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
  67 {
  68   //Destructor.
  69 }
  70
  71 //-----------------------------------------------------------------------------
  72 AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
  73   AliPHOSRawFitterv0(phosFitter),
  74   fNtimeSamples(25),
  75   fRMScut(11.)
  76 {
  77   //Copy constructor.
  78   fNtimeSamples=phosFitter.fNtimeSamples ;
  79   for(Int_t i=0; i<3;i++){
  80     fLGpar[i]=phosFitter.fLGpar[i] ;
  81     fHGpar[i]=phosFitter.fHGpar[i] ;
  82   }
  83   fRMScut=phosFitter.fRMScut ;
  84 }
  85
  86 //-----------------------------------------------------------------------------
  87 AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter)
  88 {
  89   //Assignment operator.
  90
  91   fNtimeSamples=phosFitter.fNtimeSamples ;
  92   for(Int_t i=0; i<3;i++){
  93     fLGpar[i]=phosFitter.fLGpar[i] ;
  94     fHGpar[i]=phosFitter.fHGpar[i] ;
  95   }
  96   fRMScut=phosFitter.fRMScut ;
  97   return *this;
  98 }
  99
 100 //-----------------------------------------------------------------------------
 101 Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
 102 {
 103   //Extract an energy deposited in the crystal,
 104   //crystal' position (module,column,row),
 105   //time and gain (high or low).
 106   //First collects sample, then evaluates it and if it has
 107   //reasonable shape, fits it with Gamma2 function and extracts
 108   //energy and time.
 109
 110   if (fNBunches > 1) {
 111     fQuality = 1000;
 112     return kTRUE;
 113   }
 114
 115   const Float_t kBaseLine   = 1.0;
 116   const Int_t   kPreSamples = 10;
 117
 118   Float_t  pedMean   = 0;
 119   Float_t  pedRMS    = 0;
 120   Int_t    nPed      = 0;
 121   UShort_t maxSample = 0;
 122   Int_t    nMax      = 0;
 123
 124   TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ;
 125   if(!cs)
 126     cs = new TCanvas("CSample","CSample") ;
 127
 128   TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ;
 129   if(!h) h = new TH1D("hSample","",200,0.,200.) ;
 130
 131   Double_t pedestal = 0;
 132   for (Int_t i=0; i<sigLength; i++) {
 133     if (i<kPreSamples) {
 134       nPed++;
 135       pedMean += signal[i];
 136       pedRMS  += signal[i]*signal[i] ;
 137     }
 138     if(signal[i] > maxSample) maxSample = signal[i];
 139     if(signal[i] == maxSample) nMax++;
 140     h->SetBinContent(i+1,signal[i]) ;
 141   }
 142   fEnergy = (Double_t)maxSample;
 143   if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
 144
 145   if (fPedSubtract) {
 146     if (nPed > 0) {
 147       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
 148       if(fPedestalRMS > 0.)
 149         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
 150       fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
 151     }
 152     else
 153       return kFALSE;
 154   }
 155   else {
 156     //take pedestals from DB
 157     pedestal = (Double_t) fAmpOffset ;
 158     if (fCalibData) {
 159       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
 160       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
 161       pedestal += truePed - altroSettings ;
 162     }
 163     else{
 164       AliWarning(Form("Can not read data from OCDB")) ;
 165     }
 166     fEnergy-=pedestal ;
 167   }
 168
 169   if (fEnergy < kBaseLine) {
 170     fEnergy = 0;
 171     return kTRUE;
 172   }
 173
 174   // calculate time
 175   fTime=0. ;
 176   Double_t tRMS = 0. ;
 177   Double_t tW = 0. ;
 178   Int_t cnts=0 ;
 179   Double_t a=0,b=0,c=0 ;
 180   if(fCaloFlag == 0){ // Low gain
 181     a=fLGpar[0] ;
 182     b=fLGpar[1] ;
 183     c=fLGpar[2] ;
 184   }
 185   else if(fCaloFlag == 1){ // High gain
 186     a=fHGpar[0] ;
 187     b=fHGpar[1] ;
 188     c=fHGpar[2] ;
 189   }
 190
 191
 192   fQuality = 0. ;
 193
 194   for(Int_t i=1; i<sigLength && cnts<fNtimeSamples; i++){
 195     if(signal[i] < pedestal)
 196       continue ;
 197     Double_t de = signal[i]   - pedestal ;
 198     Double_t av = signal[i-1] - pedestal + de;
 199     if(av<=0.) //this is fluctuation around pedestal, skip it
 200       continue ;
 201     Double_t ds = signal[i] - signal[i-1] ;
 202     Double_t ti = ds/av ;       // calculate log. derivative
 203     ti = a/(ti+b)-c*ti ;        // and compare with parameterization
 204     ti = i - ti ;
 205     Double_t wi = TMath::Abs(ds) ;
 206     fTime += ti*wi ;
 207     tW    += wi;
 208     tRMS  += ti*ti*wi ;
 209     cnts++ ;
 210   }
 211
 212   if(tW>0.){
 213     fTime/=tW ;
 214     fQuality = tRMS/tW-fTime*fTime ;
 215     fTime+=sigStart;
 216   }
 217   else{
 218     fTime=-999. ;
 219     fQuality=999. ;
 220   }
 221
 222   Bool_t isBad = 0 ;
 223   for(Int_t i=1; i<sigLength-1&&!isBad; i++){
 224     if(signal[i] > signal[i-1]+5 && signal[i] > signal[i+1]+5) { //single jump
 225       isBad=1 ;
 226     }
 227   }
 228   if(pedestal < 10.)
 229     isBad=1 ;
 230
 231   if(fPedestalRMS > 0.1)
 232     isBad=1 ;
 233
 234   for(Int_t i=1; i<sigLength-1&&!isBad; i++){
 235     if(signal[i] < pedestal-1)
 236       isBad=1 ;
 237   }
 238
 239   if(fEnergy>10. && !isBad ){
 240     printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ;
 241     if(fOverflow)printf(" Overflow \n") ;
 242     if(isBad)printf("bad") ;
 243     cs->cd() ;
 244     h->Draw() ;
 245     cs->Update() ;
 246     getchar() ;
 247   }
 248
 249   return kTRUE;
 250 }