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