[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 "AliLog.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() :
    AliVEventHandler(),
    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("./"),
    fExtension(""),
    fFileNumber(0),
    fEventsPerFile(0)
{
    // Default constructor
}

AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
    AliVEventHandler(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("./"),
    fExtension(""),
    fFileNumber(0),
    fEventsPerFile(0)
{
    // Constructor
}
AliMCEventHandler::~AliMCEventHandler()
{ 
    // Destructor
    delete fFileE;
    delete fFileK;
    delete fFileTR;
}

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

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

Bool_t AliMCEventHandler::GetEvent(Int_t iev)
{
    // Load the event number iev
    //
    // Calculate the file number
    Int_t inew  = iev/fEventsPerFile;
    if (inew != fFileNumber) {
	fFileNumber = inew;
	if (!OpenFile(fFileNumber)){
	    return kFALSE;
	}
    }
    // Folder name
    char folder[20];
    sprintf(folder, "Event%d", iev);

    // TreeE
    fTreeE->GetEntry(iev);
    fStack = fHeader->Stack();
    // Tree K
    TDirectoryFile* dirK  = 0;
    fFileK->GetObject(folder, dirK);
    if (!dirK) {
	AliWarning(Form("AliMCEventHandler: Event #%5d not found\n", iev));
	return kFALSE;
    }
    dirK->GetObject("TreeK", fTreeK);
    fStack->ConnectTree(fTreeK);
    fStack->GetEvent();
    //Tree TR 
    if (fFileTR) {
	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();
    AliInfo(Form("AliMCEventHandler: Event#: %5d Number of particles %5d \n", fEvent, fNparticles));
    
    return kTRUE;
}

Bool_t AliMCEventHandler::OpenFile(Int_t i)
{
    // Open file i
    Bool_t ok = kTRUE;
    if (i > 0) {
	fExtension = Form("%d", i);
    } else {
	fExtension = "";
    }
    
    
    delete fFileK;
    fFileK = new TFile(Form("%sKinematics%s.root", fPathName, fExtension));
    if (!fFileK) {
	AliFatal(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName));
	ok = kFALSE;
    }
    
    delete fFileTR;
    fFileTR = new TFile(Form("%sTrackRefs%s.root", fPathName, fExtension));
    if (!fFileTR) {
	AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName));
	ok = kFALSE;
    }
    
    return ok;
}

Bool_t AliMCEventHandler::BeginEvent()
{ 
    // Read the next event
    fEvent++;
    if (fEvent >= fNEvent) {
	AliWarning(Form("AliMCEventHandler: Event number out of range %5d\n", fEvent));
	return kFALSE;
    }
    return GetEvent(fEvent);
}

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

