[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;
    }
    if (fVertex) {
	delete fVertex;
	fVertex = 0;
    }
    
}

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