]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
add calculation and histograms for MC cross section
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 15 Jun 2013 11:17:31 +0000 (11:17 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 15 Jun 2013 11:17:31 +0000 (11:17 +0000)
PWG/CaloTrackCorrBase/AliAnalysisTaskCounter.cxx
PWG/CaloTrackCorrBase/AliAnalysisTaskCounter.h

index abfc4f3caf32b601d9775c322647a5208ed75c51..2b4860577616a4100c91414efe24b8c740d1e82f 100644 (file)
 // Author: Gustavo Conesa Balbastre (LPSC)
 //         
 //_________________________________________________________________________
-
-#include "TH2F.h"
+#include <TSystem.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TList.h>
+#include <TClonesArray.h>
 #include "AliAODHeader.h"
 #include "AliTriggerAnalysis.h"
 #include "AliESDEvent.h"
@@ -57,14 +63,17 @@ AliAnalysisTaskCounter::AliAnalysisTaskCounter(const char *name)
   fAcceptFastCluster(kTRUE),
   fZVertexCut(10.), 
   fTrackMultEtaCut(0.8),
+  fAvgTrials(-1),
   fCaloFilterPatch(kFALSE),
   fOutputContainer(0x0), 
   fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
   fTriggerAnalysis (new AliTriggerAnalysis),
-  fhNEvents(0),
+  fCurrFileName(0),
+  fhNEvents(0),   
   fhXVertex(0),    fhYVertex(0),    fhZVertex(0),
   fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0),
-  fhCentrality(0), fhEventPlaneAngle(0)
+  fhCentrality(0), fhEventPlaneAngle(0),
+  fh1Xsec(0),      fh1Trials(0)
 {
   //ctor
   DefineOutput(1, TList::Class());
@@ -76,14 +85,17 @@ AliAnalysisTaskCounter::AliAnalysisTaskCounter()
     fAcceptFastCluster(kTRUE),
     fZVertexCut(10.),
     fTrackMultEtaCut(0.8),
+    fAvgTrials(-1),
     fCaloFilterPatch(kFALSE),
     fOutputContainer(0x0), 
     fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
     fTriggerAnalysis (new AliTriggerAnalysis),
+    fCurrFileName(0),
     fhNEvents(0),
     fhXVertex(0),    fhYVertex(0),    fhZVertex(0),
     fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0),
