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