[u/mrichter/AliRoot.git] / STEER / AliMCEvent.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//-------------------------------------------------------------------------
//     Class for Kinematic Events
//     Author: Andreas Morsch, CERN
//-------------------------------------------------------------------------
#include <TArrow.h>
#include <TMarker.h>
#include <TH2F.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TRefArray.h>
#include <TList.h>
#include <TArrayF.h>

#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCVertex.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"


Int_t AliMCEvent::fgkBgLabelOffset(10000000);


AliMCEvent::AliMCEvent():
    AliVEvent(),
    fStack(0),
    fMCParticles(0),
    fMCParticleMap(0),
    fHeader(new AliHeader()),
    fTRBuffer(0),
    fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
    fTreeTR(0),
    fTmpTreeTR(0),
    fTmpFileTR(0),
    fNprimaries(-1),
    fNparticles(-1),
    fSubsidiaryEvents(0),
    fPrimaryOffset(0),
    fSecondaryOffset(0),
    fExternal(0),
    fVertex(0)
{
    // Default constructor
}

AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
    AliVEvent(mcEvnt),
    fStack(mcEvnt.fStack),
    fMCParticles(mcEvnt.fMCParticles),
    fMCParticleMap(mcEvnt.fMCParticleMap),
    fHeader(mcEvnt.fHeader),
    fTRBuffer(mcEvnt.fTRBuffer),
    fTrackReferences(mcEvnt.fTrackReferences),
    fTreeTR(mcEvnt.fTreeTR),
    fTmpTreeTR(mcEvnt.fTmpTreeTR),
    fTmpFileTR(mcEvnt.fTmpFileTR),
    fNprimaries(mcEvnt.fNprimaries),
    fNparticles(mcEvnt.fNparticles),
    fSubsidiaryEvents(0),
    fPrimaryOffset(0),
    fSecondaryOffset(0),
    fExternal(0),
    fVertex(mcEvnt.fVertex)
{ 
// Copy constructor
}


AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
{
    // assignment operator
    if (this!=&mcEvnt) { 
	AliVEvent::operator=(mcEvnt); 
    }
  
    return *this; 
}

void AliMCEvent::ConnectTreeE (TTree* tree)
{
    // Connect the event header tree
    tree->SetBranchAddress("Header", &fHeader);
}

void AliMCEvent::ConnectTreeK (TTree* tree)
{
    // Connect the kinematics tree to the stack
    if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
    //
    fStack = fHeader->Stack();
    fStack->ConnectTree(tree);
    //
    // Load the event
    fStack->GetEvent();
    fNparticles = fStack->GetNtrack();
    fNprimaries = fStack->GetNprimary();
    Int_t iev  = fHeader->GetEvent();
    Int_t ievr = fHeader->GetEventNrInRun();
    
    AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n", 
		 iev, ievr, fNparticles, fNprimaries));
 
    // This is a cache for the TParticles converted to MCParticles on user request
    if (fMCParticleMap) {
	fMCParticleMap->Clear();
	fMCParticles->Delete();
	if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
    }
    else
	fMCParticleMap = new TRefArray(fNparticles);
}

void AliMCEvent::ConnectTreeTR (TTree* tree)
{
    // Connect the track reference tree
    fTreeTR = tree;
    
    if (fTreeTR->GetBranch("AliRun")) {
	if (fTmpFileTR) {
	    fTmpFileTR->Close();
	    delete fTmpFileTR;
	}
	// This is an old format with one branch per detector not in synch with TreeK
	ReorderAndExpandTreeTR();
    } else {
	// New format 
	fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
    }
}

Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
{
    // Retrieve entry i
    if (i < 0 || i >= fNparticles) {
	AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
	particle = 0;
	trefs    = 0;
	return (-1);
    }
    particle = fStack->Particle(i);
    if (fTreeTR) {
	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	trefs    = fTRBuffer;
	return trefs->GetEntries();
    } else {
	trefs = 0;
	return -1;
    }
}


void AliMCEvent::Clean()
{
    // Clean-up before new trees are connected
    delete fStack; fStack = 0;

    // Clear TR
    if (fTRBuffer) {
	fTRBuffer->Delete();
	delete fTRBuffer;
	fTRBuffer = 0;
    }
}

#include <iostream>

void AliMCEvent::FinishEvent()
{
  // Clean-up after event
  //    
    if (fStack) fStack->Reset(0);
    fMCParticles->Delete();
    
    if (fMCParticleMap) 
      fMCParticleMap->Clear();
    if (fTRBuffer) {
      fTRBuffer->Delete();
    }
    //    fTrackReferences->Delete();
    fTrackReferences->Clear();
    fNparticles = -1;
    fNprimaries = -1;    
    fStack      =  0;
//    fSubsidiaryEvents->Clear();
    fSubsidiaryEvents = 0;
}


