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