PWG2/RESONANCES/AliRsnLoopEffPair.cxx

   1 //
   2 // Class AliRsnLoopEffPair
   3 //
   4 // Inherits from basic AliRsnLoopEff for efficiency,
   5 // and computed efficiencies for single-tracks
   6 //
   7 // author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
   8 //
   9
  10 #include "AliStack.h"
  11
  12 #include "AliRsnEvent.h"
  13 #include "AliRsnCutManager.h"
  14 #include "AliRsnDaughterDef.h"
  15 #include "AliRsnPairDef.h"
  16 #include "AliRsnLoopEffPair.h"
  17
  18 ClassImp(AliRsnLoopEffPair)
  19
  20 //_____________________________________________________________________________
  21 AliRsnLoopEffPair::AliRsnLoopEffPair(const char *name, AliRsnPairDef *def) :
  22    AliRsnLoopEff(name, 2),
  23    fDef(def),
  24    fMother()
  25 {
  26 //
  27 // Default constructor.
  28 // Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
  29 //
  30 }
  31
  32 //_____________________________________________________________________________
  33 AliRsnLoopEffPair::AliRsnLoopEffPair(const AliRsnLoopEffPair& copy) :
  34    AliRsnLoopEff(copy),
  35    fDef(copy.fDef),
  36    fMother(copy.fMother)
  37 {
  38 //
  39 // Copy constrtuctor.
  40 //
  41 }
  42
  43 //_____________________________________________________________________________
  44 AliRsnLoopEffPair& AliRsnLoopEffPair::operator=(const AliRsnLoopEffPair& copy)
  45 {
  46 //
  47 // Assignment operator.
  48 // Owned data members are meaningless for this operator.
  49 //
  50
  51    AliRsnLoopEff::operator=(copy);
  52    fDef = copy.fDef;
  53
  54    return (*this);
  55 }
  56
  57 //_____________________________________________________________________________
  58 Int_t AliRsnLoopEffPair::DoLoop(AliRsnEvent *rsn, AliRsnDaughterSelector*, AliRsnEvent*, AliRsnDaughterSelector*)
  59 {
  60 //
  61 // Loop on event and fill containers
  62 //
  63
  64    // check event cuts
  65    if (!OkEvent(rsn)) return 0;
  66
  67    // check presence of MC reference
  68    if (!rsn->GetRefMC()) {
  69       AliError("Need a MC to compute efficiency");
  70       return 0;
  71    }
  72
  73    // check event type:
  74    // must be ESD or AOD, and then use a bool to know in the rest
  75    if (!rsn->IsESD() && !rsn->IsAOD()) {
  76       AliError("Need to process ESD or AOD input");
  77       return 0;
  78    }
  79
  80    // additional coherence checks
  81    AliVEvent    *mcEvent = 0x0;
  82    AliStack     *listESD = 0x0;
  83    TClonesArray *listAOD = 0x0;
  84    Int_t         npart   = 0;
  85    if (rsn->IsESD()) {
  86       mcEvent = rsn->GetRefMCESD();
  87       listESD = ((AliMCEvent*)mcEvent)->Stack();
  88       if (!listESD) {
  89          AliError("Stack is not present");
  90          return 0;
  91       }
  92       npart = listESD->GetNprimary();
  93    } else {
  94       mcEvent = rsn->GetRefMCAOD();
  95       listAOD = (TClonesArray*)((AliAODEvent*)mcEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
  96       if (!listAOD) {
  97          AliError("Stack is not present");
  98          return 0;
  99       }
 100       npart = listAOD->GetEntries();
 101    }
 102
 103    // check number of particles
 104    if (!npart) {
 105       AliInfo("Empty event");
 106       return 0;
 107    }
 108
 109    // by default, assign daughters to mother
 110    // in the correct order (ref. pairDef)
 111    fMother.SetDaughter(0, &fDaughter[0]);
 112    fMother.SetDaughter(1, &fDaughter[1]);
 113    fMother.SetRefEvent(rsn);
 114    fDaughter[0].SetOwnerEvent(rsn);
 115    fDaughter[1].SetOwnerEvent(rsn);
 116
 117    // utility variables
 118    Bool_t            stop;
 119    Int_t             ipart, pdg, ndaughters, dindex[2], id, label, charge, imatch, index, istep, count = 0, nsteps = NStepsArray();
 120    AliVParticle     *mother = 0x0, *daughter = 0x0;
 121    AliMCParticle    *motherESD = 0x0, *daughterESD = 0x0;
 122    AliAODMCParticle *motherAOD = 0x0, *daughterAOD = 0x0;
 123
 124    // loop over particles
 125    for (ipart = 0; ipart < npart; ipart++) {
 126
 127       // get next particle and take some quantities
 128       // in different way from ESD and AOD MC
 129       if (rsn->IsESD()) {
 130          mother = mcEvent->GetTrack(ipart);
 131          motherESD = static_cast<AliMCParticle*>(mother);
 132          pdg = TMath::Abs(motherESD->Particle()->GetPdgCode());
 133          ndaughters = motherESD->Particle()->GetNDaughters();
 134       } else {
 135          mother = (AliVParticle*)listAOD->At(ipart);
 136          motherAOD = static_cast<AliAODMCParticle*>(mother);
 137          pdg = TMath::Abs(motherAOD->GetPdgCode());
 138          ndaughters = motherAOD->GetNDaughters();
 139       }
 140
 141       // skip particles with wrong PDG code,
 142       if (pdg != fDef->GetMotherPDG()) continue;
 143
 144       // fill first step
 145       fMother.SumMC().SetXYZM(mother->Px(), mother->Py(), mother->Pz(), fDef->GetMotherMass());
 146       GetOutput()->Fill(&fMother, 0);
 147       count++;
 148
 149       // reject particles with more/less than 2 daughters
 150       if (ndaughters != 2) continue;
 151
 152       // retrieve daughters
 153       // if one of them matches the definition for daughte #1 or #2 in pairDef,
 154       // assign it to the corresponding AliRsnDaughter member, and set it to 'good'
 155       // then, if both daughters are not 'good', skip the pair
 156       fDaughter[0].Reset();
 157       fDaughter[1].Reset();
 158       if (rsn->IsESD()) {
 159          dindex[0] = motherESD->Particle()->GetFirstDaughter();
 160          dindex[1] = motherESD->Particle()->GetLastDaughter();
 161       } else {
 162          dindex[0] = motherAOD->GetDaughter(0);
 163          dindex[1] = motherAOD->GetDaughter(1);
 164       }
 165
 166       // try to assign a daughter to each fDaughter[] data member
 167       stop = kFALSE;
 168       for (id = 0; id < 2; id++) {
 169          // avoid overflows
 170          if (dindex[id] < 0 || dindex[id] > npart) {
 171             AliWarning(Form("Found a stack overflow in dindex[%d]: value = %d, max = %d", id, dindex[id], npart));
 172             stop = kTRUE;
 173             break;
 174          }
 175          // retrieve daughter and copy some info
 176          if (rsn->IsESD()) {
 177             daughter = mcEvent->GetTrack(dindex[id]);
 178             daughterESD = static_cast<AliMCParticle*>(daughter);
 179             pdg = TMath::Abs(daughterESD->Particle()->GetPdgCode());
 180             charge = (Short_t)(daughterESD->Particle()->GetPDG()->Charge() / 3);
 181          } else {
 182             daughter = (AliAODMCParticle*)listAOD->At(dindex[id]);
 183             daughterAOD = static_cast<AliAODMCParticle*>(daughter);
 184             pdg = TMath::Abs(daughterAOD->GetPdgCode());
 185             charge = (Short_t)daughterAOD->Charge();
 186          }
 187          label = TMath::Abs(daughter->GetLabel());
 188          // check if can be assigned
 189          imatch = -1;
 190          if ( fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge) )
 191             imatch = 0;
 192          else if ( fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge) )
 193             imatch = 1;
 194          if (imatch < 0 || imatch > 1) continue;
 195          // assign
 196          label = daughter->GetLabel();
 197          fDaughter[imatch].SetRefMC(daughter);
 198          fDaughter[imatch].SetGood();
 199          index = FindTrack(label, mcEvent);
 200          if (index > 0) {
 201             fDaughter[imatch].SetRef(rsn->GetRef()->GetTrack(index));
 202          }
 203       }
 204
 205       // if both daughters were assigned, this means that the decay is correct
 206       // and then we fill second step
 207       if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
 208          GetOutput()->Fill(&fMother, 1);
 209          for (istep = 0; istep < nsteps; istep++) {
 210             AliRsnCutManager *cuts = (AliRsnCutManager*)fSteps[istep];
 211             if (!cuts->IsSelected(&fMother)) break;
 212             GetOutput()->Fill(&fMother, 2 + istep);
 213          }
 214       }
 215    }
 216
 217    return count;
 218 }