91e3ca334ccf404318c45e054dc4ab27ec12350c
[u/mrichter/AliRoot.git] / 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
21 #include <TMath.h>
22
23 ClassImp(AliTRDPartID)
24
25
26 AliTRDPartID::AliTRDPartID()
27 {
28 // create a TRD particle identifier
29
30   fBetheBloch = NULL;
31   fRes = 0.2;
32   fRange = 3.;
33 }
34
35 AliTRDPartID::AliTRDPartID(TF1* betheBloch, Double_t res, Double_t range)
36 {
37 // create a TRD particle identifier with custom settings
38
39   fBetheBloch = betheBloch;
40   fRes = res;
41   fRange = range;
42 }
43
44 AliTRDPartID::~AliTRDPartID()
45 {
46   if (fBetheBloch) delete fBetheBloch;
47 }
48
49
50 Bool_t AliTRDPartID::MakePID(AliESDtrack* track)
51 {
52 // This function calculates the "detector response" PID probabilities 
53
54   static const Double_t masses[]={
55     0.000511, 0.105658, 0.139570, 0.493677, 0.938272
56   };
57
58   if (((track->GetStatus()&AliESDtrack::kTRDin) == 0) &&
59       ((track->GetStatus()&AliESDtrack::kTRDout) == 0)) return kFALSE;
60   Double_t momentum = track->GetP();
61   if (momentum < 0.001) return kFALSE;
62
63   // get the probability densities
64   Double_t pSum = 0;
65   for (Int_t iSpecies = 0; iSpecies < AliESDtrack::kSPECIES; iSpecies++) {
66     Double_t expectedSignal = fBetheBloch->Eval(momentum/masses[iSpecies]);
67     Double_t expectedError = fRes * expectedSignal;
68     Double_t measuredSignal = track->GetTRDsignal();
69     if (TMath::Abs(measuredSignal - expectedSignal) > fRange * expectedError) {
70       track->SetTRDpid(iSpecies, 0.);
71     } else {
72       Double_t p = TMath::Gaus(measuredSignal, expectedSignal, expectedError);
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 < AliESDtrack::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, 
94                         0.001, 100000., 5);
95   fBetheBloch->SetParameters(1, 10, 0.00002, 2, 2);
96   fBetheBloch->SetParLimits(2, 0., 0.01);
97   fBetheBloch->SetParLimits(3, 0., 10.);
98   fBetheBloch->SetParLimits(4, 0., 10.);
99   fBetheBloch->SetFillStyle(0);
100   fBetheBloch->SetLineColor(kRed);
101   dEdxVsBetaGamma->Fit(fBetheBloch, "N");
102 }
103
104 TF1* AliTRDPartID::CreateBetheBloch(Double_t mass)
105 {
106 // create a function for expected dE/dx vs p
107
108   TF1* result = new TF1("betheBlochMass", fcnBetheBlochMass, 
109                         0.001, 100000., 6);
110   result->SetParameter(0, mass);
111   if (fBetheBloch) {
112     for (Int_t iPar = 0; iPar < 5; iPar++) {
113       result->SetParameter(iPar+1, fBetheBloch->GetParameter(iPar));
114       result->SetParError(iPar+1, fBetheBloch->GetParError(iPar));
115     }
116   }
117   result->SetFillStyle(0);
118   return result;
119 }
120
121 Double_t AliTRDPartID::fcnBetheBloch(Double_t* xx, Double_t* par)
122 {
123 // parametrized bethe bloch function (xx[0] = beta*gamma = p/m):
124 //
125 //   p0/beta^p3 * [ p1 - log(p2 + 1/(beta*gamma)^p4) - beta^p3 ]
126 //
127
128   Double_t betaGamma2 = xx[0] * xx[0];
129   Double_t beta2 = betaGamma2 / (1. + betaGamma2);
130   Double_t betaPar3 = TMath::Power(beta2, par[3]/2.);
131   return par[0]/betaPar3 * (par[1] - TMath::Log(TMath::Abs(par[2] + TMath::Power(betaGamma2, -par[4]/2.))) - betaPar3);
132 }
133
134 Double_t AliTRDPartID::fcnBetheBlochMass(Double_t* xx, Double_t* par)
135 {
136 // parametrized bethe bloch function:  xx[0] = p,  p0 = mass
137
138   Double_t betaGamma = xx[0] / par[0];
139   return fcnBetheBloch(&betaGamma, &par[1]);
140 }