[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 "AliLog.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"


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

AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
    AliVEvent(mcEvnt),
    fStack(0),
    fMCParticles(0),
    fMCParticleMap(0),
    fHeader(0),
    fTRBuffer(0),
    fTrackReferences(0),
    fTreeTR(0),
    fTmpTreeTR(0),
    fTmpFileTR(0),
    fNprimaries(-1),
    fNparticles(-1) 
{ 
// 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
    fStack = fHeader->Stack();
    fStack->ConnectTree(tree);
    //
    // Load the event
    fStack->GetEvent();
    fNparticles = fStack->GetNtrack();
    fNprimaries = fStack->GetNprimary();
    AliInfo(Form("AliMCEvent: Number of particles: %5d (all) %5d (primaries)\n", 
		 fNparticles, fNprimaries));
 
    // This is a cache for the TParticles converted to MCParticles on user request
    if (fMCParticleMap) {
	fMCParticleMap->Clear();
	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
    
    if (fHeader) {
	delete fHeader;
	fHeader = 0;
    }
    
    delete fStack;

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

void AliMCEvent::FinishEvent()
{
    // Clean-up after event
     fStack->Reset(0);
     fMCParticles->Clear();
     fTrackReferences->Clear();
}


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->Clear();
		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->Clear();
	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;
}

AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
{
    // Get MC Particle i
    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);
    }
    

    //
    // First check of 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);
	mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
	fMCParticleMap->AddAt(mcParticle, i);
    } else {
	mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
    }

    return mcParticle;
}

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

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>
415d9f5c	30
	31	#include "AliLog.h"
	32	#include "AliMCEvent.h"
	33	#include "AliStack.h"
	34	#include "AliTrackReference.h"
	35	#include "AliHeader.h"
83626b6b	36	#include "AliGenEventHeader.h"
415d9f5c	37
	38
	39	AliMCEvent::AliMCEvent():
	40	AliVEvent(),
	41	fStack(0),
	42	fMCParticles(new TClonesArray("AliMCParticle",1000)),
	43	fMCParticleMap(0),
	44	fHeader(0),
93df0e9b	45	fTRBuffer(0),
93df0e9b	46	fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
415d9f5c	47	fTreeTR(0),
	48	fTmpTreeTR(0),
	49	fTmpFileTR(0),
	50	fNprimaries(-1),
	51	fNparticles(-1)
	52	{
	53	// Default constructor
	54	}
	55
	56	AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
	57	AliVEvent(mcEvnt),
	58	fStack(0),
	59	fMCParticles(0),
	60	fMCParticleMap(0),
	61	fHeader(0),
