[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 <TString.h>
#include <TClonesArray.h>
#include <TDirectoryFile.h>
#include <TArrow.h>
#include <TMarker.h>
#include <TH2F.h>
#include <TDirectoryFile.h>


ClassImp(AliMCEventHandler)

AliMCEventHandler::AliMCEventHandler() :
    AliVEventHandler(),
    fFileE(0),
    fFileK(0),
    fFileTR(0),
    fTmpFileTR(0),
    fTreeE(0),
    fTreeK(0),
    fTreeTR(0),
    fTmpTreeTR(0),
    fDirK(0),
    fDirTR(0),
    fStack(0),
    fHeader(0),
    fTrackReferences(0),
    fNEvent(-1),
    fEvent(-1),
    fNprimaries(-1),
    fNparticles(-1),
    fPathName(new TString("./")),
    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),
    fDirK(0),
    fDirTR(0),
    fStack(0),
    fHeader(0),
    fTrackReferences(0),
    fNEvent(-1),
    fEvent(-1),
    fNprimaries(-1),
    fNparticles(-1),
    fPathName(new TString("./")),
    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->Data()));
    if (!fFileE) AliFatal(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName->Data()));

    fFileE->GetObject("TE", fTreeE);
    fTreeE->SetBranchAddress("Header", &fHeader);
    fNEvent = fTreeE->GetEntries();
    //
    // Tree K
    fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), 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->Data(), fExtension));
    if (!fFileTR) AliWarning(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName->Data()));
    //
    // 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
    fFileK->GetObject(folder, fDirK);
    if (!fDirK) {
	AliWarning(Form("AliMCEventHandler: Event #%5d not found\n", iev));
	return kFALSE;
    }
    fDirK ->GetObject("TreeK", fTreeK);
    fStack->ConnectTree(fTreeK);
    fStack->GetEvent();
    //Tree TR 
    if (fFileTR) {
	// Check which format has been read
	fFileTR->GetObject(folder, fDirTR);
	fDirTR->GetObject("TreeTR", fTreeTR);
	if (fTreeTR->GetBranch("AliRun")) {
	    if (fTmpFileTR) {
		fTmpFileTR->Close();
		delete fTmpFileTR;
	    }
	    // This is an old format with one branch per detector not in synch with TreeK
	    ReorderAndExpandTreeTR();
	} else {
	    // New format 
	    fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
	}
    }
    //
    fNparticles = fStack->GetNtrack();
    fNprimaries = fStack->GetNprimary();
    AliInfo(Form("AliMCEventHandler: Event#: %5d Number of particles: %5d (all) %5d (primaries)\n", 
		 fEvent, fNparticles, fNprimaries));
    
    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 = TFile::Open(Form("%sKinematics%s.root", fPathName->Data(), fExtension));
    if (!fFileK) {
	AliFatal(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
	ok = kFALSE;
    }
    
    delete fFileTR;
    fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName->Data(), fExtension));
    if (!fFileTR) {
	AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName->Data()));
	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 < 0 || 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, Int_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 <= search && 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);
    delete fPathName;
    fPathName = new TString(path);
    ResetIO();
    InitIO("");
    return kTRUE;
}
    
void AliMCEventHandler::ResetIO()
{
//  Clear header and stack

    if (fHeader) {
	delete fHeader;
	fHeader = 0;
    }

    delete fStack;
    delete fTreeE; fTreeE = 0;
    
// Clear TR
    
    if (fTrackReferences) {
	fTrackReferences->Clear();
	delete fTrackReferences;
	fTrackReferences = 0;
    }
    
// Reset files
    if (fFileE)  delete fFileE;
    if (fFileK)  delete fFileK;
    if (fFileTR) delete fFileTR;
}

			    
Bool_t AliMCEventHandler::FinishEvent()
{
    // Clean-up after each event
    delete fDirTR;  fDirTR = 0;
    delete fDirK;   fDirK  = 0;    
    Stack()->Reset(0);
    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

    fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
    fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
    if (!fTrackReferences)  fTrackReferences = new TClonesArray("AliTrackReference", 100);
    fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 32000, 0);
    