void AliMCEventHandler::DrawCheck(Int_t i, Bool_t search)
{
    // Retrieve entry i and draw momentum vector and hits
    if (!fFileTR) {
	AliWarning("No Track Reference information available");
	return;
    } 

    if (i > -1 && i < fNparticles) {
	fTreeTR->GetEntry(fStack->TreeKEntry(i));
    } else {
	AliWarning("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(const char *path)
{
    // Notify about directory change
    // The directory is taken from the 'path' argument
    // Reconnect trees

    printf("AliMCEventHandler::Notify() file: %s\n", path);
    fPathName = Form("%s",  path);
    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) fTmpTreeTR->Delete();
    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"
2c36081e	32	#include "AliLog.h"
5fe09262	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() :
d2f1d9ef	47	AliVEventHandler(),
5fe09262	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),
9aea8469	63	fPathName("./"),
	64	fExtension(""),
	65	fFileNumber(0),
	66	fEventsPerFile(0)
5fe09262	67	{
	68	// Default constructor
	69	}
	70
	71	AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
d2f1d9ef	72	AliVEventHandler(name, title),
5fe09262	73	fFileE(0),
	74	fFileK(0),
	75	fFileTR(0),
	76	fTmpFileTR(0),
	77	fTreeE(0),
	78	fTreeK(0),
	79	fTreeTR(0),
	80	fTmpTreeTR(0),
	81	fStack(0),
	82	fHeader(new AliHeader()),
	83	fTrackReferences(new TClonesArray("AliTrackReference", 200)),
	84	fNEvent(-1),
	85	fEvent(-1),
	86	fNprimaries(-1),
	87	fNparticles(-1),
9aea8469	88	fPathName("./"),
	89	fExtension(""),
	90	fFileNumber(0),
	91	fEventsPerFile(0)
5fe09262	92	{
	93	// Constructor
	94	}
	95	AliMCEventHandler::~AliMCEventHandler()
	96	{
	97	// Destructor
	98	delete fFileE;
	99	delete fFileK;
	100	delete fFileTR;
	101	}
	102
	103	Bool_t AliMCEventHandler::InitIO(Option_t* /opt/)
	104	{
	105	// Initialize input
	106	//
0c476ce8	107	fFileE = TFile::Open(Form("%sgalice.root", fPathName));
2c36081e	108	if (!fFileE) AliFatal(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName));
5fe09262	109
	110	fFileE->GetObject("TE", fTreeE);
	111	fTreeE->SetBranchAddress("Header", &fHeader);
	112	fNEvent = fTreeE->GetEntries();
	113	//
	114	// Tree K
0c476ce8	115	fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName, fExtension));
2c36081e	116	if (!fFileK) AliFatal(Form("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName));
9aea8469	117	fEventsPerFile = fFileK->GetNkeys() - fFileK->GetNProcessIDs();
5fe09262	118	//
5fe09262	119	// Tree TR
0c476ce8	120	fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName, fExtension));
2c36081e	121	if (!fFileTR) AliWarning(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName));
5fe09262	122	//
5fe09262	123	// Reset the event number
2c36081e	124	fEvent = -1;
2c36081e	125	fFileNumber = 0;
9aea8469	126
2c36081e	127	AliInfo(Form("AliMCEventHandler:Number of events in this directory %5d \n", fNEvent));
5fe09262	128	return kTRUE;
	129
	130	}
	131