93df0e9b	62	fTRBuffer(0),
415d9f5c	63	fTrackReferences(0),
	64	fTreeTR(0),
	65	fTmpTreeTR(0),
	66	fTmpFileTR(0),
	67	fNprimaries(-1),
	68	fNparticles(-1)
	69	{
	70	// Copy constructor
	71	}
	72
	73
	74	AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
	75	{
	76	// assignment operator
	77	if (this!=&mcEvnt) {
	78	AliVEvent::operator=(mcEvnt);
	79	}
	80
	81	return *this;
	82	}
	83
	84	void AliMCEvent::ConnectTreeE (TTree* tree)
	85	{
	86	// Connect the event header tree
	87	tree->SetBranchAddress("Header", &fHeader);
	88	}
	89
	90	void AliMCEvent::ConnectTreeK (TTree* tree)
	91	{
	92	// Connect the kinematics tree to the stack
	93	fStack = fHeader->Stack();
	94	fStack->ConnectTree(tree);
	95	//
	96	// Load the event
	97	fStack->GetEvent();
	98	fNparticles = fStack->GetNtrack();
	99	fNprimaries = fStack->GetNprimary();
	100	AliInfo(Form("AliMCEvent: Number of particles: %5d (all) %5d (primaries)\n",
	101	fNparticles, fNprimaries));
	102
	103	// This is a cache for the TParticles converted to MCParticles on user request
	104	if (fMCParticleMap) {
	105	fMCParticleMap->Clear();
	106	if (fNparticles>0) fMCParticleMap->Expand(fNparticles);}
	107	else
93df0e9b	108	fMCParticleMap = new TRefArray(fNparticles);
93df0e9b	109
415d9f5c	110	}
	111
	112	void AliMCEvent::ConnectTreeTR (TTree* tree)
	113	{
	114	// Connect the track reference tree
	115	fTreeTR = tree;
	116
	117	if (fTreeTR->GetBranch("AliRun")) {
	118	if (fTmpFileTR) {
	119	fTmpFileTR->Close();
	120	delete fTmpFileTR;
	121	}
	122	// This is an old format with one branch per detector not in synch with TreeK
	123	ReorderAndExpandTreeTR();
	124	} else {
	125	// New format
93df0e9b	126	fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
415d9f5c	127	}
	128	}
	129
	130	Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle& particle, TClonesArray& trefs)
	131	{
	132	// Retrieve entry i
	133	if (i < 0 \|\| i >= fNparticles) {
	134	AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
	135	particle = 0;
	136	trefs = 0;
	137	return (-1);
	138	}
	139	particle = fStack->Particle(i);
	140	if (fTreeTR) {
	141	fTreeTR->GetEntry(fStack->TreeKEntry(i));
93df0e9b	142	trefs = fTRBuffer;
415d9f5c	143	return trefs->GetEntries();
	144	} else {
	145	trefs = 0;
	146	return -1;
	147	}
	148	}
	149
	150
	151	void AliMCEvent::Clean()
	152	{
	153	// Clean-up before new trees are connected
	154
	155	if (fHeader) {
	156	delete fHeader;
	157	fHeader = 0;
	158	}
	159
	160	delete fStack;
	161
	162	// Clear TR
93df0e9b	163	if (fTRBuffer) {
	164	fTRBuffer->Clear();
	165	delete fTRBuffer;
	166	fTRBuffer = 0;
415d9f5c	167	}
	168	}
	169
	170	void AliMCEvent::FinishEvent()
	171	{
	172	// Clean-up after event
	173	fStack->Reset(0);
93df0e9b	174	fMCParticles->Clear();
93df0e9b	175	fTrackReferences->Clear();
415d9f5c	176	}
	177
	178
