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