MCEvent and MCEventHandler can be used as containers. This allows us to compose MC...
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Jul 2009 13:10:55 +0000 (13:10 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Jul 2009 13:10:55 +0000 (13:10 +0000)
Main use case: analysis of merged events.

ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx
STEER/AliAODHandler.cxx
STEER/AliAODMCParticle.cxx
STEER/AliAODMCParticle.h
STEER/AliMCEvent.cxx
STEER/AliMCEvent.h
STEER/AliMCEventHandler.cxx
STEER/AliMCEventHandler.h
STEER/AliMCParticle.cxx
STEER/AliMCParticle.h

index e2ca73deb1fa2acda327ef120dcc4128624e7272..8409390d5fd460af29466fc5f2f57699eca58982 100644 (file)
@@ -164,9 +164,8 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
     return;
   }
 
-  AliStack* stack = mcE->Stack();
   Int_t np    = mcE->GetNumberOfTracks();
-  Int_t nprim = stack->GetNprimary();
+  Int_t nprim = mcE->GetNumberOfPrimaries();
   // TODO ADD MC VERTEX
   
   // We take all real primaries
@@ -189,30 +188,30 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
     } else if (part->GetUniqueID() == 4) {
       // Particles from decay
       // Check that the decay chain ends at a primary particle
-      TParticle* mother = part;
-      Int_t imo = part->GetFirstMother();
+      AliMCParticle* mother = mcpart;
+      Int_t imo = mcpart->GetMother();
       while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
-       mother =  mcE->Stack()->Particle(imo);
-       imo =  mother->GetFirstMother();
+       mother = mcE->GetTrack(imo);
+       imo =  mother->GetMother();
       }
       // Select according to pseudorapidity and production point of primary ancestor
-      if (imo < nprim && Select(mcE->Stack()->Particle(imo), rv, zv))write = kTRUE;         
+      if (imo < nprim && Select(mcE->GetTrack(imo)->Particle(), rv, zv))write = kTRUE;         
     } else if (part->GetUniqueID() == 5) {
       // Now look for pair production
-      Int_t imo = part->GetFirstMother();
+      Int_t imo = mcpart->GetMother();
       if (imo < nprim) {
        // Select, if the gamma is a primary
        write = kTRUE;
       } else {
        // Check if the gamma comes from the decay chain of a primary particle
-       TParticle* mother = mcE->Stack()->Particle(imo);
-       imo = mother->GetFirstMother();
+       AliMCParticle* mother = mcE->GetTrack(imo);
+       imo = mother->GetMother();
        while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
-         mother =  mcE->Stack()->Particle(imo);
-         imo =  mother->GetFirstMother();
+         mother =  mcE->GetTrack(imo);
+         imo =  mother->GetMother();
        }
        // Select according to pseudorapidity and production point 
-       if (imo < nprim && Select(mother, rv, zv)) 
+       if (imo < nprim && Select(mother->Particle(), rv, zv)) 
          write = kTRUE;
       }
     }
@@ -251,12 +250,12 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
       // cases 
       iCharm++;
       Printf("Decay Mother %s",part->GetPDG()->GetName());
