AliRsnCutManager:
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnCutPID.cxx
CommitLineData
2dab9030 1//
2// Class AliRsnCutPID
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 "AliRsnCutPID.h"
24
25ClassImp(AliRsnCutPID)
26
27//_________________________________________________________________________________________________
28AliRsnCutPID::AliRsnCutPID() :
32992791 29 AliRsnCut(),
2dab9030 30 fPerfect(kFALSE),
31 fUseDefault(kTRUE)
32{
33//
34// Default constructor.
35// Sets the cut to realistic PID with default weights,
36// and defines the 'fMinI' value of the base class as the PID
37// to which we want to compare this object.
38//
39
40 Int_t i;
41
42 for (i = 0; i < kDetectors; i++)
43 {
44 fUseDetector[i] = kFALSE;
45 fPtThreshold[i] = 0.0;
46 fGoAboveThreshold[i] = kTRUE;
47 }
48
49 for (i = 0; i < AliPID::kSPECIES; i++)
50 {
51 fWeight[i] = 0.0;
52 fPrior[i] = 1.0;
53 }
32992791 54
55 SetTargetType(AliRsnTarget::kDaughter);
2dab9030 56}
57
58//_________________________________________________________________________________________________
59AliRsnCutPID::AliRsnCutPID(const char *name, AliPID::EParticleType pid, Double_t probMin, Bool_t perfectPID) :
60 AliRsnCut(name, AliRsnCut::kDaughter, (Int_t)pid),
61 fPerfect(perfectPID),
62 fUseDefault(kTRUE)
63{
64//
65// Default constructor.
66// Sets the cut to realistic PID with default weights,
67// and defines the 'fMinI' value of the base class as the PID
68// to which we want to compare this object.
69//
70
71 Int_t i;
72
73 for (i = 0; i < kDetectors; i++)
74 {
75 fUseDetector[i] = kFALSE;
76 fPtThreshold[i] = 0.0;
77 fGoAboveThreshold[i] = kTRUE;
78 }
79
80 for (i = 0; i < AliPID::kSPECIES; i++)
81 {
82 fWeight[i] = 0.0;
83 fPrior[i] = 1.0;
84 }
85
86 fMinD = probMin;
87 fMaxD = 1.000001;
32992791 88
89 SetTargetType(AliRsnTarget::kDaughter);
2dab9030 90}
91
92//_____________________________________________________________________________
93Bool_t AliRsnCutPID::CheckThreshold(EDetector det, Double_t value)
94{
95//
96// Checks if the passed value (track pT) stays in the
97// interval where the detector should be accepted
98//
99
100 if (!CheckBounds(det)) return kFALSE;
101
102 if (fGoAboveThreshold[det]) return (value >= fPtThreshold[det]);
103 else return (value <= fPtThreshold[det]);
104}
105
106//_________________________________________________________________________________________________
107Bool_t AliRsnCutPID::ComputeWeights(AliRsnDaughter *daughter)
108{
109//
110// Compute the PID weights using the class settings.
111// If the argument is an ESD track, customized weights can be computed
112// It the argument is a track (ESD or AOD), at least default weights
113// can be computed, otherwise, no weights can be defined.
114// The return value tells if the operation was successful.
115//
116
117 Int_t i, j;
118 Bool_t useDefault = fUseDefault;
119 Bool_t perfectPID = fPerfect;
35765294 120 if (perfectPID && !daughter->GetRefMC()) return kFALSE;
2dab9030 121 if (!daughter->GetRefESDtrack()) useDefault = kTRUE;
122 if (!daughter->GetRefESDtrack() && !daughter->GetRefAODtrack()) return kFALSE;
123
35765294 124 // if perfect PID ise required,
125 // compare the PDG code of the type stored in 'fMinI' of the cut
126 // and that of the particle which is checked, looking at its MC
2dab9030 127 if (perfectPID)
128 {
35765294 129 i = TMath::Abs(AliPID::ParticleCode(fMinI));
ec927a7d 130 j = daughter->GetPDG();
35765294 131 return (i == j);
2dab9030 132 }
133
134 // if default weights are (or need to be) used,
135 // they are taken directly and function exits
136 if (useDefault)
137 {
138 if (daughter->GetRefESDtrack())
139 daughter->GetRefESDtrack()->GetESDpid(fWeight);
140 else
141 {
142 for (i = 0; i < AliPID::kSPECIES; i++)
143 fWeight[i] = daughter->GetRefAODtrack()->PID()[i];
144 }
145 return kTRUE;
146 }
147
148 // if we arrive here, this means that we have an ESD track
149 // and we want to customize the PID
150 AliESDtrack *track = daughter->GetRefESDtrack();
151 Double_t w[kDetectors][AliPID::kSPECIES];
152 track->GetITSpid(w[kITS]);
153 track->GetTPCpid(w[kTPC]);
154 track->GetTRDpid(w[kTRD]);
155 track->GetTOFpid(w[kTOF]);
156 track->GetHMPIDpid(w[kHMPID]);
157
158 // if a detector is excluded or the track has not the right pT
159 // all related weights are set to 1 in order not to contribute
160 // to final multiplication
161 for (i = 0; i < kDetectors; i++)
162 {
163 if (!fUseDetector[i] || !CheckThreshold((EDetector)i, track->Pt()))
164 {
165 for (j = 0; j < AliPID::kSPECIES; j++) {
166 w[i][j] = 1.0;
167 }
168 }
169 }
170
171 // multiplicate all weights to compute final one
172 for (i = 0; i < AliPID::kSPECIES; i++)
173 {
174 fWeight[i] = w[kITS][i] * w[kTPC][i] * w[kTRD][i] * w[kTOF][i] * w[kHMPID][i];
175 }
176
177 return kTRUE;
178}
179
180//_________________________________________________________________________________________________
2e521c29 181Int_t AliRsnCutPID::RealisticPID(AliRsnDaughter * const daughter, Double_t &prob)
2dab9030 182{
183//
2e521c29 184// Combines the weights collected from chosen detectors with the priors
185// and gives the corresponding particle with the largest probability,
186// in terms of the AliPID particle type enumeration.
187// The argument, passed by reference, gives the corresponding probability,
188// since the cut could require that it is larger than a threshold.
2dab9030 189//
190
2e521c29 191 // try to compute the weights
192 if (!ComputeWeights(daughter))
2dab9030 193 {
2e521c29 194 prob = -1.0;
195 return AliPID::kUnknown;
2dab9030 196 }
197
2e521c29 198 // combine with priors and normalize
2dab9030 199 Int_t i;
200 Double_t sum = 0.0, w[AliPID::kSPECIES];
201 for (i = 0; i < AliPID::kSPECIES; i++)
202 {
203 w[i] = fWeight[i] * fPrior[i];
204 sum += w[i];
205 }
206 if (sum <= 0.0)
207 {
208 AliError("Sum = 0");
2e521c29 209 prob = -1.0;
210 return AliPID::kUnknown;
2dab9030 211 }
212 for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
213
214 // find the largest probability and related PID
2e521c29 215 Int_t ibest = 0;
216 prob = w[0];
2dab9030 217 for (i = 1; i < AliPID::kSPECIES; i++)
218 {
2e521c29 219 if (w[i] > prob)
2dab9030 220 {
2e521c29 221 prob = w[i];
222 ibest = i;
2dab9030 223 }
224 }
225
2e521c29 226 // return the value, while the probability
227 // will be consequentially set
228 return ibest;
229}
230
231//_________________________________________________________________________________________________
232Int_t AliRsnCutPID::PerfectPID(AliRsnDaughter * const daughter)
233{
234//
235// If MC is present, retrieve the particle corresponding to the used track
236// (using the fRefMC data member) and then return the true particle type
237// taken from the PDG code of the reference particle itself, converted
238// into the enumeration from AliPID object.
239//
240
241 // works only if the MC is present
242 if (!daughter->GetRefMC()) return AliPID::kUnknown;
243
244 // get the PDG code of the particle
ec927a7d 245 Int_t pdg = daughter->GetPDG();
2e521c29 246
247 // loop over all species listed in AliPID to find the match
248 Int_t i;
249 for (i = 0; i < AliPID::kSPECIES; i++)
250 {
251 if (AliPID::ParticleCode(i) == pdg) return i;
252 }
253
254 return AliPID::kUnknown;
255}
256
257//_________________________________________________________________________________________________
32992791 258Bool_t AliRsnCutPID::IsSelected(TObject *object)
2e521c29 259{
260//
261// Cut checker.
262//
2e521c29 263
264 // convert the object into the unique correct type
32992791 265
35e49ca5 266 if (!TargetOK(object))
96e9d35d 267 {
268 AliError(Form("[%s]: this cut works only with AliRsnDaughter objects", GetName()));
269 return kTRUE;
270 }
2dab9030 271
32992791 272 // if target is OK, do a dynamic cast
6aff5015 273 AliRsnDaughter *daughter = fDaughter;
32992791 274
2e521c29 275 // depending on the PID type, recalls the appropriate method:
276 // in case of perfect PID, checks only if the PID type is
277 // corresponding to the request in the cut (fMinI)
278 // while in case of realistic PID checks also the probability
279 // to be within the required interval
6aff5015 280 if (fPerfect && daughter)
2e521c29 281 {
282 fCutValueI = PerfectPID(daughter);
283 return OkValueI();
284 }
6aff5015 285 else if (daughter)
2e521c29 286 {
287 fCutValueI = RealisticPID(daughter, fCutValueD);
288 return OkValueI() && OkRangeD();
289 }
6aff5015 290 else
291 return kFALSE;
2dab9030 292}
293
2e521c29 294//__________________________________________________________________________________________________
2dab9030 295void AliRsnCutPID::IncludeDetector(EDetector det, Double_t threshold, Bool_t goAbove)
296{
297//
298// Include a detector for a customized weight computing
299// and specify also its eventual threshold and if the detector
300// must be used above or below the threshold.
301// By default the threshold is zero and detector is always used.
302//
303
304 if (!CheckBounds(det)) return;
305
306 fUseDetector[det] = kTRUE;
307 fPtThreshold[det] = threshold;
308 fGoAboveThreshold[det] = goAbove;
309}