STEER/STEERBase/AliITSPidParams.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 // Implementation of the class to store parameters of ITS        //
  21 // response funcions for dE/dx based PID                         //
  22 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  23 //                                                               //
  24 ///////////////////////////////////////////////////////////////////
  25
  26 #include <TFormula.h>
  27 #include <TNamed.h>
  28 #include <TMath.h>
  29 #include "AliITSPidParams.h"
  30
  31 ClassImp(AliITSPidParams)
  32
  33 //______________________________________________________________________
  34 AliITSPidParams::AliITSPidParams(Bool_t isMC):
  35 TNamed("default",""),
  36   fSDDElecLandauWidth(0),
  37   fSDDElecGaussWidth(0),
  38   fSSDElecLandauWidth(0),
  39   fSSDElecGaussWidth(0),
  40   fSDDPionLandauWidth(0),
  41   fSDDPionGaussWidth(0),
  42   fSSDPionLandauWidth(0),
  43   fSSDPionGaussWidth(0),
  44   fSDDKaonLandauWidth(0),
  45   fSDDKaonGaussWidth(0),
  46   fSSDKaonLandauWidth(0),
  47   fSSDKaonGaussWidth(0),
  48   fSDDProtLandauWidth(0),
  49   fSDDProtGaussWidth(0),
  50   fSSDProtLandauWidth(0),
  51   fSSDProtGaussWidth(0)
  52 {
  53   // default constructor
  54   if (isMC) InitMC();
  55   else InitData();
  56 }
  57 //______________________________________________________________________
  58 AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
  59   TNamed(name,""),
  60   fSDDElecLandauWidth(0),
  61   fSDDElecGaussWidth(0),
  62   fSSDElecLandauWidth(0),
  63   fSSDElecGaussWidth(0),
  64   fSDDPionLandauWidth(0),
  65   fSDDPionGaussWidth(0),
  66   fSSDPionLandauWidth(0),
  67   fSSDPionGaussWidth(0),
  68   fSDDKaonLandauWidth(0),
  69   fSDDKaonGaussWidth(0),
  70   fSSDKaonLandauWidth(0),
  71   fSSDKaonGaussWidth(0),
  72   fSDDProtLandauWidth(0),
  73   fSDDProtGaussWidth(0),
  74   fSSDProtLandauWidth(0),
  75   fSSDProtGaussWidth(0)
  76 {
  77   // standard constructor
  78   if (isMC) InitMC();
  79   else InitData();
  80 }
  81 //______________________________________________________________________
  82 AliITSPidParams::~AliITSPidParams(){
  83   //
  84   if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
  85   if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
  86
  87   if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
  88   if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
  89
  90   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  91   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
  92
  93   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  94   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  95
  96   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  97   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
  98
  99   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 100   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 101
 102   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 103   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 104
 105   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 106   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 107 }
 108 //______________________________________________________________________
 109 Double_t AliITSPidParams::BetheBloch(Double_t mom, Double_t mass, const Double_t *p) const{
 110  //
 111   // This is the empirical parameterization of the Bethe-Bloch formula.
 112   // It is normalized to 1 at the minimum.
 113   //
 114   // bg - beta*gamma
 115   //
 116   Double_t bg = mom/mass;
 117   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 118   Double_t aa = TMath::Power(beta,p[3]);
 119   Double_t bb = TMath::Power(1./bg,p[4]);
 120   bb=TMath::Log(p[2]+bb);
 121   return (p[0]-aa-bb)*p[1]/aa;
 122 }
 123 //______________________________________________________________________
 124 Double_t AliITSPidParams::ExtrapolateWidth(Double_t mom, Double_t x1, Double_t y1,Double_t x2, Double_t y2) const{
 125   //
 126   // This is a linear extrapolation of Landau width and Gaussian width
 127   // for low momentum.
 128
 129  Double_t slope = (y2-y1)/(x2-x1);
 130  return slope*mom+(y1-slope*x1);
 131 }//______________________________________________________________________
 132 void AliITSPidParams::InitMC(){
 133   // initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
 134   // parameter values from LHC10d1
 135   // MPV BetheBloch parameters;
 136
 137   //sdd MC electrons parameters
 138   fSDDElecMPVBetheParams[0] = -0.0931934;
 139   fSDDElecMPVBetheParams[1] = 77.8422;
 140   fSDDElecMPVBetheParams[2] = -0.889085;
 141   fSDDElecMPVBetheParams[3] = -154.455;
 142   fSDDElecMPVBetheParams[4] = -0.000256076;
 143
 144   //ssd MC electrons parameters
 145   fSSDElecMPVBetheParams[0] = -0.0989358;
 146   fSSDElecMPVBetheParams[1] = 77.8271;
 147   fSSDElecMPVBetheParams[2] = -0.900887;
 148   fSSDElecMPVBetheParams[3] =-1241.14;
 149   fSSDElecMPVBetheParams[4] = -0.0014204;
 150
 151    // electrons
 152   if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
 153   fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
 154   fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);
 155
 156   if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
 157   fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
 158   fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);
 159
 160   if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
 161   fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
 162   fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);
 163
 164   if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
 165   fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
 166   fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
 167
 168   //sdd MC hadrons parameters
 169   fSDDHadronMPVBetheParams[0] = 1.13531;
 170   fSDDHadronMPVBetheParams[1] = -156.651;
 171   fSDDHadronMPVBetheParams[2] = 1.87562;
 172   fSDDHadronMPVBetheParams[3] = 0.45819;
 173   fSDDHadronMPVBetheParams[4] = 2.26506;
 174
 175   //ssd MC hadrons parameters
 176   fSSDHadronMPVBetheParams[0] = -0.451908;
 177   fSSDHadronMPVBetheParams[1] = -55.4368;
 178   fSSDHadronMPVBetheParams[2] = 0.984636;
 179   fSSDHadronMPVBetheParams[3] = 0.97078;
 180   fSSDHadronMPVBetheParams[4] = 2.50883;
 181
 182   // pions
 183   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
 184   fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
 185   fSDDPionLandauWidth->SetParameters(0.08026, 5.87922);
 186
 187   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
 188   fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 189   fSDDPionGaussWidth->SetParameters(-0.090495, 7.705286);
 190
 191   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
 192   fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
 193   fSSDPionLandauWidth->SetParameters(0.083882, 5.823419);
 194
 195   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
 196   fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 197   fSSDPionGaussWidth->SetParameters(-0.105218, 5.650956);
 198
 199   // kaons
 200   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
 201   fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
 202   fSDDKaonLandauWidth->SetParameters(1.430692, 5.581389);
 203
 204   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
 205   fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 206   fSDDKaonGaussWidth->SetParameters(-1.6537, 8.071832);
 207
 208   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 209   fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
 210   fSSDKaonLandauWidth->SetParameters(1.368824, 5.639291);
 211
 212   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 213   fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 214   fSSDKaonGaussWidth->SetParameters(-1.901858, 6.0932281);
 215
 216   // protons
 217   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 218   fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
 219   fSDDProtLandauWidth->SetParameters(6.529418, 5.049098);
 220
 221   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 222   fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 223   fSDDProtGaussWidth->SetParameters(-9.360599, 11.318026);
 224
 225   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 226   fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
 227   fSSDProtLandauWidth->SetParameters(6.419493, 4.925070);
 228
 229   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 230   fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 231   fSSDProtGaussWidth->SetParameters(-9.599064, 9.358656);
 232 }
 233 //______________________________________________________________________
 234 void AliITSPidParams::InitData(){
 235   // initialize TFormulas to Real Data values (=p-p ALICE Experiment)
 236   // parameter values from LHC10b
 237   // MPV BetheBloch parameters;
 238
 239  //sdd data electrons parameters
 240   fSDDElecMPVBetheParams[0] = -0.130204;
 241   fSDDElecMPVBetheParams[1] = 75.1267;
 242   fSDDElecMPVBetheParams[2] = -0.889337;
 243   fSDDElecMPVBetheParams[3] = 388.372;
 244   fSDDElecMPVBetheParams[4] = 0.00134649;
 245
 246   //ssd data electrons parameters
 247   fSSDElecMPVBetheParams[0] = -0.162773;
 248   fSSDElecMPVBetheParams[1] = 72.9393;
 249   fSSDElecMPVBetheParams[2] = -0.896944;
 250   fSSDElecMPVBetheParams[3] = 3233.02;
 251   fSSDElecMPVBetheParams[4] = 0.00146896;
 252
 253    // electrons
 254   if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
 255   fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
 256   fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);
 257
 258   if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
 259   fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
 260   fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);
 261
 262   if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
 263   fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
 264   fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);
 265
 266   if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
 267   fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
 268   fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
 269
 270   //sdd data hadrons parameters
 271   fSDDHadronMPVBetheParams[0] = -18.1867;
 272   fSDDHadronMPVBetheParams[1] = -3.45806;
 273   fSDDHadronMPVBetheParams[2] =  2.23635;
 274   fSDDHadronMPVBetheParams[3] =  2.08328;
 275   fSDDHadronMPVBetheParams[4] = -1.92331;
 276
 277   //ssd data hadrons parameters
 278   fSSDHadronMPVBetheParams[0] = -12.6459;
 279   fSSDHadronMPVBetheParams[1] = -4.84598;
 280   fSSDHadronMPVBetheParams[2] =  1.50253;
 281   fSSDHadronMPVBetheParams[3] =  2.08328;
 282   fSSDHadronMPVBetheParams[4] = -1.3719;
 283
 284   // pions
 285   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
 286   fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
 287   fSDDPionLandauWidth->SetParameters(0.07694, 5.468728);
 288
 289   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
 290   fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 291   fSDDPionGaussWidth->SetParameters(-0.098209, 10.409441);
 292
 293   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
 294   fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
 295   fSSDPionLandauWidth->SetParameters(0.071602, 5.365442);
 296
 297   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
 298   fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 299   fSSDPionGaussWidth->SetParameters(-0.098045, 7.617583);
 300
 301   // kaons
 302   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
 303   fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
 304   fSDDKaonLandauWidth->SetParameters(0.998191, 5.461668);
 305
 306   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
 307   fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
 308   fSDDKaonGaussWidth->SetParameters(1.629308, 9.851873);
 309
 310   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 311   fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
 312   fSSDKaonLandauWidth->SetParameters(0.773113, 5.618683);
 313
 314   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 315   fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
 316   fSSDKaonGaussWidth->SetParameters(1.510713, 6.862774);
 317
 318   // protons
 319   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 320   fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
 321   fSDDProtLandauWidth->SetParameters(3.561429, 5.372105);
 322
 323   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 324   fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
 325   fSDDProtGaussWidth->SetParameters(5.395926, 10.044613);
 326
 327   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 328   fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
 329   fSSDProtLandauWidth->SetParameters(2.647428, 5.678460);
 330
 331   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 332   fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
 333   fSSDProtGaussWidth->SetParameters(4.91025, 6.779763);
 334 }
 335 //_______________________________________________________________________
 336 Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
 337   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
 338   if(TMath::Abs(pdgCode)==11) return GetLandauGausNorm(dedx,AliPID::kElectron,mom,lay);
 339   else if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
 340   else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
 341   else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
 342   else return 0.;
 343 }
 344 //_______________________________________________________________________
 345 Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
 346   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
 347
 348   Double_t par[4];
 349   Bool_t isSet=kFALSE;
 350   if(partType==AliPID::kElectron){
 351     if(lay==3 || lay==4){
 352       par[0]=GetSDDElecLandauWidth(mom);
 353       par[1]=GetSDDElecMPV(mom);
 354       par[2]=GetSDDElecGaussWidth(mom);
 355       isSet=kTRUE;
 356     }
 357     else if(lay==5 || lay==6){
 358       par[0]=GetSSDElecLandauWidth(mom);
 359       par[1]=GetSSDElecMPV(mom);
 360       par[2]=GetSSDElecGaussWidth(mom);
 361       isSet=kTRUE;
 362     }
 363   }else if(partType==AliPID::kPion){
 364     if(lay==3 || lay==4){
 365       par[0]=GetSDDPionLandauWidth(mom);
 366       par[1]=GetSDDPionMPV(mom);
 367       par[2]=GetSDDPionGaussWidth(mom);
 368       isSet=kTRUE;
 369     }
 370     else if(lay==5 || lay==6){
 371       par[0]=GetSSDPionLandauWidth(mom);
 372       par[1]=GetSSDPionMPV(mom);
 373       par[2]=GetSSDPionGaussWidth(mom);
 374       isSet=kTRUE;
 375     }
 376   }else if(partType==AliPID::kKaon){
 377     if(lay==3 || lay==4){
 378       par[0]=GetSDDKaonLandauWidth(mom);
 379       par[1]=GetSDDKaonMPV(mom);
 380       par[2]=GetSDDKaonGaussWidth(mom);
 381       isSet=kTRUE;
 382     }
 383     else if(lay==5 || lay==6){
 384       par[0]=GetSSDKaonLandauWidth(mom);
 385       par[1]=GetSSDKaonMPV(mom);
 386       par[2]=GetSSDKaonGaussWidth(mom);
 387       isSet=kTRUE;
 388     }
 389   }else if(partType==AliPID::kProton){
 390     if(lay==3 || lay==4){
 391       par[0]=GetSDDProtLandauWidth(mom);
 392       par[1]=GetSDDProtMPV(mom);
 393       par[2]=GetSDDProtGaussWidth(mom);
 394       isSet=kTRUE;
 395     }
 396     else if(lay==5 || lay==6){
 397       par[0]=GetSSDProtLandauWidth(mom);
 398       par[1]=GetSSDProtMPV(mom);
 399       par[2]=GetSSDProtGaussWidth(mom);
 400       isSet=kTRUE;
 401     }
 402   }
 403   if(!isSet) return 0.;
 404   // Numeric constants
 405   const Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 406   const Double_t mpshift  = -0.22278298;       // Landau maximum location
 407   // Control constants
 408   const Double_t np = 100.0;      // number of convolution steps
 409   const Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 410   // Variables
 411   Double_t xx;
 412   Double_t mpc;
 413   Double_t fland;
 414   Double_t sum = 0.0;
 415   Double_t xlow,xupp;
 416   Double_t step;
 417   Double_t i;
 418
 419   // MP shift correction
 420   mpc = par[1] - mpshift * par[0];
 421   // Range of convolution integral
 422   xlow = dedx - sc * par[2];
 423   xupp = dedx + sc * par[2];
 424   step = (xupp-xlow) / np;
 425
 426   // Convolution integral of Landau and Gaussian by sum
 427   for(i=1.0; i<=np/2; i++) {
 428     xx = xlow + (i-.5) * step;
 429
 430     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 431     sum += fland * TMath::Gaus(dedx,xx,par[2]);
 432
 433     xx = xupp - (i-.5) * step;
 434     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 435     sum += fland * TMath::Gaus(dedx,xx,par[2]);
 436   }
 437
 438   return (step * sum * invsq2pi / par[2]);
 439 }
 440