STEER/AliMCEvent.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //     Class for Kinematic Events
  20 //     Author: Andreas Morsch, CERN
  21 //-------------------------------------------------------------------------
  22 #include <TArrow.h>
  23 #include <TMarker.h>
  24 #include <TH2F.h>
  25 #include <TTree.h>
  26 #include <TFile.h>
  27 #include <TParticle.h>
  28 #include <TClonesArray.h>
  29 #include <TRefArray.h>
  30 #include <TList.h>
  31
  32 #include "AliLog.h"
  33 #include "AliMCEvent.h"
  34 #include "AliStack.h"
  35 #include "AliTrackReference.h"
  36 #include "AliHeader.h"
  37 #include "AliGenEventHeader.h"
  38
  39
  40 Int_t AliMCEvent::fgkBgLabelOffset(10000000);
  41
  42
  43 AliMCEvent::AliMCEvent():
  44     AliVEvent(),
  45     fStack(0),
  46     fMCParticles(0),
  47     fMCParticleMap(0),
  48     fHeader(new AliHeader()),
  49     fTRBuffer(0),
  50     fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
  51     fTreeTR(0),
  52     fTmpTreeTR(0),
  53     fTmpFileTR(0),
  54     fNprimaries(-1),
  55     fNparticles(-1),
  56     fSubsidiaryEvents(0),
  57     fPrimaryOffset(0),
  58     fSecondaryOffset(0),
  59     fExternal(0)
  60 {
  61     // Default constructor
  62 }
  63
  64 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
  65     AliVEvent(mcEvnt),
  66     fStack(mcEvnt.fStack),
  67     fMCParticles(mcEvnt.fMCParticles),
  68     fMCParticleMap(mcEvnt.fMCParticleMap),
  69     fHeader(mcEvnt.fHeader),
  70     fTRBuffer(mcEvnt.fTRBuffer),
  71     fTrackReferences(mcEvnt.fTrackReferences),
  72     fTreeTR(mcEvnt.fTreeTR),
  73     fTmpTreeTR(mcEvnt.fTmpTreeTR),
  74     fTmpFileTR(mcEvnt.fTmpFileTR),
  75     fNprimaries(mcEvnt.fNprimaries),
  76     fNparticles(mcEvnt.fNparticles),
  77     fSubsidiaryEvents(0),
  78     fPrimaryOffset(0),
  79     fSecondaryOffset(0),
  80     fExternal(0)
  81 {
  82 // Copy constructor
  83 }
  84
  85
  86 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
  87 {
  88     // assignment operator
  89     if (this!=&mcEvnt) {
  90         AliVEvent::operator=(mcEvnt);
  91     }
  92
  93     return *this;
  94 }
  95
  96 void AliMCEvent::ConnectTreeE (TTree* tree)
  97 {
  98     // Connect the event header tree
  99     tree->SetBranchAddress("Header", &fHeader);
 100 }
 101
 102 void AliMCEvent::ConnectTreeK (TTree* tree)
 103 {
 104     // Connect the kinematics tree to the stack
 105     if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
 106     //
 107     fStack = fHeader->Stack();
 108     fStack->ConnectTree(tree);
 109     //
 110     // Load the event
 111     fStack->GetEvent();
 112     fNparticles = fStack->GetNtrack();
 113     fNprimaries = fStack->GetNprimary();
 114     Int_t iev  = fHeader->GetEvent();
 115     Int_t ievr = fHeader->GetEventNrInRun();
 116
 117     AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
 118                  iev, ievr, fNparticles, fNprimaries));
 119
 120     // This is a cache for the TParticles converted to MCParticles on user request
 121     if (fMCParticleMap) {
 122         fMCParticleMap->Clear();
 123         fMCParticles->Delete();
 124         if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
 125     }
 126     else
 127         fMCParticleMap = new TRefArray(fNparticles);
 128 }
 129
 130 void AliMCEvent::ConnectTreeTR (TTree* tree)
 131 {
 132     // Connect the track reference tree
 133     fTreeTR = tree;
 134
 135     if (fTreeTR->GetBranch("AliRun")) {
 136         if (fTmpFileTR) {
 137             fTmpFileTR->Close();
 138             delete fTmpFileTR;
 139         }
 140         // This is an old format with one branch per detector not in synch with TreeK
 141         ReorderAndExpandTreeTR();
 142     } else {
 143         // New format
 144         fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
 145     }
 146 }
 147
 148 Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 149 {
 150     // Retrieve entry i
 151     if (i < 0 || i >= fNparticles) {
 152         AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
 153         particle = 0;
 154         trefs    = 0;
 155         return (-1);
 156     }
 157     particle = fStack->Particle(i);
 158     if (fTreeTR) {
 159         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 160         trefs    = fTRBuffer;
 161         return trefs->GetEntries();
 162     } else {
 163         trefs = 0;
 164         return -1;
 165     }
 166 }
 167
 168
 169 void AliMCEvent::Clean()
 170 {
 171     // Clean-up before new trees are connected
 172 //    if (fHeader) {
 173 //      delete fHeader;
 174 //      fHeader = 0;
 175 //    }
 176
 177     delete fStack; fStack = 0;
 178
 179     // Clear TR
 180     if (fTRBuffer) {
 181         fTRBuffer->Delete();
 182         delete fTRBuffer;
 183         fTRBuffer = 0;
 184     }
 185 }
 186
 187 void AliMCEvent::FinishEvent()
 188 {
 189   // Clean-up after event
 190   //
 191     fStack->Reset(0);
 192     fMCParticles->Delete();
 193     fMCParticleMap->Clear();
 194     if (fTRBuffer)
 195       fTRBuffer->Delete();
 196     fTrackReferences->Delete();
 197     fNparticles = -1;
 198     fNprimaries = -1;
 199     fStack      =  0;
 200 //    fSubsidiaryEvents->Clear();
 201     fSubsidiaryEvents = 0;
 202 }
 203
 204
 205
 206 void AliMCEvent::DrawCheck(Int_t i, Int_t search)
 207 {
 208     //
 209     // Simple event display for debugging
 210     if (!fTreeTR) {
 211         AliWarning("No Track Reference information available");
 212         return;
 213     }
 214
 215     if (i > -1 && i < fNparticles) {
 216         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 217     } else {
 218         AliWarning("AliMCEvent::GetEntry: Index out of range");
 219     }
 220
 221     Int_t nh = fTRBuffer->GetEntries();
 222
 223
 224     if (search) {
 225         while(nh <= search && i < fNparticles - 1) {
 226             i++;
 227             fTreeTR->GetEntry(fStack->TreeKEntry(i));
 228             nh =  fTRBuffer->GetEntries();
 229         }
 230         printf("Found Hits at %5d\n", i);
 231     }
 232     TParticle* particle = fStack->Particle(i);
 233
 234     TH2F*    h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
 235     Float_t x0 = particle->Vx();
 236     Float_t y0 = particle->Vy();
 237
 238     Float_t x1 = particle->Vx() + particle->Px() * 50.;
 239     Float_t y1 = particle->Vy() + particle->Py() * 50.;
 240
 241     TArrow*  a = new TArrow(x0, y0, x1, y1, 0.01);
 242     h->Draw();
 243     a->SetLineColor(2);
 244
 245     a->Draw();
 246
 247     for (Int_t ih = 0; ih < nh; ih++) {
 248         AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
 249         TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
 250         m->Draw();
 251         m->SetMarkerSize(0.4);
 252
 253     }
 254 }
 255
 256
 257 void AliMCEvent::ReorderAndExpandTreeTR()
 258 {
 259 //
 260 //  Reorder and expand the track reference tree in order to match the kinematics tree.
 261 //  Copy the information from different branches into one
 262 //
 263 //  TreeTR
 264
 265     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
 266     fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
 267     if (!fTRBuffer)  fTRBuffer = new TClonesArray("AliTrackReference", 100);
 268     fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
 269
 270
 271 //
 272 //  Activate the used branches only. Otherwisw we get a bad memory leak.
 273     fTreeTR->SetBranchStatus("*",        0);
 274     fTreeTR->SetBranchStatus("AliRun.*", 1);
 275     fTreeTR->SetBranchStatus("ITS.*",    1);
 276     fTreeTR->SetBranchStatus("TPC.*",    1);
 277     fTreeTR->SetBranchStatus("TRD.*",    1);
 278     fTreeTR->SetBranchStatus("TOF.*",    1);
 279     fTreeTR->SetBranchStatus("FRAME.*",  1);
 280     fTreeTR->SetBranchStatus("MUON.*",   1);
 281 //
 282 //  Connect the active branches
 283     TClonesArray* trefs[7];
 284     for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
 285     if (fTreeTR){
 286         // make branch for central track references
 287         if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
 288         if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
 289         if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
 290         if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
 291         if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
 292         if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
 293         if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
 294     }
 295
 296     Int_t np = fStack->GetNprimary();
 297     Int_t nt = fTreeTR->GetEntries();
 298
 299     //
 300     // Loop over tracks and find the secondaries with the help of the kine tree
 301     Int_t ifills = 0;
 302     Int_t it     = 0;
 303     Int_t itlast = 0;
 304     TParticle* part;
 305
 306     for (Int_t ip = np - 1; ip > -1; ip--) {
 307         part = fStack->Particle(ip);
 308 //      printf("Particle %5d %5d %5d %5d %5d %5d \n",
 309 //             ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
 310 //             part->GetLastDaughter(), part->TestBit(kTransportBit));
 311
 312         // Determine range of secondaries produced by this primary during transport
 313         Int_t dau1  = part->GetFirstDaughter();
 314         if (dau1 < np) continue;  // This particle has no secondaries produced during transport
 315         Int_t dau2  = -1;
 316         if (dau1 > -1) {
 317             Int_t inext = ip - 1;
 318             while (dau2 < 0) {
 319                 if (inext >= 0) {
 320                     part = fStack->Particle(inext);
 321                     dau2 =  part->GetFirstDaughter();
 322                     if (dau2 == -1 || dau2 < np) {
 323                         dau2 = -1;
 324                     } else {
 325                         dau2--;
 326                     }
 327                 } else {
 328                     dau2 = fStack->GetNtrack() - 1;
 329                 }
 330                 inext--;
 331             } // find upper bound
 332         }  // dau2 < 0
 333
 334
 335 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
 336 //
 337 // Loop over reference hits and find secondary label
 338 // First the tricky part: find the entry in treeTR than contains the hits or
 339 // make sure that no hits exist.
 340 //
 341         Bool_t hasHits   = kFALSE;
 342         Bool_t isOutside = kFALSE;
 343
 344         it = itlast;
 345         while (!hasHits && !isOutside && it < nt) {
 346             fTreeTR->GetEntry(it++);
 347             for (Int_t ib = 0; ib < 7; ib++) {
 348                 if (!trefs[ib]) continue;
 349                 Int_t nh = trefs[ib]->GetEntries();
 350                 for (Int_t ih = 0; ih < nh; ih++) {
 351                     AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 352                     Int_t label = tr->Label();
 353                     if (label >= dau1 && label <= dau2) {
 354                         hasHits = kTRUE;
 355                         itlast = it - 1;
 356                         break;
 357                     }
 358                     if (label > dau2 || label < ip) {
 359                         isOutside = kTRUE;
 360                         itlast = it - 1;
 361                         break;
 362                     }
 363                 } // hits
 364                 if (hasHits || isOutside) break;
 365             } // branches
 366         } // entries
 367
 368         if (!hasHits) {
 369             // Write empty entries
 370             for (Int_t id = dau1; (id <= dau2); id++) {
 371                 fTmpTreeTR->Fill();
 372                 ifills++;
 373             }
 374         } else {
 375             // Collect all hits
 376             fTreeTR->GetEntry(itlast);
 377             for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
 378                 for (Int_t ib = 0; ib < 7; ib++) {
 379                     if (!trefs[ib]) continue;
 380                     Int_t nh = trefs[ib]->GetEntries();
 381                     for (Int_t ih = 0; ih < nh; ih++) {
 382                         AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 383                         Int_t label = tr->Label();
 384                         // Skip primaries
 385                         if (label == ip) continue;
 386                         if (label > dau2 || label < dau1)
 387                             printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
 388                                    itlast, label, dau1, dau2);
 389                         if (label == id) {
 390                             // secondary found
 391                             tr->SetDetectorId(ib-1);
 392                             Int_t nref =  fTRBuffer->GetEntriesFast();
 393                             TClonesArray &lref = *fTRBuffer;
 394                             new(lref[nref]) AliTrackReference(*tr);
 395                         }
 396                     } // hits
 397                 } // branches
 398                 fTmpTreeTR->Fill();
 399                 fTRBuffer->Delete();
 400                 ifills++;
 401             } // daughters
 402         } // has hits
 403     } // tracks
 404
 405     //
 406     // Now loop again and write the primaries
 407     //
 408     it = nt - 1;
 409     for (Int_t ip = 0; ip < np; ip++) {
 410         Int_t labmax = -1;
 411         while (labmax < ip && it > -1) {
 412             fTreeTR->GetEntry(it--);
 413             for (Int_t ib = 0; ib < 7; ib++) {
 414                 if (!trefs[ib]) continue;
 415                 Int_t nh = trefs[ib]->GetEntries();
 416                 //
 417                 // Loop over reference hits and find primary labels
 418                 for (Int_t ih = 0; ih < nh; ih++) {
 419                     AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
 420                     Int_t label = tr->Label();
 421                     if (label < np && label > labmax) {
 422                         labmax = label;
 423                     }
 424
 425                     if (label == ip) {
 426                         tr->SetDetectorId(ib-1);
 427                         Int_t nref = fTRBuffer->GetEntriesFast();
 428                         TClonesArray &lref = *fTRBuffer;
 429                         new(lref[nref]) AliTrackReference(*tr);
 430                     }
 431                 } // hits
 432             } // branches
 433         } // entries
 434         it++;
 435         fTmpTreeTR->Fill();
 436         fTRBuffer->Delete();
 437         ifills++;
 438     } // tracks
 439     // Check
 440
 441
 442     // Clean-up
 443     delete fTreeTR; fTreeTR = 0;
 444
 445     for (Int_t ib = 0; ib < 7; ib++) {
 446         if (trefs[ib]) {
 447             trefs[ib]->Clear();
 448             delete trefs[ib];
 449             trefs[ib] = 0;
 450         }
 451     }
 452
 453     if (ifills != fStack->GetNtrack())
 454         printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
 455                ifills, fStack->GetNtrack());
 456
 457     fTmpTreeTR->Write();
 458     fTreeTR = fTmpTreeTR;
 459 }
 460
 461 AliVParticle* AliMCEvent::GetTrack(Int_t i) const
 462 {
 463     // Get MC Particle i
 464     //
 465
 466     if (fExternal) {
 467         return ((AliVParticle*) (fMCParticles->At(i)));
 468     }
 469
 470     //
 471     // Check first if this explicitely accesses the subsidiary event
 472
 473     if (i >= BgLabelOffset()) {
 474         if (fSubsidiaryEvents) {
 475             AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
 476             return (bgEvent->GetTrack(i - BgLabelOffset()));
 477         } else {
 478             return 0;
 479         }
 480     }
 481
 482     //
 483     AliMCParticle *mcParticle = 0;
 484     TParticle     *particle   = 0;
 485     TClonesArray  *trefs      = 0;
 486     Int_t          ntref      = 0;
 487     TRefArray     *rarray     = 0;
 488
 489
 490
 491     // Out of range check
 492     if (i < 0 || i >= fNparticles) {
 493         AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
 494         mcParticle = 0;
 495         return (mcParticle);
 496     }
 497
 498
 499     if (fSubsidiaryEvents) {
 500         AliMCEvent*   mc;
 501         Int_t idx = FindIndexAndEvent(i, mc);
 502         return (mc->GetTrack(idx));
 503     }
 504
 505     //
 506     // First check If the MC Particle has been already cached
 507     if(!fMCParticleMap->At(i)) {
 508         // Get particle from the stack
 509         particle   = fStack->Particle(i);
 510         // Get track references from Tree TR
 511         if (fTreeTR) {
 512             fTreeTR->GetEntry(fStack->TreeKEntry(i));
 513             trefs     = fTRBuffer;
 514             ntref     = trefs->GetEntriesFast();
 515             rarray    = new TRefArray(ntref);
 516             Int_t nen = fTrackReferences->GetEntriesFast();
 517             for (Int_t j = 0; j < ntref; j++) {
 518                 // Save the track references in a TClonesArray
 519                 AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
 520                 // Save the pointer in a TRefArray
 521                 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
 522                 rarray->AddAt((*fTrackReferences)[nen], j);
 523                 nen++;
 524             } // loop over track references for entry i
 525         } // if TreeTR available
 526         Int_t nentries = fMCParticles->GetEntriesFast();
 527         new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
 528         mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
 529         fMCParticleMap->AddAt(mcParticle, i);
 530
 531         TParticle* part = mcParticle->Particle();
 532         Int_t imo  = part->GetFirstMother();
 533         Int_t id1  = part->GetFirstDaughter();
 534         Int_t id2  = part->GetLastDaughter();
 535         if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
 536             // Remapping of the mother and daughter indices
 537             if (imo < fNprimaries) {
 538                 mcParticle->SetMother(imo + fPrimaryOffset);
 539             } else {
 540                 mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
 541             }
 542
 543             if (id1 < fNprimaries) {
 544                 mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
 545                 mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
 546             } else {
 547                 mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
 548                 mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
 549             }
 550
 551
 552             if (i > fNprimaries) {
 553                 mcParticle->SetLabel(i + fPrimaryOffset);
 554             } else {
 555                 mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
 556             }
 557
 558
 559         } else {
 560             mcParticle->SetFirstDaughter(id1);
 561             mcParticle->SetLastDaughter (id2);
 562             mcParticle->SetMother       (imo);
 563         }
 564
 565     } else {
 566         mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
 567     }
 568
 569
 570     return mcParticle;
 571 }
 572
 573     AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
 574
 575
 576 void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
 577 {
 578     // Add a subsidiary event to the list; for example merged background event.
 579     if (!fSubsidiaryEvents) {
 580         TList* events = new TList();
 581         events->Add(new AliMCEvent(*this));
 582         fSubsidiaryEvents = events;
 583     }
 584
 585     fSubsidiaryEvents->Add(event);
 586 }
 587
 588 Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
 589 {
 590     // Find the index and event in case of composed events like signal + background
 591     TIter next(fSubsidiaryEvents);
 592     next.Reset();
 593      if (oldidx < fNprimaries) {
 594         while((event = (AliMCEvent*)next())) {
 595             if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
 596         }
 597         return (oldidx - event->GetPrimaryOffset());
 598     } else {
 599         while((event = (AliMCEvent*)next())) {
 600             if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
 601         }
 602         return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
 603     }
 604 }
 605
 606 Int_t AliMCEvent::BgLabelToIndex(Int_t label)
 607 {
 608     // Convert a background label to an absolute index
 609     if (fSubsidiaryEvents) {
 610         AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
 611         label -= BgLabelOffset();
 612         if (label < bgEvent->GetNumberOfPrimaries()) {
 613             label += bgEvent->GetPrimaryOffset();
 614         } else {
 615             label += (bgEvent->GetSecondaryOffset() - fNprimaries);
 616         }
 617     }
 618     return (label);
 619 }
 620
 621
 622 Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i)
 623 {
 624 //
 625 // Delegate to subevent if necesarry
 626
 627
 628     if (!fSubsidiaryEvents) {
 629         return fStack->IsPhysicalPrimary(i);
 630     } else {
 631         AliMCEvent* evt = 0;
 632         Int_t idx = FindIndexAndEvent(i, evt);
 633         return (evt->IsPhysicalPrimary(idx));
 634     }
 635 }
 636
 637 void AliMCEvent::InitEvent()
 638 {
 639 //
 640 // Initialize the subsidiary event structure
 641     if (fSubsidiaryEvents) {
 642         TIter next(fSubsidiaryEvents);
 643         AliMCEvent* evt;
 644         fNprimaries = 0;
 645         fNparticles = 0;
 646
 647         while((evt = (AliMCEvent*)next())) {
 648             fNprimaries += evt->GetNumberOfPrimaries();
 649             fNparticles += evt->GetNumberOfTracks();
 650         }
 651
 652         Int_t ioffp = 0;
 653         Int_t ioffs = fNprimaries;
 654         next.Reset();
 655
 656         while((evt = (AliMCEvent*)next())) {
 657             evt->SetPrimaryOffset(ioffp);
 658             evt->SetSecondaryOffset(ioffs);
 659             ioffp += evt->GetNumberOfPrimaries();
 660             ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
 661         }
 662     }
 663 }
 664 void AliMCEvent::PreReadAll()
 665 {
 666     // Preread the MC information
 667     Int_t i;
 668     // secondaries
 669     for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
 670     {
 671         GetTrack(i);
 672     }
 673     // primaries
 674     for (i = 0; i < fStack->GetNprimary(); i++)
 675     {
 676         GetTrack(i);
 677     }
 678
 679
 680 }
 681
 682
 683 ClassImp(AliMCEvent)