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