X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FEMCALJetTasks%2FUserTasks%2FAliAnalysisTaskChargedJetsPA.cxx;h=f9d8d446a4903c5e09ce1687f66960278d2a84e1;hb=2d6ed1a3a88b676028c00e59c56a3a740476379b;hp=0f5b67d172aab73dc11cd3f41dcd8738c3ff3490;hpb=fa866fc819f2817e8aaead709dcfd795b98c480d;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx index 0f5b67d172a..f9d8d446a49 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx @@ -2,7 +2,6 @@ #include #include #include -#include #include #include #include @@ -29,6 +28,9 @@ #include "AliAnalysisTaskSE.h" #endif +#include +#include "TFormula.h" +#include "AliESDtrackCuts.h" #include #include #include "AliGenPythiaEventHeader.h" @@ -42,11 +44,11 @@ #include "AliAODMCParticle.h" #include "AliAnalysisUtils.h" #include "AliRhoParameter.h" +#include "TVector3.h" #include "AliAnalysisTaskChargedJetsPA.h" using std::min; -//TODO: Not accessing the particles when using MC //TODO: FillHistogram can be done better with virtual TH1(?) ClassImp(AliAnalysisTaskChargedJetsPA) @@ -57,49 +59,65 @@ void AliAnalysisTaskChargedJetsPA::Init() AliInfo("Creating histograms."); #endif + SetCurrentOutputList(0); + + // Cuts TH1D* tmpHisto = AddHistogram1D("hNumberEvents", "Number of events (0 = before cuts, 1 = after cuts)", "", 2, 0, 2, "stage","N^{Events}/cut"); tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts"); tmpHisto->GetXaxis()->SetBinLabel(2, "After cuts"); - tmpHisto = AddHistogram1D("hEventAcceptance", "Accepted events (0 = before cuts, 1 = after pile up, 2 = after vertex)", "", 3, 0, 3, "stage","N^{Events}/cut"); tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts"); tmpHisto->GetXaxis()->SetBinLabel(2, "After pile up"); tmpHisto->GetXaxis()->SetBinLabel(3, "After vertex"); - tmpHisto = AddHistogram1D("hTrackAcceptance", "Accepted tracks (0 = before cuts, 1 = after eta, 2 = after pT)", "", 3, 0, 3, "stage","N^{Tracks}/cut"); tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts"); tmpHisto->GetXaxis()->SetBinLabel(2, "After eta"); tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}"); - tmpHisto = AddHistogram1D("hJetAcceptance", "Accepted jets (0 = before cuts, 1 = after eta, 2 = after pT, 3 = after area)", "", 4, 0, 4, "stage","N^{Jets}/cut"); tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts"); tmpHisto->GetXaxis()->SetBinLabel(2, "After eta"); tmpHisto->GetXaxis()->SetBinLabel(3, "After p_{T}"); tmpHisto->GetXaxis()->SetBinLabel(4, "After area"); - - // NOTE: Jet histograms - if (fAnalyzeJets) + TH2* tmpHisto2D = AddHistogram2D("hJetPtCutStages", "Jets p_{T} distribution", "", 500, -50., 200., 4, 0, 4, "p_{T} (GeV/c)","Cut stage","dN^{Jets}/dp_{T}"); + tmpHisto2D->GetYaxis()->SetBinLabel(1, "Before cuts"); + tmpHisto2D->GetYaxis()->SetBinLabel(2, "After eta"); + tmpHisto2D->GetYaxis()->SetBinLabel(3, "After p_{T}"); + tmpHisto2D->GetYaxis()->SetBinLabel(4, "After area"); + + AddHistogram1D("hVertexX", "X distribution of the vertex", "", 2000, -1., 1., "#Delta x(cm)","dN^{Events}/dx"); + AddHistogram1D("hVertexY", "Y distribution of the vertex", "", 2000, -1., 1., "#Delta y(cm)","dN^{Events}/dy"); + AddHistogram2D("hVertexXY", "XY distribution of the vertex", "COLZ", 500, -1., 1., 500, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy"); + AddHistogram1D("hVertexZBeforeVertexCut", "Z distribution of the vertex (before std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz"); + AddHistogram1D("hVertexZAfterVertexCut", "Z distribution of the vertex (after std. vertex cut)", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz"); + AddHistogram1D("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr"); + AddHistogram1D("hCentralityV0M", "Centrality distribution V0M", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); + AddHistogram1D("hCentralityV0A", "Centrality distribution V0A", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); + AddHistogram1D("hCentralityV0C", "Centrality distribution V0C", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); + AddHistogram1D("hCentralityZNA", "Centrality distribution ZNA", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); + AddHistogram1D("hCentrality", Form("Centrality distribution %s", fCentralityType.Data()), "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); + + if(fDoJetAnalysis) { - // ######## Jet spectra - AddHistogram1D("hRawJetPt", "Raw jets p_{T} distribution (before cuts)", "", 500, 0., 250., "p_{T} (GeV/c)", "dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPt", "Jets p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPtBgrdSubtractedKTImprovedCMS", "Jets p_{T} distribution, KT background (Improved CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - + // Background corrected jet spectra + AddHistogram2D("hJetPtBgrdSubtractedExternal", "Jets p_{T} distribution, external bgrd. subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedPP", "Jets p_{T} distribution, pp background subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedExternal_Phi1", "Jets p_{T} distribution, external background (Improved CMS) subtracted (1st part of azimuth)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedExternal_Phi2", "Jets p_{T} distribution, external background (Improved CMS) subtracted (2nd part of azimuth)", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTImprovedCMS", "Jets p_{T} distribution, KT background (Improved CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background (Cone R=0.6 around jets excluded) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTPbPb", "Jets p_{T} distribution, KT background (PbPb w/o ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTPbPbWithGhosts", "Jets p_{T} distribution, KT background (PbPb w/ ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTCMS", "Jets p_{T} distribution, KT background (CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTMean", "Jets p_{T} distribution, KT background (Mean) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTTrackLike", "Jets p_{T} distribution, KT background (track-like) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + + AddHistogram2D("hJetPtSubtractedRhoExternal", "Mean subtracted KT (External) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean"); AddHistogram2D("hJetPtSubtractedRhoKTImprovedCMS", "Mean subtracted KT (CMS w/o signal) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean"); - AddHistogram2D("hJetPtSubtractedRhoKTImprovedCMS020", "Mean subtracted KT (CMS w/o signal) background from jets, 0-20", "COLZ", 600, 0, 150, 400,0.,40., "Jet p_{T} (GeV/c)", "#rho (GeV/c)", "dN^{Events}/dp_{T}#rho"); - - if(fAnalyzeDeprecatedBackgrounds) - { - AddHistogram2D("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background (Cone R=0.6 around jets excluded) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPtBgrdSubtractedKTPbPb", "Jets p_{T} distribution, KT background (PbPb w/o ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPtBgrdSubtractedKTPbPbWithGhosts", "Jets p_{T} distribution, KT background (PbPb w/ ghosts) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPtBgrdSubtractedKTCMS", "Jets p_{T} distribution, KT background (CMS) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPtBgrdSubtractedKTMean", "Jets p_{T} distribution, KT background (Mean) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPtBgrdSubtractedKTTrackLike", "Jets p_{T} distribution, KT background (track-like) subtracted", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - } + AddHistogram2D("hJetPtSubtractedRhoPP", "Mean subtracted KT (pp from Michal) background from jets", "COLZ", 600, 0, 150, fNumberOfCentralityBins, 0, 100, "Jet p_{T}", "Centrality", "#rho mean"); - // ######## Jet stuff + // Jet QA plots AddHistogram2D("hJetConstituentPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}"); + AddHistogram2D("hJetConstituentPtVsJetPt", "Jet constituents p_{T} distribution", "", 500, -50., 200., 200, 0, 200, "#it{p}_{T} (GeV/c)","#it{p}_{T}^{jet} (GeV/c)","dN^{Tracks}/dp_{T}"); AddHistogram1D("hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}"); AddHistogram1D("hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}"); AddHistogram2D("hJetCount", "Correlation jets/accepted jets", "", 200, 0., 200., 200, 0., 200., "N jets","N jets accepted", "d^{2}N^{Events}/dN^{Jets dN^{Jets, acc}}"); @@ -110,79 +128,37 @@ void AliAnalysisTaskChargedJetsPA::Init() AddHistogram1D("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)"); AddHistogram1D("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)"); - // ########## Dijet stuff - AddHistogram1D("hDijetLeadingJetPt", "Dijet leading jet p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hDijetConstituentsPt", "Dijet constituents p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram2D("hDijetPtCorrelation", "Dijet constituents p_{T} correlation", "COLZ", 500, 5., 100., 500, 5., 100., "1st leading jet p_{T} (GeV/c)","2nd leading jet p_{T} (GeV/c)","dN^{Dijets}/d^{2}p_{T}"); - } - - // NOTE: Jet background histograms - if (fAnalyzeBackground) - { - // ########## Default background estimates + // Background distributions + AddHistogram2D("hKTBackgroundExternal", "KT background density (External task)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackgroundExternal20GeV", "KT background density (External task, jet p_{T} > 20 GeV)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); AddHistogram2D("hKTBackgroundImprovedCMS", "KT background density (Improved CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundImprovedCMSExternal", "KT background density (Improved CMS approach from external task)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hDeltaPtKTImprovedCMS", "Background fluctuations #delta p_{T} (KT, Improved CMS-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTImprovedCMSPartialExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, partial jet exclusion)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTImprovedCMSPartialExclusion_Signal", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, partial jet exclusion w/ 1/N_{sig} probability)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTImprovedCMSFullExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, full leading jet exclusion)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtNoBackground", "Background fluctuations #delta p_{T} (No background)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtNoBackgroundNoEmptyCones", "Background fluctuations #delta p_{T} (No background, no empty cones)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - - AddHistogram1D("hKTMeanBackgroundImprovedCMS", "KT background mean (Improved CMS approach)", "", 100, 0, 100, "Centrality", "#rho mean"); - - AddHistogram2D("hDijetBackground", "Background density (dijets excluded)", "", 200, 0., 20., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 200, 0., 20., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - if(fAnalyzeDeprecatedBackgrounds) - { - // ########## Different background estimates - AddHistogram2D("hKTBackgroundPbPb", "KT background density (PbPb approach, no ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundPbPbWithGhosts", "KT background density (PbPb approach w/ ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundCMS", "KT background density (CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundMean", "KT background density (Mean approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundTrackLike", "KT background density (Track-like approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - AddHistogram2D("hTRBackgroundNoExcl", "TR background density (No signal excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hTRBackgroundCone02", "TR background density (Cones R=0.2 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hTRBackgroundCone04", "TR background density (Cones R=0.4 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hTRBackgroundCone06", "TR background density (Cones R=0.6 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hTRBackgroundCone08", "TR background density (Cones R=0.8 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram2D("hTRBackgroundExact", "TR background density (signal jets exactly excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - // ########## Delta Pt - AddHistogram2D("hDeltaPtKTPbPb", "Background fluctuations #delta p_{T} (KT, PbPb w/o ghosts)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTPbPbWithGhosts", "Background fluctuations #delta p_{T} (KT, PbPb w/ ghosts)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTCMS", "Background fluctuations #delta p_{T} (KT, CMS-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTMean", "Background fluctuations #delta p_{T} (KT, Mean)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtKTTrackLike", "Background fluctuations #delta p_{T} (KT, track-like)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - AddHistogram2D("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR, cone R=0.6)", "", 1201, -40.0, 40.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); - - // ########## Profiles for background means vs. centrality - AddHistogram1D("hKTMeanBackgroundPbPb", "KT background mean (PbPb approach w/o ghosts)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean"); - AddHistogram1D("hKTMeanBackgroundPbPbWithGhosts", "KT background mean (PbPb approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean"); - AddHistogram1D("hKTMeanBackgroundCMS", "KT background mean (CMS approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean"); - AddHistogram1D("hKTMeanBackgroundMean", "KT background mean (Mean approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean"); - AddHistogram1D("hKTMeanBackgroundTPC", "KT background mean (Track-like approach)", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean"); - AddHistogram1D("hTRMeanBackground", "TR background mean", "", fNumberOfCentralityBins, 0, 100, "Centrality", "#rho mean"); - } - } - - - // NOTE: Track & Cluster & QA histograms - if (fAnalyzeQA) - { - AddHistogram1D("hVertexX", "X distribution of the vertex", "", 2000, -1., 1., "#Delta x(cm)","dN^{Events}/dx"); - AddHistogram1D("hVertexY", "Y distribution of the vertex", "", 2000, -1., 1., "#Delta y(cm)","dN^{Events}/dy"); - AddHistogram2D("hVertexXY", "XY distribution of the vertex", "COLZ", 500, -1., 1., 500, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy"); - AddHistogram1D("hVertexZ", "Z distribution of the vertex", "", 200, -20., 20., "#Delta z(cm)","dN^{Events}/dz"); - AddHistogram1D("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr"); - AddHistogram1D("hCentralityV0M", "Centrality distribution V0M", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); - AddHistogram1D("hCentralityV0A", "Centrality distribution V0A", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); - AddHistogram1D("hCentralityV0C", "Centrality distribution V0C", "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); - AddHistogram1D("hCentrality", Form("Centrality distribution %s", fCentralityType.Data()), "", fNumberOfCentralityBins, 0., 100., "Centrality","dN^{Events}"); - - + AddHistogram2D("hPPBackground", "PP background density (Michals approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackgroundPbPb", "KT background density (PbPb approach, no ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackgroundPbPbWithGhosts", "KT background density (PbPb approach w/ ghosts)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackgroundCMS", "KT background density (CMS approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackgroundMean", "KT background density (Mean approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackgroundTrackLike", "KT background density (Track-like approach)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackgroundNoExcl", "TR background density (No signal excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackgroundCone02", "TR background density (Cones R=0.2 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackgroundCone04", "TR background density (Cones R=0.4 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackgroundCone06", "TR background density (Cones R=0.6 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackgroundCone08", "TR background density (Cones R=0.8 around signal jets excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackgroundExact", "TR background density (signal jets exactly excluded)", "LEGO2", 400, 0., 40., fNumberOfCentralityBins, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + + // Delta pt distributions + AddHistogram2D("hDeltaPtExternalBgrd", "Background fluctuations #delta p_{T} (KT, External)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtExternalBgrdVsPt", "Background fluctuations #delta p_{T} (KT, External, in p_{T} bins)", "", 1801, -40.0, 80.0, 200, 0, 200, "#delta p_{T} (GeV/c)","Raw jet p_{T}","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTImprovedCMS", "Background fluctuations #delta p_{T} (KT, Improved CMS-like)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTImprovedCMSFullExclusion", "Background fluctuations #delta p_{T} (KT, Improved CMS-like, full leading jet exclusion)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtNoBackground", "Background fluctuations #delta p_{T} (No background)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTPbPb", "Background fluctuations #delta p_{T} (KT, PbPb w/o ghosts)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTPbPbWithGhosts", "Background fluctuations #delta p_{T} (KT, PbPb w/ ghosts)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTCMS", "Background fluctuations #delta p_{T} (KT, CMS-like)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTMean", "Background fluctuations #delta p_{T} (KT, Mean)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTTrackLike", "Background fluctuations #delta p_{T} (KT, track-like)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR, cone R=0.6)", "", 1801, -40.0, 80.0, fNumberOfCentralityBins, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + + // Track QA plots AddHistogram2D("hTrackCountAcc", "Number of tracks in acceptance vs. centrality", "LEGO2", 750, 0., 750., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}"); AddHistogram2D("hTrackPt", "Tracks p_{T} distribution", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}"); AddHistogram2D("hTrackPtNegEta", "Tracks p_{T} distribution (negative #eta)", "", 1000, 0., 250., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Tracks}/dp_{T}"); @@ -190,41 +166,74 @@ void AliAnalysisTaskChargedJetsPA::Init() AddHistogram1D("hTrackCharge", "Charge", "", 11, -5, 5, "Charge (e)","dN^{Tracks}/dq"); AddHistogram1D("hTrackPhi", "Track #phi distribution", "", 360, 0, TMath::TwoPi(), "#phi","dN^{Tracks}/d#phi"); AddHistogram2D("hTrackPhiEta", "Track angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dN^{Tracks}/(d#phi d#eta)"); - + AddHistogram2D("hTrackPtPhiEta", "Track p_{T} angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dp_{T}^{Tracks}/(d#phi d#eta)"); AddHistogram2D("hTrackPhiPtCut", "Track #phi distribution for different pT cuts", "LEGO2", 360, 0, TMath::TwoPi(), 20, 0, 20, "#phi", "p_{T} lower cut", "dN^{Tracks}/d#phi dp_{T}"); - AddHistogram2D("hTrackPhiLabel", "Track #phi distribution for different labels", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi"); AddHistogram2D("hTrackPhiTrackType", "Track #phi distribution for different track types", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi"); - AddHistogram2D("hTrackEta", "Track #eta distribution", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Tracks}/d#eta"); - if (fAnalyzeJets) + AddHistogram2D("hTrackPtTrackType", "Track p_{T} distribution for different track types", "LEGO2", 1000, 0., 250., 3, 0, 3, "p_{T} (GeV/c)", "Label", "dN^{Tracks}/dp_{T}"); + AddHistogram2D("hTrackEta", "Track #eta distribution", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Tracks}/d#eta"); + + // Jet QA plots + AddHistogram1D("hRawJetArea", "Jets area distribution w/o area cut", "", 200, 0., 2., "Area","dN^{Jets}/dA"); + AddHistogram2D("hJetArea", "Jets area distribution", "COLZ", 200, 0., 2., 500, -50., 200, "Area","Jet p_{T}","dN^{Jets}/dA"); + AddHistogram2D("hRawJetPhiEta", "Raw Jets angular distribution w/o #eta cut", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); + AddHistogram2D("hJetEta", "Jets #eta distribution", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta"); + AddHistogram2D("hJetEta2GeVTracks", "Jets #eta distribution, track p_{T} > 2 GeV", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta"); + AddHistogram2D("hJetEta4GeVTracks", "Jets #eta distribution, track p_{T} > 4 GeV", "COLZ", 180, fMinEta, fMaxEta, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta"); + AddHistogram2D("hJetPhiEta", "Jets angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); + AddHistogram2D("hJetPtPhiEta", "Jets p_{T} angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dp_{T}^{Jets}/(d#phi d#eta)"); + AddHistogram2D("hJetPtVsConstituentCount", "Jets number of constituents vs. jet p_{T}", "COLZ", 400, 0., 200., 100, 0., 100., "p_{T}","N^{Tracks}","dN^{Jets}/(dp_{T} dN^{tracks})"); + + // ######## Jet profiles + if(fAnalyzeJetProfile) { - // ######## Jet QA - AddHistogram1D("hRawJetArea", "Jets area distribution w/o area cut", "", 200, 0., 2., "Area","dN^{Jets}/dA"); - AddHistogram1D("hJetArea", "Jets area distribution", "", 200, 0., 2., "Area","dN^{Jets}/dA"); - AddHistogram2D("hRawJetPhiEta", "Raw Jets angular distribution w/o #eta cut", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); - AddHistogram2D("hJetEta", "Jets #eta distribution", "COLZ", 180, -fTrackEtaWindow, +fTrackEtaWindow, fNumberOfCentralityBins, 0., 100., "#eta", "Centrality", "dN^{Jets}/d#eta"); - AddHistogram2D("hJetPhiEta", "Jets angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -1.0, 1.0, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); - AddHistogram2D("hJetPtVsConstituentCount", "Jets number of constituents vs. jet p_{T}", "COLZ", 400, 0., 200., 100, 0., 100., "p_{T}","N^{Tracks}","dN^{Jets}/(dp_{T} dN^{tracks})"); + SetCurrentOutputList(1); + AddHistogram2D("hJetProfile10GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 10 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + AddHistogram2D("hJetProfile20GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 20 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + AddHistogram2D("hJetProfile30GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 30 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + AddHistogram2D("hJetProfile40GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 40 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + AddHistogram2D("hJetProfile50GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 50 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + AddHistogram2D("hJetProfile60GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 60 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + AddHistogram2D("hJetProfile70GeV", "Jet profile, cone p_{T}/jet p_{T} vs. jet radius, jet p_{T} > 70 GeV", "", 12, 0, 0.6,200, 0., 2., "Cone radius","dN^{Jets}/dR", "Ratio"); + SetCurrentOutputList(0); } } - - // NOTE: Pythia histograms - if (fAnalyzePythia) + // ######## Jet track cuts + if(fAnalyzeTrackcuts) { - AddHistogram1D("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}"); + SetCurrentOutputList(2); + + AddCutHistogram("hCutsNumberClusters", "Trackcut histogram: Number of clusters", "Number of clusters", 40, 20, 160); + AddCutHistogram("hCutsChi2TPC", "Trackcut histogram: #chi^{2} per TPC cluster", "#chi^{2}", 40, 0, 8); + AddCutHistogram("hCutsChi2ITS", "Trackcut histogram: #chi^{2} per ITS cluster", "#chi^{2}", 25, 0., 50); + AddCutHistogram("hCutsChi2Constrained", "Trackcut histogram: #chi^{2} for global constrained tracks", "#chi^{2}", 60, 0, 60); + AddCutHistogram("hCutsDCAZ", "Trackcut histogram: Max. DCA z for prim. vertex", "DCA z", 20, 0, 4); + AddCutHistogram("hCutsSPDHit", "Trackcut histogram: Hit in SPD layer", "Hit or not", 2, -0.5, 1.5); + AddCutHistogram("hCutsNumberCrossedRows", "Trackcut histogram: Number of crossed rows", "Number of crossed rows", 40, 20, 160); + AddCutHistogram("hCutsNumberCrossedRowsOverFindableClusters", "Trackcut histogram: Number of crossed rows over findable clusters", "Number of crossed rows over findable clusters", 26, 0.4, 1.8); + AddCutHistogram("hCutsSharedTPC", "Trackcut histogram: Shared TPC clusters", "Shared fraction", 40, 0, 1); + AddCutHistogram("hCutsTPCRefit", "Trackcut histogram: TPC refit", "Has TPC refit", 2, -0.5, 1.5); + AddCutHistogram("hCutsTPCLength", "Trackcut histogram: TPC length", "TPC length", 40, 0, 170); + AddCutHistogram("hCutsTrackConstrained", "Trackcut histogram: Tracks constrained to vertex", "Track is constrained", 2, -0.5, 1.5); + AddCutHistogram("hCutsClustersPtDependence", "Trackcut histogram: Use pT dependence for Number of clusters cut", "Use pT dependence", 2, -0.5, 1.5); + + SetCurrentOutputList(0); } - // register Histograms - for (Int_t i = 0; i < fHistCount; i++) + PostData(1, fOutputLists[0]); + if(fAnalyzeJetProfile) + PostData(2, fOutputLists[1]); + if(fAnalyzeTrackcuts) { - fOutputList->Add(fHistList->At(i)); + if(fAnalyzeJetProfile) + PostData(3, fOutputLists[2]); + else + PostData(2, fOutputLists[1]); } - - PostData(1,fOutputList); // important for merging } //________________________________________________________________________ -AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeJets(1), fAnalyzeQA(1), fAnalyzeBackground(1), fAnalyzeDeprecatedBackgrounds(1), fAnalyzePythia(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsKinematics(0), fUseVertexCut(1), fUsePileUpCut(1), fSetCentralityToOne(0), fPartialAnalysisNParts(1), fPartialAnalysisIndex(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fUsePtHardBin(-1), fRhoTaskName(), fNcoll(6.88348), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.6), fNumberRandCones(8), fNumberExcludedJets(-1), fDijetMaxAngleDeviation(10.0), fPhysicalJetRadius(0.6), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetPt(1.0), fMinJetArea(0.5), fMinBackgroundJetPt(0.0), fMinDijetLeadingPt(10.0), fNumberOfCentralityBins(20), fCentralityType("V0A"), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0), fIsDEBUG(0), fEventCounter(0) +AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName, Bool_t analyzeJetProfile, Bool_t analyzeTrackcuts) : AliAnalysisTaskSE(name), fOutputLists(), fCurrentOutputList(0), fDoJetAnalysis(1), fAnalyzeJetProfile(0), fAnalyzeTrackcuts(0), fParticleLevel(0), fUseDefaultVertexCut(1), fUsePileUpCut(1), fSetCentralityToOne(0), fNoExternalBackground(0), fBackgroundForJetProfile(0), fPartialAnalysisNParts(1), fPartialAnalysisIndex(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(), fTrackArrayName(), fBackgroundJetArrayName(), fRhoTaskName(), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fNumberExcludedJets(-1), fMinEta(-0.9), fMaxEta(0.9), fMinJetEta(-0.5), fMaxJetEta(0.5), fMinTrackPt(0.150), fMinJetPt(5.0), fMinJetArea(0.5), fMinBackgroundJetPt(0.0), fMinNCrossedRows(70), fNumberOfCentralityBins(20), fCentralityType("V0A"), fPrimaryVertex(0), fFirstLeadingJet(0), fSecondLeadingJet(0), fFirstLeadingKTJet(0), fSecondLeadingKTJet(0), fNumberSignalJets(0), fNumberSignalJetsAbove5GeV(0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fIsDEBUG(0), fIsPA(1), fEventCounter(0), fHybridESDtrackCuts(0), fHybridESDtrackCuts_noPtDep(0) { #ifdef DEBUGMODE AliInfo("Calling constructor."); @@ -235,144 +244,364 @@ AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, con fTaskInstanceCounter = instance; instance++; - fTrackArrayName = new TString(trackArrayName); - if (fTrackArrayName->Contains("MCParticles") || fTrackArrayName->Contains("mcparticles")) - fIsKinematics = kTRUE; + fAnalyzeJetProfile = analyzeJetProfile; + fAnalyzeTrackcuts = analyzeTrackcuts; - fJetArrayName = new TString(jetArrayName); - if (strcmp(fJetArrayName->Data(),"") == 0) - fAnalyzeJets = kFALSE; - else - fAnalyzeJets = kTRUE; - - fBackgroundJetArrayName = new TString(backgroundJetArrayName); - if (strcmp(fBackgroundJetArrayName->Data(),"") == 0) - fAnalyzeBackground = kFALSE; - else - fAnalyzeBackground = kTRUE; + // Save the observables array names + fTrackArrayName = trackArrayName; + fJetArrayName = jetArrayName; + fBackgroundJetArrayName = backgroundJetArrayName; - DefineOutput(1, TList::Class()); - - fHistList = new TList(); - - for(Int_t i=0;i<1024;i++) - fSignalJets[i] = NULL; + if (fTrackArrayName.Contains("MCParticles") || fTrackArrayName.Contains("mcparticles")) + fParticleLevel = kTRUE; + DefineOutput(1, TList::Class()); + if(fAnalyzeJetProfile) + DefineOutput(2, TList::Class()); + if(fAnalyzeTrackcuts) + { + if(fAnalyzeJetProfile) + DefineOutput(3, TList::Class()); + else + DefineOutput(2, TList::Class()); + } #ifdef DEBUGMODE AliInfo("Constructor done."); #endif - } //________________________________________________________________________ -inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius) +void AliAnalysisTaskChargedJetsPA::InitializeTrackcuts() { - Double_t tmpConePt = 0.0; - - for (Int_t i = 0; i < fTrackArray->GetEntries(); i++) + AliESDtrackCuts* commonTrackCuts = new AliESDtrackCuts; + commonTrackCuts->SetMaxChi2PerClusterTPC(4); + commonTrackCuts->SetMaxChi2PerClusterITS(36); + commonTrackCuts->SetAcceptKinkDaughters(kFALSE); + commonTrackCuts->SetRequireTPCRefit(kTRUE); + commonTrackCuts->SetRequireITSRefit(kTRUE); + commonTrackCuts->SetRequireSigmaToVertex(kFALSE); + commonTrackCuts->SetMaxDCAToVertexXY(2.4); + commonTrackCuts->SetMaxDCAToVertexZ(3.2); + commonTrackCuts->SetDCAToVertex2D(kTRUE); + commonTrackCuts->SetMaxFractionSharedTPCClusters(0.4); + commonTrackCuts->SetMaxChi2TPCConstrainedGlobal(36); + commonTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); + + AliESDtrackCuts* fTrackCutsPA_global = NULL; + AliESDtrackCuts* fTrackCutsPA_complementary = NULL; + AliESDtrackCuts* fTrackCutsPP_global = NULL; + AliESDtrackCuts* fTrackCutsPP_complementary = NULL; + AliESDtrackCuts* fTrackCutsPP_global_noPtDep = NULL; + AliESDtrackCuts* fTrackCutsPP_complementary_noPtDep = NULL; + + //pPb + fTrackCutsPA_global = static_cast(commonTrackCuts->Clone("fTrackCutsPA_global")); + fTrackCutsPA_global->SetMinNCrossedRowsTPC(fMinNCrossedRows); + fTrackCutsPA_global->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8); + fTrackCutsPA_complementary = static_cast(fTrackCutsPA_global->Clone("fTrackCutsPA_complementary")); + fTrackCutsPA_complementary->SetRequireITSRefit(kFALSE); + fTrackCutsPA_complementary->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff); + + //pp + fTrackCutsPP_global = static_cast(commonTrackCuts->Clone("fTrackCutsPP_global")); + TFormula *f1NClustersTPCLinearPtDep = new TFormula("f1NClustersTPCLinearPtDep","70.+30./20.*x"); + fTrackCutsPP_global->SetMinNClustersTPCPtDep(f1NClustersTPCLinearPtDep,20.); + fTrackCutsPP_global->SetMinNClustersTPC(70); + fTrackCutsPP_global->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1 + fTrackCutsPP_global->SetEtaRange(-0.9,0.9); + fTrackCutsPP_global->SetPtRange(0.15, 1e15); + fTrackCutsPP_complementary = static_cast(fTrackCutsPP_global->Clone("fTrackCutsPP_complementary")); + fTrackCutsPP_complementary->SetRequireITSRefit(kFALSE); + fTrackCutsPP_complementary->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff); + + //pp, no pT dependence + + fTrackCutsPP_global_noPtDep = static_cast(commonTrackCuts->Clone("fTrackCutsPP_global_noPtDep")); + fTrackCutsPP_global_noPtDep->SetMinNClustersTPC(70); + fTrackCutsPP_global_noPtDep->SetRequireTPCStandAlone(kTRUE); //cut on NClustersTPC and chi2TPC Iter1 + fTrackCutsPP_global_noPtDep->SetEtaRange(-0.9,0.9); + fTrackCutsPP_global_noPtDep->SetPtRange(0.15, 1e15); + fTrackCutsPP_complementary_noPtDep = static_cast(fTrackCutsPP_global_noPtDep->Clone("fTrackCutsPP_complementary_noPtDep")); + fTrackCutsPP_complementary_noPtDep->SetRequireITSRefit(kFALSE); + fTrackCutsPP_complementary_noPtDep->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff); + + fHybridESDtrackCuts = new AliESDHybridTrackcuts(); + fHybridESDtrackCuts_noPtDep = new AliESDHybridTrackcuts(); + if(fIsPA) { - AliVTrack* tmpTrack = static_cast(fTrackArray->At(i)); - if (IsTrackInAcceptance(tmpTrack)) - if(IsTrackInCone(tmpTrack, eta, phi, radius)) - tmpConePt = tmpConePt + tmpTrack->Pt(); + fHybridESDtrackCuts->SetMainCuts(fTrackCutsPA_global); + fHybridESDtrackCuts->SetAdditionalCuts(fTrackCutsPA_complementary); + fHybridESDtrackCuts_noPtDep->SetMainCuts(fTrackCutsPA_global); + fHybridESDtrackCuts_noPtDep->SetAdditionalCuts(fTrackCutsPA_complementary); + } + else + { + fHybridESDtrackCuts->SetMainCuts(fTrackCutsPP_global); + fHybridESDtrackCuts->SetAdditionalCuts(fTrackCutsPP_complementary); + fHybridESDtrackCuts_noPtDep->SetMainCuts(fTrackCutsPP_global_noPtDep); + fHybridESDtrackCuts_noPtDep->SetAdditionalCuts(fTrackCutsPP_complementary_noPtDep); } - return tmpConePt; -} + delete commonTrackCuts; +} //________________________________________________________________________ -inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard() +void AliAnalysisTaskChargedJetsPA::CreateCutHistograms() { - #ifdef DEBUGMODE - AliInfo("Starting GetPtHard."); - #endif - AliGenPythiaEventHeader* pythiaHeader = dynamic_cast(MCEvent()->GenEventHeader()); - if (MCEvent()) - if (!pythiaHeader) - { - // Check if AOD - AliAODMCHeader* aodMCH = dynamic_cast(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); - if (aodMCH) - { - for(UInt_t i = 0;iGetNCocktailHeaders();i++) - { - pythiaHeader = dynamic_cast(aodMCH->GetCocktailHeader(i)); - if (pythiaHeader) break; - } - } + AliESDEvent* fESD = dynamic_cast( InputEvent() ); + if (!fESD) + { + AliError("For cut analysis, ESDs must be processed!"); + return; + } + + + Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut + for (Int_t i=0;i < fESD->GetNumberOfTracks(); i++) + { + AliESDtrack* track = fESD->GetTrack(i); + + // Basics kinematic variables + Double_t pT = track->Pt(); + Double_t eta = track->Eta(); + Double_t phi = track->Phi(); + + // Number of clusters + Double_t nclsTPC = track->GetTPCncls(); + Double_t nclsITS = track->GetITSclusters(0); + + // Crossed rows + Double_t ncrTPC = track->GetTPCCrossedRows(); + Double_t nCRoverFC = 0; + if(track->GetTPCNclsF()) + nCRoverFC = track->GetTPCCrossedRows()/track->GetTPCNclsF(); + + // Chi2 of tracks + Double_t chi2ITS = 999.; + if (nclsITS) + chi2ITS = track->GetITSchi2()/nclsITS; + Double_t chi2TPC = 999.; + if (nclsTPC) + chi2TPC = track->GetTPCchi2()/nclsTPC; + Double_t chi2TPCConstrained = track->GetChi2TPCConstrainedVsGlobal(static_cast(fPrimaryVertex)); + + // Misc + Double_t SharedTPCClusters = 999.; + Double_t nClustersTPC = 0; + if(fHybridESDtrackCuts->GetMainCuts()->GetRequireTPCStandAlone()) { + nClustersTPC = track->GetTPCNclsIter1(); + } + else { + nClustersTPC = track->GetTPCclusters(0); } - #ifdef DEBUGMODE - AliInfo("Ending GetPtHard."); - #endif - if (pythiaHeader) - return pythiaHeader->GetPtHard(); + if(track->GetTPCncls()) + SharedTPCClusters = static_cast(track->GetTPCnclsS())/static_cast(nClustersTPC); + Double_t tpcLength = 0.; + if (track->GetInnerParam() && track->GetESDEvent()) { + tpcLength = track->GetLengthInActiveZone(1, 1.8, 220, track->GetESDEvent()->GetMagneticField()); + } + track->GetImpactParameters(dca, cov); - AliWarning(Form("In task %s: GetPtHard() failed!", GetName())); - return -1.0; -} + // Basic kinematic cuts + if((pT<0.15) || (TMath::Abs(eta)>0.9)) + continue; + SetCurrentOutputList(2); -//________________________________________________________________________ -inline Double_t AliAnalysisTaskChargedJetsPA::GetPythiaTrials() -{ - #ifdef DEBUGMODE - AliInfo("Starting GetPythiaTrials."); - #endif - AliGenPythiaEventHeader* pythiaHeader = dynamic_cast(MCEvent()->GenEventHeader()); - if (MCEvent()) - if (!pythiaHeader) - { - // Check if AOD - AliAODMCHeader* aodMCH = dynamic_cast(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName())); + Int_t trackType = 0; - if (aodMCH) - { - for(UInt_t i = 0;iGetNCocktailHeaders();i++) - { - pythiaHeader = dynamic_cast(aodMCH->GetCocktailHeader(i)); - if (pythiaHeader) break; - } - } + // ################################################################ + // ################################################################ + Int_t minNclsTPC = fHybridESDtrackCuts->GetMainCuts()->GetMinNClusterTPC(); + Int_t minNclsTPC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinNClusterTPC(); + fHybridESDtrackCuts->GetMainCuts()->SetMinNClustersTPC(0); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNClustersTPC(0); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsNumberClusters", nclsTPC, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMinNClustersTPC(minNclsTPC); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNClustersTPC(minNclsTPC_Additional); + // ################################################################ + // ################################################################ + Float_t maxChi2 = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2PerClusterTPC(); + Float_t maxChi2_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2PerClusterTPC(); + fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterTPC(999.); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterTPC(999.); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsChi2TPC", chi2TPC, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterTPC(maxChi2); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterTPC(maxChi2_Additional); + + // ################################################################ + // ################################################################ + Float_t maxChi2TPCConstrained = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2TPCConstrainedGlobal(); + Float_t maxChi2TPCConstrained_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2TPCConstrainedGlobal(); + fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2TPCConstrainedGlobal(999.); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2TPCConstrainedGlobal(999.); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsChi2Constrained", chi2TPCConstrained, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2TPCConstrainedGlobal(maxChi2TPCConstrained); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2TPCConstrainedGlobal(maxChi2TPCConstrained_Additional); + + // ################################################################ + // ################################################################ + Float_t maxDcaZ = fHybridESDtrackCuts->GetMainCuts()->GetMaxDCAToVertexZ(); + Float_t maxDcaZ_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxDCAToVertexZ(); + fHybridESDtrackCuts->GetMainCuts()->SetMaxDCAToVertexZ(999.); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxDCAToVertexZ(999.); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsDCAZ", TMath::Abs(dca[1]), pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMaxDCAToVertexZ(maxDcaZ); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxDCAToVertexZ(maxDcaZ_Additional); + // ################################################################ + // ################################################################ + AliESDtrackCuts::ITSClusterRequirement clusterReq = fHybridESDtrackCuts->GetMainCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD); + AliESDtrackCuts::ITSClusterRequirement clusterReq_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD); + fHybridESDtrackCuts->GetMainCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff); + fHybridESDtrackCuts->GetAdditionalCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff); + + Int_t hasPoint = 0; + if (track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1)) hasPoint = 1; + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsSPDHit", hasPoint, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, clusterReq); + fHybridESDtrackCuts->GetAdditionalCuts()->SetClusterRequirementITS(AliESDtrackCuts::kSPD, clusterReq_Additional); + + // ################################################################ + // ################################################################ + Float_t minNcrTPC = fHybridESDtrackCuts->GetMainCuts()->GetMinNCrossedRowsTPC(); + Float_t minNcrTPC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinNCrossedRowsTPC(); + fHybridESDtrackCuts->GetMainCuts()->SetMinNCrossedRowsTPC(0); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNCrossedRowsTPC(0); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsNumberCrossedRows", ncrTPC, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMinNCrossedRowsTPC(minNcrTPC); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinNCrossedRowsTPC(minNcrTPC_Additional); + + // ################################################################ + // ################################################################ + Float_t minCRoverFC = fHybridESDtrackCuts->GetMainCuts()->GetMinRatioCrossedRowsOverFindableClustersTPC(); + Float_t minCRoverFC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinRatioCrossedRowsOverFindableClustersTPC(); + fHybridESDtrackCuts->GetMainCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(0.); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(0.); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsNumberCrossedRowsOverFindableClusters", nCRoverFC, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinRatioCrossedRowsOverFindableClustersTPC(minCRoverFC_Additional); + + // ################################################################ + // ################################################################ + Float_t maxSharedTPC = fHybridESDtrackCuts->GetMainCuts()->GetMaxFractionSharedTPCClusters(); + Float_t maxSharedTPC_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxFractionSharedTPCClusters(); + fHybridESDtrackCuts->GetMainCuts()->SetMaxFractionSharedTPCClusters(999.); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxFractionSharedTPCClusters(999.); + + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsSharedTPC", SharedTPCClusters, pT, eta, phi, trackType-1); + + fHybridESDtrackCuts->GetMainCuts()->SetMaxFractionSharedTPCClusters(maxSharedTPC); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxFractionSharedTPCClusters(maxSharedTPC_Additional); + + // ################################################################ + // ################################################################ + Bool_t reqTPCRefit = fHybridESDtrackCuts->GetMainCuts()->GetRequireTPCRefit(); + Bool_t reqTPCRefit_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetRequireTPCRefit(); + + fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(1); + fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(1); + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsTPCRefit", 1, pT, eta, phi, trackType-1); + else // track is not accepted as global hybrid with TPC refit requirement + { + fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(0); + fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(0); + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsTPCRefit", 0, pT, eta, phi, trackType-1); } + fHybridESDtrackCuts->GetMainCuts()->SetRequireTPCRefit(reqTPCRefit); + fHybridESDtrackCuts->GetAdditionalCuts()->SetRequireTPCRefit(reqTPCRefit_Additional); - #ifdef DEBUGMODE - AliInfo("Ending GetPythiaTrials."); - #endif - if (pythiaHeader) - return pythiaHeader->Trials(); + // ################################################################ + // ################################################################ + Float_t maxChi2ITS = fHybridESDtrackCuts->GetMainCuts()->GetMaxChi2PerClusterITS(); + Float_t maxChi2ITS_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMaxChi2PerClusterITS(); + fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterITS(999.); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterITS(999.); - AliWarning(Form("In task %s: GetPythiaTrials() failed!", GetName())); - return -1.0; -} + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsChi2ITS", chi2ITS, pT, eta, phi, trackType-1); + fHybridESDtrackCuts->GetMainCuts()->SetMaxChi2PerClusterITS(maxChi2ITS); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMaxChi2PerClusterITS(maxChi2ITS_Additional); + // ################################################################ + // ################################################################ + Float_t minTpcLength = fHybridESDtrackCuts->GetMainCuts()->GetMinLengthActiveVolumeTPC(); // Active length TPC + Float_t minTpcLength_Additional = fHybridESDtrackCuts->GetAdditionalCuts()->GetMinLengthActiveVolumeTPC(); + fHybridESDtrackCuts->GetMainCuts()->SetMinLengthActiveVolumeTPC(0); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinLengthActiveVolumeTPC(0); -//________________________________________________________________________ -inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin() -{ - #ifdef DEBUGMODE - AliInfo("Starting GetPtHardBin."); - #endif - // ########## PT HARD BIN EDGES - const Int_t kPtHardLowerEdges[] = { 0, 5,11,21,36,57, 84,117,152,191,234}; - const Int_t kPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000}; + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsTPCLength", tpcLength, pT, eta, phi, trackType-1); - Int_t tmpPtHardBin = 0; - Double_t tmpPtHard = GetPtHard(); - - for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++) - if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin]) - break; + fHybridESDtrackCuts->GetMainCuts()->SetMinLengthActiveVolumeTPC(minTpcLength); + fHybridESDtrackCuts->GetAdditionalCuts()->SetMinLengthActiveVolumeTPC(minTpcLength_Additional); - #ifdef DEBUGMODE - AliInfo("Ending GetPtHardBin."); - #endif - return tmpPtHardBin; + // ################################################################ + // ################################################################ + + if(!fIsPA) + { + trackType = fHybridESDtrackCuts->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsClustersPtDependence", 1, pT, eta, phi, trackType-1); + + trackType = fHybridESDtrackCuts_noPtDep->AcceptTrack(track); + if (trackType) + FillCutHistogram("hCutsClustersPtDependence", 0, pT, eta, phi, trackType-1); + } + // ################################################################ + // ################################################################ + if((fHybridESDtrackCuts->GetMainCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD) == AliESDtrackCuts::kOff) + || (fHybridESDtrackCuts->GetAdditionalCuts() && (fHybridESDtrackCuts->GetAdditionalCuts()->GetClusterRequirementITS(AliESDtrackCuts::kSPD) == AliESDtrackCuts::kOff))) + { + Bool_t isConstrainedWithITSRefit = static_cast(track->GetConstrainedParam()) && ((track->GetStatus())&AliESDtrack::kITSrefit); + if (trackType) + FillCutHistogram("hCutsTrackConstrained", isConstrainedWithITSRefit, pT, eta, phi, trackType-1); + } + } + + SetCurrentOutputList(0); } + //________________________________________________________________________ Double_t AliAnalysisTaskChargedJetsPA::GetExternalRho() { @@ -391,53 +620,6 @@ Double_t AliAnalysisTaskChargedJetsPA::GetExternalRho() return (rho->GetVal()); } - -//________________________________________________________________________ -inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius) -{ - // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi) - Double_t trackPhi = 0.0; - if (track->Phi() > (TMath::TwoPi() - (radius-phi))) - trackPhi = track->Phi() - TMath::TwoPi(); - else if (track->Phi() < (phi+radius - TMath::TwoPi())) - trackPhi = track->Phi() + TMath::TwoPi(); - else - trackPhi = track->Phi(); - - if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius) - return kTRUE; - - return kFALSE; -} - -//________________________________________________________________________ -inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInJet(AliEmcalJet* jet, Int_t trackIndex) -{ - for (Int_t i = 0; i < jet->GetNumberOfTracks(); ++i) - { - Int_t jetTrack = jet->TrackAt(i); - if (jetTrack == trackIndex) - return kTRUE; - } - return kFALSE; -} - -//________________________________________________________________________ -inline Bool_t AliAnalysisTaskChargedJetsPA::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2) -{ - for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i) - { - Int_t jet1Track = jet1->TrackAt(i); - for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j) - { - Int_t jet2Track = jet2->TrackAt(j); - if (jet1Track == jet2Track) - return kTRUE; - } - } - return kFALSE; -} - //________________________________________________________________________ inline Bool_t AliAnalysisTaskChargedJetsPA::IsEventInAcceptance(AliVEvent* event) { @@ -451,12 +633,24 @@ inline Bool_t AliAnalysisTaskChargedJetsPA::IsEventInAcceptance(AliVEvent* event FillHistogram("hEventAcceptance", 1.5); // number of events after pileup cuts - if(fAnalyzeQA) - FillHistogram("hVertexZ",event->GetPrimaryVertex()->GetZ()); + fPrimaryVertex = event->GetPrimaryVertex(); + + FillHistogram("hVertexZBeforeVertexCut",fPrimaryVertex->GetZ()); - if(fUseVertexCut) + if(fUseDefaultVertexCut) + { if(!fHelperClass->IsVertexSelected2013pA(event)) return kFALSE; + } + else // Failsafe vertex cut + { + if(!fPrimaryVertex || (TMath::Abs(fPrimaryVertex->GetZ()) > 10.0) || (fPrimaryVertex->GetNContributors()<1)) + return kFALSE; + } + + + FillHistogram("hVertexZAfterVertexCut",fPrimaryVertex->GetZ()); + FillHistogram("hEventAcceptance", 2.5); // number of events after vertex cut return kTRUE; @@ -468,13 +662,7 @@ inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* tr FillHistogram("hTrackAcceptance", 0.5); if (track != 0) { - if(fIsKinematics) - { - // TODO: Only working for AOD MC - if((!track->Charge()) || (!(static_cast(track))->IsPhysicalPrimary()) ) - return kFALSE; - } - if (TMath::Abs(track->Eta()) <= fTrackEtaWindow) + if ((track->Eta() < fMaxEta) && (track->Eta() >= fMinEta)) { FillHistogram("hTrackAcceptance", 1.5); if (track->Pt() >= fMinTrackPt) @@ -491,7 +679,7 @@ inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* tr inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet) { if (jet != 0) - if (TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow) + if ((jet->Eta() >= fMinJetEta) && (jet->Eta() < fMaxJetEta)) if (jet->Pt() >= fMinBackgroundJetPt) return kTRUE; @@ -499,36 +687,36 @@ inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcal } //________________________________________________________________________ -inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet) -{ +inline Bool_t AliAnalysisTaskChargedJetsPA::IsSignalJetInAcceptance(AliEmcalJet *jet, Bool_t usePtCut) +{ + Bool_t acceptedWithPtCut = kFALSE; + Bool_t acceptedWithoutPtCut = kFALSE; + FillHistogram("hJetAcceptance", 0.5); if (jet != 0) - if (TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow) + if ((jet->Eta() >= fMinJetEta) && (jet->Eta() < fMaxJetEta)) { FillHistogram("hJetAcceptance", 1.5); - if (jet->Pt() >= fMinJetPt) + if (jet->Pt() >= fMinJetPt) // jet fulfills pt cut { FillHistogram("hJetAcceptance", 2.5); if (jet->Area() >= fMinJetArea) { FillHistogram("hJetAcceptance", 3.5); - return kTRUE; + acceptedWithPtCut = kTRUE; } } + else if(!usePtCut) // jet does not fulfill pt cut + { + if (jet->Area() >= fMinJetArea) + acceptedWithoutPtCut = kTRUE; + } } - return kFALSE; -} - -//________________________________________________________________________ -inline Bool_t AliAnalysisTaskChargedJetsPA::IsDijet(AliEmcalJet *jet1, AliEmcalJet *jet2) -{ - // Output from GetDeltaPhi is < pi in any case - if ((jet1 != 0) && (jet2 != 0)) - if((TMath::Pi() - GetDeltaPhi(jet1->Phi(),jet2->Phi())) < fDijetMaxAngleDeviation) - if ((jet1->Pt() > fMinDijetLeadingPt) && (jet2->Pt() > fMinDijetLeadingPt)) //TODO: Introduce recoil jet? - return kTRUE; - return kFALSE; + if(usePtCut) + return (acceptedWithPtCut); + else + return (acceptedWithPtCut || acceptedWithoutPtCut); } //________________________________________________________________________ @@ -540,77 +728,44 @@ void AliAnalysisTaskChargedJetsPA::ExecOnce() fInitialized = kTRUE; // Check for track array - if (strcmp(fTrackArrayName->Data(), "") != 0) + if (strcmp(fTrackArrayName.Data(), "") != 0) { - fTrackArray = dynamic_cast(InputEvent()->FindListObject(fTrackArrayName->Data())); - fHasTracks = kTRUE; + fTrackArray = dynamic_cast(InputEvent()->FindListObject(fTrackArrayName.Data())); if (!fTrackArray) - { - AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data())); - fHasTracks = kFALSE; - } + AliWarning(Form("%s: Could not retrieve tracks %s!", GetName(), fTrackArrayName.Data())); else { TClass *cl = fTrackArray->GetClass(); if (!cl->GetBaseClass("AliVParticle")) { - AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data())); - fTrackArray = 0; - fHasTracks = kFALSE; + AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName.Data())); + fTrackArray = 0; } } } // Check for jet array - if (strcmp(fJetArrayName->Data(), "") != 0) + if (strcmp(fJetArrayName.Data(), "") != 0) { - fJetArray = dynamic_cast(InputEvent()->FindListObject(fJetArrayName->Data())); - fHasJets = kTRUE; - + fJetArray = dynamic_cast(InputEvent()->FindListObject(fJetArrayName.Data())); if (!fJetArray) - { - AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data())); - fHasJets = kFALSE; - } + AliWarning(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrayName.Data())); else { if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet")) { - AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data())); + AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName.Data())); fJetArray = 0; - fHasJets = kFALSE; } } } // Check for background object - if (strcmp(fBackgroundJetArrayName->Data(), "") != 0) + if (strcmp(fBackgroundJetArrayName.Data(), "") != 0) { - fBackgroundJetArray = dynamic_cast(InputEvent()->FindListObject(fBackgroundJetArrayName->Data())); - fHasBackgroundJets = kTRUE; + fBackgroundJetArray = dynamic_cast(InputEvent()->FindListObject(fBackgroundJetArrayName.Data())); if (!fBackgroundJetArray) - { - AliInfo(Form("%s: Could not retrieve background jets %s! This is OK, if background is not demanded.", GetName(), fBackgroundJetArrayName->Data())); - fHasBackgroundJets = kFALSE; - } - } - - // Look, if initialization is OK - if (!fHasTracks && fAnalyzeBackground) - { - AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName())); - fAnalyzeBackground = kFALSE; - } - if ((!fHasJets && fAnalyzeJets) || (!fHasJets && fAnalyzeBackground)) - { - AliError(Form("%s: Jets NOT successfully casted although demanded! Deactivating jet- and background analysis",GetName())); - fAnalyzeJets = kFALSE; - fAnalyzeBackground = kFALSE; - } - if (!fHasBackgroundJets && fAnalyzeBackground) - { - AliError(Form("%s: Background NOT successfully casted although demanded! Deactivating background analysis",GetName())); - fAnalyzeBackground = kFALSE; + AliInfo(Form("%s: Could not retrieve background jets %s!", GetName(), fBackgroundJetArrayName.Data())); } // Initialize helper class (for vertex selection & pile up correction) @@ -618,6 +773,8 @@ void AliAnalysisTaskChargedJetsPA::ExecOnce() fHelperClass->SetCutOnZVertexSPD(kFALSE); // Histogram init Init(); + // Trackcut initialization + InitializeTrackcuts(); #ifdef DEBUGMODE AliInfo("ExecOnce done."); @@ -626,111 +783,171 @@ void AliAnalysisTaskChargedJetsPA::ExecOnce() } //________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetSignalJets() +void AliAnalysisTaskChargedJetsPA::GetLeadingJets() { // Reset vars fFirstLeadingJet = NULL; fSecondLeadingJet = NULL; + fFirstLeadingKTJet = NULL; + fSecondLeadingKTJet = NULL; + fNumberSignalJets = 0; + fNumberSignalJetsAbove5GeV = 0; + + Int_t jetIDArray[] = {-1, -1}; + Float_t maxJetPts[] = {0, 0}; + jetIDArray[0] = -1; + jetIDArray[1] = -1; + + Int_t jetIDArrayKT[] = {-1, -1}; + Float_t maxJetPtsKT[] = {0, 0}; + jetIDArrayKT[0] = -1; + jetIDArrayKT[1] = -1; - TList tmpJets; + // Find leading signal jets for (Int_t i = 0; i < fJetArray->GetEntries(); i++) { AliEmcalJet* jet = static_cast(fJetArray->At(i)); - if (!jet) + if (!jet) { AliError(Form("%s: Could not receive jet %d", GetName(), i)); continue; } - if (!IsSignalJetInAcceptance(jet)) - continue; - for (Int_t j = 0; j <= tmpJets.GetEntries(); j++) - { - if (j>tmpJets.GetEntries()-1) // When passed last item add the jet at the end - { - tmpJets.Add(jet); - break; - } + if (!IsSignalJetInAcceptance(jet)) continue; - AliEmcalJet* listJet = static_cast(tmpJets.At(j)); - - if(jet->Pt() < listJet->Pt()) // Insert jet before that one in list if pt smaller - { - tmpJets.AddAt(jet, j); - break; - } + if (jet->Pt() > maxJetPts[0]) + { + maxJetPts[1] = maxJetPts[0]; + jetIDArray[1] = jetIDArray[0]; + maxJetPts[0] = jet->Pt(); + jetIDArray[0] = i; } + else if (jet->Pt() > maxJetPts[1]) + { + maxJetPts[1] = jet->Pt(); + jetIDArray[1] = i; + } + fNumberSignalJets++; + if(jet->Pt() >= 5.) + fNumberSignalJetsAbove5GeV++; } - for (Int_t i = 0; i < tmpJets.GetEntries(); i++) + // Find leading background jets + for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++) { - AliEmcalJet* jet = static_cast(tmpJets.At(i)); - fSignalJets[fNumberSignalJets] = jet; - fNumberSignalJets++; + AliEmcalJet* jet = static_cast(fBackgroundJetArray->At(i)); + if (!jet) + { + AliError(Form("%s: Could not receive jet %d", GetName(), i)); + continue; + } + + if (!IsBackgroundJetInAcceptance(jet)) continue; + + if (jet->Pt() > maxJetPtsKT[0]) + { + maxJetPtsKT[1] = maxJetPtsKT[0]; + jetIDArrayKT[1] = jetIDArrayKT[0]; + maxJetPtsKT[0] = jet->Pt(); + jetIDArrayKT[0] = i; + } + else if (jet->Pt() > maxJetPtsKT[1]) + { + maxJetPtsKT[1] = jet->Pt(); + jetIDArrayKT[1] = i; + } } - - if (fNumberSignalJets > 0) - fFirstLeadingJet = static_cast(tmpJets.At(0)); - if (fNumberSignalJets > 1) - fSecondLeadingJet = static_cast(tmpJets.At(1)); + if (jetIDArray[0] > -1) + fFirstLeadingJet = static_cast(fJetArray->At(jetIDArray[0])); + if (jetIDArray[1] > -1) + fSecondLeadingJet = static_cast(fJetArray->At(jetIDArray[1])); + if (jetIDArrayKT[0] > -1) + fFirstLeadingKTJet = static_cast(fBackgroundJetArray->At(jetIDArrayKT[0])); + if (jetIDArrayKT[1] > -1) + fSecondLeadingKTJet = static_cast(fBackgroundJetArray->At(jetIDArrayKT[1])); } //________________________________________________________________________ -Int_t AliAnalysisTaskChargedJetsPA::GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets) +inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t phi, Double_t radius) { -// Writes first two leading jets into already registered array jetIDArray + Double_t tmpConePt = 0.0; - if (!jetArray) + for (Int_t i = 0; i < fTrackArray->GetEntries(); i++) { - AliError("Could not get the jet array to get leading jets from it!"); - return 0; + AliVTrack* tmpTrack = static_cast(fTrackArray->At(i)); + if (IsTrackInAcceptance(tmpTrack)) + if(IsTrackInCone(tmpTrack, eta, phi, radius)) + tmpConePt = tmpConePt + tmpTrack->Pt(); } + return tmpConePt; +} - Float_t maxJetPts[] = {0, 0}; - jetIDArray[0] = -1; - jetIDArray[1] = -1; - - Int_t jetCount = jetArray->GetEntries(); - Int_t jetCountAccepted = 0; +//________________________________________________________________________ +inline Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedConePt(Double_t eta, Double_t phi, Double_t radius, Double_t background) +{ + Double_t tmpConePt = 0.0; - for (Int_t i = 0; i < jetCount; i++) + for (Int_t i = 0; i < fTrackArray->GetEntries(); i++) { - AliEmcalJet* jet = static_cast(jetArray->At(i)); - if (!jet) - { - AliError(Form("%s: Could not receive jet %d", GetName(), i)); - continue; - } + AliVTrack* tmpTrack = static_cast(fTrackArray->At(i)); + if (IsTrackInAcceptance(tmpTrack)) + if(IsTrackInCone(tmpTrack, eta, phi, radius)) + tmpConePt = tmpConePt + tmpTrack->Pt(); + } + Double_t realConeArea = (1.0*(fMaxEta-fMinEta)) * TMath::TwoPi() * MCGetOverlapCircleRectancle(eta, phi, radius, fMinEta, fMaxEta, 0., TMath::TwoPi()); + tmpConePt -= background * realConeArea; // subtract background - if(isSignalJets) - { - if (!IsSignalJetInAcceptance(jet)) continue; - } - else - { - if (!IsBackgroundJetInAcceptance(jet)) continue; - } + return tmpConePt; +} - if (jet->Pt() > maxJetPts[0]) - { - maxJetPts[1] = maxJetPts[0]; - jetIDArray[1] = jetIDArray[0]; - maxJetPts[0] = jet->Pt(); - jetIDArray[0] = i; - } - else if (jet->Pt() > maxJetPts[1]) +//________________________________________________________________________ +inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius) +{ + // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi) + Double_t trackPhi = 0.0; + if (track->Phi() > (TMath::TwoPi() - (radius-phi))) + trackPhi = track->Phi() - TMath::TwoPi(); + else if (track->Phi() < (phi+radius - TMath::TwoPi())) + trackPhi = track->Phi() + TMath::TwoPi(); + else + trackPhi = track->Phi(); + + if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius) + return kTRUE; + + return kFALSE; +} + +//________________________________________________________________________ +inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInJet(AliEmcalJet* jet, Int_t trackIndex) +{ + for (Int_t i = 0; i < jet->GetNumberOfTracks(); ++i) + { + Int_t jetTrack = jet->TrackAt(i); + if (jetTrack == trackIndex) + return kTRUE; + } + return kFALSE; +} + +//________________________________________________________________________ +inline Bool_t AliAnalysisTaskChargedJetsPA::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2) +{ + for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i) + { + Int_t jet1Track = jet1->TrackAt(i); + for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j) { - maxJetPts[1] = jet->Pt(); - jetIDArray[1] = i; + Int_t jet2Track = jet2->TrackAt(j); + if (jet1Track == jet2Track) + return kTRUE; } - jetCountAccepted++; } - return jetCountAccepted; + return kFALSE; } - //________________________________________________________________________ Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background) { @@ -738,14 +955,7 @@ Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Doubl AliInfo("Getting corrected jet spectra."); #endif - if(!jet) - { - AliError("Jet pointer passed to GetCorrectedJetPt() not valid!"); - return -1.0; - } - Double_t correctedPt = -1.0; - // if the passed background is not valid, do not subtract it if(background < 0) background = 0; @@ -774,8 +984,8 @@ Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leading // Define eta range Double_t etaMin, etaMax; - etaMin = -(fTrackEtaWindow-fRandConeRadius); - etaMax = +(fTrackEtaWindow-fRandConeRadius); + etaMin = fMinEta+fRandConeRadius; + etaMax = fMaxEta-fRandConeRadius; // Define random cone Bool_t coneValid = kTRUE; @@ -792,7 +1002,7 @@ Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leading AliEmcalJet* tmpJet = static_cast(fJetArray->At(i)); // if jet is in acceptance and higher, take as new leading if (tmpJet) - if ((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea)) + if ( ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta)) && (tmpJet->Area() >= fMinJetArea)) if((!tmpLeading) || (tmpJet->Pt() > tmpLeading->Pt())) tmpLeading = tmpJet; } @@ -815,7 +1025,6 @@ Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t rho, Double_t leading } } - // Get the cones' pt and calculate delta pt if (coneValid) deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi()); @@ -858,10 +1067,6 @@ void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensityAll(Int_t numberExclude Int_t rhoMeanJetCount = 0; - // Find 2 leading KT jets for the original PbPb approach - Int_t leadingKTJets[] = {-1, -1}; - GetLeadingJets(fBackgroundJetArray, &leadingKTJets[0], kFALSE); - // Exclude UP TO numberExcludeLeadingJets if(numberExcludeLeadingJets==-1) numberExcludeLeadingJets = fNumberSignalJets; @@ -878,19 +1083,6 @@ void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensityAll(Int_t numberExclude continue; } - // Search for overlap with signal jets - Bool_t isOverlapping = kFALSE; - for(Int_t j=0;jArea(); if(backgroundJet->Pt() > 0.150) tmpCoveredArea += backgroundJet->Area(); @@ -898,144 +1090,90 @@ void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensityAll(Int_t numberExclude if (!IsBackgroundJetInAcceptance(backgroundJet)) continue; - Double_t tmpRho = 0.0; - if(backgroundJet->Area()) - tmpRho = backgroundJet->Pt() / backgroundJet->Area(); - - // PbPb approach (take ghosts into account) - if ((i != leadingKTJets[0]) && (i != leadingKTJets[1])) - { - tmpRhoPbPbWithGhosts[rhoPbPbWithGhostsJetCount] = tmpRho; - rhoPbPbWithGhostsJetCount++; - } - - if(backgroundJet->Pt() > 0.150) - { - // CMS approach: don't take ghosts into acount - tmpRhoCMS[rhoCMSJetCount] = tmpRho; - rhoCMSJetCount++; - - // Improved CMS approach: like CMS but excluding signal - if(!isOverlapping) - { - tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho; - rhoImprovedCMSJetCount++; - } - - // PbPb w/o ghosts approach (just neglect ghosts) - if ((i != leadingKTJets[0]) && (i != leadingKTJets[1])) - { - tmpRhoPbPb[rhoPbPbJetCount] = tmpRho; - rhoPbPbJetCount++; - } - } - - // (no overlap with signal jets) - if(!isOverlapping) - { - // Mean approach - tmpRhoMean[rhoMeanJetCount] = tmpRho; - rhoMeanJetCount++; - - // Track like approach approach - tmpPtTrackLike += backgroundJet->Pt(); - tmpAreaTrackLike += backgroundJet->Area(); - } - - } - - if (tmpAreaTrackLike > 0) - rhoTrackLike = tmpPtTrackLike/tmpAreaTrackLike; - if (rhoPbPbJetCount > 0) - rhoPbPb = TMath::Median(rhoPbPbJetCount, tmpRhoPbPb); - if (rhoPbPbWithGhostsJetCount > 0) - rhoPbPbWithGhosts = TMath::Median(rhoPbPbWithGhostsJetCount, tmpRhoPbPbWithGhosts); - if (rhoCMSJetCount > 0) - rhoCMS = TMath::Median(rhoCMSJetCount, tmpRhoCMS) * tmpCoveredArea/tmpSummedArea; - if (rhoImprovedCMSJetCount > 0) - { - rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea; - } - if (rhoMeanJetCount > 0) - rhoMean = TMath::Mean(rhoMeanJetCount, tmpRhoMean); - - #ifdef DEBUGMODE - AliInfo("Got ALL KT background density."); - #endif -} - -//________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoImprovedCMS) -{ - #ifdef DEBUGMODE - AliInfo("Getting KT background density."); - #endif - - static Double_t tmpRhoImprovedCMS[1024]; - Double_t tmpCoveredArea = 0.0; - Double_t tmpSummedArea = 0.0; - - // Setting invalid values - rhoImprovedCMS = 0.0; - - Int_t rhoImprovedCMSJetCount = 0; - - // Exclude UP TO numberExcludeLeadingJets - if(numberExcludeLeadingJets==-1) - numberExcludeLeadingJets = fNumberSignalJets; - if (fNumberSignalJets < numberExcludeLeadingJets) - numberExcludeLeadingJets = fNumberSignalJets; - - for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++) - { - AliEmcalJet* backgroundJet = static_cast(fBackgroundJetArray->At(i)); - - if (!backgroundJet) - { - AliError(Form("%s: Could not receive jet %d", GetName(), i)); - continue; - } - // Search for overlap with signal jets Bool_t isOverlapping = kFALSE; for(Int_t j=0;jPt() < 5.0) + continue; + if(IsJetOverlapping(signalJet, backgroundJet)) { isOverlapping = kTRUE; break; - } - } - - tmpSummedArea += backgroundJet->Area(); - if(backgroundJet->Pt() > 0.150) - tmpCoveredArea += backgroundJet->Area(); + } + } - if (!IsBackgroundJetInAcceptance(backgroundJet)) - continue; + Double_t tmpRho = 0.0; + if(backgroundJet->Area()) + tmpRho = backgroundJet->Pt() / backgroundJet->Area(); - Double_t tmpRho = backgroundJet->Pt() / backgroundJet->Area(); + // PbPb approach (take ghosts into account) + if((backgroundJet != fFirstLeadingKTJet) || (backgroundJet != fSecondLeadingKTJet)) + { + tmpRhoPbPbWithGhosts[rhoPbPbWithGhostsJetCount] = tmpRho; + rhoPbPbWithGhostsJetCount++; + } if(backgroundJet->Pt() > 0.150) - if(!isOverlapping) + { + // CMS approach: don't take ghosts into acount + tmpRhoCMS[rhoCMSJetCount] = tmpRho; + rhoCMSJetCount++; + + // Improved CMS approach: like CMS but excluding signal + if((backgroundJet != fFirstLeadingKTJet) || (backgroundJet != fSecondLeadingKTJet)) { tmpRhoImprovedCMS[rhoImprovedCMSJetCount] = tmpRho; rhoImprovedCMSJetCount++; } + + // PbPb w/o ghosts approach (just neglect ghosts) + if((backgroundJet != fFirstLeadingKTJet) || (backgroundJet != fSecondLeadingKTJet)) + { + tmpRhoPbPb[rhoPbPbJetCount] = tmpRho; + rhoPbPbJetCount++; + } + } + + // (no overlap with signal jets) + if(!isOverlapping) + { + // Mean approach + tmpRhoMean[rhoMeanJetCount] = tmpRho; + rhoMeanJetCount++; + + // Track like approach approach + tmpPtTrackLike += backgroundJet->Pt(); + tmpAreaTrackLike += backgroundJet->Area(); + } + } + if (tmpAreaTrackLike > 0) + rhoTrackLike = tmpPtTrackLike/tmpAreaTrackLike; + if (rhoPbPbJetCount > 0) + rhoPbPb = TMath::Median(rhoPbPbJetCount, tmpRhoPbPb); + if (rhoPbPbWithGhostsJetCount > 0) + rhoPbPbWithGhosts = TMath::Median(rhoPbPbWithGhostsJetCount, tmpRhoPbPbWithGhosts); + if (rhoCMSJetCount > 0) + rhoCMS = TMath::Median(rhoCMSJetCount, tmpRhoCMS) * tmpCoveredArea/tmpSummedArea; if (rhoImprovedCMSJetCount > 0) { rhoImprovedCMS = TMath::Median(rhoImprovedCMSJetCount, tmpRhoImprovedCMS) * tmpCoveredArea/tmpSummedArea; } + if (rhoMeanJetCount > 0) + rhoMean = TMath::Mean(rhoMeanJetCount, tmpRhoMean); + #ifdef DEBUGMODE - AliInfo("Got KT background density."); + AliInfo("Got ALL KT background density."); #endif } - //________________________________________________________________________ void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoNoExclusion, Double_t& rhoConeExclusion02, Double_t& rhoConeExclusion04, Double_t& rhoConeExclusion06, Double_t& rhoConeExclusion08, Double_t& rhoExactExclusion) { @@ -1060,9 +1198,14 @@ void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLea // Exclude UP TO numberExcludeLeadingJets if(numberExcludeLeadingJets==-1) - numberExcludeLeadingJets = fNumberSignalJets; + numberExcludeLeadingJets = fNumberSignalJetsAbove5GeV; if (fNumberSignalJets < numberExcludeLeadingJets) - numberExcludeLeadingJets = fNumberSignalJets; + numberExcludeLeadingJets = fNumberSignalJetsAbove5GeV; + if(numberExcludeLeadingJets>2) + { + AliWarning("Warning: GetTRBackgroundDensity() can only exclude up to 2 leading jets!"); + numberExcludeLeadingJets = 2; + } for (Int_t i = 0; i < fTrackArray->GetEntries(); i++) { @@ -1074,7 +1217,12 @@ void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLea // Check if tracks overlaps with jet for(Int_t j=0;jPt() < 5.0) + continue; // Exact jet exclusion if (IsTrackInJet(signalJet, i)) @@ -1133,27 +1281,36 @@ void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLea // Calculate the correct area where the tracks were taking from - Double_t tmpFullTPCArea = (2.0*fTrackEtaWindow) * TMath::TwoPi(); + Double_t tmpFullTPCArea = (1.0*(fMaxEta-fMinEta)) * TMath::TwoPi(); Double_t tmpAreaCone02 = tmpFullTPCArea; Double_t tmpAreaCone04 = tmpFullTPCArea; Double_t tmpAreaCone06 = tmpFullTPCArea; Double_t tmpAreaCone08 = tmpFullTPCArea; Double_t tmpAreaWithinJets = tmpFullTPCArea; - std::vector tmpEtas(numberExcludeLeadingJets); - std::vector tmpPhis(numberExcludeLeadingJets); + std::vector tmpEtas(fNumberSignalJetsAbove5GeV); + std::vector tmpPhis(fNumberSignalJetsAbove5GeV); + Int_t iSignal = 0; for(Int_t i=0;iEta(); - tmpPhis[i] = tmpJet->Phi(); - tmpAreaWithinJets -= tmpJet->Area(); + AliEmcalJet* signalJet = fFirstLeadingJet; + if(i==1) + signalJet = fSecondLeadingJet; + + if(signalJet->Pt() < 5.0) + continue; + + tmpEtas[iSignal] = signalJet->Eta(); + tmpPhis[iSignal] = signalJet->Phi(); + tmpAreaWithinJets -= signalJet->Area(); + + iSignal++; } - tmpAreaCone02 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.2, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi()); - tmpAreaCone04 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.4, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi()); - tmpAreaCone06 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.6, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi()); - tmpAreaCone08 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(numberExcludeLeadingJets, tmpEtas, tmpPhis, 0.8, -fTrackEtaWindow, +fTrackEtaWindow, 0., TMath::TwoPi()); + tmpAreaCone02 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.2, fMinEta, fMaxEta, 0., TMath::TwoPi()); + tmpAreaCone04 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.4, fMinEta, fMaxEta, 0., TMath::TwoPi()); + tmpAreaCone06 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.6, fMinEta, fMaxEta, 0., TMath::TwoPi()); + tmpAreaCone08 -= tmpFullTPCArea * MCGetOverlapMultipleCirclesRectancle(fNumberSignalJetsAbove5GeV, tmpEtas, tmpPhis, 0.8, fMinEta, fMaxEta, 0., TMath::TwoPi()); rhoConeExclusion02 = summedTracksPtCone02/tmpAreaCone02; rhoConeExclusion04 = summedTracksPtCone04/tmpAreaCone04; @@ -1169,63 +1326,43 @@ void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLea } //________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular) +void AliAnalysisTaskChargedJetsPA::GetPPBackgroundDensity(Double_t& background) { - #ifdef DEBUGMODE - AliInfo("Getting TR background density."); - #endif + // This is the background that was used for the pp 7 TeV ALICE paper + // The background is estimated using the leading jet - // Setting invalid values - Double_t summedTracksPt = 0.0; - rhoMean = 0.0; - area = -1.0; + background = 0; - Double_t tmpRadius = 0.0; - if (doSearchPerpendicular) - tmpRadius = 0.4*TMath::Pi(); // exclude 90 degrees around jets + AliEmcalJet* jet = NULL; + if(fFirstLeadingJet) + jet = fFirstLeadingJet; else - tmpRadius = 0.8; - - numberExcludeLeadingJets = 2; // dijet is excluded here in any case - + return; + Double_t jetMom[3] = { jet->Px(), jet->Py(), jet->Pz() }; + TVector3 jet3mom1(jetMom); + TVector3 jet3mom2(jetMom); - if (!fTrackArray || !fJetArray) - { - AliError("Could not get the track/jet array to calculate track rho!"); - return; - } + jet3mom1.RotateZ(TMath::Pi()); + jet3mom2.RotateZ(-TMath::Pi()); - Int_t trackCount = fTrackArray->GetEntries(); - Int_t trackCountAccepted = 0; - for (Int_t i = 0; i < trackCount; i++) + for (int i = 0; i < fTrackArray->GetEntries(); i++) { - AliVTrack* tmpTrack = static_cast(fTrackArray->At(i)); - if (IsTrackInAcceptance(tmpTrack)) - { - if (IsTrackInCone(tmpTrack, excludeJet1->Eta(), excludeJet1->Phi(), tmpRadius)) - continue; + AliVTrack* track = static_cast(fTrackArray->At(i)); + if (!IsTrackInAcceptance(track)) + continue; - if (numberExcludeLeadingJets > 1) - if (IsTrackInCone(tmpTrack, excludeJet2->Eta(), excludeJet2->Phi(), tmpRadius)) - continue; + Double_t trackMom[3] = { track->Px(), track->Py(), track->Pz() }; + TVector3 track3mom(trackMom); - // Add track pt to array - summedTracksPt = summedTracksPt + tmpTrack->Pt(); - trackCountAccepted++; - } - } + Double_t dR1 = jet3mom1.DeltaR(track3mom); + Double_t dR2 = jet3mom2.DeltaR(track3mom); - if (trackCountAccepted > 0) - { - Double_t tmpArea = 2.0*fTrackEtaWindow*TMath::TwoPi() - 2*(tmpRadius*tmpRadius * TMath::Pi()); //TPC area - excluding jet area - rhoMean = summedTracksPt/tmpArea; - area = tmpArea; + if (dR1 <= fSignalJetRadius || dR2 <= fSignalJetRadius) + background += track3mom.Pt(); } - #ifdef DEBUGMODE - AliInfo("Got TR background density."); - #endif + background /= (2 * TMath::Pi() * fSignalJetRadius * fSignalJetRadius); } //________________________________________________________________________ @@ -1238,11 +1375,6 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) fEventCounter++; - // Check, if analysis should be done in pt hard bins - if(fUsePtHardBin != -1) - if(GetPtHardBin() != fUsePtHardBin) - return; - // This is to take only every Nth event if((fEventCounter+fPartialAnalysisIndex) % fPartialAnalysisNParts != 0) return; @@ -1258,6 +1390,11 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) AliInfo("Calculate()::Init done."); #endif + ////////////////////// NOTE: Create cut histograms + + if(fAnalyzeTrackcuts) + CreateCutHistograms(); + ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background // Get centrality @@ -1267,12 +1404,14 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) Double_t centralityPercentileV0A = 0.0; Double_t centralityPercentileV0C = 0.0; Double_t centralityPercentileV0M = 0.0; + Double_t centralityPercentileZNA = 0.0; if (tmpCentrality != NULL) { centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data()); centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A"); centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C"); centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M"); + centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA"); } if((centralityPercentile < 0.0) || (centralityPercentile > 100.0)) @@ -1281,16 +1420,26 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) if(fSetCentralityToOne) centralityPercentile = 1.0; - // Get jets - if (fAnalyzeBackground || fAnalyzeJets) - GetSignalJets(); + ////////////////////// NOTE: Get event QA histograms - // Get background estimates - Double_t backgroundKTImprovedCMS = -1.0; - Double_t backgroundKTImprovedCMSExternal = -1.0; - Double_t backgroundDijet = -1.0; - Double_t backgroundDijetPerpendicular = -1.0; + FillHistogram("hVertexX",fPrimaryVertex->GetX()); + FillHistogram("hVertexY",fPrimaryVertex->GetY()); + FillHistogram("hVertexXY",fPrimaryVertex->GetX(), fPrimaryVertex->GetY()); + FillHistogram("hVertexR",TMath::Sqrt(fPrimaryVertex->GetX()*fPrimaryVertex->GetX() + fPrimaryVertex->GetY()*fPrimaryVertex->GetY())); + FillHistogram("hCentralityV0M",centralityPercentileV0M); + FillHistogram("hCentralityV0A",centralityPercentileV0A); + FillHistogram("hCentralityV0C",centralityPercentileV0C); + FillHistogram("hCentralityZNA",centralityPercentileZNA); + FillHistogram("hCentrality",centralityPercentile); + + if(!fDoJetAnalysis) + return; + + GetLeadingJets(); + // ##################### Calculate background densities + Double_t backgroundKTImprovedCMS = -1.0; + Double_t backgroundExternal = -1.0; Double_t backgroundKTPbPb = -1.0; Double_t backgroundKTPbPbWithGhosts = -1.0; Double_t backgroundKTCMS = -1.0; @@ -1302,300 +1451,265 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) Double_t backgroundTRCone06 = -1.0; Double_t backgroundTRCone08 = -1.0; Double_t backgroundTRExact = -1.0; + Double_t backgroundPP = -1.0; + Double_t backgroundJetProfile = -1.0; - // Calculate background for different jet exclusions - - if (fAnalyzeBackground) - { - - if(fAnalyzeDeprecatedBackgrounds) - GetKTBackgroundDensityAll (fNumberExcludedJets, backgroundKTPbPb, backgroundKTPbPbWithGhosts, backgroundKTCMS, backgroundKTImprovedCMS, backgroundKTMean, backgroundKTTrackLike); - else - GetKTBackgroundDensity (fNumberExcludedJets, backgroundKTImprovedCMS); - - if(fAnalyzeDeprecatedBackgrounds) - GetTRBackgroundDensity (fNumberExcludedJets, backgroundTRNoExcl, backgroundTRCone02, backgroundTRCone04, backgroundTRCone06, backgroundTRCone08, backgroundTRExact); - - backgroundKTImprovedCMSExternal = GetExternalRho(); - } + // Get background estimates + GetKTBackgroundDensityAll (fNumberExcludedJets, backgroundKTPbPb, backgroundKTPbPbWithGhosts, backgroundKTCMS, backgroundKTImprovedCMS, backgroundKTMean, backgroundKTTrackLike); + GetTRBackgroundDensity (fNumberExcludedJets, backgroundTRNoExcl, backgroundTRCone02, backgroundTRCone04, backgroundTRCone06, backgroundTRCone08, backgroundTRExact); + GetPPBackgroundDensity(backgroundPP); + + backgroundExternal = GetExternalRho(); + if(fNoExternalBackground) + backgroundExternal = 0; + + if(fBackgroundForJetProfile==0) + backgroundJetProfile = backgroundExternal; + else if(fBackgroundForJetProfile==1) + backgroundJetProfile = backgroundKTImprovedCMS; + else if(fBackgroundForJetProfile==2) + backgroundJetProfile = backgroundKTCMS; + else if(fBackgroundForJetProfile==3) + backgroundJetProfile = backgroundPP; + else if(fBackgroundForJetProfile==4) + backgroundJetProfile = backgroundTRCone06; + else if(fBackgroundForJetProfile==5) + backgroundJetProfile = 0; #ifdef DEBUGMODE AliInfo("Calculate()::Centrality&SignalJets&Background-Calculation done."); #endif - if (fAnalyzeQA) - { - FillHistogram("hVertexX",event->GetPrimaryVertex()->GetX()); - FillHistogram("hVertexY",event->GetPrimaryVertex()->GetY()); - FillHistogram("hVertexXY",event->GetPrimaryVertex()->GetX(), event->GetPrimaryVertex()->GetY()); - FillHistogram("hVertexR",TMath::Sqrt(event->GetPrimaryVertex()->GetX()*event->GetPrimaryVertex()->GetX() + event->GetPrimaryVertex()->GetY()*event->GetPrimaryVertex()->GetY())); - FillHistogram("hCentralityV0M",centralityPercentileV0M); - FillHistogram("hCentralityV0A",centralityPercentileV0A); - FillHistogram("hCentralityV0C",centralityPercentileV0C); - FillHistogram("hCentrality",centralityPercentile); - - Int_t trackCountAcc = 0; - Int_t nTracks = fTrackArray->GetEntries(); - for (Int_t i = 0; i < nTracks; i++) - { - AliVTrack* track = static_cast(fTrackArray->At(i)); + // ##################### Fill event QA histograms - if (track != 0) - if (track->Pt() >= fMinTrackPt) - FillHistogram("hTrackPhiEta", track->Phi(),track->Eta(), 1); + Int_t trackCountAcc = 0; + Int_t nTracks = fTrackArray->GetEntries(); + for (Int_t i = 0; i < nTracks; i++) + { + AliVTrack* track = static_cast(fTrackArray->At(i)); - if (IsTrackInAcceptance(track)) + if (track != 0) + if (track->Pt() >= fMinTrackPt) { - FillHistogram("hTrackPt", track->Pt(), centralityPercentile); - if(track->Eta() >= 0) - FillHistogram("hTrackPtPosEta", track->Pt(), centralityPercentile); - else - FillHistogram("hTrackPtNegEta", track->Pt(), centralityPercentile); - - FillHistogram("hTrackEta", track->Eta(), centralityPercentile); - FillHistogram("hTrackPhi", track->Phi()); - - if(static_cast(track)) - { - FillHistogram("hTrackPhiTrackType", track->Phi(), (static_cast(track))->GetTrackType()); - FillHistogram("hTrackPhiLabel", track->Phi(), (static_cast(track))->GetLabel()); - } - for(Int_t j=0;j<20;j++) - if(track->Pt() > j) - FillHistogram("hTrackPhiPtCut", track->Phi(), track->Pt()); + FillHistogram("hTrackPhiEta", track->Phi(),track->Eta(), 1); + FillHistogram("hTrackPtPhiEta", track->Phi(),track->Eta(), track->Pt()); + } - FillHistogram("hTrackCharge", track->Charge()); - trackCountAcc++; + if (IsTrackInAcceptance(track)) + { + FillHistogram("hTrackPt", track->Pt(), centralityPercentile); + + if(track->Eta() >= 0) + FillHistogram("hTrackPtPosEta", track->Pt(), centralityPercentile); + else + FillHistogram("hTrackPtNegEta", track->Pt(), centralityPercentile); + + FillHistogram("hTrackEta", track->Eta(), centralityPercentile); + FillHistogram("hTrackPhi", track->Phi()); + + if(static_cast(track)) + { + FillHistogram("hTrackPhiTrackType", track->Phi(), (static_cast(track))->GetTrackType()); + FillHistogram("hTrackPtTrackType", track->Pt(), (static_cast(track))->GetTrackType()); } - } - FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile); + for(Int_t j=0;j<20;j++) + if(track->Pt() > j) + FillHistogram("hTrackPhiPtCut", track->Phi(), track->Pt()); + + FillHistogram("hTrackCharge", track->Charge()); + trackCountAcc++; + } } + FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile); + #ifdef DEBUGMODE AliInfo("Calculate()::QA done."); #endif - ////////////////////// NOTE: Jet analysis and calculations + // ##################### Fill jet histograms - if (fAnalyzeJets) + FillHistogram("hJetCountAll", fJetArray->GetEntries()); + FillHistogram("hJetCountAccepted", fNumberSignalJets); + FillHistogram("hJetCount", fJetArray->GetEntries(), fNumberSignalJets); + if (fFirstLeadingJet) { - for (Int_t i = 0; iGetEntries(); i++) - { - AliEmcalJet* tmpJet = static_cast(fJetArray->At(i)); - if (!tmpJet) - continue; + FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt()); + FillHistogram("hCorrectedLeadingJetPt", GetCorrectedJetPt(fFirstLeadingJet,backgroundExternal)); + } + if (fSecondLeadingJet) + { + FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt()); + FillHistogram("hCorrectedSecondLeadingJetPt", GetCorrectedJetPt(fSecondLeadingJet,backgroundExternal)); + } - FillHistogram("hRawJetPt", tmpJet->Pt()); + for (Int_t i = 0; iGetEntries(); i++) + { + AliEmcalJet* tmpJet = static_cast(fJetArray->At(i)); + if (!tmpJet) + continue; + + // ### JETS BEFORE ANY CUTS + if (tmpJet->Area() >= fMinJetArea) + FillHistogram("hRawJetPhiEta", tmpJet->Phi(), tmpJet->Eta()); + if ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta)) + FillHistogram("hRawJetArea", tmpJet->Area()); + + FillHistogram("hJetPtCutStages", tmpJet->Pt(), 0.5); + if ((tmpJet->Eta() >= fMinJetEta) && (tmpJet->Eta() < fMaxJetEta)) + { + FillHistogram("hJetPtCutStages", tmpJet->Pt(), 1.5); if (tmpJet->Pt() >= fMinJetPt) { - // ### RAW JET ANALYSIS + FillHistogram("hJetPtCutStages", tmpJet->Pt(), 2.5); if (tmpJet->Area() >= fMinJetArea) - FillHistogram("hRawJetPhiEta", tmpJet->Phi(), tmpJet->Eta()); - if (TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) - FillHistogram("hRawJetArea", tmpJet->Area()); - } - - if(IsSignalJetInAcceptance(tmpJet)) - { - // ### SIGNAL JET ANALYSIS - // Jet spectra - FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile); - FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile); - FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, backgroundKTImprovedCMS); - if(centralityPercentile<=20.0) - FillHistogram("hJetPtSubtractedRhoKTImprovedCMS020", tmpJet->Pt(), backgroundKTImprovedCMS); - - if(fAnalyzeDeprecatedBackgrounds) { - FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRCone06), centralityPercentile); - FillHistogram("hJetPtBgrdSubtractedKTPbPb", GetCorrectedJetPt(tmpJet, backgroundKTPbPb), centralityPercentile); - FillHistogram("hJetPtBgrdSubtractedKTPbPbWithGhosts", GetCorrectedJetPt(tmpJet, backgroundKTPbPbWithGhosts), centralityPercentile); - FillHistogram("hJetPtBgrdSubtractedKTCMS", GetCorrectedJetPt(tmpJet, backgroundKTCMS), centralityPercentile); - FillHistogram("hJetPtBgrdSubtractedKTMean", GetCorrectedJetPt(tmpJet, backgroundKTMean), centralityPercentile); - FillHistogram("hJetPtBgrdSubtractedKTTrackLike", GetCorrectedJetPt(tmpJet, backgroundKTTrackLike), centralityPercentile); + FillHistogram("hJetPtCutStages", tmpJet->Pt(), 3.5); } + } + } + // ### JETS AFTER CUTS + if(IsSignalJetInAcceptance(tmpJet)) + { + // Background corrected jet spectra + FillHistogram("hJetPtBgrdSubtractedExternal", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTImprovedCMS", GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedPP", GetCorrectedJetPt(tmpJet, backgroundPP), centralityPercentile); + if(tmpJet->Phi() >= TMath::Pi()) + FillHistogram("hJetPtBgrdSubtractedExternal_Phi2", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile); + else + FillHistogram("hJetPtBgrdSubtractedExternal_Phi1", GetCorrectedJetPt(tmpJet, backgroundExternal), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRCone06), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTPbPb", GetCorrectedJetPt(tmpJet, backgroundKTPbPb), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTPbPbWithGhosts", GetCorrectedJetPt(tmpJet, backgroundKTPbPbWithGhosts), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTCMS", GetCorrectedJetPt(tmpJet, backgroundKTCMS), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTMean", GetCorrectedJetPt(tmpJet, backgroundKTMean), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTTrackLike", GetCorrectedJetPt(tmpJet, backgroundKTTrackLike), centralityPercentile); + + FillHistogram("hJetPtSubtractedRhoExternal", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundExternal)); + FillHistogram("hJetPtSubtractedRhoKTImprovedCMS", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundKTImprovedCMS)); + FillHistogram("hJetPtSubtractedRhoPP", tmpJet->Pt(), centralityPercentile, tmpJet->Pt() - GetCorrectedJetPt(tmpJet, backgroundPP)); + FillHistogram("hDeltaPtExternalBgrdVsPt", GetDeltaPt(backgroundExternal), GetCorrectedJetPt(tmpJet, backgroundExternal)); + + FillHistogram("hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks()); + + for(Int_t j=0; jGetNumberOfTracks(); j++) + FillHistogram("hJetConstituentPtVsJetPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), tmpJet->Pt()); + + if(tmpJet->Pt() >= 5.0) + { + Double_t lowestTrackPt = 1e99; + Double_t highestTrackPt = 0.0; for(Int_t j=0; jGetNumberOfTracks(); j++) - FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile); - - if(fAnalyzeQA) { - FillHistogram("hJetArea", tmpJet->Area()); - FillHistogram("hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks()); - FillHistogram("hJetPhiEta", tmpJet->Phi(),tmpJet->Eta()); - FillHistogram("hJetEta", tmpJet->Eta(), centralityPercentile); + FillHistogram("hJetConstituentPt", tmpJet->TrackAt(j, fTrackArray)->Pt(), centralityPercentile); + // Find the lowest pT of a track in the jet + if (tmpJet->TrackAt(j, fTrackArray)->Pt() < lowestTrackPt) + lowestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt(); + if (tmpJet->TrackAt(j, fTrackArray)->Pt() > highestTrackPt) + highestTrackPt = tmpJet->TrackAt(j, fTrackArray)->Pt(); } + FillHistogram("hJetArea", tmpJet->Area(), tmpJet->Pt()); // Signal jet vs. signal jet - "Combinatorial" - for (Int_t j = i+1; jPhi(), fSignalJets[j]->Phi())); - } - } + for (Int_t j = 0; jGetEntries(); j++) + { + AliEmcalJet* tmpJet2 = static_cast(fJetArray->At(j)); + if (!tmpJet2) + continue; + if(tmpJet2->Pt() >= 5.0) + FillHistogram("hJetDeltaPhi", GetDeltaPhi(tmpJet->Phi(), tmpJet2->Phi())); + } - // ### DIJETS - if(fNumberSignalJets >= 2) - { - FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi())); + FillHistogram("hJetPhiEta", tmpJet->Phi(),tmpJet->Eta()); + FillHistogram("hJetPtPhiEta", tmpJet->Phi(),tmpJet->Eta(),tmpJet->Pt()); + FillHistogram("hJetEta", tmpJet->Eta(), centralityPercentile); - if (IsDijet(fFirstLeadingJet, fSecondLeadingJet)) - { - FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt()); - FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt()); - - FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt()); - FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt()); - Double_t dummyArea = 0; - GetTRBackgroundDensity (2, backgroundDijet, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE); - GetTRBackgroundDensity (2, backgroundDijetPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE); + if(lowestTrackPt>=2.0) + FillHistogram("hJetEta2GeVTracks", tmpJet->Eta(), centralityPercentile); + if(lowestTrackPt>=4.0) + FillHistogram("hJetEta4GeVTracks", tmpJet->Eta(), centralityPercentile); } } + } // end of jet loop - // ### SOME JET PLOTS - FillHistogram("hJetCountAll", fJetArray->GetEntries()); - FillHistogram("hJetCountAccepted", fNumberSignalJets); - FillHistogram("hJetCount", fJetArray->GetEntries(), fNumberSignalJets); - if (fFirstLeadingJet) - { - FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt()); - FillHistogram("hCorrectedLeadingJetPt", GetCorrectedJetPt(fFirstLeadingJet,backgroundKTImprovedCMS)); - } - if (fSecondLeadingJet) - { - FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt()); - FillHistogram("hCorrectedSecondLeadingJetPt", GetCorrectedJetPt(fSecondLeadingJet,backgroundKTImprovedCMS)); - } - - } //endif AnalyzeJets + if(fAnalyzeJetProfile) + CreateJetProfilePlots(backgroundJetProfile); #ifdef DEBUGMODE AliInfo("Calculate()::Jets done."); #endif - ////////////////////// NOTE: Background analysis - - if (fAnalyzeBackground) - { - // Calculate background in centrality classes - FillHistogram("hKTBackgroundImprovedCMS", backgroundKTImprovedCMS, centralityPercentile); - - FillHistogram("hKTBackgroundImprovedCMSExternal", backgroundKTImprovedCMSExternal, centralityPercentile); - - FillHistogram("hKTMeanBackgroundImprovedCMS", centralityPercentile, backgroundKTImprovedCMS); - - // In case of dijets -> look at the background - if (backgroundDijet >= 0) - FillHistogram("hDijetBackground", backgroundDijet, centralityPercentile); - if (backgroundDijetPerpendicular >= 0) - FillHistogram("hDijetBackgroundPerpendicular", backgroundDijetPerpendicular, centralityPercentile); - - if(fAnalyzeDeprecatedBackgrounds) - { - FillHistogram("hKTBackgroundPbPb", backgroundKTPbPb, centralityPercentile); - FillHistogram("hKTBackgroundPbPbWithGhosts", backgroundKTPbPbWithGhosts, centralityPercentile); - FillHistogram("hKTBackgroundCMS", backgroundKTCMS, centralityPercentile); - FillHistogram("hKTBackgroundMean", backgroundKTMean, centralityPercentile); - FillHistogram("hKTBackgroundTrackLike", backgroundKTTrackLike, centralityPercentile); - - FillHistogram("hTRBackgroundNoExcl", backgroundTRNoExcl, centralityPercentile); - FillHistogram("hTRBackgroundCone02", backgroundTRCone02, centralityPercentile); - FillHistogram("hTRBackgroundCone04", backgroundTRCone04, centralityPercentile); - FillHistogram("hTRBackgroundCone06", backgroundTRCone06, centralityPercentile); - FillHistogram("hTRBackgroundCone08", backgroundTRCone08, centralityPercentile); - FillHistogram("hTRBackgroundExact", backgroundTRExact, centralityPercentile); - - // Calculate background profiles in terms of centrality - FillHistogram("hKTMeanBackgroundPbPb", centralityPercentile, backgroundKTPbPb); - FillHistogram("hKTMeanBackgroundPbPbWithGhosts", centralityPercentile, backgroundKTPbPbWithGhosts); - FillHistogram("hKTMeanBackgroundCMS", centralityPercentile, backgroundKTCMS); - FillHistogram("hKTMeanBackgroundMean", centralityPercentile, backgroundKTMean); - FillHistogram("hKTMeanBackgroundTPC", centralityPercentile, backgroundKTTrackLike); - FillHistogram("hTRMeanBackground", centralityPercentile, backgroundTRCone06); - } - - - // Calculate the delta pt - Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0); - Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS); - - Double_t tmpDeltaPtKTImprovedCMSPartialExclusion = 0.0; - if(fNcoll) - tmpDeltaPtKTImprovedCMSPartialExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0/fNcoll); - else - tmpDeltaPtKTImprovedCMSPartialExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0); - - Double_t tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = 0.0; - if(fNumberSignalJets) - tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = GetDeltaPt(backgroundKTImprovedCMS, 1.0/fNumberSignalJets); - else - tmpDeltaPtKTImprovedCMSPartialExclusion_Signal = GetDeltaPt(backgroundKTImprovedCMS, 1.0); - - Double_t tmpDeltaPtKTImprovedCMSFullExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0); - - Double_t tmpDeltaPtKTPbPb = 0; - Double_t tmpDeltaPtKTPbPbWithGhosts = 0; - Double_t tmpDeltaPtKTCMS = 0; - Double_t tmpDeltaPtKTMean = 0; - Double_t tmpDeltaPtKTTrackLike = 0; - Double_t tmpDeltaPtTR = 0; - - if(fAnalyzeDeprecatedBackgrounds) - { - tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb); - tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts); - tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS); - tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean); - tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike); - tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06); - } - - // If valid, fill the delta pt histograms - - if(tmpDeltaPtKTImprovedCMS > -10000.0) - FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile); - if(tmpDeltaPtKTImprovedCMSPartialExclusion > -10000.0) - FillHistogram("hDeltaPtKTImprovedCMSPartialExclusion", tmpDeltaPtKTImprovedCMSPartialExclusion, centralityPercentile); - if(tmpDeltaPtKTImprovedCMSPartialExclusion_Signal > -10000.0) - FillHistogram("hDeltaPtKTImprovedCMSPartialExclusion_Signal", tmpDeltaPtKTImprovedCMSPartialExclusion_Signal, centralityPercentile); - if(tmpDeltaPtKTImprovedCMSFullExclusion > -10000.0) - FillHistogram("hDeltaPtKTImprovedCMSFullExclusion", tmpDeltaPtKTImprovedCMSFullExclusion, centralityPercentile); - - if(tmpDeltaPtNoBackground > 0.000001) - FillHistogram("hDeltaPtNoBackgroundNoEmptyCones", tmpDeltaPtNoBackground, centralityPercentile); - else if(tmpDeltaPtNoBackground > -10000.0) - FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile); + // ##################### Fill background plots + + FillHistogram("hKTBackgroundExternal", backgroundExternal, centralityPercentile); + if(fFirstLeadingJet && (fFirstLeadingJet->Pt()>=20.)) + FillHistogram("hKTBackgroundExternal20GeV", backgroundExternal, centralityPercentile); + + FillHistogram("hKTBackgroundImprovedCMS", backgroundKTImprovedCMS, centralityPercentile); + FillHistogram("hPPBackground", backgroundPP, centralityPercentile); + FillHistogram("hKTBackgroundPbPb", backgroundKTPbPb, centralityPercentile); + FillHistogram("hKTBackgroundPbPbWithGhosts", backgroundKTPbPbWithGhosts, centralityPercentile); + FillHistogram("hKTBackgroundCMS", backgroundKTCMS, centralityPercentile); + FillHistogram("hKTBackgroundMean", backgroundKTMean, centralityPercentile); + FillHistogram("hKTBackgroundTrackLike", backgroundKTTrackLike, centralityPercentile); + FillHistogram("hTRBackgroundNoExcl", backgroundTRNoExcl, centralityPercentile); + FillHistogram("hTRBackgroundCone02", backgroundTRCone02, centralityPercentile); + FillHistogram("hTRBackgroundCone04", backgroundTRCone04, centralityPercentile); + FillHistogram("hTRBackgroundCone06", backgroundTRCone06, centralityPercentile); + FillHistogram("hTRBackgroundCone08", backgroundTRCone08, centralityPercentile); + FillHistogram("hTRBackgroundExact", backgroundTRExact, centralityPercentile); + + // Calculate the delta pt + Double_t tmpDeltaPtNoBackground = GetDeltaPt(0.0); + Double_t tmpDeltaPtExternalBgrd = GetDeltaPt(backgroundExternal); + Double_t tmpDeltaPtKTImprovedCMS = GetDeltaPt(backgroundKTImprovedCMS); + Double_t tmpDeltaPtKTImprovedCMSFullExclusion = GetDeltaPt(backgroundKTImprovedCMS, 1.0); + + Double_t tmpDeltaPtKTPbPb = 0; + Double_t tmpDeltaPtKTPbPbWithGhosts = 0; + Double_t tmpDeltaPtKTCMS = 0; + Double_t tmpDeltaPtKTMean = 0; + Double_t tmpDeltaPtKTTrackLike = 0; + Double_t tmpDeltaPtTR = 0; + + tmpDeltaPtKTPbPb = GetDeltaPt(backgroundKTPbPb); + tmpDeltaPtKTPbPbWithGhosts = GetDeltaPt(backgroundKTPbPbWithGhosts); + tmpDeltaPtKTCMS = GetDeltaPt(backgroundKTCMS); + tmpDeltaPtKTMean = GetDeltaPt(backgroundKTMean); + tmpDeltaPtKTTrackLike = GetDeltaPt(backgroundKTTrackLike); + tmpDeltaPtTR = GetDeltaPt(backgroundTRCone06); + + + // If valid, fill the delta pt histograms + + if(tmpDeltaPtExternalBgrd > -10000.0) + FillHistogram("hDeltaPtExternalBgrd", tmpDeltaPtExternalBgrd, centralityPercentile); + if(tmpDeltaPtKTImprovedCMS > -10000.0) + FillHistogram("hDeltaPtKTImprovedCMS", tmpDeltaPtKTImprovedCMS, centralityPercentile); + if(tmpDeltaPtKTImprovedCMSFullExclusion > -10000.0) + FillHistogram("hDeltaPtKTImprovedCMSFullExclusion", tmpDeltaPtKTImprovedCMSFullExclusion, centralityPercentile); + + if(tmpDeltaPtNoBackground > -10000.0) + FillHistogram("hDeltaPtNoBackground", tmpDeltaPtNoBackground, centralityPercentile); + + if(tmpDeltaPtKTPbPb > -10000.0) + FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile); + if(tmpDeltaPtKTPbPbWithGhosts > -10000.0) + FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile); + if(tmpDeltaPtKTCMS > -10000.0) + FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile); + if(tmpDeltaPtKTMean > -10000.0) + FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile); + if(tmpDeltaPtKTTrackLike > -10000.0) + FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile); + if(tmpDeltaPtTR > -10000.0) + FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile); - if(fAnalyzeDeprecatedBackgrounds) - { - if(tmpDeltaPtKTPbPb > -10000.0) - FillHistogram("hDeltaPtKTPbPb", tmpDeltaPtKTPbPb, centralityPercentile); - if(tmpDeltaPtKTPbPbWithGhosts > -10000.0) - FillHistogram("hDeltaPtKTPbPbWithGhosts", tmpDeltaPtKTPbPbWithGhosts, centralityPercentile); - if(tmpDeltaPtKTCMS > -10000.0) - FillHistogram("hDeltaPtKTCMS", tmpDeltaPtKTCMS, centralityPercentile); - if(tmpDeltaPtKTMean > -10000.0) - FillHistogram("hDeltaPtKTMean", tmpDeltaPtKTMean, centralityPercentile); - if(tmpDeltaPtKTTrackLike > -10000.0) - FillHistogram("hDeltaPtKTTrackLike", tmpDeltaPtKTTrackLike, centralityPercentile); - - if(tmpDeltaPtTR > -10000.0) - FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile); - } - } - #ifdef DEBUGMODE AliInfo("Calculate()::Background done."); #endif - ////////////////////// NOTE: Pythia histograms - if(fAnalyzePythia) - { - FillHistogram("hPythiaPtHard", GetPtHard()); - FillHistogram("hPythiaNTrials", GetPtHardBin()+0.1, GetPythiaTrials()); - FillHistogram("hPythiaXSection", GetPtHardBin()+0.1, fCrossSection); - - #ifdef DEBUGMODE - AliInfo("Calculate()::Pythia done."); - #endif - } #ifdef DEBUGMODE AliInfo("Calculate() done."); #endif @@ -1604,78 +1718,134 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) //________________________________________________________________________ Bool_t AliAnalysisTaskChargedJetsPA::UserNotify() { - // Implemented Notify() to read the cross sections - // and number of trials from pyxsec.root - // - #ifdef DEBUGMODE - AliInfo("UserNotify started."); - #endif + return kTRUE; +} - if(fAnalyzePythia) +//________________________________________________________________________ +void AliAnalysisTaskChargedJetsPA::CreateJetProfilePlots(Double_t bgrd) +{ + for (Int_t i = 0; iGetEntries(); i++) { - TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); - TFile *currFile = tree->GetCurrentFile(); - - TString file(currFile->GetName()); + AliEmcalJet* tmpJet = static_cast(fJetArray->At(i)); + if (!tmpJet) + continue; + if(!IsSignalJetInAcceptance(tmpJet)) + continue; - 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()),""); - } - - 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; + SetCurrentOutputList(1); + // Jet profile analysis + if(TMath::Abs(tmpJet->Eta()) <= 0.3) + { + if(tmpJet->Pt()>=70.0) + { + FillHistogram("hJetProfile70GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile70GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); } - 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(key->ReadObj()); - if(!list){ - fxsec->Close(); - return kFALSE; - } - fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1); - fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1); - fxsec->Close(); + else if(GetCorrectedJetPt(tmpJet, bgrd)>=60.0) + { + FillHistogram("hJetProfile60GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile60GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); } - } // no tree pyxsec.root - else { - TTree *xtree = (TTree*)fxsec->Get("Xsection"); - if(!xtree){ - fxsec->Close(); - return kFALSE; + else if(GetCorrectedJetPt(tmpJet, bgrd)>=50.0) + { + FillHistogram("hJetProfile50GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile50GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + } + else if(GetCorrectedJetPt(tmpJet, bgrd)>=40.0) + { + FillHistogram("hJetProfile40GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile40GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + } + else if(GetCorrectedJetPt(tmpJet, bgrd)>=30.0) + { + FillHistogram("hJetProfile30GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile30GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + } + else if(GetCorrectedJetPt(tmpJet, bgrd)>=20.0) + { + FillHistogram("hJetProfile20GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile20GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + } + else if(GetCorrectedJetPt(tmpJet, bgrd)>=10.0) + { + FillHistogram("hJetProfile10GeV", 0.05-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.05, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.10-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.10, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.15-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.15, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.20-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.20, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.25-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.25, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.30-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.30, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.35-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.35, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.40-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.40, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.45-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.45, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.50-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.50, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.55-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.55, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); + FillHistogram("hJetProfile10GeV", 0.60-0.05/2, (GetCorrectedConePt(tmpJet->Eta(), tmpJet->Phi(), 0.60, bgrd))/GetCorrectedJetPt(tmpJet, bgrd)); } - UInt_t ntrials = 0; - Double_t xsection = 0; - xtree->SetBranchAddress("xsection",&xsection); - xtree->SetBranchAddress("ntrials",&ntrials); - xtree->GetEntry(0); - fTrials = ntrials; - fCrossSection = xsection; - fxsec->Close(); } + SetCurrentOutputList(0); } - #ifdef DEBUGMODE - AliInfo("UserNotify ended."); - #endif - return kTRUE; } - //________________________________________________________________________ inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg) {return 2.*atan(exp(-arg));} @@ -1752,7 +1922,7 @@ Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapMultipleCirclesRectancle(Int_ //________________________________________________________________________ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x) { - TH1* tmpHist = static_cast(fOutputList->FindObject(GetHistoName(key))); + TH1* tmpHist = static_cast(fCurrentOutputList->FindObject(GetHistoName(key))); if(!tmpHist) { AliError(Form("Cannot find histogram <%s> ",key)) ; @@ -1765,7 +1935,7 @@ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double //________________________________________________________________________ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y) { - TH1* tmpHist = static_cast(fOutputList->FindObject(GetHistoName(key))); + TH1* tmpHist = static_cast(fCurrentOutputList->FindObject(GetHistoName(key))); if(!tmpHist) { AliError(Form("Cannot find histogram <%s> ",key)); @@ -1781,7 +1951,7 @@ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double //________________________________________________________________________ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x, Double_t y, Double_t add) { - TH2* tmpHist = static_cast(fOutputList->FindObject(GetHistoName(key))); + TH2* tmpHist = static_cast(fCurrentOutputList->FindObject(GetHistoName(key))); if(!tmpHist) { AliError(Form("Cannot find histogram <%s> ",key)); @@ -1790,6 +1960,21 @@ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double tmpHist->Fill(x,y,add); } + +//________________________________________________________________________ +inline void AliAnalysisTaskChargedJetsPA::FillCutHistogram(const char * key, Double_t cut, Double_t pT, Double_t eta, Double_t phi, Int_t isAdditionalTrack) +{ + THnF* tmpHist = static_cast(fCurrentOutputList->FindObject(GetHistoName(key))); + if(!tmpHist) + { + AliError(Form("Cannot find histogram <%s> ",key)); + return; + } + + Double_t tmpVec[5] = {cut, pT, eta, phi, static_cast(isAdditionalTrack)}; + tmpHist->Fill(tmpVec); +} + //________________________________________________________________________ template T* AliAnalysisTaskChargedJetsPA::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle) { @@ -1801,9 +1986,8 @@ template T* AliAnalysisTaskChargedJetsPA::AddHistogram1D(const char* n tmpHist->SetMarkerStyle(kFullCircle); tmpHist->Sumw2(); - fHistList->Add(tmpHist); - fHistCount++; - + fCurrentOutputList->Add(tmpHist); + return tmpHist; } @@ -1818,23 +2002,81 @@ template T* AliAnalysisTaskChargedJetsPA::AddHistogram2D(const char* n tmpHist->SetMarkerStyle(kFullCircle); tmpHist->Sumw2(); - fHistList->Add(tmpHist); - fHistCount++; + fCurrentOutputList->Add(tmpHist); return tmpHist; } +//________________________________________________________________________ +THnF* AliAnalysisTaskChargedJetsPA::AddCutHistogram(const char* name, const char* title, const char* cutName, Int_t nBins, Double_t xMin, Double_t xMax) +{ + // Cut, pT, eta, phi, type + Int_t bins [5] = { nBins, 100, 20, 18, 2}; + Double_t minEdges[5] = { xMin, 0.1, -1, 0, -0.5}; + Double_t maxEdges[5] = { xMax, 40, +1, 2*TMath::Pi(), 1.5}; + + TString axisName[5] = {cutName,"#it{p}_{T}","#eta","#phi","Track type"}; + TString axisTitle[5] = {cutName,"#it{p}_{T}","#eta","#phi","Track type"}; + + THnF * histo = new THnF(name, title, 5, bins, minEdges, maxEdges); + BinLogAxis(histo, 1); + + for (Int_t iaxis=0; iaxis<5;iaxis++){ + histo->GetAxis(iaxis)->SetName(axisName[iaxis]); + histo->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]); + } + + fCurrentOutputList->Add(histo); + return histo; +} + +//________________________________________________________________________ +void AliAnalysisTaskChargedJetsPA::BinLogAxis(const THn *h, Int_t axisNumber) +{ + // Method for the correct logarithmic binning of histograms + TAxis *axis = h->GetAxis(axisNumber); + int bins = axis->GetNbins(); + + Double_t from = axis->GetXmin(); + Double_t to = axis->GetXmax(); + Double_t *newBins = new Double_t[bins + 1]; + + newBins[0] = from; + Double_t factor = pow(to/from, 1./bins); + + for (int i = 1; i <= bins; i++) { + newBins[i] = factor * newBins[i-1]; + } + axis->Set(bins, newBins); + delete [] newBins; +} + //________________________________________________________________________ void AliAnalysisTaskChargedJetsPA::Terminate(Option_t *) { - PostData(1, fOutputList); +/* + PostData(1, fOutputLists[0]); + fOutputLists[0] = dynamic_cast (GetOutputData(1)); // >1 refers to output slots - // Mandatory - fOutputList = dynamic_cast (GetOutputData(1)); // '1' refers to the output slot - if (!fOutputList) { - printf("ERROR: Output list not available\n"); - return; + if(fAnalyzeJetProfile) + { + PostData(2, fOutputLists[1]); + fOutputLists[1] = dynamic_cast (GetOutputData(2)); // >1 refers to output slots + } + if(fAnalyzeTrackcuts) + { + if(fAnalyzeJetProfile) + { + PostData(3, fOutputLists[2]); + fOutputLists[2] = dynamic_cast (GetOutputData(3)); // >1 refers to output slots + } + else + { + PostData(2, fOutputLists[1]); + fOutputLists[1] = dynamic_cast (GetOutputData(2)); // >1 refers to output slots + } } +*/ } //________________________________________________________________________ @@ -1842,9 +2084,14 @@ AliAnalysisTaskChargedJetsPA::~AliAnalysisTaskChargedJetsPA() { // Destructor. Clean-up the output list, but not the histograms that are put inside // (the list is owner and will clean-up these histograms). Protect in PROOF case. - if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { - delete fOutputList; - } +/* + delete fHybridESDtrackCuts; + delete fHybridESDtrackCuts_noPtDep; + + for(Int_t i=0; i(fOutputLists.size()); i++) + if (fOutputLists[i] && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) + delete fOutputLists[i]; +*/ } //________________________________________________________________________ @@ -1857,15 +2104,19 @@ void AliAnalysisTaskChargedJetsPA::UserCreateOutputObjects() fRandom = new TRandom3(0); - fOutputList = new TList(); - fOutputList->SetOwner(); // otherwise it produces leaks in merging - - // NOTE: Pythia histograms - AddHistogram1D("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+1, 0, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}"); - AddHistogram1D("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+1, 0, fNumPtHardBins+1, "p_{T} hard bin", "Trials"); + Int_t tmpListCount = 1; + if(fAnalyzeJetProfile) + tmpListCount++; + if(fAnalyzeTrackcuts) + tmpListCount++; - - PostData(1, fOutputList); + fOutputLists.resize(tmpListCount); + for(Int_t i=0; iSetOwner(); // otherwise it produces leaks in merging + PostData(i+1, fOutputLists[i]); + } } //________________________________________________________________________ @@ -1885,6 +2136,16 @@ void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *) ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos Calculate(InputEvent()); - - PostData(1, fOutputList); + + PostData(1, fOutputLists[0]); + if(fAnalyzeJetProfile) + PostData(2, fOutputLists[1]); + if(fAnalyzeTrackcuts) + { + if(fAnalyzeJetProfile) + PostData(3, fOutputLists[2]); + else + PostData(2, fOutputLists[1]); + } + }