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 }