void AliMCEvent::DrawCheck(Int_t i, Int_t search)
{
    //
    // Simple event display for debugging
    if (!fTreeTR) {
	AliWarning("No Track Reference information available");
	return;
    } 
    
    if (i > -1 && i < fNparticles) {
	fTreeTR->GetEntry(fStack->TreeKEntry(i));
    } else {
	AliWarning("AliMCEvent::GetEntry: Index out of range");
    }
    
    Int_t nh = fTRBuffer->GetEntries();
    
    
    if (search) {
	while(nh <= search && i < fNparticles - 1) {
	    i++;
	    fTreeTR->GetEntry(fStack->TreeKEntry(i));
	    nh =  fTRBuffer->GetEntries();
	}
	printf("Found Hits at %5d\n", i);
    }
    TParticle* particle = fStack->Particle(i);
    
    TH2F*    h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
    Float_t x0 = particle->Vx();
    Float_t y0 = particle->Vy();

    Float_t x1 = particle->Vx() + particle->Px() * 50.;
    Float_t y1 = particle->Vy() + particle->Py() * 50.;
    
    TArrow*  a = new TArrow(x0, y0, x1, y1, 0.01);
    h->Draw();
    a->SetLineColor(2);
    
    a->Draw();
    
    for (Int_t ih = 0; ih < nh; ih++) {
	AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	m->Draw();
	m->SetMarkerSize(0.4);
	
    }
}


void AliMCEvent::ReorderAndExpandTreeTR()
{
//
//  Reorder and expand the track reference tree in order to match the kinematics tree.
//  Copy the information from different branches into one
//
//  TreeTR

    fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
    fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
    if (!fTRBuffer)  fTRBuffer = new TClonesArray("AliTrackReference", 100);
    fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 0);
    

//
//  Activate the used branches only. Otherwisw we get a bad memory leak.
    fTreeTR->SetBranchStatus("*",        0);
    fTreeTR->SetBranchStatus("AliRun.*", 1);
    fTreeTR->SetBranchStatus("ITS.*",    1);
    fTreeTR->SetBranchStatus("TPC.*",    1);
    fTreeTR->SetBranchStatus("TRD.*",    1);
    fTreeTR->SetBranchStatus("TOF.*",    1);
    fTreeTR->SetBranchStatus("FRAME.*",  1);
    fTreeTR->SetBranchStatus("MUON.*",   1);
