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