PWG2/RESONANCES/AliRsnAnalysisTaskEffPair.cxx

   1 //
   2 // Class AliRsnAnalysisTaskEffPair
   3 //
   4 // Inherits from basic AliRsnAnalysisTaskEff for efficiency,
   5 // and computed efficiencies for pairs
   6 //
   7 // author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
   8 //
   9
  10 #include "AliStack.h"
  11 #include "AliMCEvent.h"
  12 #include "AliESDEvent.h"
  13 #include "AliAODEvent.h"
  14
  15 #include "AliRsnPairDef.h"
  16 #include "AliRsnCutManager.h"
  17 #include "AliRsnAnalysisTaskEffPair.h"
  18
  19 ClassImp(AliRsnAnalysisTaskEffPair)
  20
  21 //_____________________________________________________________________________
  22 AliRsnAnalysisTaskEffPair::AliRsnAnalysisTaskEffPair(const char *name) :
  23    AliRsnAnalysisTaskEff(name),
  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 AliRsnAnalysisTaskEffPair::AliRsnAnalysisTaskEffPair(const AliRsnAnalysisTaskEffPair& copy) :
  34    AliRsnAnalysisTaskEff(copy),
  35    fMother()
  36 {
  37 //
  38 // Copy constrtuctor.
  39 //
  40 }
  41
  42 //_____________________________________________________________________________
  43 AliRsnAnalysisTaskEffPair& AliRsnAnalysisTaskEffPair::operator=(const AliRsnAnalysisTaskEffPair& copy)
  44 {
  45 //
  46 // Assignment operator.
  47 // Owned data members are meaningless for this operator.
  48 //
  49
  50    AliRsnAnalysisTaskEff::operator=(copy);
  51    return (*this);
  52 }
  53
  54 //_____________________________________________________________________________
  55 Int_t AliRsnAnalysisTaskEffPair::NGoodSteps()
  56 {
  57 //
  58 // Checks how many 'reconstruction' steps are passed by current daughter
  59 //
  60
  61    Int_t istep, count = 0;
  62    Int_t nSteps = fStepsRec.GetEntries();
  63
  64    for (istep = 0; istep < nSteps; istep++) {
  65       AliRsnCutManager *cutMgr = (AliRsnCutManager*)fStepsRec[istep];
  66
  67       if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
  68       if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
  69       if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
  70       if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
  71       if (!cutMgr->PassMotherCuts(&fMother)) break;
  72
  73       count++;
  74    }
  75
  76    return count;
  77 }
  78
  79 //_____________________________________________________________________________
  80 void AliRsnAnalysisTaskEffPair::ProcessEventESD()
  81 {
  82 //
  83 // Process current event with the definitions of the specified step in MC list
  84 // and store results in the container slot defined in second argument.
  85 // It is associated with the AliCFContainer with the name of the pair.
  86 //
  87
  88    AliESDEvent   *esd   = fRsnEvent[0].GetRefESD();
  89    AliMCEvent    *mc    = fRsnEvent[0].GetRefMCESD();
  90    AliStack      *stack = mc->Stack();
  91    TArrayI        indexes[2];
  92    Int_t          i, j, istep, imax, icheck, itrack[2], ipart;
  93    Int_t          pdg, label[2];
  94    Short_t        charge, pairDefMatch[2];
  95    TParticle     *part = 0x0;
  96    AliMCParticle *mother = 0x0;
  97
  98    // set pointers
  99    fMother.SetDaughter(0, &fDaughter[0]);
 100    fMother.SetDaughter(1, &fDaughter[1]);
 101    fMother.SetRefEvent(&fRsnEvent[0]);
 102
 103    // loop on definitions
 104    AliRsnPairDef *def = 0x0;
 105    TObjArrayIter nextDef(&fDefs);
 106    while ( (def = (AliRsnPairDef*)nextDef()) ) {
 107
 108       // loop on the MC list of particles
 109       for (ipart = 0; ipart < stack->GetNprimary(); ipart++) {
 110
 111          // reset daughters
 112          fDaughter[0].Reset();
 113          fDaughter[1].Reset();
 114
 115          // taks MC particle
 116          mother = (AliMCParticle*)mc->GetTrack(ipart);
 117
 118          // check that it is a binary decay and the PDG code matches
 119          if (mother->Particle()->GetNDaughters() != 2) continue;
 120          if (mother->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
 121
 122          // store the labels of the two daughters
 123          // and check that they are in the stack
 124          label[0] = mother->Particle()->GetFirstDaughter();
 125          label[1] = mother->Particle()->GetLastDaughter();
 126          if (label[0] < 0 || label[0] > stack->GetNtrack()) continue;
 127          if (label[1] < 0 || label[1] > stack->GetNtrack()) continue;
 128
 129          // for each daughter, check what slot in the pair definition it matches
 130          // or if it does not match any of them
 131          for (i = 0; i < 2; i++) {
 132             pairDefMatch[i] = -1;
 133             part = stack->Particle(label[i]);
 134             if (part) {
 135                pdg    = TMath::Abs(part->GetPdgCode());
 136                charge = (Short_t)(part->GetPDG()->Charge() / 3);
 137                if (def->GetDef1()->MatchesPDG(pdg) && def->GetDef1()->MatchesCharge(charge)) pairDefMatch[i] = 0;
 138                if (def->GetDef2()->MatchesPDG(pdg) && def->GetDef2()->MatchesCharge(charge)) pairDefMatch[i] = 1;
 139             }
 140          }
 141
 142          // if the two label match the two definitions for the pair
 143          // and if they are in the wrong order, swap them,
 144          // otherwise, if they don't match the definition in any order, skip
 145          if (pairDefMatch[0] == 1 && pairDefMatch[1] == 0) {
 146             icheck = label[0];
 147             label[0] = label[1];
 148             label[1] = icheck;
 149          }
 150          else if (pairDefMatch[0] < 0 || pairDefMatch[1] < 0) continue;
 151
 152          // from now on, we are sure that label[0] refers to the particle
 153          // that matches definitions of first daughter, and label[1] to
 154          // the particle that matches definitions of second daughter
 155          fDaughter[0].SetRefMC(mc->GetTrack(label[0]));
 156          fDaughter[1].SetRefMC(mc->GetTrack(label[1]));
 157
 158          // assign masses and fill the MC steps,
 159          // where reconstruction is not taken into account
 160          fMother.ComputeSum(def->GetMass1(), def->GetMass2());
 161          FillContainer(kTRUE, def);
 162
 163          // search for all reconstructed tracks which have these labels
 164          for (i = 0; i < 2; i++) indexes[i] = FindTracks(label[i], esd);
 165
 166          // if not both tracks have been reconstructed, stop here
 167          if (indexes[0].GetSize() < 1 || indexes[1].GetSize() < 1) continue;
 168
 169          // if both daughters were reconstructed
 170          // search for the best combination of indexes for this pair
 171          imax = itrack[0] = itrack[1] = 0;
 172          for (i = 0; i < indexes[0].GetSize(); i++) {
 173             for (j = 0; j < indexes[1].GetSize(); j++) {
 174                fDaughter[0].SetRef(esd->GetTrack(indexes[0][i]));
 175                fDaughter[1].SetRef(esd->GetTrack(indexes[1][j]));
 176                fMother.ComputeSum(def->GetMass1(), def->GetMass2());
 177                istep = NGoodSteps();
 178                if (istep > imax) {
 179                   itrack[0] = indexes[0][i];
 180                   itrack[1] = indexes[1][j];
 181                   imax = istep;
 182                }
 183             }
 184          }
 185
 186          // then assign definitely the best combination and fill rec container
 187          fDaughter[0].SetRef(esd->GetTrack(itrack[0]));
 188          fDaughter[1].SetRef(esd->GetTrack(itrack[1]));
 189          fMother.ComputeSum(def->GetMass1(), def->GetMass2());
 190          FillContainer(kFALSE, def);
 191       }
 192    }
 193 }
 194
 195 //_____________________________________________________________________________
 196 void AliRsnAnalysisTaskEffPair::ProcessEventAOD()
 197 {
 198 //
 199 // Process current event with the definitions of the specified step in MC list
 200 // and store results in the container slot defined in second argument.
 201 // It is associated with the AliCFContainer with the name of the pair.
 202 //
 203
 204    AliAODEvent  *aod     = fRsnEvent[0].GetRefAOD();
 205    TClonesArray *mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 206    if (!mcArray) return;
 207
 208    TArrayI           indexes[2];
 209    Int_t             i, j, pdg, imax, istep, icheck, itrack[2], label[2];
 210    AliAODMCParticle *part = 0x0, *mother = 0x0;
 211    Short_t           charge = 0, pairDefMatch[2];
 212
 213    // set pointers
 214    fMother.SetDaughter(0, &fDaughter[0]);
 215    fMother.SetDaughter(1, &fDaughter[1]);
 216
 217    // loop on definitions
 218    AliRsnPairDef *def = 0x0;
 219    TObjArrayIter nextDef(&fDefs);
 220    while ( (def = (AliRsnPairDef*)nextDef()) ) {
 221
 222       // loop on the MC list of particles
 223       TObjArrayIter next(mcArray);
 224       while ((mother = (AliAODMCParticle*)next())) {
 225
 226          // check that it is a binary decay and the PDG code matches
 227          if (mother->GetNDaughters() != 2) continue;
 228          if (mother->GetPdgCode() != def->GetMotherPDG()) continue;
 229
 230          // store the labels of the two daughters
 231          label[0] = mother->GetDaughter(0);
 232          label[1] = mother->GetDaughter(1);
 233
 234          // for each daughter, check what slot in the pair definition it matches
 235          // or if it does not match any of them
 236          for (i = 0; i < 2; i++) {
 237             pairDefMatch[i] = -1;
 238             part = (AliAODMCParticle*)mcArray->At(label[i]);
 239             if (part) {
 240                pdg    = TMath::Abs(part->GetPdgCode());
 241                charge = (Short_t)part->Charge();
 242                if (def->GetDef1()->MatchesPDG(pdg) && def->GetDef1()->MatchesCharge(charge)) pairDefMatch[i] = 0;
 243                if (def->GetDef2()->MatchesPDG(pdg) && def->GetDef2()->MatchesCharge(charge)) pairDefMatch[i] = 1;
 244             }
 245          }
 246
 247          // if the two label match the two definitions for the pair
 248          // and if they are in the wrong order, swap them,
 249          // otherwise, if they don't match the definition in any order, skip
 250          if (pairDefMatch[0] == 1 && pairDefMatch[1] == 0) {
 251             icheck = label[0];
 252             label[0] = label[1];
 253             label[1] = icheck;
 254          }
 255          else if (pairDefMatch[0] < 0 || pairDefMatch[1] < 0) continue;
 256
 257          // from now on, we are sure that label[0] refers to the particle
 258          // that matches definitions of first daughter, and label[1] to
 259          // the particle that matches definitions of second daughter
 260          fDaughter[0].SetRefMC((AliAODMCParticle*)mcArray->At(label[0]));
 261          fDaughter[1].SetRefMC((AliAODMCParticle*)mcArray->At(label[1]));
 262
 263          // assign masses and fill the MC steps,
 264          // where reconstruction is not taken into account
 265          fMother.ComputeSum(def->GetMass1(), def->GetMass2());
 266          FillContainer(kTRUE, def);
 267
 268          // search for all reconstructed tracks which have these labels
 269          for (i = 0; i < 2; i++) indexes[i] = FindTracks(label[i], aod);
 270
 271          // if not both tracks have been reconstructed, stop here
 272          if (indexes[0].GetSize() < 1 || indexes[1].GetSize() < 1) continue;
 273
 274          // if both daughters were reconstructed
 275          // search for the best combination of indexes for this pair
 276          imax = itrack[0] = itrack[1] = 0;
 277          for (i = 0; i < indexes[0].GetSize(); i++) {
 278             for (j = 0; j < indexes[1].GetSize(); j++) {
 279                fDaughter[0].SetRef(aod->GetTrack(indexes[0][i]));
 280                fDaughter[1].SetRef(aod->GetTrack(indexes[1][j]));
 281                fMother.ComputeSum(def->GetMass1(), def->GetMass2());
 282                istep = NGoodSteps();
 283                if (istep > imax) {
 284                   itrack[0] = indexes[0][i];
 285                   itrack[1] = indexes[1][j];
 286                }
 287             }
 288          }
 289
 290          // then assign definitely the best combination and fill rec container
 291          fDaughter[0].SetRef(aod->GetTrack(itrack[0]));
 292          fDaughter[1].SetRef(aod->GetTrack(itrack[1]));
 293          fMother.ComputeSum(def->GetMass1(), def->GetMass2());
 294          FillContainer(kFALSE, def);
 295       }
 296    }
 297 }
 298
 299 //_____________________________________________________________________________
 300 void AliRsnAnalysisTaskEffPair::FillContainer(Bool_t mcList, TObject *defObj)
 301 {
 302 //
 303 // Fill the container corresponding to current definition.
 304 //
 305
 306    // cast the object into the def
 307    if (!defObj->InheritsFrom(AliRsnPairDef::Class())) {
 308       AliError("Def object does not inherit from AliRsnPairDef!");
 309       return;
 310    }
 311    AliRsnPairDef *def = static_cast<AliRsnPairDef*>(defObj);
 312
 313    // retrieve container
 314    AliCFContainer *cont = (AliCFContainer*)fOutList->FindObject(def->GetName());
 315    if (!cont) return;
 316
 317    TObjArray &stepList =  (mcList ? fStepsMC : fStepsRec);
 318    Int_t      firstStep = (mcList ? 0 : ((Int_t)fStepsMC.GetEntries()));
 319    Int_t      iaxis, nAxes  = fAxes.GetEntries();
 320    Int_t      istep, nSteps = stepList.GetEntries();
 321    Bool_t     computeOK;
 322
 323    // compute values for all axes
 324    for (iaxis = 0; iaxis < nAxes; iaxis++) {
 325       AliRsnValue *fcnAxis = (AliRsnValue*)fAxes.At(iaxis);
 326       fVar[iaxis] = -1E10;
 327       switch (fcnAxis->GetTargetType()) {
 328          case AliRsnTarget::kMother:
 329             fcnAxis->SetSupportObject(def);
 330             computeOK = fcnAxis->Eval(&fMother, mcList);
 331             break;
 332          case AliRsnTarget::kEvent:
 333             computeOK = fcnAxis->Eval(&fRsnEvent[0]);
 334             break;
 335          default:
 336             AliError(Form("Allowed targets are mothers and events; cannot use axis '%s' which has target '%s'", fcnAxis->GetName(), fcnAxis->GetTargetTypeName()));
 337             computeOK = kFALSE;
 338       }
 339       if (computeOK) fVar[iaxis] = ((Float_t)fcnAxis->GetComputedValue());
 340    }
 341
 342    // fill all successful steps
 343    for (istep = 0; istep < nSteps; istep++) {
 344       AliRsnCutManager *cutMgr = (AliRsnCutManager*)stepList[istep];
 345
 346       if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
 347       if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
 348       if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
 349       if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
 350       if (!cutMgr->PassMotherCuts(&fMother)) break;
 351
 352       AliDebug(AliLog::kDebug + 2, Form("DEF: %s --> filling step %d", def->GetName(), istep));
 353       cont->Fill(fVar.GetArray(), istep + firstStep);
 354    }
 355 }