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 = dynamic_cast<AliVTrack*>(fDaughter->GetRef());
 161   if (!vtrack)
 162   {
 163     AliDebug(AliLog::kDebug + 2, Form("This object is not either an ESD nor AOD track, it is an %s", fDaughter->GetRef()->ClassName()));
 164     return kFALSE;
 165   }
 166   ULong_t status   = (ULong_t)vtrack->GetStatus();
 167   ULong_t checkOn  = status & fFlagsOn;
 168   ULong_t checkOff = status & fFlagsOff;
 169   if (fFlagsOn != 0x0 && checkOn != fFlagsOn)
 170   {
 171     AliDebug(AliLog::kDebug + 2, Form("Not all required flags are present: required %lx, track has %lx", fFlagsOn, status));
 172     return kFALSE;
 173   }
 174   if (checkOff != 0)
 175   {
 176     AliDebug(AliLog::kDebug + 2, Form("Some forbidden flags are present: required %lx, track has %lx", fFlagsOff, status));
 177     return kFALSE;
 178   }
 179
 180   // retrieve real object type
 181   AliESDtrack *esdTrack = fDaughter->GetRefESDtrack();
 182   AliAODTrack *aodTrack = fDaughter->GetRefAODtrack();
 183   if (esdTrack)
 184   {
 185     AliDebug(AliLog::kDebug + 2, "Checking an ESD track");
 186     return CheckESD(esdTrack);
 187   }
 188   else if (aodTrack)
 189   {
 190     AliDebug(AliLog::kDebug + 2, "Checking an AOD track");
 191     return CheckAOD(aodTrack);
 192   }
 193   else
 194   {
 195     AliDebug(AliLog::kDebug + 2, Form("This object is not either an ESD nor AOD track, it is an %s", fDaughter->GetRef()->ClassName()));
 196     return kFALSE;
 197   }
 198 }
 199
 200 //_________________________________________________________________________________________________
 201 Bool_t AliRsnCutTrackQuality::CheckESD(AliESDtrack *track)
 202 {
 203 //
 204 // Check an ESD track.
 205 // This is done using the default track checker for ESD.
 206 // It is declared static, not to recreate it every time.
 207 //
 208
 209   static AliESDtrackCuts cuts;
 210
 211   // general acceptance/pt cuts
 212   cuts.SetPtRange (fPt[0], fPt[1]);
 213   cuts.SetEtaRange(fEta[0], fEta[1]);
 214
 215   // transverse DCA cuts
 216   if (fDCARfixed)
 217     cuts.SetMaxDCAToVertexXY(fDCARmax);
 218   else
 219     cuts.SetMaxDCAToVertexXYPtDep(fDCARptFormula.Data());
 220
 221   // longitudinal DCA cuts
 222   if (fDCAZfixed)
 223     cuts.SetMaxDCAToVertexZ(fDCAZmax);
 224   else
 225     cuts.SetMaxDCAToVertexZPtDep(fDCAZptFormula.Data());
 226
 227   // these options are always disabled in current version
 228   cuts.SetDCAToVertex2D(kFALSE);
 229   cuts.SetRequireSigmaToVertex(kFALSE);
 230
 231   // TPC related cuts for TPC+ITS tracks
 232   cuts.SetMinNClustersTPC(fTPCminNClusters);
 233   cuts.SetMaxChi2PerClusterTPC(fTPCmaxChi2);
 234   cuts.SetAcceptKinkDaughters(!fRejectKinkDaughters);
 235
 236   // ITS related cuts for TPC+ITS tracks
 237   if (fSPDminNClusters > 0)
 238     cuts.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 239
 240   // now that all is initialized, do the check
 241   return cuts.IsSelected(track);
 242 }
 243
 244 //_________________________________________________________________________________________________
 245 Bool_t AliRsnCutTrackQuality::CheckAOD(AliAODTrack *track)
 246 {
 247 //
 248 // Check an AOD track.
 249 // This is done doing directly all checks, since there is not
 250 // an equivalend checker for AOD tracks
 251 //
 252
 253   // step #0: check SPD and ITS clusters
 254   Int_t nSPD = 0;
 255   nSPD  = TESTBIT(track->GetITSClusterMap(), 0);
 256   nSPD += TESTBIT(track->GetITSClusterMap(), 1);
 257   if (nSPD < fSPDminNClusters)
 258   {
 259     AliDebug(AliLog::kDebug + 2, "Not enough SPD clusters in this track. Rejected");
 260     return kFALSE;
 261   }
 262
 263   // step #1: check number of clusters in TPC
 264   if (track->GetTPCNcls() < fTPCminNClusters)
 265   {
 266     AliDebug(AliLog::kDebug + 2, "Too few TPC clusters. Rejected");
 267     return kFALSE;
 268   }
 269   if (track->GetITSNcls() < fITSminNClusters)
 270   {
 271     AliDebug(AliLog::kDebug + 2, "Too few ITS clusters. Rejected");
 272     return kFALSE;
 273   }
 274
 275   // step #2: check chi square
 276   if (track->Chi2perNDF() > fTPCmaxChi2)
 277   {
 278     AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
 279     return kFALSE;
 280   }
 281   if (track->Chi2perNDF() > fITSmaxChi2)
 282   {
 283     AliDebug(AliLog::kDebug + 2, "Bad chi2. Rejected");
 284     return kFALSE;
 285   }
 286
 287   // step #3: reject kink daughters
 288   AliAODVertex *vertex = track->GetProdVertex();
 289   if (vertex && fRejectKinkDaughters)
 290   {
 291     if (vertex->GetType() == AliAODVertex::kKink)
 292     {
 293       AliDebug(AliLog::kDebug + 2, "Kink daughter. Rejected");
 294       return kFALSE;
 295     }
 296   }
 297
 298   // step #4: DCA cut (transverse)
 299   Double_t b[2], cov[3];
 300   vertex = AliRsnTarget::GetCurrentEvent()->GetRefAOD()->GetPrimaryVertex();
 301   if (!vertex)
 302   {
 303     AliDebug(AliLog::kDebug + 2, "NULL vertex");
 304     return kFALSE;
 305   }
 306   if (!track->PropagateToDCA(vertex, AliRsnTarget::GetCurrentEvent()->GetRefAOD()->GetMagneticField(), kVeryBig, b, cov))
 307   {
 308     AliDebug(AliLog::kDebug + 2, "Failed propagation to vertex");
 309     return kFALSE;
 310   }
 311   // if the DCA cut is not fixed, compute current value
 312   if (!fDCARfixed)
 313   {
 314     static TString str(fDCARptFormula);
 315     str.ReplaceAll("pt", "x");
 316     static const TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data());
 317     fDCARmax = dcaXY.Eval(track->Pt());
 318   }
 319   // check the cut
 320   if (TMath::Abs(b[0]) > fDCARmax)
 321   {
 322     AliDebug(AliLog::kDebug + 2, "Too large transverse DCA");
 323     return kFALSE;
 324   }
 325
 326   // step #5: DCA cut (longitudinal)
 327   // the DCA has already been computed above
 328   // if the DCA cut is not fixed, compute current value
 329   if (!fDCAZfixed)
 330   {
 331     static TString str(fDCAZptFormula);
 332     str.ReplaceAll("pt", "x");
 333     static const TFormula dcaZ(Form("%s_dcaXY", GetName()), str.Data());
 334     fDCAZmax = dcaZ.Eval(track->Pt());
 335   }
 336   // check the cut
 337   if (TMath::Abs(b[1]) > fDCAZmax)
 338   {
 339     AliDebug(AliLog::kDebug + 2, "Too large longitudinal DCA");
 340     return kFALSE;
 341   }
 342
 343   // step #6: check eta/pt range
 344   if (track->Eta() < fEta[0] || track->Eta() > fEta[1])
 345   {
 346     AliDebug(AliLog::kDebug + 2, "Outside ETA acceptance");
 347     return kFALSE;
 348   }
 349   if (track->Pt() < fPt[0] || track->Pt() > fPt[1])
 350   {
 351     AliDebug(AliLog::kDebug + 2, "Outside PT acceptance");
 352     return kFALSE;
 353   }
 354
 355   // if we are here, all cuts were passed and no exit point was got
 356   return kTRUE;
 357 }
 358
 359 //_________________________________________________________________________________________________
 360 void AliRsnCutTrackQuality::Print(const Option_t *) const
 361 {
 362 //
 363 // Print information on this cut
 364 //
 365
 366   AliInfo(Form("Cut name                : %s", GetName()));
 367   AliInfo(Form("Required flags (off, on): %lx %lx", fFlagsOn, fFlagsOff));
 368   AliInfo(Form("Ranges in eta, pt       : %.2f - %.2f, %.2f - %.2f", fEta[0], fEta[1], fPt[0], fPt[1]));
 369   AliInfo(Form("Kink daughters are      : %s", (fRejectKinkDaughters ? "rejected" : "accepted")));
 370   AliInfo(Form("TPC requirements        : min. cluster = %d, max chi2 = %f", fTPCminNClusters, fTPCmaxChi2));
 371   AliInfo(Form("ITS requirements        : min. cluster = %d (all), %d (SPD), max chi2 = %f", fITSminNClusters, fSPDminNClusters, fITSmaxChi2));
 372
 373   if (fDCARfixed)
 374   {
 375     AliInfo(Form("DCA r cut               : fixed to %f cm", fDCARmax));
 376   }
 377   else
 378   {
 379     AliInfo(Form("DCA r cut formula       : %s", fDCARptFormula.Data()));
 380   }
 381
 382   if (fDCAZfixed)
 383   {
 384     AliInfo(Form("DCA z cut               : fixed to %f cm", fDCAZmax));
 385   }
 386   else
 387   {
 388     AliInfo(Form("DCA z cut formula       : %s", fDCAZptFormula.Data()));
 389   }
 390 }