-      Int_t d0 =  part->GetFirstDaughter();
-      Int_t d1 =  part->GetLastDaughter();
+      Int_t d0 =  mcpart->GetFirstDaughter();
+      Int_t d1 =  mcpart->GetLastDaughter();
       if(d0>0&&d1>0){
        for(int id = d0;id <= d1;id++){
-         TParticle* daughter =  mcE->Stack()->Particle(id);
-         Printf("Decay Daughter %s",daughter->GetPDG()->GetName());
+         AliMCParticle* daughter =  mcE->GetTrack(id);
+         Printf("Decay Daughter %s",daughter->Particle()->GetPDG()->GetName());
          iAll++;
          if(mcH->IsParticleSelected(id))iTaken++;
        }
index a957b5ea565b82cbafd9f01599679087e85ecbd8..596d4ab2771ddaa93c24ce2a4b14544ee20ba235 100644 (file)
@@ -220,88 +220,87 @@ void AliAODHandler::StoreMCParticles(){
 
 
   // Store the AliAODParticlesMC
-
-  Int_t np    = pStack->GetNtrack();
-  Int_t nprim = pStack->GetNprimary();
+  AliMCEvent* mcEvent = fMCEventH->MCEvent();
+  
+  Int_t np    = mcEvent->GetNumberOfTracks();
+  Int_t nprim = mcEvent->GetNumberOfPrimaries();
 
 
   Int_t j = 0;
   TClonesArray& l = *mcarray;
 
-  for(int i = 0;i < np;++i){
-    if(fMCEventH->IsParticleSelected(i)){
-
-      Int_t flag = 0;
-      TParticle *part = pStack->Particle(i);
-      if(i<nprim)flag |= AliAODMCParticle::kPrimary;
-      if(pStack->IsPhysicalPrimary(i))flag |= AliAODMCParticle::kPhysicalPrim;
-
-      if(fMCEventH->GetNewLabel(i)!=j){
-       AliError(Form("MISMATCH New label %d j: %d",fMCEventH->GetNewLabel(i),j));
-      }
-      AliAODMCParticle mcpart_tmp(part,i,flag);
-
-      // 
-      Int_t d0 =  mcpart_tmp.GetDaughter(0);
-      Int_t d1 =  mcpart_tmp.GetDaughter(1);
-      Int_t m =  mcpart_tmp.GetMother();
-
-      // other than for the track labels, negative values mean
-      // no daughter/mother so preserve it
-      if(d0<0 && d1<0){
-       // no first daughter -> no second daughter
-       // nothing to be done
-       // second condition not needed just for sanity check at the end
-       mcpart_tmp.SetDaughter(0,d0);
-       mcpart_tmp.SetDaughter(1,d1);
-      }
-      else if(d1 < 0 && d0 >= 0){
-       // Only one daughter
-       // second condition not needed just for sanity check at the end
-       if(fMCEventH->IsParticleSelected(d0)){
-         mcpart_tmp.SetDaughter(0,fMCEventH->GetNewLabel(d0));
-       }
-       else{
-         mcpart_tmp.SetDaughter(0,-1);
-       }
-       mcpart_tmp.SetDaughter(1,d1);
-      }
-      else if (d0 > 0 && d1 > 0 ){
-       // we have two or more daughters loop on the stack to see if they are
-       // selected
-       Int_t d0_tmp = -1;
-       Int_t d1_tmp = -1;
-       for(int id = d0; id<=d1;++id){
-         if(fMCEventH->IsParticleSelected(id)){
-           if(d0_tmp==-1){
-             // first time
-             d0_tmp = fMCEventH->GetNewLabel(id);
-             d1_tmp = d0_tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
-           }
-           else d1_tmp = fMCEventH->GetNewLabel(id);
+  for(int i = 0; i < np; ++i){
+      if(fMCEventH->IsParticleSelected(i)){
+         Int_t flag = 0;
+         AliMCParticle* mcpart =  mcEvent->GetTrack(i);
+         TParticle *part = mcpart->Particle();
+         if(i<nprim)flag |= AliAODMCParticle::kPrimary;
+         
+         if(mcEvent->IsPhysicalPrimary(i))flag |= AliAODMCParticle::kPhysicalPrim;
+         
+         if(fMCEventH->GetNewLabel(i)!=j){
+             AliError(Form("MISMATCH New label %d j: %d",fMCEventH->GetNewLabel(i),j));
          }
-       }
-       mcpart_tmp.SetDaughter(0,d0_tmp);
-       mcpart_tmp.SetDaughter(1,d1_tmp);
-      }
-      else{
-       AliError(Form("Unxpected indices %d %d",d0,d1));
-      }
 
-      if(m<0){
-       mcpart_tmp.SetMother(m);
-      }
-      else{
-       if(fMCEventH->IsParticleSelected(m))mcpart_tmp.SetMother(fMCEventH->GetNewLabel(m));
-       else AliError("PROBLEM Mother not selected");
+         AliAODMCParticle mcpart_tmp(mcpart,i,flag);
+         
+         // 
+         Int_t d0 =  mcpart_tmp.GetDaughter(0);
+         Int_t d1 =  mcpart_tmp.GetDaughter(1);
+         Int_t m =   mcpart_tmp.GetMother();
+         
+         // other than for the track labels, negative values mean
+         // no daughter/mother so preserve it
+         
+         if(d0<0 && d1<0){
+             // no first daughter -> no second daughter
+             // nothing to be done
+             // second condition not needed just for sanity check at the end
+             mcpart_tmp.SetDaughter(0,d0);
+             mcpart_tmp.SetDaughter(1,d1);
+         } else if(d1 < 0 && d0 >= 0) {
+             // Only one daughter
+             // second condition not needed just for sanity check at the end
+             if(fMCEventH->IsParticleSelected(d0)){
+                 mcpart_tmp.SetDaughter(0,fMCEventH->GetNewLabel(d0));
+             } else {
+                 mcpart_tmp.SetDaughter(0,-1);
+             }
+             mcpart_tmp.SetDaughter(1,d1);
+         }
+         else if (d0 > 0 && d1 > 0 ){
+             // we have two or more daughters loop on the stack to see if they are
+             // selected
+             Int_t d0_tmp = -1;
+             Int_t d1_tmp = -1;
+             for(int id = d0; id<=d1;++id){
+                 if(fMCEventH->IsParticleSelected(id)){
+                     if(d0_tmp==-1){
+                         // first time
+                         d0_tmp = fMCEventH->GetNewLabel(id);
+                         d1_tmp = d0_tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
+                     }
+                     else d1_tmp = fMCEventH->GetNewLabel(id);
+                 }
+             }
+             mcpart_tmp.SetDaughter(0,d0_tmp);
+             mcpart_tmp.SetDaughter(1,d1_tmp);
+         } else {
+             AliError(Form("Unxpected indices %d %d",d0,d1));
+         }
+         
+         if(m<0){
+             mcpart_tmp.SetMother(m);
+         } else {
+             if(fMCEventH->IsParticleSelected(m))mcpart_tmp.SetMother(fMCEventH->GetNewLabel(m));
+             else AliError(Form("PROBLEM Mother not selected %d \n", m));
+         }
+         
+         new (l[j++]) AliAODMCParticle(mcpart_tmp);
+         
       }
-
-      new (l[j++]) AliAODMCParticle(mcpart_tmp);
-      
-    }
   }
-  AliInfo(Form("AliAODHandler::StoreMCParticles: Selected %d (Primaries %d / total %d) after validation",
+  AliInfo(Form("AliAODHandler::StoreMCParticles: Selected %d (Primaries %d / total %d) after validation \n",
               j,nprim,np));
   
   // Set the labels in the AOD output...
@@ -313,13 +312,17 @@ void AliAODHandler::StoreMCParticles(){
     for(int it = 0; it < fAODEvent->GetNTracks();++it){
       AliAODTrack *track = fAODEvent->GetTrack(it);
       
-      if(TMath::Abs(track->GetLabel())>np||track->GetLabel()==0){
-       AliWarning(Form("Wrong ESD track label %d",track->GetLabel()));
+      Int_t label = TMath::Abs(track->GetLabel());
+      if (label >= AliMCEvent::BgLabelOffset()) label =  mcEvent->BgLabelToIndex(label);
+      
+      
+      if(label > np || track->GetLabel() == 0){
+       AliWarning(Form("Wrong ESD track label %5d (%5d)",track->GetLabel(), label));
       }
-      if(fMCEventH->GetNewLabel(track->GetLabel())==0){
-       AliWarning(Form("New label not found for %d",track->GetLabel()));
+      if(fMCEventH->GetNewLabel(label) == 0){
+       AliWarning(Form("New label not found for %5d (%5d)",track->GetLabel(), label));
       }
-      track->SetLabel(fMCEventH->GetNewLabel(track->GetLabel()));
+      track->SetLabel(fMCEventH->GetNewLabel(label));
     }
   }
   
index bc7fcaa2175d82b887ef96754fb199542844e6bc..bbe75d763db4b88a1f06e57a5628767088aaf33b 100644 (file)
@@ -53,24 +53,24 @@ AliVParticle(),
 }
 
     
-AliAODMCParticle::AliAODMCParticle(TParticle* part, Int_t label,Int_t flag):
+AliAODMCParticle::AliAODMCParticle(AliMCParticle* mcpart, Int_t label,Int_t flag):
     AliVParticle(),
-    fPdgCode(part->GetPdgCode()),
+    fPdgCode(mcpart->Particle()->GetPdgCode()),
     fFlag(flag),
     fLabel(label),
-    fMother(part->GetMother(0)),
-    fPx(part->Px()),
-    fPy(part->Py()),
-    fPz(part->Pz()),
-    fE(part->Energy()),
-    fVx(part->Vx()),
-    fVy(part->Vy()),
-    fVz(part->Vz())
+    fMother(mcpart->GetMother()),
+    fPx(mcpart->Particle()->Px()),
+    fPy(mcpart->Particle()->Py()),
+    fPz(mcpart->Particle()->Pz()),
+    fE(mcpart->Particle()->Energy()),
+    fVx(mcpart->Particle()->Vx()),
+    fVy(mcpart->Particle()->Vy()),
+    fVz(mcpart->Particle()->Vz())
 {
-  fDaughter[0] =  part->GetDaughter(0); 
-  fDaughter[1] =  part->GetDaughter(1);
+  fDaughter[0] =  mcpart->GetFirstDaughter(); 
+  fDaughter[1] =  mcpart->GetLastDaughter();
   // Set unique id
-  TObject::SetUniqueID(part->GetUniqueID());
+  TObject::SetUniqueID(mcpart->Particle()->GetUniqueID());
 }
     
     
index 6f7d529fc89d72cf010371aee61b10ecd580f1cd..9e38e6d7c84adb65d31f4b5b6cf163d8c2f27183 100644 (file)
@@ -18,6 +18,7 @@
 
 #include "AliTrackReference.h"
 #include "AliVParticle.h"
+#include "AliMCParticle.h"
 
 class AliAODEvent;
 class TParticle;
@@ -26,7 +27,7 @@ class TClonesArray;
 class AliAODMCParticle: public AliVParticle {
  public:
   AliAODMCParticle();
-  AliAODMCParticle(TParticle* part, Int_t label=0,Int_t flag = 0);
+  AliAODMCParticle(AliMCParticle* part, Int_t label=0,Int_t flag = 0);
   virtual ~AliAODMCParticle(){};
   AliAODMCParticle(const AliAODMCParticle& mcPart); 
   AliAODMCParticle& operator=(const AliAODMCParticle& mcPart);
index e7a58faa7800b217e47d196a2e60b967fac1d688..286c9a96716bcb697aad244616a083c4b591e7e1 100644 (file)
@@ -27,6 +27,7 @@
 #include <TParticle.h>
 #include <TClonesArray.h>
 #include <TRefArray.h>
+#include <TList.h>
 
 #include "AliLog.h"
 #include "AliMCEvent.h"
@@ -36,6 +37,9 @@
 #include "AliGenEventHeader.h"
 
 
+Int_t AliMCEvent::fgkBgLabelOffset(10000000);
+
+
 AliMCEvent::AliMCEvent():
     AliVEvent(),
     fStack(0),
@@ -48,24 +52,30 @@ AliMCEvent::AliMCEvent():
     fTmpTreeTR(0),
     fTmpFileTR(0),
     fNprimaries(-1),
-    fNparticles(-1)
+    fNparticles(-1),
+    fSubsidiaryEvents(0),
+    fPrimaryOffset(0),
+    fSecondaryOffset(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)
 { 
 // Copy constructor
 }
@@ -183,7 +193,8 @@ void AliMCEvent::FinishEvent()
     fNparticles = -1;
     fNprimaries = -1;    
     fStack      =  0;
-    
+//    fSubsidiaryEvents->Clear();
+    fSubsidiaryEvents = 0;
 }
 
 
@@ -446,21 +457,43 @@ void AliMCEvent::ReorderAndExpandTreeTR()
 AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
 {
     // Get MC Particle 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);
@@ -484,12 +517,138 @@ AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
        new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
        mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
        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);
+           } 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() {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;
+       }
+       return (oldidx - event->GetPrimaryOffset());
+    } else {
+       while((event = (AliMCEvent*)next())) {
+           if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
+       }
+       return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
+    }
+}
+
+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()
+{
+    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());      
+       }
+    }
+}
+
 
 ClassImp(AliMCEvent)
index cd47b2e4edac7f433736de079b5f097b2c49c23f..57ddc248993475923d6bfcc615fae1dc8348248f 100644 (file)
@@ -26,6 +26,7 @@ class AliHeader;
 class AliGenEventHeader;
 
 class TClonesArray;
+class TList;
 
 class AliMCEvent : public AliVEvent {
 
@@ -103,24 +104,39 @@ public:
     virtual void      ConnectTreeK (TTree* tree);
     virtual void      ConnectTreeTR(TTree* tree);
     virtual void      Clean();
+    virtual void      InitEvent();
     virtual void      FinishEvent();
     virtual Int_t     GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs);
     virtual void      DrawCheck(Int_t i, Int_t search);
+    virtual void      AddSubsidiaryEvent(AliMCEvent* event);
+    virtual Int_t     GetNumberOfPrimaries() {return fNprimaries;}
+    virtual Int_t     GetPrimaryOffset()    const {return fPrimaryOffset;}
+    virtual Int_t     GetSecondaryOffset()  const {return fSecondaryOffset;}    
+    virtual void      SetPrimaryOffset(Int_t ioff)    {fPrimaryOffset = ioff;}
+    virtual void      SetSecondaryOffset(Int_t ioff)  {fSecondaryOffset = ioff;}    
+    virtual Bool_t    IsPhysicalPrimary(Int_t i);
+    virtual Int_t     BgLabelToIndex(Int_t label);
+    static  Int_t     BgLabelOffset() {return fgkBgLabelOffset;}
+
 private:
     virtual void      ReorderAndExpandTreeTR();
-    
-private:
-    AliStack         *fStack;           // Current pointer to stack
-    TClonesArray     *fMCParticles;     // Pointer to list of particles
-    TRefArray        *fMCParticleMap;   // Map of MC Particles
-    AliHeader        *fHeader;          // Current pointer to header
-    TClonesArray     *fTRBuffer;        // Track reference buffer    
-    TClonesArray     *fTrackReferences; // Array 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
+    virtual Int_t     FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const;
+private: 
+    AliStack         *fStack;            // Current pointer to stack
+    TClonesArray     *fMCParticles;      // Pointer to list of particles
+    TRefArray        *fMCParticleMap;    // Map of MC Particles
+    AliHeader        *fHeader;           // Current pointer to header
+    TClonesArray     *fTRBuffer;         // Track reference buffer    
+    TClonesArray     *fTrackReferences;  // Array 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
+    TList            *fSubsidiaryEvents; // List of possible subsidiary events (for example merged underlying event) 
+    Int_t             fPrimaryOffset;    // Offset for primaries
+    Int_t             fSecondaryOffset;  // Offset for secondaries
+    static Int_t      fgkBgLabelOffset;  // Standard branch name    
     ClassDef(AliMCEvent, 1)  // AliVEvent realisation for MC data
 };
 
index a4a9824584a8d187a187b5ea2ea216543439a500..6a67ab2ef29fa9c97fa8d70c17290f86d8b62598 100644 (file)
@@ -1,4 +1,5 @@
-/************************************************************************* * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved *
+/************************************************************************** 
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved  *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
  * Contributors are mentioned in the code where appropriate.              *
@@ -26,6 +27,7 @@
 
 #include "AliMCEventHandler.h"
 #include "AliMCEvent.h"
+#include "AliMCParticle.h"
 #include "AliPDG.h"
 #include "AliTrackReference.h"
 #include "AliHeader.h"
@@ -34,6 +36,7 @@
 
 #include <TTree.h>
 #include <TFile.h>
+#include <TList.h>
 #include <TParticle.h>
 #include <TString.h>
 #include <TClonesArray.h>
@@ -61,7 +64,9 @@ AliMCEventHandler::AliMCEventHandler() :
     fFileNumber(0),
     fEventsPerFile(0),
     fReadTR(kTRUE),
-    fInitOk(kFALSE)
+    fInitOk(kFALSE),
+    fSubsidiaryHandlers(0),
+    fEventsInContainer(0)
 {
   //
   // Default constructor
@@ -90,7 +95,9 @@ AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
     fFileNumber(0),
     fEventsPerFile(0),
     fReadTR(kTRUE),
-    fInitOk(kFALSE)
+    fInitOk(kFALSE),
+    fSubsidiaryHandlers(0),
+    fEventsInContainer(0)
 {
   //
   // Constructor
@@ -152,6 +159,17 @@ Bool_t AliMCEventHandler::Init(Option_t* opt)
     fFileNumber =  0;
     AliInfo(Form("Number of events in this directory %5d \n", fNEvent));
     fInitOk = kTRUE;
+
+
+    if (fSubsidiaryHandlers) {
+       TIter next(fSubsidiaryHandlers);
+       AliMCEventHandler *handler;
+       while((handler = (AliMCEventHandler*)next())) {
+           handler->Init(opt);
+           handler->SetNumberOfEventsInContainer(fNEvent);
+       }
+    }
+
     return kTRUE;
 }
 
@@ -193,6 +211,7 @@ Bool_t AliMCEventHandler::GetEvent(Int_t iev)
        // Connect TR to MCEvent
        fMCEvent->ConnectTreeTR(fTreeTR);
     }
+
     //
     return kTRUE;
 }
@@ -226,6 +245,7 @@ Bool_t AliMCEventHandler::OpenFile(Int_t i)
     }
     
     fInitOk = kTRUE;
+
     return kTRUE;
 }
 
@@ -234,6 +254,12 @@ Bool_t AliMCEventHandler::BeginEvent(Long64_t entry)
     fParticleSelected.Delete();
     fLabelMap.Delete();
     // Read the next event
+
+    if (fEventsInContainer != 0) {
+       entry = (Long64_t) entry * Float_t(fNEvent) / Float_t (fEventsInContainer);
+    }
+
+
     if (entry == -1) {
        fEvent++;
        entry = fEvent;
@@ -242,17 +268,34 @@ Bool_t AliMCEventHandler::BeginEvent(Long64_t entry)
     }
 
     if (entry >= fNEvent) {
-       AliWarning(Form("AliMCEventHandler: Event number out of range %5d %5d\n", entry,fNEvent));
+       AliWarning(Form("AliMCEventHandler: Event number out of range %5d %5d\n", entry, fNEvent));
        return kFALSE;
     }
-    return GetEvent(entry);
+    
+    Bool_t result = GetEvent(entry);
+
+    if (fSubsidiaryHandlers) {
+       TIter next(fSubsidiaryHandlers);
+       AliMCEventHandler *handler;
+       while((handler = (AliMCEventHandler*)next())) {
+           handler->BeginEvent(entry);
+       }
+       next.Reset();
+       while((handler = (AliMCEventHandler*)next())) {
+           fMCEvent->AddSubsidiaryEvent(handler->MCEvent());
+       }
+       fMCEvent->InitEvent();
+    }
+    return result;
+    
 }
 
 void AliMCEventHandler::SelectParticle(Int_t i){
   // taking the absolute values here, need to take care 
   // of negative daughter and mother
   // IDs when setting!
-  if(!IsParticleSelected(TMath::Abs(i)))fParticleSelected.Add(TMath::Abs(i),1);
+    if (TMath::Abs(i) >= AliMCEvent::BgLabelOffset()) i =  fMCEvent->BgLabelToIndex(TMath::Abs(i));
+    if(!IsParticleSelected(TMath::Abs(i)))fParticleSelected.Add(TMath::Abs(i),1);
 }
 
 Bool_t AliMCEventHandler::IsParticleSelected(Int_t i)  {
@@ -276,10 +319,9 @@ void AliMCEventHandler::CreateLabelMap(){
   }
 
   VerifySelectedParticles();
-  AliStack *pStack = fMCEvent->Stack();
 
   Int_t iNew = 0;
-  for(int i = 0;i < pStack->GetNtrack();++i){
+  for(int i = 0;i < fMCEvent->GetNumberOfTracks();++i){
     if(IsParticleSelected(i)){
       fLabelMap.Add(i,iNew);
       iNew++;
@@ -288,7 +330,7 @@ void AliMCEventHandler::CreateLabelMap(){
 }
 
 Int_t AliMCEventHandler::GetNewLabel(Int_t i) {
-  // Gets the labe from the new created Map
+  // Gets the label from the new created Map
   // Call CreatLabelMap before
   // otherwise only 0 returned
   return fLabelMap.GetValue(TMath::Abs(i));
@@ -302,28 +344,29 @@ void  AliMCEventHandler::VerifySelectedParticles(){
   // Private, should be only called by CreateLabelMap
 
   if(!fMCEvent){
-    fParticleSelected.Delete();
-    return;
+      fParticleSelected.Delete();
+      return;
   }
-  AliStack *pStack = fMCEvent->Stack();
 
-  Int_t nprim = pStack->GetNprimary();
-
-  for(int i = 0;i < pStack->GetNtrack();++i){
-    if(i<nprim){
-      SelectParticle(i);// take all primaries
-      continue;
-    }
-    if(!IsParticleSelected(i))continue;
-    TParticle *part = pStack->Particle(i);
-    Int_t imo = part->GetFirstMother();
-    while((imo >= nprim)&&!IsParticleSelected(imo)){
-      // Mother not yet selected
-      SelectParticle(imo);
-      TParticle *mother = pStack->Particle(imo);
-      imo = mother->GetFirstMother();
-    }
-    // after last step we may have a unselected primary
+  Int_t nprim = fMCEvent->GetNumberOfPrimaries();
+
+  for(int i = 0;i < fMCEvent->GetNumberOfTracks(); ++i){
+      if(i < nprim){
+         SelectParticle(i);// take all primaries
+         continue;
+      }
+
+      if(!IsParticleSelected(i))continue;
+
+      AliMCParticle* mcpart = fMCEvent->GetTrack(i);
+      Int_t imo = mcpart->GetMother();
+      while((imo >= nprim)&&!IsParticleSelected(imo)){
+         // Mother not yet selected
+         SelectParticle(imo);
+         AliMCParticle* mcpart = fMCEvent->GetTrack(imo);
+         imo = mcpart->GetMother();
+      }
+    // after last step we may have an unselected primary
     // mother
     if(imo>=0){
       if(!IsParticleSelected(imo))
@@ -367,12 +410,29 @@ Bool_t AliMCEventHandler::Notify(const char *path)
     else if (fileName.BeginsWith("root:")) {
       fileName.Append("?ZIP=");
     }
+
     *fPathName = fileName;
-    AliInfo(Form("Notify() Path: %s\n", fPathName->Data()));
+    AliInfo(Form("Path: -%s-\n", fPathName->Data()));
     
     ResetIO();
     InitIO("");
 
+// Handle subsidiary handlers
+    if (fSubsidiaryHandlers) {
+       TIter next(fSubsidiaryHandlers);
+       AliMCEventHandler *handler;
+       while((handler = (AliMCEventHandler*) next())) {
+           TString* spath = handler->GetInputPath();
+           if (spath->Contains("merged")) {
+               if (! fPathName->IsNull()) {
+                   handler->Notify(Form("%s/../.", fPathName->Data()));
+               } else {
+                   handler->Notify("../");
+               }
+           }
+       }
+    }
+    
     return kTRUE;
 }
 
@@ -391,6 +451,15 @@ void AliMCEventHandler::ResetIO()
     if (fFileTR) {delete fFileTR; fFileTR = 0;}
     fExtension="";
     fInitOk = kFALSE;
+
+    if (fSubsidiaryHandlers) {
+       TIter next(fSubsidiaryHandlers);
+       AliMCEventHandler *handler;
+       while((handler = (AliMCEventHandler*)next())) {
+           handler->ResetIO();
+       }
+    }
+
 }
 
                            
@@ -400,6 +469,15 @@ Bool_t AliMCEventHandler::FinishEvent()
     delete fDirTR;  fDirTR = 0;
     delete fDirK;   fDirK  = 0;    
     if (fInitOk) fMCEvent->FinishEvent();
+
+    if (fSubsidiaryHandlers) {
+       TIter next(fSubsidiaryHandlers);
+       AliMCEventHandler *handler;
+       while((handler = (AliMCEventHandler*)next())) {
+           handler->FinishEvent();
+       }
+    }
+
     return kTRUE;
 }
 
@@ -422,3 +500,11 @@ void AliMCEventHandler::SetInputPath(const char* fname)
     delete fPathName;
     fPathName = new TString(fname);
 }
+
+void AliMCEventHandler::AddSubsidiaryHandler(AliMCEventHandler* handler)
+{
+    // Add a subsidiary handler. For example for background events
+
+    if (!fSubsidiaryHandlers) fSubsidiaryHandlers = new TList();
+    fSubsidiaryHandlers->Add(handler);
+}
index b3e8f0bf4e7c12005bfa1ccaeeb6f1e15eab98fc..8fa8cc4d6d3fc1852d6072dbc4f904230a3961ae 100644 (file)
@@ -21,6 +21,8 @@
 
 class TFile;
 class TTree;
+class TList;
+
 class TParticle;
 class TString;
 class TClonesArray;
@@ -53,6 +55,8 @@ public:
     virtual void         ResetIO();
     virtual Bool_t       GetEvent(Int_t iev);
     virtual void         SetReadTR(Bool_t flag) { fReadTR = flag; }
+    virtual void         AddSubsidiaryHandler(AliMCEventHandler* handler);
+    virtual void         SetNumberOfEventsInContainer(Int_t nev) {fEventsInContainer = nev;}
     //
     AliMCEvent* MCEvent() const {return fMCEvent;} 
     TTree*      TreeTR()  const {return fTreeTR;}
@@ -72,25 +76,27 @@ private:
     AliMCEventHandler(const AliMCEventHandler& handler);             
     AliMCEventHandler& operator=(const AliMCEventHandler& handler);  
 private:
-    AliMCEvent       *fMCEvent;          //! MC Event
-    TFile            *fFileE;            //! File with TreeE
-    TFile            *fFileK;            //! File with TreeK
-    TFile            *fFileTR;           //! File with TreeTR
-    TTree            *fTreeE;            //! TreeE  (Event Headers)
-    TTree            *fTreeK;            //! TreeK  (kinematics tree)
-    TTree            *fTreeTR;           //! TreeTR (track references tree)
-    TDirectoryFile   *fDirK;             //! Directory for Kine Tree
-    TDirectoryFile   *fDirTR;            //! Directory for TR Tree
-    TExMap            fParticleSelected; //! List of selected MC particles for t
-    TExMap            fLabelMap;         //! Stores the Map of MC (ESDLabel,AODlabel)  
-    Int_t             fNEvent;           //! Number of events
-    Int_t             fEvent;            //! Current event
-    TString          *fPathName;         //! Input file path 
-    const Char_t     *fExtension;        //! File name extension 
-    Int_t             fFileNumber;       //! Input file number
-    Int_t             fEventsPerFile;    //! Number of events per file
-    Bool_t            fReadTR;           // determines if TR shall be read
-    Bool_t            fInitOk;           // Initialization ok
+    AliMCEvent            *fMCEvent;            //! MC Event
+    TFile                 *fFileE;              //! File with TreeE
+    TFile                 *fFileK;              //! File with TreeK
+    TFile                 *fFileTR;             //! File with TreeTR
+    TTree                 *fTreeE;              //! TreeE  (Event Headers)
+    TTree                 *fTreeK;              //! TreeK  (kinematics tree)
+    TTree                 *fTreeTR;             //! TreeTR (track references tree)
+    TDirectoryFile        *fDirK;               //! Directory for Kine Tree
+    TDirectoryFile        *fDirTR;              //! Directory for TR Tree
+    TExMap                 fParticleSelected;   //! List of selected MC particles for t
+    TExMap                 fLabelMap;           //! Stores the Map of MC (ESDLabel,AODlabel)  
+    Int_t                  fNEvent;             //! Number of events
+    Int_t                  fEvent;              //! Current event
+    TString               *fPathName;           //! Input file path 
+    const Char_t          *fExtension;          //! File name extension 
+    Int_t                  fFileNumber;         //! Input file number
+    Int_t                  fEventsPerFile;      //! Number of events per file
+    Bool_t                 fReadTR;             // determines if TR shall be read
+    Bool_t                 fInitOk;             // Initialization ok
+    TList                 *fSubsidiaryHandlers; //! List of subsidiary MC handlers (for example for Background)
+    Int_t                  fEventsInContainer;  //! Number of events in container class 
     ClassDef(AliMCEventHandler,1)  //MC Truth EventHandler class
 };
 #endif 
index e378eebfa7170d4e779983ce9cd5beccbc813dcd..aff6710d20d077e5784fbfefb6a5a81c2c987145 100644 (file)
@@ -34,7 +34,10 @@ AliMCParticle::AliMCParticle():
     fParticle(0),
     fTrackReferences(0),
     fNTrackRef(0),
-    fLabel(-1)
+    fLabel(-1),
+    fMother(-1),
+    fFirstDaughter(-1),
+    fLastDaughter(-1)
 {
     // Constructor
 }
@@ -45,7 +48,10 @@ AliMCParticle::AliMCParticle(TParticle* part, TRefArray* rarray, Int_t index):
     fParticle(part),
     fTrackReferences(rarray),
     fNTrackRef(0),
-    fLabel(index)
+    fLabel(index),
+    fMother(-1),
+    fFirstDaughter(-1),
+    fLastDaughter(-1)
 {
     // Constructor
     if (rarray != 0) {
@@ -59,7 +65,10 @@ AliMCParticle::AliMCParticle(const AliMCParticle& mcPart) :
     fParticle(0),    
     fTrackReferences(0),
     fNTrackRef(0),
-    fLabel(-1)
+    fLabel(-1),
+    fMother(-1),
+    fFirstDaughter(-1),
+    fLastDaughter(-1)
 {
 // Copy constructor
 }
index bca1ac3efe889184b8b751ec0d90938b9adac6c6..75d4a266f8639329bb10f6555792646f1019ce27 100644 (file)
@@ -66,12 +66,23 @@ public:
 
     // "Trackable" criteria
     Float_t  GetTPCTrackLength(Float_t bz, Float_t ptmin, Int_t &counter, Float_t deadWidth);
-    
+    // Navigation
+    Int_t GetMother()        const {return fMother;}
+    Int_t GetFirstDaughter() const {return fFirstDaughter;}
+    Int_t GetLastDaughter()  const {return fLastDaughter;}
+    void  SetMother(Int_t idx)        {fMother        = idx;}
+    void  SetFirstDaughter(Int_t idx) {fFirstDaughter = idx;}
+    void  SetLastDaughter(Int_t idx)  {fLastDaughter  = idx;}
+    void  SetLabel(Int_t label)       {fLabel         = label;}
+           
  private:
     TParticle *fParticle;             // The wrapped TParticle
     TRefArray *fTrackReferences;      // Reference array to track references
     Int_t      fNTrackRef;            // Number of track references
     Int_t      fLabel;                // fParticle Label in the Stack
+    Int_t      fMother;               // Mother particles
+    Int_t      fFirstDaughter;        // First daughter
+    Int_t      fLastDaughter;         // LastDaughter
     
   ClassDef(AliMCParticle,0)  // AliVParticle realisation for MCParticles
 };