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   fEnergy  = 0;
 142   if (fNBunches > 1) {
 143     fQuality = 1000;
 144     return kTRUE;
 145   }
 146
 147   const Float_t kBaseLine   = 1.0;
 148   const Int_t   kPreSamples = 10;
 149
 150   Float_t  pedMean   = 0;
 151   Float_t  pedRMS    = 0;
 152   Int_t    nPed      = 0;
 153   UShort_t maxSample = 0;
 154   Int_t    nMax      = 0;
 155
 156   for (Int_t i=0; i<sigLength; i++) {
 157     if (i>sigLength-kPreSamples) { //inverse signal time order
 158       nPed++;
 159       pedMean += signal[i];
 160       pedRMS  += signal[i]*signal[i] ;
 161     }
 162     if(signal[i] >  maxSample){ maxSample = signal[i]; nMax=0;}
 163     if(signal[i] == maxSample) nMax++;
 164
 165   }
 166
 167
 168   fEnergy = (Double_t)maxSample;
 169   if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
 170
 171   Double_t pedestal = 0 ;
 172   if (fPedSubtract) {
 173     if (nPed > 0) {
 174       fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
 175       if(fPedestalRMS > 0.)
 176         fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
 177       pedestal = (Double_t)(pedMean/nPed);
 178     }
 179     else
 180       return kFALSE;
 181   }
 182   else {
 183     //take pedestals from DB
 184     pedestal = (Double_t) fAmpOffset ;
 185     if (fCalibData) {
 186       Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
 187       Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
 188       pedestal += truePed - altroSettings ;
 189     }
 190     else{
 191       AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
 192     }
 193   }
 194   fEnergy-=pedestal ;
 195   if (fEnergy < kBaseLine) fEnergy = 0;
 196
 197   //Evaluate time
 198   fTime = sigStart;
 199   const Int_t nLine= 6 ;        //Parameters of fitting
 200   const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
 201   const Float_t kAmp=0.35 ;     //Result slightly depends on them, so no getters
 202   // Avoid too low peak:
 203   if(fEnergy < eMinTOF){
 204      return kTRUE;
 205   }
 206
 207   // Find index posK (kLevel is a level of "timestamp" point Tk):
 208   Int_t posK =sigLength-1 ; //last point before crossing k-level
 209   Double_t levelK = pedestal + kAmp*fEnergy;
 210   while(signal[posK] <= levelK && posK>=0){
 211      posK-- ;
 212   }
 213   posK++ ;
 214
 215   if(posK == 0 || posK==sigLength-1){
 216     return kTRUE;
 217   }
 218
 219   // Find crosing point by solving linear equation (least squares method)
 220   Int_t np = 0;
 221   Int_t iup=posK-1 ;
 222   Int_t idn=posK ;
 223   Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
 224   Double_t x,y ;
 225
 226   while(np<nLine){
 227     //point above crossing point
 228     if(iup>=0){
 229       x = sigLength-iup-1;
 230       y = signal[iup];
 231       sx += x;
 232       sy += y;
 233       sxx += (x*x);
 234       sxy += (x*y);
 235       np++ ;
 236       iup-- ;
 237     }
 238     //Point below crossing point
 239     if(idn<sigLength){
 240       if(signal[idn]<pedestal){
 241         idn=sigLength-1 ; //do not scan further
 242         idn++ ;
 243         continue ;
 244       }
 245       x = sigLength-idn-1;
 246       y = signal[idn];
 247       sx += x;
 248       sy += y;
 249       sxx += (x*x);
 250       sxy += (x*y);
 251       np++;
 252       idn++ ;
 253     }
 254     if(idn>=sigLength && iup<0){
 255       break ; //can not fit futher
 256     }
 257   }
 258
 259   Double_t det = np*sxx - sx*sx;
 260   if(det == 0){
 261     return kTRUE;
 262   }
 263   Double_t c1 = (np*sxy - sx*sy)/det;  //slope
 264   Double_t c0 = (sy-c1*sx)/np; //offset
 265   if(c1 == 0){
 266     return kTRUE;
 267   }
 268
 269   // Find where the line cross kLevel:
 270   fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
 271   return kTRUE;
 272
 273 }