9aea8469	132	Bool_t AliMCEventHandler::GetEvent(Int_t iev)
	133	{
	134	// Load the event number iev
2c36081e	135	//
2c36081e	136	// Calculate the file number
9aea8469	137	Int_t inew = iev/fEventsPerFile;
	138	if (inew != fFileNumber) {
	139	fFileNumber = inew;
	140	if (!OpenFile(fFileNumber)){
	141	return kFALSE;
	142	}
	143	}
2c36081e	144	// Folder name
5fe09262	145	char folder[20];
9aea8469	146	sprintf(folder, "Event%d", iev);
2c36081e	147
5fe09262	148	// TreeE
9aea8469	149	fTreeE->GetEntry(iev);
5fe09262	150	fStack = fHeader->Stack();
	151	// Tree K
	152	TDirectoryFile* dirK = 0;
	153	fFileK->GetObject(folder, dirK);
9aea8469	154	if (!dirK) {
2c36081e	155	AliWarning(Form("AliMCEventHandler: Event #%5d not found\n", iev));
9aea8469	156	return kFALSE;
9aea8469	157	}
5fe09262	158	dirK->GetObject("TreeK", fTreeK);
	159	fStack->ConnectTree(fTreeK);
	160	fStack->GetEvent();
	161	//Tree TR
07bd4750	162	if (fFileTR) {
2c36081e	163	TDirectoryFile* dirTR = 0;
	164	fFileTR->GetObject(folder, dirTR);
	165	dirTR->GetObject("TreeTR", fTreeTR);
	166	if (fTreeTR->GetBranch("AliRun")) {
	167	// This is an old format with one branch per detector not in synch with TreeK
	168	ReorderAndExpandTreeTR();
	169	} else {
	170	// New format
	171	fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
	172	}
5fe09262	173	}
2c36081e	174
5fe09262	175	//
	176	fNparticles = fStack->GetNtrack();
	177	fNprimaries = fStack->GetNprimary();
2c36081e	178	AliInfo(Form("AliMCEventHandler: Event#: %5d Number of particles %5d \n", fEvent, fNparticles));
5fe09262	179
5fe09262	180	return kTRUE;
9aea8469	181	}
	182
	183	Bool_t AliMCEventHandler::OpenFile(Int_t i)
	184	{
	185	// Open file i
	186	Bool_t ok = kTRUE;
	187	if (i > 0) {
	188	fExtension = Form("%d", i);
	189	} else {
	190	fExtension = "";
	191	}
	192
	193
	194	delete fFileK;
	195	fFileK = new TFile(Form("%sKinematics%s.root", fPathName, fExtension));
	196	if (!fFileK) {
2c36081e	197	AliFatal(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName));
9aea8469	198	ok = kFALSE;
	199	}
	200
	201	delete fFileTR;
	202	fFileTR = new TFile(Form("%sTrackRefs%s.root", fPathName, fExtension));
	203	if (!fFileTR) {
2c36081e	204	AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName));
9aea8469	205	ok = kFALSE;
	206	}
	207
	208	return ok;
	209	}
	210
	211	Bool_t AliMCEventHandler::BeginEvent()
	212	{
	213	// Read the next event
	214	fEvent++;
	215	if (fEvent >= fNEvent) {
2c36081e	216	AliWarning(Form("AliMCEventHandler: Event number out of range %5d\n", fEvent));
9aea8469	217	return kFALSE;
	218	}
	219	return GetEvent(fEvent);
5fe09262	220	}
	221
	222	Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle& particle, TClonesArray& trefs)
	223	{
	224	// Retrieve entry i
	225	if (i > -1 && i < fNparticles) {
2c36081e	226	AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
	227	particle = 0;
	228	trefs = 0;
	229	return (-1);
	230	}
	231	particle = fStack->Particle(i);
	232	if (fFileTR) {
5fe09262	233	fTreeTR->GetEntry(fStack->TreeKEntry(i));
2c36081e	234	trefs = fTrackReferences;
2c36081e	235	return trefs->GetEntries();
5fe09262	236	} else {
2c36081e	237	trefs = 0;
2c36081e	238	return -1;
5fe09262	239	}
5fe09262	240	}
5fe09262	241
36e4cfbb	242	void AliMCEventHandler::DrawCheck(Int_t i, Bool_t search)
5fe09262	243	{
5fe09262	244	// Retrieve entry i and draw momentum vector and hits
2c36081e	245	if (!fFileTR) {
	246	AliWarning("No Track Reference information available");
	247	return;
	248	}
	249
5fe09262	250	if (i > -1 && i < fNparticles) {
	251	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	252	} else {
2c36081e	253	AliWarning("AliMCEventHandler::GetEntry: Index out of range");
5fe09262	254	}
2c36081e	255
36e4cfbb	256	Int_t nh = fTrackReferences->GetEntries();
36e4cfbb	257
5fe09262	258
36e4cfbb	259	if (search) {
	260	while(nh == 0 && i < fNparticles - 1) {
	261	i++;
	262	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	263	nh = fTrackReferences->GetEntries();
	264	}
	265	printf("Found Hits at %5d\n", i);
	266	}
	267	TParticle* particle = fStack->Particle(i);
	268
5fe09262	269	TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
	270	Float_t x0 = particle->Vx();
	271	Float_t y0 = particle->Vy();
	272
	273	Float_t x1 = particle->Vx() + particle->Px() * 50.;
	274	Float_t y1 = particle->Vy() + particle->Py() * 50.;
	275
	276	TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
	277	h->Draw();
	278	a->SetLineColor(2);
	279
	280	a->Draw();
	281
	282	for (Int_t ih = 0; ih < nh; ih++) {
	283	AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
	284	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	285	m->Draw();
	286	m->SetMarkerSize(0.4);
	287
	288	}
	289	}
	290
