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