STEER/STEERBase/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   TObject(),\r
  59   fNorm(NULL)\r
  60 {\r
  61   //\r
  62   //  The default constructor\r
  63   //\r
  64 \r
  65   for(Int_t i = 0; i < fNptBins; i++){\r
  66 \r
  67     fPtCutMin[i] = 0.0;\r
  68 \r
  69     for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++){\r
  70 \r
  71       fMeanEoP[j][i]  = 0.0;\r
  72       fSigmaEoP[j][i] = 0.0;\r
  73       fProbLow[j][i]  = 0.0;\r
  74       fProbHigh[j][i] = 0.0;\r
  75 \r
  76     }\r
  77   } \r
  78   fPtCutMin[fNptBins] = 0.0;\r
  79 \r
  80   fNorm = new TF1("fNorm","gaus",-20,20); \r
  81 \r
  82   SetPtBoundary();\r
  83   SetParametrizations();\r
  84 }\r
  85 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
  86 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):\r
  87   TObject(other),\r
  88   fNorm(other.fNorm)\r
  89 {\r
  90   //\r
  91   //  The copy constructor\r
  92   //\r
  93     for(Int_t i = 0; i < fNptBins; i++)\r
  94     {\r
  95         fPtCutMin[i] = 0.0;\r
  96         for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++)\r
  97         {\r
  98             fMeanEoP[j][i]  = 0.0;\r
  99             fSigmaEoP[j][i] = 0.0;\r
 100             fProbLow[j][i]  = 0.0;\r
 101             fProbHigh[j][i] = 0.0;\r
 102         }\r
 103     }\r
 104   \r
 105     fPtCutMin[fNptBins] = 0.0;\r
 106     SetPtBoundary();\r
 107     SetParametrizations();\r
 108 }\r
 109 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 110 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)\r
 111 {\r
 112   //\r
 113   //  The assignment operator\r
 114   //\r
 115 \r
 116   if(this == &other) return *this;\r
 117   \r
 118   // Make copy\r
 119   TObject::operator=(other);\r
 120   fNorm = other.fNorm;\r
 121 \r
 122   for(Int_t i = 0; i < fNptBins; i++)\r
 123     {\r
 124         fPtCutMin[i] = 0.0;\r
 125         for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++)\r
 126         {\r
 127             fMeanEoP[j][i]  = 0.0;\r
 128             fSigmaEoP[j][i] = 0.0;\r
 129             fProbLow[j][i]  = 0.0;\r
 130             fProbHigh[j][i] = 0.0;\r
 131         }\r
 132     }\r
 133   \r
 134   fPtCutMin[fNptBins] = 0.0;\r
 135   SetPtBoundary();\r
 136   SetParametrizations();\r
 137 \r
 138   return *this;\r
 139 }\r
 140 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 141 AliEMCALPIDResponse::~AliEMCALPIDResponse() {\r
 142 \r
 143   delete fNorm;\r
 144 \r
 145 }\r
 146 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 147 void AliEMCALPIDResponse::SetPtBoundary(){\r
 148   //\r
 149   // Set boundaries for momentum bins\r
 150   //\r
 151   fPtCutMin[0] = 1.5;\r
 152   fPtCutMin[1] = 2.5;\r
 153   fPtCutMin[2] = 3.5;\r
 154   fPtCutMin[3] = 4.5;\r
 155   fPtCutMin[4] = 5.5;\r
 156   fPtCutMin[5] = 6.5;\r
 157 \r
 158 }\r
 159 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 160 void AliEMCALPIDResponse::SetParametrizations(){\r
 161 \r
 162 // This are the preliminary parametrizations (hard coded)\r
 163 // For electrons from LHC11a (Deepa Thomas)\r
 164 // For NON-electrons from LHC10e (TOF/TPC analysis)\r
 165   \r
 166   // Gaussian mean\r
 167   Float_t mean[4][6]    = {\r
 168     { 0.932, 0.997, 0.998, 1.001, 1.011, 1.011 },                              // electrons\r
 169     { 0.227804, 0.34839, 0.404077, -0.107795, -4.14584, 0.5 },             // NON electrons\r
 170     { -2.10377, 0.0582898, 0.0582898, 0.0582898, 0.0582898, 0.0582898  },   // protons\r
 171     { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}                                              // anti-protons\r
 172   };        \r
 173   \r
 174   // Gaussian sigma\r
 175   Float_t sigma[4][6]= {\r
 176     { 0.0866, 0.0693, 0.0664, 0.0583, 0.0488, 0.0515  },                      // electrons\r
 177     { 0.310831, 0.267586, 0.404077, 0.381968, 1.46183, 0.314687 },          // NON electrons\r
 178     { 0.603209, 0.255332, 0.255332, 0.255332, 0.255332, 0.255332},          // protons\r
 179     { 0.516837, 0.351516,0.351516,0.351516,0.351516,0.351516 }              // anti - protons\r
 180   };\r
 181 \r
 182   // lower probability\r
 183   Float_t probL[3][6] = {\r
 184     { 0.928689, 0.938455, 0.940448, 0.948496, 0.955981, 0.951923 },         // NON electrons\r
 185     { 0.974518, 0.978088, 0.975089, 0.975089, 0.975089,0.975089},           // protons\r
 186     { 0.824037, 0.861149, 0.898734, 0.898734, 0.898734, 0.898734},          // anti - protons\r
 187   };\r
 188 \r
 189   // upper probability\r
 190   Float_t probH[3][6] = {\r
 191     { 0.00030227, 4.04106e-05, 0.000147406, 0., 0.000956938, 0.00106838 },  // NON electrons\r
 192     { 0.000157945, 0., 0., 0., 0., 0. },                                    // protons\r
 193     { 0.00343237, 0., 0., 0., 0., 0.}                                       // anti - protons\r
 194   };\r
 195   \r
 196   \r
 197   // set parametrizations\r
 198   Int_t spec = 0;\r
 199   for (Int_t species = 0; species < 2*AliPID::kSPECIES; species++) {          // first negative particles and then positive\r
 200     for (Int_t pt = 0; pt < fNptBins; pt++){\r
 201 \r
 202       switch(species){\r
 203       case 0:      // electrons\r
 204         spec = 0;\r
 205         break;\r
 206       case 4:      // anti - protons\r
 207         spec = 3;\r
 208         break;\r
 209       case 5:      // positrons\r
 210         spec = 0;\r
 211         break;\r
 212       case 9:      // protons\r
 213         spec = 2;\r
 214         break;\r
 215       default:     // NON electrons\r
 216         spec = 1;\r
 217         break;\r
 218       }\r
 219     \r
 220 \r
 221       fMeanEoP[species][pt]  = mean[spec][pt];    \r
 222       fSigmaEoP[species][pt] = sigma[spec][pt];  \r
 223       if( spec == 0) {     // electrons have NO lower and upper probability thresholds --> set to 0\r
 224         fProbLow[species][pt]  = 0.;\r
 225         fProbHigh[species][pt] = 0.;    \r
 226       }\r
 227       else{\r
 228         fProbLow[species][pt]  = probL[spec-1][pt];\r
 229         fProbHigh[species][pt] = probH[spec-1][pt];     \r
 230       } \r
 231     \r
 232     }//loop pt bins\r
 233   }//loop species\r
 234 }\r
 235 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 236 Int_t AliEMCALPIDResponse::GetPtBin(Float_t pt) const {\r
 237   //\r
 238   // Returns the momentum bin index\r
 239   //\r
 240 \r
 241   Int_t i = -1;\r
 242   while(pt > fPtCutMin[i+1] && i+1 < fNptBins) i++;\r
 243 \r
 244   return i;\r
 245 }\r
 246 \r
 247 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 248 Double_t AliEMCALPIDResponse::GetExpectedSignal( Float_t pt, AliPID::EParticleType n, Int_t charge) const  {\r
 249   //\r
 250   // Calculates the expected PID signal as the function of \r
 251   // the information stored in the track, for the specified particle type \r
 252   //  \r
 253 \r
 254   Double_t signal = 0.;\r
 255 \r
 256   // Check the charge\r
 257   if( charge != -1 && charge != 1){\r
 258     \r
 259     return signal;\r
 260   }\r
 261   \r
 262   // Get the pt bin\r
 263   Int_t ptBin = GetPtBin(pt);\r
 264   \r
 265   // Get the species (first negative , then positive)\r
 266   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 267 \r
 268   // Get the signal\r
 269   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 270     signal = fMeanEoP[species][ptBin];\r
 271   }\r
 272 \r
 273   return signal;\r
 274 \r
 275 }\r
 276 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 277 Double_t AliEMCALPIDResponse::GetExpectedSigma( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
 278   //\r
 279   // Calculates the expected sigma of the PID signal as the function of \r
 280   // the information stored in the track, for the specified particle type \r
 281   //  \r
 282   //\r
 283   \r
 284   Double_t sigma = 999.;\r
 285 \r
 286   // Check the charge\r
 287   if( charge != -1 && charge != 1){\r
 288     \r
 289     return sigma;\r
 290   }  \r
 291 \r
 292   // Get the pt bin\r
 293   Int_t ptBin = GetPtBin(pt);\r
 294   \r
 295   // Get the species (first negative , then positive)\r
 296   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 297 \r
 298   // Get the sigma\r
 299   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 300     sigma = fSigmaEoP[species][ptBin];\r
 301   }\r
 302 \r
 303   return sigma;\r
 304  \r
 305 }\r
 306 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 307 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
 308   //\r
 309   // Calculates the expected sigma of the PID signal as the function of \r
 310   // the information stored in the track, for the specified particle type \r
 311   //  \r
 312   //\r
 313   \r
 314   Double_t norm = 1.;\r
 315 \r
 316   // Check the charge\r
 317   if( charge != -1 && charge != 1){\r
 318     \r
 319     return norm;\r
 320   }\r
 321 \r
 322   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )\r
 323   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
 324   norm = 1./fNorm->Integral(fLowEoP,fHighEoP)*(1-GetLowProb(pt,n,charge)-GetHighProb(pt,n,charge));\r
 325 \r
 326   return norm;\r
 327 }\r
 328 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 329 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {\r
 330       \r
 331   Double_t mean  = GetExpectedSignal(pt,n,charge);\r
 332   Double_t sigma = GetExpectedSigma(pt,n,charge);\r
 333 \r
 334   // if electron\r
 335   if(n == AliPID::kElectron){\r
 336     return (eop - mean) / sigma;\r
 337   }\r
 338 \r
 339   // if NON electron\r
 340   else{\r
 341     if ( eop < fLowEoP )\r
 342       return -999.;    // not parametrized \r
 343     else if ( eop > fHighEoP )\r
 344       return 999.;     // not parametrized \r
 345     else{\r
 346       return (eop - mean) / sigma; \r
 347     }\r
 348   }\r
 349 }\r
 350 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 351 Double_t AliEMCALPIDResponse::GetLowProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
 352   //\r
 353   //\r
 354   \r
 355   Double_t prob = 0.;\r
 356 \r
 357   // Check the charge\r
 358   if( charge != -1 && charge != 1){\r
 359     \r
 360     return prob;\r
 361   }\r
 362   \r
 363   // Get the pt bin\r
 364   Int_t ptBin = GetPtBin(pt);\r
 365   \r
 366   // Get the species (first negative , then positive)\r
 367   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 368 \r
 369   // Get the probability\r
 370   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 371     prob = fProbLow[species][ptBin];\r
 372   }\r
 373 \r
 374   return prob;\r
 375  \r
 376 }\r
 377 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 378 Double_t AliEMCALPIDResponse::GetHighProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
 379   //\r
 380   //\r
 381   \r
 382   Double_t prob = 0.;\r
 383 \r
 384   // Check the charge\r
 385   if( charge != -1 && charge != 1){\r
 386     \r
 387     return prob;\r
 388   }\r
 389   \r
 390   // Get the pt bin\r
 391   Int_t ptBin = GetPtBin(pt);\r
 392   \r
 393   // Get the species (first negative , then positive)\r
 394   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
 395 \r
 396   // Get the probability\r
 397   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
 398     prob = fProbHigh[species][ptBin];\r
 399   }\r
 400 \r
 401   return prob;\r
 402   \r
 403 }\r
 404 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
 405 Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {\r
 406   //\r
 407   //\r
 408   Double_t fRange  = 5.0;   // hardcoded \r
 409   Double_t nsigma  = 0.0;\r
 410   Double_t sigma   = 0.0;\r
 411   Double_t proba   = 999.;\r
 412   \r
 413 \r
 414   // Check the charge\r
 415   if( charge != -1 && charge != 1){\r
 416     \r
 417     return proba;\r
 418   }\r
 419 \r
 420  \r
 421   // default value (will be returned, if pt below threshold)\r
 422   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
 423     pEMCAL[species] = 999.;\r
 424   }\r
 425 \r
 426   if( GetPtBin(pt) > -1 ){\r
 427 \r
 428     // set E/p range\r
 429     if(eop < 0.05) eop = 0.05;\r
 430     if(eop > 2.00) eop = 2.00;\r
 431 \r
 432     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
 433 \r
 434       AliPID::EParticleType type = AliPID::EParticleType(species);\r
 435 \r
 436       // get nsigma value for each particle type at this E/p value\r
 437       nsigma = GetNumberOfSigmas(pt,eop,type,charge);\r
 438       sigma  = GetExpectedSigma(pt,type,charge);\r
 439 \r
 440       // electrons (standard Gaussian calculation of probabilities)\r
 441       if(type == AliPID::kElectron){\r
 442         if (TMath::Abs(nsigma) > fRange) {\r
 443           pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
 444         }\r
 445         else{\r
 446           pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
 447         }\r
 448       }\r
 449       //NON electrons\r
 450       else{\r
 451         // E/p < 0.5  -->  return probability below E/p = 0.5\r
 452         if ( nsigma == -999){\r
 453           pEMCAL[species] = GetLowProb(pt,type,charge);\r
 454         }\r
 455         // E/p > 1.5  -->  return probability above E/p = 1.5\r
 456         else if ( nsigma == 999){\r
 457           pEMCAL[species] = GetHighProb(pt,type,charge);\r
 458         }\r
 459         // in parametrized region --> calculate probability for corresponding Gauss curve\r
 460         else{\r
 461           pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
 462         \r
 463           // normalize to total probability == 1\r
 464           pEMCAL[species]*=GetExpectedNorm(pt,type,charge);\r
 465         }\r
 466       }\r
 467     }\r
 468 \r
 469     // return the electron probability\r
 470     proba = pEMCAL[AliPID::kElectron];  \r
 471 \r
 472   }\r
 473 \r
 474   return proba;\r
 475 \r
 476 }\r