The new class AliMCEvent derives from AliVEvent.
AliVEvent implements the method GetTrack which gives back MC information in form of AliMCParticles
deriving from AliVParticle.
--- /dev/null
+/**************************************************************************
+ * 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 <TObjArray.h>
+
+#include "AliLog.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliTrackReference.h"
+#include "AliHeader.h"
+
+
+AliMCEvent::AliMCEvent():
+ AliVEvent(),
+ fStack(0),
+ fMCParticles(new TClonesArray("AliMCParticle",1000)),
+ fMCParticleMap(0),
+ fHeader(0),
+ fTrackReferences(0),
+ fTreeTR(0),
+ fTmpTreeTR(0),
+ fTmpFileTR(0),
+ fNprimaries(-1),
+ fNparticles(-1)
+{
+ // Default constructor
+}
+
+AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
+ AliVEvent(mcEvnt),
+ fStack(0),
+ fMCParticles(0),
+ fMCParticleMap(0),
+ fHeader(0),
+ fTrackReferences(0),
+ fTreeTR(0),
+ fTmpTreeTR(0),
+ fTmpFileTR(0),
+ fNprimaries(-1),
+ fNparticles(-1)
+{
+// 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();
+ AliInfo(Form("AliMCEvent: Number of particles: %5d (all) %5d (primaries)\n",
+ fNparticles, fNprimaries));
+
+ // This is a cache for the TParticles converted to MCParticles on user request
+ if (fMCParticleMap) {
+ fMCParticleMap->Clear();
+ if (fNparticles>0) fMCParticleMap->Expand(fNparticles);}
+ else
+ fMCParticleMap = new TObjArray(fNparticles);
+ fMCParticles->Clear();
+}
+
+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", &fTrackReferences);
+ }
+}
+
+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 = fTrackReferences;
+ 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;
+
+ // Clear TR
+ if (fTrackReferences) {
+ fTrackReferences->Clear();
+ delete fTrackReferences;
+ fTrackReferences = 0;
+ }
+}
+
+void AliMCEvent::FinishEvent()
+{
+ // Clean-up after event
+ fStack->Reset(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 = fTrackReferences->GetEntries();
+
+
+ if (search) {
+ while(nh <= search && i < fNparticles - 1) {
+ i++;
+ fTreeTR->GetEntry(fStack->TreeKEntry(i));
+ nh = fTrackReferences->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*) fTrackReferences->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 (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
+ fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 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 = fTrackReferences->GetEntriesFast();
+ TClonesArray &lref = *fTrackReferences;
+ new(lref[nref]) AliTrackReference(*tr);
+ }
+ } // hits
+ } // branches
+ fTmpTreeTR->Fill();
+ fTrackReferences->Clear();
+ 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 = fTrackReferences->GetEntriesFast();
+ TClonesArray &lref = *fTrackReferences;
+ new(lref[nref]) AliTrackReference(*tr);
+ }
+ } // hits
+ } // branches
+ } // entries
+ it++;
+ fTmpTreeTR->Fill();
+ fTrackReferences->Clear();
+ 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
+ AliMCParticle *mcParticle;
+ TParticle *particle;
+
+ // Out of range check
+ if (i < 0 || i >= fNparticles) {
+ AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
+ mcParticle = 0;
+ return (mcParticle);
+ }
+
+
+ particle = fStack->Particle(i);
+ //
+ // First check of the MC Particle has been already cached
+ if(!fMCParticleMap->At(i)) {
+ Int_t nentries = fMCParticles->GetEntriesFast();
+ new ((*fMCParticles)[nentries]) AliMCParticle(particle);
+ mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
+ fMCParticleMap->AddAt(mcParticle, i);
+ } else {
+ mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
+ }
+
+ return mcParticle;
+}
+
+inline AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
+
+ClassImp(AliMCEvent)
--- /dev/null
+// -*- mode: C++ -*-
+#ifndef ALIMCEVENT_H
+#define ALIMCEVENT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+// Class AliMCEvent
+//
+// Origin: Andreas.Morsch, CERN, andreas.morsch@cern.ch
+//-------------------------------------------------------------------------
+
+
+#include <TTree.h>
+#include <TObjArray.h>
+
+#include <AliVEvent.h>
+#include "AliVHeader.h"
+#include "AliMCParticle.h"
+
+class AliStack;
+class AliHeader;
+class AliGenEventHeader;
+
+class TClonesArray;
+
+class AliMCEvent : public AliVEvent {
+
+public:
+
+ AliMCEvent();
+ virtual ~AliMCEvent() {;}
+ AliMCEvent(const AliMCEvent& mcEvnt);
+ AliMCEvent& operator=(const AliMCEvent& mcEvnt);
+ //
+ // Methods implementing the interface
+ //
+ // Services
+ virtual void AddObject(TObject* /*obj*/) {;}
+ virtual TObject* FindListObject(const char */*name*/) {return 0;}
+ virtual TList* GetList() const {return 0;}
+ virtual void CreateStdContent() {;}
+ virtual void GetStdContent() {;}
+ virtual void ReadFromTree(TTree * /*tree*/) {;}
+ virtual const void WriteToTree(TTree* /*tree*/) const {;}
+
+ virtual void SetStdNames() {;}
+ virtual void Print(Option_t */*option=""*/) const {;}
+
+
+ // Header
+ virtual AliVHeader* GetHeader() const {return 0;}
+
+ // Delegated methods for fESDRun or AODHeader
+
+ virtual void SetRunNumber(Int_t /*n*/) {;}
+ virtual void SetPeriodNumber(UInt_t /*n*/) {;}
+ virtual void SetMagneticField(Double_t /*mf*/) {;}
+
+
+ virtual Int_t GetRunNumber() const {return 0;}
+ virtual UInt_t GetPeriodNumber() const {return 0;}
+ virtual Double_t GetMagneticField() const {return 0.;}
+
+ // Setters not needed
+ virtual void SetOrbitNumber(UInt_t /*n*/) {;}
+ virtual void SetBunchCrossNumber(UShort_t /*n*/) {;}
+ virtual void SetEventType(UInt_t /*eventType*/) {;}
+ virtual void SetTriggerMask(ULong64_t /*n*/) {;}
+ virtual void SetTriggerCluster(UChar_t /*n*/) {;}
+
+ virtual UInt_t GetOrbitNumber() const {return 0;}
+ virtual UShort_t GetBunchCrossNumber() const {return 0;}
+
+ virtual UInt_t GetEventType() const {return 0;}
+
+ virtual ULong64_t GetTriggerMask() const {return 0;}
+ virtual UChar_t GetTriggerCluster() const {return 0;}
+ virtual Double_t GetZDCN1Energy() const {return 0.;}
+ virtual Double_t GetZDCP1Energy() const {return 0.;}
+ virtual Double_t GetZDCN2Energy() const {return 0.;}
+ virtual Double_t GetZDCP2Energy() const {return 0.;}
+ virtual Double_t GetZDCEMEnergy() const {return 0.;}
+ // Tracks
+ virtual AliMCParticle *GetTrack(Int_t i) const;
+ virtual Int_t GetNumberOfTracks() const {return fNparticles;}
+
+ //
+ // MC Specific methods
+ //
+ // Getters
+ AliStack* Stack() {return fStack;}
+ AliHeader* Header() {return fHeader;}
+ AliGenEventHeader* GenEventHeader();
+ // Services
+ virtual void ConnectTreeE (TTree* tree);
+ virtual void ConnectTreeK (TTree* tree);
+ virtual void ConnectTreeTR(TTree* tree);
+ virtual void Clean();
+ virtual void FinishEvent();
+ virtual Int_t GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs);
+ virtual void DrawCheck(Int_t i, Int_t search);
+private:
+ virtual void ReorderAndExpandTreeTR();
+
+private:
+ AliStack *fStack; //! Current pointer to stack
+ TClonesArray *fMCParticles; //! Pointer to list of particles
+ TObjArray *fMCParticleMap; //! Map of MC Particles
+ AliHeader *fHeader; //! Current pointer to header
+ TClonesArray *fTrackReferences; //! Current list of track references
+ TTree *fTreeTR; //! Pointer to Track Reference Tree
+ TTree *fTmpTreeTR; //! Temporary tree TR to read old format
+ TFile *fTmpFileTR; //! Temporary file with TreeTR to read old format
+ Int_t fNprimaries; //! Number of primaries
+ Int_t fNparticles; //! Number of particles
+ ClassDef(AliMCEvent, 0) // AliVEvent realisation for MC data
+};
+
+
+#endif
+
#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliStack.h"
#include <TString.h>
#include <TClonesArray.h>
#include <TDirectoryFile.h>
-#include <TArrow.h>
-#include <TMarker.h>
-#include <TH2F.h>
-#include <TDirectoryFile.h>
-
ClassImp(AliMCEventHandler)
AliMCEventHandler::AliMCEventHandler() :
AliVEventHandler(),
+ fMCEvent(new AliMCEvent()),
fFileE(0),
fFileK(0),
fFileTR(0),
- fTmpFileTR(0),
fTreeE(0),
fTreeK(0),
fTreeTR(0),
- fTmpTreeTR(0),
fDirK(0),
fDirTR(0),
- fStack(0),
- fHeader(0),
- fTrackReferences(0),
fNEvent(-1),
fEvent(-1),
- fNprimaries(-1),
- fNparticles(-1),
fPathName(new TString("./")),
fExtension(""),
fFileNumber(0),
AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
AliVEventHandler(name, title),
+ fMCEvent(new AliMCEvent()),
fFileE(0),
fFileK(0),
fFileTR(0),
- fTmpFileTR(0),
fTreeE(0),
fTreeK(0),
fTreeTR(0),
- fTmpTreeTR(0),
fDirK(0),
fDirTR(0),
- fStack(0),
- fHeader(0),
- fTrackReferences(0),
fNEvent(-1),
fEvent(-1),
- fNprimaries(-1),
- fNparticles(-1),
fPathName(new TString("./")),
fExtension(""),
fFileNumber(0),
AliMCEventHandler::~AliMCEventHandler()
{
// Destructor
+ delete fMCEvent;
delete fFileE;
delete fFileK;
delete fFileTR;
fFileE = TFile::Open(Form("%sgalice.root", fPathName->Data()));
if (!fFileE) AliFatal(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName->Data()));
+ //
+ // Tree E
fFileE->GetObject("TE", fTreeE);
- fTreeE->SetBranchAddress("Header", &fHeader);
+ // Connect Tree E to the MCEvent
+ fMCEvent->ConnectTreeE(fTreeE);
fNEvent = fTreeE->GetEntries();
//
// Tree K
AliInfo(Form("AliMCEventHandler:Number of events in this directory %5d \n", fNEvent));
return kTRUE;
-
}
Bool_t AliMCEventHandler::GetEvent(Int_t iev)
// Load the event number iev
//
// Calculate the file number
- Int_t inew = iev/fEventsPerFile;
+ Int_t inew = iev / fEventsPerFile;
if (inew != fFileNumber) {
fFileNumber = inew;
if (!OpenFile(fFileNumber)){
// Folder name
char folder[20];
sprintf(folder, "Event%d", iev);
-
// TreeE
fTreeE->GetEntry(iev);
- fStack = fHeader->Stack();
// Tree K
fFileK->GetObject(folder, fDirK);
if (!fDirK) {
return kFALSE;
}
fDirK ->GetObject("TreeK", fTreeK);
- fStack->ConnectTree(fTreeK);
- fStack->GetEvent();
+ // Connect TreeK to MCEvent
+ fMCEvent->ConnectTreeK(fTreeK);
//Tree TR
if (fFileTR) {
// Check which format has been read
fFileTR->GetObject(folder, fDirTR);
fDirTR->GetObject("TreeTR", fTreeTR);
- 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", &fTrackReferences);
- }
+ //
+ // Connect TR to MCEvent
+ fMCEvent->ConnectTreeTR(fTreeTR);
}
//
- fNparticles = fStack->GetNtrack();
- fNprimaries = fStack->GetNprimary();
- AliInfo(Form("AliMCEventHandler: Event#: %5d Number of particles: %5d (all) %5d (primaries)\n",
- fEvent, fNparticles, fNprimaries));
-
return kTRUE;
}
Int_t AliMCEventHandler::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 (fFileTR) {
- fTreeTR->GetEntry(fStack->TreeKEntry(i));
- trefs = fTrackReferences;
- return trefs->GetEntries();
- } else {
- trefs = 0;
- return -1;
- }
+ return (fMCEvent->GetParticleAndTR(i, particle, trefs));
}
void AliMCEventHandler::DrawCheck(Int_t i, Int_t search)
{
// Retrieve entry i and draw momentum vector and hits
- if (!fFileTR) {
- AliWarning("No Track Reference information available");
- return;
- }
-
- if (i > -1 && i < fNparticles) {
- fTreeTR->GetEntry(fStack->TreeKEntry(i));
- } else {
- AliWarning("AliMCEventHandler::GetEntry: Index out of range");
- }
-
- Int_t nh = fTrackReferences->GetEntries();
-
-
- if (search) {
- while(nh <= search && i < fNparticles - 1) {
- i++;
- fTreeTR->GetEntry(fStack->TreeKEntry(i));
- nh = fTrackReferences->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*) fTrackReferences->At(ih);
- TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
- m->Draw();
- m->SetMarkerSize(0.4);
-
- }
+ fMCEvent->DrawCheck(i, search);
}
Bool_t AliMCEventHandler::Notify(const char *path)
void AliMCEventHandler::ResetIO()
{
// Clear header and stack
-
- if (fHeader) {
- delete fHeader;
- fHeader = 0;
- }
-
- delete fStack;
- delete fTreeE; fTreeE = 0;
-
-// Clear TR
+ fMCEvent->Clean();
- if (fTrackReferences) {
- fTrackReferences->Clear();
- delete fTrackReferences;
- fTrackReferences = 0;
- }
+// Delete Tree E
+ delete fTreeE; fTreeE = 0;
// Reset files
if (fFileE) delete fFileE;
// Clean-up after each event
delete fDirTR; fDirTR = 0;
delete fDirK; fDirK = 0;
- Stack()->Reset(0);
+ fMCEvent->FinishEvent();
return kTRUE;
}
return kTRUE;
}
-void AliMCEventHandler::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 (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference", 100);
- fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 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 = fTrackReferences->GetEntriesFast();
- TClonesArray &lref = *fTrackReferences;
- new(lref[nref]) AliTrackReference(*tr);
- }
- } // hits
- } // branches
- fTmpTreeTR->Fill();
- fTrackReferences->Clear();
- 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 = fTrackReferences->GetEntriesFast();
- TClonesArray &lref = *fTrackReferences;
- new(lref[nref]) AliTrackReference(*tr);
- }
- } // hits
- } // branches
- } // entries
- it++;
- fTmpTreeTR->Fill();
- fTrackReferences->Clear();
- 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("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
- ifills, fStack->GetNtrack());
-
- fTmpTreeTR->Write();
- fTreeTR = fTmpTreeTR;
-}
void AliMCEventHandler::SetInputPath(char* fname)
{
//-------------------------------------------------------------------------
// Class AliMCEvent
// This class gives access to MC truth during the analysis.
-// Monte Carlo truth is containe in the kinematics tree (produced particles) and
+// Monte Carlo truth is contained in the kinematics tree (produced particles) and
// the tree of reference hits.
//
// Origin: Andreas Morsch, CERN, andreas.morsch@cern.ch
//-------------------------------------------------------------------------
-
#include "AliVEventHandler.h"
#include "AliHeader.h"
class TFile;
class TClonesArray;
class TDirectoryFile;
-class AliHeader;
-class AliGenEventHeader;
-class AliStack;
+class AliMCEvent;
+
class AliMCEventHandler : public AliVEventHandler
virtual void ResetIO();
virtual Bool_t GetEvent(Int_t iev);
//
- AliStack* Stack() {return fStack;}
- AliHeader* Header() {return fHeader;}
- AliGenEventHeader* GenEventHeader() {return (fHeader->GenEventHeader());}
- TTree* TreeTR() {return fTreeTR;}
- Int_t GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs);
- void DrawCheck(Int_t i, Int_t search=0);
+ AliMCEvent* MCEvent() {return fMCEvent;}
+ TTree* TreeTR() {return fTreeTR;}
+ Int_t GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs);
+ void DrawCheck(Int_t i, Int_t search=0);
private:
- Bool_t OpenFile(Int_t i);
- void ReorderAndExpandTreeTR();
-
+ Bool_t OpenFile(Int_t i);
private:
+ AliMCEvent *fMCEvent; //! MC Event
TFile *fFileE; //! File with TreeE
TFile *fFileK; //! File with TreeK
TFile *fFileTR; //! File with TreeTR
- TFile *fTmpFileTR; //! Temporary file with TreeTR to read old format
TTree *fTreeE; //! TreeE (Event Headers)
TTree *fTreeK; //! TreeK (kinematics tree)
TTree *fTreeTR; //! TreeTR (track references tree)
- TTree *fTmpTreeTR; //! Temporary tree TR to read old format
TDirectoryFile *fDirK; //! Directory for Kine Tree
TDirectoryFile *fDirTR; //! Directory for TR Tree
- AliStack *fStack; //! Current pointer to stack
- AliHeader *fHeader; //! Current pointer to header
- TClonesArray *fTrackReferences; //! Current list of track references
Int_t fNEvent; //! Number of events
Int_t fEvent; //! Current event
- Int_t fNprimaries; //! Number of primaries
- Int_t fNparticles; //! Number of particles
TString *fPathName; //! Input file path
char *fExtension; //! File name extension
Int_t fFileNumber; //! Input file number
--- /dev/null
+/**************************************************************************
+ * 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$ */
+
+//-------------------------------------------------------------------------
+// Realisation of AliVParticle for MC Particles
+// Implementation wraps a TParticle and delegates the methods
+// Author: Andreas Morsch, CERN
+//-------------------------------------------------------------------------
+
+#include "AliMCParticle.h"
+
+
+ClassImp(AliMCParticle)
+
+AliMCParticle::AliMCParticle():
+ AliVParticle(),
+ fParticle(0)
+{
+ // Constructor
+}
+
+AliMCParticle::AliMCParticle(TParticle* part):
+ AliVParticle(),
+ fParticle(part)
+{
+ // Constructor
+}
+
+
+AliMCParticle::AliMCParticle(const AliMCParticle& mcPart) :
+ AliVParticle(mcPart),
+ fParticle(0)
+{
+// Copy constructor
+}
+
+AliMCParticle& AliMCParticle::operator=(const AliMCParticle& mcPart)
+{ if (this!=&mcPart) {
+ AliVParticle::operator=(mcPart);
+ }
+
+ return *this;
+}
--- /dev/null
+#ifndef AliMCParticle_H
+#define AliMCParticle_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//-------------------------------------------------------------------------
+// AliVParticle realisation for MC Particles
+// Author: Markus Oldenburg, CERN
+//-------------------------------------------------------------------------
+
+#include <Rtypes.h>
+#include <AliVParticle.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+
+class AliMCParticle: public AliVParticle {
+
+public:
+ AliMCParticle();
+ AliMCParticle(TParticle* part);
+ virtual ~AliMCParticle() {}
+ AliMCParticle(const AliMCParticle& mcPart);
+ AliMCParticle& operator=(const AliMCParticle& mcPart);
+
+ // kinematics
+ virtual Double_t Px() const;
+ virtual Double_t Py() const;
+ virtual Double_t Pz() const;
+ virtual Double_t Pt() const;
+ virtual Double_t P() const;
+
+ virtual Double_t OneOverPt() const;
+ virtual Double_t Phi() const;
+ virtual Double_t Theta() const;
+
+
+ virtual Double_t E() const;
+ virtual Double_t M() const;
+
+ virtual Double_t Eta() const;
+ virtual Double_t Y() const;
+
+ virtual Short_t Charge() const;
+
+ // PID
+ virtual const Double_t *PID() const {return 0;} // return PID object (to be defined, still)
+
+ private:
+ TParticle *fParticle; // The TParticle
+
+ ClassDef(AliMCParticle,0) // AliVParticle realisation for MCParticles
+};
+
+inline Double_t AliMCParticle::Px() const {return fParticle->Px();}
+inline Double_t AliMCParticle::Py() const {return fParticle->Py();}
+inline Double_t AliMCParticle::Pz() const {return fParticle->Pz();}
+inline Double_t AliMCParticle::Pt() const {return fParticle->Pt();}
+inline Double_t AliMCParticle::P() const {return fParticle->P(); }
+inline Double_t AliMCParticle::OneOverPt() const {return 1. / fParticle->Pt();}
+inline Double_t AliMCParticle::Phi() const {return fParticle->Phi();}
+inline Double_t AliMCParticle::Theta() const {return fParticle->Theta();}
+inline Double_t AliMCParticle::E() const {return fParticle->Energy();}
+inline Double_t AliMCParticle::Eta() const {return fParticle->Eta();}
+
+inline Double_t AliMCParticle::M() const
+{
+ TParticlePDG* pdg = fParticle->GetPDG();
+ if (pdg) {
+ return (pdg->Mass());
+ } else {
+ return (fParticle->GetCalcMass());
+ }
+}
+
+
+inline Double_t AliMCParticle::Y() const
+{
+ Double_t e = E();
+ Double_t pz = TMath::Abs(Pz());
+ if (e != pz) {
+ return 0.5*TMath::Log((e+pz)/(e-pz));
+ } else {
+ return -999.;
+ }
+}
+
+inline Short_t AliMCParticle::Charge() const
+{
+ TParticlePDG* pdg = fParticle->GetPDG();
+ if (pdg) {
+ return (Short_t (pdg->Charge()));
+ } else {
+ return -99;
+ }
+}
+
+#endif
#pragma link C++ class AliTrackReference+;
#pragma link C++ class AliSysInfo+;
+#pragma link C++ class AliMCEvent+;
+#pragma link C++ class AliMCParticle+;
+
#endif
AliGenPythiaEventHeader.cxx AliGenCocktailEventHeader.cxx \
AliGenHijingEventHeader.cxx AliCollisionGeometry.cxx \
AliStack.cxx AliMCEventHandler.cxx \
- AliTrackReference.cxx AliSysInfo.cxx
+ AliTrackReference.cxx AliSysInfo.cxx \
+ AliMCEvent.cxx AliMCParticle.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff