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     fBBsa[0]=2.73198E7; //pure PHOBOS parameterization
  47     fBBsa[1]=6.92389;
  48     fBBsa[2]=1.90088E-6;
  49     fBBsa[3]=1.90088E-6;
  50     fBBsa[4]=3.40644E-7;
  51     fBBsaHybrid[0]=1.43505E7;  //PHOBOS+Polinomial parameterization
  52     fBBsaHybrid[1]=49.3402;
  53     fBBsaHybrid[2]=1.77741E-7;
  54     fBBsaHybrid[3]=1.77741E-7;
  55     fBBsaHybrid[4]=1.01311E-7;
  56     fBBsaHybrid[5]=77.2777;
  57     fBBsaHybrid[6]=33.4099;
  58     fBBsaHybrid[7]=46.0089;
  59     fBBsaHybrid[8]=-2.26583;
  60     fBBsaElectron[0]=4.05799E6;  //electrons in the ITS
  61     fBBsaElectron[1]=38.5713;
  62     fBBsaElectron[2]=1.46462E-7;
  63     fBBsaElectron[3]=1.46462E-7;
  64     fBBsaElectron[4]=4.40284E-7;
  65     fResolSA[0]=1.;   // 0 cluster tracks should not be used
  66     fResolSA[1]=0.25;  // rough values for tracks with 1
  67     fResolSA[2]=0.131;   // value from pp 2010 run (L. Milano, 16-Jun-11)
  68     fResolSA[3]=0.113; // value from pp 2010 run
  69     fResolSA[4]=0.104; // value from pp 2010 run
  70     for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
  71   }else{
  72     fBBtpcits[0]=1.04;
  73     fBBtpcits[1]=27.14;
  74     fBBtpcits[2]=1.00;
  75     fBBtpcits[3]=0.964;
  76     fBBtpcits[4]=2.59;
  77     fBBsa[0]=2.02078E7; //pure PHOBOS parameterization
  78     fBBsa[1]=14.0724;
  79     fBBsa[2]=3.84454E-7;
  80     fBBsa[3]=3.84454E-7;
  81     fBBsa[4]=2.43913E-7;
  82     fBBsaHybrid[0]=1.05381E7; //PHOBOS+Polinomial parameterization
  83     fBBsaHybrid[1]=89.3933;
  84     fBBsaHybrid[2]=2.4831E-7;
  85     fBBsaHybrid[3]=2.4831E-7;
  86     fBBsaHybrid[4]=7.80591E-8;
  87     fBBsaHybrid[5]=62.9214;
  88     fBBsaHybrid[6]=32.347;
  89     fBBsaHybrid[7]=58.7661;
  90     fBBsaHybrid[8]=-3.39869;
  91     fBBsaElectron[0]=2.26807E6; //electrons in the ITS
  92     fBBsaElectron[1]=99.985;
  93     fBBsaElectron[2]=0.000714841;
  94     fBBsaElectron[3]=0.000259585;
  95     fBBsaElectron[4]=1.39412E-7;
  96     fResolSA[0]=1.;   // 0 cluster tracks should not be used
  97     fResolSA[1]=0.25;  // rough values for tracks with 1
  98     fResolSA[2]=0.126;   // value from pp 2010 simulations (L. Milano, 16-Jun-11)
  99     fResolSA[3]=0.109; // value from pp 2010 simulations
 100     fResolSA[4]=0.097; // value from pp 2010 simulations
 101     for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
 102   }
 103 }
 104
 105 /*
 106 //_________________________________________________________________________
 107 AliITSPIDResponse::AliITSPIDResponse(Double_t *param):
 108   fRes(param[0]),
 109   fKp1(15.77),
 110   fKp2(4.95),
 111   fKp3(0.312),
 112   fKp4(2.14),
 113   fKp5(0.82)
 114 {
 115   //
 116   //  The main constructor
 117   //
 118   for (Int_t i=0; i<5;i++) {
 119       fBBsa[i]=0.;
 120       fBBtpcits[i]=0.;
 121       fResolSA[i]=0.;
 122       fResolTPCITS[i]=0.;
 123   }
 124 }
 125 */
 126
 127 //_________________________________________________________________________
 128 Double_t AliITSPIDResponse::BetheAleph(Double_t p, Double_t mass) const {
 129   //
 130   // returns AliExternalTrackParam::BetheBloch normalized to
 131   // fgMIP at the minimum
 132   //
 133
 134   Double_t bb=
 135     AliExternalTrackParam::BetheBlochAleph(p/mass,fKp1,fKp2,fKp3,fKp4,fKp5);
 136   return bb;
 137 }
 138
 139 //_________________________________________________________________________
 140 Double_t AliITSPIDResponse::Bethe(Double_t p, Double_t mass, Bool_t isSA) const {
 141   //
 142   // returns AliExternalTrackParam::BetheBloch normalized to
 143   // fgMIP at the minimum
 144   //
 145
 146   Double_t bg=p/mass;
 147   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 148   Double_t gamma=bg/beta;
 149   Double_t bb=1.;
 150
 151   Double_t par[5];
 152   if(isSA){
 153     if(mass>0.0005 && mass<0.00052){
 154       //if is an electron use a specific BB parameterization
 155       //To be used only between 100 and 160 MeV/c
 156       for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
 157     }else{
 158       for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsa[ip];
 159     }
 160   }else{
 161     for(Int_t ip=0; ip<5;ip++) par[ip]=fBBtpcits[ip];
 162   }
 163   Double_t eff=1.0;
 164   if(bg<par[2])
 165     eff=(bg-par[3])*(bg-par[3])+par[4];
 166   else
 167     eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 168
 169   if(gamma>=0. && beta>0.){
 170     bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 171   }
 172   return bb;
 173 }
 174
 175 //_________________________________________________________________________
 176 Double_t AliITSPIDResponse::BetheITSsaHybrid(Double_t p, Double_t mass) const {
 177   //
 178   // returns AliExternalTrackParam::BetheBloch normalized to
 179   // fgMIP at the minimum. The PHOBOS parameterization is used for beta*gamma>0.76.
 180   // For beta*gamma<0.76 a polinomial function is used
 181
 182   Double_t bg=p/mass;
 183   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 184   Double_t gamma=bg/beta;
 185   Double_t bb=1.;
 186
 187   Double_t par[9];
 188   //parameters for pi, K, p
 189   for(Int_t ip=0; ip<9;ip++) par[ip]=fBBsaHybrid[ip];
 190   //if it is an electron the PHOBOS part of the parameterization is tuned for e
 191   //in the range used for identification beta*gamma is >0.76 for electrons
 192   //To be used only between 100 and 160 MeV/c
 193   if(mass>0.0005 && mass<0.00052)for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
 194
 195   if(gamma>=0. && beta>0. && bg>0.1){
 196     if(bg>0.76){//PHOBOS
 197       Double_t eff=1.0;
 198       if(bg<par[2])
 199         eff=(bg-par[3])*(bg-par[3])+par[4];
 200       else
 201         eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 202
 203       bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 204     }else{//Polinomial
 205       bb=par[5] + par[6]/bg + par[7]/(bg*bg) + par[8]/(bg*bg*bg);
 206     }
 207   }
 208   return bb;
 209 }
 210
 211 //_________________________________________________________________________
 212 Double_t AliITSPIDResponse::GetResolution(Double_t bethe,
 213                                           Int_t nPtsForPid,
 214                                           Bool_t isSA) const {
 215   //
 216   // Calculate expected resolution for truncated mean
 217   //
 218   Float_t r;
 219   if(isSA) r=fResolSA[nPtsForPid];
 220   else r=fResolTPCITS[nPtsForPid];
 221   return r*bethe;
 222 }
 223
 224
 225
 226
 227 //_________________________________________________________________________
 228 void AliITSPIDResponse::GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES], Bool_t isMC) const {
 229   //
 230   // Method to calculate PID probabilities for a single track
 231   // using the likelihood method
 232   //
 233   const Int_t nLay = 4;
 234   const Int_t nPart= 4;
 235
 236   static AliITSPidParams pars(isMC);  // Pid parametrisation parameters
 237
 238   Double_t itsProb[nPart] = {1,1,1,1}; // e, p, K, pi
 239
 240   for (Int_t iLay = 0; iLay < nLay; iLay++) {
 241     if (qclu[iLay] <= 50.)
 242       continue;
 243
 244     Float_t dedx = qclu[iLay];
 245     Float_t layProb = pars.GetLandauGausNorm(dedx,AliPID::kProton,mom,iLay+3);
 246     itsProb[0] *= layProb;
 247
 248     layProb = pars.GetLandauGausNorm(dedx,AliPID::kKaon,mom,iLay+3);
 249     itsProb[1] *= layProb;
 250
 251     layProb = pars.GetLandauGausNorm(dedx,AliPID::kPion,mom,iLay+3);
 252     itsProb[2] *= layProb;
 253
 254     layProb = pars.GetLandauGausNorm(dedx,AliPID::kElectron,mom,iLay+3);
 255     itsProb[3] *= layProb;
 256   }
 257
 258   // Normalise probabilities
 259   Double_t sumProb = 0;
 260   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 261     sumProb += itsProb[iPart];
 262   }
 263   sumProb += itsProb[2]; // muon cannot be distinguished from pions
 264
 265   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 266     itsProb[iPart]/=sumProb;
 267   }
 268   condprobfun[AliPID::kElectron] = itsProb[3];
 269   condprobfun[AliPID::kMuon] = itsProb[2];
 270   condprobfun[AliPID::kPion] = itsProb[2];
 271   condprobfun[AliPID::kKaon] = itsProb[1];
 272   condprobfun[AliPID::kProton] = itsProb[0];
 273   return;
 274 }
 275
 276 //_________________________________________________________________________
 277 Double_t AliITSPIDResponse::GetNumberOfSigmas( const AliVTrack* track, AliPID::EParticleType type) const
 278 {
 279   //
 280   // number of sigmas
 281   //
 282   UChar_t clumap=track->GetITSClusterMap();
 283   Int_t nPointsForPid=0;
 284   for(Int_t i=2; i<6; i++){
 285     if(clumap&(1<<i)) ++nPointsForPid;
 286   }
 287   Float_t mom=track->P();
 288
 289   //check for ITS standalone tracks
 290   Bool_t isSA=kTRUE;
 291   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 292
 293   const Float_t dEdx=track->GetITSsignal();
 294
 295   //TODO: in case of the electron, use the SA parametrisation,
 296   //      this needs to be changed if ITS provides a parametrisation
 297   //      for electrons also for ITS+TPC tracks
 298   return GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
 299 }
 300
 301 //_________________________________________________________________________
 302 Double_t AliITSPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
 303 {
 304   //
 305   // Signal - expected
 306   //
 307   const Float_t mom=track->P();
 308   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
 309   Bool_t isSA=kTRUE;
 310   if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
 311
 312   const Float_t dEdx=track->GetITSsignal();
 313
 314   //TODO: in case of the electron, use the SA parametrisation,
 315   //      this needs to be changed if ITS provides a parametrisation
 316   //      for electrons also for ITS+TPC tracks
 317
 318   const Float_t bethe = Bethe(mom,AliPID::ParticleMassZ(type), isSA || (type==AliPID::kElectron))*chargeFactor;
 319
 320   Double_t delta=-9999.;
 321   if (!ratio) delta=dEdx-bethe;
 322   else if (bethe>1.e-20) delta=dEdx/bethe;
 323
 324   return delta;
 325 }
 326
 327 //_________________________________________________________________________
 328 Int_t AliITSPIDResponse::GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA) const{
 329   // method to get particle identity with simple cuts on dE/dx vs. momentum
 330
 331   Double_t massp=AliPID::ParticleMass(AliPID::kProton);
 332   Double_t massk=AliPID::ParticleMass(AliPID::kKaon);
 333   Double_t bethep=Bethe(mom,massp,isSA);
 334   Double_t bethek=Bethe(mom,massk,isSA);
 335   if(signal>(0.5*(bethep+bethek))) return AliPID::kProton;
 336   Double_t masspi=AliPID::ParticleMass(AliPID::kPion);
 337   Double_t bethepi=Bethe(mom,masspi,isSA);
 338   if(signal>(0.5*(bethepi+bethek))) return AliPID::kKaon;
 339   return AliPID::kPion;
 340
 341 }