[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>


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(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),
    fStack(0),
    fHeader(new AliHeader()),
    fTrackReferences(new TClonesArray("AliTrackReference", 200)),
    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
    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();
	    delete dirTR;
	} 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()
{
    // Reset files
    if (fFileE)  delete fFileE;
    if (fFileK)  delete fFileK;
    if (fFileTR) delete fFileTR;
}

			    
Bool_t AliMCEventHandler::FinishEvent()
{
    // Reset the stack 
    Stack()->Reset();
    
    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("all");
    }
    
    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);

//
//
//    fTreeTR->SetBranchStatus("*",      0);
    fTreeTR->SetBranchStatus("ABSO",  0);
    fTreeTR->SetBranchStatus("BODY",  0);
    fTreeTR->SetBranchStatus("DIPO",  0);
    fTreeTR->SetBranchStatus("EMCAL", 0);
    fTreeTR->SetBranchStatus("FMD",   0);
    fTreeTR->SetBranchStatus("HALL",  0);
    fTreeTR->SetBranchStatus("MAG",   0);
    fTreeTR->SetBranchStatus("PHOS",  0);
    fTreeTR->SetBranchStatus("PIPE",  0);
    fTreeTR->SetBranchStatus("PMD",   0);
    fTreeTR->SetBranchStatus("RICH",  0);
    fTreeTR->SetBranchStatus("SHIL",  0);
    fTreeTR->SetBranchStatus("START", 0);
    fTreeTR->SetBranchStatus("VZERO", 0);
    fTreeTR->SetBranchStatus("ZDC",   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);

    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
    delete fTreeTR;
    for (Int_t ib = 0; ib < 7; ib++) {
	if (trefs[ib]) delete trefs[ib];
    }

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