890126ab	291	Bool_t AliMCEventHandler::Notify(const char *path)
5fe09262	292	{
5fe09262	293	// Notify about directory change
890126ab	294	// The directory is taken from the 'path' argument
36e4cfbb	295	// Reconnect trees
890126ab	296
	297	printf("AliMCEventHandler::Notify() file: %s\n", path);
	298	fPathName = Form("%s", path);
5fe09262	299	ResetIO();
	300	InitIO("");
	301	return kTRUE;
	302	}
	303
	304	void AliMCEventHandler::ResetIO()
	305	{
36e4cfbb	306	// Reset files
5fe09262	307	if (fFileE) delete fFileE;
	308	if (fFileK) delete fFileK;
	309	if (fFileTR) delete fFileTR;
	310	}
	311
	312
	313	Bool_t AliMCEventHandler::FinishEvent()
	314	{
	315	// Dummy
	316	return kTRUE;
	317	}
	318
	319	Bool_t AliMCEventHandler::Terminate()
	320	{
	321	// Dummy
	322	return kTRUE;
	323	}
	324
	325	Bool_t AliMCEventHandler::TerminateIO()
	326	{
	327	// Dummy
	328	return kTRUE;
	329	}
	330
	331	void AliMCEventHandler::ReorderAndExpandTreeTR()
	332	{
	333	//
	334	// Reorder and expand the track reference tree in order to match the kinematics tree.
	335	// Copy the information from different branches into one
	336	//
	337	// TreeTR
777b7d9c	338	if (fTmpTreeTR) fTmpTreeTR->Delete();
5fe09262	339	if (fTmpFileTR) {
	340	fTmpFileTR->Close();
	341	delete fTmpFileTR;
	342	}
	343
	344	fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
	345	fTmpTreeTR = new TTree("TreeTR", "Track References");
	346	if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
	347	fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
	348	//
	349	TClonesArray* trefs[7];
	350	for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
	351	if (fTreeTR){
	352	// make branch for central track references
	353	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	354	if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
	355	if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
	356	if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
	357	if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
	358	if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
	359	if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
	360	}
	361
	362	Int_t np = fStack->GetNprimary();
	363	Int_t nt = fTreeTR->GetEntries();
	364	//
	365	// Loop over tracks and find the secondaries with the help of the kine tree
	366	Int_t ifills = 0;
	367	Int_t it = 0;
	368	Int_t itlast = 0;
	369
	370	for (Int_t ip = np - 1; ip > -1; ip--) {
	371	TParticle *part = fStack->Particle(ip);
	372	// printf("Particle %5d %5d %5d %5d %5d %5d \n",
	373	// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
	374	// part->GetLastDaughter(), part->TestBit(kTransportBit));
36e4cfbb	375
36e4cfbb	376	// Determine range of secondaries produced by this primary during transport
5fe09262	377	Int_t dau1 = part->GetFirstDaughter();
36e4cfbb	378	if (dau1 < np) continue; // This particle has no secondaries produced during transport
5fe09262	379	Int_t dau2 = -1;
5fe09262	380	if (dau1 > -1) {
	381	Int_t inext = ip - 1;
	382	while (dau2 < 0) {
	383	if (inext >= 0) {
	384	part = fStack->Particle(inext);
	385	dau2 = part->GetFirstDaughter();
	386	if (dau2 == -1 \|\| dau2 < np) {
	387	dau2 = -1;
	388	} else {
	389	dau2--;
	390	}
	391	} else {
	392	dau2 = fStack->GetNtrack() - 1;
	393	}
	394	inext--;
	395	} // find upper bound
	396	} // dau2 < 0
36e4cfbb	397
5fe09262	398
5fe09262	399	// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
36e4cfbb	400	//
	401	// Loop over reference hits and find secondary label
	402	// First the tricky part: find the entry in treeTR than contains the hits or
	403	// make sure that no hits exist.
	404	//
5fe09262	405	Bool_t hasHits = kFALSE;
5fe09262	406	Bool_t isOutside = kFALSE;
36e4cfbb	407
5fe09262	408	it = itlast;
5fe09262	409	while (!hasHits && !isOutside && it < nt) {
	410	fTreeTR->GetEntry(it++);
	411	for (Int_t ib = 0; ib < 7; ib++) {
	412	if (!trefs[ib]) continue;
	413	Int_t nh = trefs[ib]->GetEntries();
	414	for (Int_t ih = 0; ih < nh; ih++) {
	415	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	416	Int_t label = tr->Label();
	417	if (label >= dau1 && label <= dau2) {
	418	hasHits = kTRUE;
	419	itlast = it - 1;
36e4cfbb	420	break;
5fe09262	421	}
	422	if (label > dau2 \|\| label < ip) {
	423	isOutside = kTRUE;
	424	itlast = it - 1;
36e4cfbb	425	break;
5fe09262	426	}
5fe09262	427	} // hits
36e4cfbb	428	if (hasHits \|\| isOutside) break;
5fe09262	429	} // branches
	430	} // entries
	431
	432	if (!hasHits) {
36e4cfbb	433	// Write empty entries
5fe09262	434	for (Int_t id = dau1; (id <= dau2); id++) {
	435	fTmpTreeTR->Fill();
	436	ifills++;
	437	}
	438	} else {
36e4cfbb	439	// Collect all hits
5fe09262	440	fTreeTR->GetEntry(itlast);
	441	for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
	442	for (Int_t ib = 0; ib < 7; ib++) {
	443	if (!trefs[ib]) continue;
	444	Int_t nh = trefs[ib]->GetEntries();
	445	for (Int_t ih = 0; ih < nh; ih++) {
	446	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	447	Int_t label = tr->Label();
	448	// Skip primaries
	449	if (label == ip) continue;
	450	if (label > dau2 \|\| label < dau1)
36e4cfbb	451	printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
36e4cfbb	452	itlast, label, dau1, dau2);
5fe09262	453	if (label == id) {
	454	// secondary found
	455	tr->SetDetectorId(ib-1);
	456	Int_t nref = fTrackReferences->GetEntriesFast();
	457	TClonesArray &lref = *fTrackReferences;
	458	new(lref[nref]) AliTrackReference(*tr);
	459	}
	460	} // hits
	461	} // branches
	462	fTmpTreeTR->Fill();
	463	fTrackReferences->Clear();
	464	ifills++;
	465	} // daughters
	466	} // has hits
	467	} // tracks
