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