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