PWG2 resonances wagon included. Still problems so disabled for the moment.

[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskMCParticleFilter.cxx

diff --git a/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx b/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx
index 870452f9ddd1d86e4eba89b3167a9870d4c54a89..e2ca73deb1fa2acda327ef120dcc4128624e7272 100644 (file)
--- a/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx
+++ b/ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx
@@ -32,6 +32,7 @@
 #include "AliMCEventHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODHeader.h"
+#include "AliAODMCHeader.h"
 #include "AliAODHandler.h"
 #include "AliAODVertex.h"
 #include "AliAODMCParticle.h"
@@ -92,20 +93,38 @@ void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
   // Create the output container
     if (OutputTree()&&fTrackFilterMother) 
        OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
-    // how is this is reset cleared in the UserExec....
-    // Can this be handled by the framework?
+
+    // this part is mainly needed to set the MCEventHandler
+    // to the AODHandler, this will not be needed when
+    // AODHandler goes to ANALYSISalice
+    // setting in the steering macro will not work on proof :(
+    // so we have to do it in a task
+
+    // the branch booking can also go into the AODHandler then
+
+
+    // mcparticles
     TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
     tca->SetName(AliAODMCParticle::StdBranchName());
     AddAODBranch("TClonesArray",&tca);
-    AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
 
+    // MC header...
+    AliAODMCHeader *mcHeader = new AliAODMCHeader();
+    Printf("AODMCHeader %p",mcHeader);
+    Printf("AODMCHeader ** %p",&mcHeader);
+    mcHeader->SetName(AliAODMCHeader::StdBranchName());
+    AddAODBranch("AliAODMCHeader",&mcHeader);    
+
+    
+
+    AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
     AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
     if(!aodH){
       Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
       return;
     }
     aodH->SetMCEventHandler(mcH);
-    // TODO ADD MC VERTEX
+
     
 }
 
@@ -153,21 +172,17 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
   // We take all real primaries
 
 
-           
   Int_t j=0;
   for (Int_t ip = 0; ip < np; ip++){
     AliMCParticle* mcpart =  mcE->GetTrack(ip);
     TParticle* part = mcpart->Particle();
-    
     Float_t xv = part->Vx();
     Float_t yv = part->Vy();
     Float_t zv = part->Vz();
     Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
       
     Bool_t write = kFALSE;
-    Int_t flag = 0;
-
-
+    
     if (ip < nprim) {
       // Select the primary event
       write = kTRUE;
@@ -180,12 +195,10 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
        mother =  mcE->Stack()->Particle(imo);
        imo =  mother->GetFirstMother();
       }
-      // Select according to pseudorapidity and production point 
-       if (imo < nprim && Select(mother, rv, zv)) 
-         write = kTRUE;
+      // Select according to pseudorapidity and production point of primary ancestor
+      if (imo < nprim && Select(mcE->Stack()->Particle(imo), rv, zv))write = kTRUE;         
     } else if (part->GetUniqueID() == 5) {
       // Now look for pair production
-      continue;
       Int_t imo = part->GetFirstMother();
       if (imo < nprim) {
        // Select, if the gamma is a primary
@@ -224,6 +237,35 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
     }
   }
 
+
+    // check for charm daughters
+  static int  iTaken = 0;
+  static int  iAll = 0;
+  static int  iCharm = 0;
+  for (Int_t ip = 0; ip < np; ip++){
+    AliMCParticle* mcpart =  mcE->GetTrack(ip);
+    TParticle* part = mcpart->Particle();
+
+    //    if((TMath::Abs(part->GetPdgCode())/400)==1){
+    if((TMath::Abs(part->GetPdgCode()))==411){
+      // cases 
+      iCharm++;
+      Printf("Decay Mother %s",part->GetPDG()->GetName());
+      Int_t d0 =  part->GetFirstDaughter();
+      Int_t d1 =  part->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());
+         iAll++;
+         if(mcH->IsParticleSelected(id))iTaken++;
+       }
+      }
+    }
+  }
+  Printf("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm);
+
+
   return;
 }
 
@@ -242,8 +284,12 @@ Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Floa
   // larger for V0s in the TPC
   //  if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
 
+  if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;   
+
   // Andreas' Cuts
-  if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;   
+  //  if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;   
+
+
 
   // Muon arm
   if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms