[u/mrichter/AliRoot.git] / STEER / AliMCEvent.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, 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 for Kinematic Events
//     Author: Andreas Morsch, CERN
//-------------------------------------------------------------------------
#include <TArrow.h>
#include <TMarker.h>
#include <TH2F.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TRefArray.h>
#include <TList.h>

#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"


Int_t AliMCEvent::fgkBgLabelOffset(10000000);


AliMCEvent::AliMCEvent():
    AliVEvent(),
    fStack(0),
    fMCParticles(new TClonesArray("AliMCParticle",1000)),
    fMCParticleMap(0),
    fHeader(new AliHeader()),
    fTRBuffer(0),
    fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
    fTreeTR(0),
    fTmpTreeTR(0),
    fTmpFileTR(0),
    fNprimaries(-1),
    fNparticles(-1),
    fSubsidiaryEvents(0),
    fPrimaryOffset(0),
    fSecondaryOffset(0)
{
    // Default constructor
}

AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
    AliVEvent(mcEvnt),
    fStack(mcEvnt.fStack),
    fMCParticles(mcEvnt.fMCParticles),
    fMCParticleMap(mcEvnt.fMCParticleMap),
    fHeader(mcEvnt.fHeader),
    fTRBuffer(mcEvnt.fTRBuffer),
    fTrackReferences(mcEvnt.fTrackReferences),
    fTreeTR(mcEvnt.fTreeTR),
    fTmpTreeTR(mcEvnt.fTmpTreeTR),
    fTmpFileTR(mcEvnt.fTmpFileTR),
    fNprimaries(mcEvnt.fNprimaries),
    fNparticles(mcEvnt.fNparticles),
    fSubsidiaryEvents(0),
    fPrimaryOffset(0),
    fSecondaryOffset(0)
{ 
// Copy constructor
}


AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
{
    // assignment operator
    if (this!=&mcEvnt) { 
	AliVEvent::operator=(mcEvnt); 
    }
  
    return *this; 
}

void AliMCEvent::ConnectTreeE (TTree* tree)
{
    // Connect the event header tree
    tree->SetBranchAddress("Header", &fHeader);
}

void AliMCEvent::ConnectTreeK (TTree* tree)
{
    // Connect the kinematics tree to the stack
    fStack = fHeader->Stack();
    fStack->ConnectTree(tree);
    //
    // Load the event
    fStack->GetEvent();
    fNparticles = fStack->GetNtrack();
    fNprimaries = fStack->GetNprimary();
    Int_t iev  = fHeader->GetEvent();
    Int_t ievr = fHeader->GetEventNrInRun();
    
    AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n", 
		 iev, ievr, fNparticles, fNprimaries));
 
    // This is a cache for the TParticles converted to MCParticles on user request
    if (fMCParticleMap) {
	fMCParticleMap->Clear();
	fMCParticles->Delete();
	if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
    }
    else
	fMCParticleMap = new TRefArray(fNparticles);
}

void AliMCEvent::ConnectTreeTR (TTree* tree)
{
    // Connect the track reference tree
    fTreeTR = tree;
    
    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", &fTRBuffer);
    }
}

Int_t AliMCEvent::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 (fTreeTR) {
	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	trefs    = fTRBuffer;
	return trefs->GetEntries();
    } else {
	trefs = 0;
	return -1;
    }
}


void AliMCEvent::Clean()
{
    // Clean-up before new trees are connected
//    if (fHeader) {
//	delete fHeader;
//	fHeader = 0;
//    }

    delete fStack; fStack = 0;

    // Clear TR
    if (fTRBuffer) {
	fTRBuffer->Delete();
	delete fTRBuffer;
	fTRBuffer = 0;
    }
}

void AliMCEvent::FinishEvent()
{
  // Clean-up after event
  //    
    fStack->Reset(0);
    fMCParticles->Delete();
    fMCParticleMap->Clear();
    if (fTRBuffer)
      fTRBuffer->Delete();
    fTrackReferences->Delete();
    fNparticles = -1;
    fNprimaries = -1;    
    fStack      =  0;
//    fSubsidiaryEvents->Clear();
    fSubsidiaryEvents = 0;
}


