Added weihting with number of trials and cross section, adapted track cuts
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Oct 2010 09:10:10 +0000 (09:10 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Oct 2010 09:10:10 +0000 (09:10 +0000)
PWG4/JetTasks/AliPWG4HighPtQAMC.cxx
PWG4/JetTasks/AliPWG4HighPtQAMC.h
PWG4/JetTasks/AliPWG4HighPtQATPConly.cxx
PWG4/JetTasks/AliPWG4HighPtQATPConly.h
PWG4/JetTasks/AliPWG4HighPtSpectra.cxx
PWG4/JetTasks/AliPWG4HighPtSpectra.h
PWG4/macros/AddTaskPWG4HighPtQAMC.C
PWG4/macros/AddTaskPWG4HighPtQATPConly.C
PWG4/macros/AddTaskPWG4HighPtSpectra.C

index c5af048..19ee37b 100644 (file)
 #include "TH1.h"
 #include "TH2.h"
 #include "TH3.h"
+#include "TProfile.h"
 #include "TList.h"
+#include "TFile.h"
 #include "TChain.h"
 #include "TH3F.h"
+#include "TKey.h"
+#include "TSystem.h"
+
+#include "AliAnalysisTask.h"
 #include "AliAnalysisManager.h"
 #include "AliESDInputHandler.h"
 #include "AliMCEvent.h"
 #include "AliESDtrackCuts.h"
 #include "AliExternalTrackParam.h"
 #include "AliLog.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
+//#include "AliAnalysisHelperJetTasks.h"
 
 using namespace std; //required for resolving the 'cout' symbol
 
 ClassImp(AliPWG4HighPtQAMC)
 
-AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(): AliAnalysisTask("AliPWG4HighPtQAMC", ""), 
+AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
+: AliAnalysisTask("AliPWG4HighPtQAMC", ""), 
   fESD(0), 
   fMC(0),
   fTrackCuts(0), 
   fTrackCutsITS(0),
   fTrackType(0),
   fPtMax(100.),
+  fAvgTrials(1),
   fNEventAll(0),
   fNEventSel(0),
+  fh1Xsec(0),
+  fh1Trials(0),
+  fh1PtHard(0),
+  fh1PtHardTrials(0),
   fPtAll(0),  
   fPtSel(0),  
   fPtAllminPtMCvsPtAll(0),
@@ -86,15 +101,20 @@ AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(): AliAnalysisTask("AliPWG4HighPtQAMC", "")
 }
 //________________________________________________________________________
 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name): 
-  AliAnalysisTask(name, ""), 
+  AliAnalysisTask(name,""), 
   fESD(0),
   fMC(0),
   fTrackCuts(),
   fTrackCutsITS(),
   fTrackType(0),
   fPtMax(100.),
+  fAvgTrials(1),
   fNEventAll(0),
   fNEventSel(0),
+  fh1Xsec(0),
+  fh1Trials(0),
+  fh1PtHard(0),
+  fh1PtHardTrials(0),
   fPtAll(0),
   fPtSel(0),
   fPtAllminPtMCvsPtAll(0),
@@ -128,12 +148,9 @@ AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
   AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
   // Input slot #0 works with a TChain ESD
   DefineInput(0, TChain::Class());
-  // Output slot #0 writes into a TList
+  // Output slot #0, #1 write into a TList
   DefineOutput(0, TList::Class());
-  // Output slot #1 writes into a TList
   DefineOutput(1, TList::Class());
-  // Output slot #2 writes into a TList
-  DefineOutput(2, TList::Class());
 }
 
 //________________________________________________________________________
@@ -165,6 +182,7 @@ void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
 
 }
 
