PWG2/RESONANCES/AliRsnReader.cxx

   1 //
   2 // Class AliRsnReader
   3 //
   4 // This is the universal converter from any kind of source event
   5 // (i.e. ESD, standard AOD, MC) into the internal non-standard
   6 // AOD format used by RSN package.
   7 // ---
   8 // This class reads all tracks in an input event and converts them
   9 // into AliRsnDaughters, and computes preliminarily the PID probabilities
  10 // by doing the Bayesian combination of a set of assigned prior probabilities
  11 // with the PID weights defined in each track.
  12 // ---
  13 // When filling the output event (AliRsnEvent), some arrays of indexes
  14 // are created in order to organize tracks according to their PID and charge,
  15 // which will then be used in further analysis steps.
  16 //
  17 // author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
  18 //
  19
  20 #include <Riostream.h>
  21 #include <TString.h>
  22
  23 #include "AliLog.h"
  24
  25 #include "AliVEvent.h"
  26
  27 #include "AliESDEvent.h"
  28 #include "AliESDtrack.h"
  29 #include "AliESDVertex.h"
  30 #include "AliESDtrackCuts.h"
  31 #include "AliMultiplicity.h"
  32
  33 #include "AliAODEvent.h"
  34 #include "AliAODTrack.h"
  35 #include "AliAODVertex.h"
  36
  37 #include "AliStack.h"
  38 #include "AliMCEvent.h"
  39 #include "AliMCParticle.h"
  40 #include "AliGenEventHeader.h"
  41
  42 #include "AliRsnMCInfo.h"
  43 #include "AliRsnDaughter.h"
  44 #include "AliRsnEvent.h"
  45 #include "AliRsnPIDDefESD.h"
  46 #include "AliRsnCut.h"
  47
  48 #include "AliRsnReader.h"
  49
  50 ClassImp(AliRsnReader)
  51
  52 //_____________________________________________________________________________
  53 AliRsnReader::AliRsnReader() :
  54   fCheckSplit(kFALSE),
  55   fRejectFakes(kFALSE),
  56   fTPCOnly(kFALSE),
  57   fUseESDTrackCuts(kFALSE),
  58   fUseRsnTrackCuts(kFALSE),
  59   fCheckVertexStatus(kFALSE),
  60   fMinNContributors(0),
  61   fPID(),
  62   fPIDDef(),
  63   fPIDArraysSize(1000),
  64   fITSClusters(0),
  65   fTPCClusters(0),
  66   fTRDClusters(0),
  67   fTrackRefs(0),
  68   fTrackRefsITS(0),
  69   fTrackRefsTPC(0),
  70   fESDTrackCuts(),
  71   fRsnTrackCuts("primaryCuts")
  72 {
  73 //
  74 // Constructor.
  75 //
  76 }
  77
  78 //_____________________________________________________________________________
  79 Bool_t AliRsnReader::AreSplitted(AliESDtrack *track1, AliESDtrack *track2)
  80 {
  81 //
  82 // Checks if two tracks are splitted.
  83 // Currently, two split tracks are the ones which
  84 // have equal GEANT labels.
  85 // On real data, this criterion will have to change
  86 // into an "experimental" one.
  87 //
  88
  89   Int_t lab1 = TMath::Abs(track1->GetLabel());
  90   Int_t lab2 = TMath::Abs(track2->GetLabel());
  91
  92   return (lab1 == lab2);
  93 }
  94
  95 //_____________________________________________________________________________
  96 Bool_t AliRsnReader::ResolveSplit(AliESDtrack *track1, AliESDtrack *track2)
  97 {
  98 //
  99 // Resolves a split of two tracks.
 100 // Only two values can be returned:
 101 //  kTRUE  --> accept "track1" and reject "track2"
 102 //  kFALSE --> accept "track2" and reject "track1"
 103 //
 104
 105   Double_t chiSq1 = track1->GetConstrainedChi2();
 106   Double_t chiSq2 = track2->GetConstrainedChi2();
 107
 108   if (chiSq1 < chiSq2) return 1; else return 2;
 109 }
 110
 111 //_____________________________________________________________________________
 112 void AliRsnReader::SetTPCOnly(Bool_t doit)
 113 {
 114 //
 115 // Set the flag for TPC only.
 116 // If this is true, exclude all other detectors from PID
 117 //
 118
 119   fTPCOnly = doit;
 120
 121   if (fTPCOnly) {
 122     fPIDDef.ExcludeAll();
 123     fPIDDef.IncludeDet(AliRsnPIDDefESD::kTPC);
 124     fPIDDef.SetDivValue(AliRsnPIDDefESD::kTPC, 0.0);
 125   }
 126 }
 127
 128 //_____________________________________________________________________________
 129 Bool_t AliRsnReader::ConvertTrack(AliRsnDaughter *daughter, AliESDtrack *esdTrack)
 130 {
 131 //
 132 // Translates an ESD track into a RSN daughter
 133 //
 134
 135   if (!esdTrack) {
 136     AliError("Passed NULL object: nothing done");
 137     return kFALSE;
 138   }
 139
 140   Double_t p[3], v[3], pid[AliRsnPID::kSpecies];
 141
 142   // copy values which don't need treatment:
 143   // momentum, vertex, chi2, flags, charge, number of TPC and ITS clusters
 144   esdTrack->GetPxPyPz(p);
 145   daughter->SetP(p[0], p[1], p[2]);
 146   esdTrack->GetXYZ(v);
 147   daughter->SetV(v[0], v[1], v[2]);
 148   daughter->SetChi2(esdTrack->GetConstrainedChi2());
 149   daughter->SetFlags(esdTrack->GetStatus());
 150   daughter->SetCharge((Short_t)esdTrack->Charge());
 151   daughter->SetNumberOfITSClusters(esdTrack->GetITSclusters(0x0));
 152   daughter->SetNumberOfTPCClusters(esdTrack->GetTPCclusters(0x0));
 153
 154   // define the kink index:
 155   //  0 = no kink
 156   //  1 = kink daughter
 157   // -1 = kink mother
 158   Int_t i, ik[3];
 159   for (i = 0; i < 3; i++) ik[i] = esdTrack->GetKinkIndex(i);
 160   if (ik[0] < 0 || ik[1] < 0 || ik[2] < 0) daughter->SetKinkMother();
 161   else if (ik[0] > 0 || ik[1] > 0 || ik[2] > 0) daughter->SetKinkDaughter();
 162   else daughter->SetNoKink();
 163
 164   // store PID weights according to definition
 165   // check: if the sum of all weights is null, the adoption fails
 166   Double_t sum = 0.0;
 167   fPIDDef.ComputeWeights(esdTrack, pid);
 168   for (i = 0; i < AliRsnPID::kSpecies; i++)
 169   {
 170     daughter->SetPIDWeight(i, pid[i]);
 171     sum += pid[i];
 172   }
 173   if (sum <= 0.0) return kFALSE;
 174
 175   // compute probabilities
 176   if (!fPID.ComputeProbs(daughter)) return kFALSE;
 177
 178   // calculate N sigma to vertex
 179   AliESDtrackCuts trkCut;
 180   daughter->SetNSigmaToVertex(trkCut.GetSigmaToVertex(esdTrack));
 181
 182   return kTRUE;
 183 }
 184
 185 //_____________________________________________________________________________
 186 Bool_t AliRsnReader::ConvertTrack(AliRsnDaughter *daughter, AliAODTrack *aodTrack)
 187 {
 188 //
 189 // Translates an AOD track into a RSN daughter
 190 //
 191
 192   if (!aodTrack)
 193   {
 194     AliError("Passed NULL object: nothing done");
 195     return kFALSE;
 196   }
 197
 198   // copy momentum  and vertex
 199   daughter->SetP(aodTrack->Px(), aodTrack->Py(), aodTrack->Pz());
 200   daughter->SetV(aodTrack->Xv(), aodTrack->Yv(), aodTrack->Zv());
 201
 202   // chi2
 203   daughter->SetChi2(aodTrack->Chi2perNDF());
 204
 205   // copy PID weights
 206   Int_t i;
 207   for (i = 0; i < 5; i++) daughter->SetPIDWeight(i, aodTrack->PID()[i]);
 208
 209   // compute probabilities
 210   if (!fPID.ComputeProbs(daughter)) return kFALSE;
 211
 212   // copy flags
 213   daughter->SetFlags(aodTrack->GetStatus());
 214
 215   // copy sign
 216   daughter->SetCharge(aodTrack->Charge());
 217
 218   return kTRUE;
 219 }
 220
 221 //_____________________________________________________________________________
 222 Bool_t AliRsnReader::ConvertTrack(AliRsnDaughter *daughter, TParticle *particle)
 223 {
 224 //
 225 // Translates an MC particle into a RSN daughter
 226 //
 227
 228   if (!particle) return kFALSE;
 229
 230   // copy other MC info (mother PDG code cannot be retrieved here)
 231   daughter->InitMCInfo(particle);
 232
 233   // copy momentum  and vertex
 234   daughter->SetP(particle->Px(), particle->Py(), particle->Pz());
 235   daughter->SetV(particle->Vx(), particle->Vy(), particle->Vz());
 236
 237   // recognize charge sign from PDG code sign
 238   Int_t pdg = particle->GetPdgCode();
 239   Int_t absPDG = TMath::Abs(pdg);
 240   if (absPDG == 11 || absPDG == 13)
 241   {
 242     if (pdg > 0) daughter->SetCharge(-1); else daughter->SetCharge(1);
 243   }
 244   else if (absPDG == 211 || absPDG == 321 || absPDG == 2212)
 245   {
 246     if (pdg > 0) daughter->SetCharge(1); else daughter->SetCharge(-1);
 247   }
 248   else
 249   {
 250     // when trying to "adopt" a neutral track (photon, neutron, etc.)
 251     // for the moment a "failed" message is returned
 252     return kFALSE;
 253   }
 254
 255   // flags and PID weights make no sense with MC tracks
 256   daughter->SetFlags(0);
 257   for (pdg = 0; pdg < AliRsnPID::kSpecies; pdg++) daughter->SetPIDWeight(pdg, 0.0);
 258   daughter->SetPIDWeight(AliRsnPID::InternalType(absPDG), 1.0);
 259
 260   // compute probabilities
 261   if (!fPID.ComputeProbs(daughter)) return kFALSE;
 262
 263   return kTRUE;
 264 }
 265
 266 //_____________________________________________________________________________
 267 Bool_t AliRsnReader::Fill
 268 (AliRsnEvent *rsn, AliVEvent *event, AliMCEvent *mc)
 269 {
 270 //
 271 // According to the class type of event and the selected source
 272 // recalls one of the private reading methods to fill the RsnEvent
 273 // passed by reference as first argument.
 274 //
 275
 276   Bool_t success = kFALSE;
 277   TString str(event->ClassName());
 278
 279   if (!str.CompareTo("AliESDEvent")) {
 280     AliDebug(1, "Reading an ESD event");
 281     success = FillFromESD(rsn, (AliESDEvent*)event, mc);
 282   }
 283   else if (!str.CompareTo("AliAODEvent")) {
 284     AliDebug(1, "Reading an AOD event");
 285     success = FillFromAOD(rsn, (AliAODEvent*)event, mc);
 286   }
 287   else if (!str.CompareTo("AliMCEvent")) {
 288     AliDebug(1, "Reading an MC event");
 289     success = FillFromMC(rsn, (AliMCEvent*)event);
 290   }
 291   else {
 292     AliError(Form("ClassName '%s' not recognized as possible source data: aborting.", str.Data()));
 293     return kFALSE;
 294   }
 295
 296   return success;
 297 }
 298
 299 //_____________________________________________________________________________
 300 Bool_t AliRsnReader::FillFromESD(AliRsnEvent *rsn, AliESDEvent *esd, AliMCEvent *mc)
 301 {
 302 //
 303 // Filler from an ESD event.
 304 // Stores all tracks (if a filter is defined, it will store
 305 // only the ones which survive the cuts).
 306 // If a reference MC event is provided, it is used to store
 307 // the MC informations for each track (true PDG code,
 308 // GEANT label of mother, PDG code of mother, if any).
 309 // When this is used, the 'source' flag of the output
 310 // AliRsnEvent object will be set to 'kESD'.
 311 //
 312
 313   // retrieve stack (if possible)
 314   AliStack *stack = 0x0;
 315   if (mc) stack = mc->Stack();
 316
 317   // set MC primary vertex if present
 318   TArrayF fvertex(3);
 319   Double_t mcvx = 0., mcvy = 0., mcvz = 0.;
 320   if (mc) {
 321     mc->GenEventHeader()->PrimaryVertex(fvertex);
 322     mcvx = (Double_t)fvertex[0];
 323     mcvy = (Double_t)fvertex[1];
 324     mcvz = (Double_t)fvertex[2];
 325     rsn->SetPrimaryVertexMC(mcvx, mcvy, mcvz);
 326   }
 327
 328   // get number of tracks
 329   Int_t ntracks = esd->GetNumberOfTracks();
 330   if (!ntracks) return kFALSE;
 331
 332   // set SPD multiplicity
 333   const AliMultiplicity* SPDMult = esd->GetMultiplicity();
 334   Int_t nTracklets = SPDMult->GetNumberOfTracklets();
 335   rsn->SetTrueMultiplicity(nTracklets);
 336
 337   // if required with the flag, scans the event
 338   // and searches all split tracks (= 2 tracks with the same label);
 339   // for each pair of split tracks, only the better (best chi2) is kept
 340   // and the other is rejected: this info is stored into a Boolean array
 341   Int_t i, i1, i2;
 342   Bool_t *accept = new Bool_t[ntracks];
 343   for (i = 0; i < ntracks; i++) accept[i] = kTRUE;
 344   if (fCheckSplit) {
 345     for (i1 = 0; i1 < ntracks; i1++) {
 346       AliESDtrack *track1 = esd->GetTrack(i1);
 347       for (i2 = i1 + 1; i2 < ntracks; i2++) {
 348         AliESDtrack *track2 = esd->GetTrack(i2);
 349         if (AreSplitted(track1, track2)) {
 350           if (ResolveSplit(track1, track2)) accept[i2] = kFALSE;
 351           else accept[i1] = kFALSE;
 352         }
 353       }
 354     }
 355   }
 356
 357   // get primary vertex
 358   Double_t vertex[3];
 359   if (!fTPCOnly) {
 360     // when taking vertex from ESD event there are two options:
 361     // if a vertex with tracks was successfully reconstructed,
 362     // it is used for computing DCA;
 363     // otherwise, the one computed with SPD is used.
 364     // This is known from the "Status" parameter of the vertex itself.
 365     const AliESDVertex *v = esd->GetPrimaryVertex();
 366     if (!v->GetStatus()) v = esd->GetPrimaryVertexSPD();
 367     if (fCheckVertexStatus && (v->GetStatus() == kFALSE)) return kFALSE;
 368     if (fMinNContributors > 0 && v->GetNContributors() < fMinNContributors) return kFALSE;
 369
 370     // get primary vertex
 371     vertex[0] = (Double_t)v->GetXv();
 372     vertex[1] = (Double_t)v->GetYv();
 373     vertex[2] = (Double_t)v->GetZv();
 374   }
 375   else {
 376     const AliESDVertex *v = esd->GetPrimaryVertexTPC();
 377     if (fCheckVertexStatus && v->GetStatus() == kFALSE) return kFALSE;
 378     if (fMinNContributors > 0 && v->GetNContributors() < fMinNContributors) return kFALSE;
 379
 380     // get primary vertex
 381     vertex[0] = (Double_t)v->GetXv();
 382     vertex[1] = (Double_t)v->GetYv();
 383     vertex[2] = (Double_t)v->GetZv();
 384   }
 385   rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
 386
 387   // store tracks from ESD
 388   Float_t p[2], cov[3];
 389   Int_t   index, label, labmum;
 390   Bool_t  check;
 391   AliRsnDaughter temp;
 392   AliESDtrack *esdTrack = 0x0, esdTrackTmp;
 393   for (index = 0; index < ntracks; index++) {
 394     // skip track recognized as the worse one in a splitted pair
 395     if (!accept[index]) {
 396       AliDebug(10, Form("Rejecting split track #%d in this event", index));
 397       continue;
 398     }
 399     // get ESD track
 400     esdTrack = esd->GetTrack(index);
 401     // check for ESD track cuts (if required)
 402     if (fUseESDTrackCuts && (!fESDTrackCuts.AcceptTrack(esdTrack))) continue;
 403     // check for fake tracks
 404     label = esdTrack->GetLabel();
 405     if (fRejectFakes && (label < 0)) continue;
 406     // check for minimum number of clusters in ITS, TPC and TRD
 407     if (fITSClusters > 0 && (esdTrack->GetITSclusters(0x0) < fITSClusters)) if (!fTPCOnly) continue;
 408     if (fTPCClusters > 0 && (esdTrack->GetTPCclusters(0x0) < fTPCClusters)) continue;
 409     if (fTRDClusters > 0 && (esdTrack->GetTRDclusters(0x0) < fTRDClusters)) if (!fTPCOnly) continue;
 410     // switch to TPC data if required
 411     if (fTPCOnly) {
 412       esdTrack->GetImpactParametersTPC(p, cov);
 413       if (p[0] == 0. && p[1] == 0.) {
 414         esdTrack->RelateToVertexTPC(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), kVeryBig);
 415       }
 416       check = esdTrack->FillTPCOnlyTrack(esdTrackTmp);
 417       if (check) esdTrack = &esdTrackTmp; else continue;
 418     }
 419     // try to get information from this track
 420     // output value tells if this was successful or not
 421     check = ConvertTrack(&temp, esdTrack);
 422     if (!check) {
 423       AliDebug(10, Form("Failed adopting track #%d", index));
 424       continue;
 425     }
 426
 427     // if stack is present, copy MC info
 428     if (stack) {
 429       TParticle *part = stack->Particle(TMath::Abs(label));
 430       if (part) {
 431         temp.InitMCInfo(part);
 432         labmum = part->GetFirstMother();
 433         if (labmum >= 0) {
 434           TParticle *mum = stack->Particle(labmum);
 435           temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 436         }
 437       }
 438     }
 439
 440     // set index and label
 441     temp.SetIndex(index);
 442     temp.SetLabel(label);
 443
 444     // check this object against the Rsn cuts (if required)
 445     if (fUseRsnTrackCuts) {
 446       if (!fRsnTrackCuts.IsSelected(AliRsnCut::kParticle, &temp)) continue;
 447     }
 448
 449     // try to add track to collection and returns an error in case of problems
 450     AliRsnDaughter *ptr = rsn->AddTrack(temp);
 451     if (!ptr) AliWarning(Form("Failed storing track#%d", index));
 452   }
 453
 454   // compute total multiplicity
 455   rsn->MakeComputations();
 456   if (rsn->GetMultiplicity() <= 0) {
 457     AliDebug(1, "Zero Multiplicity in this event");
 458     return kFALSE;
 459   }
 460
 461   // correct tracks for primary vertex
 462   //rsn->CorrectTracks();
 463
 464   // fill the PID index arrays
 465   rsn->FillPIDArrays(fPIDArraysSize);
 466
 467   return kTRUE;
 468 }
 469
 470 //_____________________________________________________________________________
 471 Bool_t AliRsnReader::FillFromAOD(AliRsnEvent *rsn, AliAODEvent *aod, AliMCEvent *mc)
 472 {
 473 //
 474 // Filler from an AOD event.
 475 // Stores all tracks (if a filter is defined, it will store
 476 // only the ones which survive the cuts).
 477 // If a reference MC event is provided, it is used to store
 478 // the MC informations for each track (true PDG code,
 479 // GEANT label of mother, PDG code of mother, if any).
 480 // When this is used, the 'source' flag of the output
 481 // AliRsnEvent object will be set to 'kAOD'.
 482 //
 483
 484   // retrieve stack (if possible)
 485   AliStack *stack = 0x0;
 486   if (mc) stack = mc->Stack();
 487
 488   // get number of tracks
 489   Int_t ntracks = aod->GetNTracks();
 490   if (!ntracks) return kFALSE;
 491   rsn->SetTrueMultiplicity(ntracks);
 492
 493   // get primary vertex
 494   Double_t vertex[3];
 495   vertex[0] = aod->GetPrimaryVertex()->GetX();
 496   vertex[1] = aod->GetPrimaryVertex()->GetY();
 497   vertex[2] = aod->GetPrimaryVertex()->GetZ();
 498   rsn->SetPrimaryVertex(vertex[0], vertex[1], vertex[2]);
 499
 500   // set MC primary vertex if present
 501   TArrayF fvertex(3);
 502   Double_t mcvx = 0., mcvy = 0., mcvz = 0.;
 503   if (mc) {
 504     mc->GenEventHeader()->PrimaryVertex(fvertex);
 505     mcvx = (Double_t)fvertex[0];
 506     mcvy = (Double_t)fvertex[1];
 507     mcvz = (Double_t)fvertex[2];
 508     rsn->SetPrimaryVertexMC(mcvx, mcvy, mcvz);
 509   }
 510
 511   // store tracks from ESD
 512   Int_t  index, label, labmum;
 513   Bool_t check;
 514   AliAODTrack *aodTrack = 0;
 515   AliRsnDaughter temp;
 516   TObjArrayIter iter(aod->GetTracks());
 517   while ((aodTrack = (AliAODTrack*)iter.Next()))
 518   {
 519     // retrieve index
 520     index = aod->GetTracks()->IndexOf(aodTrack);
 521     label = aodTrack->GetLabel();
 522     if (fRejectFakes && (label < 0)) continue;
 523     // copy ESD track data into RsnDaughter
 524     // if unsuccessful, this track is skipped
 525     check = ConvertTrack(&temp, aodTrack);
 526     if (!check) continue;
 527     // if stack is present, copy MC info
 528     if (stack)
 529     {
 530       TParticle *part = stack->Particle(TMath::Abs(label));
 531       if (part)
 532       {
 533         temp.InitMCInfo(part);
 534         labmum = part->GetFirstMother();
 535         if (labmum >= 0)
 536         {
 537           TParticle *mum = stack->Particle(labmum);
 538           temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 539         }
 540       }
 541     }
 542     // set index and label and add this object to the output container
 543     temp.SetIndex(index);
 544     temp.SetLabel(label);
 545
 546     // check this object against the Rsn cuts (if required)
 547     if (fUseRsnTrackCuts) {
 548       if (!fRsnTrackCuts.IsSelected(AliRsnCut::kParticle, &temp)) continue;
 549     }
 550
 551     AliRsnDaughter *ptr = rsn->AddTrack(temp);
 552     // if problems occurred while storin, that pointer is NULL
 553     if (!ptr) AliWarning(Form("Failed storing track#%d"));
 554   }
 555
 556   // compute total multiplicity
 557   rsn->MakeComputations();
 558   if (rsn->GetMultiplicity() <= 0)
 559   {
 560     AliDebug(1, "Zero multiplicity in this event");
 561     return kFALSE;
 562   }
 563
 564   // fill the PID index arrays
 565   rsn->FillPIDArrays(fPIDArraysSize);
 566
 567   // correct tracks for primary vertex
 568   rsn->CorrectTracks();
 569
 570   return kTRUE;
 571 }
 572
 573 //_____________________________________________________________________________
 574 Bool_t AliRsnReader::FillFromMC(AliRsnEvent *rsn, AliMCEvent *mc)
 575 {
 576 //
 577 // Filler from an ESD event.
 578 // Stores all tracks which generate at least one
 579 // TrackReference (a point in a sensitive volume).
 580 // In this case, the MC info is stored by default and
 581 // perfect particle identification is the unique available.
 582 // When this is used, the 'source' flag of the output
 583 // AliRsnEvent object will be set to 'kMC'.
 584 //
 585
 586   // get number of tracks
 587   Int_t ntracks = mc->GetNumberOfTracks();
 588   if (!ntracks) return kFALSE;
 589
 590   AliStack *stack = mc->Stack();
 591
 592   // get primary vertex
 593   TArrayF fvertex(3);
 594   Double_t vx, vy, vz;
 595   mc->GenEventHeader()->PrimaryVertex(fvertex);
 596   vx = (Double_t)fvertex[0];
 597   vy = (Double_t)fvertex[1];
 598   vz = (Double_t)fvertex[2];
 599   rsn->SetPrimaryVertex(vx, vy, vz);
 600   rsn->SetPrimaryVertexMC(vx, vy, vz);
 601   rsn->SetTrueMultiplicity(stack->GetNprimary());
 602
 603   // store tracks from MC
 604   Int_t  i, index, labmum, nRef, nRefITS, nRefTPC, detectorID;
 605   AliRsnDaughter temp;
 606   for (index = 0; index < ntracks; index++)
 607   {
 608     // get MC track & take index and label
 609     AliMCParticle *mcTrack = mc->GetTrack(index);
 610     temp.SetIndex(index);
 611     temp.SetLabel(mcTrack->Label());
 612
 613     // check particle track references
 614     nRef = mcTrack->GetNumberOfTrackReferences();
 615     for (i = 0, nRefITS = 0, nRefTPC = 0; i < nRef; i++) {
 616       AliTrackReference *trackRef = mcTrack->GetTrackReference(i);
 617       if (!trackRef) continue;
 618       detectorID = trackRef->DetectorId();
 619       switch (detectorID) {
 620         case AliTrackReference::kITS  : nRefITS++; break;
 621         case AliTrackReference::kTPC  : nRefTPC++; break;
 622         default: break;
 623       }
 624     }
 625     if (fTrackRefs > 0 && nRef < fTrackRefs) continue;
 626     if (fTrackRefsITS > 0 && nRefITS < fTrackRefsITS) continue;
 627     if (fTrackRefsTPC > 0 && nRefTPC < fTrackRefsTPC) continue;
 628
 629     // try to insert in the RsnDaughter its data
 630     TParticle *mcPart = mcTrack->Particle();
 631     if (!ConvertTrack(&temp, mcPart)) continue;
 632
 633     // assign mother label and PDG code (if any)
 634     labmum = temp.GetMCInfo()->Mother();
 635     if (labmum >= 0) {
 636       TParticle *mum = stack->Particle(labmum);
 637       temp.GetMCInfo()->SetMotherPDG(mum->GetPdgCode());
 638     }
 639
 640     // check this object against the Rsn cuts (if required)
 641     if (fUseRsnTrackCuts) {
 642       if (!fRsnTrackCuts.IsSelected(AliRsnCut::kParticle, &temp)) continue;
 643     }
 644
 645     AliRsnDaughter *ptr = rsn->AddTrack(temp);
 646     if (!ptr) AliWarning(Form("Failed storing track #%d", index));
 647   }
 648
 649   // compute total multiplicity
 650   rsn->MakeComputations();
 651   if (rsn->GetMultiplicity() <= 0)
 652   {
 653     AliDebug(1, "Zero multiplicity in this event");
 654     return kFALSE;
 655   }
 656
 657   // fill the PID index arrays
 658   rsn->FillPIDArrays(fPIDArraysSize);
 659
 660   // correct tracks for primary vertex
 661   rsn->CorrectTracks();
 662
 663   return kTRUE;
 664 }