void AliMCEvent::DrawCheck(Int_t i, Int_t search)
{
    //
    // Simple event display for debugging
    if (!fTreeTR) {
	AliWarning("No Track Reference information available");
	return;
    } 
    
    if (i > -1 && i < fNparticles) {
	fTreeTR->GetEntry(fStack->TreeKEntry(i));
    } else {
	AliWarning("AliMCEvent::GetEntry: Index out of range");
    }
    
    Int_t nh = fTRBuffer->GetEntries();
    
    
    if (search) {
	while(nh <= search && i < fNparticles - 1) {
	    i++;
	    fTreeTR->GetEntry(fStack->TreeKEntry(i));
	    nh =  fTRBuffer->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*) fTRBuffer->At(ih);
	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	m->Draw();
	m->SetMarkerSize(0.4);
	
    }
}


void AliMCEvent::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 (!fTRBuffer)  fTRBuffer = new TClonesArray("AliTrackReference", 100);
    fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 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 =  fTRBuffer->GetEntriesFast();
			    TClonesArray &lref = *fTRBuffer;
			    new(lref[nref]) AliTrackReference(*tr);
			}
		    } // hits
		} // branches
		fTmpTreeTR->Fill();
		fTRBuffer->Delete();
		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 = fTRBuffer->GetEntriesFast();
			TClonesArray &lref = *fTRBuffer;
			new(lref[nref]) AliTrackReference(*tr);
		    }
		} // hits
	    } // branches
	} // entries
	it++;
	fTmpTreeTR->Fill();
	fTRBuffer->Delete();
	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("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", 
	       ifills, fStack->GetNtrack());

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

AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
{
    // Get MC Particle i
    //
    //
    // Check first if this explicitely accesses the subsidiary event
    if (i >= BgLabelOffset()) {
	if (fSubsidiaryEvents) {
	    AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	    return (bgEvent->GetTrack(i - BgLabelOffset()));
	} else {
	    return 0;
	}
    }
    
    //
    AliMCParticle *mcParticle = 0;
    TParticle     *particle   = 0;
    TClonesArray  *trefs      = 0;
    Int_t          ntref      = 0;
    TRefArray     *rarray     = 0;


    // Out of range check
    if (i < 0 || i >= fNparticles) {
	AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
	mcParticle = 0;
	return (mcParticle);
    }

    
    if (fSubsidiaryEvents) {
	AliMCEvent*   mc;
	Int_t idx = FindIndexAndEvent(i, mc);
	return (mc->GetTrack(idx));
    } 

    //
    // First check If the MC Particle has been already cached
    if(!fMCParticleMap->At(i)) {
	// Get particle from the stack
	particle   = fStack->Particle(i);
	// Get track references from Tree TR
	if (fTreeTR) {
	    fTreeTR->GetEntry(fStack->TreeKEntry(i));
	    trefs     = fTRBuffer;
	    ntref     = trefs->GetEntriesFast();
	    rarray    = new TRefArray(ntref);
	    Int_t nen = fTrackReferences->GetEntriesFast();
	    for (Int_t j = 0; j < ntref; j++) {
		// Save the track references in a TClonesArray
		AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
		// Save the pointer in a TRefArray
		new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
		rarray->AddAt((*fTrackReferences)[nen], j);
		nen++;
	    } // loop over track references for entry i
	} // if TreeTR available
	Int_t nentries = fMCParticles->GetEntriesFast();
	new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
	mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
	fMCParticleMap->AddAt(mcParticle, i);

	TParticle* part = mcParticle->Particle();
	Int_t imo  = part->GetFirstMother();
	Int_t id1  = part->GetFirstDaughter();
	Int_t id2  = part->GetLastDaughter();
	if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
	    // Remapping of the mother and daughter indices
	    if (imo < fNprimaries) {
		mcParticle->SetMother(imo + fPrimaryOffset);
	    } else {
		mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
	    }

	    if (id1 < fNprimaries) {
		mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
		mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
	    } else {
		mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
		mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
	    }

	    
	    if (i > fNprimaries) {
		mcParticle->SetLabel(i + fPrimaryOffset);
	    } else {
		mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
	    }
	    
	    
	} else {
	    mcParticle->SetFirstDaughter(id1);
	    mcParticle->SetLastDaughter (id2);
	    mcParticle->SetMother       (imo);
	}

    } else {
	mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
    }
    
    
    return mcParticle;
}

    AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}


