Updates for associating particles to generators
authormorsch <andreas.morsch@cern.ch>
Thu, 3 Jul 2014 18:16:17 +0000 (20:16 +0200)
committermorsch <andreas.morsch@cern.ch>
Thu, 3 Jul 2014 18:16:17 +0000 (20:16 +0200)
STEER/AOD/AliAODInputHandler.cxx
STEER/AOD/AliAODMCParticle.h
STEER/STEERBase/AliMCEvent.cxx
STEER/STEERBase/AliMCEvent.h
STEER/STEERBase/AliVParticle.h

index f5112f0..501ec8f 100644 (file)
@@ -132,15 +132,17 @@ Bool_t AliAODInputHandler::BeginEvent(Long64_t entry)
       fEvent->InitMagneticField();
       prevRunNumber = fEvent->GetRunNumber();
     } 
+
+    AliAODMCHeader* mcHeader = (AliAODMCHeader*) fEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+      if(mcHeader) {
+      fMCEvent->SetExternalHeader(mcHeader);}
+
     TClonesArray* mcParticles = (TClonesArray*) (fEvent->FindListObject("mcparticles"));
     if (mcParticles) {
        if (!fMCEvent) fMCEvent = new AliMCEvent();
        fMCEvent->SetParticleArray(mcParticles);
     }
 
-    AliAODMCHeader* mcHeader = (AliAODMCHeader*) fEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
-      if(mcHeader) {
-      fMCEvent->SetExternalHeader(mcHeader);}
     // When merging, get current event number from GetReadEntry(), 
     // entry gives the events in the current file
     if (fTreeToMerge) fTreeToMerge->GetEntry(GetReadEntry() + fMergeOffset);
index 15820ce..2a2d676 100644 (file)
@@ -101,7 +101,7 @@ class AliAODMCParticle: public AliVParticle {
       if(b)fFlag |= kPrimary;
       else fFlag &= ~kPrimary;
     }
-    Bool_t IsPrimary() const {return ((fFlag&kPrimary)==kPrimary);} 
+    virtual Bool_t IsPrimary() const {return ((fFlag&kPrimary)==kPrimary);} 
 
     void SetPhysicalPrimary(Bool_t b = kTRUE){
      if(b)fFlag |= kPhysicalPrim;
index d1b6787..4f55a3d 100644 (file)
@@ -793,6 +793,7 @@ void AliMCEvent::PreReadAll()
     {
        GetTrack(i);
     }
+    AssignGeneratorIndex();
 }
 
 const AliVVertex * AliMCEvent::GetPrimaryVertex() const 
@@ -826,38 +827,42 @@ Bool_t AliMCEvent::IsFromBGEvent(Int_t index)
 }
 
 
-    Int_t AliMCEvent::GetCocktailList(TList*& lh)
+    TList* AliMCEvent::GetCocktailList()
     {
       //gives the CocktailHeaders when reading ESDs/AODs (corresponding to fExteral=kFALSE/kTRUE)
       //the AODMC header (and the aodmc array) is passed as an instance to MCEvent by the AliAODInputHandler
       if(fExternal==kFALSE) { 
        AliGenCocktailEventHeader* coHeader =dynamic_cast<AliGenCocktailEventHeader*> (GenEventHeader());
-       if(!coHeader) return 0;
-       lh=coHeader->GetHeaders();}
-      if(fExternal==kTRUE){ 
-       if(!fAODMCHeader) return 0;
-       lh=fAODMCHeader->GetCocktailHeaders();}
-      return 1;  
+       if(!coHeader) {
+         return 0;
+       } else {
+         return (coHeader->GetHeaders());
+       }
+      } else {
+       if(!fAODMCHeader) { 
+         return 0;
+       } else {
+         return (fAODMCHeader->GetCocktailHeaders());
+       }
+      }
     }
 
 
