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   fkPIDParams(NULL)
  65 {
  66   //
  67   //  The default constructor
  68   //
  69
  70
  71   fNorm = new TF1("fNorm","gaus",-20,20);
  72 }
  73 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  74 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
  75   TObject(other),
  76   fNorm(other.fNorm),
  77   fkPIDParams(other.fkPIDParams)
  78 {
  79   //
  80   //  The copy constructor
  81   //
  82
  83 }
  84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  85 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
  86 {
  87   //
  88   //  The assignment operator
  89   //
  90
  91   if(this == &other) return *this;
  92
  93   // Make copy
  94   TObject::operator=(other);
  95   fNorm = other.fNorm;
  96   fkPIDParams = other.fkPIDParams;
  97
  98
  99   return *this;
 100 }
 101 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 102 AliEMCALPIDResponse::~AliEMCALPIDResponse() {
 103
 104   delete fNorm;
 105
 106 }
 107 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 108 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {
 109   //
 110   // Calculates the expected sigma of the PID signal as the function of
 111   // the information stored in the track, for the specified particle type
 112   //
 113   //
 114
 115   Double_t norm = 1.;
 116
 117   // Check the charge
 118   if( charge != -1 && charge != 1){
 119     return norm;
 120   }
 121
 122   // Get the parameters for this particle type and pt
 123   const TVectorD *params = GetParams(n, pt, charge);
 124
 125   // IF not in momentum range, NULL is returned --> return default value
 126   if(!params) return norm;
 127
 128   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
 129   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
 130   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
 131   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
 132   Double_t probLow  = (*params)[6];   // probability to be below eopMin
 133   Double_t probHigh = (*params)[7];   // probability to be above eopMax
 134
 135   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
 136   fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
 137   norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
 138
 139   return norm;
 140 }
 141 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 142 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {
 143
 144   Double_t nsigma = -999.;
 145
 146   // Check the charge
 147   if( charge != -1 && charge != 1){
 148     return nsigma;
 149   }
 150
 151   // Get the parameters for this particle type and pt
 152   const TVectorD *params = GetParams(n, pt, charge);
 153
 154   // IF not in momentum range, NULL is returned --> return default value
 155   if(!params) return nsigma;
 156
 157   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
 158   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
 159   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
 160   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
 161
 162   // if electron
 163   if(n == AliPID::kElectron){
 164     if(sigma != 0) nsigma = (eop - mean) / sigma;
 165   }
 166
 167   // if NON electron
 168   else{
 169     if ( eop < eopMin )
 170       nsigma = -99;    // not parametrized (smaller than eopMin)
 171     else if ( eop > eopMax )
 172       nsigma = 99.;     // not parametrized (bigger than eopMax)
 173     else{
 174       if(sigma != 0) nsigma = (eop - mean) / sigma;
 175     }
 176   }
 177
 178   return nsigma;
 179
 180 }
 181 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 182 Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
 183   //
 184   //
 185   Double_t fRange  = 5.0;   // hardcoded (???)
 186   Double_t nsigma  = 0.0;
 187
 188
 189   // Check the charge
 190   if( charge != -1 && charge != 1){
 191     return kFALSE;
 192   }
 193
 194
 195   // default value (will be returned, if pt below threshold)
 196   for (Int_t species = 0; species < nSpecies; species++) {
 197     pEMCAL[species] = 1./nSpecies;
 198   }
 199
 200   // set E/p range
 201   if(eop < 0.05) eop = 0.05;
 202   if(eop > 2.00) eop = 2.00;
 203
 204   for (Int_t species = 0; species < nSpecies; species++) {
 205
 206     AliPID::EParticleType type = AliPID::EParticleType(species);
 207
 208     // Get the parameters for this particle type and pt
 209     const TVectorD *params = GetParams(species, pt, charge);
 210
 211     // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
 212     if(!params) return kFALSE;
 213
 214     Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
 215     Double_t probLow  = (*params)[6];   // probability to be below eopMin
 216     Double_t probHigh = (*params)[7];   // probability to be above eopMax
 217
 218     // get nsigma value for each particle type at this E/p value
 219     nsigma = GetNumberOfSigmas(pt,eop,type,charge);
 220
 221     // electrons (standard Gaussian calculation of probabilities)
 222     if(type == AliPID::kElectron){
 223       if (TMath::Abs(nsigma) > fRange) {
 224         pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
 225       }
 226       else{
 227         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
 228       }
 229     }
 230     //NON electrons
 231     else{
 232       // E/p < eopMin  -->  return probability below E/p = eopMin
 233       if ( nsigma == -99){
 234         pEMCAL[species] = probLow;
 235       }
 236       // E/p > eopMax  -->  return probability above E/p = eopMax
 237       else if ( nsigma == 99){
 238         pEMCAL[species] = probHigh;
 239       }
 240       // in parametrized region --> calculate probability for corresponding Gauss curve
 241       else{
 242         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
 243
 244         // normalize to total probability == 1
 245         pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
 246       }
 247     }
 248   }
 249
 250   return kTRUE;
 251
 252 }
 253 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 254 const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
 255   //
 256   // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
 257   //
 258   // 0 = momMin
 259   // 1 = momMax
 260   // 2 = mean of Gaus
 261   // 3 = sigma of Gaus
 262   // 4 = eopLow
 263   // 5 = eopHig
 264   // 6 = probLow  (not used for electrons)
 265   // 7 = probHigh (not used for electrons)
 266   //
 267
 268   if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
 269   if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
 270
 271   TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
 272   if(!particlePar) return NULL;
 273
 274   TIter parIter(particlePar);
 275   const TVectorD *parameters = NULL;
 276   Double_t momMin = 0.;
 277   Double_t momMax = 0.;
 278
 279   while((parameters = static_cast<const TVectorD *>(parIter()))){
 280
 281     momMin = (*parameters)[0];
 282     momMax = (*parameters)[1];
 283
 284     if( fPt > momMin && fPt < momMax ) return parameters;
 285
 286   }
 287   AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
 288
 289   return parameters;
 290 }