]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx
end-of-line normalization
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskMCParticleFilter.cxx
index 354b3a050f9dce9c6f397329d71b63f72d1904f3..7f6bae2200d94548dfdb49a254ff58e2f26d8fc5 100644 (file)
@@ -35,6 +35,7 @@
 #include "AliStack.h"
 #include "AliMCEvent.h"
 #include "AliMCEventHandler.h"
+#include "AliESDInputHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODHeader.h"
 #include "AliAODMCHeader.h"
@@ -46,6 +47,7 @@
 #include "AliGenHijingEventHeader.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliGenCocktailEventHeader.h"
+#include "AliGenEventHeaderTunedPbPb.h"
 
 #include "AliLog.h"
 
@@ -70,7 +72,8 @@ Bool_t AliAnalysisTaskMCParticleFilter::Notify()
   // Implemented Notify() to read the cross sections
   // from pyxsec.root
   // 
-  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  TTree *tree = mgr->GetTree();
   Double_t xsection = 0;
   UInt_t   ntrials  = 0;
   if(tree){
@@ -81,7 +84,11 @@ Bool_t AliAnalysisTaskMCParticleFilter::Notify()
     }
 
     TString fileName(curfile->GetName());
-    if(fileName.Contains("AliESDs.root")){
+    TString datafile = mgr->GetInputEventHandler()->GetInputFileName();
+    if (fileName.Contains(datafile)) {
+        fileName.ReplaceAll(datafile, "pyxsec.root");
+    }
+    else if(fileName.Contains("AliESDs.root")){
         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
     }
     else if(fileName.Contains("AliAOD.root")){
@@ -164,7 +171,7 @@ AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(cons
 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
 {
   // Create the output container
-  PostData(1,fHistList);
+  
 
     if (OutputTree()&&fTrackFilterMother) 
        OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
@@ -212,17 +219,19 @@ void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
     h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
     fHistList->Add(h1Trials);
 
+    PostData(1,fHistList);
     
 }
 
 //____________________________________________________________________
 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
 {
-// Execute analysis for current event
-//
-
-// Fill AOD tracks from Kinematic stack
-    
+  // Execute analysis for current event
+  //
+  
+  // Fill AOD tracks from Kinematic stack
+  PostData(1,fHistList);
+  
   // get AliAOD Event 
   AliAODEvent* aod = AODEvent();
   if (!aod) {
@@ -269,12 +278,17 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
   AliGenHijingEventHeader *hiH  = 0;
   AliCollisionGeometry    *colG = 0;
   AliGenDPMjetEventHeader *dpmH = 0;
+  AliGenEventHeaderTunedPbPb *tunedH = 0;
+
   // it can be only one save some casts
   // assuming PYTHIA and HIJING are the most likely ones...
   if(!pyH){
       hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
       if(!hiH){
          dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
+         if(!dpmH){
+           tunedH = dynamic_cast<AliGenEventHeaderTunedPbPb*>(mcEH);
+         }
       }
   }
   
@@ -290,10 +304,10 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
     if (ccEH) {
       TList *genHeaders = ccEH->GetHeaders();
       for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
-       if(!pyH)dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
-       if(!hiH)dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
+       if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
+       if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
        if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
-       if(!dpmH)dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
+       if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
       }
     }
   }
@@ -311,13 +325,10 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
     fAODMcHeader->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
     AliInfo(Form("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle()));
   }
-
-
-
-  // check varibales for charm need all daughters
-  static int  iTaken = 0;
-  static int  iAll = 0;
-  static int  iCharm = 0;
+  else if(tunedH) {
+    fAODMcHeader->SetReactionPlaneAngle(tunedH->GetPsi2());
+    fAODMcHeader->SetCrossSection(tunedH->GetCentrality());
+  }
 
 
   Int_t j=0;
@@ -344,7 +355,8 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
        imo =  mother->GetMother();
       }
       // Select according to pseudorapidity and production point of primary ancestor
-      if (imo < nprim && Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kTRUE;         
+      if (imo < nprim)write = kTRUE;         
+      // if(!Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kFALSE; // selection on eta and phi of the mother
     } else if (part->GetUniqueID() == kPPair) {
       // Now look for pair production
       Int_t imo = mcpart->GetMother();
@@ -364,44 +376,49 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
          write = kTRUE;
       }
     }
-    /*
-    else if (part->GetUniqueID() == 13){
-      // Evaporation
-      // Check that we end at a primary particle
-      TParticle* mother = part;
-      Int_t imo = part->GetFirstMother();
-      while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
-       mother =  mcE->Stack()->Particle(imo);
-       imo =  mother->GetFirstMother();
-      }
-      // Select according to pseudorapidity and production point 
-       if (imo < nprim && Select(mother, rv, zv)) 
-         write = kTRUE;
-    }
-    */    
+
     if (write) {
       if(mcH)mcH->SelectParticle(ip);
-      j++;
-      
-      // debug info to check fro charm daugthers
-      if((TMath::Abs(part->GetPdgCode()))==411){
-       iCharm++;
-       Int_t d0 =  mcpart->GetFirstDaughter();
-       Int_t d1 =  mcpart->GetLastDaughter();
-       if(d0>0&&d1>0){
-         for(int id = d0;id <= d1;id++){
-           iAll++;
-           if(mcH->IsParticleSelected(id))iTaken++;
+      j++;      
+    }
+  }
+
+  /*
+
+  // check varibales for charm need all daughters
+  static int  iTaken = 0;
+  static int  iAll = 0;
+  static int  iCharm = 0;
+
+  for (Int_t ip = 0; ip < np; ip++){
+    AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
+    TParticle* part = mcpart->Particle();
+    // debug info to check fro charm daugthers
+    if((TMath::Abs(part->GetPdgCode()))==411){
+      iCharm++;
+      Int_t d0 =  part->GetFirstDaughter();
+      Int_t d1 =  part->GetLastDaughter();
+      if(d0>0&&d1>0){
+       for(int id = d0;id <= d1;id++){
+         iAll++;
+         if(mcH->IsParticleSelected(id)){
+           iTaken++;
+           Printf("> %d/%d Taken",id,nprim);
+           PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
+         }
+         else{
+           Printf("> %d/%d NOT Taken",id,nprim);
+           PrintMCParticle( (AliMCParticle*)mcE->GetTrack(id),id);
          }
        }
-      }// if charm
-    }
+      }
+    }// if charm
+    AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
   }
+  */
 
-  AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
 
   aodH->StoreMCParticles();
-  PostData(1,fHistList);
 
   return;
 }
@@ -461,3 +478,14 @@ Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Floa
  
 }
 
+void AliAnalysisTaskMCParticleFilter::PrintMCParticle(const AliMCParticle *mcp,Int_t np){
+  
+  const TParticle *p = mcp->Particle();
+  Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughte\
+r2 %d ",
+        np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter \
+        (0),p->GetDaughter(1));
+  Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
+  p->Print();
+  Printf("---------------------------------------");
+}