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