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];
 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   Int_t iStart = sigLength-1;
 199   while(iStart>=0 && signal[iStart]-pedestal <kBaseLine) iStart-- ;
 200   fTime = sigStart-iStart-2; //2: 1 due to oversubtraction in line above, another 1 due to signal started before amp increased
 201   const Int_t nLine= 6 ;        //Parameters of fitting
 202   const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
 203   const Float_t kAmp=0.35 ;     //Result slightly depends on them, so no getters
 204   // Avoid too low peak:
 205   if(fEnergy < eMinTOF){
 206      return kTRUE;
 207   }
 208
 209   // Find index posK (kLevel is a level of "timestamp" point Tk):
 210   Int_t posK =sigLength-1 ; //last point before crossing k-level
 211   Double_t levelK = pedestal + kAmp*fEnergy;
 212   while(signal[posK] <= levelK && posK>=0){
 213      posK-- ;
 214   }
 215   posK++ ;
 216
 217   if(posK == 0 || posK==sigLength-1){
 218     return kTRUE;
 219   }
 220
 221   // Find crosing point by solving linear equation (least squares method)
 222   Int_t np = 0;
 223   Int_t iup=posK-1 ;
 224   Int_t idn=posK ;
 225   Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
 226   Double_t x,y ;
 227
 228   while(np<nLine){
 229     //point above crossing point
 230     if(iup>=0){
 231       x = sigLength-iup-1;
 232       y = signal[iup];
 233       sx += x;
 234       sy += y;
 235       sxx += (x*x);
 236       sxy += (x*y);
 237       np++ ;
 238       iup-- ;
 239     }
 240     //Point below crossing point
 241     if(idn<sigLength){
 242       if(signal[idn]<pedestal){
 243         idn=sigLength-1 ; //do not scan further
 244         idn++ ;
 245         continue ;
 246       }
 247       x = sigLength-idn-1;
 248       y = signal[idn];
 249       sx += x;
 250       sy += y;
 251       sxx += (x*x);
 252       sxy += (x*y);
 253       np++;
 254       idn++ ;
 255     }
 256     if(idn>=sigLength && iup<0){
 257       break ; //can not fit futher
 258     }
 259   }
 260
 261   Double_t det = np*sxx - sx*sx;
 262   if(det == 0){
 263     return kTRUE;
 264   }
 265   Double_t c1 = (np*sxy - sx*sy)/det;  //slope
 266   Double_t c0 = (sy-c1*sx)/np; //offset
 267   if(c1 == 0){
 268     return kTRUE;
 269   }
 270
 271   // Find where the line cross kLevel:
 272   fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
 273   return kTRUE;
 274
 275 }