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     fBBsa[0]=2.02078E7; //pure PHOBOS parameterization
 100     fBBsa[1]=14.0724;
 101     fBBsa[2]=3.84454E-7;
 102     fBBsa[3]=3.84454E-7;
 103     fBBsa[4]=2.43913E-7;
 104     fBBsaHybrid[0]=1.05381E7; //PHOBOS+Polinomial parameterization
 105     fBBsaHybrid[1]=89.3933;
 106     fBBsaHybrid[2]=2.4831E-7;
 107     fBBsaHybrid[3]=2.4831E-7;
 108     fBBsaHybrid[4]=7.80591E-8;
 109     fBBsaHybrid[5]=62.9214;
 110     fBBsaHybrid[6]=32.347;
 111     fBBsaHybrid[7]=58.7661;
 112     fBBsaHybrid[8]=-3.39869;
 113     fBBsaElectron[0]=2.26807E6; //electrons in the ITS
 114     fBBsaElectron[1]=99.985;
 115     fBBsaElectron[2]=0.000714841;
 116     fBBsaElectron[3]=0.000259585;
 117     fBBsaElectron[4]=1.39412E-7;
 118     fResolSA[0]=1.;   // 0 cluster tracks should not be used
 119     fResolSA[1]=0.25;  // rough values for tracks with 1
 120     fResolSA[2]=0.126;   // value from pp 2010 simulations (L. Milano, 16-Jun-11)
 121     fResolSA[3]=0.109; // value from pp 2010 simulations
 122     fResolSA[4]=0.097; // value from pp 2010 simulations
 123     for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
 124   }
 125 }
 126
 127 /*
 128 //_________________________________________________________________________
 129 AliITSPIDResponse::AliITSPIDResponse(Double_t *param):
 130   fRes(param[0]),
 131   fKp1(15.77),
 132   fKp2(4.95),
 133   fKp3(0.312),
 134   fKp4(2.14),
 135   fKp5(0.82)
 136 {
 137   //
 138   //  The main constructor
 139   //
 140   for (Int_t i=0; i<5;i++) {
 141       fBBsa[i]=0.;
 142       fBBtpcits[i]=0.;
 143       fResolSA[i]=0.;
 144       fResolTPCITS[i]=0.;
 145   }
 146 }
 147 */
 148
 149 //_________________________________________________________________________
 150 Double_t AliITSPIDResponse::BetheAleph(Double_t p, Double_t mass) const {
 151   //
 152   // returns AliExternalTrackParam::BetheBloch normalized to
 153   // fgMIP at the minimum
 154   //
 155
 156   Double_t bb=
 157     AliExternalTrackParam::BetheBlochAleph(p/mass,fKp1,fKp2,fKp3,fKp4,fKp5);
 158   return bb;
 159 }
 160
 161 //_________________________________________________________________________
 162 Double_t AliITSPIDResponse::Bethe(Double_t bg, const Double_t * const par, Bool_t isNuclei) const
 163 {
 164
 165   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 166   const Double_t gamma=bg/beta;
 167   Double_t bb=1.;
 168
 169   Double_t eff=1.0;
 170   if(bg<par[2])
 171     eff=(bg-par[3])*(bg-par[3])+par[4];
 172   else
 173     eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 174
 175   if(gamma>=0. && beta>0.){
 176     if(isNuclei){
 177       //Parameterization for deuteron between 0.4 - 1.5 GeV/c; triton between 0.58 - 1.65 GeV/c
 178       bb=par[0] + par[1]/bg + par[2]/(bg*bg) + par[3]/(bg*bg*bg) + par[4]/(bg*bg*bg*bg);
 179     }else{ //Parameterization for pion, kaon, proton, electron
 180       bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 181     }
 182   }
 183
 184   return bb;
 185 }
 186
 187 //_________________________________________________________________________
 188 Double_t AliITSPIDResponse::Bethe(Double_t p, Double_t mass, Bool_t isSA) const {
 189
 190   //OLD - Mantained for backward compatibility
 191   //from the mass check --> Set the Particle Type
 192   //at the end use the method Bethe(Double_t p, AliPID::EParticleType species, Bool_t isSA) const to set the right parameter
 193
 194   //
 195   // returns AliExternalTrackParam::BetheBloch normalized to
 196   // fgMIP at the minimum
 197   //
 198
 199   // NEW: Parameterization for Deuteron and Triton energy loss, reproduced with a polynomial in fixed p range
 200   // fBBdeu --> parameters for deuteron
 201   // fBBtri --> parameters for triton
 202
 203   //NOTE
 204   //NOTE: if changes are made here, please also check the alternative function below
 205   //NOTE
 206
 207   AliPID::EParticleType species = AliPID::kPion;
 208
 209     if(TMath::AreEqualAbs(mass,AliPID::ParticleMass(0),0.00001)){
 210       //if is an electron use a specific BB parameterization
 211       //To be used only between 100 and 160 MeV/c
 212       species=AliPID::kElectron;
 213     }
 214
 215     if(TMath::AreEqualAbs(mass,AliPID::ParticleMass(5),0.002)) species=AliPID::kDeuteron;
 216     if(TMath::AreEqualAbs(mass,AliPID::ParticleMass(6),0.001)) species=AliPID::kTriton;
 217
 218     return Bethe(p,species,isSA);
 219 }
 220
 221 //_________________________________________________________________________
 222 Double_t AliITSPIDResponse::Bethe(Double_t p, AliPID::EParticleType species, Bool_t isSA) const
 223 {
 224   // NEW - to be used
 225   // Alternative bethe function assuming a particle type not a mass
 226   // should be slightly faster
 227   //
 228
 229   const Double_t m=AliPID::ParticleMassZ(species);
 230   const Double_t bg=p/m;
 231   Bool_t isNuclei=kFALSE;
 232
 233   //NOTE
 234   //NOTE: if changes are made here, please also check the alternative function above
 235   //NOTE
 236   const Double_t *par=fBBtpcits;
 237   if(isSA){
 238     if(species == AliPID::kElectron){
 239       //if is an electron use a specific BB parameterization
 240       //To be used only between 100 and 160 MeV/c
 241       par=fBBsaElectron;
 242     }else{
 243       par=fBBsa;
 244     }
 245   }else{
 246     if(species == AliPID::kDeuteron) {
 247       par=fBBdeu;
 248       isNuclei=kTRUE;
 249     }
 250     if(species == AliPID::kTriton  ) {
 251       par=fBBtri;
 252       isNuclei=kTRUE;
 253     }
 254   }
 255
 256   return Bethe(bg, par, isNuclei);
 257 }
 258
 259 //_________________________________________________________________________
 260 Double_t AliITSPIDResponse::BetheITSsaHybrid(Double_t p, Double_t mass) const {
 261   //
 262   // returns AliExternalTrackParam::BetheBloch normalized to
 263   // fgMIP at the minimum. The PHOBOS parameterization is used for beta*gamma>0.76.
 264   // For beta*gamma<0.76 a polinomial function is used
 265
 266   Double_t bg=p/mass;
 267   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 268   Double_t gamma=bg/beta;
 269   Double_t bb=1.;
 270
 271   Double_t par[9];
 272   //parameters for pi, K, p
 273   for(Int_t ip=0; ip<9;ip++) par[ip]=fBBsaHybrid[ip];
 274   //if it is an electron the PHOBOS part of the parameterization is tuned for e
 275   //in the range used for identification beta*gamma is >0.76 for electrons
 276   //To be used only between 100 and 160 MeV/c
 277   if(mass>0.0005 && mass<0.00052)for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
 278
 279   if(gamma>=0. && beta>0. && bg>0.1){
 280     if(bg>0.76){//PHOBOS
 281       Double_t eff=1.0;
 282       if(bg<par[2])
 283         eff=(bg-par[3])*(bg-par[3])+par[4];
 284       else
 285         eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 286
 287       bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 288     }else{//Polinomial
 289       bb=par[5] + par[6]/bg + par[7]/(bg*bg) + par[8]/(bg*bg*bg);
 290     }
 291   }
 292   return bb;
 293 }
 294
 295 //_________________________________________________________________________
 296 Double_t AliITSPIDResponse::GetResolution(Double_t bethe,
 297                                           Int_t nPtsForPid,
 298                                          Bool_t isSA,
 299                                          Double_t p,
 300                                          AliPID::EParticleType type) const {
 301   //
 302   // Calculate expected resolution for truncated mean
 303   //
 304   // NEW: Added new variables which are Double_t p and AliPID::EParticleType type
 305   // AliPID::EParticleType type is used to set the correct resolution for the different particles
 306   // default -> AliPID::EParticleType type = AliPID::kPion
 307   // Double_t p is used for the resolution of deuteron and triton, because they are function of the momentum
 308   // default -> Double_t p=0.
 309
 310   Float_t r=0.f;
 311   Double_t c=1.; //this is a correction factor used for the nuclei resolution, while for pion/kaon/proton/electron is 1.
 312
 313   if(isSA) r=fResolSA[nPtsForPid];
 314   else{
 315     const Double_t *par=0x0;
 316     if(type==AliPID::kDeuteron){
 317       if(nPtsForPid==3) par = fResolTPCITSDeu3;
 318       if(nPtsForPid==4) par = fResolTPCITSDeu4;
 319       c=par[2];
 320       r=par[0]+par[1]*p;
 321     } else if(type==AliPID::kTriton){
 322       if(nPtsForPid==3) par = fResolTPCITSTri3;
 323       if(nPtsForPid==4) par = fResolTPCITSTri4;
 324       c=par[2];
 325       r=par[0]+par[1]*p;
 326     } else{
 327       r=fResolTPCITS[nPtsForPid];
 328     }
 329   }
 330
 331   return r*bethe*c;
 332 }
 333
 334
 335 //_________________________________________________________________________
 336 void AliITSPIDResponse::GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES], Bool_t isMC) const {
 337   //
 338   // Method to calculate PID probabilities for a single track
 339   // using the likelihood method
 340   //
 341   const Int_t nLay = 4;
 342   const Int_t nPart= 4;
 343
 344   static AliITSPidParams pars(isMC);  // Pid parametrisation parameters
 345
 346   Double_t itsProb[nPart] = {1,1,1,1}; // e, p, K, pi
 347
 348   for (Int_t iLay = 0; iLay < nLay; iLay++) {
 349     if (qclu[iLay] <= 50.)
 350       continue;
 351
 352     Float_t dedx = qclu[iLay];
 353     Float_t layProb = pars.GetLandauGausNorm(dedx,AliPID::kProton,mom,iLay+3);
 354     itsProb[0] *= layProb;
 355
 356     layProb = pars.GetLandauGausNorm(dedx,AliPID::kKaon,mom,iLay+3);
 357     itsProb[1] *= layProb;
 358
 359     layProb = pars.GetLandauGausNorm(dedx,AliPID::kPion,mom,iLay+3);
 360     itsProb[2] *= layProb;
 361
 362     layProb = pars.GetLandauGausNorm(dedx,AliPID::kElectron,mom,iLay+3);
 363     itsProb[3] *= layProb;
 364   }
 365
 366   // Normalise probabilities
 367   Double_t sumProb = 0;
 368   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 369     sumProb += itsProb[iPart];
 370   }
 371   sumProb += itsProb[2]; // muon cannot be distinguished from pions
 372
 373   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 374     itsProb[iPart]/=sumProb;
 375   }
 376   condprobfun[AliPID::kElectron] = itsProb[3];
 377   condprobfun[AliPID::kMuon] = itsProb[2];
 378   condprobfun[AliPID::kPion] = itsProb[2];
 379   condprobfun[AliPID::kKaon] = itsProb[1];
 380   condprobfun[AliPID::kProton] = itsProb[0];
 381   return;
 382 }
 383
 384 //_________________________________________________________________________
 385 Double_t AliITSPIDResponse::GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType type) const
 386 {
 387   //
 388   // number of sigmas
 389   //
 390   UChar_t clumap=track->GetITSClusterMap();
 391   Int_t nPointsForPid=0;
 392   for(Int_t i=2; i<6; i++){
 393     if(clumap&(1<<i)) ++nPointsForPid;
 394   }
 395   Float_t mom=track->P();
 396
 397   //check for ITS standalone tracks
 398   Bool_t isSA=kTRUE;
 399   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 400
 401   const Float_t dEdx=track->GetITSsignal();
 402
 403   //TODO: in case of the electron, use the SA parametrisation,
 404   //      this needs to be changed if ITS provides a parametrisation
 405   //      for electrons also for ITS+TPC tracks
 406   return GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
 407 }
 408
 409 //_________________________________________________________________________
 410 Double_t AliITSPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
 411 {
 412   //
 413   // Signal - expected
 414   //
 415   const Float_t mom=track->P();
 416   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
 417   Bool_t isSA=kTRUE;
 418   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 419
 420   const Float_t dEdx=track->GetITSsignal();
 421
 422   //TODO: in case of the electron, use the SA parametrisation,
 423   //      this needs to be changed if ITS provides a parametrisation
 424   //      for electrons also for ITS+TPC tracks
 425
 426   const Float_t bethe = Bethe(mom,AliPID::ParticleMassZ(type), isSA || (type==AliPID::kElectron))*chargeFactor;
 427
 428   Double_t delta=-9999.;
 429   if (!ratio) delta=dEdx-bethe;
 430   else if (bethe>1.e-20) delta=dEdx/bethe;
 431
 432   return delta;
 433 }
 434
 435 //_________________________________________________________________________
 436 Int_t AliITSPIDResponse::GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA) const{
 437   // method to get particle identity with simple cuts on dE/dx vs. momentum
 438
 439   Double_t massp=AliPID::ParticleMass(AliPID::kProton);
 440   Double_t massk=AliPID::ParticleMass(AliPID::kKaon);
 441   Double_t bethep=Bethe(mom,massp,isSA);
 442   Double_t bethek=Bethe(mom,massk,isSA);
 443   if(signal>(0.5*(bethep+bethek))) return AliPID::kProton;
 444   Double_t masspi=AliPID::ParticleMass(AliPID::kPion);
 445   Double_t bethepi=Bethe(mom,masspi,isSA);
 446   if(signal>(0.5*(bethepi+bethek))) return AliPID::kKaon;
 447   return AliPID::kPion;
 448
 449 }