b0bf074436fa49fa6be317e74dfdfcb7e9501fd8
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.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 <TGraph.h>
26 #include <TObjArray.h>
27 #include <TSpline.h>
28 #include <TMath.h>
29
30 #include "AliExternalTrackParam.h"
31
32 #include "AliTPCPIDResponse.h"
33
34 ClassImp(AliTPCPIDResponse)
35
36 //_________________________________________________________________________
37 AliTPCPIDResponse::AliTPCPIDResponse():
38   fMIP(50.),
39   fRes0(0.07),
40   fResN2(0.),
41   fKp1(0.0283086),
42   fKp2(2.63394e+01),
43   fKp3(5.04114e-11),
44   fKp4(2.12543),
45   fKp5(4.88663),
46   fUseDatabase(kFALSE),
47   fResponseFunctions(AliPID::kUnknown+1)
48 {
49   //
50   //  The default constructor
51   //
52 }
53
54 //_________________________________________________________________________
55 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
56   fMIP(param[0]),
57   fRes0(param[1]),
58   fResN2(param[2]),
59   fKp1(0.0283086),
60   fKp2(2.63394e+01),
61   fKp3(5.04114e-11),
62   fKp4(2.12543),
63   fKp5(4.88663),
64   fUseDatabase(kFALSE),
65   fResponseFunctions(AliPID::kUnknown+1)
66 {
67   //
68   //  The main constructor
69   //
70 }
71
72 //_________________________________________________________________________
73 Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
74   //
75   // This is the Bethe-Bloch function normalised to 1 at the minimum
76   // WARNING
77   // Simulated and reconstructed Bethe-Bloch differs
78   //           Simulated  curve is the dNprim/dx
79   //           Reconstructed is proportianal dNtot/dx
80   // Temporary fix for production -  Simple linear correction function
81   // Future    2 Bethe Bloch formulas needed
82   //           1. for simulation
83   //           2. for reconstructed PID
84   //
85   
86 //   const Float_t kmeanCorrection =0.1;
87   Double_t bb=
88     AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
89   return bb*fMIP;
90 }
91
92 //_________________________________________________________________________
93 void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
94                              Double_t kp2,
95                              Double_t kp3,
96                              Double_t kp4,
97                              Double_t kp5) {
98   //
99   // Set the parameters of the ALEPH Bethe-Bloch formula
100   //
101   fKp1=kp1;
102   fKp2=kp2;
103   fKp3=kp3;
104   fKp4=kp4;
105   fKp5=kp5;
106 }
107 //_________________________________________________________________________
108 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
109   //
110   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
111   //
112   fRes0 = res0;
113   fResN2 = resN2;
114 }
115
116 //_________________________________________________________________________
117 Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
118                                               AliPID::EParticleType n) const {
119   //
120   // Calculates the expected PID signal as the function of 
121   // the information stored in the track, for the specified particle type 
122   //  
123   // At the moment, these signals are just the results of calling the 
124   // Bethe-Bloch formula. 
125   // This can be improved. By taking into account the number of
126   // assigned clusters and/or the track dip angle, for example.  
127   //
128   
129   Double_t mass=AliPID::ParticleMassZ(n);
130   if (!fUseDatabase||fResponseFunctions.GetEntriesFast()>AliPID::kUnknown) return Bethe(mom/mass);
131   //
132   TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
133
134   if (!responseFunction) return Bethe(mom/mass);
135   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
136   // !!! Splines for light nuclei need to be normalised to this factor !!!
137   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
138   return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
139
140 }
141
142 //_________________________________________________________________________
143 Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom, 
144                                              const Int_t nPoints,
145                                          AliPID::EParticleType n) const {
146   //
147   // Calculates the expected sigma of the PID signal as the function of 
148   // the information stored in the track, for the specified particle type 
149   //  
150   //
151   
152   if (nPoints != 0) 
153     return GetExpectedSignal(mom,n)*fRes0*sqrt(1. + fResN2/nPoints);
154   else
155     return GetExpectedSignal(mom,n)*fRes0;
156 }