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