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