-    fhCentrality(0), fhEventPlaneAngle(0)
+    fhCentrality(0), fhEventPlaneAngle(0),
+    fh1Xsec(0),      fh1Trials(0)
 {
   // ctor
   DefineOutput(1, TList::Class());
@@ -115,6 +127,14 @@ void AliAnalysisTaskCounter::UserCreateOutputObjects()
   
   fOutputContainer = new TList();
   
+  fh1Xsec = new TH1F("hXsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+  fOutputContainer->Add(fh1Xsec);
+  
+  fh1Trials = new TH1F("hTrials","trials root file",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+  fOutputContainer->Add(fh1Trials);
+
   fhZVertex     = new TH1F("hZVertex", " Z vertex distribution"   , 200 , -50 , 50  ) ;
   fhZVertex->SetXTitle("v_{z} (cm)");
   fOutputContainer->Add(fhZVertex);
@@ -181,6 +201,7 @@ void AliAnalysisTaskCounter::UserCreateOutputObjects()
   
 }
 
+
 //_______________________________________________
 void AliAnalysisTaskCounter::UserExec(Option_t *) 
 {
@@ -189,6 +210,8 @@ void AliAnalysisTaskCounter::UserExec(Option_t *)
   
   //printf("___ Event __ %d __\n",(Int_t)Entry());
   
+  Notify();
+  
   fhNEvents->Fill(0.5);  
   
   AliVEvent * event = InputEvent();
@@ -419,6 +442,7 @@ void AliAnalysisTaskCounter::UserExec(Option_t *)
 
 }
 
+
 //____________________________________________________
 Bool_t AliAnalysisTaskCounter::CheckForPrimaryVertex()
 {
@@ -454,7 +478,6 @@ Bool_t AliAnalysisTaskCounter::CheckForPrimaryVertex()
 }
 
 
-
 //_____________________________________________________
 void AliAnalysisTaskCounter::FinishTaskOutput()
 {
@@ -475,3 +498,133 @@ void AliAnalysisTaskCounter::FinishTaskOutput()
     fOutputContainer->Add(histBin0); 
   
 }
+
+
+//_____________________________________
+Bool_t AliAnalysisTaskCounter::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  //
+  
+  // Fetch the aod also from the input in,
+  // have todo it in notify
+  
+  Float_t xsection = 0;
+  Float_t trials   = 1;
+  fAvgTrials = -1;
+  
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  if(!tree) return kFALSE;
+  
+  TFile *curfile = tree->GetCurrentFile();
+  
+  if(!curfile) return kFALSE;
+  
+  if(fCurrFileName == curfile->GetName()) return kFALSE;
+  
+  fCurrFileName = TString(curfile->GetName());
+  
+  if(!fh1Xsec||!fh1Trials)
+  {
+    Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+    return kFALSE;
+  }
+  
+  PythiaInfoFromFile(fCurrFileName,xsection,trials);
+  
+  fh1Xsec->Fill("<#sigma>",xsection);
+  
+  // construct a poor man average trials
+  Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+  
+  if(trials >= nEntries && nEntries > 0.) fAvgTrials = trials/nEntries;
+  
+  fh1Trials->Fill("#sum{ntrials}",trials);
+  
+  printf("AliAnalysisTaskCounter::Notify() - xs %f, trial %f, avg trials %f\n",xsection,trials, fAvgTrials);
+  
+  if(fDebug) Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
+  
+  return kTRUE;
+}
+
+//_____________________________________________________________________________________________________
+Bool_t AliAnalysisTaskCounter::PythiaInfoFromFile(TString file,Float_t & xsec,Float_t & trials)
+{
+  //
+  // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
+  // This is to called in Notify and should provide the path to the AOD/ESD file
+  
+  xsec   = 0;
+  trials = 1;
+  
+  if(file.Contains("root_archive.zip#"))
+  {
+    Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
+    Ssiz_t pos  = file.Index("#",1,pos1,TString::kExact);
+    Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
+    file.Replace(pos+1,pos2-pos1,"");
+  }
+  else
+  {
+    // not an archive take the basename....
+    file.ReplaceAll(gSystem->BaseName(file.Data()),"");
+  }
+  
+  //Printf("%s",file.Data());
+  
+  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
+  if(!fxsec)
+  {
+    // next trial fetch the histgram file
+    fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
+    if(!fxsec)
+    {
+      // not a severe condition but inciate that we have no information
+      return kFALSE;
+    }
+    else
+    {
+      // find the tlist we want to be independtent of the name so use the Tkey
+      TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
+      if(!key)
+      {
+        fxsec->Close();
+        return kFALSE;
+      }
+      
+      TList *list = dynamic_cast<TList*>(key->ReadObj());
+      if(!list)
+      {
+        fxsec->Close();
+        return kFALSE;
+      }
+      
+      xsec    = ((TProfile*)list->FindObject("h1Xsec"))  ->GetBinContent(1);
+      trials  = ((TH1F*)    list->FindObject("h1Trials"))->GetBinContent(1);
+      fxsec->Close();
+    }
+  } // no tree pyxsec.root
+  else
+  {
+    TTree *xtree = (TTree*)fxsec->Get("Xsection");
+    if(!xtree)
+    {
+      fxsec->Close();
+      return kFALSE;
+    }
+    
+    UInt_t   ntrials  = 0;
+    Double_t  xsection  = 0;
+    xtree->SetBranchAddress("xsection",&xsection);
+    xtree->SetBranchAddress("ntrials",&ntrials);
+    xtree->GetEntry(0);
+    trials = ntrials;
+    xsec = xsection;
+    fxsec->Close();
+  }
+  
+  return kTRUE;
+}
\ No newline at end of file
index 8a5fcf8d09b85df2ec550e087bbce386fc8c9109..fa360aaed384121a990a36f799de1805d6164a27 100644 (file)
@@ -23,10 +23,13 @@ class AliAnalysisTaskCounter : public AliAnalysisTaskSE {
   AliAnalysisTaskCounter(const char *name);  
   virtual ~AliAnalysisTaskCounter() ;
   
-  virtual void UserCreateOutputObjects();  
-  virtual void UserExec(Option_t *option);  
-  virtual void FinishTaskOutput();  
+  virtual void   UserCreateOutputObjects();  
+  virtual void   UserExec(Option_t *option);
+  virtual void   FinishTaskOutput();
   
+  static  Bool_t PythiaInfoFromFile(TString currFile, Float_t & xsec, Float_t & trials) ;
+  virtual Bool_t Notify();
+
   void    SetTrackMultiplicityEtaCut(Float_t eta) { fTrackMultEtaCut   = eta    ; }  
   void    SetZVertexCut(Float_t vcut)             { fZVertexCut        = vcut   ; }  
 
@@ -39,15 +42,17 @@ class AliAnalysisTaskCounter : public AliAnalysisTaskSE {
   Bool_t  IsFastClusterAccepted()       const     { return fAcceptFastCluster   ; }   
   
   Bool_t  CheckForPrimaryVertex() ;
-   
+
  private: 
   Bool_t               fAcceptFastCluster; // Accept events from fast cluster, exclude thiese events for LHC11a
   Float_t              fZVertexCut;        // Z vertex cut  
   Float_t              fTrackMultEtaCut;   // Track multiplicity eta cut  
+  Float_t              fAvgTrials;         // avg trials
   Bool_t               fCaloFilterPatch;   // CaloFilter patch  
   TList*               fOutputContainer;   //! Histogram container  
   AliESDtrackCuts    * fESDtrackCuts;      // Track cut    
   AliTriggerAnalysis * fTriggerAnalysis;   // Trigger algorithm 
+  TString              fCurrFileName;      // current file path name
   
   //Histograms
   TH1I *  fhNEvents;      //! Events that delivers the analysis frame after different assumptions  
@@ -60,10 +65,13 @@ class AliAnalysisTaskCounter : public AliAnalysisTaskSE {
   TH1F *  fhCentrality;   //! centrality
   TH1F *  fhEventPlaneAngle; //! Histogram with Event plane angle
 
+  TH1F *  fh1Xsec ;                       //! Xsec pythia
+  TH1F *  fh1Trials ;                     //! trials pythia
+  
   AliAnalysisTaskCounter(           const AliAnalysisTaskCounter&); // not implemented  
   AliAnalysisTaskCounter& operator=(const AliAnalysisTaskCounter&); // not implemented
   
-  ClassDef(AliAnalysisTaskCounter, 3);
+  ClassDef(AliAnalysisTaskCounter, 4);
 
 };