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 <Riostream.h>
  11
  12 #include "AliStack.h"
  13 #include "AliGenEventHeader.h"
  14 #include "AliAODMCHeader.h"
  15
  16 #include "AliRsnEvent.h"
  17 #include "AliRsnCutManager.h"
  18 #include "AliRsnDaughterDef.h"
  19 #include "AliRsnPairDef.h"
  20 #include "AliRsnLoopEffPair.h"
  21
  22 ClassImp(AliRsnLoopEffPair)
  23
  24 //_____________________________________________________________________________
  25 AliRsnLoopEffPair::AliRsnLoopEffPair(const char *name, AliRsnPairDef *def) :
  26    AliRsnLoopEff(name, 1),
  27    fDef(def),
  28    fMother()
  29 {
  30 //
  31 // Default constructor.
  32 // Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
  33 //
  34 }
  35
  36 //_____________________________________________________________________________
  37 AliRsnLoopEffPair::AliRsnLoopEffPair(const AliRsnLoopEffPair& copy) :
  38    AliRsnLoopEff(copy),
  39    fDef(copy.fDef),
  40    fMother(copy.fMother)
  41 {
  42 //
  43 // Copy constrtuctor.
  44 //
  45 }
  46
  47 //_____________________________________________________________________________
  48 AliRsnLoopEffPair& AliRsnLoopEffPair::operator=(const AliRsnLoopEffPair& copy)
  49 {
  50 //
  51 // Assignment operator.
  52 // Owned data members are meaningless for this operator.
  53 //
  54
  55    AliRsnLoopEff::operator=(copy);
  56    fDef = copy.fDef;
  57
  58    return (*this);
  59 }
  60
  61 //__________________________________________________________________________________________________
  62 Bool_t AliRsnLoopEffPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
  63 {
  64 //
  65 // Calls the appropriate assignment method
  66 //
  67
  68    // setup pointers
  69    fMother.SetDaughter(0, &fDaughter[0]);
  70    fMother.SetDaughter(1, &fDaughter[1]);
  71    fMother.SetRefEvent(rsnEvent);
  72    fDaughter[0].SetOwnerEvent(rsnEvent);
  73    fDaughter[1].SetOwnerEvent(rsnEvent);
  74
  75    if (rsnEvent->IsESD())
  76       return AssignMotherAndDaughtersESD(rsnEvent, ipart);
  77    else if (rsnEvent->IsAOD())
  78       return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
  79    else {
  80       AliError("Unrecognized input event");
  81       return kFALSE;
  82    }
  83 }
  84
  85 //__________________________________________________________________________________________________
  86 Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
  87 {
  88 //
  89 // Gets a particle in the MC event and try to assign it to the mother.
  90 // If it has two daughters, retrieve them and assign also them.
  91 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
  92 //       for ESD and AOD, if available
  93 // ---
  94 // Implementation for ESD inputs
  95 //
  96
  97    AliMCEvent    *mc      = rsnEvent->GetRefMCESD();
  98    AliStack      *stack   = mc->Stack();
  99    AliMCParticle *mother  = (AliMCParticle*)mc->GetTrack(ipart);
 100    TParticle     *motherP = mother->Particle();
 101    Int_t          ntracks = stack->GetNtrack();
 102
 103    // check PDG code and exit if it is wrong
 104    if (TMath::Abs(motherP->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
 105
 106    // check number of daughters and exit if it is not 2
 107    if (motherP->GetNDaughters() < 2) return kFALSE;
 108
 109    // check distance from primary vertex
 110    TLorentzVector vprod;
 111    motherP->ProductionVertex(vprod);
 112    if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
 113       AliDebugClass(1, "Distant production vertex");
 114       return kFALSE;
 115    }
 116
 117    // get the daughters and check their PDG code and charge:
 118    // if they match one of the pair daughter definitions,
 119    // assign them as MC reference of the 'fDaughter' objects
 120    fDaughter[0].Reset();
 121    fDaughter[1].Reset();
 122    Int_t index[2] = {motherP->GetFirstDaughter(), motherP->GetLastDaughter()};
 123    Int_t i, pdg;
 124    Short_t charge;
 125    AliMCParticle *daughter = 0x0;
 126    for (i = 0; i < 2; i++) {
 127       // check index for stack
 128       if (index[i] < 0 || index[i] > ntracks) {
 129          AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
 130          return kFALSE;
 131       }
 132       // get daughter and its PDG and charge
 133       daughter = (AliMCParticle*)mc->GetTrack(index[i]);
 134       pdg      = TMath::Abs(daughter->Particle()->GetPdgCode());
 135       charge   = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
 136       // check if it matches one definition
 137       if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
 138          fDaughter[0].SetGood();
 139          fDaughter[0].SetRefMC(daughter);
 140          fDaughter[0].SetLabel(index[i]);
 141       } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
 142          fDaughter[1].SetGood();
 143          fDaughter[1].SetRefMC(daughter);
 144          fDaughter[1].SetLabel(index[i]);
 145       }
 146    }
 147
 148    // return success if both daughters were assigned
 149    if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
 150       return kTRUE;
 151    } else {
 152       fDaughter[0].Reset();
 153       fDaughter[1].Reset();
 154       return kFALSE;
 155    }
 156 }
 157
 158 //__________________________________________________________________________________________________
 159 Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
 160 {
 161 //
 162 // Gets a particle in the MC event and try to assign it to the mother.
 163 // If it has two daughters, retrieve them and assign also them.
 164 // NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
 165 //       for ESD and AOD, if available
 166 // ---
 167 // Implementation for AOD inputs
 168 //
 169
 170    AliAODEvent      *aod     = rsnEvent->GetRefAOD();
 171    TClonesArray     *listAOD = (TClonesArray*)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 172    AliAODMCParticle *mother  = (AliAODMCParticle*)listAOD->At(ipart);
 173    Int_t             ntracks = listAOD->GetEntries();
 174
 175    // check PDG code and exit if it is wrong
 176    if (TMath::Abs(mother->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
 177
 178    // check number of daughters and exit if it is not 2
 179    if (mother->GetNDaughters() < 2) return kFALSE;
 180
 181    // check distance from primary vertex
 182    Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
 183    Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
 184    if (dv > fMaxDistPV) {
 185       AliDebugClass(1, "Distant production vertex");
 186       return kFALSE;
 187    }
 188
 189    // get the daughters and check their PDG code and charge:
 190    // if they match one of the pair daughter definitions,
 191    // assign them as MC reference of the 'fDaughter' objects
 192    fDaughter[0].Reset();
 193    fDaughter[1].Reset();
 194    Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
 195    Int_t i, pdg;
 196    Short_t charge;
 197    AliAODMCParticle *daughter = 0x0;
 198    for (i = 0; i < 2; i++) {
 199       // check index for stack
 200       if (index[i] < 0 || index[i] > ntracks) {
 201          AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
 202          return kFALSE;
 203       }
 204       // get daughter and its PDG and charge
 205       daughter = (AliAODMCParticle*)listAOD->At(index[i]);
 206       pdg      = TMath::Abs(daughter->GetPdgCode());
 207       charge   = (Short_t)(daughter->Charge() / 3);
 208       // check if it matches one definition
 209       if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
 210          fDaughter[0].SetGood();
 211          fDaughter[0].SetRefMC(daughter);
 212          fDaughter[0].SetLabel(index[i]);
 213       } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
 214          fDaughter[1].SetGood();
 215          fDaughter[1].SetRefMC(daughter);
 216          fDaughter[1].SetLabel(index[i]);
 217       }
 218    }
 219
 220    // return success if both daughters were assigned
 221    if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
 222       return kTRUE;
 223    } else {
 224       fDaughter[0].Reset();
 225       fDaughter[1].Reset();
 226       return kFALSE;
 227    }
 228 }
 229
 230 //_____________________________________________________________________________
 231 Int_t AliRsnLoopEffPair::DoLoop(AliRsnEvent *rsn, AliRsnDaughterSelector*, AliRsnEvent*, AliRsnDaughterSelector*)
 232 {
 233 //
 234 // Loop on event and fill containers
 235 //
 236
 237    // check event cuts
 238    if (!OkEvent(rsn)) return 0;
 239
 240    // retrieve output
 241    fOutput = (AliRsnListOutput*)fOutputs[0];
 242
 243    // check presence of MC reference
 244    if (!rsn->GetRefMC()) {
 245       AliError("Need a MC to compute efficiency");
 246       return 0;
 247    }
 248
 249    // check presence of event
 250    if (!rsn->GetRef()) {
 251       AliError("Need an event to compute efficiency");
 252       return 0;
 253    }
 254
 255    // check event type:
 256    // must be ESD or AOD, and then use a bool to know in the rest
 257    if (!rsn->IsESD() && !rsn->IsAOD()) {
 258       AliError("Need to process ESD or AOD input");
 259       return 0;
 260    }
 261
 262    // retrieve the MC primary vertex position
 263    // and do some additional coherence checks
 264    Int_t i, npart = 0;
 265    if (rsn->IsESD()) {
 266       TArrayF mcVertex(3);
 267       AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
 268       genEH->PrimaryVertex(mcVertex);
 269       for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
 270       npart = rsn->GetRefMCESD()->GetNumberOfTracks();
 271    } else {
 272       for (i = 0; i < 3; i++) fVertex[i] = 0.0;
 273       AliAODEvent *aod = rsn->GetRefMCAOD();
 274       TClonesArray *listAOD = (TClonesArray*)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 275       if (listAOD) npart = listAOD->GetEntries();
 276       AliAODMCHeader *mcH = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
 277       if (mcH) mcH->GetVertex(fVertex);
 278    }
 279
 280    // check number of particles
 281    if (!npart) {
 282       AliInfo("Empty event");
 283       return 0;
 284    }
 285
 286    // utility variables
 287    Int_t ipart, istep, count = 0, nsteps = fSteps.GetEntries();
 288    Int_t ntracks = rsn->GetAbsoluteSum();
 289    AliRsnDaughter check;
 290
 291    // loop over particles
 292    for (ipart = 0; ipart < npart; ipart++) {
 293       // check i-th particle
 294       if (!AssignMotherAndDaughters(rsn, ipart)) continue;
 295       // if assignment was successful, for first step, use MC info
 296       fDaughter[0].SetRef(fDaughter[0].GetRefMC());
 297       fDaughter[1].SetRef(fDaughter[1].GetRefMC());
 298       fMother.SetRefEvent(rsn);
 299       fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
 300       fOutput->Fill(&fMother, 0);
 301       count++;
 302       // for each further step, try to find two tracks which pass the related cuts
 303       for (istep = 0; istep < nsteps; istep++) {
 304          AliRsnCutManager *cuts = (AliRsnCutManager*)fSteps[istep];
 305          fDaughter[0].SetBad();
 306          fDaughter[1].SetBad();
 307          for (i = 0; i < ntracks; i++) {
 308             rsn->SetDaughter(check, i);
 309             if (!cuts->PassCommonDaughterCuts(&check)) continue;
 310             if (TMath::Abs(check.GetLabel()) == fDaughter[0].GetLabel()) {
 311                if (!cuts->PassDaughter1Cuts(&check)) continue;
 312                fDaughter[0].SetRef(check.GetRef());
 313                fDaughter[0].SetGood();
 314             } else if (TMath::Abs(check.GetLabel()) == fDaughter[1].GetLabel()) {
 315                if (!cuts->PassDaughter2Cuts(&check)) continue;
 316                fDaughter[1].SetRef(check.GetRef());
 317                fDaughter[1].SetGood();
 318             }
 319             if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
 320                fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
 321                if (cuts->PassMotherCuts(&fMother)) {
 322                   fOutput->Fill(&fMother, 1 + istep);
 323                   break;
 324                }
 325             }
 326          }
 327       }
 328    }
 329
 330    return count;
 331 }