[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;
    }
}

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