5fe09262	468	//
5fe09262	469	// Now loop again and write the primaries
36e4cfbb	470	//
5fe09262	471	it = nt - 1;
	472	for (Int_t ip = 0; ip < np; ip++) {
	473	Int_t labmax = -1;
	474	while (labmax < ip && it > -1) {
	475	fTreeTR->GetEntry(it--);
	476	for (Int_t ib = 0; ib < 7; ib++) {
	477	if (!trefs[ib]) continue;
	478	Int_t nh = trefs[ib]->GetEntries();
	479	//
	480	// Loop over reference hits and find primary labels
	481	for (Int_t ih = 0; ih < nh; ih++) {
	482	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	483	Int_t label = tr->Label();
	484	if (label < np && label > labmax) {
	485	labmax = label;
	486	}
	487
	488	if (label == ip) {
	489	tr->SetDetectorId(ib-1);
	490	Int_t nref = fTrackReferences->GetEntriesFast();
	491	TClonesArray &lref = *fTrackReferences;
	492	new(lref[nref]) AliTrackReference(*tr);
	493	}
	494	} // hits
	495	} // branches
	496	} // entries
	497	it++;
	498	fTmpTreeTR->Fill();
	499	fTrackReferences->Clear();
	500	ifills++;
	501	} // tracks
	502	// Check
5fe09262	503	if (ifills != fStack->GetNtrack())
36e4cfbb	504	printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
36e4cfbb	505	ifills, fStack->GetNtrack());
5fe09262	506	fTreeTR = fTmpTreeTR;
5fe09262	507	}