void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event) 
{
    // Add a subsidiary event to the list; for example merged background event.
    if (!fSubsidiaryEvents) {
	TList* events = new TList();
	events->Add(new AliMCEvent(*this));
	fSubsidiaryEvents = events;
    }
    
    fSubsidiaryEvents->Add(event);
}

Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
{
    // Find the index and event in case of composed events like signal + background
    TIter next(fSubsidiaryEvents);
    next.Reset();
     if (oldidx < fNprimaries) {
	while((event = (AliMCEvent*)next())) {
	    if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
	}
	return (oldidx - event->GetPrimaryOffset());
    } else {
	while((event = (AliMCEvent*)next())) {
	    if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
	}
	return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
    }
}

Int_t AliMCEvent::BgLabelToIndex(Int_t label)
{
    // Convert a background label to an absolute index
    if (fSubsidiaryEvents) {
	AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	label -= BgLabelOffset();
	if (label < bgEvent->GetNumberOfPrimaries()) {
	    label += bgEvent->GetPrimaryOffset();
	} else {
	    label += (bgEvent->GetSecondaryOffset() - fNprimaries);
	}
    }
    return (label);
}


Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) 
{
//
// Delegate to subevent if necesarry 

    
    if (!fSubsidiaryEvents) {
	return fStack->IsPhysicalPrimary(i);
    } else {
	AliMCEvent* evt = 0;
	Int_t idx = FindIndexAndEvent(i, evt);
	return (evt->IsPhysicalPrimary(idx));
    }
}

void AliMCEvent::InitEvent()
{
    if (fSubsidiaryEvents) {
	TIter next(fSubsidiaryEvents);
	AliMCEvent* evt;
	fNprimaries = 0;
	fNparticles = 0;
	
	while((evt = (AliMCEvent*)next())) {
	    fNprimaries += evt->GetNumberOfPrimaries();	
	    fNparticles += evt->GetNumberOfTracks();    
	}
	
	Int_t ioffp = 0;
	Int_t ioffs = fNprimaries;
	next.Reset();
	
	while((evt = (AliMCEvent*)next())) {
	    evt->SetPrimaryOffset(ioffp);
	    evt->SetSecondaryOffset(ioffs);
	    ioffp += evt->GetNumberOfPrimaries();
	    ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());	    
	}
    }
}


ClassImp(AliMCEvent)
Commit	Line	Data
415d9f5c	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, 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	//-------------------------------------------------------------------------
	19	// Class for Kinematic Events
	20	// Author: Andreas Morsch, CERN
	21	//-------------------------------------------------------------------------
	22	#include <TArrow.h>
	23	#include <TMarker.h>
	24	#include <TH2F.h>
	25	#include <TTree.h>
	26	#include <TFile.h>
	27	#include <TParticle.h>
	28	#include <TClonesArray.h>
93df0e9b	29	#include <TRefArray.h>
93836e1b	30	#include <TList.h>
415d9f5c	31
	32	#include "AliLog.h"
	33	#include "AliMCEvent.h"
	34	#include "AliStack.h"
	35	#include "AliTrackReference.h"
	36	#include "AliHeader.h"
83626b6b	37	#include "AliGenEventHeader.h"
415d9f5c	38
415d9f5c	39
93836e1b	40	Int_t AliMCEvent::fgkBgLabelOffset(10000000);
	41
	42
415d9f5c	43	AliMCEvent::AliMCEvent():
	44	AliVEvent(),
	45	fStack(0),
	46	fMCParticles(new TClonesArray("AliMCParticle",1000)),
	47	fMCParticleMap(0),
fb9aeda5	48	fHeader(new AliHeader()),
93df0e9b	49	fTRBuffer(0),
93df0e9b	50	fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
415d9f5c	51	fTreeTR(0),
	52	fTmpTreeTR(0),
	53	fTmpFileTR(0),
	54	fNprimaries(-1),
93836e1b	55	fNparticles(-1),
	56	fSubsidiaryEvents(0),
	57	fPrimaryOffset(0),
	58	fSecondaryOffset(0)
415d9f5c	59	{
	60	// Default constructor
	61	}
	62
	63	AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
	64	AliVEvent(mcEvnt),
