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

/**************************************************************************
 * Copyright(c) 1998-1999, 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 AliMCEventHandler
// This class gives access to MC truth during the analysis.
// Monte Carlo truth is containe in the kinematics tree (produced particles) and 
// the tree of reference hits.
//      
// Origin: Andreas Morsch, CERN, andreas.morsch@cern.ch 
//---------------------------------------------------------------------------------


#include "AliMCEventHandler.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliStack.h"
#include "AliAnalysisManager.h"

#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TDirectoryFile.h>
#include <TArrow.h>
#include <TMarker.h>
#include <TH2F.h>


ClassImp(AliMCEventHandler)

AliMCEventHandler::AliMCEventHandler() :
    AliVirtualEventHandler(),
    fFileE(0),
    fFileK(0),
    fFileTR(0),
    fTmpFileTR(0),
    fTreeE(0),
    fTreeK(0),
    fTreeTR(0),
    fTmpTreeTR(0),
    fStack(0),
    fHeader(0),
    fTrackReferences(0),
    fNEvent(-1),
    fEvent(-1),
    fNprimaries(-1),
    fNparticles(-1),
    fPathName("./")
{
    // Default constructor
}

AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
    AliVirtualEventHandler(name, title),
    fFileE(0),
    fFileK(0),
    fFileTR(0),
    fTmpFileTR(0),
    fTreeE(0),
    fTreeK(0),
    fTreeTR(0),
    fTmpTreeTR(0),
    fStack(0),
    fHeader(new AliHeader()),
    fTrackReferences(new TClonesArray("AliTrackReference", 200)),
    fNEvent(-1),
    fEvent(-1),
    fNprimaries(-1),
    fNparticles(-1),
    fPathName("./")
{
    // Constructor
}
AliMCEventHandler::~AliMCEventHandler()
{ 
    // Destructor
    delete fFileE;
    delete fFileK;
    delete fFileTR;
}

Bool_t AliMCEventHandler::InitIO(Option_t* /*opt*/) 
{ 
    // Initialize input
    //
    fFileE = new TFile(Form("%sgalice.root", fPathName));
    if (!fFileE) printf("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName);

    fFileE->GetObject("TE", fTreeE);
    fTreeE->SetBranchAddress("Header", &fHeader);
    fNEvent = fTreeE->GetEntries();
    //
    // Tree K
    fFileK = new TFile(Form("%sKinematics.root", fPathName));
    if (!fFileK) printf("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName);
    //
    // Tree TR
    fFileTR = new TFile(Form("%sTrackRefs.root", fPathName));
    if (!fFileTR) printf("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName);
    //
    // Reset the event number
    fEvent = -1;
    printf("AliMCEventHandler:Number of events in this directory %5d \n", fNEvent);
    return kTRUE;
    
}

Bool_t AliMCEventHandler::BeginEvent()
{ 
    // Read the next event
    fEvent++;
    char folder[20];
    sprintf(folder, "Event%d", fEvent);
    // TreeE
    fTreeE->GetEntry(fEvent);
    fStack = fHeader->Stack();
    // Tree K
    TDirectoryFile* dirK  = 0;
    fFileK->GetObject(folder, dirK);
    dirK->GetObject("TreeK", fTreeK);
    fStack->ConnectTree(fTreeK);
    fStack->GetEvent();
    //Tree TR 
    TDirectoryFile* dirTR = 0;
    fFileTR->GetObject(folder, dirTR);
    dirTR->GetObject("TreeTR", fTreeTR);
    if (fTreeTR->GetBranch("AliRun")) {
	// This is an old format with one branch per detector not in synch with TreeK
	ReorderAndExpandTreeTR();
    } else {
	// New format 
	fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
    }
	
    //
    fNparticles = fStack->GetNtrack();
    fNprimaries = fStack->GetNprimary();
    printf("AliMCEventHandler: Event#: %5d Number of particles %5d \n", fEvent, fNparticles);
    
    return kTRUE;
}

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

void AliMCEventHandler::DrawCheck(Int_t i, Bool_t search)
{
    // Retrieve entry i and draw momentum vector and hits
    if (i > -1 && i < fNparticles) {
	fTreeTR->GetEntry(fStack->TreeKEntry(i));
    } else {
	printf("AliMCEventHandler::GetEntry: Index out of range");
    }

    Int_t nh = fTrackReferences->GetEntries();
    
    
    if (search) {
	while(nh == 0 && i < fNparticles - 1) {
	    i++;
	    fTreeTR->GetEntry(fStack->TreeKEntry(i));
	    nh =  fTrackReferences->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*) fTrackReferences->At(ih);
	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	m->Draw();
	m->SetMarkerSize(0.4);
	
    }
}

Bool_t AliMCEventHandler::Notify()
{
    // Notify about directory change
    //
    // Reconnect trees
    TTree* tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
    TFile *curfile = tree->GetCurrentFile();
    TString fileName(curfile->GetName());
    fileName.ReplaceAll("AliESDs.root", "");
    printf("AliMCEventHandler::Notify() file: %s\n", fileName.Data());
    fPathName = Form("%s",  fileName.Data());
    ResetIO();
    InitIO("");
    return kTRUE;
}
    
void AliMCEventHandler::ResetIO()
{
    // Reset files
    if (fFileE)  delete fFileE;
    if (fFileK)  delete fFileK;
    if (fFileTR) delete fFileTR;
}

			    
Bool_t AliMCEventHandler::FinishEvent()
{
    // Dummy 
    return kTRUE;
}

Bool_t AliMCEventHandler::Terminate()
{ 
    // Dummy 
    return kTRUE;
}

Bool_t AliMCEventHandler::TerminateIO()
{ 
    // Dummy
    return kTRUE;
}
    
void AliMCEventHandler::ReorderAndExpandTreeTR()
{
//
//  Reorder and expand the track reference tree in order to match the kinematics tree.
//  Copy the information from different branches into one
//
//  TreeTR
    if (fTmpTreeTR) delete fTmpTreeTR;
    if (fTmpFileTR) {
	fTmpFileTR->Close();
	delete fTmpFileTR;
    }

    fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
    fTmpTreeTR = new TTree("TreeTR", "Track References");
    if (!fTrackReferences)  fTrackReferences = new TClonesArray("AliTrackReference", 100);
    fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
//
    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;
    
    for (Int_t ip = np - 1; ip > -1; ip--) {
	TParticle *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 =  fTrackReferences->GetEntriesFast();
			    TClonesArray &lref = *fTrackReferences;
			    new(lref[nref]) AliTrackReference(*tr);
			}
		    } // hits
		} // branches
		fTmpTreeTR->Fill();
		fTrackReferences->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 = fTrackReferences->GetEntriesFast();
			TClonesArray &lref = *fTrackReferences;
			new(lref[nref]) AliTrackReference(*tr);
		    }
		} // hits
	    } // branches
	} // entries
	it++;
	fTmpTreeTR->Fill();
	fTrackReferences->Clear();
	ifills++;
    } // tracks
    // Check
    if (ifills != fStack->GetNtrack()) 
	printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", 
	       ifills, fStack->GetNtrack());
    fTreeTR = fTmpTreeTR;
}
Commit	Line	Data
5fe09262	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	// Class AliMCEventHandler
	19	// This class gives access to MC truth during the analysis.
	20	// Monte Carlo truth is containe in the kinematics tree (produced particles) and
	21	// the tree of reference hits.
	22	//
	23	// Origin: Andreas Morsch, CERN, andreas.morsch@cern.ch
	24	//---------------------------------------------------------------------------------
	25
	26
	27
	28	#include "AliMCEventHandler.h"
	29	#include "AliTrackReference.h"
	30	#include "AliHeader.h"
	31	#include "AliStack.h"
	32	#include "AliAnalysisManager.h"
	33
	34	#include <TTree.h>
	35	#include <TFile.h>
	36	#include <TParticle.h>
	37	#include <TClonesArray.h>
	38	#include <TDirectoryFile.h>
	39	#include <TArrow.h>
	40	#include <TMarker.h>
	41	#include <TH2F.h>
	42
	43
	44	ClassImp(AliMCEventHandler)
	45
	46	AliMCEventHandler::AliMCEventHandler() :
	47	AliVirtualEventHandler(),
	48	fFileE(0),
	49	fFileK(0),
	50	fFileTR(0),
	51	fTmpFileTR(0),
	52	fTreeE(0),
	53	fTreeK(0),
	54	fTreeTR(0),
	55	fTmpTreeTR(0),
	56	fStack(0),
	57	fHeader(0),
	58	fTrackReferences(0),
	59	fNEvent(-1),
	60	fEvent(-1),
	61	fNprimaries(-1),
	62	fNparticles(-1),
	63	fPathName("./")
	64	{
65	// Default constructor
66	}
67
68	AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
69	AliVirtualEventHandler(name, title),
70	fFileE(0),
71	fFileK(0),
72	fFileTR(0),
73	fTmpFileTR(0),
74	fTreeE(0),
75	fTreeK(0),
76	fTreeTR(0),
77	fTmpTreeTR(0),
78	fStack(0),
79	fHeader(new AliHeader()),
80	fTrackReferences(new TClonesArray("AliTrackReference", 200)),
81	fNEvent(-1),
82	fEvent(-1),
83	fNprimaries(-1),
84	fNparticles(-1),
85	fPathName("./")
86	{
87	// Constructor
88	}
89	AliMCEventHandler::~AliMCEventHandler()
90	{
91	// Destructor
92	delete fFileE;
93	delete fFileK;
94	delete fFileTR;
95	}
96
97	Bool_t AliMCEventHandler::InitIO(Option_t* /opt/)
98	{
99	// Initialize input
100	//
101	fFileE = new TFile(Form("%sgalice.root", fPathName));
36e4cfbb	102	if (!fFileE) printf("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName);
5fe09262	103
	104	fFileE->GetObject("TE", fTreeE);
	105	fTreeE->SetBranchAddress("Header", &fHeader);
	106	fNEvent = fTreeE->GetEntries();
	107	//
	108	// Tree K
	109	fFileK = new TFile(Form("%sKinematics.root", fPathName));
36e4cfbb	110	if (!fFileK) printf("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName);
5fe09262	111	//
	112	// Tree TR
	113	fFileTR = new TFile(Form("%sTrackRefs.root", fPathName));
36e4cfbb	114	if (!fFileTR) printf("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName);
5fe09262	115	//
	116	// Reset the event number
	117	fEvent = -1;
36e4cfbb	118	printf("AliMCEventHandler:Number of events in this directory %5d \n", fNEvent);
5fe09262	119	return kTRUE;
	120
	121	}
	122
	123	Bool_t AliMCEventHandler::BeginEvent()
	124	{
	125	// Read the next event
5fe09262	126	fEvent++;
	127	char folder[20];
	128	sprintf(folder, "Event%d", fEvent);
	129	// TreeE
	130	fTreeE->GetEntry(fEvent);
	131	fStack = fHeader->Stack();
	132	// Tree K
	133	TDirectoryFile* dirK = 0;
	134	fFileK->GetObject(folder, dirK);
	135	dirK->GetObject("TreeK", fTreeK);
	136	fStack->ConnectTree(fTreeK);
	137	fStack->GetEvent();
	138	//Tree TR
	139	TDirectoryFile* dirTR = 0;
	140	fFileTR->GetObject(folder, dirTR);
	141	dirTR->GetObject("TreeTR", fTreeTR);
	142	if (fTreeTR->GetBranch("AliRun")) {
	143	// This is an old format with one branch per detector not in synch with TreeK
	144	ReorderAndExpandTreeTR();
	145	} else {
	146	// New format
	147	fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
	148	}
	149
	150	//
	151	fNparticles = fStack->GetNtrack();
	152	fNprimaries = fStack->GetNprimary();
36e4cfbb	153	printf("AliMCEventHandler: Event#: %5d Number of particles %5d \n", fEvent, fNparticles);
5fe09262	154
	155	return kTRUE;
	156	}
	157
	158	Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle& particle, TClonesArray& trefs)
	159	{
	160	// Retrieve entry i
	161	if (i > -1 && i < fNparticles) {
	162	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	163	} else {
	164	printf("AliMCEventHandler::GetEntry: Index out of range");
	165	}
	166	particle = fStack->Particle(i);
	167	trefs = fTrackReferences;
5fe09262	168	return trefs->GetEntries();
5fe09262	169	}
5fe09262	170
36e4cfbb	171	void AliMCEventHandler::DrawCheck(Int_t i, Bool_t search)
5fe09262	172	{
	173	// Retrieve entry i and draw momentum vector and hits
	174	if (i > -1 && i < fNparticles) {
	175	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	176	} else {
	177	printf("AliMCEventHandler::GetEntry: Index out of range");
	178	}
	179
36e4cfbb	180	Int_t nh = fTrackReferences->GetEntries();
36e4cfbb	181
5fe09262	182
36e4cfbb	183	if (search) {
	184	while(nh == 0 && i < fNparticles - 1) {
	185	i++;
	186	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	187	nh = fTrackReferences->GetEntries();
	188	}
	189	printf("Found Hits at %5d\n", i);
	190	}
	191	TParticle* particle = fStack->Particle(i);
	192
5fe09262	193	TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
	194	Float_t x0 = particle->Vx();
	195	Float_t y0 = particle->Vy();
	196
	197	Float_t x1 = particle->Vx() + particle->Px() * 50.;
	198	Float_t y1 = particle->Vy() + particle->Py() * 50.;
	199
	200	TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
	201	h->Draw();
	202	a->SetLineColor(2);
	203
	204	a->Draw();
	205
	206	for (Int_t ih = 0; ih < nh; ih++) {
	207	AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
	208	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	209	m->Draw();
	210	m->SetMarkerSize(0.4);
	211
	212	}
	213	}
	214
	215	Bool_t AliMCEventHandler::Notify()
	216	{
	217	// Notify about directory change
36e4cfbb	218	//
36e4cfbb	219	// Reconnect trees
5fe09262	220	TTree* tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
	221	TFile *curfile = tree->GetCurrentFile();
	222	TString fileName(curfile->GetName());
	223	fileName.ReplaceAll("AliESDs.root", "");
	224	printf("AliMCEventHandler::Notify() file: %s\n", fileName.Data());
	225	fPathName = Form("%s", fileName.Data());
	226	ResetIO();
	227	InitIO("");
	228	return kTRUE;
	229	}
	230
	231	void AliMCEventHandler::ResetIO()
	232	{
36e4cfbb	233	// Reset files
5fe09262	234	if (fFileE) delete fFileE;
	235	if (fFileK) delete fFileK;
	236	if (fFileTR) delete fFileTR;
	237	}
	238
	239
	240	Bool_t AliMCEventHandler::FinishEvent()
	241	{
	242	// Dummy
	243	return kTRUE;
	244	}
	245
	246	Bool_t AliMCEventHandler::Terminate()
	247	{
	248	// Dummy
	249	return kTRUE;
	250	}
	251
	252	Bool_t AliMCEventHandler::TerminateIO()
	253	{
	254	// Dummy
	255	return kTRUE;
	256	}
	257
	258	void AliMCEventHandler::ReorderAndExpandTreeTR()
	259	{
	260	//
	261	// Reorder and expand the track reference tree in order to match the kinematics tree.
	262	// Copy the information from different branches into one
	263	//
	264	// TreeTR
	265	if (fTmpTreeTR) delete fTmpTreeTR;
	266	if (fTmpFileTR) {
	267	fTmpFileTR->Close();
	268	delete fTmpFileTR;
	269	}
	270
	271	fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
	272	fTmpTreeTR = new TTree("TreeTR", "Track References");
	273	if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
	274	fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
	275	//
	276	TClonesArray* trefs[7];
	277	for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
	278	if (fTreeTR){
	279	// make branch for central track references
	280	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	281	if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
	282	if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
	283	if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
	284	if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
	285	if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
	286	if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
	287	}
	288
	289	Int_t np = fStack->GetNprimary();
	290	Int_t nt = fTreeTR->GetEntries();
	291	//
	292	// Loop over tracks and find the secondaries with the help of the kine tree
	293	Int_t ifills = 0;
	294	Int_t it = 0;
	295	Int_t itlast = 0;
	296
	297	for (Int_t ip = np - 1; ip > -1; ip--) {
298	TParticle *part = fStack->Particle(ip);
299	// printf("Particle %5d %5d %5d %5d %5d %5d \n",
300	// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
301	// part->GetLastDaughter(), part->TestBit(kTransportBit));
36e4cfbb	302
36e4cfbb	303	// Determine range of secondaries produced by this primary during transport
5fe09262	304	Int_t dau1 = part->GetFirstDaughter();
36e4cfbb	305	if (dau1 < np) continue; // This particle has no secondaries produced during transport
5fe09262	306	Int_t dau2 = -1;
5fe09262	307	if (dau1 > -1) {
	308	Int_t inext = ip - 1;
	309	while (dau2 < 0) {
	310	if (inext >= 0) {
	311	part = fStack->Particle(inext);
	312	dau2 = part->GetFirstDaughter();
	313	if (dau2 == -1 \|\| dau2 < np) {
	314	dau2 = -1;
	315	} else {
	316	dau2--;
	317	}
	318	} else {
	319	dau2 = fStack->GetNtrack() - 1;
	320	}
	321	inext--;
	322	} // find upper bound
	323	} // dau2 < 0
36e4cfbb	324
5fe09262	325
5fe09262	326	// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
36e4cfbb	327	//
	328	// Loop over reference hits and find secondary label
	329	// First the tricky part: find the entry in treeTR than contains the hits or
	330	// make sure that no hits exist.
	331	//
5fe09262	332	Bool_t hasHits = kFALSE;
5fe09262	333	Bool_t isOutside = kFALSE;
36e4cfbb	334
5fe09262	335	it = itlast;
5fe09262	336	while (!hasHits && !isOutside && it < nt) {
	337	fTreeTR->GetEntry(it++);
	338	for (Int_t ib = 0; ib < 7; ib++) {
	339	if (!trefs[ib]) continue;
	340	Int_t nh = trefs[ib]->GetEntries();
	341	for (Int_t ih = 0; ih < nh; ih++) {
	342	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	343	Int_t label = tr->Label();
	344	if (label >= dau1 && label <= dau2) {
	345	hasHits = kTRUE;
	346	itlast = it - 1;
36e4cfbb	347	break;
5fe09262	348	}
	349	if (label > dau2 \|\| label < ip) {
	350	isOutside = kTRUE;
	351	itlast = it - 1;
36e4cfbb	352	break;
5fe09262	353	}
5fe09262	354	} // hits
36e4cfbb	355	if (hasHits \|\| isOutside) break;
5fe09262	356	} // branches
	357	} // entries
	358
	359	if (!hasHits) {
36e4cfbb	360	// Write empty entries
5fe09262	361	for (Int_t id = dau1; (id <= dau2); id++) {
	362	fTmpTreeTR->Fill();
	363	ifills++;
	364	}
	365	} else {
36e4cfbb	366	// Collect all hits
5fe09262	367	fTreeTR->GetEntry(itlast);
	368	for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
	369	for (Int_t ib = 0; ib < 7; ib++) {
	370	if (!trefs[ib]) continue;
	371	Int_t nh = trefs[ib]->GetEntries();
	372	for (Int_t ih = 0; ih < nh; ih++) {
	373	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	374	Int_t label = tr->Label();
	375	// Skip primaries
	376	if (label == ip) continue;
	377	if (label > dau2 \|\| label < dau1)
36e4cfbb	378	printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
36e4cfbb	379	itlast, label, dau1, dau2);
5fe09262	380	if (label == id) {
	381	// secondary found
	382	tr->SetDetectorId(ib-1);
	383	Int_t nref = fTrackReferences->GetEntriesFast();
	384	TClonesArray &lref = *fTrackReferences;
	385	new(lref[nref]) AliTrackReference(*tr);
	386	}
	387	} // hits
	388	} // branches
	389	fTmpTreeTR->Fill();
	390	fTrackReferences->Clear();
	391	ifills++;
	392	} // daughters
	393	} // has hits
	394	} // tracks
5fe09262	395	//
5fe09262	396	// Now loop again and write the primaries
36e4cfbb	397	//
5fe09262	398	it = nt - 1;
	399	for (Int_t ip = 0; ip < np; ip++) {
	400	Int_t labmax = -1;
	401	while (labmax < ip && it > -1) {
	402	fTreeTR->GetEntry(it--);
	403	for (Int_t ib = 0; ib < 7; ib++) {
	404	if (!trefs[ib]) continue;
	405	Int_t nh = trefs[ib]->GetEntries();
	406	//
	407	// Loop over reference hits and find primary labels
	408	for (Int_t ih = 0; ih < nh; ih++) {
	409	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	410	Int_t label = tr->Label();
	411	if (label < np && label > labmax) {
	412	labmax = label;
	413	}
	414
	415	if (label == ip) {
	416	tr->SetDetectorId(ib-1);
	417	Int_t nref = fTrackReferences->GetEntriesFast();
	418	TClonesArray &lref = *fTrackReferences;
	419	new(lref[nref]) AliTrackReference(*tr);
	420	}
	421	} // hits
	422	} // branches
	423	} // entries
	424	it++;
	425	fTmpTreeTR->Fill();
	426	fTrackReferences->Clear();
	427	ifills++;
	428	} // tracks
	429	// Check
5fe09262	430	if (ifills != fStack->GetNtrack())
36e4cfbb	431	printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
36e4cfbb	432	ifills, fStack->GetNtrack());
5fe09262	433	fTreeTR = fTmpTreeTR;
5fe09262	434	}