PHOS/AliPHOSRawFitterv0.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 // A coarse algorithm takes the energy as the maximum
  21 // sample, time is the first sample index, and the quality is the number
  22 // of bunches in the signal.
  23 //
  24 //     AliPHOSRawFitterv0 *fitterv0=new AliPHOSRawFitterv0();
  25 //     fitterv0->SetChannelGeo(module,cellX,cellZ,caloFlag);
  26 //     fitterv0->SetCalibData(fgCalibData) ;
  27 //     fitterv0->Eval(sig,sigStart,sigLength);
  28 //     Double_t amplitude = fitterv0.GetEnergy();
  29 //     Double_t time      = fitterv0.GetTime();
  30 //     Bool_t   isLowGain = fitterv0.GetCaloFlag()==0;
  31
  32 // Author: Yuri Kharlov
  33
  34 // --- ROOT system ---
  35 #include "TArrayI.h"
  36 #include "TMath.h"
  37 #include "TObject.h"
  38
  39 // --- AliRoot header files ---
  40 #include "AliPHOSRawFitterv0.h"
  41 #include "AliPHOSCalibData.h"
  42 #include "AliLog.h"
  43
  44 ClassImp(AliPHOSRawFitterv0)
  45
  46 //-----------------------------------------------------------------------------
  47 AliPHOSRawFitterv0::AliPHOSRawFitterv0():
  48   TObject(),
  49   fModule(0),
  50   fCellX(0),
  51   fCellZ(0),
  52   fCaloFlag(0),
  53   fNBunches(0),
  54   fPedSubtract(kFALSE),
  55   fEnergy(-111),
  56   fTime(-111),
  57   fQuality(0.),
  58   fPedestalRMS(0.),
  59   fAmpOffset(0),
  60   fAmpThreshold(0),
  61   fOverflow(kFALSE),
  62   fCalibData(0)
  63 {
  64   //Default constructor
  65 }
  66
  67 //-----------------------------------------------------------------------------
  68 AliPHOSRawFitterv0::~AliPHOSRawFitterv0()
  69 {
  70   //Destructor
  71 }
  72
  73 //-----------------------------------------------------------------------------
  74 AliPHOSRawFitterv0::AliPHOSRawFitterv0(const AliPHOSRawFitterv0 &phosFitter ):
  75   TObject(),
  76   fModule      (phosFitter.fModule),
  77   fCellX       (phosFitter.fCellX),
  78   fCellZ       (phosFitter.fCellZ),
  79   fCaloFlag    (phosFitter.fCaloFlag),
  80   fNBunches    (phosFitter.fNBunches),
  81   fPedSubtract (phosFitter.fPedSubtract),
  82   fEnergy      (phosFitter.fEnergy),
  83   fTime        (phosFitter.fTime),
  84   fQuality     (phosFitter.fQuality),
  85   fPedestalRMS (phosFitter.fPedestalRMS),
  86   fAmpOffset   (phosFitter.fAmpOffset),
  87   fAmpThreshold(phosFitter.fAmpThreshold),
  88   fOverflow    (phosFitter.fOverflow),
  89   fCalibData   (phosFitter.fCalibData)
  90 {
  91   //Copy constructor
  92 }
  93
  94 //-----------------------------------------------------------------------------
  95 AliPHOSRawFitterv0& AliPHOSRawFitterv0::operator = (const AliPHOSRawFitterv0 &phosFitter)
  96 {
  97   //Assignment operator.
  98
  99   if(this != &phosFitter) {
 100     fModule       = phosFitter.fModule;
 101     fCellX        = phosFitter.fCellX;
 102     fCellZ        = phosFitter.fCellZ;
 103     fCaloFlag     = phosFitter.fCaloFlag;
 104     fNBunches     = phosFitter.fNBunches;
 105     fPedSubtract  = phosFitter.fPedSubtract;
 106     fEnergy       = phosFitter.fEnergy;
 107     fTime         = phosFitter.fTime;
 108     fQuality      = phosFitter.fQuality;
 109     fPedestalRMS  = phosFitter.fPedestalRMS;
 110     fAmpOffset    = phosFitter.fAmpOffset;
 111     fAmpThreshold = phosFitter.fAmpThreshold;
 112     fOverflow     = phosFitter.fOverflow;
 113     fCalibData    = phosFitter.fCalibData;
 114   }
 115
 116   return *this;
 117 }
 118
 119 //-----------------------------------------------------------------------------
 120
 121 void AliPHOSRawFitterv0::SetChannelGeo(const Int_t module, const Int_t cellX,
 122                                      const Int_t cellZ,  const Int_t caloFlag)
 123 {
 124   // Set geometry address of the channel
 125   // (for a case if fitting parameters are different for different channels)
 126
 127   fModule   = module;
 128   fCellX    = cellX;
 129   fCellZ    = cellZ;
 130   fCaloFlag = caloFlag;
 131 }
 132 //-----------------------------------------------------------------------------
 133
 134 Bool_t AliPHOSRawFitterv0::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
 135 {
 136   // Calculate signal parameters (energy, time, quality) from array of samples
 137   // Energy is a maximum sample minus pedestal 9
 138   // Time is the first time bin
 139   // Signal overflows is there are at least 3 samples of the same amplitude above 900
 140
 141   fOverflow= kFALSE ;
 142   fEnergy  = 0;
 143   if (fNBunches > 1) {
 144     fQuality = 1000;
 145     return kTRUE;
 146   }
 147
 148   const Float_t kBaseLine   = 1.0;
 149   const Int_t   kPreSamples = 10;
 150
 151   Float_t  pedMean   = 0;
 152   Float_t  pedRMS    = 0;
 153   Int_t    nPed      = 0;
 154   UShort_t maxSample = 0;
 155   Int_t    nMax      = 0;
 156
 157   for (Int_t i=0; i<sigLength; i++) {
 158     if (i>sigLength-kPreSamples) { //inverse signal time order
 159       nPed++;
 160       pedMean += signal[i];
 161       pedRMS  += signal[i]*signal[i] ;
 162     }
 163     if(signal[i] >  maxSample){ maxSample = signal[i]; nMax=0;}
 164     if(signal[i] == maxSample) nMax++;
 165
 166   }
 167
 168
 169   fEnergy = (Double_t)maxSample;
 170   if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
 171
 172   Double_t pedestal = 0 ;
 173   if (fPedSubtract) {
 174     if (nPed > 0) {
 175       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
 176       if(fPedestalRMS > 0.)
 177         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
 178       pedestal = (Double_t)(pedMean/nPed);
 179     }
 180     else
 181       return kFALSE;
 182   }
 183   else {
 184     //take pedestals from DB
 185     pedestal = (Double_t) fAmpOffset ;
 186     if (fCalibData) {
 187       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
 188       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
 189       pedestal += truePed - altroSettings ;
 190     }
 191     else{
 192       AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
 193     }
 194   }
 195   fEnergy-=pedestal ;
 196   if (fEnergy < kBaseLine) fEnergy = 0;
 197
 198   //Evaluate time
 199   fTime = sigStart-sigLength-3;
 200   const Int_t nLine= 6 ;        //Parameters of fitting
 201   const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
 202   const Float_t kAmp=0.35 ;     //Result slightly depends on them, so no getters
 203   // Avoid too low peak:
 204   if(fEnergy < eMinTOF){
 205      return kTRUE;
 206   }
 207
 208   // Find index posK (kLevel is a level of "timestamp" point Tk):
 209   Int_t posK =sigLength-1 ; //last point before crossing k-level
 210   Double_t levelK = pedestal + kAmp*fEnergy;
 211   while(signal[posK] <= levelK && posK>=0){
 212      posK-- ;
 213   }
 214   posK++ ;
 215
 216   if(posK == 0 || posK==sigLength-1){
 217     return kTRUE;
 218   }
 219
 220   // Find crosing point by solving linear equation (least squares method)
 221   Int_t np = 0;
 222   Int_t iup=posK-1 ;
 223   Int_t idn=posK ;
 224   Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
 225   Double_t x,y ;
 226
 227   while(np<nLine){
 228     //point above crossing point
 229     if(iup>=0){
 230       x = sigLength-iup-1;
 231       y = signal[iup];
 232       sx += x;
 233       sy += y;
 234       sxx += (x*x);
 235       sxy += (x*y);
 236       np++ ;
 237       iup-- ;
 238     }
 239     //Point below crossing point
 240     if(idn<sigLength){
 241       if(signal[idn]<pedestal){
 242         idn=sigLength-1 ; //do not scan further
 243         idn++ ;
 244         continue ;
 245       }
 246       x = sigLength-idn-1;
 247       y = signal[idn];
 248       sx += x;
 249       sy += y;
 250       sxx += (x*x);
 251       sxy += (x*y);
 252       np++;
 253       idn++ ;
 254     }
 255     if(idn>=sigLength && iup<0){
 256       break ; //can not fit futher
 257     }
 258   }
 259
 260   Double_t det = np*sxx - sx*sx;
 261   if(det == 0){
 262     return kTRUE;
 263   }
 264   Double_t c1 = (np*sxy - sx*sy)/det;  //slope
 265   Double_t c0 = (sy-c1*sx)/np; //offset
 266   if(c1 == 0){
 267     return kTRUE;
 268   }
 269
 270   // Find where the line cross kLevel:
 271   fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
 272   return kTRUE;
 273
 274 }