93836e1b	65	fStack(mcEvnt.fStack),
	66	fMCParticles(mcEvnt.fMCParticles),
	67	fMCParticleMap(mcEvnt.fMCParticleMap),
	68	fHeader(mcEvnt.fHeader),
	69	fTRBuffer(mcEvnt.fTRBuffer),
	70	fTrackReferences(mcEvnt.fTrackReferences),
	71	fTreeTR(mcEvnt.fTreeTR),
	72	fTmpTreeTR(mcEvnt.fTmpTreeTR),
	73	fTmpFileTR(mcEvnt.fTmpFileTR),
	74	fNprimaries(mcEvnt.fNprimaries),
	75	fNparticles(mcEvnt.fNparticles),
	76	fSubsidiaryEvents(0),
	77	fPrimaryOffset(0),
	78	fSecondaryOffset(0)
415d9f5c	79	{
	80	// Copy constructor
	81	}
	82
	83
	84	AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
	85	{
	86	// assignment operator
	87	if (this!=&mcEvnt) {
	88	AliVEvent::operator=(mcEvnt);
	89	}
	90
	91	return *this;
	92	}
	93
	94	void AliMCEvent::ConnectTreeE (TTree* tree)
	95	{
	96	// Connect the event header tree
	97	tree->SetBranchAddress("Header", &fHeader);
	98	}
	99
	100	void AliMCEvent::ConnectTreeK (TTree* tree)
	101	{
	102	// Connect the kinematics tree to the stack
	103	fStack = fHeader->Stack();
	104	fStack->ConnectTree(tree);
	105	//
	106	// Load the event
	107	fStack->GetEvent();
	108	fNparticles = fStack->GetNtrack();
	109	fNprimaries = fStack->GetNprimary();
ed97dc98	110	Int_t iev = fHeader->GetEvent();
	111	Int_t ievr = fHeader->GetEventNrInRun();
	112
	113	AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
	114	iev, ievr, fNparticles, fNprimaries));
415d9f5c	115
	116	// This is a cache for the TParticles converted to MCParticles on user request
	117	if (fMCParticleMap) {
678a0069	118	fMCParticleMap->Clear();
ee2774b7	119	fMCParticles->Delete();
678a0069	120	if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
678a0069	121	}
415d9f5c	122	else
93df0e9b	123	fMCParticleMap = new TRefArray(fNparticles);
415d9f5c	124	}
	125
	126	void AliMCEvent::ConnectTreeTR (TTree* tree)
	127	{
	128	// Connect the track reference tree
	129	fTreeTR = tree;
	130
	131	if (fTreeTR->GetBranch("AliRun")) {
	132	if (fTmpFileTR) {
	133	fTmpFileTR->Close();
	134	delete fTmpFileTR;
	135	}
	136	// This is an old format with one branch per detector not in synch with TreeK
	137	ReorderAndExpandTreeTR();
	138	} else {
	139	// New format
93df0e9b	140	fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
415d9f5c	141	}
	142	}
	143
	144	Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle& particle, TClonesArray& trefs)
	145	{
	146	// Retrieve entry i
	147	if (i < 0 \|\| i >= fNparticles) {
	148	AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
	149	particle = 0;
	150	trefs = 0;
	151	return (-1);
	152	}
	153	particle = fStack->Particle(i);
	154	if (fTreeTR) {
	155	fTreeTR->GetEntry(fStack->TreeKEntry(i));
93df0e9b	156	trefs = fTRBuffer;
415d9f5c	157	return trefs->GetEntries();
	158	} else {
	159	trefs = 0;
	160	return -1;
	161	}
	162	}
	163
	164
	165	void AliMCEvent::Clean()
	166	{
	167	// Clean-up before new trees are connected
fb9aeda5	168	// if (fHeader) {
	169	// delete fHeader;
	170	// fHeader = 0;
	171	// }
c6360f33	172
c6360f33	173	delete fStack; fStack = 0;
415d9f5c	174
415d9f5c	175	// Clear TR
93df0e9b	176	if (fTRBuffer) {
ee2774b7	177	fTRBuffer->Delete();
93df0e9b	178	delete fTRBuffer;
93df0e9b	179	fTRBuffer = 0;
415d9f5c	180	}
	181	}
	182
	183	void AliMCEvent::FinishEvent()
	184	{
f1ac8a97	185	// Clean-up after event
c6360f33	186	//
	187	fStack->Reset(0);
	188	fMCParticles->Delete();
678a0069	189	fMCParticleMap->Clear();
568d5aed	190	if (fTRBuffer)
568d5aed	191	fTRBuffer->Delete();
c6360f33	192	fTrackReferences->Delete();
2bb794ba	193	fNparticles = -1;
	194	fNprimaries = -1;
	195	fStack = 0;
93836e1b	196	// fSubsidiaryEvents->Clear();
93836e1b	197	fSubsidiaryEvents = 0;
415d9f5c	198	}
	199
	200
