#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),
}
//________________________________________________________________________
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),
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());
}
//________________________________________________________________________
}
+
//________________________________________________________________________
void AliPWG4HighPtQAMC::CreateOutputObjects() {
//Create output objects
OpenFile(0);
fHistList = new TList();
+ fHistList->SetOwner(kTRUE);
OpenFile(1);
fHistListITS = new TList();
+ fHistListITS->SetOwner(kTRUE);
Int_t fgkNPhiBins=18;
Float_t kMinPhi = 0.;
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);
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;
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());
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));
PostData(1, fHistListITS);
}
+//________________________________________________________________________
+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 *)
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff