PWGLF/RESONANCES/AliRsnLoopPair.cxx

   1 //
   2 // *** Class AliRsnLoopPair ***
   3 //
   4 // "Core" method for defining the work on a pari of particles.
   5 // For one analysis, one must setup one of this for each pair he wants to analyze,
   6 // adding to it all analysis which he desires to do.
   7 // Here he defines the cuts, and the particle types and charges, and can add
   8 // functions which do different operations on the same pair, and some binning
   9 // with respect to some kinematic variables (eta, momentum)
  10 //
  11 // authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
  12 //          M. Vala (email: martin.vala@cern.ch)
  13 //
  14
  15 #include <Riostream.h>
  16 #include <TList.h>
  17 #include <TEntryList.h>
  18
  19 #include "AliLog.h"
  20
  21 #include "AliRsnEvent.h"
  22 #include "AliRsnPairDef.h"
  23 #include "AliRsnMother.h"
  24 #include "AliRsnCutSet.h"
  25 #include "AliRsnDaughterSelector.h"
  26
  27 #include "AliRsnLoopPair.h"
  28
  29 ClassImp(AliRsnLoopPair)
  30
  31 //_____________________________________________________________________________
  32 AliRsnLoopPair::AliRsnLoopPair(const char *name, AliRsnPairDef *def, Bool_t isMixed) :
  33    AliRsnLoop(name, isMixed),
  34    fTrueMC(kFALSE),
  35    fOnlyTrue(kFALSE),
  36    fUseMCRef(kFALSE),
  37    fCheckDecay(kFALSE),
  38    fRangeY(1E20),
  39    fRangeDCAproduct(1E20),
  40    fPairDef(def),
  41    fPairCuts(0x0),
  42    fMother()
  43 {
  44 //
  45 // Default constructor
  46 //
  47
  48    fListID[0] = -1;
  49    fListID[1] = -1;
  50 }
  51
  52 //_____________________________________________________________________________
  53 AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair &copy) :
  54    AliRsnLoop(copy),
  55    fTrueMC(copy.fTrueMC),
  56    fOnlyTrue(copy.fOnlyTrue),
  57    fUseMCRef(copy.fUseMCRef),
  58    fCheckDecay(copy.fCheckDecay),
  59    fRangeY(copy.fRangeY),
  60    fRangeDCAproduct(copy.fRangeDCAproduct),
  61    fPairDef(copy.fPairDef),
  62    fPairCuts(copy.fPairCuts),
  63    fMother(copy.fMother)
  64 {
  65 //
  66 // Copy constructor
  67 //
  68
  69    fListID[0] = copy.fListID[0];
  70    fListID[1] = copy.fListID[1];
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliRsnLoopPair &AliRsnLoopPair::operator=(const AliRsnLoopPair &copy)
  75 {
  76 //
  77 // Assignment operator
  78 //
  79
  80    AliRsnLoop::operator=(copy);
  81    if (this == &copy)
  82       return *this;
  83    fTrueMC = copy.fTrueMC;
  84    fOnlyTrue = copy.fOnlyTrue;
  85    fUseMCRef = copy.fUseMCRef;
  86    fCheckDecay = copy.fCheckDecay;
  87    fRangeY = copy.fRangeY;
  88    fRangeDCAproduct = copy.fRangeDCAproduct;
  89    fPairDef = copy.fPairDef;
  90    fPairCuts = copy.fPairCuts;
  91    fMother = copy.fMother;
  92    fListID[0] = copy.fListID[0];
  93    fListID[1] = copy.fListID[1];
  94
  95    return (*this);
  96 }
  97
  98 //_____________________________________________________________________________
  99 AliRsnLoopPair::~AliRsnLoopPair()
 100 {
 101 //
 102 // Destructor
 103 //
 104 }
 105
 106 //_____________________________________________________________________________
 107 void AliRsnLoopPair::Print(Option_t * /*option*/) const
 108 {
 109 //
 110 // Prints info about pair
 111 //
 112
 113    AliRsnLoop::Print();
 114 }
 115
 116 //_____________________________________________________________________________
 117 Bool_t AliRsnLoopPair::Init(const char *prefix, TList *list)
 118 {
 119 //
 120 // Initialization function.
 121 // Loops on all functions and eventual the ntuple, to initialize output objects.
 122 //
 123
 124    // assign data members relations
 125    fMother.SetDaughter(0, &fDaughter[0]);
 126    fMother.SetDaughter(1, &fDaughter[1]);
 127    AliInfo(Form("[%s] Initialization", GetName()));
 128
 129    TString name(prefix);
 130    name += '.';
 131    name += GetName();
 132 //    if (IsMixed()) name.Append("_mix");
 133
 134    return AliRsnLoop::Init(name.Data(), list);
 135 }
 136
 137 //_____________________________________________________________________________
 138 Int_t AliRsnLoopPair::DoLoop
 139 (AliRsnEvent *evMain, AliRsnDaughterSelector *selMain, AliRsnEvent *evMix, AliRsnDaughterSelector *selMix)
 140 {
 141 //
 142 // Loop function.
 143 // Computes what is needed from passed events.
 144 // Returns the number of pairs successfully processed.
 145 //
 146
 147    if (fIsMixed) {
 148       AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
 149       if (!evMix || !selMix) {
 150          AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
 151          return 0;
 152       }
 153    } else {
 154       AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
 155       evMix = evMain;
 156       selMix = selMain;
 157    }
 158    fMother.SetRefEvent(evMain);
 159
 160    // check cuts
 161    if (!OkEvent(evMain)) {
 162       AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
 163       return 0;
 164    }
 165    if (!OkEvent(evMix)) {
 166       AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
 167       return 0;
 168    }
 169
 170    // if it is required to loop over True MC, do this here and skip the rest of the method
 171    if (fTrueMC) return LoopTrueMC(evMain);
 172
 173    Int_t i0, i1, start, npairs = 0;
 174
 175    TEntryList *list0 = selMain->GetSelected(fListID[0], fPairDef->GetDef1().GetChargeC());
 176    TEntryList *list1 = selMix ->GetSelected(fListID[1], fPairDef->GetDef2().GetChargeC());
 177    if (!list0 || !list1) {
 178       AliError("Can't process NULL lists");
 179       return 0;
 180    }
 181    AliDebugClass(3, Form("[%s]: list counts: %lld, %lld", GetName(), list0->GetN(), list1->GetN()));
 182    if (!list0->GetN() || !list1->GetN()) {
 183       AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
 184       return 0;
 185    }
 186
 187    TObjArrayIter next(&fOutputs);
 188    AliRsnListOutput *out = 0x0;
 189    Long64_t iEntry1,iEntry2;
 190    for (i0 = 0; i0 < list0->GetN(); i0++) {
 191       iEntry1 = list0->GetEntry(i0);
 192       evMain->SetDaughter(fDaughter[0], iEntry1,fUseMCRef);
 193       fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
 194       start = 0;
 195       if (!fIsMixed && list0 == list1) start = i0 + 1;
 196       for (i1 = start; i1 < list1->GetN(); i1++) {
 197          iEntry2 = list1->GetEntry(i1);
 198          if (iEntry1 == iEntry2) continue;
 199          AliDebugClass(4, Form("Checking entries pair: %d (%lld) with %d (%lld)", i0, iEntry1, i1, iEntry2));
 200          evMix->SetDaughter(fDaughter[1], iEntry2,fUseMCRef);
 201          fDaughter[1].FillP(fPairDef->GetDef2().GetMass());
 202          fMother.Sum(0) = fDaughter[0].Prec() + fDaughter[1].Prec();
 203          fMother.Sum(1) = fDaughter[0].Psim() + fDaughter[1].Psim();
 204          fMother.Ref(0).SetXYZM(fMother.Sum(0).X(), fMother.Sum(0).Y(), fMother.Sum(0).Z(), fPairDef->GetMotherMass());
 205          fMother.Ref(1).SetXYZM(fMother.Sum(1).X(), fMother.Sum(1).Y(), fMother.Sum(1).Z(), fPairDef->GetMotherMass());
 206          // check rapidity range
 207          if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
 208             AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
 209             continue;
 210          }
 211          // check daughter's DCA product range
 212          if (TMath::Abs(fMother.DCAproduct()) > fRangeDCAproduct) {
 213             AliDebugClass(2, Form("[%s]: Outside daughter's DCA products range", GetName()));
 214             continue;
 215          }
 216          // check mother
 217          if (fOnlyTrue) {
 218             if (!IsTrueMother()) {
 219                AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
 220                continue;
 221             }
 222          }
 223          // check cuts
 224          if (fPairCuts) {
 225             if (!fPairCuts->IsSelected(&fMother)) {
 226                AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
 227                continue;
 228             }
 229          }
 230          // fill outputs
 231          next.Reset();
 232          while ( (out = (AliRsnListOutput *)next()) ) {
 233             if (out->Fill(&fMother)) npairs++;
 234             else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
 235          }
 236       }
 237    }
 238
 239    return npairs;
 240 }
 241
 242 //_____________________________________________________________________________
 243 Bool_t AliRsnLoopPair::IsTrueMother()
 244 {
 245 //
 246 // Checks to see if the mother comes from a true resonance.
 247 // It is triggered by the 'SetOnlyTrue()' function
 248 //
 249
 250    // check #1:
 251    // daughters have same mother with the right PDG code
 252    Int_t commonPDG = fMother.CommonMother();
 253    if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
 254    AliDebugClass(1, "Found a true mother");
 255
 256    // check #2:
 257    // checks if daughter have the right particle type
 258    // (activated by fCheckDecay)
 259    if (fCheckDecay) {
 260       AliRsnDaughterDef &def1 = fPairDef->GetDef1();
 261       AliRsnDaughterDef &def2 = fPairDef->GetDef2();
 262       if (!def1.MatchesPID(&fDaughter[0])) return kFALSE;
 263       if (!def2.MatchesPID(&fDaughter[1])) return kFALSE;
 264    }
 265    AliDebugClass(1, "Decay products match");
 266
 267    return kTRUE;
 268 }
 269
 270 //__________________________________________________________________________________________________
 271 Bool_t AliRsnLoopPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
 272 {
 273 //
 274 // Calls the appropriate assignment method
 275 //
 276
 277    // setup pointers
 278    fMother.SetDaughter(0, &fDaughter[0]);
 279    fMother.SetDaughter(1, &fDaughter[1]);
 280    fMother.SetRefEvent(rsnEvent);
 281    fDaughter[0].SetOwnerEvent(rsnEvent);
 282    fDaughter[1].SetOwnerEvent(rsnEvent);
 283
 284    if (rsnEvent->IsESD())
 285       return AssignMotherAndDaughtersESD(rsnEvent, ipart);
 286    else if (rsnEvent->IsAOD())
 287       return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
 288    else {
 289       AliError("Unrecognized input event");
 290       return kFALSE;
 291    }
 292 }
 293
 294 //__________________________________________________________________________________________________
 295 Bool_t AliRsnLoopPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
 296 {
 297 //
 298 // Gets a particle in the MC event and try to assign it to the mother.
 299 // If it has two daughters, retrieve them and assign also them.
 300 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
 301 //       for ESD and AOD, if available
 302 // ---
 303 // Implementation for ESD inputs
 304 //
 305
 306    AliMCEvent    *mc      = rsnEvent->GetRefMCESD();
 307    AliStack      *stack   = mc->Stack();
 308    AliMCParticle *mother  = (AliMCParticle *)mc->GetTrack(ipart);
 309    TParticle     *motherP = mother->Particle();
 310    Int_t          ntracks = stack->GetNtrack();
 311
 312    // check PDG code and exit if it is wrong
 313    if (TMath::Abs(motherP->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
 314
 315    // check number of daughters and exit if it is not 2
 316    if (motherP->GetNDaughters() < 2) return kFALSE;
 317    //if (!stack->IsPhysicalPrimary(ipart)) return kFALSE;
 318
 319    /*
 320    // check distance from primary vertex
 321    TLorentzVector vprod;
 322    motherP->ProductionVertex(vprod);
 323    if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
 324       AliDebugClass(1, "Distant production vertex");
 325       return kFALSE;
 326    }
 327    */
 328
 329    // get the daughters and check their PDG code and charge:
 330    // if they match one of the pair daughter definitions,
 331    // assign them as MC reference of the 'fDaughter' objects
 332    fDaughter[0].Reset();
 333    fDaughter[1].Reset();
 334    Int_t index[2] = {motherP->GetDaughter(0), motherP->GetDaughter(1)};
 335    Int_t i, pdg;
 336    Short_t charge;
 337    AliMCParticle *daughter = 0x0;
 338    for (i = 0; i < 2; i++) {
 339       // check index for stack
 340       if (index[i] < 0 || index[i] > ntracks) {
 341          AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
 342          return kFALSE;
 343       }
 344       // get daughter and its PDG and charge
 345       daughter = (AliMCParticle *)mc->GetTrack(index[i]);
 346       pdg      = TMath::Abs(daughter->Particle()->GetPdgCode());
 347       charge   = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
 348       // check if it matches one definition
 349       if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
 350          fDaughter[0].SetGood();
 351          fDaughter[0].SetRefMC(daughter);
 352          fDaughter[0].SetLabel(index[i]);
 353       } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
 354          fDaughter[1].SetGood();
 355          fDaughter[1].SetRefMC(daughter);
 356          fDaughter[1].SetLabel(index[i]);
 357       }
 358    }
 359
 360    // return success if both daughters were assigned
 361    if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
 362       return kTRUE;
 363    } else {
 364       fDaughter[0].Reset();
 365       fDaughter[1].Reset();
 366       return kFALSE;
 367    }
 368 }
 369
 370 //__________________________________________________________________________________________________
 371 Bool_t AliRsnLoopPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
 372 {
 373 //
 374 // Gets a particle in the MC event and try to assign it to the mother.
 375 // If it has two daughters, retrieve them and assign also them.
 376 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
 377 //       for ESD and AOD, if available
 378 // ---
 379 // Implementation for AOD inputs
 380 //
 381
 382    AliAODEvent      *aod     = rsnEvent->GetRefAOD();
 383    TClonesArray     *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 384    AliAODMCParticle *mother  = (AliAODMCParticle *)listAOD->At(ipart);
 385
 386    // check PDG code and exit if it is wrong
 387    if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
 388
 389    // check number of daughters and exit if it is not 2
 390    if (mother->GetNDaughters() < 2) return kFALSE;
 391    if (!mother->IsPrimary()) return kFALSE;
 392
 393    /*
 394    // check distance from primary vertex
 395    Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
 396    Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
 397    if (dv > fMaxDistPV) {
 398       AliDebugClass(1, "Distant production vertex");
 399       return kFALSE;
 400    }
 401    */
 402
 403    // get the daughters and check their PDG code and charge:
 404    // if they match one of the pair daughter definitions,
 405    // assign them as MC reference of the 'fDaughter' objects
 406    fDaughter[0].Reset();
 407    fDaughter[1].Reset();
 408    Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
 409    Int_t ntracks = listAOD->GetEntriesFast();
 410    Int_t i, pdg;
 411    Short_t charge;
 412    AliAODMCParticle *daughter = 0x0;
 413    for (i = 0; i < 2; i++) {
 414       // check index for stack
 415       if (index[i] < 0 || index[i] > ntracks) {
 416          AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
 417          return kFALSE;
 418       }
 419       // get daughter and its PDG and charge
 420       daughter = (AliAODMCParticle *)listAOD->At(index[i]);
 421       pdg      = TMath::Abs(daughter->GetPdgCode());
 422       charge   = (Short_t)(daughter->Charge() / 3);
 423       // check if it matches one definition
 424       if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
 425          fDaughter[0].SetGood();
 426          fDaughter[0].SetRefMC(daughter);
 427          fDaughter[0].SetLabel(index[i]);
 428       } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
 429          fDaughter[1].SetGood();
 430          fDaughter[1].SetRefMC(daughter);
 431          fDaughter[1].SetLabel(index[i]);
 432       }
 433    }
 434
 435    // return success if both daughters were assigned
 436    if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
 437       return kTRUE;
 438    } else {
 439       fDaughter[0].Reset();
 440       fDaughter[1].Reset();
 441       return kFALSE;
 442    }
 443 }
 444
 445 //_____________________________________________________________________________
 446 Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
 447 {
 448 //
 449 // Loop on event and fill containers
 450 //
 451
 452    // check presence of MC reference
 453    if (!rsn->GetRefMC()) {
 454       AliError("Need a MC to compute efficiency");
 455       return 0;
 456    }
 457
 458    // check event type:
 459    // must be ESD or AOD, and then use a bool to know in the rest
 460    if (!rsn->IsESD() && !rsn->IsAOD()) {
 461       AliError("Need to process ESD or AOD input");
 462       return 0;
 463    }
 464
 465    // retrieve the MC primary vertex position
 466    // and do some additional coherence checks
 467    //Double_t fVertex[3] = {0.0, 0.0, 0.0};
 468    Int_t npart = 0;
 469    if (rsn->IsESD()) {
 470       //TArrayF mcVertex(3);
 471       //AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
 472       //genEH->PrimaryVertex(mcVertex);
 473       //for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
 474       npart = rsn->GetRefMCESD()->GetNumberOfTracks();
 475    } else {
 476       //for (i = 0; i < 3; i++) fVertex[i] = 0.0;
 477       AliAODEvent *aod = rsn->GetRefMCAOD();
 478       TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 479       if (listAOD) npart = listAOD->GetEntries();
 480       //AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
 481       //if (mcH) mcH->GetVertex(fVertex);
 482    }
 483
 484    // check number of particles
 485    if (!npart) {
 486       AliInfo("Empty event");
 487       return 0;
 488    }
 489
 490    // utility variables
 491    Int_t ipart, count = 0;
 492    AliRsnDaughter check;
 493
 494    TObjArrayIter next(&fOutputs);
 495    AliRsnListOutput *out = 0x0;
 496
 497    // loop over particles
 498    for (ipart = 0; ipart < npart; ipart++) {
 499       // check i-th particle
 500       if (!AssignMotherAndDaughters(rsn, ipart)) continue;
 501       // if assignment was successful, for first step, use MC info
 502       fDaughter[0].SetRef(fDaughter[0].GetRefMC());
 503       fDaughter[1].SetRef(fDaughter[1].GetRefMC());
 504       fMother.SetRefEvent(rsn);
 505       fMother.ComputeSum(fPairDef->GetDef1().GetMass(), fPairDef->GetDef2().GetMass(), fPairDef->GetMotherMass());
 506       // check rapidity range
 507       if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
 508          AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
 509          continue;
 510       }
 511       if (TMath::Abs(fMother.DCAproduct()) > fRangeDCAproduct) {
 512          AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
 513          continue;
 514       }
 515       // fill outputs
 516       next.Reset();
 517       while ( (out = (AliRsnListOutput *)next()) ) {
 518          if (out->Fill(&fMother)) count++;
 519          else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
 520       }
 521    }
 522
 523    return count;
 524 }
 525