93df0e9b	201
415d9f5c	202	void AliMCEvent::DrawCheck(Int_t i, Int_t search)
	203	{
	204	//
	205	// Simple event display for debugging
	206	if (!fTreeTR) {
	207	AliWarning("No Track Reference information available");
	208	return;
	209	}
	210
	211	if (i > -1 && i < fNparticles) {
	212	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	213	} else {
	214	AliWarning("AliMCEvent::GetEntry: Index out of range");
	215	}
	216
93df0e9b	217	Int_t nh = fTRBuffer->GetEntries();
415d9f5c	218
	219
	220	if (search) {
	221	while(nh <= search && i < fNparticles - 1) {
	222	i++;
	223	fTreeTR->GetEntry(fStack->TreeKEntry(i));
93df0e9b	224	nh = fTRBuffer->GetEntries();
415d9f5c	225	}
	226	printf("Found Hits at %5d\n", i);
	227	}
	228	TParticle* particle = fStack->Particle(i);
	229
	230	TH2F* h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
	231	Float_t x0 = particle->Vx();
	232	Float_t y0 = particle->Vy();
	233
	234	Float_t x1 = particle->Vx() + particle->Px() * 50.;
	235	Float_t y1 = particle->Vy() + particle->Py() * 50.;
	236
	237	TArrow* a = new TArrow(x0, y0, x1, y1, 0.01);
	238	h->Draw();
	239	a->SetLineColor(2);
	240
	241	a->Draw();
	242
	243	for (Int_t ih = 0; ih < nh; ih++) {
93df0e9b	244	AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
415d9f5c	245	TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
	246	m->Draw();
	247	m->SetMarkerSize(0.4);
	248
	249	}
	250	}
	251
	252
	253	void AliMCEvent::ReorderAndExpandTreeTR()
	254	{
	255	//
	256	// Reorder and expand the track reference tree in order to match the kinematics tree.
	257	// Copy the information from different branches into one
	258	//
	259	// TreeTR
	260
	261	fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
	262	fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
93df0e9b	263	if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
93df0e9b	264	fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
415d9f5c	265
	266
	267	//
	268	// Activate the used branches only. Otherwisw we get a bad memory leak.
	269	fTreeTR->SetBranchStatus("*", 0);
	270	fTreeTR->SetBranchStatus("AliRun.*", 1);
	271	fTreeTR->SetBranchStatus("ITS.*", 1);
	272	fTreeTR->SetBranchStatus("TPC.*", 1);
	273	fTreeTR->SetBranchStatus("TRD.*", 1);
	274	fTreeTR->SetBranchStatus("TOF.*", 1);
	275	fTreeTR->SetBranchStatus("FRAME.*", 1);
	276	fTreeTR->SetBranchStatus("MUON.*", 1);
	277	//
	278	// Connect the active branches
	279	TClonesArray* trefs[7];
	280	for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
	281	if (fTreeTR){
	282	// make branch for central track references
	283	if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
	284	if (fTreeTR->GetBranch("ITS")) fTreeTR->SetBranchAddress("ITS", &trefs[1]);
	285	if (fTreeTR->GetBranch("TPC")) fTreeTR->SetBranchAddress("TPC", &trefs[2]);
	286	if (fTreeTR->GetBranch("TRD")) fTreeTR->SetBranchAddress("TRD", &trefs[3]);
	287	if (fTreeTR->GetBranch("TOF")) fTreeTR->SetBranchAddress("TOF", &trefs[4]);
	288	if (fTreeTR->GetBranch("FRAME")) fTreeTR->SetBranchAddress("FRAME", &trefs[5]);
	289	if (fTreeTR->GetBranch("MUON")) fTreeTR->SetBranchAddress("MUON", &trefs[6]);
	290	}
	291
	292	Int_t np = fStack->GetNprimary();
	293	Int_t nt = fTreeTR->GetEntries();
	294
	295	//
	296	// Loop over tracks and find the secondaries with the help of the kine tree
	297	Int_t ifills = 0;
	298	Int_t it = 0;
	299	Int_t itlast = 0;
	300	TParticle* part;
	301
	302	for (Int_t ip = np - 1; ip > -1; ip--) {
	303	part = fStack->Particle(ip);
	304	// printf("Particle %5d %5d %5d %5d %5d %5d \n",
	305	// ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
	306	// part->GetLastDaughter(), part->TestBit(kTransportBit));
	307
	308	// Determine range of secondaries produced by this primary during transport
	309	Int_t dau1 = part->GetFirstDaughter();
	310	if (dau1 < np) continue; // This particle has no secondaries produced during transport
	311	Int_t dau2 = -1;
	312	if (dau1 > -1) {
	313	Int_t inext = ip - 1;
	314	while (dau2 < 0) {
	315	if (inext >= 0) {
	316	part = fStack->Particle(inext);
	317	dau2 = part->GetFirstDaughter();
	318	if (dau2 == -1 \|\| dau2 < np) {
	319	dau2 = -1;
	320	} else {
	321	dau2--;
	322	}
	323	} else {
	324	dau2 = fStack->GetNtrack() - 1;
	325	}
	326	inext--;
	327	} // find upper bound
	328	} // dau2 < 0
329
330
331	// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
332	//
333	// Loop over reference hits and find secondary label
334	// First the tricky part: find the entry in treeTR than contains the hits or
335	// make sure that no hits exist.
336	//
337	Bool_t hasHits = kFALSE;
338	Bool_t isOutside = kFALSE;
339
340	it = itlast;
341	while (!hasHits && !isOutside && it < nt) {
342	fTreeTR->GetEntry(it++);
343	for (Int_t ib = 0; ib < 7; ib++) {
344	if (!trefs[ib]) continue;
345	Int_t nh = trefs[ib]->GetEntries();
346	for (Int_t ih = 0; ih < nh; ih++) {
347	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
348	Int_t label = tr->Label();
349	if (label >= dau1 && label <= dau2) {
350	hasHits = kTRUE;
351	itlast = it - 1;
352	break;
353	}
354	if (label > dau2 \|\| label < ip) {
355	isOutside = kTRUE;
356	itlast = it - 1;
357	break;
358	}
359	} // hits
360	if (hasHits \|\| isOutside) break;
361	} // branches
362	} // entries
363
364	if (!hasHits) {
365	// Write empty entries
366	for (Int_t id = dau1; (id <= dau2); id++) {
367	fTmpTreeTR->Fill();
368	ifills++;
369	}
370	} else {
371	// Collect all hits
372	fTreeTR->GetEntry(itlast);
373	for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
374	for (Int_t ib = 0; ib < 7; ib++) {
375	if (!trefs[ib]) continue;
376	Int_t nh = trefs[ib]->GetEntries();
377	for (Int_t ih = 0; ih < nh; ih++) {
378	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
379	Int_t label = tr->Label();
380	// Skip primaries
381	if (label == ip) continue;
382	if (label > dau2 \|\| label < dau1)
383	printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
384	itlast, label, dau1, dau2);
385	if (label == id) {
386	// secondary found
387	tr->SetDetectorId(ib-1);
93df0e9b	388	Int_t nref = fTRBuffer->GetEntriesFast();
93df0e9b	389	TClonesArray &lref = *fTRBuffer;
415d9f5c	390	new(lref[nref]) AliTrackReference(*tr);
	391	}
	392	} // hits
	393	} // branches
	394	fTmpTreeTR->Fill();
ee2774b7	395	fTRBuffer->Delete();
415d9f5c	396	ifills++;
	397	} // daughters
	398	} // has hits
	399	} // tracks
	400
	401	//
	402	// Now loop again and write the primaries
	403	//
	404	it = nt - 1;
	405	for (Int_t ip = 0; ip < np; ip++) {
	406	Int_t labmax = -1;
	407	while (labmax < ip && it > -1) {
	408	fTreeTR->GetEntry(it--);
	409	for (Int_t ib = 0; ib < 7; ib++) {
	410	if (!trefs[ib]) continue;
	411	Int_t nh = trefs[ib]->GetEntries();
	412	//
	413	// Loop over reference hits and find primary labels
	414	for (Int_t ih = 0; ih < nh; ih++) {
	415	AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
	416	Int_t label = tr->Label();
	417	if (label < np && label > labmax) {
	418	labmax = label;
	419	}
	420
	421	if (label == ip) {
	422	tr->SetDetectorId(ib-1);
93df0e9b	423	Int_t nref = fTRBuffer->GetEntriesFast();
93df0e9b	424	TClonesArray &lref = *fTRBuffer;
415d9f5c	425	new(lref[nref]) AliTrackReference(*tr);
	426	}
	427	} // hits
	428	} // branches
	429	} // entries
	430	it++;
	431	fTmpTreeTR->Fill();
ee2774b7	432	fTRBuffer->Delete();
415d9f5c	433	ifills++;
	434	} // tracks
	435	// Check
	436
	437
	438	// Clean-up
	439	delete fTreeTR; fTreeTR = 0;
	440
	441	for (Int_t ib = 0; ib < 7; ib++) {
	442	if (trefs[ib]) {
	443	trefs[ib]->Clear();
	444	delete trefs[ib];
	445	trefs[ib] = 0;
	446	}
	447	}
	448
	449	if (ifills != fStack->GetNtrack())
	450	printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
	451	ifills, fStack->GetNtrack());
	452
	453	fTmpTreeTR->Write();
	454	fTreeTR = fTmpTreeTR;
	455	}
	456
	457	AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
	458	{
	459	// Get MC Particle i
93836e1b	460	//
	461	//
	462	// Check first if this explicitely accesses the subsidiary event
	463	if (i >= BgLabelOffset()) {
	464	if (fSubsidiaryEvents) {
	465	AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	466	return (bgEvent->GetTrack(i - BgLabelOffset()));
	467	} else {
	468	return 0;
	469	}
	470	}
	471
	472	//
93df0e9b	473	AliMCParticle *mcParticle = 0;
	474	TParticle *particle = 0;
	475	TClonesArray *trefs = 0;
	476	Int_t ntref = 0;
	477	TRefArray *rarray = 0;
93836e1b	478
	479
	480
415d9f5c	481	// Out of range check
	482	if (i < 0 \|\| i >= fNparticles) {
	483	AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
	484	mcParticle = 0;
	485	return (mcParticle);
	486	}
93836e1b	487
415d9f5c	488
93836e1b	489	if (fSubsidiaryEvents) {
	490	AliMCEvent* mc;
	491	Int_t idx = FindIndexAndEvent(i, mc);
	492	return (mc->GetTrack(idx));
	493	}
415d9f5c	494
415d9f5c	495	//
93836e1b	496	// First check If the MC Particle has been already cached
415d9f5c	497	if(!fMCParticleMap->At(i)) {
93df0e9b	498	// Get particle from the stack
	499	particle = fStack->Particle(i);
	500	// Get track references from Tree TR
	501	if (fTreeTR) {
	502	fTreeTR->GetEntry(fStack->TreeKEntry(i));
	503	trefs = fTRBuffer;
	504	ntref = trefs->GetEntriesFast();
	505	rarray = new TRefArray(ntref);
	506	Int_t nen = fTrackReferences->GetEntriesFast();
	507	for (Int_t j = 0; j < ntref; j++) {
	508	// Save the track references in a TClonesArray
	509	AliTrackReference* ref = dynamic_cast<AliTrackReference>((fTRBuffer)[j]);
	510	// Save the pointer in a TRefArray
	511	new ((fTrackReferences)[nen]) AliTrackReference(ref);
	512	rarray->AddAt((*fTrackReferences)[nen], j);
	513	nen++;
	514	} // loop over track references for entry i
	515	} // if TreeTR available
415d9f5c	516	Int_t nentries = fMCParticles->GetEntriesFast();
eee13e8d	517	new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
415d9f5c	518	mcParticle = dynamic_cast<AliMCParticle>((fMCParticles)[nentries]);
415d9f5c	519	fMCParticleMap->AddAt(mcParticle, i);
93836e1b	520
	521	TParticle* part = mcParticle->Particle();
	522	Int_t imo = part->GetFirstMother();
	523	Int_t id1 = part->GetFirstDaughter();
	524	Int_t id2 = part->GetLastDaughter();
	525	if (fPrimaryOffset > 0 \|\| fSecondaryOffset > 0) {
	526	// Remapping of the mother and daughter indices
	527	if (imo < fNprimaries) {
	528	mcParticle->SetMother(imo + fPrimaryOffset);
	529	} else {
	530	mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
	531	}
	532
	533	if (id1 < fNprimaries) {
	534	mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
	535	mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
	536	} else {
	537	mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
	538	mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
	539	}
	540
	541
	542	if (i > fNprimaries) {
	543	mcParticle->SetLabel(i + fPrimaryOffset);
	544	} else {
	545	mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
	546	}
	547
	548
	549	} else {
	550	mcParticle->SetFirstDaughter(id1);
	551	mcParticle->SetLastDaughter (id2);
	552	mcParticle->SetMother (imo);
	553	}
	554
415d9f5c	555	} else {
	556	mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
	557	}
93836e1b	558
93836e1b	559
415d9f5c	560	return mcParticle;
	561	}
	562
