Fixed a bug on retrieving TPC inner param
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnCutBetheBloch.cxx
1 //
2 // Class AliRsnCutBetheBloch
3 //
4 // General implementation of a single cut strategy, which can be:
5 // - a value contained in a given interval  [--> IsBetween()   ]
6 // - a value equal to a given reference     [--> MatchesValue()]
7 //
8 // In all cases, the reference value(s) is (are) given as data members
9 // and each kind of cut requires a given value type (Int, UInt, Double),
10 // but the cut check procedure is then automatized and chosen thanks to
11 // an enumeration of the implemented cut types.
12 // At the end, the user (or any other point which uses this object) has
13 // to use the method IsSelected() to check if this cut has been passed.
14 //
15 // authors: Martin Vala (martin.vala@cern.ch)
16 //          Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
17 //
18 #include "TMath.h"
19
20 #include "AliExternalTrackParam.h"
21
22 #include "AliRsnDaughter.h"
23 #include "AliRsnCutBetheBloch.h"
24
25 ClassImp(AliRsnCutBetheBloch)
26
27 //_________________________________________________________________________________________________
28 AliRsnCutBetheBloch::AliRsnCutBetheBloch() :
29     AliRsnCut(),
30     fCorrect(kTRUE),
31     fMIP(50.0),
32     fType(AliPID::kUnknown)
33 {
34 //
35 // Default constructor.
36 //
37
38   fConst[0] = fConst[1] = fConst[2] = fConst[3] = fConst[4] = 0.0;
39 }
40
41 //_________________________________________________________________________________________________
42 AliRsnCutBetheBloch::AliRsnCutBetheBloch
43 (const char *name, Double_t fractionRange, AliPID::EParticleType type, Double_t mip, Bool_t correct) :
44     AliRsnCut(name, 0.0, fractionRange),
45     fCorrect(correct),
46     fMIP(mip),
47     fType(type)
48 {
49 //
50 // Main constructor.
51 // the cut range is the relative fraction of the value:
52 // BB*(1-fraction) < TPC < BB*(1+fraction)
53 // which means:
54 // -fraction < (TPC - BB)/BB < fraction
55 //
56
57   fConst[0] = fConst[1] = fConst[2] = fConst[3] = fConst[4] = 0.0;
58 }
59
60 //_____________________________________________________________________________
61 Double_t AliRsnCutBetheBloch::BetheBloch(AliRsnDaughter * const trackRef)
62 {
63 //
64 // Computes the theoretical dE/dx according to
65 // a given mass hypothesis, from which betaGamma is computed
66 //
67 // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
68 // It is normalized to 1 at the minimum.
69 //
70 // The default values for the kp* parameters are for ALICE TPC.
71 // The value is computed in MIP units, multiplied by 50 to have it in energy.
72 //
73
74   AliPID pid;
75   Double_t mass = pid.ParticleMass(fType);
76
77   // get the track momentum at the inner wall of TPC: if absent cut is not passed
78   if (!trackRef->GetRefESD()->GetInnerParam()) return 1000000.0;
79   AliExternalTrackParam track(*trackRef->GetRefESD()->GetInnerParam());
80
81   Double_t betaGamma = track.P() / mass;
82   Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);
83   Double_t aa = TMath::Power(beta, fConst[3]);
84   Double_t bb = TMath::Power(1.0 / betaGamma, fConst[4]);
85
86   bb = TMath::Log(fConst[2] + bb);
87
88   Double_t out = (fConst[1] - aa - bb) * fConst[0] / aa;
89
90   if (fCorrect) {
91     Double_t kMeanCorr = 0.1;
92     Double_t meanCorr = (1 + (out - 1) * kMeanCorr);
93     out *= meanCorr;
94   }
95
96   return out * fMIP;
97 }
98
99 //_____________________________________________________________________________
100 Double_t AliRsnCutBetheBloch::RelDiff(AliRsnDaughter *track)
101 {
102 //
103 // Relative difference between BB value and TPC signal
104 //
105
106   if (!track->GetRefESD()) return -99999.9;
107
108   // compute Bethe-Bloch with the given mass hypothesis
109   Double_t bb = BetheBloch(track);
110   return TMath::Abs((track->GetRefESD()->GetTPCsignal() - bb) / bb);
111 }
112
113 //_________________________________________________________________________________________________
114 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget tgt, AliRsnDaughter *track)
115 {
116 //
117 // Cut checker.
118 //
119
120   // coherence check
121   if (tgt != AliRsnCut::kParticle) {
122     AliError(Form("Wrong target. Skipping cut", GetName()));
123     return kTRUE;
124   }
125
126   // if the required PID of the track is not the same as the
127   // reference of the cut, the cut is automatically skipped
128   if (track->RequiredPID() != fType) return kTRUE;
129
130   // retrieve the TPC signal
131   AliESDtrack *esd = track->GetRefESD();
132   if (!esd) {
133     AliError("ESD information unavailable");
134     return kTRUE;
135   }
136
137   // the cut range is the relative fraction of the value:
138   // BB*(1-fraction) < TPC < BB*(1+fraction)
139   // which means:
140   // -fraction < (TPC - BB)/BB < fraction
141   // so we must compute the cut value accordingly
142   fCutValueD = RelDiff(track);
143
144   // then, this cut is checked inside the range
145   return OkRange();
146 }
147
148 //_________________________________________________________________________________________________
149 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget /*tgt*/, AliRsnPairParticle* /*pair*/)
150 {
151 //
152 // Cut checker
153 //
154
155   AliWarning("Cannot apply this cut to pairs");
156   return kTRUE;
157 }
158
159 //_________________________________________________________________________________________________
160 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget /*tgt*/, AliRsnEvent* /*event*/)
161 {
162 //
163 // Cut checker
164 //
165
166   AliWarning("Cannot apply this cut to events");
167   return kTRUE;
168 }
169
170 //_________________________________________________________________________________________________
171 Bool_t AliRsnCutBetheBloch::IsSelected(ETarget /*tgt*/, AliRsnEvent* /*ev1*/, AliRsnEvent* /*ev2*/)
172 {
173 //
174 // Cut checker
175 //
176
177   AliWarning("Cannot apply this cut to event mixing");
178   return kTRUE;
179 }