X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMCEvent.cxx;h=bccbfb39337302fb7eb4a4c7bea956a8426c3f9f;hb=cd995490fc657af53fda843eac38103f1929f980;hp=fe1775601d4024d6a9a45aedc71f1c1408e74270;hpb=568d5aedc444df316622b028b9853c0f9e550cda;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMCEvent.cxx b/STEER/AliMCEvent.cxx index fe1775601d4..bccbfb39337 100644 --- a/STEER/AliMCEvent.cxx +++ b/STEER/AliMCEvent.cxx @@ -26,46 +26,62 @@ #include #include #include -#include +//#include +#include +#include #include "AliLog.h" #include "AliMCEvent.h" +#include "AliMCVertex.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)), + fMCParticles(0), fMCParticleMap(0), - fHeader(0), + fHeader(new AliHeader()), fTRBuffer(0), fTrackReferences(new TClonesArray("AliTrackReference", 1000)), fTreeTR(0), fTmpTreeTR(0), fTmpFileTR(0), fNprimaries(-1), - fNparticles(-1) + fNparticles(-1), + fSubsidiaryEvents(0), + fPrimaryOffset(0), + fSecondaryOffset(0), + fExternal(0), + fVertex(0) { // Default constructor } AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) : AliVEvent(mcEvnt), - fStack(0), - fMCParticles(0), - fMCParticleMap(0), - fHeader(0), - fTRBuffer(0), - fTrackReferences(0), - fTreeTR(0), - fTmpTreeTR(0), - fTmpFileTR(0), - fNprimaries(-1), - fNparticles(-1) + 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), + fVertex(mcEvnt.fVertex) { // Copy constructor } @@ -90,6 +106,8 @@ void AliMCEvent::ConnectTreeE (TTree* tree) 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); // @@ -97,10 +115,10 @@ void AliMCEvent::ConnectTreeK (TTree* tree) 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", + AliDebug(1, 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 @@ -110,7 +128,7 @@ void AliMCEvent::ConnectTreeK (TTree* tree) if (fNparticles>0) fMCParticleMap->Expand(fNparticles); } else - fMCParticleMap = new TRefArray(fNparticles); + fMCParticleMap = new TObjArray(fNparticles); } void AliMCEvent::ConnectTreeTR (TTree* tree) @@ -155,11 +173,6 @@ Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& void AliMCEvent::Clean() { // Clean-up before new trees are connected - if (fHeader) { - delete fHeader; - fHeader = 0; - } - delete fStack; fStack = 0; // Clear TR @@ -170,16 +183,27 @@ void AliMCEvent::Clean() } } +#include + void AliMCEvent::FinishEvent() { // Clean-up after event // - fStack->Reset(0); + if (fStack) fStack->Reset(0); fMCParticles->Delete(); - fMCParticleMap->Clear(); - if (fTRBuffer) + + if (fMCParticleMap) + fMCParticleMap->Clear(); + if (fTRBuffer) { fTRBuffer->Delete(); - fTrackReferences->Delete(); + } + // fTrackReferences->Delete(); + fTrackReferences->Clear(); + fNparticles = -1; + fNprimaries = -1; + fStack = 0; +// fSubsidiaryEvents->Clear(); + fSubsidiaryEvents = 0; } @@ -246,19 +270,22 @@ 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, 32000, 0); + fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 64000, 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); + 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); + } + // // Connect the active branches TClonesArray* trefs[7]; @@ -439,24 +466,52 @@ void AliMCEvent::ReorderAndExpandTreeTR() fTreeTR = fTmpTreeTR; } -AliMCParticle* AliMCEvent::GetTrack(Int_t i) const +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 of the MC Particle has been already cached + // First check If the MC Particle has been already cached if(!fMCParticleMap->At(i)) { // Get particle from the stack particle = fStack->Particle(i); @@ -471,21 +526,183 @@ AliMCParticle* AliMCEvent::GetTrack(Int_t i) const // Save the track references in a TClonesArray AliTrackReference* ref = dynamic_cast((*fTRBuffer)[j]); // Save the pointer in a TRefArray - new ((*fTrackReferences)[nen]) AliTrackReference(*ref); - rarray->AddAt((*fTrackReferences)[nen], j); - nen++; + if (ref) { + 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((*fMCParticles)[nentries]); + 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) { + // 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(fMCParticleMap->At(i)); } return mcParticle; } - AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());} +AliGenEventHeader* AliMCEvent::GenEventHeader() const {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; + } + 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) +{ + // 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() +{ +// +// Initialize the subsidiary event structure + 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() +{ + // Preread the MC information + Int_t i; + // secondaries + for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++) + { + GetTrack(i); + } + // primaries + for (i = 0; i < fStack->GetNprimary(); i++) + { + GetTrack(i); + } +} + +const AliVVertex * AliMCEvent::GetPrimaryVertex() const +{ + // Create a MCVertex object from the MCHeader information + 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; +} + ClassImp(AliMCEvent)