[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();
    AliDebug(1, 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 TObjArray(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.
    if (fTreeTR) {
	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
		if (ref) {
		    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();
	mcParticle = new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
	fMCParticleMap->AddAt(mcParticle, i);
	if (mcParticle) {
	    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;
	}
	if (event) {
	    return (oldidx - event->GetPrimaryOffset());
	} else {
	    return (-1);
	}
    } else {
	while((event = (AliMCEvent*)next())) {
	    if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
	}
	if (event) {
	    return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
	} else {
	    return (-1);
	}
    }
}

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>
091fa189	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();
477fb828	118
ed97dc98	119	Int_t iev = fHeader->GetEvent();
ed97dc98	120	Int_t ievr = fHeader->GetEventNrInRun();
477fb828	121	AliDebug(1, Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
ed97dc98	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
091fa189	131	fMCParticleMap = new TObjArray(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.
bb278e48	278	if (fTreeTR) {
	279	fTreeTR->SetBranchStatus("*", 0);
	280	fTreeTR->SetBranchStatus("AliRun.*", 1);
	281	fTreeTR->SetBranchStatus("ITS.*", 1);
	282	fTreeTR->SetBranchStatus("TPC.*", 1);
	283	fTreeTR->SetBranchStatus("TRD.*", 1);
	284	fTreeTR->SetBranchStatus("TOF.*", 1);
	285	fTreeTR->SetBranchStatus("FRAME.*", 1);
	286	fTreeTR->SetBranchStatus("MUON.*", 1);
	287	}
	288
415d9f5c	289	//
	290	// Connect the active branches
	291	TClonesArray* trefs[7];
	292	for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
	293	if (fTreeTR){
	294	// make branch for central track references
	295	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	296	if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
	297	if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
	298	if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
	299	if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
	300	if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
	301	if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
	302	}
	303
	304	Int_t np = fStack->GetNprimary();
	305	Int_t nt = fTreeTR->GetEntries();
	306
	307	//
	308	// Loop over tracks and find the secondaries with the help of the kine tree
	309	Int_t ifills = 0;
	310	Int_t it = 0;
	311	Int_t itlast = 0;
	312	TParticle* part;
	313
	314	for (Int_t ip = np - 1; ip > -1; ip--) {
	315	part = fStack->Particle(ip);
	316	// printf("Particle %5d %5d %5d %5d %5d %5d \n",
	317	// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
	318	// part->GetLastDaughter(), part->TestBit(kTransportBit));
	319
	320	// Determine range of secondaries produced by this primary during transport
	321	Int_t dau1 = part->GetFirstDaughter();
	322	if (dau1 < np) continue; // This particle has no secondaries produced during transport
	323	Int_t dau2 = -1;
	324	if (dau1 > -1) {
	325	Int_t inext = ip - 1;
	326	while (dau2 < 0) {
	327	if (inext >= 0) {
	328	part = fStack->Particle(inext);
	329	dau2 = part->GetFirstDaughter();
	330	if (dau2 == -1 \|\| dau2 < np) {
	331	dau2 = -1;
	332	} else {
	333	dau2--;
	334	}
	335	} else {
	336	dau2 = fStack->GetNtrack() - 1;
	337	}
	338	inext--;
	339	} // find upper bound
	340	} // dau2 < 0
	341
	342
	343	// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
	344	//
	345	// Loop over reference hits and find secondary label
	346	// First the tricky part: find the entry in treeTR than contains the hits or
	347	// make sure that no hits exist.
	348	//
	349	Bool_t hasHits = kFALSE;
	350	Bool_t isOutside = kFALSE;
	351
	352	it = itlast;
353	while (!hasHits && !isOutside && it < nt) {
354	fTreeTR->GetEntry(it++);
355	for (Int_t ib = 0; ib < 7; ib++) {
356	if (!trefs[ib]) continue;
357	Int_t nh = trefs[ib]->GetEntries();
358	for (Int_t ih = 0; ih < nh; ih++) {
359	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
360	Int_t label = tr->Label();
361	if (label >= dau1 && label <= dau2) {
362	hasHits = kTRUE;
363	itlast = it - 1;
364	break;
365	}
366	if (label > dau2 \|\| label < ip) {
367	isOutside = kTRUE;
368	itlast = it - 1;
369	break;
370	}
371	} // hits
372	if (hasHits \|\| isOutside) break;
373	} // branches
374	} // entries
375
376	if (!hasHits) {
377	// Write empty entries
378	for (Int_t id = dau1; (id <= dau2); id++) {
379	fTmpTreeTR->Fill();
380	ifills++;
381	}
382	} else {
383	// Collect all hits
384	fTreeTR->GetEntry(itlast);
385	for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
386	for (Int_t ib = 0; ib < 7; ib++) {
387	if (!trefs[ib]) continue;
388	Int_t nh = trefs[ib]->GetEntries();
389	for (Int_t ih = 0; ih < nh; ih++) {
390	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
391	Int_t label = tr->Label();
392	// Skip primaries
393	if (label == ip) continue;
394	if (label > dau2 \|\| label < dau1)
395	printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
396	itlast, label, dau1, dau2);
397	if (label == id) {
398	// secondary found
399	tr->SetDetectorId(ib-1);
93df0e9b	400	Int_t nref = fTRBuffer->GetEntriesFast();
93df0e9b	401	TClonesArray &lref = *fTRBuffer;
415d9f5c	402	new(lref[nref]) AliTrackReference(*tr);
	403	}
	404	} // hits
	405	} // branches
	406	fTmpTreeTR->Fill();
ee2774b7	407	fTRBuffer->Delete();
415d9f5c	408	ifills++;
	409	} // daughters
	410	} // has hits
	411	} // tracks
	412
	413	//
	414	// Now loop again and write the primaries
	415	//
	416	it = nt - 1;
	417	for (Int_t ip = 0; ip < np; ip++) {
	418	Int_t labmax = -1;
	419	while (labmax < ip && it > -1) {
	420	fTreeTR->GetEntry(it--);
	421	for (Int_t ib = 0; ib < 7; ib++) {
	422	if (!trefs[ib]) continue;
	423	Int_t nh = trefs[ib]->GetEntries();
	424	//
	425	// Loop over reference hits and find primary labels
	426	for (Int_t ih = 0; ih < nh; ih++) {
	427	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	428	Int_t label = tr->Label();
	429	if (label < np && label > labmax) {
	430	labmax = label;
	431	}
	432
	433	if (label == ip) {
	434	tr->SetDetectorId(ib-1);
93df0e9b	435	Int_t nref = fTRBuffer->GetEntriesFast();
93df0e9b	436	TClonesArray &lref = *fTRBuffer;
415d9f5c	437	new(lref[nref]) AliTrackReference(*tr);
	438	}
	439	} // hits
	440	} // branches
	441	} // entries
	442	it++;
	443	fTmpTreeTR->Fill();
ee2774b7	444	fTRBuffer->Delete();
415d9f5c	445	ifills++;
	446	} // tracks
	447	// Check
	448
	449
	450	// Clean-up
	451	delete fTreeTR; fTreeTR = 0;
	452
	453	for (Int_t ib = 0; ib < 7; ib++) {
	454	if (trefs[ib]) {
	455	trefs[ib]->Clear();
	456	delete trefs[ib];
	457	trefs[ib] = 0;
	458	}
	459	}
	460
	461	if (ifills != fStack->GetNtrack())
	462	printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
	463	ifills, fStack->GetNtrack());
	464
	465	fTmpTreeTR->Write();
	466	fTreeTR = fTmpTreeTR;
	467	}
	468
7aad0c47	469	AliVParticle* AliMCEvent::GetTrack(Int_t i) const
415d9f5c	470	{
415d9f5c	471	// Get MC Particle i
93836e1b	472	//
7aad0c47	473
	474	if (fExternal) {
	475	return ((AliVParticle*) (fMCParticles->At(i)));
	476	}
	477
93836e1b	478	//
93836e1b	479	// Check first if this explicitely accesses the subsidiary event
7aad0c47	480
93836e1b	481	if (i >= BgLabelOffset()) {
	482	if (fSubsidiaryEvents) {
	483	AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	484	return (bgEvent->GetTrack(i - BgLabelOffset()));
	485	} else {
	486	return 0;
	487	}
	488	}
	489
	490	//
93df0e9b	491	AliMCParticle *mcParticle = 0;
	492	TParticle *particle = 0;
	493	TClonesArray *trefs = 0;
	494	Int_t ntref = 0;
	495	TRefArray *rarray = 0;
93836e1b	496
	497
	498
415d9f5c	499	// Out of range check
	500	if (i < 0 \|\| i >= fNparticles) {
	501	AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
	502	mcParticle = 0;
	503	return (mcParticle);
	504	}
93836e1b	505
415d9f5c	506
93836e1b	507	if (fSubsidiaryEvents) {
	508	AliMCEvent* mc;
	509	Int_t idx = FindIndexAndEvent(i, mc);
	510	return (mc->GetTrack(idx));
	511	}
415d9f5c	512
415d9f5c	513	//
93836e1b	514	// First check If the MC Particle has been already cached
415d9f5c	515	if(!fMCParticleMap->At(i)) {
93df0e9b	516	// Get particle from the stack
	517	particle = fStack->Particle(i);
	518	// Get track references from Tree TR
	519	if (fTreeTR) {
	520	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	521	trefs = fTRBuffer;
	522	ntref = trefs->GetEntriesFast();
	523	rarray = new TRefArray(ntref);
	524	Int_t nen = fTrackReferences->GetEntriesFast();
	525	for (Int_t j = 0; j < ntref; j++) {
	526	// Save the track references in a TClonesArray
	527	AliTrackReference* ref = dynamic_cast<AliTrackReference>((fTRBuffer)[j]);
	528	// Save the pointer in a TRefArray
58c338d2	529	if (ref) {
	530	new ((fTrackReferences)[nen]) AliTrackReference(ref);
	531	rarray->AddAt((*fTrackReferences)[nen], j);
	532	nen++;
	533	}
93df0e9b	534	} // loop over track references for entry i
93df0e9b	535	} // if TreeTR available
415d9f5c	536	Int_t nentries = fMCParticles->GetEntriesFast();
091fa189	537	mcParticle = new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
415d9f5c	538	fMCParticleMap->AddAt(mcParticle, i);
58c338d2	539	if (mcParticle) {
	540	TParticle* part = mcParticle->Particle();
	541	Int_t imo = part->GetFirstMother();
	542	Int_t id1 = part->GetFirstDaughter();
	543	Int_t id2 = part->GetLastDaughter();
	544	if (fPrimaryOffset > 0 \|\| fSecondaryOffset > 0) {
	545	// Remapping of the mother and daughter indices
	546	if (imo < fNprimaries) {
	547	mcParticle->SetMother(imo + fPrimaryOffset);
	548	} else {
	549	mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
	550	}
	551
	552	if (id1 < fNprimaries) {
	553	mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
	554	mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
	555	} else {
	556	mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
	557	mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
	558	}
	559
	560
	561	if (i > fNprimaries) {
	562	mcParticle->SetLabel(i + fPrimaryOffset);
	563	} else {
	564	mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
	565	}
93836e1b	566	} else {
58c338d2	567	mcParticle->SetFirstDaughter(id1);
	568	mcParticle->SetLastDaughter (id2);
	569	mcParticle->SetMother (imo);
93836e1b	570	}
93836e1b	571	}
415d9f5c	572	} else {
	573	mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
	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	}
58c338d2	602	if (event) {
	603	return (oldidx - event->GetPrimaryOffset());
	604	} else {
	605	return (-1);
	606	}
93836e1b	607	} else {
	608	while((event = (AliMCEvent*)next())) {
	609	if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
	610	}
58c338d2	611	if (event) {
	612	return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
	613	} else {
	614	return (-1);
	615	}
93836e1b	616	}
	617	}
	618
	619	Int_t AliMCEvent::BgLabelToIndex(Int_t label)
	620	{
	621	// Convert a background label to an absolute index
	622	if (fSubsidiaryEvents) {
	623	AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	624	label -= BgLabelOffset();
	625	if (label < bgEvent->GetNumberOfPrimaries()) {
	626	label += bgEvent->GetPrimaryOffset();
	627	} else {
	628	label += (bgEvent->GetSecondaryOffset() - fNprimaries);
	629	}
	630	}
	631	return (label);
	632	}
	633
	634
	635	Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i)
	636	{
	637	//
	638	// Delegate to subevent if necesarry
	639
	640
	641	if (!fSubsidiaryEvents) {
	642	return fStack->IsPhysicalPrimary(i);
	643	} else {
	644	AliMCEvent* evt = 0;
	645	Int_t idx = FindIndexAndEvent(i, evt);
	646	return (evt->IsPhysicalPrimary(idx));
	647	}
	648	}
	649
	650	void AliMCEvent::InitEvent()
	651	{
35152e4b	652	//
35152e4b	653	// Initialize the subsidiary event structure
93836e1b	654	if (fSubsidiaryEvents) {
	655	TIter next(fSubsidiaryEvents);
	656	AliMCEvent* evt;
	657	fNprimaries = 0;
	658	fNparticles = 0;
	659
	660	while((evt = (AliMCEvent*)next())) {
	661	fNprimaries += evt->GetNumberOfPrimaries();
	662	fNparticles += evt->GetNumberOfTracks();
	663	}
	664
	665	Int_t ioffp = 0;
	666	Int_t ioffs = fNprimaries;
	667	next.Reset();
	668
	669	while((evt = (AliMCEvent*)next())) {
	670	evt->SetPrimaryOffset(ioffp);
	671	evt->SetSecondaryOffset(ioffs);
	672	ioffp += evt->GetNumberOfPrimaries();
	673	ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
	674	}
	675	}
	676	}
091fa189	677
35152e4b	678	void AliMCEvent::PreReadAll()
	679	{
	680	// Preread the MC information
	681	Int_t i;
	682	// secondaries
	683	for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
	684	{
	685	GetTrack(i);
	686	}
	687	// primaries
	688	for (i = 0; i < fStack->GetNprimary(); i++)
	689	{
	690	GetTrack(i);
	691	}
35152e4b	692	}
93836e1b	693
5df3496b	694	const AliVVertex * AliMCEvent::GetPrimaryVertex() const
	695	{
	696	// Create a MCVertex object from the MCHeader information
570308e4	697	TArrayF v;
570308e4	698	GenEventHeader()->PrimaryVertex(v) ;
5df3496b	699	if (!fVertex) {
5df3496b	700	fVertex = new AliMCVertex(v[0], v[1], v[2]);
570308e4	701	} else {
570308e4	702	((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
5df3496b	703	}
5df3496b	704	return fVertex;
	705	}
	706
	707
415d9f5c	708	ClassImp(AliMCEvent)