PWG2/RESONANCES/AliRsnAnalysisEffSE.cxx

   1 //
   2 // Class AliRsnAnalysisEffSE
   3 //
   4 // TODO
   5 // TODO
   6 //
   7 // authors: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
   8 //          Martin Vala (martin.vala@cern.ch)
   9 //
  10 #include <Riostream.h>
  11 #include "AliStack.h"
  12 #include "AliESDEvent.h"
  13 #include "AliAODEvent.h"
  14 #include "AliMCEvent.h"
  15
  16 #include "AliCFContainer.h"
  17 #include "AliRsnCutManager.h"
  18 #include "AliRsnValue.h"
  19 #include "AliRsnAnalysisEffSE.h"
  20 #include "AliRsnPairDef.h"
  21 #include "AliRsnCutSet.h"
  22
  23 ClassImp(AliRsnAnalysisEffSE)
  24
  25 //_____________________________________________________________________________
  26 AliRsnAnalysisEffSE::AliRsnAnalysisEffSE(const char *name) :
  27   AliRsnVAnalysisTaskSE(name),
  28   fUseITSSA(kTRUE),
  29   fUseGlobal(kTRUE),
  30   fStepListMC(0),
  31   fStepListESD(0),
  32   fAxisList("AliRsnValue", 0),
  33   fPairDefList(0),
  34   fContainerList(0x0),
  35   fOutList(0x0),
  36   fVar(0),
  37   fPair(),
  38   fEventCuts("eventCuts", AliRsnCut::kEvent)
  39 {
  40 //
  41 // Default constructor.
  42 //
  43
  44   AliDebug(AliLog::kDebug+2,"<-");
  45
  46   DefineOutput(2, TList::Class());
  47
  48   AliDebug(AliLog::kDebug+2,"->");
  49 }
  50
  51 //_____________________________________________________________________________
  52 AliRsnAnalysisEffSE::AliRsnAnalysisEffSE(const AliRsnAnalysisEffSE& copy) :
  53   AliRsnVAnalysisTaskSE(copy),
  54   fUseITSSA(copy.fUseITSSA),
  55   fUseGlobal(copy.fUseGlobal),
  56   fStepListMC(copy.fStepListMC),
  57   fStepListESD(copy.fStepListESD),
  58   fAxisList(copy.fAxisList),
  59   fPairDefList(copy.fPairDefList),
  60   fContainerList(copy.fContainerList),
  61   fOutList(0x0),
  62   fVar(0),
  63   fPair(),
  64   fEventCuts(copy.fEventCuts)
  65 {
  66 //
  67 // Copy constrtuctor
  68 //
  69 }
  70
  71 //_____________________________________________________________________________
  72 void AliRsnAnalysisEffSE::AddStepMC(AliRsnCutManager *mgr)
  73 {
  74 //
  75 // Add a step on montecarlo
  76 //
  77
  78   fStepListMC.AddLast(mgr);
  79 }
  80
  81 //_____________________________________________________________________________
  82 void AliRsnAnalysisEffSE::AddStepESD(AliRsnCutManager *mgr)
  83 {
  84 //
  85 // Add a step on ESD
  86 //
  87
  88   fStepListESD.AddLast(mgr);
  89 }
  90
  91 //_____________________________________________________________________________
  92 void AliRsnAnalysisEffSE::AddAxis(AliRsnValue *axis)
  93 {
  94 //
  95 // Add a new axis
  96 //
  97
  98   Int_t n = fAxisList.GetEntries();
  99   new (fAxisList[n]) AliRsnValue(*axis);
 100 }
 101
 102 //_____________________________________________________________________________
 103 void AliRsnAnalysisEffSE::RsnUserCreateOutputObjects()
 104 {
 105 //
 106 // Creation of output objects.
 107 // These are created through the utility methods in the analysis manager,
 108 // which produces a list of histograms for each specified set of pairs.
 109 // Each of these lists is added to the main list of this task.
 110 //
 111
 112   AliDebug(AliLog::kDebug+2,"<-");
 113
 114   // get number of steps and axes
 115   Int_t iaxis  = 0;
 116   Int_t nAxes  = fAxisList.GetEntries();
 117   Int_t nSteps = (Int_t)fStepListMC.GetEntries() + (Int_t)fStepListESD.GetEntries();
 118
 119   if (!nSteps) {
 120     AliError("No steps defined");
 121     return;
 122   }
 123   if (!nAxes) {
 124     AliError("No axes defined");
 125     return;
 126   }
 127
 128   // initialize variable list
 129   fVar.Set(nAxes);
 130
 131   // retrieve number of bins for each axis
 132   Int_t *nBins = new Int_t[nAxes];
 133   for (iaxis = 0; iaxis < nAxes; iaxis++)
 134   {
 135     AliRsnValue *fcnAxis = (AliRsnValue*)fAxisList.At(iaxis);
 136     nBins[iaxis] = fcnAxis->GetArray().GetSize() - 1;
 137   }
 138
 139   // create ouput list of containers
 140   fContainerList = new TList;
 141   fContainerList->SetOwner();
 142   fContainerList->SetName(Form("%s_containers", GetName()));
 143
 144   // initialize output list
 145   OpenFile(2);
 146   fOutList = new TList();
 147   fOutList->SetOwner();
 148
 149   // create the containers
 150   Int_t i = 0, nDef = (Int_t)fPairDefList.GetEntries();
 151   for (i = 0; i < nDef; i++)
 152   {
 153     AliRsnPairDef *def = (AliRsnPairDef*)fPairDefList[i];
 154     AliCFContainer *cont = new AliCFContainer(Form("%s", def->GetPairName()), "", nSteps, nAxes, nBins);
 155     // set the bin limits for each axis
 156     for (iaxis = 0; iaxis < nAxes; iaxis++)
 157     {
 158       AliRsnValue *fcnAxis = (AliRsnValue*)fAxisList.At(iaxis);
 159       cont->SetBinLimits(iaxis, fcnAxis->GetArray().GetArray());
 160     }
 161     // add the container to output list
 162     fContainerList->Add(cont);
 163   }
 164
 165   fOutList->Add(fContainerList);
 166   fOutList->Print();
 167
 168   PostData(2, fOutList);
 169
 170   // clear heap
 171   delete [] nBins;
 172
 173   AliDebug(AliLog::kDebug+2,"->");
 174 }
 175
 176 //_____________________________________________________________________________
 177 void AliRsnAnalysisEffSE::RsnUserExec(Option_t*)
 178 {
 179 //
 180 // Execution of the analysis task.
 181 // Recovers the input event and processes it with all included pair objects.
 182 // In this case, we NEED to have ESD and MC, otherwise we cannod do anything.
 183 //ULTIMO UNDO
 184 /*
 185   AliRsnDaughter trk;
 186   for (Int_t i = 0; i <= AliPID::kSPECIES; i++)
 187   {
 188     cout << AliPID::ParticleName((AliPID::EParticleType)i) << ": " << endl;
 189     for (Int_t m = 0; m < AliRsnDaughter::kMethods; m++)
 190     {
 191       cout << "-- method: " << AliRsnDaughter::MethodName((AliRsnDaughter::EPIDMethod)m) << endl;
 192       Char_t   sign[2] = {'+', '-'};
 193       for (Int_t s = 0; s < 2; s++)
 194       {
 195         TArrayI *a = fRsnPIDIndex.GetTracksArray((AliRsnDaughter::EPIDMethod)m, sign[s], (AliPID::EParticleType)i);
 196         Int_t n = a->GetSize();
 197         for (Int_t j = 0; j < n; j++)
 198         {
 199           Int_t k = a->At(j);
 200           cout << "-- -- track " << Form("%4d ", k) << ": ";
 201           fRsnEvent.SetDaughter(trk, k);
 202           cout << "charge = " << (trk.IsPos() ? "+ " : (trk.IsNeg() ? "- " : "0 "));
 203           cout << "truePID = " << Form("%10s ", AliPID::ParticleName(trk.PerfectPID()));
 204           cout << "realPID = " << Form("%10s ", AliPID::ParticleName(trk.RealisticPID()));
 205           cout << endl;
 206           cout << "-- -- -- weights (computed): ";
 207           for (Int_t q = 0; q < AliPID::kSPECIES; q++)
 208             cout << Form("%15.12f", trk.ComputedWeights()[q]) << ' ';
 209           cout << endl;
 210           cout << "-- -- -- weights (original): ";
 211           for (Int_t q = 0; q < AliPID::kSPECIES; q++)
 212             cout << Form("%15.12f", trk.GetRef()->PID()[q]) << ' ';
 213           cout << endl;
 214         }
 215       }
 216     }
 217   } return;
 218   */
 219
 220   AliDebug(AliLog::kDebug+2,"<-");
 221   if (fMCOnly && fMCEvent)
 222   {
 223     fRsnEvent.SetRef  (fMCEvent);
 224     fRsnEvent.SetRefMC(fMCEvent);
 225   }
 226   else if (fESDEvent)
 227   {
 228     fRsnEvent.SetRef  (fESDEvent);
 229     fRsnEvent.SetRefMC(fMCEvent);
 230   }
 231   else if (fAODEventOut)
 232   {
 233     fRsnEvent.SetRef  (fAODEventOut);
 234     fRsnEvent.SetRefMC(fAODEventOut);
 235   }
 236   else if (fAODEventIn)
 237   {
 238     fRsnEvent.SetRef  (fAODEventIn);
 239     fRsnEvent.SetRefMC(fAODEventIn);
 240   }
 241   else {
 242     AliError("NO ESD or AOD object!!! Skipping ...");
 243     return;
 244   }
 245
 246   // if general event cuts are added to the task (recommended)
 247   // they are checked here on the RSN event interface and,
 248   // if the event does not pass them, it is skipped and ProcessInfo
 249   // is updated accordingly
 250   if (!fEventCuts.IsSelected(&fRsnEvent))
 251   {
 252     fTaskInfo.SetEventUsed(kFALSE);
 253     return;
 254   }
 255
 256   // if cuts are passed or not cuts were defined,
 257   // update the task info before processing the event
 258   fTaskInfo.SetEventUsed(kTRUE);
 259
 260   // process first MC steps and then ESD steps
 261   AliRsnPairDef *pairDef = 0;
 262   TObjArrayIter iter(&fPairDefList);
 263   while ( (pairDef = (AliRsnPairDef*)iter.Next()) )
 264   {
 265     if (fRsnEvent.IsESD()) ProcessEventESD(pairDef);
 266     if (fRsnEvent.IsAOD()) ProcessEventAOD(pairDef);
 267   }
 268
 269   // Post the data
 270   PostData(2, fOutList);
 271
 272   AliDebug(AliLog::kDebug+2,"->");
 273 }
 274
 275 //_____________________________________________________________________________
 276 void AliRsnAnalysisEffSE::ProcessEventESD(AliRsnPairDef *pairDef)
 277 {
 278 //
 279 // Process current event with the definitions of the specified step in MC list
 280 // and store results in the container slot defined in second argument.
 281 // It is associated with the AliCFContainer with the name of the pair.
 282 //
 283
 284   AliStack      *stack = fRsnEvent.GetRefMCESD()->Stack();
 285   AliESDEvent   *esd   = fRsnEvent.GetRefESD();
 286   AliMCEvent    *mc    = fRsnEvent.GetRefMCESD();
 287   AliMCParticle *mother;
 288
 289   if (!pairDef) return;
 290   AliCFContainer *cont = (AliCFContainer*)fContainerList->FindObject(pairDef->GetPairName());
 291   if (!cont) return;
 292
 293   // get informations from pairDef
 294   Int_t   pdgM = 0;
 295   Int_t   pdgD[2] = {0, 0};
 296   Short_t chargeD[2] = {0, 0};
 297   pdgM       = pairDef->GetMotherPDG();
 298   pdgD   [0] = AliPID::ParticleCode(pairDef->GetPID(0));
 299   pdgD   [1] = AliPID::ParticleCode(pairDef->GetPID(1));
 300   chargeD[0] = pairDef->GetChargeShort(0);
 301   chargeD[1] = pairDef->GetChargeShort(1);
 302
 303   // other utility variables
 304   Int_t      first, j, ipart;
 305   Int_t      label[2] = {-1, -1};
 306   Short_t    charge[2] = {0, 0};
 307   Short_t    pairDefMatch[2] = {-1, -1};
 308   Int_t      esdIndex[2] = {-1, -1};
 309   TParticle *part[2] = {0, 0};
 310
 311   // in this case, we first take the resonance from MC
 312   // and then we find its daughters and compute cuts on them
 313   for (ipart = 0; ipart < stack->GetNprimary(); ipart++)
 314   {
 315     // take a track from the MC event
 316     mother = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 317
 318     // check that it is a binary decay and the PDG code matches
 319     if (mother->Particle()->GetNDaughters() != 2) continue;
 320     if (mother->Particle()->GetPdgCode() != pdgM) continue;
 321
 322     // store the labels of the two daughters
 323     label[0] = mother->Particle()->GetFirstDaughter();
 324     label[1] = mother->Particle()->GetLastDaughter();
 325
 326     // retrieve the particles and other stuff
 327     // check if they match the order in the pairDef
 328     for (j = 0; j < 2; j++)
 329     {
 330       if (label[j] < 0) continue;
 331       part[j]   = stack->Particle(label[j]);
 332       pdgD[j]   = TMath::Abs(part[j]->GetPdgCode());
 333       charge[j] = (Short_t)(part[j]->GetPDG()->Charge() / 3);
 334       if (pdgD[j] == AliPID::ParticleCode(pairDef->GetPID(0)) && charge[j] == pairDef->GetChargeShort(0))
 335         pairDefMatch[j] = 0;
 336       else if (pdgD[j] == AliPID::ParticleCode(pairDef->GetPID(1)) && charge[j] == pairDef->GetChargeShort(1))
 337         pairDefMatch[j] = 1;
 338       else
 339         pairDefMatch[j] = -1;
 340
 341       // find corresponding ESD particle: first try rejecting fakes,
 342       // and in case of failure, try accepting fakes
 343       esdIndex[j] = FindESDtrack(label[j], esd, kTRUE);
 344       //TArrayI idx = FindESDtracks(label[j], esd);
 345       //for (Int_t kk = 0; kk < idx.GetSize(); kk++) cout << "DAUGHTER " << j << " --> FOUND INDEX: " << idx[kk] << endl;
 346       if (esdIndex[j] < 0) esdIndex[j] = FindESDtrack(label[j], esd, kFALSE);
 347       //cout << "DAUGHTER " << j << " SINGLE FOUND INDEX = " << esdIndex[j] << endl;
 348     }
 349
 350     // since each candidate good resonance is taken once, we must check
 351     // that it matches the pair definition in any order, and reject in case
 352     // in none of them the pair is OK
 353     // anyway, we must associate the correct daughter to the correct data member
 354     if (pairDefMatch[0] == 0 && pairDefMatch[1] == 1)
 355     {
 356       // 1st track --> 1st member of PairDef
 357       fDaughter[0].SetRef(mc->GetTrack(label[0]));
 358       fDaughter[0].SetRefMC((AliMCParticle*)mc->GetTrack(label[0]));
 359       fDaughter[0].SetGood();
 360       // 2nd track --> 2nd member of PairDef
 361       fDaughter[1].SetRef(mc->GetTrack(label[1]));
 362       fDaughter[1].SetRefMC((AliMCParticle*)mc->GetTrack(label[1]));
 363       fDaughter[1].SetGood();
 364     }
 365     else if ((pairDefMatch[0] == 1 && pairDefMatch[1] == 0))
 366     {
 367       // 1st track --> 2nd member of PairDef
 368       fDaughter[0].SetRef(mc->GetTrack(label[1]));
 369       fDaughter[0].SetRefMC((AliMCParticle*)mc->GetTrack(label[1]));
 370       fDaughter[0].SetGood();
 371       // 2nd track --> 1st member of PairDef
 372       fDaughter[1].SetRef(mc->GetTrack(label[0]));
 373       fDaughter[1].SetRefMC((AliMCParticle*)mc->GetTrack(label[0]));
 374       fDaughter[1].SetGood();
 375     }
 376     else
 377     {
 378       fDaughter[0].SetBad();
 379       fDaughter[1].SetBad();
 380     }
 381
 382     // reject the pair if the matching was unsuccessful
 383     if (!fDaughter[0].IsOK() || !fDaughter[1].IsOK()) continue;
 384
 385     // first, we set the internal AliRsnMother object to
 386     // the MC particles and then operate the selections on MC
 387     fPair.SetDaughters(&fDaughter[0], pairDef->GetMass(0), &fDaughter[1], pairDef->GetMass(1));
 388     FillContainer(cont, &fStepListMC, pairDef, 0);
 389
 390     // then, if both particles found a good match in the ESD
 391     // reassociate the ESD tracks to the pair and fill ESD containers
 392     if (esdIndex[0] < 0 || esdIndex[1] < 0) continue;
 393     if (pairDefMatch[0] == 0 && pairDefMatch[1] == 1)
 394     {
 395       // 1st track --> 1st member of PairDef
 396       fDaughter[0].SetRef(esd->GetTrack(esdIndex[0]));
 397       // 2nd track --> 2nd member of PairDef
 398       fDaughter[1].SetRef(esd->GetTrack(esdIndex[1]));
 399       //cout << "****** MATCHING SCHEME 1" << endl;
 400     }
 401     else if ((pairDefMatch[0] == 1 && pairDefMatch[1] == 0))
 402     {
 403       // 1st track --> 2nd member of PairDef
 404       fDaughter[0].SetRef(esd->GetTrack(esdIndex[1]));
 405       // 2nd track --> 1st member of PairDef
 406       fDaughter[1].SetRef(esd->GetTrack(esdIndex[0]));
 407       //cout << "****** MATCHING SCHEME 2" << endl;
 408     }
 409     //cout << "****** IDs = " << fDaughter[0].GetID() << ' ' << fDaughter[1].GetID() << endl;
 410     // here we must remember how many steps were already filled
 411     first = (Int_t)fStepListMC.GetEntries();
 412     FillContainer(cont, &fStepListESD, pairDef, first);
 413   }
 414 }
 415
 416 //_____________________________________________________________________________
 417 void AliRsnAnalysisEffSE::ProcessEventAOD(AliRsnPairDef *pairDef)
 418 {
 419 //
 420 // Process current event with the definitions of the specified step in MC list
 421 // and store results in the container slot defined in second argument.
 422 // It is associated with the AliCFContainer with the name of the pair.
 423 //
 424
 425   AliAODEvent      *aod   = fRsnEvent.GetRefAOD();
 426   AliAODMCParticle *mother;
 427   TClonesArray     *mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 428   if (!mcArray) return;
 429
 430   if (!pairDef) return;
 431   AliCFContainer *cont = (AliCFContainer*)fContainerList->FindObject(pairDef->GetPairName());
 432   if (!cont) return;
 433
 434   // get informations from pairDef
 435   Int_t   pdgM = 0;
 436   Int_t   pdgD[2] = {0, 0};
 437   Short_t chargeD[2] = {0, 0};
 438   pdgM       = pairDef->GetMotherPDG();
 439   pdgD   [0] = AliPID::ParticleCode(pairDef->GetPID(0));
 440   pdgD   [1] = AliPID::ParticleCode(pairDef->GetPID(1));
 441   chargeD[0] = pairDef->GetChargeShort(0);
 442   chargeD[1] = pairDef->GetChargeShort(1);
 443
 444   // other utility variables
 445   Int_t             first, j;
 446   Int_t             label [2] = {-1, -1};
 447   Short_t           charge[2] = {0, 0};
 448   Short_t           pairDefMatch[2] = {-1, -1};
 449   Int_t             aodIndex[2] = {-1, -1};
 450   AliAODMCParticle *part[2] = {0, 0};
 451
 452   // loop on MC particles
 453   TObjArrayIter next(mcArray);
 454   while ( (mother = (AliAODMCParticle*)next()) )
 455   {
 456     // check that it is a binary decay and the PDG code matches
 457     if (mother->GetNDaughters() != 2) continue;
 458     if (mother->GetPdgCode() != pdgM) continue;
 459
 460     // store the labels of the two daughters
 461     label[0] = mother->GetDaughter(0);
 462     label[1] = mother->GetDaughter(1);
 463
 464     // retrieve the particles and other stuff
 465     // check if they match the order in the pairDef
 466     for (j = 0; j < 2; j++)
 467     {
 468       if (label[j] < 0) continue;
 469       part[j]   = (AliAODMCParticle*)mcArray->At(label[j]);
 470       pdgD[j]   = TMath::Abs(part[j]->GetPdgCode());
 471       charge[j] = (Short_t)(part[j]->Charge());
 472       if (pdgD[j] == AliPID::ParticleCode(pairDef->GetPID(0)) && charge[j] == pairDef->GetChargeShort(0))
 473         pairDefMatch[j] = 0;
 474       else if (pdgD[j] == AliPID::ParticleCode(pairDef->GetPID(1)) && charge[j] == pairDef->GetChargeShort(1))
 475         pairDefMatch[j] = 1;
 476       else
 477         pairDefMatch[j] = -1;
 478
 479       // find corresponding ESD particle: first try rejecting fakes,
 480       // and in case of failure, try accepting fakes
 481       aodIndex[j] = FindAODtrack(label[j], aod, kTRUE);
 482       if (aodIndex[j] < 0) aodIndex[j] = FindAODtrack(label[j], aod, kFALSE);
 483     }
 484
 485     // since each candidate good resonance is taken once, we must check
 486     // that it matches the pair definition in any order, and reject in case
 487     // in none of them the pair is OK
 488     // anyway, we must associate the correct daughter to the correct data member
 489     if (pairDefMatch[0] == 0 && pairDefMatch[1] == 1)
 490     {
 491       // 1st track --> 1st member of PairDef
 492       fDaughter[0].SetRef((AliAODMCParticle*)mcArray->At(label[0]));
 493       fDaughter[0].SetRefMC((AliAODMCParticle*)mcArray->At(label[0]));
 494       fDaughter[0].SetGood();
 495       // 2nd track --> 2nd member of PairDef
 496       fDaughter[1].SetRef((AliAODMCParticle*)mcArray->At(label[1]));
 497       fDaughter[1].SetRefMC((AliAODMCParticle*)mcArray->At(label[1]));
 498       fDaughter[1].SetGood();
 499     }
 500     else if ((pairDefMatch[0] == 1 && pairDefMatch[1] == 0))
 501     {
 502       // 1st track --> 2nd member of PairDef
 503       fDaughter[0].SetRef((AliAODMCParticle*)mcArray->At(label[1]));
 504       fDaughter[0].SetRefMC((AliAODMCParticle*)mcArray->At(label[1]));
 505       fDaughter[0].SetGood();
 506       // 2nd track --> 1st member of PairDef
 507       fDaughter[1].SetRef((AliAODMCParticle*)mcArray->At(label[0]));
 508       fDaughter[1].SetRefMC((AliAODMCParticle*)mcArray->At(label[0]));
 509       fDaughter[1].SetGood();
 510     }
 511     else
 512     {
 513       fDaughter[0].SetBad();
 514       fDaughter[1].SetBad();
 515     }
 516
 517     // reject the pair if the matching was unsuccessful
 518     if (!fDaughter[0].IsOK() || !fDaughter[1].IsOK()) continue;
 519
 520     // first, we set the internal AliRsnMother object to
 521     // the MC particles and then operate the selections on MC
 522     fPair.SetDaughters(&fDaughter[0], pairDef->GetMass(0), &fDaughter[1], pairDef->GetMass(1));
 523     FillContainer(cont, &fStepListMC, pairDef, 0);
 524
 525     // then, if both particles found a good match in the AOD
 526     // reassociate the AOD tracks to the pair and fill AOD containers
 527     if (aodIndex[0] < 0 || aodIndex[1] < 0) continue;
 528     if (pairDefMatch[0] == 0 && pairDefMatch[1] == 1)
 529     {
 530       // 1st track --> 1st member of PairDef
 531       fDaughter[0].SetRef(aod->GetTrack(aodIndex[0]));
 532       // 2nd track --> 2nd member of PairDef
 533       fDaughter[1].SetRef(aod->GetTrack(aodIndex[1]));
 534       //cout << "****** MATCHING SCHEME 1" << endl;
 535     }
 536     else if ((pairDefMatch[0] == 1 && pairDefMatch[1] == 0))
 537     {
 538       // 1st track --> 2nd member of PairDef
 539       fDaughter[0].SetRef(aod->GetTrack(aodIndex[1]));
 540       // 2nd track --> 1st member of PairDef
 541       fDaughter[1].SetRef(aod->GetTrack(aodIndex[0]));
 542       //cout << "****** MATCHING SCHEME 2" << endl;
 543     }
 544     //cout << "****** IDs = " << fDaughter[0].GetID() << ' ' << fDaughter[1].GetID() << endl;
 545     // here we must remember how many steps were already filled
 546     first = (Int_t)fStepListMC.GetEntries();
 547     FillContainer(cont, &fStepListESD, pairDef, first);
 548   }
 549 }
 550
 551 //_____________________________________________________________________________
 552 void AliRsnAnalysisEffSE::FillContainer(AliCFContainer *cont, const TObjArray *stepList, AliRsnPairDef *pd, Int_t firstOutStep)
 553 {
 554 //
 555 // Fill the containers
 556 //
 557
 558   Int_t iaxis, nAxes  = fAxisList.GetEntries();
 559   Int_t istep, nSteps = stepList->GetEntries();
 560
 561   // set daughters to pair
 562   fPair.SetDaughters(&fDaughter[0], pd->GetMass(0), &fDaughter[1], pd->GetMass(1));
 563
 564   // compute values for all axes
 565   for (iaxis = 0; iaxis < nAxes; iaxis++)
 566   {
 567     AliRsnValue *fcnAxis = (AliRsnValue*)fAxisList.At(iaxis);
 568     fVar[iaxis] = -1E10;
 569     if (fcnAxis->Eval(&fPair, pd, &fRsnEvent)) fVar[iaxis] = (Float_t)fcnAxis->GetValue();
 570   }
 571
 572   // fill all steps
 573   for (istep = 0; istep < nSteps; istep++)
 574   {
 575     AliRsnCutManager *cutMgr = (AliRsnCutManager*)stepList->At(istep);
 576     cutMgr->SetEvent(&fRsnEvent);
 577     if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
 578     if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
 579     if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
 580     if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
 581     if (!cutMgr->PassMotherCuts(&fPair)) break;
 582     //cout << "**************************************** FILLING STEP " << istep << endl;
 583     cont->Fill(fVar.GetArray(), istep + firstOutStep);
 584   }
 585 }
 586
 587 //_____________________________________________________________________________
 588 void AliRsnAnalysisEffSE::RsnTerminate(Option_t*)
 589 {
 590 //
 591 // Termination.
 592 // Could be added some monitor histograms here.
 593 //
 594
 595   AliDebug(AliLog::kDebug+2,"<-");
 596   AliDebug(AliLog::kDebug+2,"->");
 597 }
 598
 599 //_____________________________________________________________________________
 600 void AliRsnAnalysisEffSE::AddPairDef(AliRsnPairDef* pairDef)
 601 {
 602 //
 603 //  Adds pair definition
 604 //
 605   fPairDefList.AddLast(pairDef);
 606 }
 607
 608 //_____________________________________________________________________________
 609 Int_t AliRsnAnalysisEffSE::FindESDtrack(Int_t label, AliESDEvent *esd, Bool_t rejectFakes)
 610 {
 611 //
 612 // Finds in the ESD a track whose label corresponds to that in argument.
 613 // When global tracks are enabled, tries first to find a global track
 614 // satisfying that requirement.
 615 // If no global tracks are found, if ITS-SA are enable, tries to search among them
 616 // otherwise return a negative number.
 617 // If global tracks are disabled, search only among ITS SA
 618 //
 619
 620   Int_t   i = 0;
 621   Int_t   ntracks = esd->GetNumberOfTracks();
 622   ULong_t status;
 623   Bool_t  isTPC;
 624   Bool_t  isITSSA;
 625
 626   // loop for global tracks
 627   if (fUseGlobal)
 628   {
 629     for (i = 0; i < ntracks; i++)
 630     {
 631       AliESDtrack *track = esd->GetTrack(i);
 632       status  = (ULong_t)track->GetStatus();
 633       isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
 634       if (!isTPC) continue;
 635
 636       // check that label match
 637       if (TMath::Abs(track->GetLabel()) != label) continue;
 638
 639       // if required, reject fakes
 640       if (rejectFakes && track->GetLabel() < 0) continue;
 641
 642       // if all checks are passed and we are searching among global
 643       // this means that thie track is a global one with the right label
 644       // then, the return value is set to this, and returned
 645       return i;
 646     }
 647   }
 648
 649   // loop for ITS-SA tracks (this happens only if no global tracks were found
 650   // or searching among globals is disabled)
 651   if (fUseITSSA)
 652   {
 653     for (i = 0; i < ntracks; i++)
 654     {
 655       AliESDtrack *track = esd->GetTrack(i);
 656       status  = (ULong_t)track->GetStatus();
 657       isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 658       if (!isITSSA) continue;
 659
 660       // check that label match
 661       if (TMath::Abs(track->GetLabel()) != label) continue;
 662
 663       // if required, reject fakes
 664       if (rejectFakes && track->GetLabel() < 0) continue;
 665
 666       // if all checks are passed and we are searching among global
 667       // this means that thie track is a global one with the right label
 668       // then, the return value is set to this, and returned
 669       return i;
 670     }
 671   }
 672
 673   // if we reach this point, no match were found
 674   return -1;
 675 }
 676
 677 //_____________________________________________________________________________
 678 TArrayI AliRsnAnalysisEffSE::FindESDtracks(Int_t label, AliESDEvent *esd)
 679 {
 680 //
 681 // Finds in the ESD a track whose label corresponds to that in argument.
 682 // When global tracks are enabled, tries first to find a global track
 683 // satisfying that requirement.
 684 // If no global tracks are found, if ITS-SA are enable, tries to search among them
 685 // otherwise return a negative number.
 686 // If global tracks are disabled, search only among ITS SA
 687 //
 688
 689   Int_t   i = 0;
 690   Int_t   ntracks = esd->GetNumberOfTracks();
 691   ULong_t status;
 692   Bool_t  isTPC;
 693   Bool_t  isITSSA;
 694   TArrayI array(100);
 695   Int_t   nfound = 0;
 696
 697   // loop for global tracks
 698   if (fUseGlobal)
 699   {
 700     for (i = 0; i < ntracks; i++)
 701     {
 702       AliESDtrack *track = esd->GetTrack(i);
 703       status  = (ULong_t)track->GetStatus();
 704       isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
 705       if (!isTPC) continue;
 706
 707       // check that label match
 708       if (TMath::Abs(track->GetLabel()) != label) continue;
 709
 710       array[nfound++] = i;
 711     }
 712   }
 713
 714   // loop for ITS-SA tracks (this happens only if no global tracks were found
 715   // or searching among globals is disabled)
 716   if (fUseITSSA)
 717   {
 718     for (i = 0; i < ntracks; i++)
 719     {
 720       AliESDtrack *track = esd->GetTrack(i);
 721       status  = (ULong_t)track->GetStatus();
 722       isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 723       if (!isITSSA) continue;
 724
 725       // check that label match
 726       if (TMath::Abs(track->GetLabel()) != label) continue;
 727
 728       array[nfound++] = i;
 729     }
 730   }
 731
 732   array.Set(nfound);
 733   return array;
 734 }
 735
 736 //_____________________________________________________________________________
 737 Int_t AliRsnAnalysisEffSE::FindAODtrack(Int_t label, AliAODEvent *aod, Bool_t rejectFakes)
 738 {
 739 //
 740 // Finds in the ESD a track whose label corresponds to that in argument.
 741 // When global tracks are enabled, tries first to find a global track
 742 // satisfying that requirement.
 743 // If no global tracks are found, if ITS-SA are enable, tries to search among them
 744 // otherwise return a negative number.
 745 // If global tracks are disabled, search only among ITS SA
 746 //
 747
 748   Int_t   i = 0;
 749   Int_t   ntracks = aod->GetNumberOfTracks();
 750   ULong_t status;
 751   Bool_t  isTPC;
 752   Bool_t  isITSSA;
 753
 754   // loop for global tracks
 755   if (fUseGlobal)
 756   {
 757     for (i = 0; i < ntracks; i++)
 758     {
 759       AliAODTrack *track = aod->GetTrack(i);
 760       status  = (ULong_t)track->GetStatus();
 761       isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
 762       if (!isTPC) continue;
 763
 764       // check that label match
 765       if (TMath::Abs(track->GetLabel()) != label) continue;
 766
 767       // if required, reject fakes
 768       if (rejectFakes && track->GetLabel() < 0) continue;
 769
 770       // if all checks are passed and we are searching among global
 771       // this means that thie track is a global one with the right label
 772       // then, the return value is set to this, and returned
 773       return i;
 774     }
 775   }
 776
 777   // loop for ITS-SA tracks (this happens only if no global tracks were found
 778   // or searching among globals is disabled)
 779   if (fUseITSSA)
 780   {
 781     for (i = 0; i < ntracks; i++)
 782     {
 783       AliAODTrack *track = aod->GetTrack(i);
 784       status  = (ULong_t)track->GetStatus();
 785       isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 786       if (!isITSSA) continue;
 787
 788       // check that label match
 789       if (TMath::Abs(track->GetLabel()) != label) continue;
 790
 791       // if required, reject fakes
 792       if (rejectFakes && track->GetLabel() < 0) continue;
 793
 794       // if all checks are passed and we are searching among global
 795       // this means that thie track is a global one with the right label
 796       // then, the return value is set to this, and returned
 797       return i;
 798     }
 799   }
 800
 801   // if we reach this point, no match were found
 802   return -1;
 803 }
 804
 805 //_____________________________________________________________________________
 806 TArrayI AliRsnAnalysisEffSE::FindAODtracks(Int_t label, AliAODEvent *aod)
 807 {
 808 //
 809 // Finds in the ESD a track whose label corresponds to that in argument.
 810 // When global tracks are enabled, tries first to find a global track
 811 // satisfying that requirement.
 812 // If no global tracks are found, if ITS-SA are enable, tries to search among them
 813 // otherwise return a negative number.
 814 // If global tracks are disabled, search only among ITS SA
 815 //
 816
 817   Int_t   i = 0;
 818   Int_t   ntracks = aod->GetNumberOfTracks();
 819   ULong_t status;
 820   Bool_t  isTPC;
 821   Bool_t  isITSSA;
 822   TArrayI array(100);
 823   Int_t   nfound = 0;
 824
 825   // loop for global tracks
 826   if (fUseGlobal)
 827   {
 828     for (i = 0; i < ntracks; i++)
 829     {
 830       AliAODTrack *track = aod->GetTrack(i);
 831       status  = (ULong_t)track->GetStatus();
 832       isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
 833       if (!isTPC) continue;
 834
 835       // check that label match
 836       if (TMath::Abs(track->GetLabel()) != label) continue;
 837
 838       array[nfound++] = i;
 839     }
 840   }
 841
 842   // loop for ITS-SA tracks (this happens only if no global tracks were found
 843   // or searching among globals is disabled)
 844   if (fUseITSSA)
 845   {
 846     for (i = 0; i < ntracks; i++)
 847     {
 848       AliAODTrack *track = aod->GetTrack(i);
 849       status  = (ULong_t)track->GetStatus();
 850       isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 851       if (!isITSSA) continue;
 852
 853       // check that label match
 854       if (TMath::Abs(track->GetLabel()) != label) continue;
 855
 856       array[nfound++] = i;
 857     }
 858   }
 859
 860   array.Set(nfound);
 861   return array;
 862 }
 863