+
 //________________________________________________________________________
 void AliPWG4HighPtQAMC::CreateOutputObjects() {
   //Create output objects
@@ -175,8 +193,10 @@ void AliPWG4HighPtQAMC::CreateOutputObjects() {
 
   OpenFile(0);
   fHistList = new TList();
+  fHistList->SetOwner(kTRUE);
   OpenFile(1);
   fHistListITS = new TList();
+  fHistListITS->SetOwner(kTRUE);
 
   Int_t fgkNPhiBins=18;
   Float_t kMinPhi = 0.;
@@ -191,6 +211,20 @@ void AliPWG4HighPtQAMC::CreateOutputObjects() {
   fHistList->Add(fNEventAll);
   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
   fHistList->Add(fNEventSel);
+
+  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+  fHistList->Add(fh1Xsec);
+
+  fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+  fHistList->Add(fh1Trials);
+
+  fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
+  fHistList->Add(fh1PtHard);
+  fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
+  fHistList->Add(fh1PtHardTrials);
+
   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, fgkPtMin, fgkPtMax);
   fHistList->Add(fPtAll);
   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
@@ -325,7 +359,7 @@ void AliPWG4HighPtQAMC::Exec(Option_t *) {
   fNEventAll->Fill(0.);
 
   if (!fESD) {
-    AliDebug(2,Form("ERROR: fESD not available\n"));
+    AliDebug(2,Form("ERROR: fInputEvent not available\n"));
     PostData(0, fHistList);
     PostData(1, fHistListITS);
     return;
@@ -353,6 +387,29 @@ void AliPWG4HighPtQAMC::Exec(Option_t *) {
     return;
   }
 
+  //___ get MC information __________________________________________________________________
+
+  Double_t ptHard = 0.;
+  Double_t nTrials = 1; // trials for MC trigger weight for real data
+  
+  if(fMC){
+    AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
+     if(!pythiaGenHeader){
+       AliDebug(2,Form("ERROR: Could not retrieve AliGenPythiaEventHeader"));
+       PostData(0, fHistList);
+       PostData(1, fHistListITS);
+       return;
+     } else {
+        nTrials = pythiaGenHeader->Trials();
+        ptHard  = pythiaGenHeader->GetPtHard();
+
+        fh1PtHard->Fill(ptHard);
+        fh1PtHardTrials->Fill(ptHard,nTrials);
+
+        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
+     }
+   }
+
   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
   // Need vertex cut
   TString vtxName(vtx->GetName());
@@ -379,15 +436,15 @@ void AliPWG4HighPtQAMC::Exec(Option_t *) {
   
   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 
-  // Need to keep track of evts without vertex
-  fNEventSel->Fill(0.);
-
   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
     PostData(0, fHistList);
     PostData(1, fHistListITS);
     return;
   }
 
+  //Need to keep track of selected events
+  fNEventSel->Fill(0.);
+
   Int_t nTracks = fESD->GetNumberOfTracks();
   AliDebug(2,Form("nTracks ESD%d", nTracks));
 
@@ -492,6 +549,135 @@ void AliPWG4HighPtQAMC::Exec(Option_t *) {
 
 }
 //________________________________________________________________________
+Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
+  //
+  // 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
+  // Copied from AliAnalysisTaskJetSpectrum2
+  //
+
+  TString file(currFile);  
+  fXsec = 0;
+  fTrials = 1;
+
+  if(file.Contains("root_archive.zip#")){
+    Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
+    Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
+    file.Replace(pos+1,20,"");
+  }
+  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;
+      }
+      fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+      fTrials  = ((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);
+    fTrials = ntrials;
+    fXsec = xsection;
+    fxsec->Close();
+  }
+  return kTRUE;
+}
+//________________________________________________________________________
+Bool_t AliPWG4HighPtQAMC::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // Copied from AliAnalysisTaskJetSpectrum2
+  // 
+
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Float_t xsection = 0;
+  Float_t ftrials  = 1;
+
+  fAvgTrials = 1;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+    if(!fh1Xsec||!fh1Trials){
+      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+      return kFALSE;
+    }
+    PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
+    fh1Xsec->Fill("<#sigma>",xsection);
+    // construct a poor man average trials 
+    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+    if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
+  }  
+  return kTRUE;
+}
+
+//________________________________________________________________________
+AliGenPythiaEventHeader*  AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
+  
+  if(!mcEvent)return 0;
+  AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
+  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
+  if(!pythiaGenHeader){
+    // cocktail ??
+    AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
+    
+    if (!genCocktailHeader) {
+      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
+      //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
+      return 0;
+    }
+    TList* headerList = genCocktailHeader->GetHeaders();
+    for (Int_t i=0; i<headerList->GetEntries(); i++) {
+      pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
+      if (pythiaGenHeader)
+        break;
+    }
+    if(!pythiaGenHeader){
+      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
+      return 0;
+    }
+  }
+  return pythiaGenHeader;
+
+}
+
+//________________________________________________________________________
 void AliPWG4HighPtQAMC::Terminate(Option_t *)
 {
   // The Terminate() function is the last function to be called during
index 7defa80..e302560 100644 (file)
 class TH1F;
 class TH2F;
 class TH3F;
+class TProfile;
 class TList;
 class AliESDEvent;
 class AliESDtrackCuts;
 class AliMCEvent;
+class AliGenPythiaEventHeader;
+//class AliAnalysisHelperJetTasks;
 
 class AliPWG4HighPtQAMC: public AliAnalysisTask {
 
@@ -37,11 +40,12 @@ class AliPWG4HighPtQAMC: public AliAnalysisTask {
   AliPWG4HighPtQAMC();
   AliPWG4HighPtQAMC(const char *name);
   ~AliPWG4HighPtQAMC() {;}
-
   virtual void   ConnectInputData(Option_t *);
   virtual void   CreateOutputObjects();
   virtual void   Exec(Option_t *option);
   virtual void   Terminate(Option_t *);
+  virtual Bool_t Notify(); //Copied from AliAnalysisTaskJetSpectrum2
 
   void SetCuts(AliESDtrackCuts* trackCuts) {fTrackCuts = trackCuts;}
   void SetCutsITS(AliESDtrackCuts* trackCutsITS) {fTrackCutsITS = trackCutsITS;}
@@ -50,10 +54,12 @@ class AliPWG4HighPtQAMC: public AliAnalysisTask {
   void SetPtMax(Float_t ptmax) {fPtMax = ptmax;}
   Float_t GetPtMax()           {return fPtMax;}
 
+  static AliGenPythiaEventHeader*  GetPythiaEventHeader(AliMCEvent *mcEvent);
+  static Bool_t PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials);// get the cross section and the trails either from pyxsec.root or from pysec_hists.root
+
  protected:
 
  private:
-
   AliPWG4HighPtQAMC(const AliPWG4HighPtQAMC&);
   AliPWG4HighPtQAMC& operator=(const AliPWG4HighPtQAMC&);
 
@@ -63,13 +69,21 @@ class AliPWG4HighPtQAMC: public AliAnalysisTask {
   AliESDtrackCuts *fTrackCuts;    // TrackCuts for global reconstructed vs MC comparison
   AliESDtrackCuts *fTrackCutsITS; // TrackCuts including ITSrefit
 
-  Int_t fTrackType;               // 0: global track; 1:TPConly track
+  Int_t   fTrackType;             // 0: global track; 1:TPConly track
 
   Float_t fPtMax;                 // Maximum pT for histograms
 
+  Float_t fAvgTrials;             // Average number of trials
   
   TH1F *fNEventAll;                            //! Event counter
   TH1F *fNEventSel;                            //! Event counter
+  TProfile*     fh1Xsec;                       //! pythia cross section and trials
+  TH1F*         fh1Trials;                     //! trials which are added
+  TH1F*         fh1PtHard;                     //! pt hard of the event
+  TH1F*         fh1PtHardTrials;               //! pt hard of the event
+
+
   TH1F *fPtAll;                                //! Pt spectrum all charged particles
   TH1F *fPtSel;                                //! Pt spectrum all selected charged particles by fTrackCuts
   TH2F *fPtAllminPtMCvsPtAll;                  //! Momentum resolution (global vs MC)
index 0ab2b87..a2b7eac 100644 (file)
@@ -289,10 +289,13 @@ void AliPWG4HighPtQATPConly::CreateOutputObjects() {
 
   OpenFile(0);
   fHistList = new TList();
+  fHistList->SetOwner(kTRUE);
   OpenFile(1);
   fHistListTPC = new TList();
+  fHistListTPC->SetOwner(kTRUE);
   OpenFile(2);
   fHistListITS = new TList();
+  fHistListITS->SetOwner(kTRUE);
 
   Int_t fgkNPhiBins=18;
   Float_t kMinPhi = 0.;
index d0c91f0..5a92e2f 100644 (file)
@@ -1,4 +1,4 @@
-/**************************************************************************
+ /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
@@ -27,6 +27,7 @@ class TH1F;
 class TH2F;
 class TH3F;
 class TList;
+
 class AliESDEvent;
 class AliESDfriend;
 class AliESDfriendTrack;
@@ -34,6 +35,7 @@ class AliVEvent;
 class AliESDtrackCuts;
 class AliESDtrack;
 
+
 class AliPWG4HighPtQATPConly: public AliAnalysisTask {
 
  public:
index dab1d9f..25c7bab 100644 (file)
 #include "TH1.h"
 #include "TH2.h"
 #include "TH3.h"
+#include "TProfile.h"
 #include "TList.h"
 #include "TChain.h"
+#include "TKey.h"
+#include "TSystem.h"
+#include "TFile.h"
 
 #include "AliAnalysisManager.h"
 #include "AliESDInputHandler.h"
@@ -50,7 +54,8 @@
 #include "AliMCEvent.h"
 #include "AliMCEventHandler.h"
 #include "AliCFContainer.h"
-
+#include "AliGenPythiaEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
 
 //#include <iostream>
@@ -66,9 +71,14 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpe
   fESD(0),
   fTrackCuts(0),
   fTrackCutsTPConly(0),
+  fAvgTrials(1),
   fHistList(0),
   fNEventAll(0),
-  fNEventSel(0)
+  fNEventSel(0),
+  fh1Xsec(0),
+  fh1Trials(0),
+  fh1PtHard(0),
+  fh1PtHardTrials(0)
 {
   //
   //Default ctor
@@ -83,9 +93,14 @@ AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
   fESD(0),
   fTrackCuts(),
   fTrackCutsTPConly(0),
+  fAvgTrials(1),
   fHistList(0),
   fNEventAll(0),
-  fNEventSel(0)
+  fNEventSel(0),
+  fh1Xsec(0),
+  fh1Trials(0),
+  fh1PtHard(0),
+  fh1PtHardTrials(0)
 {
   //
   // Constructor. Initialization of Inputs and Outputs
@@ -189,6 +204,28 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
     
     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
   }
+
+  //___ get MC information __________________________________________________________________
+
+  Double_t ptHard = 0.;
+  Double_t nTrials = 1; // trials for MC trigger weight for real data
+  
+  if(mcEvent){
+    AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(mcEvent);
+     if(!pythiaGenHeader){
+       AliDebug(2,Form("ERROR: Could not retrieve AliGenPythiaEventHeader"));
+       PostData(0, fHistList);
+       return;
+     } else {
+        nTrials = pythiaGenHeader->Trials();
+        ptHard  = pythiaGenHeader->GetPtHard();
+
+        fh1PtHard->Fill(ptHard);
+        fh1PtHardTrials->Fill(ptHard,nTrials);
+
+        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
+     }
+  }
   
   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
@@ -358,6 +395,133 @@ void AliPWG4HighPtSpectra::Exec(Option_t *)
   PostData(2,fCFManagerNeg->GetParticleContainer());
   
 }
+//________________________________________________________________________
+Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
+  //
+  // 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
+  // Copied from AliAnalysisTaskJetSpectrum2
+  //
+
+  TString file(currFile);  
+  fXsec = 0;
+  fTrials = 1;
+
+  if(file.Contains("root_archive.zip#")){
+    Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
+    Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
+    file.Replace(pos+1,20,"");
+  }
+  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;
+      }
+      fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+      fTrials  = ((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);
+    fTrials = ntrials;
+    fXsec = xsection;
+    fxsec->Close();
+  }
+  return kTRUE;
+}
+//________________________________________________________________________
+Bool_t AliPWG4HighPtSpectra::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // Copied from AliAnalysisTaskJetSpectrum2
+  // 
+
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Float_t xsection = 0;
+  Float_t ftrials  = 1;
+
+  fAvgTrials = 1;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+    if(!fh1Xsec||!fh1Trials){
+      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+      return kFALSE;
+    }
+     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
+    fh1Xsec->Fill("<#sigma>",xsection);
+    // construct a poor man average trials 
+    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+    if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
+  }  
+  return kTRUE;
+}
+
+//________________________________________________________________________
+AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
+  
+  if(!mcEvent)return 0;
+  AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
+  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
+  if(!pythiaGenHeader){
+    // cocktail ??
+    AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
+    
+    if (!genCocktailHeader) {
+      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
+      //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
+      return 0;
+    }
+    TList* headerList = genCocktailHeader->GetHeaders();
+    for (Int_t i=0; i<headerList->GetEntries(); i++) {
+      pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
+      if (pythiaGenHeader)
+        break;
+    }
+    if(!pythiaGenHeader){
+      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
+      return 0;
+    }
+  }
+  return pythiaGenHeader;
+
+}
 
 
 //___________________________________________________________________________
@@ -382,11 +546,25 @@ void AliPWG4HighPtSpectra::CreateOutputObjects() {
   //slot #1
   OpenFile(0);
   fHistList = new TList();
+  fHistList->SetOwner(kTRUE);
   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
   fHistList->Add(fNEventAll);
   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
   fHistList->Add(fNEventSel);
 
+  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+  fHistList->Add(fh1Xsec);
+
+  fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+  fHistList->Add(fh1Trials);
+
+  fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
+  fHistList->Add(fh1PtHard);
+  fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
+  fHistList->Add(fh1PtHardTrials);
+
   TH1::AddDirectory(oldStatus);   
 
 }
index 4c64acc..71a0754 100644 (file)
 class TH1I;
 class TH1F;
 class TH1D;
-class TFile ;
+class TProfile;
+class TFile;
+class TList;
+
 //class AliCFManager;
 class AliESDtrackCuts;
 class AliESDEvent;
+class AliMCEvent;
+class AliGenPythiaEventHeader;
 
 class AliPWG4HighPtSpectra : public AliAnalysisTask {
  public:
@@ -55,6 +60,7 @@ class AliPWG4HighPtSpectra : public AliAnalysisTask {
   virtual void   CreateOutputObjects();
   virtual void   Exec(Option_t *option);
   virtual void   Terminate(Option_t *);
+  virtual Bool_t Notify(); //Copied from AliAnalysisTaskJetSpectrum2
 
   // CORRECTION FRAMEWORK RELATED FUNCTIONS
   void     SetCFManagerPos(const AliCFManager* io1) {fCFManagerPos = io1;}   // global correction manager 
@@ -69,6 +75,9 @@ class AliPWG4HighPtSpectra : public AliAnalysisTask {
   // Data types
   Bool_t IsReadAODData()   const {return fReadAODData;}
   void   SetReadAODData(Bool_t flag=kTRUE) {fReadAODData=flag;}
+
+  static AliGenPythiaEventHeader*  GetPythiaEventHeader(AliMCEvent *mcEvent);
+  static Bool_t PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials);// get the cross section and the trails either from pyxsec.root or from pysec_hists.root
   
  protected:
   Bool_t              fReadAODData ;       // flag for AOD/ESD input files
@@ -84,12 +93,19 @@ class AliPWG4HighPtSpectra : public AliAnalysisTask {
   AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra&);
   AliPWG4HighPtSpectra& operator=(const AliPWG4HighPtSpectra&);
 
+  Float_t fAvgTrials;             // Average number of trials
+
   // Histograms
   //Number of events
-  TList *fHistList;            //! List of output histograms
+  TList *fHistList;             //! List of output histograms
   TH1F  *fNEventAll;            //! Event counter
   TH1F  *fNEventSel;            //! Event counter: Selected events for analysis
 
+  TProfile*     fh1Xsec;                       //! pythia cross section and trials
+  TH1F*         fh1Trials;                     //! trials which are added
+  TH1F*         fh1PtHard;                     //! pt hard of the event
+  TH1F*         fh1PtHardTrials;               //! pt hard of the event
+
   ClassDef(AliPWG4HighPtSpectra,2);
 };
 
index a7346be..170dafe 100644 (file)
@@ -31,16 +31,18 @@ AliPWG4HighPtQAMC* AddTaskPWG4HighPtQAMC(char *prodType = "LHC10e14", int trackT
   AliESDtrackCuts *trackCuts = new AliESDtrackCuts("AliESDtrackCuts","Standard Cuts");
   //Standard Cuts
   //Set track cuts for global tracks
-  if(trackType==0) trackCuts = trackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
+  if(trackType==0) trackCuts = trackCuts->GetStandardITSTPCTrackCuts2010(kTRUE);//Primary Track Selection
   //Set track cuts for TPConly tracks
-  if(trackType==1) trackCuts = trackCuts->GetStandardTPCOnlyTrackCuts(); 
+  if(trackType==1) { 
+    trackCuts = trackCuts->GetStandardTPCOnlyTrackCuts(); 
+    trackCuts->SetMinNClustersTPC(70);
+  }
   trackCuts->SetEtaRange(-0.9,0.9);
   trackCuts->SetPtRange(0.15, 1e10);
-  trackCuts->SetRequireITSRefit(kFALSE);
   
   AliESDtrackCuts *trackCutsITS = new AliESDtrackCuts("AliESDtrackCuts","Standard Cuts with ITSrefit");
   //Standard Cuts
-  trackCutsITS=trackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
+  trackCutsITS=trackCuts->GetStandardITSTPCTrackCuts2010(kTRUE);//Primary Track Selection
   trackCutsITS->SetEtaRange(-0.9,0.9);
   trackCutsITS->SetPtRange(0.15, 1e10);
   trackCutsITS->SetRequireITSRefit(kTRUE);
@@ -50,6 +52,7 @@ AliPWG4HighPtQAMC* AddTaskPWG4HighPtQAMC(char *prodType = "LHC10e14", int trackT
   taskPWG4QAMC->SetCuts(trackCuts);
   taskPWG4QAMC->SetCutsITS(trackCutsITS);
   taskPWG4QAMC->SetTrackType(trackType);
+  
   if(!strcmp(prodType, "LHC10e14")  || !strcmp(prodType, "PbPb")) taskPWG4QAMC->SetPtMax(500.);
   else taskPWG4QAMC->SetPtMax(100.);
  
@@ -63,14 +66,14 @@ AliPWG4HighPtQAMC* AddTaskPWG4HighPtQAMC(char *prodType = "LHC10e14", int trackT
   TString outputfile = AliAnalysisManager::GetCommonFileName();
   outputfile += Form(":PWG4_HighPtQAMC%d",trackType);
   
-  AliAnalysisDataContainer *cout_hist0 = mgr->CreateContainer(Form("qa_histsMC%d",trackType), TList::Class(), AliAnalysisManager::kOutputContainer,outputfile);
+  AliAnalysisDataContainer *cout_hist1 = mgr->CreateContainer(Form("qa_histsMC%d",trackType), TList::Class(), AliAnalysisManager::kOutputContainer,outputfile);
   AliAnalysisDataContainer *cout_hist2 = mgr->CreateContainer(Form("qa_histsMCITS%d",trackType), TList::Class(), AliAnalysisManager::kOutputContainer,outputfile);  
 
   mgr->AddTask(taskPWG4QAMC);
   mgr->ConnectInput(taskPWG4QAMC,0,mgr->GetCommonInputContainer());
-  mgr->ConnectOutput(taskPWG4QAMC,0,cout_hist0);
+  mgr->ConnectOutput(taskPWG4QAMC,0,cout_hist1);
   mgr->ConnectOutput(taskPWG4QAMC,1,cout_hist2);
 
-  // Return task pointer at the end
+    // Return task pointer at the end
   return taskPWG4QAMC;
 }
index 0dd8f16..5805522 100644 (file)
@@ -29,21 +29,13 @@ AliPWG4HighPtQATPConly* AddTaskPWG4HighPtQATPConly(char *prodType = "LHC10e14",i
   //Use AliESDtrackCuts
   AliESDtrackCuts *trackCuts = new AliESDtrackCuts("AliESDtrackCuts","Standard Cuts");
   if(cuts==1) {
-    //Standard Cuts
-    trackCuts->SetAcceptKinkDaughters(kFALSE);
-    trackCuts->SetRequireTPCStandAlone(kTRUE); 
-    trackCuts->SetRequireTPCRefit(kTRUE);
-    trackCuts->SetMinNClustersTPC(70);
+    trackCuts=trackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
     trackCuts->SetEtaRange(-0.9,0.9);
-    trackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
     trackCuts->SetPtRange(0.15, 1e10);
-    trackCuts->SetMaxChi2PerClusterTPC(3.5);
-    trackCuts->SetMaxDCAToVertexXY(2.4);
-    trackCuts->SetMaxDCAToVertexZ(3.2);
-    trackCuts->SetDCAToVertex2D(kTRUE);
+    trackCuts->SetRequireITSRefit(kFALSE);
   }
   else if(cuts==2) {
-    trackCuts=trackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
+    trackCuts=trackCuts->GetStandardITSTPCTrackCuts2010(kTRUE);//Primary Track Selection
     trackCuts->SetEtaRange(-0.9,0.9);
     trackCuts->SetPtRange(0.15, 1e10);
     trackCuts->SetRequireITSRefit(kFALSE);
@@ -51,20 +43,12 @@ AliPWG4HighPtQATPConly* AddTaskPWG4HighPtQATPConly(char *prodType = "LHC10e14",i
 
   AliESDtrackCuts *trackCutsITS = new AliESDtrackCuts("AliESDtrackCuts","Standard Cuts with ITSrefit");
   if(cuts==1) {
-    trackCutsITS->SetAcceptKinkDaughters(kFALSE);
-    trackCutsITS->SetRequireTPCRefit(kTRUE);
-    trackCutsITS->SetEtaRange(-0.9,0.9);
-    trackCutsITS->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
-    trackCutsITS->SetPtRange(0.15, 1e10);
-    trackCutsITS->SetMinNClustersTPC(70);
-    trackCutsITS->SetMaxChi2PerClusterTPC(3.5);
-    trackCutsITS->SetRequireITSRefit(kTRUE);
-    trackCutsITS->SetMaxDCAToVertexXY(2.4);
-    trackCutsITS->SetMaxDCAToVertexZ(3.2);
-    trackCutsITS->SetDCAToVertex2D(kTRUE); 
+   trackCutsITS=trackCutsITS->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
+   trackCutsITS->SetEtaRange(-0.9,0.9);
+   trackCutsITS->SetPtRange(0.15, 1e10); 
   }
  else if(cuts==2) {
-   trackCutsITS=trackCutsITS->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Selection
+   trackCutsITS=trackCutsITS->GetStandardITSTPCTrackCuts2010(kTRUE);//Primary Track Selection
    trackCutsITS->SetEtaRange(-0.9,0.9);
    trackCutsITS->SetPtRange(0.15, 1e10);
  }
index 27b4493..32273dd 100644 (file)
@@ -117,7 +117,7 @@ AliPWG4HighPtSpectra* AddTaskPWG4HighPtSpectra(char *prodType = "LHC10e14")
   //CREATE THE  CUTS -----------------------------------------------
   //Use AliESDtrackCuts
   AliESDtrackCuts *trackCuts = new AliESDtrackCuts("AliESDtrackCuts","Standard Cuts");
-  trackCuts = trackCuts->GetStandardITSTPCTrackCuts2009(kTRUE);//Primary Track Quality Selection for Global tracks
+  trackCuts = trackCuts->GetStandardITSTPCTrackCuts2010(kTRUE);//Primary Track Quality Selection for Global tracks
   trackCuts->SetEtaRange(-0.9,0.9);
   trackCuts->SetPtRange(0.15, 1e10);