-
-
 TString AliMCEvent::GetGenerator(Int_t index)
 {
-  Int_t nsumpart=0;
-  
-  TList* lh;
-  Int_t nt= GetCocktailList(lh);
-  if(nt==0){ TString noheader="nococktailheader";
+  Int_t nsumpart=fNprimaries;
+  TList* lh = GetCocktailList();
+  if(!lh){ TString noheader="nococktailheader";
     return noheader;}
   Int_t nh=lh->GetEntries();
-  for(Int_t i=0;i<nh;i++){
+  for (Int_t i = nh-1; i >= 0; i--){
     AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
-     TString genname=gh->GetName();
-     Int_t npart=gh->NProduced();
-     if(index>=nsumpart && index<(nsumpart+npart)) return genname;
-     nsumpart+=npart;}
+    TString genname=gh->GetName();
+    Int_t npart=gh->NProduced();
+    if (i == 0) npart = nsumpart;
+    if(index < nsumpart && index >= (nsumpart-npart)) return genname;
+    nsumpart-=npart;
+  }
   TString empty="";
   return empty;
 }
@@ -866,34 +871,56 @@ void AliMCEvent::AssignGeneratorIndex() {
   //
   // Assign the generator index to each particle
   //
-  TList* list;
-  Int_t nt = GetCocktailList(list);
-  if (nt == 0) {
+  TList* list = GetCocktailList();
+  if (!list) {
+    return;
   } else {
     Int_t nh = list->GetEntries();
-    Int_t nsumpart = 0;
-    for(Int_t i = 0; i < nh; i++){
+    Int_t nsumpart = fNprimaries;
+    for(Int_t i = nh-1; i >= 0; i--){
       AliGenEventHeader* gh = (AliGenEventHeader*)list->At(i);
       Int_t npart = gh->NProduced();
-      for (Int_t j = nsumpart; j < npart; j++) {
+      if (i==0) npart = nsumpart;
+      //
+      // Loop over primary particles for generator i
+      for (Int_t j = nsumpart-1; j >= nsumpart-npart; j--) {
        AliVParticle* part = GetTrack(j);
        part->SetGeneratorIndex(i);
        Int_t dmin = part->GetFirstDaughter();
        Int_t dmax = part->GetLastDaughter();
-       for (Int_t k = dmin; k <= dmax; k++) {
-         AliVParticle* dpart = GetTrack(k);
-         dpart->SetGeneratorIndex(i);
-       }
+       if (dmin == -1) continue;
+       AssignGeneratorIndex(i, dmin, dmax);
       } 
-      nsumpart += npart;
+      nsumpart -= npart;
+    }
+  }
+}
+void AliMCEvent::AssignGeneratorIndex(Int_t index, Int_t dmin, Int_t dmax) {
+  for (Int_t k = dmin; k <= dmax; k++) {
+    AliVParticle* dpart = GetTrack(k);
+    dpart->SetGeneratorIndex(index);
+    Int_t dmin = dpart->GetFirstDaughter();
+    Int_t dmax = dpart->GetLastDaughter();
+    if (dmin > -1) {
+      AssignGeneratorIndex(index, dmin, dmax);
     }
   }
 }
-
 
    Bool_t  AliMCEvent::GetCocktailGenerator(Int_t index,TString &nameGen){
-  //method that gives the generator for a given particle with label index (or that of the corresponding primary)
-  
+     //method that gives the generator for a given particle with label index (or that of the corresponding primary)
+     AliVParticle* mcpart0 = (AliVParticle*) (GetTrack(index));
+     if(!mcpart0){
+       printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",index);
+       return 0;
+     }
+     /*
+     Int_t ig = mcpart0->GetGeneratorIndex();
+     if (ig != -1) {
+       nameGen = ((AliGenEventHeader*)GetCocktailList()->At(ig))->GetName();
+       return 1;
+     }
+     */
     nameGen=GetGenerator(index);
     if(nameGen.Contains("nococktailheader") )return 0;
     Int_t lab=index;
@@ -926,6 +953,16 @@ void AliMCEvent::AssignGeneratorIndex() {
    return 1;
    }
 
+void  AliMCEvent::SetParticleArray(TClonesArray* mcParticles) 
+  {
+    fMCParticles = mcParticles; 
+    fNparticles = fMCParticles->GetEntries(); 
+    fExternal = kTRUE; 
+    fNprimaries = 0;
+    for (Int_t i = 0; i < mcParticles->GetEntries(); i++) 
+      if  (((AliVParticle*) mcParticles->At(i))->IsPrimary()) fNprimaries++;
+    AssignGeneratorIndex();
+  }
 
 
 
index 15d3571..6d1a077 100644 (file)
@@ -129,19 +129,17 @@ public:
     virtual Int_t     BgLabelToIndex(Int_t label);
     static  Int_t     BgLabelOffset() {return fgkBgLabelOffset;}
     virtual Bool_t    IsFromBGEvent(Int_t index);
-    Int_t  GetCocktailList(TList*& lista);
+    TList*  GetCocktailList();
     TString  GetGenerator(Int_t index); 
     Bool_t GetCocktailGenerator(Int_t index,TString &nameGen);
     virtual Bool_t    IsSecondaryFromWeakDecay(Int_t index);
     virtual Bool_t    IsSecondaryFromMaterial(Int_t index);
     // External particle array
-    virtual void      SetParticleArray(TClonesArray* mcParticles) 
-       {fMCParticles = mcParticles; fNparticles = fMCParticles->GetEntries(); fExternal = kTRUE;}
+  virtual void      SetParticleArray(TClonesArray* mcParticles); 
     //External Header 
      virtual void SetExternalHeader(AliVHeader* aodmcHeader)
        {fAODMCHeader=aodmcHeader;}  
   virtual AliGenEventHeader *FindHeader(Int_t ipart);
-  virtual void AssignGeneratorIndex();    
     //Following needed only for mixed event
   virtual Int_t        EventIndex(Int_t)       const {return 0;}
   virtual Int_t        EventIndexForCaloCluster(Int_t) const {return 0;}
@@ -158,6 +156,8 @@ private:
     virtual void      ReorderAndExpandTreeTR();
     virtual Int_t     FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const;
     void             UpdateEventInformation();
+    virtual void      AssignGeneratorIndex();    
+    virtual void      AssignGeneratorIndex(Int_t index, Int_t dmin, Int_t dmax);    
     
 private: 
     // Stanndard implementation for ESD production
index 4298d6e..3f488c6 100644 (file)
@@ -74,6 +74,7 @@ public:
    *  @return     always kTRUE;
    */
   Bool_t IsSortable() const  { return kTRUE; }
+  virtual Bool_t IsPrimary()  const  { return kFALSE; }
   
 
   // coordinate system conversions