//
//  Activate the used branches only. Otherwisw we get a bad memory leak.
    fTreeTR->SetBranchStatus("*",        0);
    fTreeTR->SetBranchStatus("AliRun.*", 1);
    fTreeTR->SetBranchStatus("ITS.*",    1);
    fTreeTR->SetBranchStatus("TPC.*",    1);
    fTreeTR->SetBranchStatus("TRD.*",    1);
    fTreeTR->SetBranchStatus("TOF.*",    1);
    fTreeTR->SetBranchStatus("FRAME.*",  1);
    fTreeTR->SetBranchStatus("MUON.*",   1);
//
//  Connect the active branches
    TClonesArray* trefs[7];
    for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
    if (fTreeTR){
	// make branch for central track references
	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
	if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
	if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
	if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
	if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
	if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
    }

    Int_t np = fStack->GetNprimary();
    Int_t nt = fTreeTR->GetEntries();
    
    //
    // Loop over tracks and find the secondaries with the help of the kine tree
    Int_t ifills = 0;
    Int_t it     = 0;
    Int_t itlast = 0;
    TParticle* part;

    for (Int_t ip = np - 1; ip > -1; ip--) {
	part = fStack->Particle(ip);
//	printf("Particle %5d %5d %5d %5d %5d %5d \n", 
//	       ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), 
//	       part->GetLastDaughter(), part->TestBit(kTransportBit));

	// Determine range of secondaries produced by this primary during transport	
	Int_t dau1  = part->GetFirstDaughter();
	if (dau1 < np) continue;  // This particle has no secondaries produced during transport
	Int_t dau2  = -1;
	if (dau1 > -1) {
	    Int_t inext = ip - 1;
	    while (dau2 < 0) {
		if (inext >= 0) {
		    part = fStack->Particle(inext);
		    dau2 =  part->GetFirstDaughter();
		    if (dau2 == -1 || dau2 < np) {
			dau2 = -1;
		    } else {
			dau2--;
		    }
		} else {
		    dau2 = fStack->GetNtrack() - 1;
		}
		inext--;
	    } // find upper bound
	}  // dau2 < 0
	

//	printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
//
// Loop over reference hits and find secondary label
// First the tricky part: find the entry in treeTR than contains the hits or
// make sure that no hits exist.
//
	Bool_t hasHits   = kFALSE;
	Bool_t isOutside = kFALSE;

	it = itlast;
	while (!hasHits && !isOutside && it < nt) {
	    fTreeTR->GetEntry(it++);
	    for (Int_t ib = 0; ib < 7; ib++) {
		if (!trefs[ib]) continue;
		Int_t nh = trefs[ib]->GetEntries();
		for (Int_t ih = 0; ih < nh; ih++) {
		    AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
		    Int_t label = tr->Label();
		    if (label >= dau1 && label <= dau2) {
			hasHits = kTRUE;
			itlast = it - 1;
			break;
		    }
		    if (label > dau2 || label < ip) {
			isOutside = kTRUE;
			itlast = it - 1;
			break;
		    }
		} // hits
		if (hasHits || isOutside) break;
	    } // branches
	} // entries

	if (!hasHits) {
	    // Write empty entries
	    for (Int_t id = dau1; (id <= dau2); id++) {
		fTmpTreeTR->Fill();
		ifills++;
	    } 
	} else {
	    // Collect all hits
	    fTreeTR->GetEntry(itlast);
	    for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
		for (Int_t ib = 0; ib < 7; ib++) {
		    if (!trefs[ib]) continue;
		    Int_t nh = trefs[ib]->GetEntries();
		    for (Int_t ih = 0; ih < nh; ih++) {
			AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
			Int_t label = tr->Label();
			// Skip primaries
			if (label == ip) continue;
			if (label > dau2 || label < dau1) 
			    printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n", 
				   itlast, label, dau1, dau2);
			if (label == id) {
			    // secondary found
			    tr->SetDetectorId(ib-1);
			    Int_t nref =  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


    // Clean-up
    delete fTreeTR; fTreeTR = 0;
    
    for (Int_t ib = 0; ib < 7; ib++) {
	if (trefs[ib]) {
	    trefs[ib]->Clear();
	    delete trefs[ib];
	    trefs[ib] = 0;
	}
    }

    if (ifills != fStack->GetNtrack()) 
	printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", 
	       ifills, fStack->GetNtrack());

    fTmpTreeTR->Write();
    fTreeTR = fTmpTreeTR;
}

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