#include "TH2F.h"
#include "TString.h"
#include "THnSparse.h"
+#include "TProfile.h"
+#include "TFile.h"
+#include "TKey.h"
#include "AliAODInputHandler.h"
#include "AliAODHandler.h"
,fBranchGenJets("")
,fTrackTypeGen(0)
,fJetTypeGen(0)
+ ,fJetTypeRecEff(0)
,fFilterMask(0)
+ ,fUsePhysicsSelection(kTRUE)
,fTrackPtCut(0)
,fTrackEtaMin(0)
,fTrackEtaMax(0)
,fTracksRec(0)
,fTracksRecCuts(0)
,fTracksGen(0)
+ ,fTracksAODMCCharged(0)
+ ,fTracksRecQualityCuts(0)
,fJetsRec(0)
,fJetsRecCuts(0)
,fJetsGen(0)
+ ,fJetsRecEff(0)
,fQATrackHistosRec(0)
,fQATrackHistosRecCuts(0)
,fQATrackHistosGen(0)
,fQAJetHistosRecCutsLeading(0)
,fQAJetHistosGen(0)
,fQAJetHistosGenLeading(0)
+ ,fQAJetHistosRecEffLeading(0)
,fFFHistosRecCuts(0)
,fFFHistosRecLeading(0)
,fFFHistosRecLeadingTrack(0)
,fh1VertexNContributors(0)
,fh1VertexZ(0)
,fh1EvtMult(0)
+ ,fh1Xsec(0)
+ ,fh1Trials(0)
+ ,fh1PtHard(0)
+ ,fh1PtHardTrials(0)
,fh1nRecJetsCuts(0)
,fh1nGenJets(0)
+ ,fh1nRecEffJets(0)
+ ,fhnSingleTrackRecEffHisto(0)
+ ,fhnJetTrackRecEffHisto(0)
{
// default constructor
}
,fBranchGenJets("")
,fTrackTypeGen(0)
,fJetTypeGen(0)
+ ,fJetTypeRecEff(0)
,fFilterMask(0)
+ ,fUsePhysicsSelection(kTRUE)
,fTrackPtCut(0)
,fTrackEtaMin(0)
,fTrackEtaMax(0)
,fTracksRec(0)
,fTracksRecCuts(0)
,fTracksGen(0)
+ ,fTracksAODMCCharged(0)
+ ,fTracksRecQualityCuts(0)
,fJetsRec(0)
,fJetsRecCuts(0)
,fJetsGen(0)
+ ,fJetsRecEff(0)
,fQATrackHistosRec(0)
,fQATrackHistosRecCuts(0)
,fQATrackHistosGen(0)
,fQAJetHistosRecCutsLeading(0)
,fQAJetHistosGen(0)
,fQAJetHistosGenLeading(0)
+ ,fQAJetHistosRecEffLeading(0)
,fFFHistosRecCuts(0)
,fFFHistosRecLeading(0)
,fFFHistosRecLeadingTrack(0)
,fh1VertexNContributors(0)
,fh1VertexZ(0)
,fh1EvtMult(0)
+ ,fh1Xsec(0)
+ ,fh1Trials(0)
+ ,fh1PtHard(0)
+ ,fh1PtHardTrials(0)
,fh1nRecJetsCuts(0)
,fh1nGenJets(0)
+ ,fh1nRecEffJets(0)
+ ,fhnSingleTrackRecEffHisto(0)
+ ,fhnJetTrackRecEffHisto(0)
{
// constructor
,fBranchGenJets(copy.fBranchGenJets)
,fTrackTypeGen(copy.fTrackTypeGen)
,fJetTypeGen(copy.fJetTypeGen)
+ ,fJetTypeRecEff(copy.fJetTypeRecEff)
,fFilterMask(copy.fFilterMask)
+ ,fUsePhysicsSelection(copy.fUsePhysicsSelection)
,fTrackPtCut(copy.fTrackPtCut)
,fTrackEtaMin(copy.fTrackEtaMin)
,fTrackEtaMax(copy.fTrackEtaMax)
,fJetEtaMin(copy.fJetEtaMin)
,fJetEtaMax(copy.fJetEtaMax)
,fJetPhiMin(copy.fJetPhiMin)
- ,fJetPhiMax(copy.fJetPhiMin)
+ ,fJetPhiMax(copy.fJetPhiMax)
,fDiJetCut(copy.fDiJetCut)
,fDiJetDeltaPhiCut(copy.fDiJetDeltaPhiCut)
,fDiJetPtFractionCut(copy.fDiJetPtFractionCut)
,fTracksRec(copy.fTracksRec)
,fTracksRecCuts(copy.fTracksRecCuts)
,fTracksGen(copy.fTracksGen)
+ ,fTracksAODMCCharged(copy.fTracksAODMCCharged)
+ ,fTracksRecQualityCuts(copy.fTracksRecQualityCuts)
,fJetsRec(copy.fJetsRec)
,fJetsRecCuts(copy.fJetsRecCuts)
,fJetsGen(copy.fJetsGen)
+ ,fJetsRecEff(copy.fJetsRecEff)
,fQATrackHistosRec(copy.fQATrackHistosRec)
,fQATrackHistosRecCuts(copy.fQATrackHistosRecCuts)
,fQATrackHistosGen(copy.fQATrackHistosGen)
,fQAJetHistosRecCutsLeading(copy.fQAJetHistosRecCutsLeading)
,fQAJetHistosGen(copy.fQAJetHistosGen)
,fQAJetHistosGenLeading(copy.fQAJetHistosGenLeading)
+ ,fQAJetHistosRecEffLeading(copy.fQAJetHistosRecEffLeading)
,fFFHistosRecCuts(copy.fFFHistosRecCuts)
,fFFHistosRecLeading(copy.fFFHistosRecLeading)
,fFFHistosRecLeadingTrack(copy.fFFHistosRecLeadingTrack)
,fh1VertexNContributors(copy.fh1VertexNContributors)
,fh1VertexZ(copy.fh1VertexZ)
,fh1EvtMult(copy.fh1EvtMult)
+ ,fh1Xsec(copy.fh1Xsec)
+ ,fh1Trials(copy.fh1Trials)
+ ,fh1PtHard(copy.fh1PtHard)
+ ,fh1PtHardTrials(copy.fh1PtHardTrials)
,fh1nRecJetsCuts(copy.fh1nRecJetsCuts)
,fh1nGenJets(copy.fh1nGenJets)
+ ,fh1nRecEffJets(copy.fh1nRecEffJets)
+ ,fhnSingleTrackRecEffHisto(copy.fhnSingleTrackRecEffHisto)
+ ,fhnJetTrackRecEffHisto(copy.fhnJetTrackRecEffHisto)
{
// copy constructor
fBranchGenJets = o.fBranchGenJets;
fTrackTypeGen = o.fTrackTypeGen;
fJetTypeGen = o.fJetTypeGen;
+ fJetTypeRecEff = o.fJetTypeRecEff;
fFilterMask = o.fFilterMask;
+ fUsePhysicsSelection = o.fUsePhysicsSelection;
fTrackPtCut = o.fTrackPtCut;
fTrackEtaMin = o.fTrackEtaMin;
fTrackEtaMax = o.fTrackEtaMax;
fTracksRec = o.fTracksRec;
fTracksRecCuts = o.fTracksRecCuts;
fTracksGen = o.fTracksGen;
+ fTracksAODMCCharged = o.fTracksAODMCCharged;
+ fTracksRecQualityCuts = o.fTracksRecQualityCuts;
fJetsRec = o.fJetsRec;
fJetsRecCuts = o.fJetsRecCuts;
fJetsGen = o.fJetsGen;
+ fJetsRecEff = o.fJetsRecEff;
fQATrackHistosRec = o.fQATrackHistosRec;
fQATrackHistosRecCuts = o.fQATrackHistosRecCuts;
fQATrackHistosGen = o.fQATrackHistosGen;
fQAJetHistosRecCutsLeading = o.fQAJetHistosRecCutsLeading;
fQAJetHistosGen = o.fQAJetHistosGen;
fQAJetHistosGenLeading = o.fQAJetHistosGenLeading;
+ fQAJetHistosRecEffLeading = o.fQAJetHistosRecEffLeading;
fFFHistosRecCuts = o.fFFHistosRecCuts;
fFFHistosRecLeading = o.fFFHistosRecLeading;
fFFHistosRecLeadingTrack = o.fFFHistosRecLeadingTrack;
fh1VertexNContributors = o.fh1VertexNContributors;
fh1VertexZ = o.fh1VertexZ;
fh1EvtMult = o.fh1EvtMult;
+ fh1Xsec = o.fh1Xsec;
+ fh1Trials = o.fh1Trials;
+ fh1PtHard = o.fh1PtHard;
+ fh1PtHardTrials = o.fh1PtHardTrials;
fh1nRecJetsCuts = o.fh1nRecJetsCuts;
fh1nGenJets = o.fh1nGenJets;
+ fh1nRecEffJets = o.fh1nRecEffJets;
+ fhnSingleTrackRecEffHisto = o.fhnSingleTrackRecEffHisto;
+ fhnJetTrackRecEffHisto = o.fhnJetTrackRecEffHisto;
}
return *this;
{
// destructor
- if(fTracksRec) delete fTracksRec;
- if(fTracksRecCuts) delete fTracksRecCuts;
- if(fTracksGen) delete fTracksGen;
- if(fJetsRec) delete fJetsRec;
- if(fJetsRecCuts) delete fJetsRecCuts;
- if(fJetsGen) delete fJetsGen;
+ if(fTracksRec) delete fTracksRec;
+ if(fTracksRecCuts) delete fTracksRecCuts;
+ if(fTracksGen) delete fTracksGen;
+ if(fTracksAODMCCharged) delete fTracksAODMCCharged;
+ if(fTracksRecQualityCuts) delete fTracksRecQualityCuts;
+ if(fJetsRec) delete fJetsRec;
+ if(fJetsRecCuts) delete fJetsRecCuts;
+ if(fJetsGen) delete fJetsGen;
+ if(fJetsRecEff) delete fJetsRecEff;
// if(fDiJetBins) delete fDiJetBins;
if(incrementJetPt) fh1JetPt->Fill(jetPt);
fh2TrackPt->Fill(jetPt,trackPt);
- // Added by me
Double_t z = 0.;
if(jetPt>0) z = trackPt / jetPt;
Double_t xi = 0;
}
+
+Bool_t AliAnalysisTaskFragmentationFunction::Notify()
+{
+ //
+ // Implemented Notify() to read the cross sections
+ // and number of trials from pyxsec.root
+ // (taken from AliAnalysisTaskJetSpectrum)
+ //
+ TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+ Double_t xsection = 0;
+ UInt_t ntrials = 0;
+ Float_t ftrials = 0;
+ 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;
+ }
+
+ TString fileName(curfile->GetName());
+ if(fileName.Contains("AliESDs.root")){
+ fileName.ReplaceAll("AliESDs.root", "");
+ }
+ else if(fileName.Contains("AliAOD.root")){
+ fileName.ReplaceAll("AliAOD.root", "");
+ }
+ else if(fileName.Contains("AliAODs.root")){
+ fileName.ReplaceAll("AliAODs.root", "");
+ }
+ else if(fileName.Contains("galice.root")){
+ // for running with galice and kinematics alone...
+ fileName.ReplaceAll("galice.root", "");
+ }
+ TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
+ if(!fxsec){
+ if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
+ // next trial fetch the histgram file
+ fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+ if(!fxsec){
+ // not a severe condition
+ if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+ return kTRUE;
+ }
+ 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){
+ if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);
+ return kTRUE;
+ }
+ TList *list = dynamic_cast<TList*>(key->ReadObj());
+ if(!list){
+ if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);
+ return kTRUE;
+ }
+ xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+ ftrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+ }
+ }
+ else{
+ TTree *xtree = (TTree*)fxsec->Get("Xsection");
+ if(!xtree){
+ Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+ }
+ xtree->SetBranchAddress("xsection",&xsection);
+ xtree->SetBranchAddress("ntrials",&ntrials);
+ ftrials = ntrials;
+ xtree->GetEntry(0);
+ }
+ fh1Xsec->Fill("<#sigma>",xsection);
+ fh1Trials->Fill("#sum{ntrials}",ftrials);
+ }
+
+ return kTRUE;
+}
+
+
+
//__________________________________________________________________
void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
{
fTracksGen = new TList();
fTracksGen->SetOwner(kFALSE);
+ fTracksAODMCCharged = new TList();
+ fTracksAODMCCharged->SetOwner(kFALSE);
+
+ fTracksRecQualityCuts = new TList();
+ fTracksRecQualityCuts->SetOwner(kFALSE);
+
fJetsRec = new TList();
fJetsRec->SetOwner(kFALSE);
fJetsGen = new TList();
fJetsGen->SetOwner(kFALSE);
+ fJetsRecEff = new TList();
+ fJetsRecEff->SetOwner(kFALSE);
+
// fJetsKine = new TList();
// fJetsKine->SetOwner(kTRUE); // delete AOD jets using mom from Kine Tree via TList::Clear()
fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
fh1VertexNContributors = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
fh1VertexZ = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
- fh1EvtMult = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",100,0.,100.);
+ fh1EvtMult = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",120,0.,120.);
+
+ fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+ fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+ fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
+ fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+ fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
+ fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
+
fh1nRecJetsCuts = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
fh1nGenJets = new TH1F("fh1nGenJets","generated jets per event",10,-0.5,9.5);
+ fh1nRecEffJets = new TH1F("fh1nRecEffJets","reconstruction effiency: jets per event",10,-0.5,9.5);
+
+ // 5D single track eff histo: phi:eta:gen pt:rec pt:isReconstructed - use binning as for track QA
+ Int_t nBinsSingleTrackEffHisto[5] = { fQATrackNBinsPhi, fQATrackNBinsEta, fQATrackNBinsPt, fQATrackNBinsPt, 2 };
+ Double_t binMinSingleTrackEffHisto[5] = { fQATrackPhiMin, fQATrackEtaMin, fQATrackPtMin, fQATrackPtMin, 0 };
+ Double_t binMaxSingleTrackEffHisto[5] = { fQATrackPhiMax, fQATrackEtaMax, fQATrackPtMax, fQATrackPtMax, 2 };
+ const char* labelsSingleTrackEffHisto[5] = {"#phi","#eta","gen p_{T} [GeV/c]", "rec p_{T} [GeV/c]", "isRec"};
+
+ fhnSingleTrackRecEffHisto = new THnSparseF("fhnSingleTrackRecEffHisto","generated tracks phi:eta:pt:isReconstructed",5,
+ nBinsSingleTrackEffHisto,binMinSingleTrackEffHisto,binMaxSingleTrackEffHisto);
+
+ AliAnalysisTaskFragmentationFunction::SetProperties(fhnSingleTrackRecEffHisto,5,labelsSingleTrackEffHisto);
+
+
+ // 7D jets track eff histo: jet phi:eta:pt:track pt:z:xi:isReconstructed - use binning as for track/jet QA
+ Int_t nBinsJetTrackEffHisto[7] = { fQAJetNBinsPhi, fQAJetNBinsEta, fFFNBinsJetPt, fFFNBinsPt, fFFNBinsZ, fFFNBinsXi, 2};
+ Double_t binMinJetTrackEffHisto[7] = { fQAJetPhiMin, fQAJetEtaMin, fFFJetPtMin , fFFPtMin, fFFZMin , fFFXiMin, 0 };
+ Double_t binMaxJetTrackEffHisto[7] = { fQAJetPhiMax, fQAJetEtaMax, fFFJetPtMax , fFFPtMax, fFFZMax , fFFXiMax, 2 };
+ const char* labelsJetTrackEffHisto[7] = {"jet #phi","jet #eta","jet p_{T} [GeV/c]","track p_{T} [GeV/c]","z","#xi","isRec"};
+
+ fhnJetTrackRecEffHisto = new THnSparseF("fhnJetTrackRecEffHisto","generated tracks - jet phi:jet eta:jet pt:track pt:z:xi:isReconstructed",7,
+ nBinsJetTrackEffHisto,binMinJetTrackEffHisto,binMaxJetTrackEffHisto);
+
+ AliAnalysisTaskFragmentationFunction::SetProperties(fhnJetTrackRecEffHisto,7,labelsJetTrackEffHisto);
fQATrackHistosRec = new AliFragFuncQATrackHistos("Rec", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
fQAJetHistosGenLeading = new AliFragFuncQAJetHistos("GenLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
- fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
-
+ fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+ fQAJetHistosRecEffLeading = new AliFragFuncQAJetHistos("RecEffLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax,
+ fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
+
fFFHistosRecCuts = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
fFFNBinsPt, fFFPtMin, fFFPtMax,
fQAJetHistosRecCutsLeading->DefineHistos();
fQAJetHistosGen->DefineHistos();
fQAJetHistosGenLeading->DefineHistos();
+ fQAJetHistosRecEffLeading->DefineHistos();
fFFHistosRecCuts->DefineHistos();
fFFHistosRecLeading->DefineHistos();
fQADiJetHistosRecCuts->DefineQADiJetHistos();
fQADiJetHistosGen->DefineQADiJetHistos();
+ Bool_t genJets = (fJetTypeGen != kJetsUndef) ? kTRUE : kFALSE;
+ Bool_t genTracks = (fTrackTypeGen != kTrackUndef) ? kTRUE : kFALSE;
+ Bool_t recJetsEff = (fJetTypeRecEff != kJetsUndef) ? kTRUE : kFALSE;
+
+
const Int_t saveLevel = 5;
if(saveLevel>0){
fCommonHistList->Add(fh1EvtSelection);
fFFHistosRecCuts->AddToOutput(fCommonHistList);
fFFHistosRecLeading->AddToOutput(fCommonHistList);
fFFHistosRecLeadingTrack->AddToOutput(fCommonHistList);
- fFFHistosGen->AddToOutput(fCommonHistList);
- fFFHistosGenLeading->AddToOutput(fCommonHistList);
- fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+ if(genJets && genTracks){
+ fFFHistosGen->AddToOutput(fCommonHistList);
+ fFFHistosGenLeading->AddToOutput(fCommonHistList);
+ fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+ fCommonHistList->Add(fh1Xsec);
+ fCommonHistList->Add(fh1Trials);
+ fCommonHistList->Add(fh1PtHard);
+ fCommonHistList->Add(fh1PtHardTrials);
+ }
}
if(saveLevel>1){
fQATrackHistosRec->AddToOutput(fCommonHistList);
fQATrackHistosRecCuts->AddToOutput(fCommonHistList);
- fQATrackHistosGen->AddToOutput(fCommonHistList);
+ if(genTracks) fQATrackHistosGen->AddToOutput(fCommonHistList);
fQAJetHistosRec->AddToOutput(fCommonHistList);
fQAJetHistosRecCuts->AddToOutput(fCommonHistList);
fQAJetHistosRecCutsLeading->AddToOutput(fCommonHistList);
- fQAJetHistosGen->AddToOutput(fCommonHistList);
- fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
-
+ if(recJetsEff) fQAJetHistosRecEffLeading->AddToOutput(fCommonHistList);
+
+ if(genJets){
+ fQAJetHistosGen->AddToOutput(fCommonHistList);
+ fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
+ }
+
fCommonHistList->Add(fh1EvtMult);
fCommonHistList->Add(fh1nRecJetsCuts);
- fCommonHistList->Add(fh1nGenJets);
+ if(genJets) fCommonHistList->Add(fh1nGenJets);
}
if(saveLevel>2){
fCommonHistList->Add(fh1VertexNContributors);
fIJHistosRecCuts->AddToOutput(fCommonHistList);
fIJHistosRecLeading->AddToOutput(fCommonHistList);
fIJHistosRecLeadingTrack->AddToOutput(fCommonHistList);
- fIJHistosGen->AddToOutput(fCommonHistList);
- fIJHistosGenLeading->AddToOutput(fCommonHistList);
- fIJHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+ if(genJets && genTracks){
+ fIJHistosGen->AddToOutput(fCommonHistList);
+ fIJHistosGenLeading->AddToOutput(fCommonHistList);
+ fIJHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+ }
}
if(saveLevel>4){
fFFDiJetHistosRecCuts->AddToOutput(fCommonHistList);
fFFDiJetHistosRecLeading->AddToOutput(fCommonHistList);
fFFDiJetHistosRecLeadingTrack->AddToOutput(fCommonHistList);
- fFFDiJetHistosGen->AddToOutput(fCommonHistList);
- fFFDiJetHistosGenLeading->AddToOutput(fCommonHistList);
- fFFDiJetHistosGenLeadingTrack->AddToOutput(fCommonHistList);
fQADiJetHistosRecCuts->AddToOutput(fCommonHistList);
- fQADiJetHistosGen->AddToOutput(fCommonHistList);
+
+ if(genJets && genTracks){
+ fFFDiJetHistosGen->AddToOutput(fCommonHistList);
+ fFFDiJetHistosGenLeading->AddToOutput(fCommonHistList);
+ fFFDiJetHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+ fQADiJetHistosGen->AddToOutput(fCommonHistList);
+ }
+
+ if(recJetsEff && genTracks){
+ fCommonHistList->Add(fhnSingleTrackRecEffHisto);
+ fCommonHistList->Add(fhnJetTrackRecEffHisto);
+ fCommonHistList->Add(fh1nRecEffJets);
+ }
}
// =========== Switch on Sumw2 for all histos ===========
AliInputEventHandler* inputHandler = (AliInputEventHandler*)
((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
- if(inputHandler->IsEventSelected()&AliVEvent::kMB){
+ if(inputHandler->IsEventSelected() & AliVEvent::kMB){
if(fDebug > 1) Printf(" Trigger Selection: event ACCEPTED ... ");
fh1EvtSelection->Fill(1.);
} else {
fh1EvtSelection->Fill(0.);
- if(inputHandler->InheritsFrom("AliESDInputHandler")){ // PhysicsSelection only with ESD input
+ if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
PostData(1, fCommonHistList);
return;
}
if (fDebug > 1) Printf("%s:%d primary vertex selection: event ACCEPTED ...",(char*)__FILE__,__LINE__);
fh1EvtSelection->Fill(5.);
+
+
+ //___ get MC information __________________________________________________________________
+
+ Double_t ptHard = 0.;
+ Double_t nTrials = 1; // trials for MC trigger weight for real data
+
+ AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMCEvent);
+ if(!pythiaGenHeader){
+ if(fJetTypeGen != kJetsUndef && fTrackTypeGen != kTrackUndef){
+ Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
+ return;
+ }
+ } else {
+ nTrials = pythiaGenHeader->Trials();
+ ptHard = pythiaGenHeader->GetPtHard();
+
+ fh1PtHard->Fill(ptHard);
+ fh1PtHardTrials->Fill(ptHard,nTrials);
+ }
//___ fetch jets __________________________________________________________________________
fh1nGenJets->Fill(nGenJets);
+ if(fJetTypeRecEff==kJetsKine || fJetTypeRecEff == kJetsKineAcceptance) fJetsRecEff->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear()
+ Int_t nJRecEff = GetListOfJets(fJetsRecEff, fJetTypeRecEff);
+ Int_t nRecEffJets = 0;
+ if(nJRecEff>=0) nRecEffJets = fJetsRecEff->GetEntries();
+ if(fDebug>2)Printf("%s:%d Selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
+ if(nJRecEff != nRecEffJets) Printf("%s:%d Mismatch selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
+ fh1nRecEffJets->Fill(nRecEffJets);
+
+
//____ fetch particles __________________________________________________________
Int_t nT = GetListOfTracks(fTracksRec, kTrackAOD);
//_______ DiJet part _____________________________________________________
- if (nRecJets > 1)
+ if (nRecJetsCuts > 1) // OB: debugged this
{
AliAODJet* jet1 = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(0));
} // End if(nGenJets > 1)
+ // ____ efficiency _______________________________
+
+ // arrays for generated particles: reconstructed AOD track index, isPrimary flag
+ TArrayI indexAODTr;
+ TArrayS isGenPrim;
+
+ // array for reconcstructed AOD tracks: generated particle index
+ TArrayI indexMCTr;
+
+ Int_t nTracksAODMCCharged = GetListOfTracks(fTracksAODMCCharged, kTrackAODMCCharged);
+ if(fDebug>2)Printf("%s:%d selected AODMC tracks: %d ",(char*)__FILE__,__LINE__,nTracksAODMCCharged);
+
+ Int_t nTracksRecQualityCuts = GetListOfTracks(fTracksRecQualityCuts, kTrackAODQualityCuts);
+ if(fDebug>2)Printf("%s:%d selected rec tracks quality after cuts, full acceptance/pt : %d ",(char*)__FILE__,__LINE__,nTracksRecQualityCuts);
+
+ // associate gen and rec tracks, store indices in TArrays
+ AssociateGenRec(fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,indexMCTr,isGenPrim);
+
+ // single track eff
+ FillSingleTrackRecEffHisto(fhnSingleTrackRecEffHisto,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim);
+
+ // jet track eff
+ for(Int_t ij=0; ij<nRecEffJets; ++ij){
+
+ AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecEff->At(ij));
+
+ if(ij==0){ // leading jet
+
+ TList* jettracklist = new TList();
+ Double_t sumPt = 0.;
+
+ GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt); // for efficiency: gen tracks from pointing with gen/rec jet
+
+ Double_t jetEta = jet->Eta();
+ Double_t jetPhi = TVector2::Phi_0_2pi(jet->Phi());
+ Double_t jetPt = sumPt;
+
+ fQAJetHistosRecEffLeading->FillJetQA( jetEta, jetPhi, jetPt );
+ FillJetTrackRecEffHisto(fhnJetTrackRecEffHisto,jetPhi,jetEta,jetPt,jettracklist,fTracksAODMCCharged,indexAODTr,isGenPrim);
+
+ delete jettracklist;
+ }
+ }
+
+ //___________________
+
fTracksRec->Clear();
fTracksRecCuts->Clear();
fTracksGen->Clear();
+ fTracksAODMCCharged->Clear();
+ fTracksRecQualityCuts->Clear();
+
fJetsRec->Clear();
fJetsRecCuts->Clear();
fJetsGen->Clear();
+ fJetsRecEff->Clear();
//Post output data.
PostData(1, fCommonHistList);
{
jetBin = fDiJetJetInvMassMin + i*stepInvMass/2.;
if(((fDiJetJetInvMassMin+i*stepInvMass) <= invMass) &&
- (fDiJetJetInvMassMin + (i+1)*stepInvMass) > invMass) jetBinOk = jetBin;
+ (fDiJetJetInvMassMin + (i+1)*stepInvMass) > invMass) {jetBinOk = jetBin; break;}
+ else jetBinOk = -1.;
}
- jetBinOk = -1.;
}
else if (kindBins == 3)
{
{
jetBin = fDiJetJetPtMin + i*stepPt/2.;
if(((fDiJetJetPtMin+i*stepPt) <= EtMean) &&
- (fDiJetJetPtMin + (i+1)*stepPt) > EtMean) jetBinOk = jetBin;
+ (fDiJetJetPtMin + (i+1)*stepPt) > EtMean) {jetBinOk = jetBin; break;}
+ else jetBinOk = -1.;
}
- jetBinOk = -1.;
}
else if (kindBins == 2)
{
{
jetBin = fDiJetJetPtMin + i*stepPt/2.;
if(((fDiJetJetPtMin+i*stepPt) <= leadingJetPt) &&
- (fDiJetJetPtMin + (i+1)*stepPt) > leadingJetPt) jetBinOk = jetBin;
+ (fDiJetJetPtMin + (i+1)*stepPt) > leadingJetPt) {jetBinOk = jetBin; break;}
+ else jetBinOk = -1.;
}
- jetBinOk = -1.;
}
else {Printf("WARNING: kindBins wrongly set ! Please make sure to call SetKindSlices() and set the kind parameter to 1, 2 or 3.\n");}
if(type==kTrackUndef) return 0;
Int_t iCount = 0;
- if(type==kTrackAODCuts || type==kTrackAOD){
+ if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD){
// all rec. tracks, esd filter mask, eta range
if(!fAOD) return -1;
for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
AliAODTrack *tr = fAOD->GetTrack(it);
- if(type == kTrackAODCuts){
+ if(type == kTrackAODCuts || type==kTrackAODQualityCuts ){
if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue;
- if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
- if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
- if(tr->Pt() < fTrackPtCut) continue;
+ if(type == kTrackAODCuts){
+ if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
+ if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
+ if(tr->Pt() < fTrackPtCut) continue;
+ }
}
list->Add(tr);
iCount++;
list->Sort();
}
+
+// _ ________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::AssociateGenRec(TList* tracksAODMCCharged,TList* tracksRec, TArrayI& indexAODTr,TArrayI& indexMCTr,TArrayS& isGenPrim)
+{
+ // associate generated and reconstructed tracks, fill TArrays of list indices
+
+
+ Int_t nTracksRec = tracksRec->GetSize();
+ Int_t nTracksGen = tracksAODMCCharged->GetSize();
+ TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
+
+ if(!nTracksGen) return;
+ if(!tca) return;
+
+ // set size
+ indexAODTr.Set(nTracksGen);
+ indexMCTr.Set(nTracksRec);
+ isGenPrim.Set(nTracksGen);
+
+ indexAODTr.Reset(-1);
+ indexMCTr.Reset(-1);
+ isGenPrim.Reset(0);
+
+ // loop over reconstructed tracks, get generated track
+
+ for(Int_t iRec=0; iRec<nTracksRec; iRec++){
+
+ AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec));
+
+ Int_t label = rectrack->GetLabel();
+
+ // find MC track in our list
+ AliAODMCParticle* gentrack = 0x0;
+ if(label>=0) gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(label));
+
+ Int_t listIndex = -1;
+ if(gentrack) listIndex = tracksAODMCCharged->IndexOf(gentrack);
+
+ if(listIndex>=0){
+
+ indexAODTr[listIndex] = iRec;
+ indexMCTr[iRec] = listIndex;
+ }
+ }
+
+
+ // define primary sample for reconstruction efficiency
+
+ for(Int_t iGen=0; iGen<nTracksGen; iGen++){
+
+ AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tracksAODMCCharged->At(iGen));
+
+ Int_t pdg = gentrack->GetPdgCode();
+
+ // 211 - pi, 2212 - proton, 321 - Kaon, 11 - electron, 13 - muon
+ if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 ||
+ TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
+
+ isGenPrim[iGen] = kTRUE;
+ }
+ }
+}
+
+// _____________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::FillSingleTrackRecEffHisto(THnSparse* histo, TList* tracksGen, TList* tracksRec, TArrayI& indexAODTr, TArrayS& isGenPrim){
+
+ // fill THnSparse for single track reconstruction efficiency
+
+ Int_t nTracksGen = tracksGen->GetSize();
+
+ if(!nTracksGen) return;
+
+ for(Int_t iGen=0; iGen<nTracksGen; iGen++){
+
+ if(isGenPrim[iGen] != 1) continue; // select primaries
+
+ AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tracksGen->At(iGen));
+
+ Double_t ptGen = gentrack->Pt();
+ Double_t etaGen = gentrack->Eta();
+ Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
+
+ // apply same acc & pt cuts as for FF
+ // could in principle also be done setting THNsparse axis limits before projecting,
+ // but then the binning needs to be fine grained enough
+
+ if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
+ if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
+ if(ptGen < fTrackPtCut) continue;
+
+ Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track
+ Double_t isRec = 0;
+ Double_t ptRec = -1;
+
+ if(iRec>=0){
+
+ AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec));
+ ptRec = rectrack->Pt();
+ isRec = 1;
+ }
+
+ Double_t entries[5] = {phiGen,etaGen,ptGen,ptRec,isRec};
+ histo->Fill(entries);
+ }
+}
+
+// ______________________________________________________________________________________________________________________________________________________
+void AliAnalysisTaskFragmentationFunction::FillJetTrackRecEffHisto(THnSparse* histo,Double_t jetPhi, Double_t jetEta, Double_t jetPt, TList* jetTrackList,
+ TList* tracksGen, TArrayI& indexAODTr, TArrayS& isGenPrim)
+{
+ // fill THnSparse for jet track reconstruction efficiency
+
+ Int_t nTracksJet = jetTrackList->GetSize(); // list with AODMC tracks
+
+ if(!nTracksJet) return;
+
+ for(Int_t iTr=0; iTr<nTracksJet; iTr++){
+
+ AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (jetTrackList->At(iTr));
+
+ // find jet track in gen tracks list
+ Int_t iGen = tracksGen->IndexOf(gentrack);
+
+ if(iGen<0){
+ if(fDebug>0) Printf("%s:%d gen jet track not found ",(char*)__FILE__,__LINE__);
+ continue;
+ }
+
+ if(isGenPrim[iGen] != 1) continue; // select primaries
+
+ Double_t ptGen = gentrack->Pt();
+ Double_t etaGen = gentrack->Eta();
+ Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
+
+ // apply same acc & pt cuts as for FF
+ // could in principle also be done setting THNsparse axis limits before projecting,
+ // but then the binning needs to be fine grained enough
+
+ if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
+ if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
+ if(ptGen < fTrackPtCut) continue;
+
+ Double_t z = ptGen / jetPt;
+ Double_t xi = 0;
+ if(z>0) xi = TMath::Log(1/z);
+
+ Double_t isRec = 0;
+ Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track
+ if(iRec>=0) isRec = 1;
+
+ Double_t entries[7] = {jetPhi,jetEta,jetPt,ptGen,z,xi,isRec};
+ histo->Fill(entries);
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff