[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(0),
    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),
    fExternal(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),
    fExternal(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
    if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
    //
    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;
}

AliVParticle* AliMCEvent::GetTrack(Int_t i) const
{
    // Get MC Particle i
    //

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