STEER/AliTPCpidESD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //-----------------------------------------------------------------
  17 //           Implementation of the TPC PID class
  18 // Very naive one... Should be made better by the detector experts...
  19 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  20 // With many additions and modifications suggested by
  21 //      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
  22 //      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
  23 //-----------------------------------------------------------------
  24
  25 #include "AliTPCpidESD.h"
  26 #include "AliESDEvent.h"
  27 #include "AliESDtrack.h"
  28
  29 ClassImp(AliTPCpidESD)
  30
  31 //_________________________________________________________________________
  32 AliTPCpidESD::AliTPCpidESD():
  33     fMIP(50.),
  34     fRes(0.07),
  35     fRange(5.),
  36     fKp1(0.76176e-1),
  37     fKp2(10.632),
  38     fKp3(0.13279e-4),
  39     fKp4(1.8631),
  40     fKp5(1.9479)
  41 {
  42   //
  43   //  The default constructor
  44   //
  45 }
  46
  47 //_________________________________________________________________________
  48 AliTPCpidESD::AliTPCpidESD(Double_t *param):
  49     fMIP(param[0]),
  50     fRes(param[1]),
  51     fRange(param[2]),
  52     fKp1(0.76176e-1),
  53     fKp2(10.632),
  54     fKp3(0.13279e-4),
  55     fKp4(1.8631),
  56     fKp5(1.9479)
  57 {
  58   //
  59   //  The main constructor
  60   //
  61 }
  62
  63 Double_t AliTPCpidESD::Bethe(Double_t betaGamma) const {
  64   //
  65   // This is the Bethe-Bloch function normalised to 1 at the minimum
  66   // WARNING
  67   // Simulated and reconstructed Bethe-Bloch differs
  68   //           Simulated  curve is the dNprim/dx
  69   //           Reconstructed is proportianal dNtot/dx
  70   // Temporary fix for production -  Simple linear correction function
  71   // Future    2 Bethe Bloch formulas needed
  72   //           1. for simulation
  73   //           2. for reconstructed PID
  74   //
  75   const Float_t kmeanCorrection =0.1;
  76   Double_t bb=
  77     AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
  78   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
  79   bb *= meanCorrection;
  80   return bb;
  81 }
  82
  83 //_________________________________________________________________________
  84 void AliTPCpidESD::SetBetheBlochParameters(Double_t kp1,
  85                              Double_t kp2,
  86                              Double_t kp3,
  87                              Double_t kp4,
  88                              Double_t kp5) {
  89   //
  90   // Set the parameters of the ALEPH Bethe-Bloch formula
  91   //
  92   fKp1=kp1;
  93   fKp2=kp2;
  94   fKp3=kp3;
  95   fKp4=kp4;
  96   fKp5=kp5;
  97 }
  98
  99 //_________________________________________________________________________
 100 Bool_t AliTPCpidESD::ExpectedSigmas(const AliESDtrack *t,
 101                                      Double_t s[],
 102                                      Int_t n) const {
 103   //
 104   // Calculate the expected dE/dx resolution as the function of
 105   // the information stored in the track.
 106   //
 107   // At the moment, this resolution is just proportional to the expected
 108   // signal. This can be improved. By taking into account the number of
 109   // assigned clusters, for example.
 110   //
 111   Bool_t ok=kFALSE;
 112   Double_t signals[AliPID::kSPECIESN];
 113   if (ExpectedSignals(t,signals,n)) {
 114      for (Int_t i=0; i<n; i++) s[i] = fRes*signals[i];
 115      ok=kTRUE;
 116   }
 117   return ok;
 118 }
 119
 120 //_________________________________________________________________________
 121 Bool_t AliTPCpidESD::NumberOfSigmas(const AliESDtrack *t,
 122                                      Double_t s[],
 123                                      Int_t n) const {
 124   //
 125   // Calculate the deviation of the actual PID signal from the expected
 126   // signal, in units of expected sigmas.
 127   //
 128   Bool_t ok=kFALSE;
 129
 130   Double_t dedx=t->GetTPCsignal()/fMIP;
 131   Double_t sigmas[AliPID::kSPECIESN];
 132   if (ExpectedSigmas(t,sigmas,n)) {
 133      Double_t signals[AliPID::kSPECIESN];
 134      if (ExpectedSignals(t,signals,n)) {
 135        for (Int_t i=0; i<n; i++) s[i] = (signals[i] - dedx)/sigmas[i];
 136        ok=kTRUE;
 137      }
 138   }
 139   return ok;
 140 }
 141
 142 //_________________________________________________________________________
 143 Bool_t AliTPCpidESD::ExpectedSignals(const AliESDtrack *t,
 144                                       Double_t s[],
 145                                       Int_t n) const {
 146   //
 147   // Calculates the expected PID signals as the function of
 148   // the information stored in the track.
 149   //
 150   // At the moment, these signals are just the results of calling the
 151   // Bethe-Bloch formula.
 152   // This can be improved. By taking into account the number of
 153   // assigned clusters and/or the track dip angle, for example.
 154   //
 155
 156   Double_t mom=t->GetP();
 157   const AliExternalTrackParam *in=t->GetInnerParam();
 158   if (in) mom=in->GetP();
 159
 160   for (Int_t i=0; i<n; i++) {
 161       Double_t mass=AliPID::ParticleMass(i);
 162       s[i]=Bethe(mom/mass);
 163   }
 164
 165   return kTRUE;
 166 }
 167
 168 //_________________________________________________________________________
 169 Double_t AliTPCpidESD::GetExpectedSignal(const AliESDtrack *t,
 170                                          AliPID::EParticleType n) const {
 171   //
 172   // Calculates the expected PID signal as the function of
 173   // the information stored in the track, for the specified particle type
 174   //
 175   // At the moment, these signals are just the results of calling the
 176   // Bethe-Bloch formula.
 177   // This can be improved. By taking into account the number of
 178   // assigned clusters and/or the track dip angle, for example.
 179   //
 180
 181   Double_t mom=t->GetP();
 182   const AliExternalTrackParam *in=t->GetInnerParam();
 183   if (in) mom=in->GetP();
 184
 185   Double_t mass=AliPID::ParticleMass(n);
 186   return Bethe(mom/mass);
 187 }
 188
 189 //_________________________________________________________________________
 190 Double_t AliTPCpidESD::GetExpectedSigma(const AliESDtrack *t,
 191                                          AliPID::EParticleType n) const {
 192   //
 193   // Calculates the expected sigma of the PID signal as the function of
 194   // the information stored in the track, for the specified particle type
 195   //
 196   //
 197   // At the moment, this sigma is just proportional to the expected
 198   // signal. This can be improved. By taking into account the number of
 199   // assigned clusters, for example.
 200   //
 201
 202   return fRes*GetExpectedSignal(t,n);
 203 }
 204
 205 //_________________________________________________________________________
 206 Double_t AliTPCpidESD::GetNumberOfSigmas(const AliESDtrack *t,
 207                                          AliPID::EParticleType n) const {
 208   //
 209   // Calculate the deviation of the actual PID signal from the expected
 210   // signal, in units of expected sigmas, for the specified particle type
 211   //
 212
 213     Double_t dedx=t->GetTPCsignal()/fMIP;
 214     return (dedx - GetExpectedSignal(t,n))/GetExpectedSigma(t,n);
 215 }
 216
 217 //_________________________________________________________________________
 218 Int_t AliTPCpidESD::MakePID(AliESDEvent *event)
 219 {
 220   //
 221   //  This function calculates the "detector response" PID probabilities
 222   //
 223   Int_t ntrk=event->GetNumberOfTracks();
 224   for (Int_t i=0; i<ntrk; i++) {
 225     AliESDtrack *t=event->GetTrack(i);
 226     if ((t->GetStatus()&AliESDtrack::kTPCin )==0)
 227       if ((t->GetStatus()&AliESDtrack::kTPCout)==0) continue;
 228     Double_t p[10];
 229     Double_t dedx=t->GetTPCsignal()/fMIP;
 230     Bool_t mismatch=kTRUE, heavy=kTRUE;
 231     for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 232       AliPID::EParticleType type=AliPID::EParticleType(j);
 233       Double_t bethe=GetExpectedSignal(t,type);
 234       Double_t sigma=GetExpectedSigma(t,type);
 235       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 236         p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 237       } else {
 238         p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 239         mismatch=kFALSE;
 240       }
 241
 242       // Check for particles heavier than (AliPID::kSPECIES - 1)
 243       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
 244
 245     }
 246
 247     if (mismatch)
 248        for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
 249
 250     t->SetTPCpid(p);
 251
 252     if (heavy) t->ResetStatus(AliESDtrack::kTPCpid);
 253
 254   }
 255   return 0;
 256 }