STEER/STEERBase/AliITSPIDResponse.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2005-2007, 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 // ITS PID method # 1
  20 //           Implementation of the ITS PID class
  21 // Very naive one... Should be made better by the detector experts...
  22 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  23 //-----------------------------------------------------------------
  24 #include "TMath.h"
  25 #include "AliVTrack.h"
  26 #include "AliITSPIDResponse.h"
  27 #include "AliITSPidParams.h"
  28 #include "AliExternalTrackParam.h"
  29
  30 ClassImp(AliITSPIDResponse)
  31
  32 AliITSPIDResponse::AliITSPIDResponse(Bool_t isMC):
  33   fRes(0.13),
  34   fKp1(15.77),
  35   fKp2(4.95),
  36   fKp3(0.312),
  37   fKp4(2.14),
  38   fKp5(0.82)
  39 {
  40   if(!isMC){
  41     fBBtpcits[0]=0.73;
  42     fBBtpcits[1]=14.68;
  43     fBBtpcits[2]=0.905;
  44     fBBtpcits[3]=1.2;
  45     fBBtpcits[4]=6.6;
  46     fBBdeu[0]=76.43; // parameters for the deuteron - tpcits - value from PbPb 2010 run (S.Trogolo - July 2014)
  47     fBBdeu[1]=-34.21;
  48     fBBdeu[2]=113.2;
  49     fBBdeu[3]=-18.12;
  50     fBBdeu[4]=0.6019;
  51     fBBtri[0]=13.34; // parameters for the triton - tpcits - value from PbPb 2010 run (S.Trogolo - July 2014)
  52     fBBtri[1]=55.17;
  53     fBBtri[2]=66.41;
  54     fBBtri[3]=-6.601;
  55     fBBtri[4]=-0.4134;
  56     fBBsa[0]=2.73198E7; //pure PHOBOS parameterization
  57     fBBsa[1]=6.92389;
  58     fBBsa[2]=1.90088E-6;
  59     fBBsa[3]=1.90088E-6;
  60     fBBsa[4]=3.40644E-7;
  61     fBBsaHybrid[0]=1.43505E7;  //PHOBOS+Polinomial parameterization
  62     fBBsaHybrid[1]=49.3402;
  63     fBBsaHybrid[2]=1.77741E-7;
  64     fBBsaHybrid[3]=1.77741E-7;
  65     fBBsaHybrid[4]=1.01311E-7;
  66     fBBsaHybrid[5]=77.2777;
  67     fBBsaHybrid[6]=33.4099;
  68     fBBsaHybrid[7]=46.0089;
  69     fBBsaHybrid[8]=-2.26583;
  70     fBBsaElectron[0]=4.05799E6;  //electrons in the ITS
  71     fBBsaElectron[1]=38.5713;
  72     fBBsaElectron[2]=1.46462E-7;
  73     fBBsaElectron[3]=1.46462E-7;
  74     fBBsaElectron[4]=4.40284E-7;
  75     fResolSA[0]=1.;   // 0 cluster tracks should not be used
  76     fResolSA[1]=0.25;  // rough values for tracks with 1
  77     fResolSA[2]=0.131;   // value from pp 2010 run (L. Milano, 16-Jun-11)
  78     fResolSA[3]=0.113; // value from pp 2010 run
  79     fResolSA[4]=0.104; // value from pp 2010 run
  80     for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
  81     fResolTPCITSDeu3[0]=0.06918; // deuteron resolution vs p
  82     fResolTPCITSDeu3[1]=0.02498; // 3 ITS clusters for PId
  83     fResolTPCITSDeu3[2]=1.1; // value from PbPb 2010 run (July 2014)
  84     fResolTPCITSDeu4[0]=0.06756;// deuteron resolution vs p
  85     fResolTPCITSDeu4[1]=0.02078; // 4 ITS clusters for PId
  86     fResolTPCITSDeu4[2]=1.05; // value from PbPb 2010 run (July 2014)
  87     fResolTPCITSTri3[0]=0.07239; // triton resolution vs p
  88     fResolTPCITSTri3[1]=0.0192; // 3 ITS clusters for PId
  89     fResolTPCITSTri3[2]=1.1; // value from PbPb 2010 run (July 2014)
  90     fResolTPCITSTri4[0]=0.06083; // triton resolution
  91     fResolTPCITSTri4[1]=0.02579; // 4 ITS clusters for PId
  92     fResolTPCITSTri4[2]=1.15; // value from PbPb 2010 run (July 2014)
  93   }else{
  94     fBBtpcits[0]=1.04;
  95     fBBtpcits[1]=27.14;
  96     fBBtpcits[2]=1.00;
  97     fBBtpcits[3]=0.964;
  98     fBBtpcits[4]=2.59;
  99     fBBdeu[0]=88.22; // parameters for the deuteron - MC (LHC14a6)
 100     fBBdeu[1]=-40.74;
 101     fBBdeu[2]=107.2;
 102     fBBdeu[3]=-8.962;
 103     fBBdeu[4]=-0.766;
 104     fBBtri[0]=100.7; //parameters for the triton - MC (LHC14a6)
 105     fBBtri[1]=-68.56;
 106     fBBtri[2]=128.2;
 107     fBBtri[3]=-15.5;
 108     fBBtri[4]=0.1833;
 109     fBBsa[0]=2.02078E7; //pure PHOBOS parameterization
 110     fBBsa[1]=14.0724;
 111     fBBsa[2]=3.84454E-7;
 112     fBBsa[3]=3.84454E-7;
 113     fBBsa[4]=2.43913E-7;
 114     fBBsaHybrid[0]=1.05381E7; //PHOBOS+Polinomial parameterization
 115     fBBsaHybrid[1]=89.3933;
 116     fBBsaHybrid[2]=2.4831E-7;
 117     fBBsaHybrid[3]=2.4831E-7;
 118     fBBsaHybrid[4]=7.80591E-8;
 119     fBBsaHybrid[5]=62.9214;
 120     fBBsaHybrid[6]=32.347;
 121     fBBsaHybrid[7]=58.7661;
 122     fBBsaHybrid[8]=-3.39869;
 123     fBBsaElectron[0]=2.26807E6; //electrons in the ITS
 124     fBBsaElectron[1]=99.985;
 125     fBBsaElectron[2]=0.000714841;
 126     fBBsaElectron[3]=0.000259585;
 127     fBBsaElectron[4]=1.39412E-7;
 128     fResolSA[0]=1.;   // 0 cluster tracks should not be used
 129     fResolSA[1]=0.25;  // rough values for tracks with 1
 130     fResolSA[2]=0.126;   // value from pp 2010 simulations (L. Milano, 16-Jun-11)
 131     fResolSA[3]=0.109; // value from pp 2010 simulations
 132     fResolSA[4]=0.097; // value from pp 2010 simulations
 133     for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
 134     fResolTPCITSDeu3[0]=0.06853; // deuteron resolution vs p
 135     fResolTPCITSDeu3[1]=0.01607; // 3 ITS clusters for PId
 136     fResolTPCITSDeu3[2]=1.08; // value from PbPb 2010 run (July 2014)
 137     fResolTPCITSDeu4[0]=0.06853;
 138     fResolTPCITSDeu4[1]=0.01607;
 139     fResolTPCITSDeu4[2]=1.08;
 140     fResolTPCITSTri3[0]=0.07239; // triton resolution vs p
 141     fResolTPCITSTri3[1]=0.0192; // 3 ITS clusters for PId
 142     fResolTPCITSTri3[2]=1.12; // value from PbPb 2010 run (July 2014)
 143     fResolTPCITSTri4[0]=0.07239; // triton resolution vs p
 144     fResolTPCITSTri4[1]=0.0192; // 3 ITS clusters for PId
 145     fResolTPCITSTri4[2]=1.12;
 146   }
 147 }
 148
 149 /*
 150 //_________________________________________________________________________
 151 AliITSPIDResponse::AliITSPIDResponse(Double_t *param):
 152   fRes(param[0]),
 153   fKp1(15.77),
 154   fKp2(4.95),
 155   fKp3(0.312),
 156   fKp4(2.14),
 157   fKp5(0.82)
 158 {
 159   //
 160   //  The main constructor
 161   //
 162   for (Int_t i=0; i<5;i++) {
 163       fBBsa[i]=0.;
 164       fBBtpcits[i]=0.;
 165       fResolSA[i]=0.;
 166       fResolTPCITS[i]=0.;
 167   }
 168 }
 169 */
 170
 171 //_________________________________________________________________________
 172 Double_t AliITSPIDResponse::BetheAleph(Double_t p, Double_t mass) const {
 173   //
 174   // returns AliExternalTrackParam::BetheBloch normalized to
 175   // fgMIP at the minimum
 176   //
 177
 178   Double_t bb=
 179     AliExternalTrackParam::BetheBlochAleph(p/mass,fKp1,fKp2,fKp3,fKp4,fKp5);
 180   return bb;
 181 }
 182
 183 //_________________________________________________________________________
 184 Double_t AliITSPIDResponse::Bethe(Double_t bg, const Double_t * const par, Bool_t isNuclei) const
 185 {
 186
 187   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 188   const Double_t gamma=bg/beta;
 189   Double_t bb=1.;
 190
 191   Double_t eff=1.0;
 192   if(bg<par[2])
 193     eff=(bg-par[3])*(bg-par[3])+par[4];
 194   else
 195     eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 196
 197   if(gamma>=0. && beta>0.){
 198     if(isNuclei){
 199       //Parameterization for deuteron between 0.4 - 1.5 GeV/c; triton between 0.58 - 1.65 GeV/c
 200       bb=par[0] + par[1]/bg + par[2]/(bg*bg) + par[3]/(bg*bg*bg) + par[4]/(bg*bg*bg*bg);
 201     }else{ //Parameterization for pion, kaon, proton, electron
 202       bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 203     }
 204   }
 205
 206   return bb;
 207 }
 208
 209 //_________________________________________________________________________
 210 Double_t AliITSPIDResponse::Bethe(Double_t p, Double_t mass, Bool_t isSA) const {
 211
 212   //OLD - Mantained for backward compatibility
 213   //from the MASS check --> Set the Particle Type
 214   //at the end use the method Bethe(Double_t p, AliPID::EParticleType species, Bool_t isSA) const to set the right parameter
 215
 216   //
 217   // returns AliExternalTrackParam::BetheBloch normalized to
 218   // fgMIP at the minimum
 219   //
 220
 221   // NEW: Parameterization for Deuteron and Triton energy loss, reproduced with a polynomial in fixed p range
 222   // fBBdeu --> parameters for deuteron
 223   // fBBtri --> parameters for triton
 224
 225   //NOTE
 226   //NOTE: if changes are made here, please also check the alternative function below
 227   //NOTE
 228
 229   AliPID::EParticleType species = AliPID::kPion;
 230   Bool_t foundMatchingSpecies = kFALSE;
 231   for (Int_t spec = 0; spec < AliPID::kSPECIESC; spec++) {
 232     if (TMath::AreEqualAbs(mass,AliPID::ParticleMassZ(spec),0.001)){
 233       species = (AliPID::EParticleType)spec;
 234       foundMatchingSpecies = kTRUE;
 235       break;
 236     }
 237   }
 238   if (!foundMatchingSpecies)
 239     printf("Error AliITSPIDResponse::Bethe: Mass does not match any species. Assuming pion! Note that this function is deprecated!\n");
 240
 241     return Bethe(p,species,isSA);
 242 }
 243
 244 //_________________________________________________________________________
 245 Double_t AliITSPIDResponse::Bethe(Double_t p, AliPID::EParticleType species, Bool_t isSA) const
 246 {
 247   // NEW - to be used
 248   // **** ATTENTION: the second parameter must be the PARTICLE TYPE you want to identify ****
 249   // Alternative bethe function assuming a particle type not a mass
 250   // should be slightly faster
 251   //
 252
 253   const Double_t m=AliPID::ParticleMassZ(species);
 254   const Double_t bg=p/m;
 255   Bool_t isNuclei=kFALSE;
 256
 257   //NOTE
 258   //NOTE: if changes are made here, please also check the alternative function above
 259   //NOTE
 260   const Double_t *par=fBBtpcits;
 261   if(isSA){
 262     if(species == AliPID::kElectron){
 263       //if is an electron use a specific BB parameterization
 264       //To be used only between 100 and 160 MeV/c
 265       par=fBBsaElectron;
 266     }else{
 267       par=fBBsa;
 268     }
 269   }else{
 270     if(species == AliPID::kDeuteron) {
 271       par=fBBdeu;
 272       isNuclei=kTRUE;
 273     }
 274     if(species == AliPID::kTriton  ) {
 275       par=fBBtri;
 276       isNuclei=kTRUE;
 277     }
 278   }
 279
 280   return Bethe(bg, par, isNuclei);
 281 }
 282
 283 //_________________________________________________________________________
 284 Double_t AliITSPIDResponse::BetheITSsaHybrid(Double_t p, Double_t mass) const {
 285   //
 286   // returns AliExternalTrackParam::BetheBloch normalized to
 287   // fgMIP at the minimum. The PHOBOS parameterization is used for beta*gamma>0.76.
 288   // For beta*gamma<0.76 a polinomial function is used
 289
 290   Double_t bg=p/mass;
 291   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 292   Double_t gamma=bg/beta;
 293   Double_t bb=1.;
 294
 295   Double_t par[9];
 296   //parameters for pi, K, p
 297   for(Int_t ip=0; ip<9;ip++) par[ip]=fBBsaHybrid[ip];
 298   //if it is an electron the PHOBOS part of the parameterization is tuned for e
 299   //in the range used for identification beta*gamma is >0.76 for electrons
 300   //To be used only between 100 and 160 MeV/c
 301   if(mass>0.0005 && mass<0.00052)for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
 302
 303   if(gamma>=0. && beta>0. && bg>0.1){
 304     if(bg>0.76){//PHOBOS
 305       Double_t eff=1.0;
 306       if(bg<par[2])
 307         eff=(bg-par[3])*(bg-par[3])+par[4];
 308       else
 309         eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 310
 311       bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 312     }else{//Polinomial
 313       bb=par[5] + par[6]/bg + par[7]/(bg*bg) + par[8]/(bg*bg*bg);
 314     }
 315   }
 316   return bb;
 317 }
 318
 319 //_________________________________________________________________________
 320 Double_t AliITSPIDResponse::GetResolution(Double_t bethe,
 321                                           Int_t nPtsForPid,
 322                                          Bool_t isSA,
 323                                          Double_t p,
 324                                          AliPID::EParticleType type) const {
 325   //
 326   // Calculate expected resolution for truncated mean
 327   //
 328   // NEW: Added new variables which are Double_t p and AliPID::EParticleType type
 329   // AliPID::EParticleType type is used to set the correct resolution for the different particles
 330   // default -> AliPID::EParticleType type = AliPID::kPion
 331   // Double_t p is used for the resolution of deuteron and triton, because they are function of the momentum
 332   // default -> Double_t p=0.
 333
 334   Float_t r=0.f;
 335   Double_t c=1.; //this is a correction factor used for the nuclei resolution, while for pion/kaon/proton/electron is 1.
 336
 337   if(isSA) r=fResolSA[nPtsForPid];
 338   else{
 339     const Double_t *par=0x0;
 340     if(type==AliPID::kDeuteron){
 341       if(nPtsForPid==4) par = fResolTPCITSDeu4;
 342       else par = fResolTPCITSDeu3;
 343       c=par[2];
 344       r=par[0]+par[1]*p;
 345     } else if(type==AliPID::kTriton){
 346       if(nPtsForPid==4) par = fResolTPCITSTri4;
 347       else par = fResolTPCITSTri3;
 348       c=par[2];
 349       r=par[0]+par[1]*p;
 350     } else{
 351       r=fResolTPCITS[nPtsForPid];
 352     }
 353   }
 354
 355   return r*bethe*c;
 356 }
 357
 358
 359 //_________________________________________________________________________
 360 void AliITSPIDResponse::GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES], Bool_t isMC) const {
 361   //
 362   // Method to calculate PID probabilities for a single track
 363   // using the likelihood method
 364   //
 365   const Int_t nLay = 4;
 366   const Int_t nPart= 4;
 367
 368   static AliITSPidParams pars(isMC);  // Pid parametrisation parameters
 369
 370   Double_t itsProb[nPart] = {1,1,1,1}; // e, p, K, pi
 371
 372   for (Int_t iLay = 0; iLay < nLay; iLay++) {
 373     if (qclu[iLay] <= 50.)
 374       continue;
 375
 376     Float_t dedx = qclu[iLay];
 377     Float_t layProb = pars.GetLandauGausNorm(dedx,AliPID::kProton,mom,iLay+3);
 378     itsProb[0] *= layProb;
 379
 380     layProb = pars.GetLandauGausNorm(dedx,AliPID::kKaon,mom,iLay+3);
 381     itsProb[1] *= layProb;
 382
 383     layProb = pars.GetLandauGausNorm(dedx,AliPID::kPion,mom,iLay+3);
 384     itsProb[2] *= layProb;
 385
 386     layProb = pars.GetLandauGausNorm(dedx,AliPID::kElectron,mom,iLay+3);
 387     itsProb[3] *= layProb;
 388   }
 389
 390   // Normalise probabilities
 391   Double_t sumProb = 0;
 392   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 393     sumProb += itsProb[iPart];
 394   }
 395   sumProb += itsProb[2]; // muon cannot be distinguished from pions
 396
 397   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 398     itsProb[iPart]/=sumProb;
 399   }
 400   condprobfun[AliPID::kElectron] = itsProb[3];
 401   condprobfun[AliPID::kMuon] = itsProb[2];
 402   condprobfun[AliPID::kPion] = itsProb[2];
 403   condprobfun[AliPID::kKaon] = itsProb[1];
 404   condprobfun[AliPID::kProton] = itsProb[0];
 405   return;
 406 }
 407
 408 //_________________________________________________________________________
 409 Double_t AliITSPIDResponse::GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType type) const
 410 {
 411   //
 412   // number of sigmas
 413   //
 414   UChar_t clumap=track->GetITSClusterMap();
 415   Int_t nPointsForPid=0;
 416   for(Int_t i=2; i<6; i++){
 417     if(clumap&(1<<i)) ++nPointsForPid;
 418   }
 419   Float_t mom=track->P();
 420
 421   //check for ITS standalone tracks
 422   Bool_t isSA=kTRUE;
 423   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 424
 425   const Float_t dEdx=track->GetITSsignal();
 426
 427   //TODO: in case of the electron, use the SA parametrisation,
 428   //      this needs to be changed if ITS provides a parametrisation
 429   //      for electrons also for ITS+TPC tracks
 430   return GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
 431 }
 432
 433 //_________________________________________________________________________
 434 Double_t AliITSPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
 435 {
 436   //
 437   // Signal - expected
 438   //
 439   const Float_t mom=track->P();
 440   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
 441   Bool_t isSA=kTRUE;
 442   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 443
 444   const Float_t dEdx=track->GetITSsignal();
 445
 446   //TODO: in case of the electron, use the SA parametrisation,
 447   //      this needs to be changed if ITS provides a parametrisation
 448   //      for electrons also for ITS+TPC tracks
 449
 450   const Float_t bethe = Bethe(mom,type, isSA || (type==AliPID::kElectron))*chargeFactor;
 451
 452   Double_t delta=-9999.;
 453   if (!ratio) delta=dEdx-bethe;
 454   else if (bethe>1.e-20) delta=dEdx/bethe;
 455
 456   return delta;
 457 }
 458
 459 //_________________________________________________________________________
 460 Int_t AliITSPIDResponse::GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA) const{
 461   // method to get particle identity with simple cuts on dE/dx vs. momentum
 462
 463   Double_t massp=AliPID::ParticleMass(AliPID::kProton);
 464   Double_t massk=AliPID::ParticleMass(AliPID::kKaon);
 465   Double_t bethep=Bethe(mom,massp,isSA);
 466   Double_t bethek=Bethe(mom,massk,isSA);
 467   if(signal>(0.5*(bethep+bethek))) return AliPID::kProton;
 468   Double_t masspi=AliPID::ParticleMass(AliPID::kPion);
 469   Double_t bethepi=Bethe(mom,masspi,isSA);
 470   if(signal>(0.5*(bethepi+bethek))) return AliPID::kKaon;
 471   return AliPID::kPion;
 472
 473 }