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