#include <TArrow.h>
#include <TMarker.h>
#include <TH2F.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <TList.h>
#include <TArrayF.h>
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCVertex.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenCocktailEventHeader.h"
Int_t AliMCEvent::fgkBgLabelOffset(10000000);
AliMCEvent::AliMCEvent():
AliVEvent(),
fStack(0),
fMCParticles(0),
fMCParticleMap(0),
fHeader(new AliHeader()),
fAODMCHeader(0),
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),
fVertex(0),
fNBG(-1)
{
}
AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
AliVEvent(mcEvnt),
fStack(mcEvnt.fStack),
fMCParticles(mcEvnt.fMCParticles),
fMCParticleMap(mcEvnt.fMCParticleMap),
fHeader(mcEvnt.fHeader),
fAODMCHeader(mcEvnt.fAODMCHeader),
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),
fVertex(mcEvnt.fVertex),
fNBG(mcEvnt.fNBG)
{
}
AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
{
if (this!=&mcEvnt) {
AliVEvent::operator=(mcEvnt);
}
return *this;
}
void AliMCEvent::ConnectTreeE (TTree* tree)
{
tree->SetBranchAddress("Header", &fHeader);
}
void AliMCEvent::ConnectTreeK (TTree* tree)
{
fStack = fHeader->Stack();
fStack->ConnectTree(tree);
fStack->GetEvent();
UpdateEventInformation();
}
void AliMCEvent::ConnectHeaderAndStack(AliHeader* header)
{
fHeader = header;
fStack = fHeader->Stack();
UpdateEventInformation();
}
void AliMCEvent::UpdateEventInformation()
{
if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
fNparticles = fStack->GetNtrack();
fNprimaries = fStack->GetNprimary();
Int_t iev = fHeader->GetEvent();
Int_t ievr = fHeader->GetEventNrInRun();
AliDebug(1, Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
iev, ievr, fNparticles, fNprimaries));
if (fMCParticleMap) {
fMCParticleMap->Clear();
fMCParticles->Delete();
if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
}
else
fMCParticleMap = new TObjArray(fNparticles);
}
void AliMCEvent::ConnectTreeTR (TTree* tree)
{
fTreeTR = tree;
if (fTreeTR->GetBranch("AliRun")) {
if (fTmpFileTR) {
fTmpFileTR->Close();
delete fTmpFileTR;
}
ReorderAndExpandTreeTR();
} else {
fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
}
}
Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
{
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()
{
delete fStack; fStack = 0;
if (fTRBuffer) {
fTRBuffer->Delete();
delete fTRBuffer;
fTRBuffer = 0;
}
}
#include <iostream>
void AliMCEvent::FinishEvent()
{
if (fStack) fStack->Reset(0);
fMCParticles->Delete();
if (fMCParticleMap)
fMCParticleMap->Clear();
if (fTRBuffer) {
fTRBuffer->Delete();
}
fTrackReferences->Clear();
fNparticles = -1;
fNprimaries = -1;
fStack = 0;
fSubsidiaryEvents = 0;
fNBG = -1;
}
void AliMCEvent::DrawCheck(Int_t i, Int_t search)
{
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()
{
fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
if (!fTRBuffer) fTRBuffer = new TClonesArray("AliTrackReference", 100);
fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 0);
if (fTreeTR) {
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);
}
TClonesArray* trefs[7];
for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
if (fTreeTR){
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();
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);
Int_t dau1 = part->GetFirstDaughter();
if (dau1 < np) continue;
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--;
}
}
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;
}
}
if (hasHits || isOutside) break;
}
}
if (!hasHits) {
for (Int_t id = dau1; (id <= dau2); id++) {
fTmpTreeTR->Fill();
ifills++;
}
} else {
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();
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) {
tr->SetDetectorId(ib-1);
Int_t nref = fTRBuffer->GetEntriesFast();
TClonesArray &lref = *fTRBuffer;
new(lref[nref]) AliTrackReference(*tr);
}
}
}
fTmpTreeTR->Fill();
fTRBuffer->Delete();
ifills++;
}
}
}
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();
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);
}
}
}
}
it++;
fTmpTreeTR->Fill();
fTRBuffer->Delete();
ifills++;
}
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
{
if (fExternal) {
return ((AliVParticle*) (fMCParticles->At(i)));
}
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;
TObjArray *rarray = 0;
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));
}
if(!fMCParticleMap->At(i)) {
particle = fStack->Particle(i);
if (fTreeTR) {
fTreeTR->GetEntry(fStack->TreeKEntry(i));
trefs = fTRBuffer;
ntref = trefs->GetEntriesFast();
rarray = new TObjArray(ntref);
Int_t nen = fTrackReferences->GetEntriesFast();
for (Int_t j = 0; j < ntref; j++) {
AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
if (ref) {
new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
rarray->AddAt((*fTrackReferences)[nen], j);
nen++;
}
}
}
Int_t nentries = fMCParticles->GetEntriesFast();
mcParticle = new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
fMCParticleMap->AddAt(mcParticle, i);
if (mcParticle) {
TParticle* part = mcParticle->Particle();
Int_t imo = part->GetFirstMother();
Int_t id1 = part->GetFirstDaughter();
Int_t id2 = part->GetLastDaughter();
if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
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() const
{
if (!fExternal) {
return (fHeader->GenEventHeader());
} else {
if (fAODMCHeader) {
TList * lh = fAODMCHeader->GetCocktailHeaders();
if (lh) {return ((AliGenEventHeader*) lh->At(0));}
}
}
return 0;
}
void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event)
{
if (!fSubsidiaryEvents) {
TList* events = new TList();
events->Add(new AliMCEvent(*this));
fSubsidiaryEvents = events;
}
fSubsidiaryEvents->Add(event);
}
AliGenEventHeader *AliMCEvent::FindHeader(Int_t ipart) {
AliMCEvent *event = this;
if (fSubsidiaryEvents) {
ipart = FindIndexAndEvent(ipart,event);
}
AliGenEventHeader* header = event->GenEventHeader();
if (ipart >= header->NProduced()) {
AliWarning(Form("Not a primary -- returning 0 (idx %d, nPrimary %d)",ipart,header->NProduced()));
return 0;
}
AliGenCocktailEventHeader *coHeader = dynamic_cast<AliGenCocktailEventHeader*>(header);
if (coHeader) {
TList* headerList = coHeader->GetHeaders();
TIter headIt(headerList);
Int_t nproduced = 0;
do {
header = (AliGenEventHeader*) headIt();
if (header) nproduced += header->NProduced();
} while (header && ipart >= nproduced);
}
return header;
}
Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
{
TIter next(fSubsidiaryEvents);
next.Reset();
if (oldidx < fNprimaries) {
while((event = (AliMCEvent*)next())) {
if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
}
if (event) {
return (oldidx - event->GetPrimaryOffset());
} else {
return (-1);
}
} else {
while((event = (AliMCEvent*)next())) {
if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
}
if (event) {
return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
} else {
return (-1);
}
}
}
Int_t AliMCEvent::BgLabelToIndex(Int_t label)
{
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) const
{
if (!fSubsidiaryEvents) {
return fStack->IsPhysicalPrimary(i);
} else {
AliMCEvent* evt = 0;
Int_t idx = FindIndexAndEvent(i, evt);
return (evt->IsPhysicalPrimary(idx));
}
}
Bool_t AliMCEvent::IsSecondaryFromWeakDecay(Int_t i)
{
if (!fSubsidiaryEvents) {
return fStack->IsSecondaryFromWeakDecay(i);
} else {
AliMCEvent* evt = 0;
Int_t idx = FindIndexAndEvent(i, evt);
return (evt->IsSecondaryFromWeakDecay(idx));
}
}
Bool_t AliMCEvent::IsSecondaryFromMaterial(Int_t i)
{
if (!fSubsidiaryEvents) {
return fStack->IsSecondaryFromMaterial(i);
} else {
AliMCEvent* evt = 0;
Int_t idx = FindIndexAndEvent(i, evt);
return (evt->IsSecondaryFromMaterial(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());
}
}
}
void AliMCEvent::PreReadAll()
{
Int_t i;
for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++)
{
GetTrack(i);
}
for (i = 0; i < fStack->GetNprimary(); i++)
{
GetTrack(i);
}
AssignGeneratorIndex();
}
const AliVVertex * AliMCEvent::GetPrimaryVertex() const
{
TArrayF v;
GenEventHeader()->PrimaryVertex(v) ;
if (!fVertex) {
fVertex = new AliMCVertex(v[0], v[1], v[2]);
} else {
((AliMCVertex*) fVertex)->SetPosition(v[0], v[1], v[2]);
}
return fVertex;
}
Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
{
if (fNBG == -1) {
AliGenCocktailEventHeader* coHeader =
dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
if (!coHeader) return (0);
TList* list = coHeader->GetHeaders();
AliGenHijingEventHeader* hijingH = dynamic_cast<AliGenHijingEventHeader*>(list->FindObject("Hijing"));
if (!hijingH) return (0);
fNBG = hijingH->NProduced();
}
return (index < fNBG);
}
TList* AliMCEvent::GetCocktailList()
{
if(fExternal==kFALSE) {
AliGenCocktailEventHeader* coHeader =dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
if(!coHeader) {
return 0;
} else {
return (coHeader->GetHeaders());
}
} else {
if(!fAODMCHeader) {
return 0;
} else {
return (fAODMCHeader->GetCocktailHeaders());
}
}
}
TString AliMCEvent::GetGenerator(Int_t index)
{
Int_t nsumpart=fNprimaries;
TList* lh = GetCocktailList();
if(!lh){ TString noheader="nococktailheader";
return noheader;}
Int_t nh=lh->GetEntries();
for (Int_t i = nh-1; i >= 0; i--){
AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
TString genname=gh->GetName();
Int_t npart=gh->NProduced();
if (i == 0) npart = nsumpart;
if(index < nsumpart && index >= (nsumpart-npart)) return genname;
nsumpart-=npart;
}
TString empty="";
return empty;
}
void AliMCEvent::AssignGeneratorIndex() {
TList* list = GetCocktailList();
if (fNprimaries <= 0) {
AliWarning(Form("AliMCEvent::AssignGeneratorIndex: no primaries %10d\n", fNprimaries));
return;
}
if (!list) {
return;
} else {
Int_t nh = list->GetEntries();
Int_t nsumpart = fNprimaries;
for(Int_t i = nh-1; i >= 0; i--){
AliGenEventHeader* gh = (AliGenEventHeader*)list->At(i);
Int_t npart = gh->NProduced();
if (i==0) {
if (npart != nsumpart) {
}
npart = nsumpart;
}
for (Int_t j = nsumpart-1; j >= nsumpart-npart; j--) {
AliVParticle* part = GetTrack(j);
if (!part) {
AliWarning(Form("AliMCEvent::AssignGeneratorIndex: 0-pointer to particle j %8d npart %8d nsumpart %8d Nprimaries %8d\n",
j, npart, nsumpart, fNprimaries));
break;
}
part->SetGeneratorIndex(i);
Int_t dmin = part->GetFirstDaughter();
Int_t dmax = part->GetLastDaughter();
if (dmin == -1) continue;
AssignGeneratorIndex(i, dmin, dmax);
}
nsumpart -= npart;
}
}
}
void AliMCEvent::AssignGeneratorIndex(Int_t index, Int_t dmin, Int_t dmax) {
for (Int_t k = dmin; k <= dmax; k++) {
AliVParticle* dpart = GetTrack(k);
dpart->SetGeneratorIndex(index);
Int_t d1 = dpart->GetFirstDaughter();
Int_t d2 = dpart->GetLastDaughter();
if (d1 > -1) {
AssignGeneratorIndex(index, d1, d2);
}
}
}
Bool_t AliMCEvent::GetCocktailGenerator(Int_t index,TString &nameGen){
AliVParticle* mcpart0 = (AliVParticle*) (GetTrack(index));
if(!mcpart0){
printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",index);
return 0;
}
nameGen=GetGenerator(index);
if(nameGen.Contains("nococktailheader") )return 0;
Int_t lab=index;
while(nameGen.IsWhitespace()){
AliVParticle* mcpart = (AliVParticle*) (GetTrack(lab));
if(!mcpart){
printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",lab);
break;}
Int_t mother=0;
mother = mcpart->GetMother();
if(mother<0){
printf("AliMCEvent - BREAK: Reached primary particle without valid mother\n");
break;
}
AliVParticle* mcmom = (AliVParticle*) (GetTrack(mother));
if(!mcmom){
printf("AliMCEvent-BREAK: No valid AliMCParticle mother at label %i\n",mother);
break;
}
lab=mother;
nameGen=GetGenerator(mother);
}
return 1;
}
void AliMCEvent::SetParticleArray(TClonesArray* mcParticles)
{
fMCParticles = mcParticles;
fNparticles = fMCParticles->GetEntries();
fExternal = kTRUE;
fNprimaries = 0;
struct Local {
static Int_t binaryfirst(TClonesArray* a, Int_t low, Int_t high)
{
Int_t mid = low + (high - low)/2;
if (low > a->GetEntries()-1) return (a->GetEntries()-1);
if (!((AliVParticle*) a->At(mid))->IsPrimary()) {
if (mid > 1 && !((AliVParticle*) a->At(mid-1))->IsPrimary()) {
return binaryfirst(a, low, mid-1);
} else {
return mid;
}
} else {
return binaryfirst(a, mid+1, high);
}
}
};
fNprimaries = Local::binaryfirst(mcParticles, 0, mcParticles->GetEntries()-1);
AssignGeneratorIndex();
}
AliVEvent::EDataLayoutType AliMCEvent::GetDataLayoutType() const
{
return AliVEvent::kMC;
}
ClassImp(AliMCEvent)