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    fCheckOnlyFilterBit(kTRUE),
  52    fESDtrackCuts(0x0)
  53 {
  54 //
  55 // Default constructor.
  56 // Initializes all cuts in such a way that all of them are disabled.
  57 //
  58
  59    SetPtRange(0.0, 1E20);
  60    SetEtaRange(-1E20, 1E20);
  61 }
  62
  63 //_________________________________________________________________________________________________
  64 AliRsnCutTrackQuality::AliRsnCutTrackQuality(const AliRsnCutTrackQuality &copy) :
  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),
  79    fTPCmaxChi2(copy.fTPCmaxChi2),
  80    fCutMaxChi2TPCConstrainedVsGlobal(copy.fCutMaxChi2TPCConstrainedVsGlobal),
  81    fAODTestFilterBit(copy.fAODTestFilterBit),
  82    fCheckOnlyFilterBit(copy.fCheckOnlyFilterBit),
  83    fESDtrackCuts(copy.fESDtrackCuts)
  84 {
  85 //
  86 // Copy constructor.
  87 // Just copy all data member values.
  88 //
  89
  90    SetPtRange(copy.fPt[0], copy.fPt[1]);
  91    SetEtaRange(copy.fEta[0], copy.fEta[1]);
  92 }
  93
  94 //_________________________________________________________________________________________________
  95 AliRsnCutTrackQuality &AliRsnCutTrackQuality::operator=(const AliRsnCutTrackQuality &copy)
  96 {
  97 //
  98 // Assignment operator.
  99 // Just copy all data member values.
 100 //
 101
 102    if (this == &copy)
 103       return *this;
 104
 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;
 119    fAODTestFilterBit = copy.fAODTestFilterBit;
 120    fCheckOnlyFilterBit = copy.fCheckOnlyFilterBit;
 121    fESDtrackCuts = copy.fESDtrackCuts;
 122    SetPtRange(copy.fPt[0], copy.fPt[1]);
 123    SetEtaRange(copy.fEta[0], copy.fEta[1]);
 124
 125    return (*this);
 126 }
 127
 128 //_________________________________________________________________________________________________
 129 void AliRsnCutTrackQuality::DisableAll()
 130 {
 131 //
 132 // Disable all cuts
 133 //
 134
 135    fFlagsOn = 0x0;
 136    fFlagsOff = 0x0;
 137    fRejectKinkDaughters = kFALSE;
 138    fDCARfixed = kTRUE;
 139    fDCARptFormula = "";
 140    fDCARmax = 1E20;
 141    fDCAZfixed = kTRUE;
 142    fDCAZptFormula = "";
 143    fDCAZmax = 1E20;
 144    fSPDminNClusters = 0;
 145    fITSminNClusters = 0;
 146    fITSmaxChi2 = 1E20;
 147    fTPCminNClusters = 0;
 148    fTPCmaxChi2 = 1E20;
 149    fAODTestFilterBit = -1;
 150    if (fESDtrackCuts) {
 151       const char *cutsName = fESDtrackCuts->GetName();
 152       const char *cutsTitle = fESDtrackCuts->GetTitle();
 153       delete fESDtrackCuts;
 154       fESDtrackCuts = new AliESDtrackCuts(cutsName,cutsTitle);
 155    }
 156    SetPtRange(0.0, 1E20);
 157    SetEtaRange(-1E20, 1E20);
 158 }
 159
 160 //_________________________________________________________________________________________________
 161 Bool_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
 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
 176    AliVTrack *vtrack = fDaughter->Ref2Vtrack();
 177    if (!vtrack) {
 178       AliDebug(AliLog::kDebug + 2, "This object is not either an ESD nor AOD track");
 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) {
 185       AliDebug(AliLog::kDebug + 2, Form("Failed flag check: required  %s", Binary(fFlagsOn)));
 186       AliDebug(AliLog::kDebug + 2, Form("                   track has %s", Binary(status  )));
 187       return kFALSE;
 188    }
 189    if (checkOff != 0) {
 190       AliDebug(AliLog::kDebug + 2, Form("Failed flag check: forbidden %s", Binary(fFlagsOff)));
 191       AliDebug(AliLog::kDebug + 2, Form("                   track has %s", Binary(status  )));
 192       return kFALSE;
 193    }
 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  )));
 197
 198    // retrieve real object type
 199    AliESDtrack *esdTrack = fDaughter->Ref2ESDtrack();
 200    AliAODTrack *aodTrack = fDaughter->Ref2AODtrack();
 201    if (esdTrack) {
 202       AliDebug(AliLog::kDebug + 2, "Checking an ESD track");
 203       if (fESDtrackCuts)
 204          return fESDtrackCuts->IsSelected(esdTrack);
 205       else
 206          return CheckESD(esdTrack);
 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    }
 214 }
 215
 216 //_________________________________________________________________________________________________
 217 Bool_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   //static AliESDtrackCuts cuts;
 225   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);
 251    cuts.SetMaxChi2TPCConstrainedGlobal(fCutMaxChi2TPCConstrainedVsGlobal);
 252
 253    // ITS related cuts for TPC+ITS tracks
 254    if (fSPDminNClusters > 0)
 255       cuts.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 256    cuts.SetMaxChi2PerClusterITS(fITSmaxChi2);
 257
 258    // now that all is initialized, do the check
 259    return cuts.IsSelected(track);
 260 }
 261
 262 //_________________________________________________________________________________________________
 263 Bool_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 //
 270
 271    // if a test bit is used, check it and skip the following
 272    if (fAODTestFilterBit >= 0) {
 273       UInt_t bit = 1 << fAODTestFilterBit;
 274       AliDebugClass(2, Form("Required a test filter bit for AOD check: %u (result: %s)", bit, (track->TestFilterBit(bit) ? "accept" : "reject")));
 275       if (!track->TestFilterBit(bit))
 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;
 280          if (fCheckOnlyFilterBit) return kTRUE;
 281       }
 282    }
 283
 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
 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");
 298       return kFALSE;
 299    }
 300
 301
 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)
 332    // --> reject all tracks not ITS refitted
 333    Double_t b[2], cov[3];
 334    vertex = aodEvent->GetPrimaryVertex();
 335    if (!vertex) {
 336       AliDebug(AliLog::kDebug + 2, "NULL vertex");
 337       return kFALSE;
 338    }
 339    if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
 340       AliDebug(AliLog::kDebug + 2, "Not ITS refitted");
 341       return kFALSE;
 342    }
 343    if (!track->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b, cov)) {
 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) {
 349       TString str(fDCARptFormula);
 350       str.ReplaceAll("pt", "x");
 351       TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data());
 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) {
 364       TString str(fDCAZptFormula);
 365       str.ReplaceAll("pt", "x");
 366       TFormula dcaZ(Form("%s_dcaXY", GetName()), str.Data());
 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
 386    AliDebug(AliLog::kDebug + 2, "============================= ACCEPTED TRACK =====================================================");
 387    return kTRUE;
 388 }
 389
 390 //_________________________________________________________________________________________________
 391 void AliRsnCutTrackQuality::Print(const Option_t *) const
 392 {
 393 //
 394 // Print information on this cut
 395 //
 396
 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    }
 415
 416    AliInfo(Form("fAODTestFilterBit       : filter bit %i",fAODTestFilterBit));
 417    AliInfo(Form("fCheckOnlyFilterBit     : %i",((int) fCheckOnlyFilterBit)));
 418 }
 419 //__________________________________________________________________________________________________
 420 void 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);
 434    // SetITSmaxChi2(36);
 435    // SetMaxChi2TPCConstrainedGlobal(36);
 436    SetTPCminNClusters(70);
 437    SetTPCmaxChi2(4.0);
 438    SetRejectKinkDaughters();
 439    SetAODTestFilterBit(5);
 440 }
 441
 442 //__________________________________________________________________________________________________
 443 void 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 //__________________________________________________________________________________________________
 464 const 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 }