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 #include "AliPID.h"
  31
  32 ClassImp(AliITSPidParams)
  33
  34 //______________________________________________________________________
  35 AliITSPidParams::AliITSPidParams(Bool_t isMC):
  36 TNamed("default",""),
  37   fSDDPionMPV(0),
  38   fSDDPionLandauWidth(0),
  39   fSDDPionGaussWidth(0),
  40   fSSDPionMPV(0),
  41   fSSDPionLandauWidth(0),
  42   fSSDPionGaussWidth(0),
  43   fSDDKaonMPV(0),
  44   fSDDKaonLandauWidth(0),
  45   fSDDKaonGaussWidth(0),
  46   fSSDKaonMPV(0),
  47   fSSDKaonLandauWidth(0),
  48   fSSDKaonGaussWidth(0),
  49   fSDDProtMPV(0),
  50   fSDDProtLandauWidth(0),
  51   fSDDProtGaussWidth(0),
  52   fSSDProtMPV(0),
  53   fSSDProtLandauWidth(0),
  54   fSSDProtGaussWidth(0)
  55 {
  56   // default constructor
  57   if (isMC) InitMC();
  58   else InitData();
  59 }
  60 //______________________________________________________________________
  61 AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
  62   TNamed(name,""),
  63   fSDDPionMPV(0),
  64   fSDDPionLandauWidth(0),
  65   fSDDPionGaussWidth(0),
  66   fSSDPionMPV(0),
  67   fSSDPionLandauWidth(0),
  68   fSSDPionGaussWidth(0),
  69   fSDDKaonMPV(0),
  70   fSDDKaonLandauWidth(0),
  71   fSDDKaonGaussWidth(0),
  72   fSSDKaonMPV(0),
  73   fSSDKaonLandauWidth(0),
  74   fSSDKaonGaussWidth(0),
  75   fSDDProtMPV(0),
  76   fSDDProtLandauWidth(0),
  77   fSDDProtGaussWidth(0),
  78   fSSDProtMPV(0),
  79   fSSDProtLandauWidth(0),
  80   fSSDProtGaussWidth(0)
  81 {
  82   // standard constructor
  83   if (isMC) InitMC();
  84   else InitData();
  85 }
  86 //______________________________________________________________________
  87 AliITSPidParams::~AliITSPidParams(){
  88   //
  89   if(fSDDPionMPV) delete fSDDPionMPV;
  90   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  91   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
  92
  93   if(fSSDPionMPV) delete fSSDPionMPV;
  94   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  95   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  96
  97   if(fSDDKaonMPV) delete fSDDKaonMPV;
  98   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  99   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
 100
 101   if(fSSDKaonMPV) delete fSSDKaonMPV;
 102   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 103   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 104
 105   if(fSDDProtMPV) delete fSDDProtMPV;
 106   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 107   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 108
 109   if(fSSDProtMPV) delete fSSDProtMPV;
 110   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 111   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 112 }
 113
 114 //______________________________________________________________________
 115 void AliITSPidParams::InitMC(){
 116   // initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
 117   // parameter values from LHC10d1
 118
 119   // pions
 120   if(fSDDPionMPV) delete fSDDPionMPV;
 121   fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
 122   fSDDPionMPV->SetParameters(-0.892291, 0.003630, 1.866484, 78.378179);
 123
 124   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
 125   fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
 126   fSDDPionLandauWidth->SetParameters(0.080999, 5.917715);
 127
 128   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
 129   fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 130   fSDDPionGaussWidth->SetParameters(-0.092822, 7.839729);
 131
 132   if(fSSDPionMPV) delete fSSDPionMPV;
 133   fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
 134   fSSDPionMPV->SetParameters(-0.896507, 0.003173, 2.017155, 80.682567);
 135
 136   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
 137   fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
 138   fSSDPionLandauWidth->SetParameters(0.087182, 5.843610);
 139
 140   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
 141   fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 142   fSSDPionGaussWidth->SetParameters(-0.110444, 5.837737);
 143
 144   // kaons
 145   if(fSDDKaonMPV) delete fSDDKaonMPV;
 146   fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 147   fSDDKaonMPV->SetParameters(17.581590, -0.120134, 72.550701);
 148
 149   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
 150   fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
 151   fSDDKaonLandauWidth->SetParameters(1.271756, 5.778888);
 152
 153   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
 154   fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 155   fSDDKaonGaussWidth->SetParameters(-1.650298, 8.322084);
 156
 157   if(fSSDKaonMPV) delete fSSDKaonMPV;
 158   fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 159   fSSDKaonMPV->SetParameters(16.238778, 0.039318, 75.863719);
 160
 161   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 162   fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
 163   fSSDKaonLandauWidth->SetParameters(1.179541, 5.961353);
 164
 165   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 166   fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 167   fSSDKaonGaussWidth->SetParameters(-2.019126, 6.155977);
 168
 169   // protons
 170   if(fSDDProtMPV) delete fSDDProtMPV;
 171   fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 172   fSDDProtMPV->SetParameters(64.482762, -1.667823, 71.850731);
 173
 174   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 175   fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
 176   fSDDProtLandauWidth->SetParameters(6.948997, 3.928018);
 177
 178   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 179   fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 180   fSDDProtGaussWidth->SetParameters(-6.522760, 12.673959);
 181
 182   if(fSSDProtMPV) delete fSSDProtMPV;
 183   fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 184   fSSDProtMPV->SetParameters(63.817375, -1.221779, 73.233644);
 185
 186   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 187   fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
 188   fSSDProtLandauWidth->SetParameters(7.286942, 3.581451);
 189
 190   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 191   fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 192   fSSDProtGaussWidth->SetParameters(-8.327867, 9.723422);
 193
 194 }
 195 //______________________________________________________________________
 196 void AliITSPidParams::InitData(){
 197   // initialize TFormulas to Real Data values (=p-p simulations PYTHIA+GEANT)
 198   // parameter values from LHC10b
 199
 200   // pions
 201   if(fSDDPionMPV) delete fSDDPionMPV;
 202   fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
 203   fSDDPionMPV->SetParameters(-0.907609, 0.006521, 3.340347, 81.297942);
 204
 205   if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
 206   fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
 207   fSDDPionLandauWidth->SetParameters(0.077272, 5.478557);
 208
 209   if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
 210   fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 211   fSDDPionGaussWidth->SetParameters(-0.098529, 10.265711);
 212
 213   if(fSSDPionMPV) delete fSSDPionMPV;
 214   fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
 215   fSSDPionMPV->SetParameters(-0.920046, 0.006061, 3.428578, 81.401816);
 216
 217   if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
 218   fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
 219   fSSDPionLandauWidth->SetParameters(0.071243, 5.388830);
 220
 221   if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
 222   fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
 223   fSSDPionGaussWidth->SetParameters(-0.099189, 7.412309);
 224
 225   // kaons
 226   if(fSDDKaonMPV) delete fSDDKaonMPV;
 227   fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 228   fSDDKaonMPV->SetParameters(15.429146, 0.178251, 74.640293);
 229
 230   if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
 231   fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
 232   fSDDKaonLandauWidth->SetParameters(0.975202, 5.699311);
 233
 234   if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
 235   fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
 236   fSDDKaonGaussWidth->SetParameters(1.660840, 9.389192);
 237
 238   if(fSSDKaonMPV) delete fSSDKaonMPV;
 239   fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 240   fSSDKaonMPV->SetParameters(15.170715, 0.181379, 74.951884);
 241
 242   if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
 243   fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
 244   fSSDKaonLandauWidth->SetParameters(0.756466, 5.818274);
 245
 246   if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
 247   fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
 248   fSSDKaonGaussWidth->SetParameters(1.546693, 6.389872);
 249
 250   // protons
 251   if(fSDDProtMPV) delete fSDDProtMPV;
 252   fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 253   fSDDProtMPV->SetParameters(61.452534, 0.372908, 71.668352);
 254
 255   if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
 256   fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
 257   fSDDProtLandauWidth->SetParameters(3.667023, 5.430721);
 258
 259   if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
 260   fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
 261   fSDDProtGaussWidth->SetParameters(5.503814, 9.657439);
 262
 263   if(fSSDProtMPV) delete fSSDProtMPV;
 264   fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
 265   fSSDProtMPV->SetParameters(60.246538, 0.000323, 71.992031);
 266
 267   if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
 268   fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
 269   fSSDProtLandauWidth->SetParameters(2.568323, 5.939774);
 270
 271   if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
 272   fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
 273   fSSDProtGaussWidth->SetParameters(5.050541, 6.290964);
 274
 275 }
 276 //_______________________________________________________________________
 277 Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
 278   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
 279   if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
 280   else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
 281   else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
 282   else return 0.;
 283 }
 284 //_______________________________________________________________________
 285 Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
 286   // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
 287
 288   Double_t par[3];
 289   Bool_t isSet=kFALSE;
 290   if(partType==AliPID::kPion){
 291     if(lay==3 || lay==4){
 292       par[0]=GetSDDPionLandauWidth(mom);
 293       par[1]=GetSDDPionMPV(mom);
 294       par[2]=GetSDDPionGaussWidth(mom);
 295       isSet=kTRUE;
 296     }
 297     else if(lay==5 || lay==6){
 298       par[0]=GetSSDPionLandauWidth(mom);
 299       par[1]=GetSSDPionMPV(mom);
 300       par[2]=GetSSDPionGaussWidth(mom);
 301       isSet=kTRUE;
 302     }
 303   }else if(partType==AliPID::kKaon){
 304     if(lay==3 || lay==4){
 305       par[0]=GetSDDKaonLandauWidth(mom);
 306       par[1]=GetSDDKaonMPV(mom);
 307       par[2]=GetSDDKaonGaussWidth(mom);
 308       isSet=kTRUE;
 309     }
 310     else if(lay==5 || lay==6){
 311       par[0]=GetSSDKaonLandauWidth(mom);
 312       par[1]=GetSSDKaonMPV(mom);
 313       par[2]=GetSSDKaonGaussWidth(mom);
 314       isSet=kTRUE;
 315     }
 316   }else if(partType==AliPID::kProton){
 317     if(lay==3 || lay==4){
 318       par[0]=GetSDDProtLandauWidth(mom);
 319       par[1]=GetSDDProtMPV(mom);
 320       par[2]=GetSDDProtGaussWidth(mom);
 321       isSet=kTRUE;
 322     }
 323     else if(lay==5 || lay==6){
 324       par[0]=GetSSDProtLandauWidth(mom);
 325       par[1]=GetSSDProtMPV(mom);
 326       par[2]=GetSSDProtGaussWidth(mom);
 327       isSet=kTRUE;
 328     }
 329   }
 330   if(!isSet) return 0.;
 331   // Numeric constants
 332   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 333   Double_t mpshift  = -0.22278298;       // Landau maximum location
 334   // Control constants
 335   Double_t np = 100.0;      // number of convolution steps
 336   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 337   // Variables
 338   Double_t xx;
 339   Double_t mpc;
 340   Double_t fland;
 341   Double_t sum = 0.0;
 342   Double_t xlow,xupp;
 343   Double_t step = 0.;
 344   Double_t i;
 345
 346   // MP shift correction
 347   mpc = par[1] - mpshift * par[0];
 348   // Range of convolution integral
 349   xlow = dedx - sc * par[2];
 350   xupp = dedx + sc * par[2];
 351   if(np!=0) step = (xupp-xlow) / np;
 352
 353   // Convolution integral of Landau and Gaussian by sum
 354   for(i=1.0; i<=np/2; i++) {
 355     xx = xlow + (i-.5) * step;
 356
 357     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 358     sum += fland * TMath::Gaus(dedx,xx,par[2]);
 359
 360     xx = xupp - (i-.5) * step;
 361     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 362     sum += fland * TMath::Gaus(dedx,xx,par[2]);
 363   }
 364
 365   return (step * sum * invsq2pi / par[2]);
 366 }
 367