//
//  Connect the active branches
    TClonesArray* trefs[7];
    for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
    if (fTreeTR){
	// make branch for central track references
	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
	if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
	if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
	if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
	if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
	if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
    }

    Int_t np = fStack->GetNprimary();
    Int_t nt = fTreeTR->GetEntries();
    
    //
    // Loop over tracks and find the secondaries with the help of the kine tree
    Int_t ifills = 0;
    Int_t it     = 0;
    Int_t itlast = 0;
    TParticle* part;

    for (Int_t ip = np - 1; ip > -1; ip--) {
	part = fStack->Particle(ip);
//	printf("Particle %5d %5d %5d %5d %5d %5d \n", 
//	       ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), 
//	       part->GetLastDaughter(), part->TestBit(kTransportBit));

	// Determine range of secondaries produced by this primary during transport	
	Int_t dau1  = part->GetFirstDaughter();
	if (dau1 < np) continue;  // This particle has no secondaries produced during transport
	Int_t dau2  = -1;
	if (dau1 > -1) {
	    Int_t inext = ip - 1;
	    while (dau2 < 0) {
		if (inext >= 0) {
		    part = fStack->Particle(inext);
		    dau2 =  part->GetFirstDaughter();
		    if (dau2 == -1 || dau2 < np) {
			dau2 = -1;
		    } else {
			dau2--;
		    }
		} else {
		    dau2 = fStack->GetNtrack() - 1;
		}
		inext--;
	    } // find upper bound
	}  // dau2 < 0
	

//	printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
//
// Loop over reference hits and find secondary label
// First the tricky part: find the entry in treeTR than contains the hits or
// make sure that no hits exist.
//
	Bool_t hasHits   = kFALSE;
	Bool_t isOutside = kFALSE;

	it = itlast;
	while (!hasHits && !isOutside && it < nt) {
	    fTreeTR->GetEntry(it++);
	    for (Int_t ib = 0; ib < 7; ib++) {
		if (!trefs[ib]) continue;
		Int_t nh = trefs[ib]->GetEntries();
		for (Int_t ih = 0; ih < nh; ih++) {
		    AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
		    Int_t label = tr->Label();
		    if (label >= dau1 && label <= dau2) {
			hasHits = kTRUE;
			itlast = it - 1;
			break;
		    }
		    if (label > dau2 || label < ip) {
			isOutside = kTRUE;
			itlast = it - 1;
			break;
		    }
		} // hits
		if (hasHits || isOutside) break;
	    } // branches
	} // entries

	if (!hasHits) {
	    // Write empty entries
	    for (Int_t id = dau1; (id <= dau2); id++) {
		fTmpTreeTR->Fill();
		ifills++;
	    } 
	} else {
	    // Collect all hits
	    fTreeTR->GetEntry(itlast);
	    for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
		for (Int_t ib = 0; ib < 7; ib++) {
		    if (!trefs[ib]) continue;
		    Int_t nh = trefs[ib]->GetEntries();
		    for (Int_t ih = 0; ih < nh; ih++) {
			AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
			Int_t label = tr->Label();
			// Skip primaries
			if (label == ip) continue;
			if (label > dau2 || label < dau1) 
			    printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n", 
				   itlast, label, dau1, dau2);
			if (label == id) {
			    // secondary found
			    tr->SetDetectorId(ib-1);
			    Int_t nref =  fTRBuffer->GetEntriesFast();
			    TClonesArray &lref = *fTRBuffer;
			    new(lref[nref]) AliTrackReference(*tr);
			}
		    } // hits
		} // branches
		fTmpTreeTR->Fill();
		fTRBuffer->Delete();
		ifills++;
	    } // daughters
	} // has hits
    } // tracks

    //
    // Now loop again and write the primaries
    //
    it = nt - 1;
    for (Int_t ip = 0; ip < np; ip++) {
	Int_t labmax = -1;
	while (labmax < ip && it > -1) {
	    fTreeTR->GetEntry(it--);
	    for (Int_t ib = 0; ib < 7; ib++) {
		if (!trefs[ib]) continue;
		Int_t nh = trefs[ib]->GetEntries();
		// 
		// Loop over reference hits and find primary labels
		for (Int_t ih = 0; ih < nh; ih++) {
		    AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
		    Int_t label = tr->Label();
		    if (label < np && label > labmax) {
			labmax = label;
		    }
		    
		    if (label == ip) {
			tr->SetDetectorId(ib-1);
			Int_t nref = fTRBuffer->GetEntriesFast();
			TClonesArray &lref = *fTRBuffer;
			new(lref[nref]) AliTrackReference(*tr);
		    }
		} // hits
	    } // branches
	} // entries
	it++;
	fTmpTreeTR->Fill();
	fTRBuffer->Delete();
	ifills++;
    } // tracks
    // Check


    // Clean-up
    delete fTreeTR; fTreeTR = 0;
    
    for (Int_t ib = 0; ib < 7; ib++) {
	if (trefs[ib]) {
	    trefs[ib]->Clear();
	    delete trefs[ib];
	    trefs[ib] = 0;
	}
    }

    if (ifills != fStack->GetNtrack()) 
	printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", 
	       ifills, fStack->GetNtrack());

    fTmpTreeTR->Write();
    fTreeTR = fTmpTreeTR;
}

