From: loizides Date: Wed, 20 Feb 2013 22:51:28 +0000 (+0000) Subject: update from Ruediger X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=12eb4ac10fee0c348f3f5ab0f02897ee1637419a;p=u%2Fmrichter%2FAliRoot.git update from Ruediger --- diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx index 18d2e8cc0d4..441f6c6df7c 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx @@ -1,3 +1,43 @@ +#ifndef ALIANALYSISTASKSE_H +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliAnalysisTask.h" +#include "AliCentrality.h" +#include "AliStack.h" +#include "AliESDEvent.h" +#include "AliESDInputHandler.h" +#include "AliAODEvent.h" +#include "AliAODHandler.h" +#include "AliAnalysisManager.h" +#include "AliAnalysisTaskSE.h" +#endif + +#include +#include "AliGenPythiaEventHeader.h" +#include "AliMCEvent.h" +#include "AliLog.h" +#include +#include "AliVEventHandler.h" +#include "AliVParticle.h" +#include "AliAnalysisUtils.h" + + #include "AliAnalysisTaskChargedJetsPA.h" //TODO: Not accessing the particles when using MC @@ -10,141 +50,90 @@ void AliAnalysisTaskChargedJetsPA::Init() #ifdef DEBUGMODE AliInfo("Creating histograms."); #endif - // NOTE: Track & Cluster & QA histograms - if (fAnalyzeQA) - { - - AddHistogram1D("hNumberEvents", "Number of events (0 = before, 1 = after vertex cuts)", "", 2, 0, 2, "#Delta z(cm)","N^{Events}/cut"); - - AddHistogram1D("hAppliedEtaCorrectionFactor", "Applied #eta correction factor for the k_{T} background", "", 500, 0.5, 1.5, "Correction factor","dN^{Jets}/df"); - AddHistogram1D("hAppliedEtaCorrectionFactor2", "Applied #eta correction factor for the k_{T} background 2", "", 500, 0.5, 1.5, "Correction factor","dN^{Jets}/df"); - AddHistogram1D("hVertexZ", "Z distribution of the vertex", "", 400, -40., 40., "#Delta z(cm)","dN^{Events}/dz"); - - AddHistogram1D("hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr"); - AddHistogram1D("hCentrality", "Centrality distribution", "", 5, 0., 100., "Centrality (classes)","dN^{Events}"); - AddHistogram2D("hTrackCountAcc", "Number of tracks in acceptance vs. centrality", "LEGO2", 750, 0., 750., 5, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}"); - AddHistogram2D("hTrackPhiEta", "Track angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dN^{Tracks}/(d#phi d#eta)"); - AddHistogram1D("hTrackPt", "Tracks p_{T} distribution", "", 20000, 0., 200., "p_{T} (GeV/c)","dN^{Tracks}/dp_{T}"); - AddHistogram1D("hTrackCharge", "Charge", "", 11, -5, 5, "Charge (e)","dN^{Tracks}/dq"); - AddHistogram1D("hTrackEta", "Track #eta distribution", "", 180, -fTrackEtaWindow, +fTrackEtaWindow, "#eta","dN^{Tracks}/d#eta"); - - AddHistogram2D("hClusterCountAcc", "Number of clusters in acceptance vs. centrality", "LEGO2", 750, 0., 750., 5, 0, 100, "N clusters","Centrality", "dN^{Events}/dN^{Clusters}"); - AddHistogram1D("hClusterE", "Clusters energy distribution", "", 20000, 0., 200., "p_{T} (GeV/c)","dN^{Cluster}/dp_{T}"); - } - - // NOTE: Pythia histograms - if (fAnalyzePythia) - { - AddHistogram1D("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}"); - AddHistogram1D("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}"); - AddHistogram1D("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin", "Trials"); - } + AddHistogram1D("hNumberEvents", "Number of events (0 = before, 1 = after vertex cuts)", "", 2, 0, 2, "#Delta z(cm)","N^{Events}/cut"); // NOTE: Jet histograms if (fAnalyzeJets) { // ######## Jet spectra - AddHistogram1D("hJetPt", "Jets p_{T} distribution", "", 1000, 0., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hJetPtBgrdSubtractedRC", "Jets p_{T} distribution, RC background subtracted", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hJetPtBgrdSubtractedKT", "Jets p_{T} distribution, KT background subtracted, corrected for eta dependence)", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hJetPtBgrdSubtractedKTNoEtaCorr", "Jets p_{T} distribution, KT background subtracted", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hJetPtBgrdSubtractedKT2", "Jets p_{T} distribution, KT background 2 subtracted, corrected for eta dependence)", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hJetPtBgrdSubtractedKT2NoEtaCorr", "Jets p_{T} distribution, KT background 2 subtracted", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - - AddHistogram1D("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, Track background subtracted", "", 500, -50., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - - - AddHistogram1D("hJetArea", "Jets area distribution", "", 200, 0., 2., "Area","dN^{Jets}/dA"); - AddHistogram2D("hJetPtArea", "Jets p_{T} distribution", "LEGO2", 1000, 0., 200.,100, 0., 1., "p_{T} (GeV/c)","Area","dN^{Jets}/(dp_{T}dA)"); - AddHistogram1D("hJetDeltaPhi", "Jets combinatorial #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)"); - AddHistogram2D("hJetDeltaPhiPt", "Jets combinatorial #Delta #phi vs. p_{T}", "LEGO2", 250, 0., TMath::Pi(), 20, 0.,100., "#Delta #phi","max(p_{T,1},p_{T,2}) (GeV/c)","dN^{Jets}/d(#Delta #phi)dp_{T}"); - AddHistogram1D("hLeadingJetDeltaPhi", "1st and 2nd leading jet #Delta #phi", "", 250, 0., TMath::Pi(), "#Delta #phi","dN^{Jets}/d(#Delta #phi)"); - AddHistogram2D("hLeadingJetDeltaPhiPt", "1st and 2nd leading jet #Delta #phi vs. p_{T}", "LEGO2", 250, 0., TMath::Pi(),20, 0.,100., "#Delta #phi","1st leading p_{T} (GeV/c)","dN^{Jets}/d(#Delta #phi)dp_{T}"); - AddHistogram2D("hJetPtEta", "Jets p_{T} distribution", "LEGO2", 1000, 0., 200.,100, -0.6, 0.6, "p_{T} (GeV/c)","#eta","dN^{Jets}/(dp_{T}d#eta)"); - AddHistogram2D("hJetPtPhi", "Jets p_{T} #phi distribution", "LEGO2", 1000, 0., 200.,100, 0.0, TMath::TwoPi(), "p_{T} (GeV/c)","#phi","dN^{Jets}/(dp_{T}d#phi)"); - AddHistogram2D("hJetPtCentrality", "Jets p_{T} distribution", "LEGO2", 1000, 0., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); - AddHistogram2D("hJetPhiEta", "Jets angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -0.6, 0.6, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); + AddHistogram2D("hJetPt", "Jets p_{T} distribution", "", 1000, 0., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedRC", "Jets p_{T} distribution, RC background subtracted", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKT", "Jets p_{T} distribution, KT background subtracted", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedTR", "Jets p_{T} distribution, TR background subtracted", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedRCNoEta", "Jets p_{T} distribution, RC background subtracted (no #eta correction)", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedKTNoEta", "Jets p_{T} distribution, KT background subtracted (no #eta correction)", "", 500, -50., 200., 5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + AddHistogram2D("hJetPtBgrdSubtractedTRNoEta", "Jets p_{T} distribution, TR background subtracted (no #eta correction)", "", 500, -50., 200.,5, 0, 100, "p_{T} (GeV/c)","Centrality","dN^{Jets}/dp_{T}"); + + // ######## Jet stuff 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}"); - AddHistogram1D("hLeadingJetPt", "Leading jet p_{T}", "", 500, 0, 100, "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); AddHistogram1D("hSecondLeadingJetPt", "Second Leading jet p_{T}", "", 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}"); - AddHistogram1D("hDijetLeadingJetPt", "Dijet leading jet 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}"); - - AddHistogram1D("hDijet2ConstituentsPt", "Dijet2 constituents p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram1D("hDijet2LeadingJetPt", "Dijet2 leading jet p_{T} distribution", "", 500, 0., 100., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}"); - AddHistogram2D("hDijet2PtCorrelation", "Dijet2 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}"); - + 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)"); } + // NOTE: Jet background histograms if (fAnalyzeBackground) { - // ########## Delta Pt - AddHistogram1D("hDeltaPtKT", "Background fluctuations #delta p_{T} (KT, 0 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtKT1Excl", "Background fluctuations #delta p_{T} (KT, 1 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtKT2Excl", "Background fluctuations #delta p_{T} (KT, 2 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - - Double_t dptEtaMin = -(fTrackEtaWindow-fRandConeRadius) + 2*(fTrackEtaWindow-fRandConeRadius)/fBackgroundEtaBins * fKTDeltaPtEtaBin; - Double_t dptEtaMax = -(fTrackEtaWindow-fRandConeRadius) + 2*(fTrackEtaWindow-fRandConeRadius)/fBackgroundEtaBins * (fKTDeltaPtEtaBin+1); - - AddHistogram1D("hDeltaPtKTEta", Form("Background fluctuations #delta p_{T} (KT, 0 jets excluded, #eta=%1.3f to %1.3f)", dptEtaMin,dptEtaMax), "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtKTEta1Excl", Form("Background fluctuations #delta p_{T} (KT, 1 jets excluded, #eta=%1.3f to %1.3f)", dptEtaMin,dptEtaMax), "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtKTEta2Excl", Form("Background fluctuations #delta p_{T} (KT, 2 jets excluded, #eta=%1.3f to %1.3f)", dptEtaMin,dptEtaMax), "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - - AddHistogram1D("hDeltaPtKT2Eta2Excl", Form("Background fluctuations #delta p_{T} (KT 2, 2 jets excluded, #eta=%1.3f to %1.3f)", dptEtaMin,dptEtaMax), "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - - AddHistogram1D("hDeltaPtRC", "Background fluctuations #delta p_{T} (RC, 0 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtRC1Excl", "Background fluctuations #delta p_{T} (RC, 1 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtRC2Excl", "Background fluctuations #delta p_{T} (RC, 2 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - - AddHistogram1D("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR, 0 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtTR1Excl", "Background fluctuations #delta p_{T} (TR, 1 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - AddHistogram1D("hDeltaPtTR2Excl", "Background fluctuations #delta p_{T} (TR, 2 jets excluded)", "", 500, -20., 80., "#delta p_{T} (GeV/c)","dN^{Jets}/d#delta p_{T}"); - - - - AddHistogram2D("hKTJetPhiEta", "KT Jets angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -0.6, 0.6, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); - AddHistogram2D("hKTLeadingJetPhiEta", "KT Leading jets angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -0.6, 0.6, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); - - AddHistogram1D("hDijetBackground", "Background density (dijets excluded)", "", 400, 0., 40., "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram1D("hDijetBackgroundMostCentral", "Background density (0-20 centrality, dijets excluded)", "", 400, 0., 40., "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hDijetBackgroundVsCentrality", "Background density vs. centrality (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - AddHistogram1D("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 400, 0., 40., "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram1D("hDijetBackgroundPerpendicularMostCentral", "Background density (0-20 centrality, dijets excluded)", "", 400, 0., 40., "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hDijetBackgroundPerpendicularVsCentrality", "Background density vs. centrality (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - AddHistogram2D("hRCBackground", "RC background density (2 leading jets excluded, mean(8 RCs))", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); + // ########## Different background estimates + AddHistogram2D("hRCBackground", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hKTBackground", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hTRBackground", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - AddHistogram1D("hAccConesInRCBackground", Form("Number of cones used for RC background (|#eta| < %1.1f)", fSignalJetEtaWindow), "", 8, 0, 8, "Used cones", "dN^{Events}/dN^{Cones}"); - - AddHistogram2D("hRCBackgroundMostCentral", "RC background density (0-20 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hRCBackgroundMostPeripheral", "RC background density (80-100 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hRCBackgroundVsCentrality", "RC background density vs centrality (2 leading jets excluded)", "LEGO2", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - - AddHistogram2D("hKTBackground", "KT background density (2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundMostCentral", "KT background density (0-20 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundMostPeripheral", "KT background density (80-100 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hKTBackgroundVsCentrality", "KT background density vs centrality (2 leading jets excluded)", "LEGO2", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - AddHistogram2D("hKTBackground2", "KT background 2 density (2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hKTBackground2MostCentral", "KT background 2 density (0-20 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hKTBackground2MostPeripheral", "KT background 2 density (80-100 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hKTBackground2VsCentrality", "KT background 2 density vs centrality (2 leading jets excluded)", "LEGO2", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); - - - AddHistogram2D("hTrackBackground", "Track background density (2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hTrackBackgroundMostCentral", "Track background density (0-20 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hTrackBackgroundMostPeripheral", "Track background density (80-100 centrality, 2 leading jets excluded)", "LEGO2", fBackgroundEtaBins, -(fTrackEtaWindow-fRandConeRadius), +(fTrackEtaWindow-fRandConeRadius), 400, 0., 40., "#eta", "#rho (GeV/c)","dN^{Events}/d#rho"); - AddHistogram2D("hTrackBackgroundVsCentrality", "Track background density vs centrality (2 leading jets excluded)", "LEGO2", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + // ########## Delta Pt + AddHistogram2D("hDeltaPtKT", "Background fluctuations #delta p_{T} (KT)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtRC", "Background fluctuations #delta p_{T} (RC)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtTR", "Background fluctuations #delta p_{T} (TR)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTNoEta", "Background fluctuations #delta p_{T} (KT, no #eta correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtRCNoEta", "Background fluctuations #delta p_{T} (RC, no #eta correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtTRNoEta", "Background fluctuations #delta p_{T} (TR, no #eta correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtKTNoEtaNoExcl", "Background fluctuations #delta p_{T} (KT, no #eta correction, no leading jet correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtRCNoEtaNoExcl", "Background fluctuations #delta p_{T} (RC, no #eta correction, no leading jet correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + AddHistogram2D("hDeltaPtTRNoEtaNoExcl", "Background fluctuations #delta p_{T} (TR, no #eta correction, no leading jet correction)", "", 600, -40., 80., 5, 0, 100, "#delta p_{T} (GeV/c)","Centrality","dN^{Jets}/d#delta p_{T}"); + + // ########## Min bias background in eta bins + AddHistogram2D("hRCBackgroundEta", "RC background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta"); + AddHistogram2D("hKTBackgroundEta", "KT background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta"); + AddHistogram2D("hTRBackgroundEta", "TR background density (2 leading jets excluded)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta"); + + AddHistogram2D("hRCBackgroundEtaCorrected", "RC background density (2 leading jets excluded, #eta-corrected)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta"); + AddHistogram2D("hKTBackgroundEtaCorrected", "KT background density (2 leading jets excluded, #eta-corrected)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta"); + AddHistogram2D("hTRBackgroundEtaCorrected", "TR background density (2 leading jets excluded, #eta-corrected)", "LEGO2", 400, 0., 40., 5, -0.5, +0.5, "#rho (GeV/c)","#eta", "dN^{Events}/d#rho d#eta"); + + // ########## 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}"); + AddHistogram2D("hDijetBackground", "Background density (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + AddHistogram2D("hDijetBackgroundPerpendicular", "Background density (dijets excluded)", "", 200, 0., 20., 5, 0, 100, "#rho (GeV/c)","Centrality", "dN^{Events}/d#rho"); + } + // NOTE: Pythia histograms + if (fAnalyzePythia) + { + AddHistogram1D("hPythiaPtHard", "Pythia p_{T} hard distribution", "", 2000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}"); + AddHistogram1D("hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin","dN^{Events}/dp_{T,hard}"); + AddHistogram1D("hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins+2, -1, fNumPtHardBins+1, "p_{T} hard bin", "Trials"); } + // Save eta correction histograms to output when available + if(fJetKTEtaCorrection) + { + fHistList->Add(fJetKTEtaCorrection); + fHistCount++; + } + if(fJetRCEtaCorrection) + { + fHistList->Add(fJetRCEtaCorrection); + fHistCount++; + } + if(fJetTREtaCorrection) + { + fHistList->Add(fJetTREtaCorrection); + fHistCount++; + } // register Histograms @@ -158,33 +147,26 @@ void AliAnalysisTaskChargedJetsPA::Init() } //________________________________________________________________________ -AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* clusterArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeQA(1), fAnalyzeJets(1), fAnalyzeBackground(1), fAnalyzePythia(0), fHasTracks(0), fHasClusters(0), fHasJets(0), fHasBackgroundJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fClusterArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fClusterArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fKTDeltaPtEtaBin(3), fTrackBackgroundConeRadius(0.4), fNumberRandCones(8), fNumberExcludedJets(2), fDijetMaxAngleDeviation(10.0), fBackgroundEtaBins(5), fJetBgrdCorrectionFactors(0), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fClusterEtaWindow(0.7), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinClusterPt(0.300), fMinJetPt(1.0), fMinJetArea(0.4), fMinBackgroundJetPt(0.15), fMinDijetLeadingPt(10.0), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0) +AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA() : AliAnalysisTaskSE("AliAnalysisTaskChargedJetsPA"), fOutputList(0), fAnalyzeJets(1), fAnalyzeBackground(1), fAnalyzePythia(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.4), fNumberRandCones(8), fNumberExcludedJets(2), fDijetMaxAngleDeviation(10.0), fJetKTEtaCorrection(0), fJetRCEtaCorrection(0), fJetTREtaCorrection(0), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinJetPt(1.0), fMinJetArea(0.4), fMinBackgroundJetPt(0.15), fMinDijetLeadingPt(10.0), 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) +{ +// default constructor +} + +//________________________________________________________________________ +AliAnalysisTaskChargedJetsPA::AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeJets(1), fAnalyzeBackground(1), fAnalyzePythia(0), fHasTracks(0), fHasJets(0), fHasBackgroundJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fTRBackgroundConeRadius(0.4), fNumberRandCones(8), fNumberExcludedJets(2), fDijetMaxAngleDeviation(10.0), fJetKTEtaCorrection(0), fJetRCEtaCorrection(0), fJetTREtaCorrection(0), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinJetPt(1.0), fMinJetArea(0.4), fMinBackgroundJetPt(0.15), fMinDijetLeadingPt(10.0), 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) { #ifdef DEBUGMODE AliInfo("Calling constructor."); #endif - // Constructor - // Define input and output slots here (never in the dummy constructor) - // Input slot #0 works with a TChain - it is connected to the default input container - // Output slot #1 writes into a TH1 container - // Constructor - // Every instance of this task gets his own number static Int_t instance = 0; fTaskInstanceCounter = instance; instance++; fTrackArrayName = new TString(trackArrayName); - fClusterArrayName = new TString(clusterArrayName); - if (strcmp(fTrackArrayName->Data(),"") == 0) - fAnalyzeQA = kFALSE; - else - { - fAnalyzeQA = kTRUE; - if (fTrackArrayName->Contains("MCParticles")) //TODO: Hardcoded for now - fIsMC = kTRUE; - } + if (fTrackArrayName->Contains("MCParticles")) //TODO: Not working for now + fIsMC = kTRUE; fJetArrayName = new TString(jetArrayName); if (strcmp(fJetArrayName->Data(),"") == 0) @@ -223,6 +205,7 @@ inline Double_t AliAnalysisTaskChargedJetsPA::GetConePt(Double_t eta, Double_t p return tmpConePt; } + //________________________________________________________________________ inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard() { @@ -241,18 +224,19 @@ inline Double_t AliAnalysisTaskChargedJetsPA::GetPtHard() return tmpPtHard; } + //________________________________________________________________________ inline Int_t AliAnalysisTaskChargedJetsPA::GetPtHardBin() { // ########## PT HARD BIN EDGES - const Int_t localkPtHardLowerEdges[] = { 0, 5,11,21,36,57, 84,117,152,191,234}; - const Int_t localkPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000}; + 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}; Int_t tmpPtHardBin = 0; Double_t tmpPtHard = GetPtHard(); for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++) - if (tmpPtHard >= localkPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < localkPtHardHigherEdges[tmpPtHardBin]) + if (tmpPtHard >= kPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < kPtHardHigherEdges[tmpPtHardBin]) break; return tmpPtHardBin; @@ -288,20 +272,6 @@ inline Bool_t AliAnalysisTaskChargedJetsPA::IsTrackInAcceptance(AliVParticle* tr return kFALSE; } -//________________________________________________________________________ -inline Bool_t AliAnalysisTaskChargedJetsPA::IsClusterInAcceptance(AliVCluster* cluster) -{ - if (cluster != 0) - // if (TMath::Abs(cluster->Eta()) <= fClusterEtaWindow) -// if (cluster->Phi() <= 187.0/360.0 * TMath::TwoPi()); -// if (cluster->Phi() >= 80.0/360.0 * TMath::TwoPi()); - if (cluster->E() >= fMinClusterPt) - return kTRUE; - - return kFALSE; -} - - //________________________________________________________________________ inline Bool_t AliAnalysisTaskChargedJetsPA::IsBackgroundJetInAcceptance(AliEmcalJet *jet) { @@ -366,27 +336,6 @@ void AliAnalysisTaskChargedJetsPA::ExecOnce() } } } - // Check for cluster array - if (strcmp(fClusterArrayName->Data(), "") != 0) - { - fClusterArray = dynamic_cast(InputEvent()->FindListObject(fClusterArrayName->Data())); - fHasClusters = kTRUE; - if (!fClusterArray) - { - AliInfo(Form("%s: Could not retrieve clusters %s! This is OK, if clusters are not demanded.", GetName(), fClusterArrayName->Data())); - fHasClusters = kFALSE; - } - else - { - TClass *cl = fClusterArray->GetClass(); - if (!cl->GetBaseClass("AliVCluster")) - { - AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fClusterArrayName->Data())); - fClusterArray = 0; - fHasClusters = kFALSE; - } - } - } // Check for jet array if (strcmp(fJetArrayName->Data(), "") != 0) @@ -423,10 +372,9 @@ void AliAnalysisTaskChargedJetsPA::ExecOnce() } // Look, if initialization is OK - if ((!fHasTracks && fAnalyzeQA) || (!fHasTracks && fAnalyzeBackground)) + if (!fHasTracks && fAnalyzeBackground) { - AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating QA and background analysis",GetName())); - fAnalyzeQA = kFALSE; + AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating background analysis",GetName())); fAnalyzeBackground = kFALSE; } if ((!fHasJets && fAnalyzeJets) || (!fHasJets && fAnalyzeBackground)) @@ -552,22 +500,54 @@ Int_t AliAnalysisTaskChargedJetsPA::GetLeadingJets(TClonesArray* jetArray, Int_t } return jetCountAccepted; } + //________________________________________________________________________ -Double_t AliAnalysisTaskChargedJetsPA::GetJetBackgroundCorrFactor(Double_t eta, Double_t background) +Double_t AliAnalysisTaskChargedJetsPA::GetBackgroundEtaCorrFactor(EtaCorrectionMode mode, Double_t eta, Double_t background) { - Double_t tmpCorrFactor = 1.0; + if ((eta>=-0.5) && (eta<-0.3)) + return GetBackgroundEtaBinCorrFactor(mode, 1, background); + else if ((eta>=-0.3) && (eta<-0.1)) + return GetBackgroundEtaBinCorrFactor(mode, 2, background); + else if ((eta>=-0.1) && (eta<+0.1)) + return GetBackgroundEtaBinCorrFactor(mode, 3, background); + else if ((eta>=+0.1) && (eta<+0.3)) + return GetBackgroundEtaBinCorrFactor(mode, 4, background); + else if ((eta>=+0.3) && (eta<=+0.5)) + return GetBackgroundEtaBinCorrFactor(mode, 5, background); + else + AliError(Form("Wrong eta value! Eta=%1.4f", eta)); - if(fJetBgrdCorrectionFactors) - tmpCorrFactor = fJetBgrdCorrectionFactors->GetBinContent - ( - fJetBgrdCorrectionFactors->GetXaxis()->FindBin(eta), - fJetBgrdCorrectionFactors->GetYaxis()->FindBin(background) - ); + return 1.0; +} + +//________________________________________________________________________ +Double_t AliAnalysisTaskChargedJetsPA::GetBackgroundEtaBinCorrFactor(EtaCorrectionMode mode, Int_t eta, Double_t background) +{ + Double_t corrFactor = 1.0; - return tmpCorrFactor; + if((eta < 1) || (eta>5)) + { + AliError("Wrong eta bin!"); + return corrFactor; + } + // if background not valid, do not correct + if(background<0) + return corrFactor; + + if(mode == kKTEtaCorrection) + corrFactor = fJetKTEtaCorrection->GetBinContent(eta, fJetKTEtaCorrection->GetYaxis()->FindBin(background)); + else if(mode == kRCEtaCorrection) + corrFactor = fJetRCEtaCorrection->GetBinContent(eta, fJetRCEtaCorrection->GetYaxis()->FindBin(background)); + else if(mode == kTREtaCorrection) + corrFactor = fJetTREtaCorrection->GetBinContent(eta, fJetTREtaCorrection->GetYaxis()->FindBin(background)); + else if(mode == kNoEtaCorrection) + corrFactor = 1.0; + + return corrFactor; } + //________________________________________________________________________ -Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background, Bool_t useEtaCorrection) +Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Double_t background, EtaCorrectionMode mode) { #ifdef DEBUGMODE AliInfo("Getting corrected jet spectra."); @@ -581,12 +561,9 @@ Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Doubl Double_t correctedPt = -1.0; - // Get correction factor from saved histo or similar in dependence of jet eta and background density + // Get correction factor from saved histo in dependence of jet eta and background density Double_t corrfactor = 1.0; - if(useEtaCorrection) - { - corrfactor = GetJetBackgroundCorrFactor(jet->Eta(), background); - } + corrfactor = GetBackgroundEtaCorrFactor(mode, jet->Eta(), background); // Get Eta corrected background Double_t tmpCorrectedBackground = background * corrfactor; @@ -602,58 +579,36 @@ Double_t AliAnalysisTaskChargedJetsPA::GetCorrectedJetPt(AliEmcalJet* jet, Doubl } + //________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t& deltaPt, Double_t rho, Int_t numberExcludeLeadingJets, Int_t usedEtaBin, Bool_t useEtaCorrection) +void AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t& deltaPt, Double_t rho, EtaCorrectionMode mode, Bool_t leadingJetExclusion) { #ifdef DEBUGMODE AliInfo("Getting Delta Pt."); #endif - // Define the tmp delta pt + // Define an invalid delta pt deltaPt = -10000.0; - // Exclude UP TO numberExcludeLeadingJets - if (fNumberSignalJets < 2) - numberExcludeLeadingJets = fNumberSignalJets; - + // Define eta range Double_t etaMin, etaMax; - if (usedEtaBin==-1) - { - etaMin = -(fTrackEtaWindow-fRandConeRadius); - etaMax = +(fTrackEtaWindow-fRandConeRadius); - } - else - { - etaMin = -(fTrackEtaWindow-fRandConeRadius) + 2*(fTrackEtaWindow-fRandConeRadius)/fBackgroundEtaBins * usedEtaBin; - etaMax = -(fTrackEtaWindow-fRandConeRadius) + 2*(fTrackEtaWindow-fRandConeRadius)/fBackgroundEtaBins * (usedEtaBin+1); - } - + etaMin = -(fTrackEtaWindow-fRandConeRadius); + etaMax = +(fTrackEtaWindow-fRandConeRadius); + // Define random cone + Bool_t coneValid = kTRUE; Double_t tmpRandConeEta = 0.0; Double_t tmpRandConePhi = 0.0; - - Bool_t coneValid = kTRUE; - - tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin); tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi(); - // Apply eta correction on demand - if(useEtaCorrection) - rho = GetJetBackgroundCorrFactor(tmpRandConeEta, rho)*rho; + // Apply eta correction on background if demanded + rho *= GetBackgroundEtaCorrFactor(mode, tmpRandConeEta, rho); - // Go through all excluded leading jets and check if there's an overlap - for(Int_t j=0;jPhi(); Double_t excludedJetEta = tmpJet->Eta(); Double_t tmpDeltaPhi = GetDeltaPhi(tmpRandConePhi, excludedJetPhi); @@ -661,14 +616,19 @@ void AliAnalysisTaskChargedJetsPA::GetDeltaPt(Double_t& deltaPt, Double_t rho, I // Check, if cone has overlap with jet if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(tmpRandConeEta-excludedJetEta)*TMath::Abs(tmpRandConeEta-excludedJetEta) <= fRandConeRadius*fRandConeRadius) { - coneValid = kFALSE; - break; + // Define probability to exclude the RC + Double_t probability = (fNumberSignalJets-1)/fNumberSignalJets; + + // Only exclude cone with a given probability + if (fRandom->Rndm()<=probability) + coneValid = kFALSE; } } // Get the cones' pt and calculate delta pt if (coneValid) deltaPt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius) - (rho*fRandConeRadius*fRandConeRadius*TMath::Pi()); + #ifdef DEBUGMODE AliInfo("Got Delta Pt."); #endif @@ -743,86 +703,6 @@ void AliAnalysisTaskChargedJetsPA::GetKTBackgroundDensity(Int_t numberExcludeLea #endif } -//________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetKTBackground2Density(Int_t numberExcludeLeadingJets, Double_t& rhoMedian, Double_t& areaMean, Double_t etaMin, Double_t etaMax) -{ - #ifdef DEBUGMODE - AliInfo("Getting KT background 2 density."); - #endif - - // static declaration. Advantage: more speed. Disadvantage: Problematic for events with more than 1024 jets :) - static Double_t tmpRhos[1024]; - static Double_t tmpAreas[1024]; - - // Setting invalid values - rhoMedian = -1.0; - areaMean= -1.0; - - if ((etaMin == 0) && (etaMax == 0)) - { - etaMin = -fBackgroundJetEtaWindow; - etaMax = +fBackgroundJetEtaWindow; - } - - Int_t jetCountAccepted = 0; - Int_t jetCount = fBackgroundJetArray->GetEntries(); - - for (Int_t i = 0; i < jetCount; i++) - { - Bool_t jetValid = kTRUE; - AliEmcalJet* jet = static_cast(fBackgroundJetArray->At(i)); - if (!jet) - { - AliError(Form("%s: Could not receive jet %d", GetName(), i)); - continue; - } - - if (!((jet->Eta() >= etaMin) && (jet->Eta() < etaMax))) - continue; - if (!IsBackgroundJetInAcceptance(jet)) - continue; - - // Look, if theres an overlap of leading jets/ kT jet. If yes, exclude this jet - for(Int_t j=0;jPhi(), tmpLeadingJet->Phi()); - if ( tmpDeltaPhi*tmpDeltaPhi + TMath::Abs(jet->Eta()-tmpLeadingJet->Eta())*TMath::Abs(jet->Eta()-tmpLeadingJet->Eta()) <= fBackgroundJetRadius*fBackgroundJetRadius) - { - jetValid = kFALSE; - break; - } - } - } - - if(!jetValid) - continue; - - tmpRhos[jetCountAccepted] = jet->Pt() / jet->Area(); - tmpAreas[jetCountAccepted] = jet->Area(); - jetCountAccepted++; - } - - if (jetCountAccepted > 0) - { - rhoMedian = TMath::Median(jetCountAccepted, tmpRhos); - areaMean = TMath::Mean(jetCountAccepted, tmpAreas); - } - #ifdef DEBUGMODE - AliInfo("Got KT background 2 density."); - #endif -} - //________________________________________________________________________ Int_t AliAnalysisTaskChargedJetsPA::GetRCBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& rhoMedian, Double_t etaMin, Double_t etaMax, Int_t numberRandCones) @@ -915,10 +795,10 @@ Int_t AliAnalysisTaskChargedJetsPA::GetRCBackgroundDensity(Int_t numberExcludeLe } //________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetTrackBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, Double_t etaMin, Double_t etaMax) +void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, Double_t etaMin, Double_t etaMax) { #ifdef DEBUGMODE - AliInfo("Getting track background density."); + AliInfo("Getting TR background density."); #endif Double_t summedTracksPt = 0.0; @@ -956,7 +836,7 @@ void AliAnalysisTaskChargedJetsPA::GetTrackBackgroundDensity(Int_t numberExclude else AliFatal("Trying to exclude more than 2 jets in track background -- not implemented."); - if (IsTrackInCone(tmpTrack, tmpJet->Eta(), tmpJet->Phi(), fTrackBackgroundConeRadius)) + if (IsTrackInCone(tmpTrack, tmpJet->Eta(), tmpJet->Phi(), fTRBackgroundConeRadius)) { trackValid = kFALSE; break; @@ -979,24 +859,24 @@ void AliAnalysisTaskChargedJetsPA::GetTrackBackgroundDensity(Int_t numberExclude // Now: exclude the part of the leading jet that is in the strip if (numberExcludeLeadingJets == 2) - tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTrackBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()) -MCGetOverlapCircleRectancle(fSecondLeadingJet->Eta(), fSecondLeadingJet->Phi(), fTrackBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi())); + tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi()) -MCGetOverlapCircleRectancle(fSecondLeadingJet->Eta(), fSecondLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi())); else if (numberExcludeLeadingJets == 1) - tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTrackBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi())); + tmpArea = tmpArea*(1.0-MCGetOverlapCircleRectancle(fFirstLeadingJet->Eta(), fFirstLeadingJet->Phi(), fTRBackgroundConeRadius, etaMin, etaMax, 0., TMath::TwoPi())); rhoMean = summedTracksPt/tmpArea; area = tmpArea; } #ifdef DEBUGMODE - AliInfo("Got track background density."); + AliInfo("Got TR background density."); #endif } //________________________________________________________________________ -void AliAnalysisTaskChargedJetsPA::GetTrackBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular) +void AliAnalysisTaskChargedJetsPA::GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular) { #ifdef DEBUGMODE - AliInfo("Getting track background density."); + AliInfo("Getting TR background density."); #endif // Setting invalid values @@ -1048,7 +928,7 @@ void AliAnalysisTaskChargedJetsPA::GetTrackBackgroundDensity(Int_t numberExclude } #ifdef DEBUGMODE - AliInfo("Got track background density."); + AliInfo("Got TR background density."); #endif } @@ -1060,14 +940,6 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) #endif ////////////////////// NOTE: initialization & casting - if (!event) { - AliError("??? Event pointer == 0 ???"); - return; - } - - if (!fInitialized) - ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos - // Additional cuts FillHistogram("hNumberEvents", 0.5); // number of events before manual cuts @@ -1082,188 +954,88 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) ////////////////////// NOTE: Get Centrality, (Leading)Signal jets and Background - // Get centrality (V0A) + // Get centrality AliCentrality* tmpCentrality = NULL; tmpCentrality = event->GetCentrality(); Double_t centralityPercentile = 0.0; if (tmpCentrality != NULL) - centralityPercentile = tmpCentrality->GetCentralityPercentile("V0A"); + centralityPercentile = tmpCentrality->GetCentralityPercentile(fCentralityType.Data()); // Get jets if (fAnalyzeBackground || fAnalyzeJets) GetSignalJets(); - // Get background - - // Background with N excluded leading jets - std::vector ktBackgroundRhoMedian(fBackgroundEtaBins+1); - std::vector ktBackgroundAreaMean(fBackgroundEtaBins+1); - std::vector ktBackground2RhoMedian(fBackgroundEtaBins+1); - std::vector ktBackground2AreaMean(fBackgroundEtaBins+1); - std::vector rcBackgroundRhoMean(fBackgroundEtaBins+1); - std::vector rcBackgroundRhoMedian(fBackgroundEtaBins+1); - std::vector trackBackgroundRhoMean(fBackgroundEtaBins+1); - std::vector trackBackgroundArea(fBackgroundEtaBins+1); - Double_t dijetBackground = -1.0; // calculation only done in events with dijets I! - Double_t dijetBackgroundPerpendicular = -1.0; // calculation only done in events with dijets I! + // Get background estimates + Double_t backgroundKTMedian = -1.0; + Double_t backgroundRCMean = -1.0; + Double_t backgroundRCMedian = -1.0; + Double_t backgroundTRMean = -1.0; + Double_t backgroundKTAreaMean = -1.0; + Double_t backgroundTRAreaMean = -1.0; + Double_t dijetBackground = -1.0; + Double_t dijetBackgroundPerpendicular = -1.0; + if (fAnalyzeBackground) { - - // Get backgrounds in bins of eta - for(Int_t i = 0; iGetPrimaryVertex()->GetZ()); - FillHistogram("hVertexR",TMath::Sqrt(event->GetPrimaryVertex()->GetX()*event->GetPrimaryVertex()->GetX() + event->GetPrimaryVertex()->GetY()*event->GetPrimaryVertex()->GetY())); - FillHistogram("hCentrality",centralityPercentile); - - Int_t trackCountAcc = 0; - Int_t nTracks = fTrackArray->GetEntries(); - for (Int_t i = 0; i < nTracks; i++) + // ### SIGNAL JET ANALYSIS + for (Int_t i = 0; i(fTrackArray->At(i)); - if (IsTrackInAcceptance(track)) - { - FillHistogram("hTrackPhiEta", track->Phi(),track->Eta(), 1); - FillHistogram("hTrackPt", track->Pt()); - FillHistogram("hTrackEta", track->Eta()); - FillHistogram("hTrackCharge", track->Charge()); - trackCountAcc++; - } - } - FillHistogram("hTrackCountAcc", trackCountAcc, centralityPercentile); + AliEmcalJet* tmpJet = fSignalJets[i]; - if (fHasClusters) - { - Int_t clusterCountAcc = 0; - Int_t nClusters = fClusterArray->GetEntries(); - for (Int_t i = 0; i < nClusters; i++) - { - AliVCluster* cluster = static_cast(fClusterArray->At(i)); - if (IsClusterInAcceptance(cluster)) - { - FillHistogram("hClusterE", cluster->E()); - clusterCountAcc++; - } - } - FillHistogram("hClusterCountAcc", clusterCountAcc, centralityPercentile); - } - } - #ifdef DEBUGMODE - AliInfo("Calculate()::QA done."); - #endif + // Jet spectra + FillHistogram("hJetPt", tmpJet->Pt(), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedRC", GetCorrectedJetPt(tmpJet, backgroundRCMean, kRCEtaCorrection), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKT", GetCorrectedJetPt(tmpJet, backgroundKTMedian, kKTEtaCorrection), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, backgroundTRMean, kTREtaCorrection), centralityPercentile); - ////////////////////// NOTE: Jet analysis and calculations + FillHistogram("hJetPtBgrdSubtractedRCNoEta", GetCorrectedJetPt(tmpJet, backgroundRCMean, kNoEtaCorrection), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedKTNoEta", GetCorrectedJetPt(tmpJet, backgroundKTMedian, kNoEtaCorrection), centralityPercentile); + FillHistogram("hJetPtBgrdSubtractedTRNoEta", GetCorrectedJetPt(tmpJet, backgroundTRMean, kNoEtaCorrection), centralityPercentile); - if (fAnalyzeJets) - { - FillHistogram("hJetCountAll", fJetArray->GetEntries()); - FillHistogram("hJetCountAccepted", fNumberSignalJets); - if (fFirstLeadingJet) - FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt()); - if (fSecondLeadingJet) - FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt()); + // Signal jet vs. signal jet - "Combinatorial" + for (Int_t j = i+1; jPhi(), fSignalJets[j]->Phi())); + } - // ### Dijets I ### + // ### DIJETS if(fNumberSignalJets >= 2) { FillHistogram("hLeadingJetDeltaPhi", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi())); - FillHistogram("hLeadingJetDeltaPhiPt", GetDeltaPhi(fFirstLeadingJet->Phi(), fSecondLeadingJet->Phi()), fFirstLeadingJet->Pt()); - if (IsDijet(fFirstLeadingJet, fSecondLeadingJet)) // Gettin' the money + if (IsDijet(fFirstLeadingJet, fSecondLeadingJet)) { - FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt()); FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt()); + FillHistogram("hDijetConstituentsPt", fFirstLeadingJet->Pt()); + FillHistogram("hDijetConstituentsPt", fSecondLeadingJet->Pt()); + FillHistogram("hDijetLeadingJetPt", fFirstLeadingJet->Pt()); FillHistogram("hDijetPtCorrelation", fFirstLeadingJet->Pt(), fSecondLeadingJet->Pt()); Double_t dummyArea = 0; - GetTrackBackgroundDensity (2, dijetBackground, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE); - GetTrackBackgroundDensity (2, dijetBackgroundPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE); + GetTRBackgroundDensity (2, dijetBackground, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kFALSE); + GetTRBackgroundDensity (2, dijetBackgroundPerpendicular, dummyArea, fFirstLeadingJet, fSecondLeadingJet, kTRUE); } } - // SIGNAL JET ANALYSIS - for (Int_t i = 0; iPt(), tmpJet->Area()); - FillHistogram("hJetPtEta", tmpJet->Pt(), tmpJet->Eta()); - FillHistogram("hJetPtPhi", tmpJet->Pt(), tmpJet->Phi()); - FillHistogram("hJetPtCentrality", tmpJet->Pt(), centralityPercentile); - FillHistogram("hJetArea", tmpJet->Area()); - FillHistogram("hJetPt", tmpJet->Pt()); - FillHistogram("hJetPhiEta", tmpJet->Phi(),tmpJet->Eta()); - - // Background subtracted spectra - - FillHistogram("hJetPtBgrdSubtractedRC", GetCorrectedJetPt(tmpJet, rcBackgroundRhoMean[fBackgroundEtaBins])); - FillHistogram("hJetPtBgrdSubtractedKT", GetCorrectedJetPt(tmpJet, ktBackgroundRhoMedian[fBackgroundEtaBins], kTRUE)); - FillHistogram("hJetPtBgrdSubtractedKT2", GetCorrectedJetPt(tmpJet, ktBackground2RhoMedian[fBackgroundEtaBins], kTRUE)); - FillHistogram("hJetPtBgrdSubtractedKTNoEtaCorr", GetCorrectedJetPt(tmpJet, ktBackgroundRhoMedian[fBackgroundEtaBins])); - FillHistogram("hJetPtBgrdSubtractedKT2NoEtaCorr", GetCorrectedJetPt(tmpJet, ktBackground2RhoMedian[fBackgroundEtaBins])); - FillHistogram("hJetPtBgrdSubtractedTR", GetCorrectedJetPt(tmpJet, trackBackgroundRhoMean[fBackgroundEtaBins])); - - - Double_t tmpCorrFactor = GetJetBackgroundCorrFactor(tmpJet->Eta(), ktBackgroundRhoMedian[fBackgroundEtaBins]); - FillHistogram("hAppliedEtaCorrectionFactor", tmpCorrFactor); - tmpCorrFactor = GetJetBackgroundCorrFactor(tmpJet->Eta(), ktBackground2RhoMedian[fBackgroundEtaBins]); - FillHistogram("hAppliedEtaCorrectionFactor2", tmpCorrFactor); - - // Signal jet vs. signal jet - for (Int_t j = i+1; jPhi(), tmpJet2->Phi())); - FillHistogram("hJetDeltaPhiPt", GetDeltaPhi(tmpJet->Phi(), tmpJet2->Phi()), max(tmpJet->Pt(), tmpJet2->Pt())); + // ### SOME JET PLOTS + FillHistogram("hJetCountAll", fJetArray->GetEntries()); + FillHistogram("hJetCountAccepted", fNumberSignalJets); + if (fFirstLeadingJet) + FillHistogram("hLeadingJetPt", fFirstLeadingJet->Pt()); + if (fSecondLeadingJet) + FillHistogram("hSecondLeadingJetPt", fSecondLeadingJet->Pt()); - // ### Dijets II ### - if (IsDijet(tmpJet, tmpJet2)) // Gettin' the money - { - FillHistogram("hDijet2ConstituentsPt", tmpJet->Pt()); FillHistogram("hDijet2ConstituentsPt", tmpJet2->Pt()); - FillHistogram("hDijet2LeadingJetPt", fFirstLeadingJet->Pt()); - FillHistogram("hDijet2PtCorrelation", tmpJet->Pt(), tmpJet2->Pt()); - } - } - } } //endif AnalyzeJets #ifdef DEBUGMODE @@ -1273,139 +1045,108 @@ void AliAnalysisTaskChargedJetsPA::Calculate(AliVEvent* event) if (fAnalyzeBackground) { + // Calculate (non-eta corrected) background in centrality classes + FillHistogram("hRCBackground", backgroundRCMean, centralityPercentile); + FillHistogram("hTRBackground", backgroundTRMean, centralityPercentile); + FillHistogram("hKTBackground", backgroundKTMedian, centralityPercentile); - Int_t leadingJetIds[] = {-1, -1}; - GetLeadingJets(fBackgroundJetArray, &leadingJetIds[0], kFALSE); - - for (Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++) + // Calculate backgrounds in eta bins + for(Int_t i=0;i<5;i++) { - AliEmcalJet* jet = static_cast(fBackgroundJetArray->At(i)); - if (!jet) - { - AliError(Form("%s: Could not receive kt jet %d", GetName(), i)); - continue; - } - if (!IsBackgroundJetInAcceptance(jet)) - continue; - if (!((jet->Eta() >= -fBackgroundJetEtaWindow) && (jet->Eta() < fBackgroundJetEtaWindow))) - continue; - - FillHistogram("hKTJetPhiEta", jet->Phi(),jet->Eta()); - if(i==leadingJetIds[0]) - FillHistogram("hKTLeadingJetPhiEta", jet->Phi(),jet->Eta()); - + Double_t dummy = 0.0; + Double_t tmpKTRho = 0.0; + Double_t tmpRCRho = 0.0; + Double_t tmpTRRho = 0.0; + Double_t tmpKTRhoCorr = 0.0; + Double_t tmpRCRhoCorr = 0.0; + Double_t tmpTRRhoCorr = 0.0; + + Double_t etaMin = -(fTrackEtaWindow-fSignalJetRadius) + 2*(fTrackEtaWindow-fSignalJetRadius)/5 * i; + Double_t etaMax = -(fTrackEtaWindow-fSignalJetRadius) + 2*(fTrackEtaWindow-fSignalJetRadius)/5 * (i+1); + + // Calculate backgrounds + GetKTBackgroundDensity (fNumberExcludedJets, tmpKTRho, dummy, etaMin, etaMax); + GetTRBackgroundDensity (fNumberExcludedJets, tmpTRRho, dummy, etaMin, etaMax); + GetRCBackgroundDensity (fNumberExcludedJets, tmpRCRho, dummy, etaMin, etaMax); + // Add eta-correction + tmpKTRhoCorr = tmpKTRho * GetBackgroundEtaCorrFactor(kKTEtaCorrection, (etaMin+etaMax)/2.0, tmpKTRho); + tmpTRRhoCorr = tmpTRRho * GetBackgroundEtaCorrFactor(kTREtaCorrection, (etaMin+etaMax)/2.0, tmpTRRho); + tmpRCRhoCorr = tmpRCRho * GetBackgroundEtaCorrFactor(kRCEtaCorrection, (etaMin+etaMax)/2.0, tmpRCRho); + + FillHistogram("hRCBackgroundEta", tmpRCRho, (etaMin+etaMax)/2.0); + FillHistogram("hTRBackgroundEta", tmpTRRho, (etaMin+etaMax)/2.0); + FillHistogram("hKTBackgroundEta", tmpKTRho, (etaMin+etaMax)/2.0); + FillHistogram("hRCBackgroundEtaCorrected", tmpRCRhoCorr, (etaMin+etaMax)/2.0); + FillHistogram("hTRBackgroundEtaCorrected", tmpTRRhoCorr, (etaMin+etaMax)/2.0); + FillHistogram("hKTBackgroundEtaCorrected", tmpKTRhoCorr, (etaMin+etaMax)/2.0); } - // ############# RC, Track, and KT background calculations - Double_t etaMin = 0; - for (Int_t i=0;i= 80.) - { - FillHistogram("hRCBackgroundMostPeripheral", etaMin, rcBackgroundRhoMean[i]); - FillHistogram("hTrackBackgroundMostPeripheral", etaMin, trackBackgroundRhoMean[i]); - FillHistogram("hKTBackgroundMostPeripheral", etaMin, ktBackgroundRhoMedian[i]); - FillHistogram("hKTBackground2MostPeripheral", etaMin, ktBackground2RhoMedian[i]); - } - } - - FillHistogram("hRCBackgroundVsCentrality", rcBackgroundRhoMean[fBackgroundEtaBins], centralityPercentile); - FillHistogram("hTrackBackgroundVsCentrality", trackBackgroundRhoMean[fBackgroundEtaBins], centralityPercentile); - FillHistogram("hKTBackgroundVsCentrality", ktBackgroundRhoMedian[fBackgroundEtaBins], centralityPercentile); - FillHistogram("hKTBackground2VsCentrality", ktBackground2RhoMedian[fBackgroundEtaBins], centralityPercentile); - + // In case of dijets -> look at the background if (dijetBackground >= 0) - { - // Background in Dijet events - FillHistogram("hDijetBackground", dijetBackground); - if(centralityPercentile <= 20.) - FillHistogram("hDijetBackgroundMostCentral", dijetBackground); - FillHistogram("hDijetBackgroundVsCentrality", dijetBackground, centralityPercentile); - } + FillHistogram("hDijetBackground", dijetBackground, centralityPercentile); if (dijetBackgroundPerpendicular >= 0) - { - // Background in Dijet events - FillHistogram("hDijetBackgroundPerpendicular", dijetBackgroundPerpendicular); - if(centralityPercentile <= 20.) - FillHistogram("hDijetBackgroundPerpendicularMostCentral", dijetBackgroundPerpendicular); - FillHistogram("hDijetBackgroundPerpendicularVsCentrality", dijetBackgroundPerpendicular, centralityPercentile); - } + FillHistogram("hDijetBackgroundPerpendicular", dijetBackgroundPerpendicular, centralityPercentile); + + + // Calculate the delta pt + Double_t tmpDeltaPtKT = 0.0; + Double_t tmpDeltaPtRC = 0.0; + Double_t tmpDeltaPtTR = 0.0; + Double_t tmpDeltaPtKTNoEta = 0.0; + Double_t tmpDeltaPtRCNoEta = 0.0; + Double_t tmpDeltaPtTRNoEta = 0.0; + Double_t tmpDeltaPtKTNoEtaNoExcl = 0.0; + Double_t tmpDeltaPtRCNoEtaNoExcl = 0.0; + Double_t tmpDeltaPtTRNoEtaNoExcl = 0.0; + GetDeltaPt(tmpDeltaPtKT, backgroundKTMedian, kKTEtaCorrection); + GetDeltaPt(tmpDeltaPtRC, backgroundRCMean, kRCEtaCorrection); + GetDeltaPt(tmpDeltaPtTR, backgroundTRMean, kTREtaCorrection); + GetDeltaPt(tmpDeltaPtKTNoEta, backgroundKTMedian, kNoEtaCorrection); + GetDeltaPt(tmpDeltaPtRCNoEta, backgroundRCMean, kNoEtaCorrection); + GetDeltaPt(tmpDeltaPtTRNoEta, backgroundTRMean, kNoEtaCorrection); + GetDeltaPt(tmpDeltaPtKTNoEtaNoExcl, backgroundKTMedian, kNoEtaCorrection, kFALSE); + GetDeltaPt(tmpDeltaPtRCNoEtaNoExcl, backgroundRCMean, kNoEtaCorrection, kFALSE); + GetDeltaPt(tmpDeltaPtTRNoEtaNoExcl, backgroundTRMean, kNoEtaCorrection, kFALSE); + + // If valid, fill the delta pt histograms + if(tmpDeltaPtKT > -10000.0) + FillHistogram("hDeltaPtKT", tmpDeltaPtKT, centralityPercentile); + if(tmpDeltaPtRC > -10000.0) + FillHistogram("hDeltaPtRC", tmpDeltaPtRC, centralityPercentile); + if(tmpDeltaPtTR > -10000.0) + FillHistogram("hDeltaPtTR", tmpDeltaPtTR, centralityPercentile); + if(tmpDeltaPtKTNoEta > -10000.0) + FillHistogram("hDeltaPtKTNoEta", tmpDeltaPtKTNoEta, centralityPercentile); + if(tmpDeltaPtRCNoEta > -10000.0) + FillHistogram("hDeltaPtRCNoEta", tmpDeltaPtRCNoEta, centralityPercentile); + if(tmpDeltaPtTRNoEta > -10000.0) + FillHistogram("hDeltaPtTRNoEta", tmpDeltaPtTRNoEta, centralityPercentile); + if(tmpDeltaPtKTNoEtaNoExcl > -10000.0) + FillHistogram("hDeltaPtKTNoEtaNoExcl", tmpDeltaPtKTNoEtaNoExcl, centralityPercentile); + if(tmpDeltaPtRCNoEtaNoExcl > -10000.0) + FillHistogram("hDeltaPtRCNoEtaNoExcl", tmpDeltaPtRCNoEtaNoExcl, centralityPercentile); + if(tmpDeltaPtTRNoEtaNoExcl > -10000.0) + FillHistogram("hDeltaPtTRNoEtaNoExcl", tmpDeltaPtTRNoEtaNoExcl, centralityPercentile); - // ########## Delta pT calculations (most central, kt is eta corrected) - if (centralityPercentile <= 20.) - { - Double_t tmpDeltaPtKT, tmpDeltaPtKT2Excl, tmpDeltaPtKT1Excl; - Double_t tmpDeltaPtKTEta, tmpDeltaPtKTEta2Excl, tmpDeltaPtKTEta1Excl, tmpDeltaPtKT2Eta2Excl; - Double_t tmpDeltaPtRC, tmpDeltaPtRC2Excl, tmpDeltaPtRC1Excl; - Double_t tmpDeltaPtTR, tmpDeltaPtTR2Excl, tmpDeltaPtTR1Excl; - - GetDeltaPt(tmpDeltaPtKT, ktBackgroundRhoMedian[fBackgroundEtaBins], 0, -1, kTRUE); - GetDeltaPt(tmpDeltaPtKTEta, ktBackgroundRhoMedian[fKTDeltaPtEtaBin], 0, fKTDeltaPtEtaBin); - GetDeltaPt(tmpDeltaPtRC, rcBackgroundRhoMean[fBackgroundEtaBins], 0); - GetDeltaPt(tmpDeltaPtTR, trackBackgroundRhoMean[fBackgroundEtaBins], 0); - - GetDeltaPt(tmpDeltaPtKT1Excl, ktBackgroundRhoMedian[fBackgroundEtaBins], 1, -1, kTRUE); - GetDeltaPt(tmpDeltaPtKTEta1Excl, ktBackgroundRhoMedian[fKTDeltaPtEtaBin], 1, fKTDeltaPtEtaBin); - GetDeltaPt(tmpDeltaPtRC1Excl, rcBackgroundRhoMean[fBackgroundEtaBins], 1); - GetDeltaPt(tmpDeltaPtTR1Excl, trackBackgroundRhoMean[fBackgroundEtaBins], 1); - - GetDeltaPt(tmpDeltaPtKT2Excl, ktBackgroundRhoMedian[fBackgroundEtaBins], 2, -1, kTRUE); - GetDeltaPt(tmpDeltaPtKTEta2Excl, ktBackgroundRhoMedian[fKTDeltaPtEtaBin], 2, fKTDeltaPtEtaBin); - GetDeltaPt(tmpDeltaPtRC2Excl, rcBackgroundRhoMean[fBackgroundEtaBins], 2); - GetDeltaPt(tmpDeltaPtTR2Excl, trackBackgroundRhoMean[fBackgroundEtaBins], 2); - - GetDeltaPt(tmpDeltaPtKT2Eta2Excl, ktBackground2RhoMedian[fKTDeltaPtEtaBin], 2, fKTDeltaPtEtaBin); - - // kT Background - if(tmpDeltaPtKT > -10000.0) - FillHistogram("hDeltaPtKT", tmpDeltaPtKT); - if(tmpDeltaPtKT1Excl > -10000.0) - FillHistogram("hDeltaPtKT1Excl", tmpDeltaPtKT1Excl); - if(tmpDeltaPtKT2Excl > -10000.0) - FillHistogram("hDeltaPtKT2Excl", tmpDeltaPtKT2Excl); - - if(tmpDeltaPtKT > -10000.0) - FillHistogram("hDeltaPtKTEta", tmpDeltaPtKTEta); - if(tmpDeltaPtKTEta1Excl > -10000.0) - FillHistogram("hDeltaPtKTEta1Excl", tmpDeltaPtKTEta1Excl); - if(tmpDeltaPtKTEta2Excl > -10000.0) - FillHistogram("hDeltaPtKTEta2Excl", tmpDeltaPtKTEta2Excl); - if(tmpDeltaPtKT2Eta2Excl > -10000.0) - FillHistogram("hDeltaPtKT2Eta2Excl", tmpDeltaPtKT2Eta2Excl); - - // RC Background - if(tmpDeltaPtRC > -10000.0) - FillHistogram("hDeltaPtRC", tmpDeltaPtRC); - if(tmpDeltaPtRC1Excl > -10000.0) - FillHistogram("hDeltaPtRC1Excl", tmpDeltaPtRC1Excl); - if(tmpDeltaPtRC2Excl > -10000.0) - FillHistogram("hDeltaPtRC2Excl", tmpDeltaPtRC2Excl); - // TR Background - if(tmpDeltaPtTR > -10000.0) - FillHistogram("hDeltaPtTR", tmpDeltaPtTR); - if(tmpDeltaPtTR1Excl > -10000.0) - FillHistogram("hDeltaPtTR1Excl", tmpDeltaPtTR1Excl); - if(tmpDeltaPtTR2Excl > -10000.0) - FillHistogram("hDeltaPtTR2Excl", tmpDeltaPtTR2Excl); - } } #ifdef DEBUGMODE AliInfo("Calculate()::Background done."); #endif -} + + ////////////////////// NOTE: Pythia histograms + if(fAnalyzePythia) + { + FillHistogram("hPythiaPtHard", GetPtHard()); + FillHistogram("hPythiaNTrials", GetPtHardBin()-0.1, fTrials); + FillHistogram("hPythiaXSection", GetPtHardBin()-0.1, fCrossSection); + + #ifdef DEBUGMODE + AliInfo("Calculate()::Pythia done."); + #endif + } + +} //________________________________________________________________________ Bool_t AliAnalysisTaskChargedJetsPA::Notify() @@ -1481,13 +1222,96 @@ Bool_t AliAnalysisTaskChargedJetsPA::Notify() return kTRUE; } +//________________________________________________________________________ +void AliAnalysisTaskChargedJetsPA::SetKTEtaCorrectionFactors(TH2D* histo) +{ + // COPY given histogram + fJetKTEtaCorrection = new TH2D(*histo); + + if (!fJetKTEtaCorrection) + AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName())); + + // Look, if given histogram is compatible with given code + if (fJetKTEtaCorrection->GetXaxis()->GetNbins() != 5) + AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetKTEtaCorrection->GetXaxis()->GetNbins())); +} + +//________________________________________________________________________ +void AliAnalysisTaskChargedJetsPA::SetRCEtaCorrectionFactors(TH2D* histo) +{ + // COPY given histogram + fJetRCEtaCorrection = new TH2D(*histo); + + if (!fJetRCEtaCorrection) + AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName())); + + // Look, if given histogram is compatible with given code + if (fJetRCEtaCorrection->GetXaxis()->GetNbins() != 5) + AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetRCEtaCorrection->GetXaxis()->GetNbins())); +} + +//________________________________________________________________________ +void AliAnalysisTaskChargedJetsPA::SetTREtaCorrectionFactors(TH2D* histo) +{ + // COPY given histogram + fJetTREtaCorrection = new TH2D(*histo); + + if (!fJetTREtaCorrection) + AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName())); + + // Look, if given histogram is compatible with given code + if (fJetTREtaCorrection->GetXaxis()->GetNbins() != 5) + AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetTREtaCorrection->GetXaxis()->GetNbins())); +} + +//________________________________________________________________________ +inline Double_t AliAnalysisTaskChargedJetsPA::EtaToTheta(Double_t arg) + {return 2.*atan(exp(-arg));} +//________________________________________________________________________ +inline Double_t AliAnalysisTaskChargedJetsPA::ThetaToEta(Double_t arg) +{ + if ((arg > TMath::Pi()) || (arg < 0.0)) + { + AliError(Form("ThetaToEta got wrong input! (%f)", arg)); + return 0.0; + } + return -log(tan(arg/2.)); +} +//________________________________________________________________________ +inline Double_t AliAnalysisTaskChargedJetsPA::GetDeltaPhi(Double_t phi1, Double_t phi2) + {return min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));} + +//________________________________________________________________________ +Double_t AliAnalysisTaskChargedJetsPA::MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax) +{ + const Int_t kTests = 1000; + Int_t hits = 0; + TRandom3 randomGen(0); + + // Loop over kTests-many tests + for (Int_t i=0; i(hits)/static_cast(kTests)); +} + //________________________________________________________________________ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double_t x) { TH1* tmpHist = static_cast(fOutputList->FindObject(GetHistoName(key))); if(!tmpHist) { - AliInfo(Form("Cannot find histogram <%s> ",key)) ; + AliWarning(Form("Cannot find histogram <%s> ",key)) ; return; } @@ -1500,7 +1324,7 @@ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double TH1* tmpHist = static_cast(fOutputList->FindObject(GetHistoName(key))); if(!tmpHist) { - AliInfo(Form("Cannot find histogram <%s> ",key)); + AliWarning(Form("Cannot find histogram <%s> ",key)); return; } @@ -1516,7 +1340,7 @@ inline void AliAnalysisTaskChargedJetsPA::FillHistogram(const char * key, Double TH2* tmpHist = static_cast(fOutputList->FindObject(GetHistoName(key))); if(!tmpHist) { - AliInfo(Form("Cannot find histogram <%s> ",key)); + AliWarning(Form("Cannot find histogram <%s> ",key)); return; } @@ -1601,6 +1425,15 @@ void AliAnalysisTaskChargedJetsPA::UserExec(Option_t *) AliInfo("UserExec() started."); #endif + if (!InputEvent()) + { + AliError("??? Event pointer == 0 ???"); + return; + } + + if (!fInitialized) + ExecOnce(); // Get tracks, jets, background from arrays if not already given + Init Histos + Calculate(InputEvent()); PostData(1, fOutputList); diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.h index 7bd2a6db57d..1a1cecc1755 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.h +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.h @@ -1,8 +1,7 @@ #ifndef ALIANALYSISTASKCHARGEDJETSPA_H #define ALIANALYSISTASKCHARGEDJETSPA_H -// #define DEBUGMODE - +//#define DEBUGMODE class TH1F; class TH2F; @@ -15,215 +14,110 @@ class AliVParticle; class AliLog; class AliAnalysisUtils; -#ifndef ALIANALYSISTASKSE_H -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "AliAnalysisTask.h" -#include "AliCentrality.h" -#include "AliStack.h" -#include "AliESDEvent.h" -#include "AliESDInputHandler.h" -#include "AliAODEvent.h" -#include "AliAODHandler.h" -#include "AliAnalysisManager.h" -#include "AliAnalysisTaskSE.h" -#endif - -#include -#include "AliGenPythiaEventHeader.h" -#include "AliMCEvent.h" -#include "AliLog.h" -#include -#include -#include "AliVEventHandler.h" -#include "AliVParticle.h" -#include "AliAnalysisUtils.h" - - class AliAnalysisTaskChargedJetsPA : public AliAnalysisTaskSE { public: - - AliAnalysisTaskChargedJetsPA() : AliAnalysisTaskSE(), fOutputList(0), fAnalyzeQA(1), fAnalyzeJets(1), fAnalyzeBackground(1), fAnalyzePythia(0), fHasTracks(0), fHasClusters(0), fHasJets(0), fHasBackgroundJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fClusterArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fClusterArrayName(0), fBackgroundJetArrayName(0), fNumPtHardBins(11), fRandConeRadius(0.4), fSignalJetRadius(0.4), fBackgroundJetRadius(0.4), fKTDeltaPtEtaBin(3), fTrackBackgroundConeRadius(0.4), fNumberRandCones(8), fNumberExcludedJets(2), fDijetMaxAngleDeviation(10.0), fBackgroundEtaBins(5), fJetBgrdCorrectionFactors(0), fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fClusterEtaWindow(0.7), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinClusterPt(0.300), fMinJetPt(1.0), fMinJetArea(0.4), fMinBackgroundJetPt(0.15), fMinDijetLeadingPt(10.0), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0) {} - - AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* clusterArrayName, const char* jetArrayName, const char* backgroundJetArrayName); - - // Standard functions + // ######### CONTRUCTORS/DESTRUCTORS AND STD FUNCTIONS + AliAnalysisTaskChargedJetsPA(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName); + AliAnalysisTaskChargedJetsPA(); virtual ~AliAnalysisTaskChargedJetsPA(); virtual void UserCreateOutputObjects(); virtual void UserExec(Option_t *option); virtual void Terminate(Option_t *); - // Setters - void SetAnalyzeTracks(Bool_t val) {fAnalyzeQA = val;} - void SetAnalyzeJets(Bool_t val) {fAnalyzeJets = val;} - void SetAnalyzeBackground(Bool_t val) {fAnalyzeBackground = val;} - void SetAnalyzePythia(Bool_t val) {fAnalyzePythia = val;} - - void SetTrackMinPt(Double_t minPt) {fMinJetPt = minPt;} - void SetSignalJetMinPt(Double_t minPt) {fMinJetPt = minPt;} - void SetSignalJetMinArea(Double_t minArea) {fMinJetArea = minArea;} - void SetBackgroundJetMinPt(Double_t minPt) {fMinBackgroundJetPt = minPt;} - void SetDijetLeadingMinPt(Double_t minPt) {fMinDijetLeadingPt = minPt;} - void SetNumberOfPtHardBins(Int_t count) {fNumPtHardBins = count;} - - void SetNumberOfRandConesPerEvent(Int_t count) {fNumberRandCones = count;} - void SetRandConeRadius(Double_t radius) {fRandConeRadius = radius;} - void SetSignalJetRadius(Double_t radius) {fSignalJetRadius = radius;} - void SetBackgroundJetRadius(Double_t radius) {fBackgroundJetRadius = radius;} - void SetKTDeltaPtEtaBin(Int_t bin) {fKTDeltaPtEtaBin = bin;} - void SetTrackBackgroundConeRadius(Double_t radius) {fTrackBackgroundConeRadius = radius;} - - void SetDijetMaxAngleDeviation(Double_t degrees) {fDijetMaxAngleDeviation = degrees/360.0 * TMath::TwoPi();} // degrees are more comfortable - void SetAcceptanceWindows(Double_t trackEta, Double_t vertexZ, Double_t vertexMaxR, Double_t signalJetRadius, Double_t bgrdJetRadius) - { - fVertexWindow = vertexZ; - fVertexMaxR = vertexMaxR; - fTrackEtaWindow = trackEta; - fSignalJetRadius = signalJetRadius; - fBackgroundJetRadius = bgrdJetRadius; - fSignalJetEtaWindow = fTrackEtaWindow-fSignalJetRadius; - - fBackgroundJetEtaWindow = fTrackEtaWindow-fBackgroundJetRadius; - } - - void SetCorrectionFactors(TH2D* histo) - { - // COPY given histogram - fJetBgrdCorrectionFactors = new TH2D(*histo); - - if (!fJetBgrdCorrectionFactors) - AliError(Form("Setting the correction factors with %s (%s) failed! You won't get eta-corrected spectra!", histo->GetName(), histo->IsA()->GetName())); - - // Look, if given histogram is compatible with given code - if (fJetBgrdCorrectionFactors->GetXaxis()->GetNbins() != 5) - AliError(Form("Setting the correction factors failed, because the given histogram is not compatible! You need nbinX=5 (currently:%d)",fJetBgrdCorrectionFactors->GetXaxis()->GetNbins())); - } - - // Getters - Int_t GetInstanceCounter() {return fTaskInstanceCounter;} + // ######### SETTERS/GETTERS + void SetAnalyzeJets(Bool_t val) {fAnalyzeJets = val;} + void SetAnalyzeBackground(Bool_t val) {fAnalyzeBackground = val;} + void SetAnalyzePythia(Bool_t val) {fAnalyzePythia = val;} + + void SetTrackMinPt(Double_t minPt) {fMinJetPt = minPt;} + void SetSignalJetMinPt(Double_t minPt) {fMinJetPt = minPt;} + void SetSignalJetMinArea(Double_t minArea) {fMinJetArea = minArea;} + void SetBackgroundJetMinPt(Double_t minPt) {fMinBackgroundJetPt = minPt;} + void SetDijetLeadingMinPt(Double_t minPt) {fMinDijetLeadingPt = minPt;} + void SetNumberOfPtHardBins(Int_t count) {fNumPtHardBins = count;} + void SetNumberOfRandConesPerEvent(Int_t count) {fNumberRandCones = count;} + void SetRandConeRadius(Double_t radius) {fRandConeRadius = radius;} + void SetSignalJetRadius(Double_t radius) {fSignalJetRadius = radius;} + void SetBackgroundJetRadius(Double_t radius) {fBackgroundJetRadius = radius;} + void SetTRBackgroundConeRadius(Double_t radius) {fTRBackgroundConeRadius = radius;} + void SetCentralityType(const char* type) {fCentralityType = type;} + + void SetDijetMaxAngleDeviation(Double_t degrees) {fDijetMaxAngleDeviation = degrees/360.0 * TMath::TwoPi();} // degrees are more comfortable + void SetAcceptanceWindows(Double_t trackEta, Double_t vertexZ, Double_t vertexMaxR, Double_t signalJetRadius, Double_t bgrdJetRadius){fVertexWindow = vertexZ; fVertexMaxR = vertexMaxR; fTrackEtaWindow = trackEta; fSignalJetRadius = signalJetRadius; fBackgroundJetRadius = bgrdJetRadius; fSignalJetEtaWindow = fTrackEtaWindow-fSignalJetRadius; fBackgroundJetEtaWindow = fTrackEtaWindow-fBackgroundJetRadius;} + void SetKTEtaCorrectionFactors(TH2D* histo); + void SetRCEtaCorrectionFactors(TH2D* histo); + void SetTREtaCorrectionFactors(TH2D* histo); + Int_t GetInstanceCounter() {return fTaskInstanceCounter;} private: - // Calculation functions - void GetSignalJets(); - Int_t GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets); - Double_t GetJetBackgroundCorrFactor(Double_t eta, Double_t background); - Double_t GetCorrectedJetPt(AliEmcalJet* jet, Double_t background, Bool_t useEtaCorrection = kFALSE); - void GetDeltaPt(Double_t& deltaPt, Double_t rho, Int_t numberExcludeLeadingJets = 0, Int_t usedEtaBin = -1, Bool_t useEtaCorrection = kFALSE); - - void GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMedian, Double_t& areaMean, Double_t etaMin = 0, Double_t etaMax = 0); - void GetKTBackground2Density(Int_t numberExcludeLeadingJets, Double_t& rhoMedian, Double_t& areaMean, Double_t etaMin = 0, Double_t etaMax = 0); - Int_t GetRCBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& rhoMedian, Double_t etaMin = 0, Double_t etaMax = 0, Int_t numberRandCones = 0); - void GetTrackBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, Double_t etaMin = 0, Double_t etaMax = 0); - void GetTrackBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular); - Double_t GetConePt(Double_t eta, Double_t phi, Double_t radius); - Double_t GetPtHard(); - Int_t GetPtHardBin(); - void GetPerpendicularCone(Double_t vecPhi, Double_t vecTheta, Double_t& conePt); - - // Cut checks - Bool_t IsTrackInAcceptance(AliVParticle* track); - Bool_t IsClusterInAcceptance(AliVCluster* cluster); - Bool_t IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius); - - Bool_t IsBackgroundJetInAcceptance(AliEmcalJet* jet); - Bool_t IsSignalJetInAcceptance(AliEmcalJet* jet); - Bool_t IsDijet(AliEmcalJet* jet1, AliEmcalJet* jet2); + enum EtaCorrectionMode {kNoEtaCorrection, kKTEtaCorrection, kRCEtaCorrection, kTREtaCorrection}; + + // ######### MAIN CALCULATION FUNCTIONS + void GetSignalJets(); + Int_t GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray, Bool_t isSignalJets); + Double_t GetBackgroundEtaCorrFactor(EtaCorrectionMode mode, Double_t eta, Double_t background); + Double_t GetBackgroundEtaBinCorrFactor(EtaCorrectionMode mode, Int_t eta, Double_t background); + Double_t GetCorrectedJetPt(AliEmcalJet* jet, Double_t background, EtaCorrectionMode mode); + void GetDeltaPt(Double_t& deltaPt, Double_t rho, EtaCorrectionMode mode, Bool_t leadingJetExclusion = kTRUE); + + void GetKTBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMedian, Double_t& areaMean, Double_t etaMin = 0, Double_t etaMax = 0); + Int_t GetRCBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& rhoMedian, Double_t etaMin = 0, Double_t etaMax = 0, Int_t numberRandCones = 0); + void GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, Double_t etaMin = 0, Double_t etaMax = 0); + void GetTRBackgroundDensity(Int_t numberExcludeLeadingJets, Double_t& rhoMean, Double_t& area, AliEmcalJet* excludeJet1, AliEmcalJet* excludeJet2, Bool_t doSearchPerpendicular); + Double_t GetConePt(Double_t eta, Double_t phi, Double_t radius); + Double_t GetPtHard(); + Int_t GetPtHardBin(); + void GetPerpendicularCone(Double_t vecPhi, Double_t vecTheta, Double_t& conePt); + + // ######### CHECK FUNCTIONS + Bool_t IsTrackInAcceptance(AliVParticle* track); + Bool_t IsTrackInCone(AliVTrack* track, Double_t eta, Double_t phi, Double_t radius); + + Bool_t IsBackgroundJetInAcceptance(AliEmcalJet* jet); + Bool_t IsSignalJetInAcceptance(AliEmcalJet* jet); + Bool_t IsDijet(AliEmcalJet* jet1, AliEmcalJet* jet2); - // Some helpers - Double_t EtaToTheta(Double_t arg){return 2.*atan(exp(-arg));} - Double_t ThetaToEta(Double_t arg) - { - if ((arg > TMath::Pi()) || (arg < 0.0)) - { - AliError(Form("ThetaToEta got wrong input! (%f)", arg)); - return 0.0; - } - return -log(tan(arg/2.)); - } - Double_t GetDeltaPhi(Double_t phi1, Double_t phi2) {return min(TMath::Abs(phi1-phi2),TMath::TwoPi()- TMath::Abs(phi1-phi2));} - - // #### This functions return the ratio of a rectangle that is covered by a circle - Double_t MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax) - { - const Int_t kTests = 1000; - Int_t hits = 0; - TRandom3 randomGen(0); - - // Loop over kTests-many tests - for (Int_t i=0; i(hits)/static_cast(kTests)); - } - - void FillHistogram(const char * key, Double_t x); - void FillHistogram(const char * key, Double_t x, Double_t y); - void FillHistogram(const char * key, Double_t x, Double_t y, Double_t add); + // ######### HELPER FUNCTIONS + Double_t EtaToTheta(Double_t arg); + Double_t ThetaToEta(Double_t arg); + Double_t GetDeltaPhi(Double_t phi1, Double_t phi2); + Double_t MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax); + + // ######### HISTOGRAM FUNCTIONS + void FillHistogram(const char * key, Double_t x); + void FillHistogram(const char * key, Double_t x, Double_t y); + void FillHistogram(const char * key, Double_t x, Double_t y, Double_t add); const char* GetHistoName(const char* name) { if (fIsMC) return Form("H%d_%s_MC", fTaskInstanceCounter, name); - return Form("H%d_%s", fTaskInstanceCounter, name); } - template T* AddHistogram1D(const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, const char* xTitle = "x axis", const char* yTitle = "y axis"); - template T* AddHistogram2D(const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, Int_t yBins = 100, Double_t yMin = 0.0, Double_t yMax = 20.0, const char* xTitle = "x axis", const char* yTitle = "y axis", const char* zTitle = "z axis"); - - - // Standard functions + // ######### STANDARD FUNCTIONS Bool_t Notify(); void Calculate(AliVEvent* event); void ExecOnce(); void Init (); TList* fOutputList; //! Output list - // ########## TRIGGERS - Bool_t fAnalyzeQA; // trigger if tracks should be analyzed + // ########## USAGE TRIGGERS Bool_t fAnalyzeJets; // trigger if jets should be processed Bool_t fAnalyzeBackground; // trigger if background should be processed Bool_t fAnalyzePythia; // trigger if pythia properties should be processed Bool_t fHasTracks; // trigger if tracks are actually valid - Bool_t fHasClusters; // trigger if clusters are actually valid Bool_t fHasJets; // trigger if jets are actually valid Bool_t fHasBackgroundJets; // trigger if background is actually valid Bool_t fIsMC; // trigger if data is MC (for naming reasons) + // ########## SOURCE INFORMATION TClonesArray* fJetArray; //! object containing the jets TClonesArray* fTrackArray; //! object containing the tracks - TClonesArray* fClusterArray; //! object containing the clusters TClonesArray* fBackgroundJetArray; //! object containing background jets TString* fJetArrayName; // name of object containing the jets TString* fTrackArrayName; // name of object containing the tracks - TString* fClusterArrayName; // name of object containing the tracks TString* fBackgroundJetArrayName;// name of object containing event wise bckgrds Int_t fNumPtHardBins; // Number of used pt hard bins @@ -231,27 +125,26 @@ class AliAnalysisTaskChargedJetsPA : public AliAnalysisTaskSE { Double_t fRandConeRadius; // Radius for the random cones Double_t fSignalJetRadius; // Radius for the signal jets Double_t fBackgroundJetRadius; // Radius for the KT background jets - Int_t fKTDeltaPtEtaBin; // Bin, in which the KT delta pt is calculate in case of eta correction (-1) - Double_t fTrackBackgroundConeRadius; // Radius for the jets excluded in track background + Double_t fTRBackgroundConeRadius;// Radius for the jets excluded in track background Int_t fNumberRandCones; // Number of random cones to be put into one event Int_t fNumberExcludedJets; // Number of jets to be excluded from backgrounds Double_t fDijetMaxAngleDeviation;// Max angle deviation from pi between two jets to be accept. as dijet - Int_t fBackgroundEtaBins; // Number of eta bins for the RC/track background - TH2D* fJetBgrdCorrectionFactors;// Correction factors in bins of rho and eta to correct the eta dependence of the jet background + TH2D* fJetKTEtaCorrection; // Correction factors in bins of rho and eta to correct the eta dependence of the jet background + TH2D* fJetRCEtaCorrection; // Correction factors in bins of rho and eta to correct the eta dependence of the jet background + TH2D* fJetTREtaCorrection; // Correction factors in bins of rho and eta to correct the eta dependence of the jet background // ########## CUTS Double_t fSignalJetEtaWindow; // +- window in eta for signal jets Double_t fBackgroundJetEtaWindow;// +- window in eta for background jets Double_t fTrackEtaWindow; // +- window in eta for tracks - Double_t fClusterEtaWindow; // +- window in eta for clusters Double_t fVertexWindow; // +- window in Z for the vertex Double_t fVertexMaxR; // +- window in R for the vertex (distance in xy-plane) Double_t fMinTrackPt; // Min track pt to be accepted - Double_t fMinClusterPt; // Min track pt to be accepted Double_t fMinJetPt; // Min jet pt to be accepted Double_t fMinJetArea; // Min jet area to be accepted Double_t fMinBackgroundJetPt; // Min jet pt to be accepted as background jet Double_t fMinDijetLeadingPt; // Min jet pt to be accepted as constituent of dijet + TString fCentralityType; // Used centrality estimate (V0A, V0C, V0M, ...) // ########## EVENT PROPERTIES AliEmcalJet* fFirstLeadingJet; //! leading jet in event @@ -269,10 +162,10 @@ class AliAnalysisTaskChargedJetsPA : public AliAnalysisTaskSE { TList* fHistList; // Histogram list Int_t fHistCount; // Histogram count - AliAnalysisTaskChargedJetsPA(const AliAnalysisTaskChargedJetsPA&); // not implemented - AliAnalysisTaskChargedJetsPA& operator=(const AliAnalysisTaskChargedJetsPA&); // not implemented + AliAnalysisTaskChargedJetsPA(const AliAnalysisTaskChargedJetsPA&); + AliAnalysisTaskChargedJetsPA& operator=(const AliAnalysisTaskChargedJetsPA&); - ClassDef(AliAnalysisTaskChargedJetsPA, 1); // Charged jet analysis for pA + ClassDef(AliAnalysisTaskChargedJetsPA, 2); // Charged jet analysis for pA }; #endif diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.cxx index dc0df0dd5a4..2dabe6229b2 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.cxx +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.cxx @@ -1,3 +1,44 @@ +#ifndef ALIANALYSISTASKSE_H +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliAnalysisTask.h" +#include "AliCentrality.h" +#include "AliStack.h" +#include "AliESDEvent.h" +#include "AliESDInputHandler.h" +#include "AliAODEvent.h" +#include "AliAODHandler.h" +#include "AliAnalysisManager.h" +#include "AliAnalysisTaskSE.h" +#endif + +#include +#include "AliGenPythiaEventHeader.h" +#include "AliMCEvent.h" +#include "AliLog.h" +#include +#include +#include "AliVEventHandler.h" +#include "AliVParticle.h" +#include "AliAnalysisUtils.h" + #include "AliAnalysisTaskQualityAssurancePA.h" @@ -60,7 +101,7 @@ void AliAnalysisTaskQualityAssurancePA::Init() AddHistogram1D(tmpRunNum, "hJetArea", "Jets area distribution", "", 200, 0., 2., "Area","dN^{Jets}/dA"); AddHistogram2D(tmpRunNum, "hJetAreaVsPt", "Jets area vs. p_{T} distribution", "COLZ", 200, 0., 2., 400, 0., 200., "Area", "p_{T}", "dN^{Jets}/dA dp_{T}"); - AddHistogram2D(tmpRunNum, "hJetPhiEta", "Jets angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -0.6, 0.6, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); + AddHistogram2D(tmpRunNum, "hJetPhiEta", "Jets angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -0.6, 0.6, "#phi","#eta","dN^{Jets}/(d#phi d#eta)"); AddHistogram2D(tmpRunNum, "hJetPtVsConstituentCount", "Jets number of constituents vs. jet p_{T}", "COLZ", 800, 0., 400., 100, 0., 100., "p_{T}","N^{Tracks}","dN^{Jets}/(dp_{T} dN^{tracks})"); AddHistogram1D(tmpRunNum, "hJetCountAll", "Number of Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}"); AddHistogram1D(tmpRunNum, "hJetCountAccepted", "Number of accepted Jets", "", 200, 0., 200., "N jets","dN^{Events}/dN^{Jets}"); @@ -77,6 +118,14 @@ void AliAnalysisTaskQualityAssurancePA::Init() } + +//________________________________________________________________________ +AliAnalysisTaskQualityAssurancePA::AliAnalysisTaskQualityAssurancePA() : AliAnalysisTaskSE("AliAnalysisTaskQualityAssurancePA"), fOutputList(0), fAnalyzeQA(1), fAnalyzeJets(1), fAnalyzePythia(0), fHasTracks(0), fHasClusters(0), fHasJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fClusterArray(0), fJetArrayName(0), fTrackArrayName(0), fClusterArrayName(0), fRunNumbers(0), fNumPtHardBins(11), fSignalJetRadius(0.4), fNumberExcludedJets(2), fSignalJetEtaWindow(0.5), fTrackEtaWindow(0.9), fClusterEtaWindow(0.7), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinClusterPt(0.300), fMinJetPt(1.0), fMinJetArea(0.4), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0) +{ + // default constructor +} + + //________________________________________________________________________ AliAnalysisTaskQualityAssurancePA::AliAnalysisTaskQualityAssurancePA(const char *name, const char* trackArrayName, const char* clusterArrayName, const char* jetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeQA(1), fAnalyzeJets(1), fAnalyzePythia(0), fHasTracks(0), fHasClusters(0), fHasJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fClusterArray(0), fJetArrayName(0), fTrackArrayName(0), fClusterArrayName(0), fRunNumbers(0), fNumPtHardBins(11), fSignalJetRadius(0.4), fNumberExcludedJets(2), fSignalJetEtaWindow(0.5), fTrackEtaWindow(0.9), fClusterEtaWindow(0.7), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinClusterPt(0.300), fMinJetPt(1.0), fMinJetArea(0.4), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0) { diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.h index 54f9d88ddbe..97c63ff67b5 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.h +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.h @@ -1,8 +1,7 @@ #ifndef ALIANALYSISTASKQUALITYASSURANCEPA_H #define ALIANALYSISTASKQUALITYASSURANCEPA_H -// #define DEBUGMODE - +//#define DEBUGMODE class TH1F; class TH2F; @@ -15,159 +14,71 @@ class AliVParticle; class AliLog; class AliAnalysisUtils; -#ifndef ALIANALYSISTASKSE_H -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "AliAnalysisTask.h" -#include "AliCentrality.h" -#include "AliStack.h" -#include "AliESDEvent.h" -#include "AliESDInputHandler.h" -#include "AliAODEvent.h" -#include "AliAODHandler.h" -#include "AliAnalysisManager.h" -#include "AliAnalysisTaskSE.h" -#endif - -#include -#include "AliGenPythiaEventHeader.h" -#include "AliMCEvent.h" -#include "AliLog.h" -#include -#include -#include "AliVEventHandler.h" -#include "AliVParticle.h" -#include "AliAnalysisUtils.h" - - -class AliAnalysisTaskQualityAssurancePA : public AliAnalysisTaskSE { +class AliAnalysisTaskQualityAssurancePA : public AliAnalysisTaskSE +{ public: - AliAnalysisTaskQualityAssurancePA() : AliAnalysisTaskSE(), fOutputList(0), fAnalyzeQA(1), fAnalyzeJets(1), fAnalyzePythia(0), fHasTracks(0), fHasClusters(0), fHasJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fClusterArray(0), fJetArrayName(0), fTrackArrayName(0), fClusterArrayName(0), fRunNumbers(0), fNumPtHardBins(11), fSignalJetRadius(0.4), fNumberExcludedJets(2), fSignalJetEtaWindow(0.5), fTrackEtaWindow(0.9), fClusterEtaWindow(0.7), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinClusterPt(0.300), fMinJetPt(1.0), fMinJetArea(0.4), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0) {} - AliAnalysisTaskQualityAssurancePA(const char *name, const char* trackArrayName, const char* clusterArrayName, const char* jetArrayName); + AliAnalysisTaskQualityAssurancePA(); // Standard functions - virtual ~AliAnalysisTaskQualityAssurancePA(); - virtual void UserCreateOutputObjects(); - virtual void UserExec(Option_t *option); - virtual void Terminate(Option_t *); + virtual ~AliAnalysisTaskQualityAssurancePA(); + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *); // Setters - void SetAnalyzeTracks(Bool_t val) {fAnalyzeQA = val;} - void SetAnalyzeJets(Bool_t val) {fAnalyzeJets = val;} - void SetAnalyzePythia(Bool_t val) {fAnalyzePythia = val;} - - void SetTrackMinPt(Double_t minPt) {fMinJetPt = minPt;} - void SetSignalJetMinPt(Double_t minPt) {fMinJetPt = minPt;} - void SetSignalJetMinArea(Double_t minArea) {fMinJetArea = minArea;} - void SetRunNumbers(const char* runNumbers) {*fRunNumbers = runNumbers;} - void SetNumberOfPtHardBins(Int_t count) {fNumPtHardBins = count;} - - void SetSignalJetRadius(Double_t radius) {fSignalJetRadius = radius;} - void SetAcceptanceWindows(Double_t trackEta, Double_t vertexZ, Double_t vertexMaxR, Double_t signalJetRadius) - { - fVertexWindow = vertexZ; - fVertexMaxR = vertexMaxR; - fTrackEtaWindow = trackEta; - fSignalJetRadius = signalJetRadius; - fSignalJetEtaWindow = fTrackEtaWindow-fSignalJetRadius; - } + void SetAnalyzeTracks(Bool_t val) {fAnalyzeQA = val;} + void SetAnalyzeJets(Bool_t val) {fAnalyzeJets = val;} + void SetAnalyzePythia(Bool_t val) {fAnalyzePythia = val;} - // Getters - Int_t GetInstanceCounter() {return fTaskInstanceCounter;} + void SetTrackMinPt(Double_t minPt) {fMinJetPt = minPt;} + void SetSignalJetMinPt(Double_t minPt) {fMinJetPt = minPt;} + void SetSignalJetMinArea(Double_t minArea) {fMinJetArea = minArea;} + void SetRunNumbers(const char* runNumbers) {*fRunNumbers = runNumbers;} + void SetNumberOfPtHardBins(Int_t count) {fNumPtHardBins = count;} - private: - // Calculation functions - void GetSignalJets(); - Int_t GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray); - Double_t GetPtHard(); - Int_t GetPtHardBin(); - - // Cut checks - Bool_t IsTrackInAcceptance(AliVParticle* track); - Bool_t IsClusterInAcceptance(AliVCluster* cluster); - Bool_t IsSignalJetInAcceptance(AliEmcalJet* jet); - - // Some helpers - Double_t EtaToTheta(Double_t arg){return 2.*atan(exp(-arg));} - Double_t ThetaToEta(Double_t arg) - { - if ((arg > TMath::Pi()) || (arg < 0.0)) - { - AliError(Form("ThetaToEta got wrong input! (%f)", arg)); - return 0.0; - } - return -log(tan(arg/2.)); - } - Double_t GetDeltaPhi(Double_t phi1, Double_t phi2) {return min(TMath::Abs(phi1-phi2),TMath::TwoPi()- TMath::Abs(phi1-phi2));} + void SetSignalJetRadius(Double_t radius) {fSignalJetRadius = radius;} + void SetAcceptanceWindows(Double_t trackEta, Double_t vertexZ, Double_t vertexMaxR, Double_t signalJetRadius) {fVertexWindow = vertexZ; fVertexMaxR = vertexMaxR; fTrackEtaWindow = trackEta; fSignalJetRadius = signalJetRadius; fSignalJetEtaWindow = fTrackEtaWindow-fSignalJetRadius;} - // #### This functions return the ratio of a rectangle that is covered by a circle - Double_t MCGetOverlapCircleRectancle(Double_t cPosX, Double_t cPosY, Double_t cRadius, Double_t rPosXmin, Double_t rPosXmax, Double_t rPosYmin, Double_t rPosYmax) - { - const Int_t kTests = 1000; - Int_t hits = 0; - TRandom3 randomGen(0); - - // Loop over kTests-many tests - for (Int_t i=0; i(hits)/static_cast(kTests)); - } + // Getters + Int_t GetInstanceCounter() {return fTaskInstanceCounter;} - void FillHistogram(const char* runNumber, const char * key, Double_t x); - void FillHistogram(const char* runNumber, const char * key, Double_t x, Double_t y); - void FillHistogram(const char* runNumber, const char * key, Double_t x, Double_t y, Double_t add); - const char* GetHistoName(const char* runNumber, const char* name) + private: + // ######### MAIN CALCULATION FUNCTIONS + void GetSignalJets(); + Int_t GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray); + Double_t GetPtHard(); + Int_t GetPtHardBin(); + + // ######### CHECK FUNCTIONS + Bool_t IsTrackInAcceptance(AliVParticle* track); + Bool_t IsClusterInAcceptance(AliVCluster* cluster); + Bool_t IsSignalJetInAcceptance(AliEmcalJet* jet); + + // ######### HISTOGRAM FUNCTIONS + void FillHistogram(const char* runNumber, const char * key, Double_t x); + void FillHistogram(const char* runNumber, const char * key, Double_t x, Double_t y); + void FillHistogram(const char* runNumber, const char * key, Double_t x, Double_t y, Double_t add); + const char* GetHistoName(const char* runNumber, const char* name) { if (fIsMC) return Form("H%d_%s_%s_MC", fTaskInstanceCounter, name, runNumber); - return Form("H%d_%s_%s", fTaskInstanceCounter, name, runNumber); } - template T* AddHistogram1D(const char* runNumber, const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, const char* xTitle = "x axis", const char* yTitle = "y axis"); - template T* AddHistogram2D(const char* runNumber, const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, Int_t yBins = 100, Double_t yMin = 0.0, Double_t yMax = 20.0, const char* xTitle = "x axis", const char* yTitle = "y axis", const char* zTitle = "z axis"); - - // Standard functions - Bool_t Notify(); - void Calculate(AliVEvent* event); - void ExecOnce(); - void Init (); + // ######### STANDARD FUNCTIONS + Bool_t Notify(); + void Calculate(AliVEvent* event); + void ExecOnce(); + void Init (); TList* fOutputList; //! Output list - // ########## TRIGGERS + + // ########## USAGE TRIGGERS Bool_t fAnalyzeQA; // trigger if tracks should be analyzed Bool_t fAnalyzeJets; // trigger if jets should be processed Bool_t fAnalyzePythia; // trigger if pythia properties should be processed @@ -175,6 +86,7 @@ class AliAnalysisTaskQualityAssurancePA : public AliAnalysisTaskSE { Bool_t fHasClusters; // trigger if clusters are actually valid Bool_t fHasJets; // trigger if jets are actually valid Bool_t fIsMC; // trigger if data is MC (for naming reasons) + // ########## SOURCE INFORMATION TClonesArray* fJetArray; //! object containing the jets TClonesArray* fTrackArray; //! object containing the tracks @@ -210,7 +122,7 @@ class AliAnalysisTaskQualityAssurancePA : public AliAnalysisTaskSE { // ########## GENERAL VARS TRandom3* fRandom; //! A random number - AliAnalysisUtils* fHelperClass; //! Vertex selection helper + AliAnalysisUtils* fHelperClass; //! Vertex selection helper Bool_t fInitialized; //! trigger if tracks/jets are loaded Int_t fTaskInstanceCounter; // for naming reasons TList* fHistList; // Histogram list @@ -219,7 +131,7 @@ class AliAnalysisTaskQualityAssurancePA : public AliAnalysisTaskSE { AliAnalysisTaskQualityAssurancePA(const AliAnalysisTaskQualityAssurancePA&); // not implemented AliAnalysisTaskQualityAssurancePA& operator=(const AliAnalysisTaskQualityAssurancePA&); // not implemented - ClassDef(AliAnalysisTaskQualityAssurancePA, 1); // QA helper task for pA + ClassDef(AliAnalysisTaskQualityAssurancePA, 2); // QA helper task for pA }; #endif diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskChargedJetsPA.C b/PWGJE/EMCALJetTasks/macros/AddTaskChargedJetsPA.C index 0f3e5beb0d5..5bda573c74a 100644 --- a/PWGJE/EMCALJetTasks/macros/AddTaskChargedJetsPA.C +++ b/PWGJE/EMCALJetTasks/macros/AddTaskChargedJetsPA.C @@ -5,14 +5,16 @@ AliAnalysisTaskChargedJetsPA* AddTaskChargedJetsPA( Double_t randomConeR = 0.4, Double_t trackBgrdConeR = 0.6, const char* usedTracks = "PicoTracks", - const char* usedClusters = "CaloClustersCorr", + const char* centralityType = "V0A", Double_t trackEtaWindow = 0.9, Double_t vertexWindow = 10.0, Double_t vertexMaxR = 1.0, Double_t minJetPt = 5.0, // signal jet min pt Double_t dijetLeadingMinPt = 10.0, Double_t dijetMaxAngleDev = 10.0, - Int_t numberOfPtHardBins = 0 + Int_t numberOfPtHardBins = 0, + const char* fileEtaCorrectionFactors= "alien:///alice/cern.ch/user/r/rhaake/pA/EtaCorrectionFactors.root", + const char* externalMacro = NULL ) { // #### Detect the demanded trigger with its readable name @@ -49,10 +51,35 @@ AliAnalysisTaskChargedJetsPA* AddTaskChargedJetsPA( AliEmcalJetTask* jetFinderTask = AddTaskEmcalJet(usedTracks,"",1,jetRadius,1,0.150,0.300);// anti-kt AliEmcalJetTask* jetFinderTaskKT = AddTaskEmcalJet(usedTracks,"",0,jetRadius,1,0.150,0.300); // kt + // #### Load correction factors from alien + TH2D* corrFactorsKT = NULL; + TH2D* corrFactorsRC = NULL; + TH2D* corrFactorsTR = NULL; + if (fileEtaCorrectionFactors) + { + // trying to connect to alien + if (!TGrid::Connect("alien://")) + ::Warning("AddTaskChargedJetsPA", "AliEn connection failed!"); + else + { + ::Info("AddTaskChargedJetsPA", "AliEn connection successful!"); + // Copy eta correction file + Bool_t copied = TFile::Cp(fileEtaCorrectionFactors,"file:EtaCorrectionFactors.root"); + if(copied) + { + TFile* tmpFile= new TFile("EtaCorrectionFactors.root","READ"); + corrFactorsKT = static_cast(tmpFile->Get("EtaCorrectionFactorsKT")); + corrFactorsRC = static_cast(tmpFile->Get("EtaCorrectionFactorsRC")); + corrFactorsTR = static_cast(tmpFile->Get("EtaCorrectionFactorsTR")); + } + else + ::Warning("AddTaskChargedJetsPA", "AliEn copying failed!"); + } + } // #### Define analysis task AliAnalysisTaskChargedJetsPA *task = NULL; contHistos = manager->CreateContainer(myContName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:ChargedJetsPA", AliAnalysisManager::GetCommonFileName())); - task = new AliAnalysisTaskChargedJetsPA(Form("AnalysisPA_%s_%s", jetFinderTask->GetName(), triggerName.Data()), usedTracks, usedClusters, jetFinderTask->GetName(),jetFinderTaskKT->GetName()); + task = new AliAnalysisTaskChargedJetsPA(Form("AnalysisPA_%s_%s", jetFinderTask->GetName(), triggerName.Data()), usedTracks, jetFinderTask->GetName(),jetFinderTaskKT->GetName()); // #### Task preferences task->SetAcceptanceWindows(trackEtaWindow, vertexWindow, vertexMaxR, jetRadius, jetRadius); @@ -61,14 +88,28 @@ AliAnalysisTaskChargedJetsPA* AddTaskChargedJetsPA( task->SetDijetLeadingMinPt(dijetLeadingMinPt); task->SetDijetMaxAngleDeviation(dijetMaxAngleDev); task->SetRandConeRadius(randomConeR); - task->SetTrackBackgroundConeRadius(trackBgrdConeR); + task->SetTRBackgroundConeRadius(trackBgrdConeR); task->SelectCollisionCandidates(trigger); + task->SetCentralityType(centralityType); if(numberOfPtHardBins) task->SetNumberOfPtHardBins(numberOfPtHardBins); + if(corrFactorsKT) + task->SetKTEtaCorrectionFactors(corrFactorsKT); + if(corrFactorsRC) + task->SetRCEtaCorrectionFactors(corrFactorsRC); + if(corrFactorsTR) + task->SetTREtaCorrectionFactors(corrFactorsTR); + // #### Add analysis task manager->AddTask(task); manager->ConnectInput(task, 0, manager->GetCommonInputContainer()); manager->ConnectOutput(task, 1, contHistos); + + // #### Do some nasty piggybacking on demand + if (externalMacro) + gROOT->LoadMacro(externalMacro); + + return task; } diff --git a/PWGJE/EMCALJetTasks/macros/runEMCalJetAnalysis.C b/PWGJE/EMCALJetTasks/macros/runEMCalJetAnalysis.C index ef82cfe859e..78cf1c07cba 100644 --- a/PWGJE/EMCALJetTasks/macros/runEMCalJetAnalysis.C +++ b/PWGJE/EMCALJetTasks/macros/runEMCalJetAnalysis.C @@ -110,7 +110,7 @@ void runEMCalJetAnalysis( cout << "Using " << localFiles.Data() << " as input file list.\n"; // Create MC handler, if MC is demanded - if (isMC) + if (isMC && (usedData != "AOD")) { AliMCEventHandler* mcH = new AliMCEventHandler(); mcH->SetPreReadMode(AliMCEventHandler::kLmPreRead);