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