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