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, Bool_t isNuclei) const {
 189
 190   //
 191   // returns AliExternalTrackParam::BetheBloch normalized to
 192   // fgMIP at the minimum
 193   //
 194
 195   // NEW: Parameterization for Deuteron and Triton energy loss, reproduced with a polynomial in fixed p range
 196   // fBBdeu --> parameters for deuteron
 197   // fBBtri --> parameters for triton
 198
 199
 200   const Double_t bg=p/mass;
 201
 202   //NOTE
 203   //NOTE: if changes are made here, please also check the alternative function below
 204   //NOTE
 205   const Double_t *par=fBBtpcits;
 206   if(isSA){
 207     if(TMath::AreEqualAbs(mass,AliPID::ParticleMass(0),0.00001)){
 208       //if is an electron use a specific BB parameterization
 209       //To be used only between 100 and 160 MeV/c
 210       par=fBBsaElectron;
 211     }else{
 212       par=fBBsa;
 213     }
 214   }else{
 215     if(isNuclei){
 216       if(TMath::AreEqualAbs(mass,AliPID::ParticleMass(5),0.002)) par=fBBdeu;
 217       if(TMath::AreEqualAbs(mass,AliPID::ParticleMass(6),0.001)) par=fBBtri;
 218     }
 219   }
 220
 221   return Bethe(bg, par, isNuclei);
 222 }
 223
 224 //_________________________________________________________________________
 225 Double_t AliITSPIDResponse::Bethe(Double_t p, AliPID::EParticleType species, Bool_t isSA) const
 226 {
 227   //
 228   // Aliternative bethe function assuming a particle type not a mass
 229   // should be slightly faster
 230   //
 231
 232   const Double_t m=AliPID::ParticleMassZ(species);
 233   const Double_t bg=p/m;
 234   Bool_t isNuclei=kFALSE;
 235
 236   //NOTE
 237   //NOTE: if changes are made here, please also check the alternative function above
 238   //NOTE
 239   const Double_t *par=fBBtpcits;
 240   if(isSA){
 241     if(species == AliPID::kElectron){
 242       //if is an electron use a specific BB parameterization
 243       //To be used only between 100 and 160 MeV/c
 244       par=fBBsaElectron;
 245     }else{
 246       par=fBBsa;
 247     }
 248   }else{
 249     if(species == AliPID::kDeuteron) {
 250       par=fBBdeu;
 251       isNuclei=kTRUE;
 252     }
 253     if(species == AliPID::kTriton  ) {
 254       par=fBBtri;
 255       isNuclei=kTRUE;
 256     }
 257   }
 258
 259   return Bethe(bg, par, isNuclei);
 260 }
 261
 262 //_________________________________________________________________________
 263 Double_t AliITSPIDResponse::BetheITSsaHybrid(Double_t p, Double_t mass) const {
 264   //
 265   // returns AliExternalTrackParam::BetheBloch normalized to
 266   // fgMIP at the minimum. The PHOBOS parameterization is used for beta*gamma>0.76.
 267   // For beta*gamma<0.76 a polinomial function is used
 268
 269   Double_t bg=p/mass;
 270   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 271   Double_t gamma=bg/beta;
 272   Double_t bb=1.;
 273
 274   Double_t par[9];
 275   //parameters for pi, K, p
 276   for(Int_t ip=0; ip<9;ip++) par[ip]=fBBsaHybrid[ip];
 277   //if it is an electron the PHOBOS part of the parameterization is tuned for e
 278   //in the range used for identification beta*gamma is >0.76 for electrons
 279   //To be used only between 100 and 160 MeV/c
 280   if(mass>0.0005 && mass<0.00052)for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
 281
 282   if(gamma>=0. && beta>0. && bg>0.1){
 283     if(bg>0.76){//PHOBOS
 284       Double_t eff=1.0;
 285       if(bg<par[2])
 286         eff=(bg-par[3])*(bg-par[3])+par[4];
 287       else
 288         eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 289
 290       bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 291     }else{//Polinomial
 292       bb=par[5] + par[6]/bg + par[7]/(bg*bg) + par[8]/(bg*bg*bg);
 293     }
 294   }
 295   return bb;
 296 }
 297
 298 //_________________________________________________________________________
 299 Double_t AliITSPIDResponse::GetResolution(Double_t bethe,
 300                                           Int_t nPtsForPid,
 301                                          Bool_t isSA,
 302                                          Double_t p,
 303                                          AliPID::EParticleType type) const {
 304   //
 305   // Calculate expected resolution for truncated mean
 306   //
 307   // NEW: Added new variables which are Double_t p and AliPID::EParticleType type
 308   // AliPID::EParticleType type is used to set the correct resolution for the different particles
 309   // default -> AliPID::EParticleType type = AliPID::kPion
 310   // Double_t p is used for the resolution of deuteron and triton, because they are function of the momentum
 311   // default -> Double_t p=0.
 312
 313   Float_t r=0.f;
 314   Double_t c=1.; //this is a correction factor used for the nuclei resolution, while for pion/kaon/proton/electron is 1.
 315
 316   if(isSA) r=fResolSA[nPtsForPid];
 317   else{
 318     const Double_t *par=0x0;
 319     if(type==AliPID::kDeuteron){
 320       if(nPtsForPid==3) par = fResolTPCITSDeu3;
 321       if(nPtsForPid==4) par = fResolTPCITSDeu4;
 322       c=par[2];
 323       r=par[0]+par[1]*p;
 324     } else if(type==AliPID::kTriton){
 325       if(nPtsForPid==3) par = fResolTPCITSTri3;
 326       if(nPtsForPid==4) par = fResolTPCITSTri4;
 327       c=par[2];
 328       r=par[0]+par[1]*p;
 329     } else{
 330       r=fResolTPCITS[nPtsForPid];
 331     }
 332   }
 333
 334   return r*bethe*c;
 335 }
 336
 337
 338 //_________________________________________________________________________
 339 void AliITSPIDResponse::GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES], Bool_t isMC) const {
 340   //
 341   // Method to calculate PID probabilities for a single track
 342   // using the likelihood method
 343   //
 344   const Int_t nLay = 4;
 345   const Int_t nPart= 4;
 346
 347   static AliITSPidParams pars(isMC);  // Pid parametrisation parameters
 348
 349   Double_t itsProb[nPart] = {1,1,1,1}; // e, p, K, pi
 350
 351   for (Int_t iLay = 0; iLay < nLay; iLay++) {
 352     if (qclu[iLay] <= 50.)
 353       continue;
 354
 355     Float_t dedx = qclu[iLay];
 356     Float_t layProb = pars.GetLandauGausNorm(dedx,AliPID::kProton,mom,iLay+3);
 357     itsProb[0] *= layProb;
 358
 359     layProb = pars.GetLandauGausNorm(dedx,AliPID::kKaon,mom,iLay+3);
 360     itsProb[1] *= layProb;
 361
 362     layProb = pars.GetLandauGausNorm(dedx,AliPID::kPion,mom,iLay+3);
 363     itsProb[2] *= layProb;
 364
 365     layProb = pars.GetLandauGausNorm(dedx,AliPID::kElectron,mom,iLay+3);
 366     itsProb[3] *= layProb;
 367   }
 368
 369   // Normalise probabilities
 370   Double_t sumProb = 0;
 371   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 372     sumProb += itsProb[iPart];
 373   }
 374   sumProb += itsProb[2]; // muon cannot be distinguished from pions
 375
 376   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 377     itsProb[iPart]/=sumProb;
 378   }
 379   condprobfun[AliPID::kElectron] = itsProb[3];
 380   condprobfun[AliPID::kMuon] = itsProb[2];
 381   condprobfun[AliPID::kPion] = itsProb[2];
 382   condprobfun[AliPID::kKaon] = itsProb[1];
 383   condprobfun[AliPID::kProton] = itsProb[0];
 384   return;
 385 }
 386
 387 //_________________________________________________________________________
 388 Double_t AliITSPIDResponse::GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType type) const
 389 {
 390   //
 391   // number of sigmas
 392   //
 393   UChar_t clumap=track->GetITSClusterMap();
 394   Int_t nPointsForPid=0;
 395   for(Int_t i=2; i<6; i++){
 396     if(clumap&(1<<i)) ++nPointsForPid;
 397   }
 398   Float_t mom=track->P();
 399
 400   //check for ITS standalone tracks
 401   Bool_t isSA=kTRUE;
 402   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 403
 404   const Float_t dEdx=track->GetITSsignal();
 405
 406   //TODO: in case of the electron, use the SA parametrisation,
 407   //      this needs to be changed if ITS provides a parametrisation
 408   //      for electrons also for ITS+TPC tracks
 409   return GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
 410 }
 411
 412 //_________________________________________________________________________
 413 Double_t AliITSPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
 414 {
 415   //
 416   // Signal - expected
 417   //
 418   const Float_t mom=track->P();
 419   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
 420   Bool_t isSA=kTRUE;
 421   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 422
 423   const Float_t dEdx=track->GetITSsignal();
 424
 425   //TODO: in case of the electron, use the SA parametrisation,
 426   //      this needs to be changed if ITS provides a parametrisation
 427   //      for electrons also for ITS+TPC tracks
 428
 429   const Float_t bethe = Bethe(mom,AliPID::ParticleMassZ(type), isSA || (type==AliPID::kElectron))*chargeFactor;
 430
 431   Double_t delta=-9999.;
 432   if (!ratio) delta=dEdx-bethe;
 433   else if (bethe>1.e-20) delta=dEdx/bethe;
 434
 435   return delta;
 436 }
 437
 438 //_________________________________________________________________________
 439 Int_t AliITSPIDResponse::GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA) const{
 440   // method to get particle identity with simple cuts on dE/dx vs. momentum
 441
 442   Double_t massp=AliPID::ParticleMass(AliPID::kProton);
 443   Double_t massk=AliPID::ParticleMass(AliPID::kKaon);
 444   Double_t bethep=Bethe(mom,massp,isSA);
 445   Double_t bethek=Bethe(mom,massk,isSA);
 446   if(signal>(0.5*(bethep+bethek))) return AliPID::kProton;
 447   Double_t masspi=AliPID::ParticleMass(AliPID::kPion);
 448   Double_t bethepi=Bethe(mom,masspi,isSA);
 449   if(signal>(0.5*(bethepi+bethek))) return AliPID::kKaon;
 450   return AliPID::kPion;
 451
 452 }