93df0e9b	179
415d9f5c	180	void AliMCEvent::DrawCheck(Int_t i, Int_t search)
	181	{
	182	//
	183	// Simple event display for debugging
	184	if (!fTreeTR) {
	185	AliWarning("No Track Reference information available");
	186	return;
	187	}
	188
	189	if (i > -1 && i < fNparticles) {
	190	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	191	} else {
	192	AliWarning("AliMCEvent::GetEntry: Index out of range");
	193	}
	194
93df0e9b	195	Int_t nh = fTRBuffer->GetEntries();
415d9f5c	196
	197
	198	if (search) {
	199	while(nh <= search && i < fNparticles - 1) {
	200	i++;
	201	fTreeTR->GetEntry(fStack->TreeKEntry(i));
93df0e9b	202	nh = fTRBuffer->GetEntries();
415d9f5c	203	}
	204	printf("Found Hits at %5d\n", i);
	205	}
	206	TParticle* particle = fStack->Particle(i);
	207
	208	TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
	209	Float_t x0 = particle->Vx();
	210	Float_t y0 = particle->Vy();
	211
	212	Float_t x1 = particle->Vx() + particle->Px() * 50.;
	213	Float_t y1 = particle->Vy() + particle->Py() * 50.;
	214
	215	TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
	216	h->Draw();
	217	a->SetLineColor(2);
	218
	219	a->Draw();
	220
	221	for (Int_t ih = 0; ih < nh; ih++) {
93df0e9b	222	AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
415d9f5c	223	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	224	m->Draw();
	225	m->SetMarkerSize(0.4);
	226
	227	}
	228	}
	229
	230
	231	void AliMCEvent::ReorderAndExpandTreeTR()
	232	{
	233	//
	234	// Reorder and expand the track reference tree in order to match the kinematics tree.
	235	// Copy the information from different branches into one
	236	//
	237	// TreeTR
	238
	239	fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
	240	fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
93df0e9b	241	if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
93df0e9b	242	fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
415d9f5c	243
	244
	245	//
	246	// Activate the used branches only. Otherwisw we get a bad memory leak.
	247	fTreeTR->SetBranchStatus("*", 0);
	248	fTreeTR->SetBranchStatus("AliRun.*", 1);
	249	fTreeTR->SetBranchStatus("ITS.*", 1);
	250	fTreeTR->SetBranchStatus("TPC.*", 1);
	251	fTreeTR->SetBranchStatus("TRD.*", 1);
	252	fTreeTR->SetBranchStatus("TOF.*", 1);
	253	fTreeTR->SetBranchStatus("FRAME.*", 1);
	254	fTreeTR->SetBranchStatus("MUON.*", 1);
	255	//
	256	// Connect the active branches
	257	TClonesArray* trefs[7];
	258	for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
	259	if (fTreeTR){
	260	// make branch for central track references
	261	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	262	if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
	263	if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
	264	if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
	265	if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
	266	if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
	267	if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
	268	}
	269
	270	Int_t np = fStack->GetNprimary();
	271	Int_t nt = fTreeTR->GetEntries();
	272
	273	//
	274	// Loop over tracks and find the secondaries with the help of the kine tree
	275	Int_t ifills = 0;
	276	Int_t it = 0;
	277	Int_t itlast = 0;
	278	TParticle* part;
	279
	280	for (Int_t ip = np - 1; ip > -1; ip--) {
	281	part = fStack->Particle(ip);
	282	// printf("Particle %5d %5d %5d %5d %5d %5d \n",
	283	// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
	284	// part->GetLastDaughter(), part->TestBit(kTransportBit));
	285
	286	// Determine range of secondaries produced by this primary during transport
	287	Int_t dau1 = part->GetFirstDaughter();
	288	if (dau1 < np) continue; // This particle has no secondaries produced during transport
	289	Int_t dau2 = -1;
	290	if (dau1 > -1) {
	291	Int_t inext = ip - 1;
	292	while (dau2 < 0) {
	293	if (inext >= 0) {
	294	part = fStack->Particle(inext);
	295	dau2 = part->GetFirstDaughter();
	296	if (dau2 == -1 \|\| dau2 < np) {
	297	dau2 = -1;
	298	} else {
	299	dau2--;
	300	}
	301	} else {
	302	dau2 = fStack->GetNtrack() - 1;
	303	}
	304	inext--;
	305	} // find upper bound
	306	} // dau2 < 0
307
308
309	// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
310	//
311	// Loop over reference hits and find secondary label
312	// First the tricky part: find the entry in treeTR than contains the hits or
313	// make sure that no hits exist.
314	//
315	Bool_t hasHits = kFALSE;
316	Bool_t isOutside = kFALSE;
317
318	it = itlast;
319	while (!hasHits && !isOutside && it < nt) {
320	fTreeTR->GetEntry(it++);
321	for (Int_t ib = 0; ib < 7; ib++) {
322	if (!trefs[ib]) continue;
323	Int_t nh = trefs[ib]->GetEntries();
324	for (Int_t ih = 0; ih < nh; ih++) {
325	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
326	Int_t label = tr->Label();
327	if (label >= dau1 && label <= dau2) {
328	hasHits = kTRUE;
329	itlast = it - 1;
330	break;
331	}
332	if (label > dau2 \|\| label < ip) {
333	isOutside = kTRUE;
334	itlast = it - 1;
335	break;
336	}
337	} // hits
338	if (hasHits \|\| isOutside) break;
339	} // branches
340	} // entries
341
342	if (!hasHits) {
343	// Write empty entries
344	for (Int_t id = dau1; (id <= dau2); id++) {
345	fTmpTreeTR->Fill();
346	ifills++;
347	}
348	} else {
349	// Collect all hits
350	fTreeTR->GetEntry(itlast);
351	for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
352	for (Int_t ib = 0; ib < 7; ib++) {
353	if (!trefs[ib]) continue;
354	Int_t nh = trefs[ib]->GetEntries();
355	for (Int_t ih = 0; ih < nh; ih++) {
356	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
357	Int_t label = tr->Label();
358	// Skip primaries
359	if (label == ip) continue;
360	if (label > dau2 \|\| label < dau1)
361	printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
362	itlast, label, dau1, dau2);
363	if (label == id) {
364	// secondary found
365	tr->SetDetectorId(ib-1);
93df0e9b	366	Int_t nref = fTRBuffer->GetEntriesFast();
93df0e9b	367	TClonesArray &lref = *fTRBuffer;
415d9f5c	368	new(lref[nref]) AliTrackReference(*tr);
	369	}
	370	} // hits
	371	} // branches
	372	fTmpTreeTR->Fill();
93df0e9b	373	fTRBuffer->Clear();
415d9f5c	374	ifills++;
	375	} // daughters
	376	} // has hits
	377	} // tracks
	378
	379	//
	380	// Now loop again and write the primaries
	381	//
	382	it = nt - 1;
	383	for (Int_t ip = 0; ip < np; ip++) {
	384	Int_t labmax = -1;
	385	while (labmax < ip && it > -1) {
	386	fTreeTR->GetEntry(it--);
	387	for (Int_t ib = 0; ib < 7; ib++) {
	388	if (!trefs[ib]) continue;
	389	Int_t nh = trefs[ib]->GetEntries();
	390	//
	391	// Loop over reference hits and find primary labels
	392	for (Int_t ih = 0; ih < nh; ih++) {
	393	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	394	Int_t label = tr->Label();
	395	if (label < np && label > labmax) {
	396	labmax = label;
	397	}
	398
	399	if (label == ip) {
	400	tr->SetDetectorId(ib-1);
93df0e9b	401	Int_t nref = fTRBuffer->GetEntriesFast();
93df0e9b	402	TClonesArray &lref = *fTRBuffer;
415d9f5c	403	new(lref[nref]) AliTrackReference(*tr);
	404	}
	405	} // hits
	406	} // branches
	407	} // entries
	408	it++;
	409	fTmpTreeTR->Fill();
93df0e9b	410	fTRBuffer->Clear();
415d9f5c	411	ifills++;
	412	} // tracks
	413	// Check
	414
	415
	416	// Clean-up
	417	delete fTreeTR; fTreeTR = 0;
	418
	419	for (Int_t ib = 0; ib < 7; ib++) {
	420	if (trefs[ib]) {
	421	trefs[ib]->Clear();
	422	delete trefs[ib];
	423	trefs[ib] = 0;
	424	}
	425	}
	426
	427	if (ifills != fStack->GetNtrack())
	428	printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
	429	ifills, fStack->GetNtrack());
	430
	431	fTmpTreeTR->Write();
	432	fTreeTR = fTmpTreeTR;
	433	}
	434
	435	AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
	436	{
	437	// Get MC Particle i
93df0e9b	438	AliMCParticle *mcParticle = 0;
	439	TParticle *particle = 0;
	440	TClonesArray *trefs = 0;
	441	Int_t ntref = 0;
	442	TRefArray *rarray = 0;
	443
415d9f5c	444	// Out of range check
	445	if (i < 0 \|\| i >= fNparticles) {
	446	AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
	447	mcParticle = 0;
	448	return (mcParticle);
	449	}
	450
	451
93df0e9b	452
415d9f5c	453	//
	454	// First check of the MC Particle has been already cached
	455	if(!fMCParticleMap->At(i)) {
93df0e9b	456	// Get particle from the stack
	457	particle = fStack->Particle(i);
	458	// Get track references from Tree TR
	459	if (fTreeTR) {
	460	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	461	trefs = fTRBuffer;
	462	ntref = trefs->GetEntriesFast();
	463	rarray = new TRefArray(ntref);
	464	Int_t nen = fTrackReferences->GetEntriesFast();
	465	for (Int_t j = 0; j < ntref; j++) {
	466	// Save the track references in a TClonesArray
	467	AliTrackReference* ref = dynamic_cast<AliTrackReference>((fTRBuffer)[j]);
	468	// Save the pointer in a TRefArray
	469	new ((fTrackReferences)[nen]) AliTrackReference(ref);
	470	rarray->AddAt((*fTrackReferences)[nen], j);
	471	nen++;
	472	} // loop over track references for entry i
	473	} // if TreeTR available
	474
415d9f5c	475	Int_t nentries = fMCParticles->GetEntriesFast();
93df0e9b	476	new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray);
415d9f5c	477	mcParticle = dynamic_cast<AliMCParticle>((fMCParticles)[nentries]);
	478	fMCParticleMap->AddAt(mcParticle, i);
	479	} else {
	480	mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
	481	}
	482
	483	return mcParticle;
	484	}
	485
83626b6b	486	AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
415d9f5c	487
415d9f5c	488	ClassImp(AliMCEvent)