PWGLF/RESONANCES/AliRsnCutV0.cxx

   1 //
   2 // Class AliRsnCutV0
   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: Massimo Venaruzzo (massimo.venaruzzo@ts.infn.it)
  16 // modified: Enrico Fragiacomo (enrico.fragiacomo@ts.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 "AliRsnCutV0.h"
  29
  30 ClassImp(AliRsnCutV0)
  31
  32 //_________________________________________________________________________________________________
  33 //AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis) :
  34 AliRsnCutV0::AliRsnCutV0(const char *name, Int_t hypothesis, AliPID::EParticleType pid, AliPID::EParticleType pid2) :
  35    AliRsnCut(name, AliRsnTarget::kDaughter),
  36    fHypothesis(0),
  37    fMass(0.0),
  38    fTolerance(0.01),
  39    fMaxDCAVertex(0.3),
  40    fMinCosPointAngle(0.95),
  41    fMaxDaughtersDCA(0.5),
  42    fMaxRapidity(0.8),
  43    fPID(pid),
  44    fPID2(pid2),
  45    fPIDCut1(0),
  46    fPIDCut2(0),
  47    fPIDCut3(0),
  48    fESDtrackCuts(0x0),
  49    fCutQuality(Form("%sDaughtersQuality", name)),
  50    fAODTestFilterBit(5)
  51 {
  52 //
  53 // Default constructor.
  54 // Initializes all cuts in such a way that all of them are disabled.
  55 //
  56
  57    SetHypothesis(hypothesis);
  58 }
  59
  60 //_________________________________________________________________________________________________
  61 AliRsnCutV0::AliRsnCutV0(const AliRsnCutV0 &copy) :
  62    AliRsnCut(copy),
  63    fHypothesis(copy.fHypothesis),
  64    fMass(copy.fMass),
  65    fTolerance(copy.fTolerance),
  66    fMaxDCAVertex(copy.fMaxDCAVertex),
  67    fMinCosPointAngle(copy.fMinCosPointAngle),
  68    fMaxDaughtersDCA(copy.fMaxDaughtersDCA),
  69    fMaxRapidity(copy.fMaxRapidity),
  70    fPID(copy.fPID),
  71    fPID2(copy.fPID2),
  72    fPIDCut1(copy.fPIDCut1),
  73    fPIDCut2(copy.fPIDCut2),
  74    fPIDCut3(copy.fPIDCut3),
  75    fESDtrackCuts(copy.fESDtrackCuts),
  76    fCutQuality(copy.fCutQuality),
  77    fAODTestFilterBit(copy.fAODTestFilterBit)
  78 {
  79 //
  80 // Copy constructor.
  81 // Just copy all data member values.:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
  82
  83 //
  84    fCutQuality.SetPtRange(0.15, 1E+20);
  85    fCutQuality.SetEtaRange(-0.8, 0.8);
  86    fCutQuality.SetDCARmax(0.05);
  87    //fCutQuality.SetDCARPtFormula("0.0182+0.0350/pt^1.01");
  88    fCutQuality.SetDCAZmax(2.0);
  89    fCutQuality.SetSPDminNClusters(1);
  90    fCutQuality.SetITSminNClusters(0);
  91    fCutQuality.SetITSmaxChi2(1E+20);
  92    fCutQuality.SetTPCminNClusters(70);
  93    fCutQuality.SetTPCmaxChi2(4.0);
  94    fCutQuality.SetRejectKinkDaughters();
  95    fCutQuality.SetAODTestFilterBit(5);
  96
  97 }
  98
  99 //_________________________________________________________________________________________________
 100 AliRsnCutV0 &AliRsnCutV0::operator=(const AliRsnCutV0 &copy)
 101 {
 102 //
 103 // Assignment operator.
 104 // Just copy all data member values.
 105 //
 106    if (this == &copy)
 107      return *this;
 108    fHypothesis = copy.fHypothesis;
 109    fMass = copy.fMass;
 110    fTolerance = copy.fTolerance;
 111    fMaxDCAVertex = copy.fMaxDCAVertex;
 112    fMinCosPointAngle = copy.fMinCosPointAngle;
 113    fMaxDaughtersDCA = copy.fMaxDaughtersDCA;
 114    fMaxRapidity = copy.fMaxRapidity;
 115    fPID = copy.fPID;
 116    fPID2 = copy.fPID2;
 117    fPIDCut1 = copy.fPIDCut1;
 118    fPIDCut2 = copy.fPIDCut2;
 119    fPIDCut3 = copy.fPIDCut3;
 120    fESDtrackCuts = copy.fESDtrackCuts;
 121    fCutQuality = copy.fCutQuality;
 122    fAODTestFilterBit = copy.fAODTestFilterBit;
 123
 124    return (*this);
 125 }
 126
 127 //_________________________________________________________________________________________________
 128 Bool_t AliRsnCutV0::IsSelected(TObject *object)
 129 {
 130 //:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
 131
 132 // Cut checker.
 133 // Checks the type of object being evaluated
 134 // and then calls the appropriate sub-function (for ESD or AOD)
 135 //
 136
 137    // coherence check
 138    if (!TargetOK(object)) return kFALSE;
 139
 140    // check cast
 141    AliESDv0 *v0esd = fDaughter->Ref2ESDv0();
 142    AliAODv0 *v0aod = fDaughter->Ref2AODv0();
 143    //cout << fDaughter->GetRef()->ClassName() << ' ' << v0esd << ' ' << v0aod << endl;
 144
 145    // operate depending on cast:IsSelected: Object is not a V0 (RESONANCES/AliRsnCutV0.cxx:149)
 146
 147    if (v0esd) {
 148       return CheckESD(v0esd);
 149    } else if (v0aod) {
 150       return CheckAOD(v0aod);
 151    } else {
 152       AliDebugClass(1, "Object is not a V0");
 153       return kFALSE;
 154    }
 155 }
 156
 157 //_________________________________________________________________________________________________
 158 Bool_t AliRsnCutV0::CheckESD(AliESDv0 *v0)
 159 {
 160 //
 161 // Check an ESD V0.
 162 // This is done using the default track checker for ESD.
 163 // It is declared static, not to recreate it every time.
 164 //
 165
 166    AliDebugClass(1, "Check ESD");
 167    if (v0->GetOnFlyStatus()) {
 168       AliDebugClass(1, "Rejecting V0 in 'on fly' status");
 169       return kFALSE; // if kTRUE, then this V0 is recontructed
 170    }
 171
 172    // retrieve pointer to owner event
 173    AliESDEvent *lESDEvent = fEvent->GetRefESD();
 174    Double_t xPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetX();
 175    Double_t yPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetY();
 176    Double_t zPrimaryVertex = lESDEvent->GetPrimaryVertex()->GetZ();
 177    AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));
 178
 179    // retrieve the V0 daughters
 180    UInt_t lIdxPos      = (UInt_t) TMath::Abs(v0->GetPindex());
 181    UInt_t lIdxNeg      = (UInt_t) TMath::Abs(v0->GetNindex());
 182    AliESDtrack *pTrack = lESDEvent->GetTrack(lIdxPos);
 183    AliESDtrack *nTrack = lESDEvent->GetTrack(lIdxNeg);
 184
 185    // check quality cuts
 186    if (fESDtrackCuts) {
 187       AliDebugClass(2, "Checking quality cuts");
 188       if (!fESDtrackCuts->IsSelected(pTrack)) {
 189          AliDebugClass(2, "Positive daughter failed quality cuts");
 190          return kFALSE;
 191       }
 192       if (!fESDtrackCuts->IsSelected(nTrack)) {
 193          AliDebugClass(2, "Negative daughter failed quality cuts");
 194          return kFALSE;
 195       }
 196    }
 197
 198    // filter like-sign V0
 199    if ( TMath::Abs( ((pTrack->GetSign()) - (nTrack->GetSign())) ) < 0.1) {
 200       AliDebugClass(2, "Failed like-sign V0 check");
 201       return kFALSE;
 202    }
 203
 204
 205    // check compatibility with expected species hypothesis
 206    v0->ChangeMassHypothesis(fHypothesis);
 207    if ((TMath::Abs(v0->GetEffMass() - fMass)) > fTolerance) {
 208       AliDebugClass(2, "V0 is not in the expected inv mass range");
 209       return kFALSE;
 210    }
 211
 212    // topological checks
 213    if (TMath::Abs(v0->GetD(xPrimaryVertex, yPrimaryVertex, zPrimaryVertex)) > fMaxDCAVertex) {
 214       AliDebugClass(2, "Failed check on DCA to primary vertes");
 215       return kFALSE;
 216    }
 217    if (TMath::Abs(v0->GetV0CosineOfPointingAngle()) < fMinCosPointAngle) {
 218       AliDebugClass(2, "Failed check on cosine of pointing angle");
 219       return kFALSE;
 220    }
 221    if (TMath::Abs(v0->GetDcaV0Daughters()) > fMaxDaughtersDCA) {
 222       AliDebugClass(2, "Failed check on DCA between daughters");
 223       return kFALSE;
 224    }
 225    if (TMath::Abs(v0->Y(fHypothesis)) > fMaxRapidity) {
 226       AliDebugClass(2, "Failed check on V0 rapidity");
 227       return kFALSE;
 228    }
 229
 230
 231    // check PID on proton or antiproton from V0
 232
 233    // check initialization of PID object
 234    AliPIDResponse *pid = fEvent->GetPIDResponse();
 235    if (!pid) {
 236       AliFatal("NULL PID response");
 237       return kFALSE;
 238    }
 239
 240    // check if TOF is matched
 241    // and computes all values used in the PID cut
 242    //Bool_t   isTOFpos  = MatchTOF(ptrack);
 243    //Bool_t   isTOFneg  = MatchTOF(ntrack);
 244    Double_t pospTPC   = pTrack->GetTPCmomentum();
 245    Double_t negpTPC   = nTrack->GetTPCmomentum();
 246    //Double_t posp      = pTrack->P();
 247    //Double_t negp      = nTrack->P();
 248    Double_t posnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
 249    Double_t posnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
 250    //Double_t posnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
 251    Double_t negnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
 252    Double_t negnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
 253    //Double_t negnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
 254    Double_t maxTPC = 1E20;
 255    Double_t maxTPC2 = 1E20;
 256    //Double_t maxTOF = 1E20;
 257
 258    // applies the cut differently depending on the PID and the momentum
 259
 260    if(fHypothesis==kLambda0) {
 261       //if (isTOFpos) {
 262       // TPC: 5sigma cut for all
 263       //if (posnsTPC > 5.0) return kFALSE;
 264       // TOF: 3sigma
 265       // maxTOF = 3.0;
 266       //return (posnsTOF <= maxTOF);
 267       //} else {
 268       // TPC:
 269       // below 600 MeV: 4sigma
 270       // above 600 MeV: 3sigma
 271
 272       if (pospTPC <= 0.6 && fPID==AliPID::kProton)
 273          maxTPC = fPIDCut1;
 274       else if (pospTPC > 0.6 && fPID==AliPID::kProton)
 275          maxTPC = fPIDCut2;
 276       else
 277          return kFALSE;
 278
 279       maxTPC2 = fPIDCut3;
 280
 281       if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
 282          AliDebugClass(2, "Failed check on V0 PID");
 283          return kFALSE;
 284       }
 285    }
 286
 287    //}
 288
 289    if(fHypothesis==kLambda0Bar) {
 290       //if (isTOFneg) {
 291       // TPC: 5sigma cut for all
 292       //if (negnsTPC > 5.0) return kFALSE;
 293       // TOF: 3sigma
 294       // maxTOF = 3.0;
 295       //return (negnsTOF <= maxTOF);
 296       //} else {
 297       // TPC:
 298       // below 600 MeV: 4sigma
 299       // above 600 MeV: 3sigma
 300
 301       if (negpTPC <= 0.6 && fPID==AliPID::kProton)
 302          maxTPC = fPIDCut1;
 303       else if (negpTPC > 0.6 && fPID==AliPID::kProton)
 304          maxTPC = fPIDCut2;
 305       else
 306          return kFALSE;
 307
 308       maxTPC2 = fPIDCut3;
 309
 310       if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
 311          AliDebugClass(2, "Failed check on V0 PID");
 312          return kFALSE;
 313       }
 314    }
 315    //}
 316
 317
 318    // if we reach this point, all checks were successful
 319    AliDebugClass(2, "Good V0 (hallelujah)");
 320    return kTRUE;
 321 }
 322
 323 //_________________________________________________________________________________________________
 324 Bool_t AliRsnCutV0::CheckAOD(AliAODv0 *v0)
 325 {
 326 //
 327 // Check an AOD V0.
 328 // This is done doing directly all checks, since there is not
 329 // an equivalend checker for AOD tracks
 330 //
 331
 332    AliDebugClass(2, "Check AOD");
 333    if (v0->GetOnFlyStatus()) {
 334       AliDebugClass(2, "Rejecting V0 in 'on fly' status");
 335       return kFALSE; // if kTRUE, then this V0 is recontructed
 336    }
 337
 338
 339
 340    // retrieve pointer to owner event
 341    AliAODEvent *lAODEvent = fEvent->GetRefAOD();
 342    Double_t xPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetX();
 343    Double_t yPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetY();
 344    Double_t zPrimaryVertex = lAODEvent->GetPrimaryVertex()->GetZ();
 345    AliDebugClass(2, Form("Primary vertex: %f %f %f", xPrimaryVertex, yPrimaryVertex, zPrimaryVertex));
 346
 347    // retrieve the V0 daughters
 348    AliAODTrack *pTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
 349    AliAODTrack *nTrack = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
 350
 351
 352    // check quality cuts
 353    UInt_t  filtermapP = 9999;
 354    UInt_t  filtermapN = 9999;
 355    filtermapP = pTrack->GetFilterMap();
 356    filtermapN = nTrack->GetFilterMap();
 357
 358
 359    if ( !pTrack->TestFilterBit(fAODTestFilterBit)   ) {
 360       AliDebugClass(2, Form("Positive daughter failed quality cuts filtermapP=%d",filtermapP));
 361       return kFALSE;
 362    }
 363    if ( !nTrack->TestFilterBit(fAODTestFilterBit)   ) {
 364       AliDebugClass(2, Form("Negative daughter failed quality cuts filtermapN=%d",filtermapN));
 365       return kFALSE;
 366    }
 367
 368
 369
 370    /*AliDebugClass(1, Form("fESDtrackCuts=%p",fESDtrackCuts));
 371    if (fESDtrackCuts) { // use fESDtrackCuts to retrieve cuts values
 372
 373
 374       AliDebugClass(2, "Checking quality cuts");
 375
 376
 377
 378       const Bool_t  CutAcceptKinkDaughters  = fESDtrackCuts->GetAcceptKinkDaughters(); // 0 = kFalse
 379       const Float_t CutMaxChi2PerClusterTPC = fESDtrackCuts->GetMaxChi2PerClusterTPC();
 380       const Int_t   CutMinNClusterTPC       = fESDtrackCuts->GetMinNClusterTPC();
 381       const Bool_t  CutRequireTPCRefit      = fESDtrackCuts->GetRequireTPCRefit();
 382       //AliDebugClass(2, Form("accept kink=%d maxchi2=%f minnclSTPC=%d requireTPCrefit=%d", CutAcceptKinkDaughters,
 383       //        CutMaxChi2PerClusterTPC, CutMinNClusterTPC, CutRequireTPCRefit));
 384       //AliDebugClass(2, Form("pTrack->TPCNcls=%d", pTrack->GetTPCNcls() ));
 385       //AliDebugClass(2, Form("nTrack->TPCNcls=%d", nTrack->GetTPCNcls() ));
 386
 387
 388       if(pTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Positive daughter not MinNclsTPC"); return kFALSE;}
 389       if(nTrack->GetTPCNcls() < CutMinNClusterTPC) {AliDebugClass(2, "Negative daughter not MinNclsTPC"); return kFALSE;}
 390
 391    }*/
 392
 393    // filter like-sign V0
 394    if ( TMath::Abs( ((pTrack->Charge()) - (nTrack->Charge())) ) < 0.1) {
 395       AliDebugClass(2, "Failed like-sign V0 check");
 396       return kFALSE;
 397    }
 398
 399    // check compatibility with expected species hypothesis
 400    Double_t mass = 0.0;
 401    if(fHypothesis==kLambda0) {
 402       mass = v0->MassLambda();
 403    }
 404    else if (fHypothesis==kLambda0Bar) {
 405       mass = v0->MassAntiLambda();
 406    }
 407    if ((TMath::Abs(mass - fMass)) > fTolerance) {
 408       AliDebugClass(2, Form("V0 is not in the expected inv mass range  Mass: %d %f %f", fHypothesis, fMass, mass));
 409       return kFALSE;
 410    }
 411    AliDebugClass(2, Form("Mass: %d %f %f", fHypothesis, fMass, mass));
 412
 413
 414    // topological checks
 415    if (TMath::Abs(v0->DcaV0ToPrimVertex()) > fMaxDCAVertex) {
 416       AliDebugClass(2, Form("Failed check on DCA to primary vertes dca=%f maxdca=%f",TMath::Abs(v0->DcaV0ToPrimVertex()),fMaxDCAVertex));
 417       return kFALSE;
 418    }
 419    if (TMath::Abs( TMath::Cos(v0->OpenAngleV0()) ) < fMinCosPointAngle) {
 420       AliDebugClass(2, "Failed check on cosine of pointing angle");
 421       return kFALSE;
 422    }
 423    if (TMath::Abs(v0->DcaV0Daughters()) > fMaxDaughtersDCA) {
 424       AliDebugClass(2, "Failed check on DCA between daughters");
 425       return kFALSE;
 426    }
 427    if (TMath::Abs(v0->RapLambda()) > fMaxRapidity) {
 428       AliDebugClass(2, "Failed check on V0 rapidity");
 429       return kFALSE;
 430    }
 431
 432
 433
 434
 435
 436    // check initialization of PID object
 437    AliPIDResponse *pid = fEvent->GetPIDResponse();
 438    if (!pid) {
 439       AliFatal("NULL PID response");
 440       return kFALSE;
 441    }
 442
 443    // check if TOF is matched
 444    // and computes all values used in the PID cut
 445    //Bool_t   isTOFpos  = MatchTOF(ptrack);
 446    //Bool_t   isTOFneg  = MatchTOF(ntrack);
 447    Double_t pospTPC   = pTrack->GetTPCmomentum();
 448    Double_t negpTPC   = nTrack->GetTPCmomentum();
 449    //Double_t posp      = pTrack->P();
 450    //Double_t negp      = nTrack->P();
 451    Double_t posnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID));
 452    Double_t posnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(pTrack, fPID2));
 453    //Double_t posnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ptrack, fPID));
 454    Double_t negnsTPC   = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID));
 455    Double_t negnsTPC2  = TMath::Abs(pid->NumberOfSigmasTPC(nTrack, fPID2));
 456    //Double_t negnsTOF  = TMath::Abs(pid->NumberOfSigmasTOF(ntrack, fPID));
 457    Double_t maxTPC = 1E20;
 458    Double_t maxTPC2 = 1E20;
 459    //Double_t maxTOF = 1E20;
 460
 461    // applies the cut differently depending on the PID and the momentum
 462
 463    if(fHypothesis==kLambda0) {
 464
 465
 466       //if (isTOFpos) {
 467       // TPC: 5sigma cut for all
 468       //if (posnsTPC > 5.0) return kFALSE;
 469       // TOF: 3sigma
 470       // maxTOF = 3.0;
 471       //return (posnsTOF <= maxTOF);
 472       //} else {
 473       // TPC:
 474       // below 600 MeV: 4sigma
 475       // above 600 MeV: 3sigma
 476
 477       if (pospTPC <= 0.6 && fPID==AliPID::kProton)
 478          maxTPC = fPIDCut1; // EF safer value to run on MC
 479       else if (pospTPC > 0.6 && fPID==AliPID::kProton)
 480          maxTPC = fPIDCut2; // EF safer value to run on MC
 481       else
 482          return kFALSE;
 483
 484       maxTPC2 = fPIDCut3; // EF safer value to run on MC
 485
 486       if (! ((posnsTPC <= maxTPC) && (negnsTPC2 <= maxTPC2)) ) {
 487          AliDebugClass(2, "Failed check on V0 PID");
 488          return kFALSE;
 489       }
 490    }
 491
 492    if(fHypothesis==kLambda0Bar) {
 493
 494       //if (isTOFneg) {
 495       // TPC: 5sigma cut for all
 496       //if (negnsTPC > 5.0) return kFALSE;
 497       // TOF: 3sigma
 498       // maxTOF = 3.0;
 499       //return (negnsTOF <= maxTOF);
 500       //} else {
 501       // TPC:
 502       // below 600 MeV: 4sigma
 503       // above 600 MeV: 3sigma
 504
 505       if (negpTPC <= 0.6 && fPID==AliPID::kProton)
 506          maxTPC = fPIDCut1; // EF safer value to run on MC
 507       else if (negpTPC > 0.6 && fPID==AliPID::kProton)
 508          maxTPC = fPIDCut2; // EF safer value to run on MC
 509       else
 510          return kFALSE;
 511
 512       maxTPC2 = fPIDCut3; // EF safer value to run on MC
 513
 514       if(! ((negnsTPC <= maxTPC) && (posnsTPC2 <= maxTPC2)) ) {
 515          AliDebugClass(2, "Failed check on V0 PID");
 516          return kFALSE;
 517       }
 518    }
 519
 520
 521
 522    // if we reach this point, all checks were successful
 523    AliDebugClass(1, "Good AOD V0 (hallelujah)");
 524    AliDebugClass(1, Form("Mass: %d %f %f %d %d", fHypothesis, fMass, mass, filtermapP, filtermapN));
 525    return kTRUE;
 526
 527 }
 528
 529 //_________________________________________________________________________________________________
 530 void AliRsnCutV0::Print(const Option_t *) const
 531 {
 532 //
 533 // Print information on this cut
 534 //
 535 }