STEER/AliEMCALPIDResponse.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 //////////////////////////////////////////////////////////////////////////\r
  16 //                                                                      //\r
  17 // AliEMCALPIDResponse                                                  //\r
  18 //                                                                      //\r
  19 // EMCAL class to perfom PID                                            //\r
  20 // This is a prototype and still under development                      //\r
  21 //                                                                      //\r
  22 // ---------------------------------------------------------------------//\r
  23 // GetNumberOfSigmas():                                                 //\r
  24 //                                                                      //\r
  25 // Electrons:  Number of Sigmas for E/p value                           //\r
  26 //             Parametrization of LHC11a (after recalibration)          //   \r
  27 //                                                                      //\r
  28 // NON electrons:                                                       //\r
  29 //             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //\r
  30 //             --> return +/- 999                                       //\r
  31 //             Otherwise                                                //\r
  32 //             --> return nsigma (parametrization of LHC10e)            //   \r
  33 //                                                                      //\r
  34 // ---------------------------------------------------------------------//\r
  35 // ComputeEMCALProbability():                                           //\r
  36 //                                                                      //\r
  37 // Electrons:  Probability from Gaussian distribution                    //\r
  38 //                                                                      //\r
  39 // NON electrons:                                                       //\r
  40 //             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //\r
  41 //             --> probability to find particles below or above thr.    //\r
  42 //             Otherwise                                                //\r
  43 //             -->  Probability from Gaussian distribution              //  \r
  44 //                  (proper normalization to each other?)               //\r
  45 //                                                                      //\r
  46 //////////////////////////////////////////////////////////////////////////\r
  47 \r
  48 #include <TF1.h>\r
  49 #include <TMath.h>\r
  50 \r
  51 #include "AliEMCALPIDResponse.h"       //class header\r
  52 \r
  53 #include "AliLog.h"   \r
  54 \r
  55 ClassImp(AliEMCALPIDResponse)\r
  56 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
  57 AliEMCALPIDResponse::AliEMCALPIDResponse():\r
  58 fNorm(0)\r
  59 {\r
  60   //\r
  61   //  The default constructor\r
  62   //\r
  63 \r
  64   for(Int_t i = 0; i < fNptBins; i++){\r
  65 \r
  66     fPtCutMin[i] = 0.0;\r
  67 \r
  68     for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++){\r
  69 \r
  70       fMeanEoP[j][i]  = 0.0;\r
  71       fSigmaEoP[j][i] = 0.0;\r
  72       fProbLow[j][i]  = 0.0;\r
  73       fProbHigh[j][i] = 0.0;\r
  74 \r
  75     }\r
  76   } \r
  77   fPtCutMin[fNptBins] = 0.0;\r
  78 \r
  79   fNorm = new TF1("fNorm","gaus",-20,20); \r
  80 \r
  81   SetPtBoundary();\r
  82   SetParametrizations();\r
  83 }\r
  84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
  85 void AliEMCALPIDResponse::SetPtBoundary(){\r
  86   //\r
  87   // Set boundaries for momentum bins\r
  88   //\r
  89   fPtCutMin[0] = 1.5;\r
  90   fPtCutMin[1] = 2.5;\r
  91   fPtCutMin[2] = 3.5;\r
  92   fPtCutMin[3] = 4.5;\r
  93   fPtCutMin[4] = 5.5;\r
  94   fPtCutMin[5] = 6.5;\r
  95 \r
  96 }\r
  97 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
  98 void AliEMCALPIDResponse::SetParametrizations(){\r
  99 \r
 100 // This are the preliminary parametrizations (hard coded)\r
 101 // For electrons from LHC11a (Deepa Thomas)\r
 102 // For NON-electrons from LHC10e (TOF/TPC analysis)\r
 103   \r
 104   // Gaussian mean\r
 105   Float_t mean[4][6]    = {\r
 106     { 0.932, 0.997, 0.998, 1.001, 1.011, 1.011 },                              // electrons\r
 107     { 0.227804, 0.34839, 0.404077, -0.107795, -4.14584, 0.5 },             // NON electrons\r
 108     { -2.10377, 0.0582898, 0.0582898, 0.0582898, 0.0582898, 0.0582898  },   // protons\r
 109     { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}                                              // anti-protons\r
 110   };        \r
 111   \r
 112   // Gaussian sigma\r
 113   Float_t sigma[4][6]= {\r
 114     { 0.0866, 0.0693, 0.0664, 0.0583, 0.0488, 0.0515  },                      // electrons\r
 115     { 0.310831, 0.267586, 0.404077, 0.381968, 1.46183, 0.314687 },          // NON electrons\r
 116     { 0.603209, 0.255332, 0.255332, 0.255332, 0.255332, 0.255332},          // protons\r
 117     { 0.516837, 0.351516,0.351516,0.351516,0.351516,0.351516 }              // anti - protons\r
 118   };\r
 119 \r
 120   // lower probability\r
 121   Float_t probL[3][6] = {\r
 122     { 0.928689, 0.938455, 0.940448, 0.948496, 0.955981, 0.951923 },         // NON electrons\r
 123     { 0.974518, 0.978088, 0.975089, 0.975089, 0.975089,0.975089},           // protons\r
 124     { 0.824037, 0.861149, 0.898734, 0.898734, 0.898734, 0.898734},          // anti - protons\r
 125   };\r
 126 \r
 127   // upper probability\r
 128   Float_t probH[3][6] = {\r
 129     { 0.00030227, 4.04106e-05, 0.000147406, 0., 0.000956938, 0.00106838 },  // NON electrons\r
 130     { 0.000157945, 0., 0., 0., 0., 0. },                                    // protons\r
 131     { 0.00343237, 0., 0., 0., 0., 0.}                                       // anti - protons\r
 132   };\r
 133   \r
 134   \r
 135   // set parametrizations\r
 136   Int_t spec = 0;\r
 137   for (Int_t species = 0; species < 2*AliPID::kSPECIES; species++) {          // first negative particles and then positive\r
 138     for (Int_t pt = 0; pt < fNptBins; pt++){\r
 139 \r
 140       switch(species){\r
 141       case 0:      // electrons\r
 142         spec = 0;\r
 143         break;\r
 144       case 4:      // anti - protons\r
 145         spec = 3;\r
 146         break;\r
 147       case 5:      // positrons\r
 148         spec = 0;\r
 149         break;\r
 150       case 9:      // protons\r
 151         spec = 2;\r
 152         break;\r
 153       default:     // NON electrons\r
 154         spec = 1;\r
 155         break;\r
 156       }\r
 157     \r
 158 \r
 159       fMeanEoP[species][pt]  = mean[spec][pt];    \r
 160       fSigmaEoP[species][pt] = sigma[spec][pt];  \r
 161       if( spec == 0) {     // electrons have NO lower and upper probability thresholds --> set to 0\r
 162         fProbLow[species][pt]  = 0.;\r
 163         fProbHigh[species][pt] = 0.;    \r
 164       }\r
 165       else{\r
 166         fProbLow[species][pt]  = probL[spec-1][pt];\r
 167         fProbHigh[species][pt] = probH[spec-1][pt];     \r
 168       } \r
 169     \r
 170     }//loop pt bins\r
 171   }//loop species\r
 172 }\r
 173 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 174 Int_t AliEMCALPIDResponse::GetPtBin(Float_t pt) const {\r
 175   //\r
 176   // Returns the momentum bin index\r
 177   //\r
 178 \r
 179   Int_t i = -1;\r
 180   while(pt > fPtCutMin[i+1] && i+1 < fNptBins) i++;\r
 181 \r
 182   return i;\r
 183 }\r
 184 \r
 185 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 186 Double_t AliEMCALPIDResponse::GetExpectedSignal( Float_t pt, AliPID::EParticleType n, Int_t charge) const  {\r
 187   //\r
 188   // Calculates the expected PID signal as the function of \r
 189   // the information stored in the track, for the specified particle type \r
 190   //  \r
 191 \r
 192   Double_t signal = 0.;\r
 193 \r
 194   // Check the charge\r
 195   if( charge != -1 && charge != 1){\r
 196     \r
 197     return signal;\r
 198   }\r
 199   \r
 200   // Get the pt bin\r
 201   Int_t ptBin = GetPtBin(pt);\r
 202   \r
 203   // Get the species (first negative , then positive)\r
 204   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 205 \r
 206   // Get the signal\r
 207   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 208     signal = fMeanEoP[species][ptBin];\r
 209   }\r
 210 \r
 211   return signal;\r
 212 \r
 213 }\r
 214 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 215 Double_t AliEMCALPIDResponse::GetExpectedSigma( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
 216   //\r
 217   // Calculates the expected sigma of the PID signal as the function of \r
 218   // the information stored in the track, for the specified particle type \r
 219   //  \r
 220   //\r
 221   \r
 222   Double_t sigma = 999.;\r
 223 \r
 224   // Check the charge\r
 225   if( charge != -1 && charge != 1){\r
 226     \r
 227     return sigma;\r
 228   }  \r
 229 \r
 230   // Get the pt bin\r
 231   Int_t ptBin = GetPtBin(pt);\r
 232   \r
 233   // Get the species (first negative , then positive)\r
 234   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 235 \r
 236   // Get the sigma\r
 237   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 238     sigma = fSigmaEoP[species][ptBin];\r
 239   }\r
 240 \r
 241   return sigma;\r
 242  \r
 243 }\r
 244 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 245 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
 246   //\r
 247   // Calculates the expected sigma of the PID signal as the function of \r
 248   // the information stored in the track, for the specified particle type \r
 249   //  \r
 250   //\r
 251   \r
 252   Double_t norm = 1.;\r
 253 \r
 254   // Check the charge\r
 255   if( charge != -1 && charge != 1){\r
 256     \r
 257     return norm;\r
 258   }\r
 259 \r
 260   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )\r
 261   fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*GetExpectedSigma(pt,n,charge)*GetExpectedSigma(pt,n,charge)),GetExpectedSignal(pt,n,charge),GetExpectedSigma(pt,n,charge));\r
 262   norm = 1./fNorm->Integral(fLowEoP,fHighEoP)*(1-GetLowProb(pt,n,charge)-GetHighProb(pt,n,charge));\r
 263 \r
 264   return norm;\r
 265 }\r
 266 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 267 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {\r
 268       \r
 269   Double_t mean  = GetExpectedSignal(pt,n,charge);\r
 270   Double_t sigma = GetExpectedSigma(pt,n,charge);\r
 271 \r
 272   // if electron\r
 273   if(n == AliPID::kElectron){\r
 274     return (eop - mean) / sigma;\r
 275   }\r
 276 \r
 277   // if NON electron\r
 278   else{\r
 279     if ( eop < fLowEoP )\r
 280       return -999.;    // not parametrized \r
 281     else if ( eop > fHighEoP )\r
 282       return 999.;     // not parametrized \r
 283     else{\r
 284       return (eop - mean) / sigma; \r
 285     }\r
 286   }\r
 287 }\r
 288 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 289 Double_t AliEMCALPIDResponse::GetLowProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
 290   //\r
 291   //\r
 292   \r
 293   Double_t prob = 0.;\r
 294 \r
 295   // Check the charge\r
 296   if( charge != -1 && charge != 1){\r
 297     \r
 298     return prob;\r
 299   }\r
 300   \r
 301   // Get the pt bin\r
 302   Int_t ptBin = GetPtBin(pt);\r
 303   \r
 304   // Get the species (first negative , then positive)\r
 305   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 306 \r
 307   // Get the probability\r
 308   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 309     prob = fProbLow[species][ptBin];\r
 310   }\r
 311 \r
 312   return prob;\r
 313  \r
 314 }\r
 315 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 316 Double_t AliEMCALPIDResponse::GetHighProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
 317   //\r
 318   //\r
 319   \r
 320   Double_t prob = 0.;\r
 321 \r
 322   // Check the charge\r
 323   if( charge != -1 && charge != 1){\r
 324     \r
 325     return prob;\r
 326   }\r
 327   \r
 328   // Get the pt bin\r
 329   Int_t ptBin = GetPtBin(pt);\r
 330   \r
 331   // Get the species (first negative , then positive)\r
 332   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 333 \r
 334   // Get the probability\r
 335   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 336     prob = fProbHigh[species][ptBin];\r
 337   }\r
 338 \r
 339   return prob;\r
 340   \r
 341 }\r
 342 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 343 Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {\r
 344   //\r
 345   //\r
 346   Double_t fRange  = 5.0;   // hardcoded \r
 347   Double_t nsigma  = 0.0;\r
 348   Double_t sigma   = 0.0;\r
 349   Double_t proba   = -1.;\r
 350 \r
 351 \r
 352   // Check the charge\r
 353   if( charge != -1 && charge != 1){\r
 354     \r
 355     return proba;\r
 356   }\r
 357 \r
 358  \r
 359   // default value (will be returned, if pt below threshold)\r
 360   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
 361     pEMCAL[species] = -1;\r
 362   }\r
 363 \r
 364   if( GetPtBin(pt) > -1 ){\r
 365 \r
 366     // set E/p range\r
 367     if(eop < 0.05) eop = 0.05;\r
 368     if(eop > 2.00) eop = 2.00;\r
 369 \r
 370     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
 371 \r
 372       AliPID::EParticleType type = AliPID::EParticleType(species);\r
 373 \r
 374       // get nsigma value for each particle type at this E/p value\r
 375       nsigma = GetNumberOfSigmas(pt,eop,type,charge);\r
 376       sigma  = GetExpectedSigma(pt,type,charge);\r
 377 \r
 378       // electrons (standard Gaussian calculation of probabilities)\r
 379       if(type == AliPID::kElectron){\r
 380         if (TMath::Abs(nsigma) > fRange) {\r
 381           pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
 382         }\r
 383         else{\r
 384           pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
 385         }\r
 386       }\r
 387       //NON electrons\r
 388       else{\r
 389         // E/p < 0.5  -->  return probability below E/p = 0.5\r
 390         if ( nsigma == -999){\r
 391           pEMCAL[species] = GetLowProb(pt,type,charge);\r
 392         }\r
 393         // E/p > 1.5  -->  return probability above E/p = 1.5\r
 394         else if ( nsigma == 999){\r
 395           pEMCAL[species] = GetHighProb(pt,type,charge);\r
 396         }\r
 397         // in parametrized region --> calculate probability for corresponding Gauss curve\r
 398         else{\r
 399           pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
 400         \r
 401           // normalize to total probability == 1\r
 402           pEMCAL[species]*=GetExpectedNorm(pt,type,charge);\r
 403         }\r
 404       }\r
 405     }\r
 406 \r
 407     // return the electron probability\r
 408     proba = pEMCAL[AliPID::kElectron];  \r
 409 \r
 410   }\r
 411 \r
 412   return proba;\r
 413 \r
 414 }\r