Changes for #82873: Module debugging broken (Christian)
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnCutTrackQuality.cxx
1 //
2 // Class AliRsnCutTrackQuality
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
19 #include <Riostream.h>
20 #include <TFormula.h>
21 #include <TBits.h>
22
23 #include "AliLog.h"
24 #include "AliESDtrackCuts.h"
25
26 #include "AliRsnEvent.h"
27 #include "AliRsnDaughter.h"
28 #include "AliRsnCutTrackQuality.h"
29
30 ClassImp(AliRsnCutTrackQuality)
31
32 //_________________________________________________________________________________________________
33 AliRsnCutTrackQuality::AliRsnCutTrackQuality(const char *name) :
34    AliRsnCut(name, AliRsnTarget::kDaughter, 0.0, 0.0),
35    fFlagsOn(0x0),
36    fFlagsOff(0x0),
37    fRejectKinkDaughters(kTRUE),
38    fDCARfixed(kTRUE),
39    fDCARptFormula(""),
40    fDCARmax(1E20),
41    fDCAZfixed(kTRUE),
42    fDCAZptFormula(""),
43    fDCAZmax(1E20),
44    fSPDminNClusters(0),
45    fITSminNClusters(0),
46    fITSmaxChi2(1E20),
47    fTPCminNClusters(0),
48    fTPCmaxChi2(1E20),
49    fAODTestFilterBit(-1)
50 {
51 //
52 // Default constructor.
53 // Initializes all cuts in such a way that all of them are disabled.
54 //
55
56    SetPtRange(0.0, 1E20);
57    SetEtaRange(-1E20, 1E20);
58 }
59
60 //_________________________________________________________________________________________________
61 AliRsnCutTrackQuality::AliRsnCutTrackQuality(const AliRsnCutTrackQuality &copy) :
62    AliRsnCut(copy),
63    fFlagsOn(copy.fFlagsOn),
64    fFlagsOff(copy.fFlagsOff),
65    fRejectKinkDaughters(copy.fRejectKinkDaughters),
66    fDCARfixed(copy.fDCARfixed),
67    fDCARptFormula(copy.fDCARptFormula),
68    fDCARmax(copy.fDCARmax),
69    fDCAZfixed(copy.fDCAZfixed),
70    fDCAZptFormula(copy.fDCAZptFormula),
71    fDCAZmax(copy.fDCAZmax),
72    fSPDminNClusters(copy.fSPDminNClusters),
73    fITSminNClusters(copy.fITSminNClusters),
74    fITSmaxChi2(copy.fITSmaxChi2),
75    fTPCminNClusters(copy.fTPCminNClusters),
76    fTPCmaxChi2(copy.fTPCmaxChi2),
77    fAODTestFilterBit(copy.fAODTestFilterBit)
78 {
79 //
80 // Copy constructor.
81 // Just copy all data member values.
82 //
83
84    SetPtRange(copy.fPt[0], copy.fPt[1]);
85    SetEtaRange(copy.fEta[0], copy.fEta[1]);
86 }
87
88 //_________________________________________________________________________________________________
89 AliRsnCutTrackQuality& AliRsnCutTrackQuality::operator=(const AliRsnCutTrackQuality &copy)
90 {
91 //
92 // Assignment operator.
93 // Just copy all data member values.
94 //
95
96
97    fFlagsOn = copy.fFlagsOn;
98    fFlagsOff = copy.fFlagsOff;
99    fRejectKinkDaughters = copy.fRejectKinkDaughters;
100    fDCARfixed = copy.fDCARfixed;
101    fDCARptFormula = copy.fDCARptFormula;
102    fDCARmax = copy.fDCARmax;
103    fDCAZfixed = copy.fDCAZfixed;
104    fDCAZptFormula = copy.fDCAZptFormula;
105    fDCAZmax = copy.fDCAZmax;
106    fSPDminNClusters = copy.fSPDminNClusters;
107    fITSminNClusters = copy.fITSminNClusters;
108    fITSmaxChi2 = copy.fITSmaxChi2;
109    fTPCminNClusters = copy.fTPCminNClusters;
110    fTPCmaxChi2 = copy.fTPCmaxChi2;
111    fAODTestFilterBit = copy.fAODTestFilterBit;
112
113    SetPtRange(copy.fPt[0], copy.fPt[1]);
114    SetEtaRange(copy.fEta[0], copy.fEta[1]);
115
116    return (*this);
117 }
118
119 //_________________________________________________________________________________________________
120 void AliRsnCutTrackQuality::DisableAll()
121 {
122 //
123 // Disable all cuts
124 //
125
126    fFlagsOn = 0x0;
127    fFlagsOff = 0x0;
128    fRejectKinkDaughters = kFALSE;
129    fDCARfixed = kTRUE;
130    fDCARptFormula = "";
131    fDCARmax = 1E20;
132    fDCAZfixed = kTRUE;
133    fDCAZptFormula = "";
134    fDCAZmax = 1E20;
135    fSPDminNClusters = 0;
136    fITSminNClusters = 0;
137    fITSmaxChi2 = 1E20;
138    fTPCminNClusters = 0;
139    fTPCmaxChi2 = 1E20;
140    fAODTestFilterBit = -1;
141
142    SetPtRange(0.0, 1E20);
143    SetEtaRange(-1E20, 1E20);
144 }
145
146 //_________________________________________________________________________________________________
147 Bool_t AliRsnCutTrackQuality::IsSelected(TObject *object)
148 {
149 //
150 // Cut checker.
151 // Checks the type of object being evaluated
152 // and then calls the appropriate sub-function (for ESD or AOD)
153 //
154
155    // coherence check
156    if (!TargetOK(object)) return kFALSE;
157
158    // status is checked in the same way for all tracks, using AliVTrack
159    // as a convention, if a the collection of 'on' flags is '0x0', it
160    // is assumed that no flags are required, and this check is skipped;
161    // for the collection of 'off' flags this is not needed
162    AliVTrack *vtrack = fDaughter->Ref2Vtrack();
163    if (!vtrack) {
164       AliDebug(AliLog::kDebug + 2, "This object is not either an ESD nor AOD track");
165       return kFALSE;
166    }
167    ULong_t status   = (ULong_t)vtrack->GetStatus();
168    ULong_t checkOn  = status & fFlagsOn;
169    ULong_t checkOff = status & fFlagsOff;
170    if (fFlagsOn != 0x0 && checkOn != fFlagsOn) {
171       AliDebug(AliLog::kDebug + 2, Form("Failed flag check: required  %s", Binary(fFlagsOn)));
172       AliDebug(AliLog::kDebug + 2, Form("                   track has %s", Binary(status  )));
173       return kFALSE;
174    }
175    if (checkOff != 0) {
176       AliDebug(AliLog::kDebug + 2, Form("Failed flag check: forbidden %s", Binary(fFlagsOff)));
177       AliDebug(AliLog::kDebug + 2, Form("                   track has %s", Binary(status  )));
178       return kFALSE;
179    }
180    AliDebug(AliLog::kDebug + 3, Form("Flag check OK: required  %s", Binary(fFlagsOn)));
181    AliDebug(AliLog::kDebug + 3, Form("               forbidden %s", Binary(fFlagsOff)));
182    AliDebug(AliLog::kDebug + 3, Form("               track has %s", Binary(status  )));
183
184    // retrieve real object type
185    AliESDtrack *esdTrack = fDaughter->Ref2ESDtrack();
186    AliAODTrack *aodTrack = fDaughter->Ref2AODtrack();
187    if (esdTrack) {
188       AliDebug(AliLog::kDebug + 2, "Checking an ESD track");
189       return CheckESD(esdTrack);
190    } else if (aodTrack) {
191       AliDebug(AliLog::kDebug + 2, "Checking an AOD track");
192       return CheckAOD(aodTrack);
193    } else {
194       AliDebug(AliLog::kDebug + 2, Form("This object is not either an ESD nor AOD track, it is an %s", fDaughter->GetRef()->ClassName()));
195       return kFALSE;
196    }
197 }
198
199 //_________________________________________________________________________________________________
200 Bool_t AliRsnCutTrackQuality::CheckESD(AliESDtrack *track)
201 {
202 //
203 // Check an ESD track.
204 // This is done using the default track checker for ESD.
205 // It is declared static, not to recreate it every time.
206 //
207
208    static AliESDtrackCuts cuts;
209
210    // general acceptance/pt cuts
211    cuts.SetPtRange(fPt[0], fPt[1]);
212    cuts.SetEtaRange(fEta[0], fEta[1]);
213
214    // transverse DCA cuts
215    if (fDCARfixed)
216       cuts.SetMaxDCAToVertexXY(fDCARmax);
217    else
218       cuts.SetMaxDCAToVertexXYPtDep(fDCARptFormula.Data());
219
220    // longitudinal DCA cuts
221    if (fDCAZfixed)
222       cuts.SetMaxDCAToVertexZ(fDCAZmax);
223    else
224       cuts.SetMaxDCAToVertexZPtDep(fDCAZptFormula.Data());
225
226    // these options are always disabled in current version
227    cuts.SetDCAToVertex2D(kFALSE);
228    cuts.SetRequireSigmaToVertex(kFALSE);
229
230    // TPC related cuts for TPC+ITS tracks
231    cuts.SetMinNClustersTPC(fTPCminNClusters);
232    cuts.SetMaxChi2PerClusterTPC(fTPCmaxChi2);
233    cuts.SetAcceptKinkDaughters(!fRejectKinkDaughters);
234
235    // ITS related cuts for TPC+ITS tracks
236    if (fSPDminNClusters > 0)
237       cuts.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
238
239    // now that all is initialized, do the check
240    return cuts.IsSelected(track);
241 }
242
243 //_________________________________________________________________________________________________
244 Bool_t AliRsnCutTrackQuality::CheckAOD(AliAODTrack *track)
245 {
246 //
247 // Check an AOD track.
248 // This is done doing directly all checks, since there is not
249 // an equivalend checker for AOD tracks
250 //
251
252    // if a test bit is used, check it and skip the following
253    if (fAODTestFilterBit >= 0) {
254       UInt_t bit = (UInt_t)fAODTestFilterBit;
255       AliDebugClass(2, Form("Required a test filter bit for AOD check: %u (result: %s)", bit, (track->TestFilterBit(bit) ? "accept" : "reject")));
256       if (!track->TestFilterBit(bit)) 
257          return kFALSE;
258       else {
259          if (track->Pt() < fPt[0] || track->Pt() > fPt[1]) return kFALSE;
260          if (track->Eta() < fEta[0] || track->Eta() > fEta[1]) return kFALSE;
261          return kTRUE;
262       }
263    }
264
265    // try to retrieve the reference AOD event
266    AliAODEvent *aodEvent = 0x0;
267    if (fEvent) aodEvent = fEvent->GetRefAOD();
268    if (!aodEvent) {
269       AliError("AOD reference event is not initialized!");
270       return kFALSE;
271    }
272
273    // step #0: check SPD and ITS clusters
274    Int_t nSPD = 0;
275    nSPD  = TESTBIT(track->GetITSClusterMap(), 0);
276    nSPD += TESTBIT(track->GetITSClusterMap(), 1);
277    if (nSPD < fSPDminNClusters) {
278       AliDebug(AliLog::kDebug + 2, "Not enough SPD clusters in this track. Rejected");
279       return kFALSE;
280    }
281
282    // step #1: check number of clusters in TPC
283    if (track->GetTPCNcls() < fTPCminNClusters) {
284       AliDebug(AliLog::kDebug + 2, "Too few TPC clusters. Rejected");
285       return kFALSE;
286    }
287    if (track->GetITSNcls() < fITSminNClusters) {
288       AliDebug(AliLog::kDebug + 2, "Too few ITS clusters. Rejected");
289       return kFALSE;
290    }
291
292    // step #2: check chi square
293    if (track->Chi2perNDF() > fTPCmaxChi2) {
294       AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
295       return kFALSE;
296    }
297    if (track->Chi2perNDF() > fITSmaxChi2) {
298       AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
299       return kFALSE;
300    }
301
302    // step #3: reject kink daughters
303    AliAODVertex *vertex = track->GetProdVertex();
304    if (vertex && fRejectKinkDaughters) {
305       if (vertex->GetType() == AliAODVertex::kKink) {
306          AliDebug(AliLog::kDebug + 2, "Kink daughter. Rejected");
307          return kFALSE;
308       }
309    }
310
311    // step #4: DCA cut (transverse)
312    // --> reject all tracks not ITS refitted
313    Double_t b[2], cov[3];
314    vertex = aodEvent->GetPrimaryVertex();
315    if (!vertex) {
316       AliDebug(AliLog::kDebug + 2, "NULL vertex");
317       return kFALSE;
318    }
319    if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
320       AliDebug(AliLog::kDebug + 2, "Not ITS refitted");
321       return kFALSE;
322    }
323    if (!track->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b, cov)) {
324       AliDebug(AliLog::kDebug + 2, "Failed propagation to vertex");
325       return kFALSE;
326    }
327    // if the DCA cut is not fixed, compute current value
328    if (!fDCARfixed) {
329       static TString str(fDCARptFormula);
330       str.ReplaceAll("pt", "x");
331       static const TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data());
332       fDCARmax = dcaXY.Eval(track->Pt());
333    }
334    // check the cut
335    if (TMath::Abs(b[0]) > fDCARmax) {
336       AliDebug(AliLog::kDebug + 2, "Too large transverse DCA");
337       return kFALSE;
338    }
339
340    // step #5: DCA cut (longitudinal)
341    // the DCA has already been computed above
342    // if the DCA cut is not fixed, compute current value
343    if (!fDCAZfixed) {
344       static TString str(fDCAZptFormula);
345       str.ReplaceAll("pt", "x");
346       static const TFormula dcaZ(Form("%s_dcaXY", GetName()), str.Data());
347       fDCAZmax = dcaZ.Eval(track->Pt());
348    }
349    // check the cut
350    if (TMath::Abs(b[1]) > fDCAZmax) {
351       AliDebug(AliLog::kDebug + 2, "Too large longitudinal DCA");
352       return kFALSE;
353    }
354
355    // step #6: check eta/pt range
356    if (track->Eta() < fEta[0] || track->Eta() > fEta[1]) {
357       AliDebug(AliLog::kDebug + 2, "Outside ETA acceptance");
358       return kFALSE;
359    }
360    if (track->Pt() < fPt[0] || track->Pt() > fPt[1]) {
361       AliDebug(AliLog::kDebug + 2, "Outside PT acceptance");
362       return kFALSE;
363    }
364
365    // if we are here, all cuts were passed and no exit point was got
366    AliDebug(AliLog::kDebug + 2, "============================= ACCEPTED TRACK =====================================================");
367    return kTRUE;
368 }
369
370 //_________________________________________________________________________________________________
371 void AliRsnCutTrackQuality::Print(const Option_t *) const
372 {
373 //
374 // Print information on this cut
375 //
376
377    AliInfo(Form("Cut name                : %s", GetName()));
378    AliInfo(Form("Required flags (off, on): %lx %lx", fFlagsOn, fFlagsOff));
379    AliInfo(Form("Ranges in eta, pt       : %.2f - %.2f, %.2f - %.2f", fEta[0], fEta[1], fPt[0], fPt[1]));
380    AliInfo(Form("Kink daughters are      : %s", (fRejectKinkDaughters ? "rejected" : "accepted")));
381    AliInfo(Form("TPC requirements        : min. cluster = %d, max chi2 = %f", fTPCminNClusters, fTPCmaxChi2));
382    AliInfo(Form("ITS requirements        : min. cluster = %d (all), %d (SPD), max chi2 = %f", fITSminNClusters, fSPDminNClusters, fITSmaxChi2));
383
384    if (fDCARfixed) {
385       AliInfo(Form("DCA r cut               : fixed to %f cm", fDCARmax));
386    } else {
387       AliInfo(Form("DCA r cut formula       : %s", fDCARptFormula.Data()));
388    }
389
390    if (fDCAZfixed) {
391       AliInfo(Form("DCA z cut               : fixed to %f cm", fDCAZmax));
392    } else {
393       AliInfo(Form("DCA z cut formula       : %s", fDCAZptFormula.Data()));
394    }
395 }