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