93836e1b	563	AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
	564
	565
	566	void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
	567	{
	568	// Add a subsidiary event to the list; for example merged background event.
	569	if (!fSubsidiaryEvents) {
	570	TList* events = new TList();
	571	events->Add(new AliMCEvent(*this));
	572	fSubsidiaryEvents = events;
	573	}
	574
	575	fSubsidiaryEvents->Add(event);
	576	}
	577
	578	Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
	579	{
	580	// Find the index and event in case of composed events like signal + background
	581	TIter next(fSubsidiaryEvents);
	582	next.Reset();
	583	if (oldidx < fNprimaries) {
	584	while((event = (AliMCEvent*)next())) {
	585	if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
	586	}
	587	return (oldidx - event->GetPrimaryOffset());
	588	} else {
	589	while((event = (AliMCEvent*)next())) {
	590	if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
	591	}
	592	return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
	593	}
	594	}
	595
	596	Int_t AliMCEvent::BgLabelToIndex(Int_t label)
	597	{
	598	// Convert a background label to an absolute index
	599	if (fSubsidiaryEvents) {
	600	AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
	601	label -= BgLabelOffset();
	602	if (label < bgEvent->GetNumberOfPrimaries()) {
	603	label += bgEvent->GetPrimaryOffset();
	604	} else {
	605	label += (bgEvent->GetSecondaryOffset() - fNprimaries);
	606	}
	607	}
	608	return (label);
	609	}
	610
	611
	612	Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i)
	613	{
	614	//
	615	// Delegate to subevent if necesarry
	616
	617
	618	if (!fSubsidiaryEvents) {
	619	return fStack->IsPhysicalPrimary(i);
	620	} else {
	621	AliMCEvent* evt = 0;
	622	Int_t idx = FindIndexAndEvent(i, evt);
	623	return (evt->IsPhysicalPrimary(idx));
	624	}
	625	}
	626
627	void AliMCEvent::InitEvent()
628	{
629	if (fSubsidiaryEvents) {
630	TIter next(fSubsidiaryEvents);
631	AliMCEvent* evt;
632	fNprimaries = 0;
633	fNparticles = 0;
634
635	while((evt = (AliMCEvent*)next())) {
636	fNprimaries += evt->GetNumberOfPrimaries();
637	fNparticles += evt->GetNumberOfTracks();
638	}
639
640	Int_t ioffp = 0;
641	Int_t ioffs = fNprimaries;
642	next.Reset();
643
644	while((evt = (AliMCEvent*)next())) {
645	evt->SetPrimaryOffset(ioffp);
646	evt->SetSecondaryOffset(ioffs);
647	ioffp += evt->GetNumberOfPrimaries();
648	ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());
649	}
650	}
651	}
652
415d9f5c	653
415d9f5c	654	ClassImp(AliMCEvent)