STEER/STEERBase/AliEMCALPIDResponse.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //                                                                      //
  17 // AliEMCALPIDResponse                                                  //
  18 //                                                                      //
  19 // EMCAL class to perfom PID                                            //
  20 // This is a prototype and still under development                      //
  21 //                                                                      //
  22 // ---------------------------------------------------------------------//
  23 // GetNumberOfSigmas():                                                 //
  24 //                                                                      //
  25 // Electrons:  Number of Sigmas for E/p value                           //
  26 //             Parametrization of LHC11a (after recalibration)          //
  27 //                                                                      //
  28 // NON electrons:                                                       //
  29 //             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
  30 //             --> return +/- 99                                        //
  31 //             Otherwise                                                //
  32 //             --> return nsigma (parametrization of LHC10e)            //
  33 //                                                                      //
  34 // NO Parametrization (outside pT range): --> return -999               //
  35 //                                                                      //
  36 // ---------------------------------------------------------------------//
  37 // ComputeEMCALProbability():                                           //
  38 //                                                                      //
  39 // Electrons:  Probability from Gaussian distribution                   //
  40 //                                                                      //
  41 // NON electrons:                                                       //
  42 //             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
  43 //             --> probability to find particles below or above thr.    //
  44 //             Otherwise                                                //
  45 //             -->  Probability from Gaussian distribution              //
  46 //                  (proper normalization to each other?)               //
  47 //                                                                      //
  48 // NO Parametrization (outside pT range): --> return kFALSE             //
  49 //////////////////////////////////////////////////////////////////////////
  50
  51 #include <TF1.h>
  52 #include <TMath.h>
  53
  54 #include "AliEMCALPIDResponse.h"       //class header
  55
  56 #include "AliLog.h"
  57
  58 ClassImp(AliEMCALPIDResponse)
  59
  60 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  61 AliEMCALPIDResponse::AliEMCALPIDResponse():
  62   TObject(),
  63   fNorm(NULL),
  64   fCurrCentrality(-1.),
  65   fkPIDParams(NULL)
  66 {
  67   //
  68   //  The default constructor
  69   //
  70
  71
  72   fNorm = new TF1("fNorm","gaus",-20,20);
  73 }
  74 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  75 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
  76   TObject(other),
  77   fNorm(other.fNorm),
  78   fCurrCentrality(other.fCurrCentrality),
  79   fkPIDParams(other.fkPIDParams)
  80 {
  81   //
  82   //  The copy constructor
  83   //
  84
  85 }
  86 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  87 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
  88 {
  89   //
  90   //  The assignment operator
  91   //
  92
  93   if(this == &other) return *this;
  94
  95   // Make copy
  96   TObject::operator=(other);
  97   fNorm = other.fNorm;
  98   fCurrCentrality = other.fCurrCentrality;
  99   fkPIDParams = other.fkPIDParams;
 100
 101
 102   return *this;
 103 }
 104 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 105 AliEMCALPIDResponse::~AliEMCALPIDResponse() {
 106
 107   delete fNorm;
 108
 109 }
 110 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 111 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {
 112   //
 113   // Calculates the expected sigma of the PID signal as the function of
 114   // the information stored in the track, for the specified particle type
 115   //
 116   //
 117
 118   Double_t norm = 1.;
 119
 120   // Check the charge
 121   if( charge != -1 && charge != 1){
 122     return norm;
 123   }
 124
 125   // Get the parameters for this particle type and pt
 126   const TVectorD *params = GetParams(n, pt, charge);
 127
 128   // IF not in momentum range, NULL is returned --> return default value
 129   if(!params) return norm;
 130
 131   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
 132   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
 133   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
 134   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
 135   Double_t probLow  = (*params)[6];   // probability to be below eopMin
 136   Double_t probHigh = (*params)[7];   // probability to be above eopMax
 137
 138   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
 139   fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
 140   norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
 141
 142   return norm;
 143 }
 144 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 145 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {
 146
 147   Double_t nsigma = -999.;
 148
 149   // Check the charge
 150   if( charge != -1 && charge != 1){
 151     return nsigma;
 152   }
 153
 154   // Get the parameters for this particle type and pt
 155   const TVectorD *params = GetParams(n, pt, charge);
 156
 157   // IF not in momentum range, NULL is returned --> return default value
 158   if(!params) return nsigma;
 159
 160   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
 161   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
 162   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
 163   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
 164
 165   // if electron
 166   if(n == AliPID::kElectron){
 167     if(sigma != 0) nsigma = (eop - mean) / sigma;
 168   }
 169
 170   // if NON electron
 171   else{
 172     if ( eop < eopMin )
 173       nsigma = -99;    // not parametrized (smaller than eopMin)
 174     else if ( eop > eopMax )
 175       nsigma = 99.;     // not parametrized (bigger than eopMax)
 176     else{
 177       if(sigma != 0) nsigma = (eop - mean) / sigma;
 178     }
 179   }
 180
 181   return nsigma;
 182
 183 }
 184 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 185 Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
 186   //
 187   //
 188   Double_t fRange  = 5.0;   // hardcoded (???)
 189   Double_t nsigma  = 0.0;
 190
 191
 192   // Check the charge
 193   if( charge != -1 && charge != 1){
 194     return kFALSE;
 195   }
 196
 197
 198   // default value (will be returned, if pt below threshold)
 199   for (Int_t species = 0; species < nSpecies; species++) {
 200     pEMCAL[species] = 1./nSpecies;
 201   }
 202
 203   // set E/p range
 204   if(eop < 0.05) eop = 0.05;
 205   if(eop > 2.00) eop = 2.00;
 206
 207   for (Int_t species = 0; species < nSpecies; species++) {
 208
 209     AliPID::EParticleType type = AliPID::EParticleType(species);
 210
 211     // Get the parameters for this particle type and pt
 212     const TVectorD *params = GetParams(species, pt, charge);
 213
 214     // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
 215     if(!params) return kFALSE;
 216
 217     Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
 218     Double_t probLow  = (*params)[6];   // probability to be below eopMin
 219     Double_t probHigh = (*params)[7];   // probability to be above eopMax
 220
 221     // get nsigma value for each particle type at this E/p value
 222     nsigma = GetNumberOfSigmas(pt,eop,type,charge);
 223
 224     // electrons (standard Gaussian calculation of probabilities)
 225     if(type == AliPID::kElectron){
 226       if (TMath::Abs(nsigma) > fRange) {
 227         pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
 228       }
 229       else{
 230         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
 231       }
 232     }
 233     //NON electrons
 234     else{
 235       // E/p < eopMin  -->  return probability below E/p = eopMin
 236       if ( nsigma == -99){
 237         pEMCAL[species] = probLow;
 238       }
 239       // E/p > eopMax  -->  return probability above E/p = eopMax
 240       else if ( nsigma == 99){
 241         pEMCAL[species] = probHigh;
 242       }
 243       // in parametrized region --> calculate probability for corresponding Gauss curve
 244       else{
 245         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
 246
 247         // normalize to total probability == 1
 248         pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
 249       }
 250     }
 251   }
 252
 253   return kTRUE;
 254
 255 }
 256 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 257 const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
 258   //
 259   // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
 260   //
 261   // 0 = momMin
 262   // 1 = momMax
 263   // 2 = mean of Gaus
 264   // 3 = sigma of Gaus
 265   // 4 = eopLow
 266   // 5 = eopHig
 267   // 6 = probLow  (not used for electrons)
 268   // 7 = probHigh (not used for electrons)
 269   //
 270   // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
 271   // so first the correct centrality bin has to be found
 272
 273   // **** Centrality bins (hard coded for the moment)
 274   const Int_t nCent = 7;
 275   Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
 276
 277   if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
 278   if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
 279
 280   TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
 281   if(!particlePar) return NULL;
 282
 283   const TVectorD *parameters = NULL;
 284   Double_t momMin = 0.;
 285   Double_t momMax = 0.;
 286
 287   // is the centrality dependent parametrization used
 288   TString arrayName = particlePar->GetName();
 289
 290   // centrality dependent parametrization
 291   if(arrayName.Contains("Centrality")){
 292
 293     for(Int_t iCent = 0; iCent < nCent; iCent++){
 294
 295       if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
 296
 297         TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
 298         if(!centPar) return NULL;
 299
 300         TIter centIter(centPar);
 301         parameters = NULL;
 302         momMin = 0.;
 303         momMax = 0.;
 304
 305         while((parameters = static_cast<const TVectorD *>(centIter()))){
 306
 307           momMin = (*parameters)[0];
 308           momMax = (*parameters)[1];
 309
 310           if( fPt > momMin && fPt < momMax ) return parameters;
 311
 312         }
 313       }
 314     }
 315   }
 316
 317   // NO centrality dependent parametrization
 318   else{
 319
 320     TIter parIter(particlePar);
 321     while((parameters = static_cast<const TVectorD *>(parIter()))){
 322
 323       momMin = (*parameters)[0];
 324       momMax = (*parameters)[1];
 325
 326       if( fPt > momMin && fPt < momMax ) return parameters;
 327
 328     }
 329   }
 330   AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
 331
 332   return parameters;
 333 }