[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("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("Kinematics.root not found in directory %s ! \n", fPathName);
    //
    // Tree TR
    fFileTR = new TFile(Form("%sTrackRefs.root", fPathName));
    if (!fFileTR) printf("TrackRefs.root not found in directory %s ! \n", fPathName);
    //
    // Reset the event number
    fEvent = -1;
    printf("Number of events in this directory %5d \n", fNEvent);
    return kTRUE;
    
}

Bool_t AliMCEventHandler::BeginEvent()
{ 
    // Read the next event
    printf("AliMCEventHandler::BeginEvent called \n");
    
    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("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;
    printf("Returning %5d \n", particle->GetPdgCode());
    
    return trefs->GetEntries();
    
}

void AliMCEventHandler::DrawCheck(Int_t i)
{
    // 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");
    }

    TParticle* particle = fStack->Particle(i);
    Int_t nh =  fTrackReferences->GetEntries();
    
    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
    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()
{
    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));
	
	// Skip primaries that have not been transported
	Int_t dau1  = part->GetFirstDaughter();
	Int_t dau2  = -1;
	// Determine range of secondaries produced by this primary
	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
	Bool_t hasHits   = kFALSE;
	Bool_t isOutside = kFALSE;
	it = itlast;
	if (dau1 < np) continue;
	
	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;
		    }
		    if (label > dau2 || label < ip) {
			isOutside = kTRUE;
			itlast = it - 1;
		    }
		} // hits
	    } // branches
	} // entries

	if (!hasHits) {
	    for (Int_t id = dau1; (id <= dau2); id++) {
		fTmpTreeTR->Fill();
		ifills++;
	    } 
	} else {
	    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("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
    printf("Secondaries %5d %5d \n", ifills, fStack->GetNtrack() - fStack->GetNprimary());
    //
    // 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
    printf("All %5d %5d \n", ifills, fStack->GetNtrack());
    if (ifills != fStack->GetNtrack()) 
	printf("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));
102	if (!fFileE) printf("galice.root not found in directory %s ! \n", fPathName);
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));
110	if (!fFileK) printf("Kinematics.root not found in directory %s ! \n", fPathName);
111	//
112	// Tree TR
113	fFileTR = new TFile(Form("%sTrackRefs.root", fPathName));
114	if (!fFileTR) printf("TrackRefs.root not found in directory %s ! \n", fPathName);
115	//
116	// Reset the event number
117	fEvent = -1;
118	printf("Number of events in this directory %5d \n", fNEvent);
119	return kTRUE;
120
121	}
122
123	Bool_t AliMCEventHandler::BeginEvent()
124	{
125	// Read the next event
126	printf("AliMCEventHandler::BeginEvent called \n");
127
128	fEvent++;
129	char folder[20];
130	sprintf(folder, "Event%d", fEvent);
131	// TreeE
132	fTreeE->GetEntry(fEvent);
133	fStack = fHeader->Stack();
134	// Tree K
135	TDirectoryFile* dirK = 0;
136	fFileK->GetObject(folder, dirK);
137	dirK->GetObject("TreeK", fTreeK);
138	fStack->ConnectTree(fTreeK);
139	fStack->GetEvent();
140	//Tree TR
141	TDirectoryFile* dirTR = 0;
142	fFileTR->GetObject(folder, dirTR);
143	dirTR->GetObject("TreeTR", fTreeTR);
144	if (fTreeTR->GetBranch("AliRun")) {
145	// This is an old format with one branch per detector not in synch with TreeK
146	ReorderAndExpandTreeTR();
147	} else {
148	// New format
149	fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
150	}
151
152	//
153	fNparticles = fStack->GetNtrack();
154	fNprimaries = fStack->GetNprimary();
155	printf("Event#: %5d Number of particles %5d \n", fEvent, fNparticles);
156
157	return kTRUE;
158	}
159
160	Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle& particle, TClonesArray& trefs)
161	{
162	// Retrieve entry i
163	if (i > -1 && i < fNparticles) {
164	fTreeTR->GetEntry(fStack->TreeKEntry(i));
165	} else {
166	printf("AliMCEventHandler::GetEntry: Index out of range");
167	}
168	particle = fStack->Particle(i);
169	trefs = fTrackReferences;
170	printf("Returning %5d \n", particle->GetPdgCode());
171
172	return trefs->GetEntries();
173
174	}
175
176	void AliMCEventHandler::DrawCheck(Int_t i)
177	{
178	// Retrieve entry i and draw momentum vector and hits
179	if (i > -1 && i < fNparticles) {
180	fTreeTR->GetEntry(fStack->TreeKEntry(i));
181	} else {
182	printf("AliMCEventHandler::GetEntry: Index out of range");
183	}
184
185	TParticle* particle = fStack->Particle(i);
186	Int_t nh = fTrackReferences->GetEntries();
187
188	TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
189	Float_t x0 = particle->Vx();
190	Float_t y0 = particle->Vy();
191
192	Float_t x1 = particle->Vx() + particle->Px() * 50.;
193	Float_t y1 = particle->Vy() + particle->Py() * 50.;
194
195	TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
196	h->Draw();
197	a->SetLineColor(2);
198
199	a->Draw();
200
201	for (Int_t ih = 0; ih < nh; ih++) {
202	AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
203	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
204	m->Draw();
205	m->SetMarkerSize(0.4);
206
207	}
208	}
209
210	Bool_t AliMCEventHandler::Notify()
211	{
212	// Notify about directory change
213	TTree* tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
214	TFile *curfile = tree->GetCurrentFile();
215	TString fileName(curfile->GetName());
216	fileName.ReplaceAll("AliESDs.root", "");
217	printf("AliMCEventHandler::Notify() file: %s\n", fileName.Data());
218	fPathName = Form("%s", fileName.Data());
219	ResetIO();
220	InitIO("");
221	return kTRUE;
222	}
223
224	void AliMCEventHandler::ResetIO()
225	{
226	if (fFileE) delete fFileE;
227	if (fFileK) delete fFileK;
228	if (fFileTR) delete fFileTR;
229	}
230
231
232	Bool_t AliMCEventHandler::FinishEvent()
233	{
234	// Dummy
235	return kTRUE;
236	}
237
238	Bool_t AliMCEventHandler::Terminate()
239	{
240	// Dummy
241	return kTRUE;
242	}
243
244	Bool_t AliMCEventHandler::TerminateIO()
245	{
246	// Dummy
247	return kTRUE;
248	}
249
250	void AliMCEventHandler::ReorderAndExpandTreeTR()
251	{
252	//
253	// Reorder and expand the track reference tree in order to match the kinematics tree.
254	// Copy the information from different branches into one
255	//
256	// TreeTR
257	if (fTmpTreeTR) delete fTmpTreeTR;
258	if (fTmpFileTR) {
259	fTmpFileTR->Close();
260	delete fTmpFileTR;
261	}
262
263	fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
264	fTmpTreeTR = new TTree("TreeTR", "Track References");
265	if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
266	fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
267	//
268	TClonesArray* trefs[7];
269	for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
270	if (fTreeTR){
271	// make branch for central track references
272	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
273	if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
274	if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
275	if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
276	if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
277	if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
278	if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
279	}
280
281	Int_t np = fStack->GetNprimary();
282	Int_t nt = fTreeTR->GetEntries();
283	//
284	// Loop over tracks and find the secondaries with the help of the kine tree
285	Int_t ifills = 0;
286	Int_t it = 0;
287	Int_t itlast = 0;
288
289	for (Int_t ip = np - 1; ip > -1; ip--) {
290	TParticle *part = fStack->Particle(ip);
291	// printf("Particle %5d %5d %5d %5d %5d %5d \n",
292	// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
293	// part->GetLastDaughter(), part->TestBit(kTransportBit));
294
295	// Skip primaries that have not been transported
296	Int_t dau1 = part->GetFirstDaughter();
297	Int_t dau2 = -1;
298	// Determine range of secondaries produced by this primary
299	if (dau1 > -1) {
300	Int_t inext = ip - 1;
301	while (dau2 < 0) {
302	if (inext >= 0) {
303	part = fStack->Particle(inext);
304	dau2 = part->GetFirstDaughter();
305	if (dau2 == -1 \|\| dau2 < np) {
306	dau2 = -1;
307	} else {
308	dau2--;
309	}
310	} else {
311	dau2 = fStack->GetNtrack() - 1;
312	}
313	inext--;
314	} // find upper bound
315	} // dau2 < 0
316
317	// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
318
319	//
320	// Loop over reference hits and find secondary label
321	Bool_t hasHits = kFALSE;
322	Bool_t isOutside = kFALSE;
323	it = itlast;
324	if (dau1 < np) continue;
325
326	while (!hasHits && !isOutside && it < nt) {
327	fTreeTR->GetEntry(it++);
328	for (Int_t ib = 0; ib < 7; ib++) {
329	if (!trefs[ib]) continue;
330	Int_t nh = trefs[ib]->GetEntries();
331	for (Int_t ih = 0; ih < nh; ih++) {
332	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
333	Int_t label = tr->Label();
334	if (label >= dau1 && label <= dau2) {
335	hasHits = kTRUE;
336	itlast = it - 1;
337	}
338	if (label > dau2 \|\| label < ip) {
339	isOutside = kTRUE;
340	itlast = it - 1;
341	}
342	} // hits
343	} // branches
344	} // entries
345
346	if (!hasHits) {
347	for (Int_t id = dau1; (id <= dau2); id++) {
348	fTmpTreeTR->Fill();
349	ifills++;
350	}
351	} else {
352	fTreeTR->GetEntry(itlast);
353	for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
354	for (Int_t ib = 0; ib < 7; ib++) {
355	if (!trefs[ib]) continue;
356	Int_t nh = trefs[ib]->GetEntries();
357	for (Int_t ih = 0; ih < nh; ih++) {
358	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
359	Int_t label = tr->Label();
360	// Skip primaries
361	if (label == ip) continue;
362	if (label > dau2 \|\| label < dau1)
363	printf("Track Reference Label out of range !: %5d %5d %5d %5d \n", itlast, label, dau1, dau2);
364	if (label == id) {
365	// secondary found
366	tr->SetDetectorId(ib-1);
367	Int_t nref = fTrackReferences->GetEntriesFast();
368	TClonesArray &lref = *fTrackReferences;
369	new(lref[nref]) AliTrackReference(*tr);
370	}
371	} // hits
372	} // branches
373	fTmpTreeTR->Fill();
374	fTrackReferences->Clear();
375	ifills++;
376	} // daughters
377	} // has hits
378	} // tracks
379	printf("Secondaries %5d %5d \n", ifills, fStack->GetNtrack() - fStack->GetNprimary());
380	//
381	// Now loop again and write the primaries
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);
401	Int_t nref = fTrackReferences->GetEntriesFast();
402	TClonesArray &lref = *fTrackReferences;
403	new(lref[nref]) AliTrackReference(*tr);
404	}
405	} // hits
406	} // branches
407	} // entries
408	it++;
409	fTmpTreeTR->Fill();
410	fTrackReferences->Clear();
411	ifills++;
412	} // tracks
413	// Check
414	printf("All %5d %5d \n", ifills, fStack->GetNtrack());
415	if (ifills != fStack->GetNtrack())
416	printf("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, fStack->GetNtrack());
417	fTreeTR = fTmpTreeTR;
418
419	}