AliVParticle* AliMCEvent::GetTrack(Int_t i) const
{
    // Get MC Particle i
    //

    if (fExternal) {
	return ((AliVParticle*) (fMCParticles->At(i)));
    }
    
    //
    // Check first if this explicitely accesses the subsidiary event
    
    if (i >= BgLabelOffset()) {
	if (fSubsidiaryEvents) {
	    AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	    return (bgEvent->GetTrack(i - BgLabelOffset()));
	} else {
	    return 0;
	}
    }
    
    //
    AliMCParticle *mcParticle = 0;
    TParticle     *particle   = 0;
    TClonesArray  *trefs      = 0;
    Int_t          ntref      = 0;
    TRefArray     *rarray     = 0;


    // Out of range check
    if (i < 0 || i >= fNparticles) {
	AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
	mcParticle = 0;
	return (mcParticle);
    }

    
    if (fSubsidiaryEvents) {
	AliMCEvent*   mc;
	Int_t idx = FindIndexAndEvent(i, mc);
	return (mc->GetTrack(idx));
    } 

    //
    // First check If the MC Particle has been already cached
    if(!fMCParticleMap->At(i)) {
	// Get particle from the stack
	particle   = fStack->Particle(i);
	// Get track references from Tree TR
	if (fTreeTR) {
	    fTreeTR->GetEntry(fStack->TreeKEntry(i));
	    trefs     = fTRBuffer;
	    ntref     = trefs->GetEntriesFast();
	    rarray    = new TRefArray(ntref);
	    Int_t nen = fTrackReferences->GetEntriesFast();
	    for (Int_t j = 0; j < ntref; j++) {
		// Save the track references in a TClonesArray
		AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
		// Save the pointer in a TRefArray
		new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
		rarray->AddAt((*fTrackReferences)[nen], j);
		nen++;
	    } // loop over track references for entry i
	} // if TreeTR available
	Int_t nentries = fMCParticles->GetEntriesFast();
	new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
	mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
	fMCParticleMap->AddAt(mcParticle, i);

	TParticle* part = mcParticle->Particle();
	Int_t imo  = part->GetFirstMother();
	Int_t id1  = part->GetFirstDaughter();
	Int_t id2  = part->GetLastDaughter();
	if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
	    // Remapping of the mother and daughter indices
	    if (imo < fNprimaries) {
		mcParticle->SetMother(imo + fPrimaryOffset);
	    } else {
		mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
	    }

	    if (id1 < fNprimaries) {
		mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
		mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
	    } else {
		mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
		mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
	    }

	    
	    if (i > fNprimaries) {
		mcParticle->SetLabel(i + fPrimaryOffset);
	    } else {
		mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
	    }
	    
	    
	} else {
	    mcParticle->SetFirstDaughter(id1);
	    mcParticle->SetLastDaughter (id2);
	    mcParticle->SetMother       (imo);
	}

    } else {
	mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
    }
    
    
    return mcParticle;
}

AliGenEventHeader* AliMCEvent::GenEventHeader() const {return (fHeader->GenEventHeader());}


void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event) 
{
    // Add a subsidiary event to the list; for example merged background event.
    if (!fSubsidiaryEvents) {
	TList* events = new TList();
	events->Add(new AliMCEvent(*this));
	fSubsidiaryEvents = events;
    }
    
    fSubsidiaryEvents->Add(event);
}

Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
{
    // Find the index and event in case of composed events like signal + background
    TIter next(fSubsidiaryEvents);
    next.Reset();
     if (oldidx < fNprimaries) {
	while((event = (AliMCEvent*)next())) {
	    if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
	}
	return (oldidx - event->GetPrimaryOffset());
    } else {
	while((event = (AliMCEvent*)next())) {
	    if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
	}
	return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
    }
}

Int_t AliMCEvent::BgLabelToIndex(Int_t label)
{
    // Convert a background label to an absolute index
    if (fSubsidiaryEvents) {
	AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	label -= BgLabelOffset();
	if (label < bgEvent->GetNumberOfPrimaries()) {
	    label += bgEvent->GetPrimaryOffset();
	} else {
	    label += (bgEvent->GetSecondaryOffset() - fNprimaries);
	}
    }
    return (label);
}


Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) 
{
//
// Delegate to subevent if necesarry 

    
    if (!fSubsidiaryEvents) {
	return fStack->IsPhysicalPrimary(i);
    } else {
	AliMCEvent* evt = 0;
	Int_t idx = FindIndexAndEvent(i, evt);
	return (evt->IsPhysicalPrimary(idx));
    }
}

void AliMCEvent::InitEvent()
{
//
// Initialize the subsidiary event structure
    if (fSubsidiaryEvents) {
	TIter next(fSubsidiaryEvents);
	AliMCEvent* evt;
	fNprimaries = 0;
	fNparticles = 0;
	
	while((evt = (AliMCEvent*)next())) {
	    fNprimaries += evt->GetNumberOfPrimaries();	
	    fNparticles += evt->GetNumberOfTracks();    
	}
	
	Int_t ioffp = 0;
	Int_t ioffs = fNprimaries;
	next.Reset();
	
	while((evt = (AliMCEvent*)next())) {
	    evt->SetPrimaryOffset(ioffp);
	    evt->SetSecondaryOffset(ioffs);
	    ioffp += evt->GetNumberOfPrimaries();
	    ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());	    
	}
    }
}
void AliMCEvent::PreReadAll()                              
{
    // Preread the MC information
    Int_t i;
    // secondaries
    for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++) 
    {
	GetTrack(i);
    }
    // primaries
    for (i = 0; i < fStack->GetNprimary(); i++) 
    {
	GetTrack(i);
    }
    
    
}


const AliVVertex * AliMCEvent::GetPrimaryVertex() const 
{
    // Create a MCVertex object from the MCHeader information
    TArrayF v;
    GenEventHeader()->PrimaryVertex(v) ;
    if (!fVertex) {
	fVertex = new AliMCVertex(v[0], v[1], v[2]);
    } else {
	((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
    }
    return fVertex;
}


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