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