AliRsnCutManager:
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnCutPID.cxx
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
25 ClassImp(AliRsnCutPID)
26
27 //_________________________________________________________________________________________________
28 AliRsnCutPID::AliRsnCutPID() :
29   AliRsnCut(),
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   }
54   
55   SetTargetType(AliRsnTarget::kDaughter);
56 }
57
58 //_________________________________________________________________________________________________
59 AliRsnCutPID::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;
88   
89   SetTargetType(AliRsnTarget::kDaughter);
90 }
91
92 //_____________________________________________________________________________
93 Bool_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 //_________________________________________________________________________________________________
107 Bool_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;
120   if (perfectPID && !daughter->GetRefMC()) return kFALSE;
121   if (!daughter->GetRefESDtrack()) useDefault = kTRUE;
122   if (!daughter->GetRefESDtrack() && !daughter->GetRefAODtrack()) return kFALSE;
123   
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
127   if (perfectPID)
128   {
129     i = TMath::Abs(AliPID::ParticleCode(fMinI));
130     j = daughter->GetPDG();
131     return (i == j);
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 //_________________________________________________________________________________________________
181 Int_t AliRsnCutPID::RealisticPID(AliRsnDaughter * const daughter, Double_t &prob)
182 {
183 //
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.
189 //
190
191   // try to compute the weights
192   if (!ComputeWeights(daughter)) 
193   {
194     prob = -1.0;
195     return AliPID::kUnknown;
196   }
197   
198   // combine with priors and normalize
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");
209     prob = -1.0;
210     return AliPID::kUnknown;
211   }
212   for (i = 0; i < AliPID::kSPECIES; i++) w[i] /= sum;
213   
214   // find the largest probability and related PID
215   Int_t ibest = 0;
216   prob = w[0];
217   for (i = 1; i < AliPID::kSPECIES; i++)
218   {
219     if (w[i] > prob) 
220     {
221       prob = w[i];
222       ibest = i;
223     }
224   }
225   
226   // return the value, while the probability
227   // will be consequentially set
228   return ibest;
229 }
230   
231 //_________________________________________________________________________________________________
232 Int_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
245   Int_t pdg = daughter->GetPDG();
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 //_________________________________________________________________________________________________
258 Bool_t AliRsnCutPID::IsSelected(TObject *object)
259 {
260 //
261 // Cut checker.
262 //
263   
264   // convert the object into the unique correct type
265   
266   if (!TargetOK(object))
267   {
268     AliError(Form("[%s]: this cut works only with AliRsnDaughter objects", GetName()));
269     return kTRUE;
270   }
271   
272   // if target is OK, do a dynamic cast
273   AliRsnDaughter *daughter = fDaughter;
274   
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
280   if (fPerfect && daughter)
281   {
282     fCutValueI = PerfectPID(daughter);
283     return OkValueI();
284   }
285   else if (daughter)
286   {
287     fCutValueI = RealisticPID(daughter, fCutValueD);
288     return OkValueI() && OkRangeD();
289   }
290   else
291     return kFALSE;
292 }
293
294 //__________________________________________________________________________________________________
295 void 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 }