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

Bool_t AliMCEventHandler::GetEvent(Int_t iev)
{
    // Load the event number iev
    Int_t inew  = iev/fEventsPerFile;
    if (inew != fFileNumber) {
	fFileNumber = inew;
	if (!OpenFile(fFileNumber)){
	    return kFALSE;
	}
    }
    
    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) {
	printf("AliMCEventHandler: Event #%5d not found\n", iev);
	return kFALSE;
    }
    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;

}

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