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