#include <TFile.h>
#include <TParticle.h>
#include <TClonesArray.h>
-#include <TRefArray.h>
+//#include <TRefArray.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"
+Int_t AliMCEvent::fgkBgLabelOffset(10000000);
+
+
AliMCEvent::AliMCEvent():
AliVEvent(),
fStack(0),
- fMCParticles(new TClonesArray("AliMCParticle",1000)),
+ fMCParticles(0),
fMCParticleMap(0),
fHeader(new AliHeader()),
fTRBuffer(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(new AliHeader()),
- 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
}
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);
//
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
if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
}
else
- fMCParticleMap = new TRefArray(fNparticles);
+ fMCParticleMap = new TObjArray(fNparticles);
}
void AliMCEvent::ConnectTreeTR (TTree* tree)
void AliMCEvent::Clean()
{
// Clean-up before new trees are connected
-// if (fHeader) {
-// delete fHeader;
-// fHeader = 0;
-// }
-
delete fStack; fStack = 0;
// Clear TR
}
}
+#include <iostream>
+
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;
}
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];
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);
// 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++;
+ 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<AliMCParticle*>((*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<AliMCParticle*>(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)