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