]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMCEvent.cxx
Added selection on FiredTriggerClass
[u/mrichter/AliRoot.git] / STEER / AliMCEvent.cxx
index e7b687aec4a1e4d7f65c6175f3d6f8f394c555bf..bccbfb39337302fb7eb4a4c7bea956a8426c3f9f 100644 (file)
 #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"
@@ -56,7 +58,8 @@ AliMCEvent::AliMCEvent():
     fSubsidiaryEvents(0),
     fPrimaryOffset(0),
     fSecondaryOffset(0),
-    fExternal(0)
+    fExternal(0),
+    fVertex(0)
 {
     // Default constructor
 }
@@ -77,7 +80,8 @@ AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
     fSubsidiaryEvents(0),
     fPrimaryOffset(0),
     fSecondaryOffset(0),
-    fExternal(0)
+    fExternal(0),
+    fVertex(mcEvnt.fVertex)
 { 
 // Copy constructor
 }
@@ -111,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
@@ -124,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)
@@ -169,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
@@ -184,16 +183,22 @@ void AliMCEvent::Clean()
     }
 }
 
+#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;
@@ -265,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];
@@ -518,59 +526,56 @@ AliVParticle* AliMCEvent::GetTrack(Int_t i) const
                // 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);
-
-       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);
+       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->SetLabel(i + fSecondaryOffset - fNprimaries);
+               mcParticle->SetFirstDaughter(id1);
+               mcParticle->SetLastDaughter (id2);
+               mcParticle->SetMother       (imo);
            }
-           
-           
-       } 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) 
@@ -594,12 +599,20 @@ Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
        while((event = (AliMCEvent*)next())) {
            if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
        }
-       return (oldidx - event->GetPrimaryOffset());
+       if (event) {
+           return (oldidx - event->GetPrimaryOffset());
+       } else {
+           return (-1);
+       }
     } else {
        while((event = (AliMCEvent*)next())) {
            if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
        }
-       return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
+       if (event) {
+           return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
+       } else {
+           return (-1);
+       }
     }
 }
 
@@ -636,6 +649,8 @@ Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i)
 
 void AliMCEvent::InitEvent()
 {
+//
+// Initialize the subsidiary event structure
     if (fSubsidiaryEvents) {
        TIter next(fSubsidiaryEvents);
        AliMCEvent* evt;
@@ -660,5 +675,34 @@ void AliMCEvent::InitEvent()
     }
 }
 
+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)