TRD/AliTRDPartID.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 #include "AliTRDPartID.h"
  17 #include "AliESDtrack.h"
  18 #include <TProfile.h>
  19 #include <TF1.h>
  20 #include <TFile.h>
  21 #include <TROOT.h>
  22 #include <TMath.h>
  23 #include <Riostream.h>
  24
  25 ClassImp(AliTRDPartID)
  26
  27
  28 AliTRDPartID::AliTRDPartID()
  29 {
  30 // create a TRD particle identifier
  31
  32   fBetheBloch = NULL;
  33   fRes = 0.2;
  34   fRange = 3.;
  35 }
  36
  37 AliTRDPartID::AliTRDPartID(TF1* betheBloch, Double_t res, Double_t range)
  38 {
  39 // create a TRD particle identifier with custom settings
  40
  41   fBetheBloch = betheBloch;
  42   fRes = res;
  43   fRange = range;
  44 }
  45
  46 AliTRDPartID::~AliTRDPartID()
  47 {
  48   if (fBetheBloch) delete fBetheBloch;
  49 }
  50
  51
  52 Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
  53 {
  54 // This function calculates the "detector response" PID probabilities
  55
  56   if (((track->GetStatus()&AliESDtrack::kTRDin) == 0) &&
  57       ((track->GetStatus()&AliESDtrack::kTRDout) == 0)) return kFALSE;
  58   Double_t momentum = track->GetP();
  59   if (momentum < 0.001) return kFALSE;
  60
  61   // get the probability densities
  62   Double_t pSum = 0;
  63   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
  64     Double_t expectedSignal = fBetheBloch->Eval(momentum/AliPID::ParticleMass(iSpecies));
  65     Double_t expectedError = fRes * expectedSignal;
  66     Double_t measuredSignal = track->GetTRDsignal();
  67     if (TMath::Abs(measuredSignal - expectedSignal) > fRange * expectedError) {
  68       track->SetTRDpid(iSpecies, 0.);
  69     } else {
  70       Double_t delta = (measuredSignal-expectedSignal) / expectedError;
  71       const Double_t kInvSqr2Pi = 0.398942280401432703;
  72       Double_t p = kInvSqr2Pi / expectedError * TMath::Exp(-delta*delta / 2.);
  73       pSum += p;
  74       track->SetTRDpid(iSpecies, p);
  75     }
  76   }
  77
  78   // "normalize" the probability densities
  79   if (pSum <= 0) return kFALSE;
  80   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
  81     track->SetTRDpid(iSpecies, track->GetTRDpid(iSpecies) / pSum);
  82   }
  83
  84   return kTRUE;
  85 }
  86
  87
  88 void AliTRDPartID::FitBetheBloch(TProfile* dEdxVsBetaGamma)
  89 {
  90 // determine the parameters of the bethe bloch function
  91
  92   if (fBetheBloch) delete fBetheBloch;
  93   fBetheBloch = new TF1("fBetheBlochTRD", fcnBetheBloch, 0.001, 100000., 5);
  94   fBetheBloch->SetParameters(1, 10, 0.00002, 2, 2);
  95   fBetheBloch->SetParLimits(2, 0., 0.01);
  96   fBetheBloch->SetParLimits(3, 0., 10.);
  97   fBetheBloch->SetParLimits(4, 0., 10.);
  98   fBetheBloch->SetFillStyle(0);
  99   fBetheBloch->SetLineColor(kRed);
 100   dEdxVsBetaGamma->Fit(fBetheBloch, "NIR", "goff", 0.6, dEdxVsBetaGamma->GetXaxis()->GetXmax());
 101 }
 102
 103 TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
 104 {
 105 // create a function for expected dE/dx vs p
 106
 107   TF1* result = new TF1("betheBlochMass", fcnBetheBlochMass,
 108                         0.001, 100000., 6);
 109   result->SetParameter(0, mass);
 110   if (fBetheBloch) {
 111     for (Int_t iPar = 0; iPar < 5; iPar++) {
 112       result->SetParameter(iPar+1, fBetheBloch->GetParameter(iPar));
 113       result->SetParError(iPar+1, fBetheBloch->GetParError(iPar));
 114     }
 115   }
 116   result->SetFillStyle(0);
 117   return result;
 118 }
 119
 120 AliTRDPartID* AliTRDPartID::GetFromFile(const char* fileName)
 121 {
 122 // read an AliTRDPartID object from a file
 123
 124   TFile* pidFile = (TFile*) gROOT->GetListOfFiles()->FindObject(fileName);
 125   Bool_t fileOpened = kFALSE;
 126   if (!pidFile) {
 127     pidFile = TFile::Open(fileName);
 128     fileOpened = kTRUE;
 129   }
 130   if (!pidFile->IsOpen()) {
 131     cerr << "Can't open " << fileName << " !\n";
 132     if (fileOpened) delete pidFile;
 133     return NULL;
 134   }
 135   gROOT->cd();
 136
 137   AliTRDPartID* trdPID = (AliTRDPartID*) pidFile->Get("AliTRDPartID")->Clone();
 138   if (!trdPID) {
 139     cerr << "Can't get PID object !\n";
 140   } else {
 141     trdPID->GetBetheBloch()->SetFunction(fcnBetheBloch);
 142   }
 143
 144   if (fileOpened) {
 145     pidFile->Close();
 146     delete pidFile;
 147   }
 148
 149   return trdPID;
 150 }
 151
 152 Double_t AliTRDPartID::fcnBetheBloch(Double_t* xx, Double_t* par)
 153 {
 154 // parametrized bethe bloch function (xx[0] = beta*gamma = p/m):
 155 //
 156 //   p0/beta^p3 * [ p1 - log(p2 + 1/(beta*gamma)^p4) - beta^p3 ]
 157 //
 158
 159   Double_t betaGamma2 = xx[0] * xx[0];
 160   Double_t beta2 = betaGamma2 / (1. + betaGamma2);
 161   Double_t betaPar3 = TMath::Power(beta2, par[3]/2.);
 162   return par[0]/betaPar3 * (par[1] - TMath::Log(TMath::Abs(par[2] + TMath::Power(betaGamma2, -par[4]/2.))) - betaPar3);
 163 }
 164
 165 Double_t AliTRDPartID::fcnBetheBlochMass(Double_t* xx, Double_t* par)
 166 {
 167 // parametrized bethe bloch function:  xx[0] = p,  p0 = mass
 168
 169   Double_t betaGamma = xx[0] / par[0];
 170   return fcnBetheBloch(&betaGamma, &par[1]);
 171 }