Added an output file for average cross section and number of trials from pyxsec.root...
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 1 Sep 2009 13:29:59 +0000 (13:29 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 1 Sep 2009 13:29:59 +0000 (13:29 +0000)
ANALYSIS/AliAnalysisTaskMCParticleFilter.cxx
ANALYSIS/AliAnalysisTaskMCParticleFilter.h
ANALYSIS/macros/AddTaskESDFilter.C
ANALYSIS/runAODFilterMC.C

index 0aed4ab..5d5233a 100644 (file)
 #include <TChain.h>
 #include <TFile.h>
 #include "TParticle.h"
+#include "TParameter.h"
 #include "TString.h"
+#include "TList.h"
+#include "TProfile.h"
+#include "TH1F.h"
+
 
 #include "AliAnalysisTaskMCParticleFilter.h"
 #include "AliAnalysisManager.h"
@@ -47,20 +52,75 @@ ClassImp(AliAnalysisTaskMCParticleFilter)
 //____________________________________________________________________
 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
 AliAnalysisTaskSE(),
-  fTrackFilterMother(0x0)
+  fTrackFilterMother(0x0),
+  fHistList(0x0)
 {
   // Default constructor
 }
 
+Bool_t AliAnalysisTaskMCParticleFilter::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // 
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Double_t xsection = 0;
+  UInt_t   ntrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "pyxsec.root");
+    }
+
+
+    TFile *fxsec = TFile::Open(fileName.Data());
+    if(!fxsec){
+      AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
+      // not a severe condition we just do not have the information...
+      return kTRUE;
+    }
+    TTree *xtree = (TTree*)fxsec->Get("Xsection");
+    if(!xtree){
+      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+      return kTRUE;
+    }
+    xtree->SetBranchAddress("xsection",&xsection);
+    xtree->SetBranchAddress("ntrials",&ntrials);
+    xtree->GetEntry(0);
+    ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
+    ((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials); 
+
+  }
+  return kTRUE;
+}
+
+
+
 //____________________________________________________________________
 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
     AliAnalysisTaskSE(name),
-    fTrackFilterMother(0x0)
+    fTrackFilterMother(0x0),
+    fHistList(0x0)
 {
   // Default constructor
+  DefineOutput(1,TList::Class());  
 }
 
-
+/*
 //____________________________________________________________________
 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
     AliAnalysisTaskSE(obj),
@@ -68,14 +128,14 @@ AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalys
 {
   // Copy constructor
 }
-
+*/
 //____________________________________________________________________
 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
 {
   //  if( fTrackFilterMother ) delete fTrackFilterMother;
 }
 
-
+/*
 //____________________________________________________________________
 AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
 {
@@ -86,7 +146,7 @@ AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(cons
   }
   return *this;
 }
-
+*/
 //____________________________________________________________________
 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
 {
@@ -123,6 +183,22 @@ void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
     }
     aodH->SetMCEventHandler(mcH);
 
+
+    // these are histograms, for averaging and summing
+    // do it via histograms to be PROOF-proof
+    // which merges the results from different workers
+    // histos are not saved reside only in memory
+    
+    
+
+    fHistList = new TList();
+    TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
+    h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+    fHistList->Add(h1Xsec);
+    TH1F* h1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
+    h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+    fHistList->Add(h1Trials);
+
     
 }
 
@@ -251,7 +327,7 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
       Int_t d1 =  mcpart->GetLastDaughter();
       if(d0>0&&d1>0){
        for(int id = d0;id <= d1;id++){
-         AliMCParticle* daughter =  mcE->GetTrack(id);
+         // AliMCParticle* daughter =  mcE->GetTrack(id);
          iAll++;
          if(mcH->IsParticleSelected(id))iTaken++;
        }
@@ -261,10 +337,50 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
 
   AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
 
+  PostData(1,fHistList);
 
   return;
 }
 
+void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
+  //
+  // Terminate the execution save the PYTHIA x_section to the UserInfo()
+  //
+
+
+  TTree *tree = (TTree*)GetOutputData(0);
+  fHistList = (TList*)GetOutputData(1);
+  if(!fHistList){
+    Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
+    return;
+  }
+
+  if(!tree){
+    Printf("%s:%d Output tree not found",(char*)__FILE__,__LINE__);
+    return;
+  }
+
+  TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
+  TH1F* hTrials  = ((TH1F*)(fHistList->FindObject("h1Trials")));
+  if(!hXsec)return;
+  if(!hTrials)return;
+
+  Float_t xsec = hXsec->GetBinContent(1);
+  Float_t trials = hTrials->Integral();
+  AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
+  
+  // This only works if the tree is not a special output
+  // in this case it is already written
+  /*
+  TParameter<float> *x = new TParameter<float>("xsection",xsec);
+  tree->GetUserInfo()->Add(x);
+  TParameter<float> *t = new TParameter<float>("trials",trials);
+  tree->GetUserInfo()->Add(t);
+  */
+  
+
+}
+
 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
 {
   // Selection accoring to eta of the mother and production point
index 4a9442a..31f5b3c 100644 (file)
@@ -9,34 +9,40 @@
 //  Fill AOD tracks from Kinematic stack
 //
 
-#include <TList.h> 
 #include "AliAnalysisTaskSE.h"
 
 class AliAnalysisFilter;
 class TString;
+class TList;
 
 class AliAnalysisTaskMCParticleFilter : public AliAnalysisTaskSE
 {
  public:
                                   AliAnalysisTaskMCParticleFilter();
                                   AliAnalysisTaskMCParticleFilter( const char* name );
-                                  AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj);
     virtual                      ~AliAnalysisTaskMCParticleFilter();
-     AliAnalysisTaskMCParticleFilter&   operator=(const AliAnalysisTaskMCParticleFilter& other);
     
     // Implementation of interface methods
     virtual                void   UserCreateOutputObjects();
     virtual                void   UserExec( Option_t *option );
-    
+    virtual                Bool_t Notify();
+    virtual                void   Terminate( Option_t *option );
     // Setters
     virtual                void   SetTrackFilterMother(AliAnalysisFilter* trackF) { fTrackFilterMother = trackF; }
     
  private:
     Bool_t Select(TParticle* part, Float_t rv, Float_t zv);                
+    
+    // pivate c'tors to prevent misuse
+    AliAnalysisTaskMCParticleFilter&   operator=(const AliAnalysisTaskMCParticleFilter& other);
+    AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj);
+
+
 
     AliAnalysisFilter*  fTrackFilterMother;   //  Track Filter
-                
-    ClassDef( AliAnalysisTaskMCParticleFilter, 1 ); // Analysis task for Kinematic filtering
+    TList *fHistList;                         // list to store e histograms, only as exchange
+
+    ClassDef( AliAnalysisTaskMCParticleFilter, 2 ); // Analysis task for Kinematic filtering
 };
  
 #endif
index bd8ddc4..2c0de06 100644 (file)
@@ -91,6 +91,8 @@ AliAnalysisTaskESDfilter *AddTaskESDFilter(Bool_t useKineFilter=kTRUE)
    if (useKineFilter) {
       mgr->ConnectInput  (kinefilter,  0, mgr->GetCommonInputContainer());
       mgr->ConnectOutput (kinefilter,  0, mgr->GetCommonOutputContainer());
+      AliAnalysisDataContainer *coutputEx = mgr->CreateContainer("cFilterList", TList::Class(),
+                                                                  AliAnalysisManager::kOutputContainer,"pyxsec_hists.root");
    }   
    return esdfilter;
 }   
index 513482c..20a686d 100644 (file)
@@ -96,6 +96,9 @@ void runAODFilterMC()
       if(bKineFilter){
        mgr->ConnectInput  (kinefilter,     0, cinput1  );
        mgr->ConnectOutput (kinefilter,     0, coutput1 );
+       AliAnalysisDataContainer *coutputEx = mgr->CreateContainer("cFilterList", TList::Class(),
+                                                                  AliAnalysisManager::kOutputContainer,"pyxsec_hists.root");
+       mgr->ConnectOutput (kinefilter,     1, coutputEx );
       }
 
        mgr->ConnectInput  (esdfilter,     0, cinput1  );