X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FEMCALJetTasks%2FAliJetResponseMaker.cxx;h=b8731048df988fe2dc7cac64baeec3d584b56f62;hb=21b312002d48f0af6f47b2ace15b0d2eec941018;hp=aafdd5e4dfdfbc7d033460f14d939ae4be5aabbe;hpb=b16bb00196553a8ad7361d75b928d5acd27f7e4a;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx index aafdd5e4dfd..b8731048df9 100644 --- a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx +++ b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx @@ -7,125 +7,182 @@ #include "AliJetResponseMaker.h" #include -#include #include -#include +#include #include +#include "AliAnalysisManager.h" #include "AliVCluster.h" #include "AliVTrack.h" #include "AliEmcalJet.h" -#include "AliGenPythiaEventHeader.h" #include "AliLog.h" -#include "AliMCEvent.h" +#include "AliRhoParameter.h" +#include "AliNamedArrayI.h" +#include "AliJetContainer.h" +#include "AliParticleContainer.h" +#include "AliClusterContainer.h" ClassImp(AliJetResponseMaker) //________________________________________________________________________ AliJetResponseMaker::AliJetResponseMaker() : AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE), - fMCTracksName("MCParticles"), - fMCJetsName("MCJets"), - fMaxDistance(0.25), - fDoWeighting(kFALSE), - fEventWeightHist(kFALSE), - fMCFiducial(kFALSE), - fMCminEta(0), - fMCmaxEta(0), - fMCminPhi(0), - fMCmaxPhi(0), - fSelectPtHardBin(-999), - fDoMatching(kTRUE), - fPythiaHeader(0), - fEventWeight(0), - fPtHardBin(0), - fNTrials(0), - fMCTracks(0), - fMCJets(0), - fHistNTrials(0), - fHistEvents(0), - fHistMCJetsPhiEta(0), - fHistMCJetsPtArea(0), - fHistMCJetsPhiEtaFiducial(0), - fHistMCJetsPtAreaFiducial(0), - fHistMCJetsNEFvsPt(0), - fHistMCJetsZvsPt(0), - fHistJetsPhiEta(0), - fHistJetsPtArea(0), - fHistJetsCorrPtArea(0), - fHistJetsNEFvsPt(0), - fHistJetsZvsPt(0), - fHistMatchingLevelMCPt(0), - fHistClosestDeltaEtaPhiMCPt(0), - fHistClosestDeltaPtMCPt(0), - fHistClosestDeltaCorrPtMCPt(0), - fHistNonMatchedMCJetsPtArea(0), - fHistNonMatchedJetsPtArea(0), - fHistNonMatchedJetsCorrPtArea(0), - fHistPartvsDetecPt(0), - fHistPartvsDetecCorrPt(0), - fHistMissedMCJetsPtArea(0) + fMatching(kNoMatching), + fMatchingPar1(0), + fMatchingPar2(0), + fUseCellsToMatch(kFALSE), + fMinJetMCPt(1), + fHistoType(0), + fDeltaPtAxis(0), + fDeltaEtaDeltaPhiAxis(0), + fNEFAxis(0), + fZAxis(0), + fIsJet1Rho(kFALSE), + fIsJet2Rho(kFALSE), + fHistRejectionReason1(0), + fHistRejectionReason2(0), + fHistJets1(0), + fHistJets2(0), + fHistMatching(0), + fHistJets1PhiEta(0), + fHistJets1PtArea(0), + fHistJets1CorrPtArea(0), + fHistJets1NEFvsPt(0), + fHistJets1CEFvsCEFPt(0), + fHistJets1ZvsPt(0), + fHistJets2PhiEta(0), + fHistJets2PtArea(0), + fHistJets2CorrPtArea(0), + fHistJets2NEFvsPt(0), + fHistJets2CEFvsCEFPt(0), + fHistJets2ZvsPt(0), + fHistCommonEnergy1vsJet1Pt(0), + fHistCommonEnergy2vsJet2Pt(0), + fHistDistancevsJet1Pt(0), + fHistDistancevsJet2Pt(0), + fHistDistancevsCommonEnergy1(0), + fHistDistancevsCommonEnergy2(0), + fHistCommonEnergy1vsCommonEnergy2(0), + fHistDeltaEtaDeltaPhi(0), + fHistJet2PtOverJet1PtvsJet2Pt(0), + fHistJet1PtOverJet2PtvsJet1Pt(0), + fHistDeltaPtvsJet1Pt(0), + fHistDeltaPtvsJet2Pt(0), + fHistDeltaPtOverJet1PtvsJet1Pt(0), + fHistDeltaPtOverJet2PtvsJet2Pt(0), + fHistDeltaPtvsDistance(0), + fHistDeltaPtvsCommonEnergy1(0), + fHistDeltaPtvsCommonEnergy2(0), + fHistDeltaPtvsArea1(0), + fHistDeltaPtvsArea2(0), + fHistDeltaPtvsDeltaArea(0), + fHistJet1PtvsJet2Pt(0), + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0), + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0), + fHistDeltaCorrPtvsJet1CorrPt(0), + fHistDeltaCorrPtvsJet2CorrPt(0), + fHistDeltaCorrPtvsDistance(0), + fHistDeltaCorrPtvsCommonEnergy1(0), + fHistDeltaCorrPtvsCommonEnergy2(0), + fHistDeltaCorrPtvsArea1(0), + fHistDeltaCorrPtvsArea2(0), + fHistDeltaCorrPtvsDeltaArea(0), + fHistJet1CorrPtvsJet2CorrPt(0), + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0), + fHistDeltaMCPtOverJet2PtvsJet2Pt(0), + fHistDeltaMCPtvsJet1MCPt(0), + fHistDeltaMCPtvsJet2Pt(0), + fHistDeltaMCPtvsDistance(0), + fHistDeltaMCPtvsCommonEnergy1(0), + fHistDeltaMCPtvsCommonEnergy2(0), + fHistDeltaMCPtvsArea1(0), + fHistDeltaMCPtvsArea2(0), + fHistDeltaMCPtvsDeltaArea(0), + fHistJet1MCPtvsJet2Pt(0) { // Default constructor. - for (Int_t i = 0; i < 11; i++) { - fHistEventWeight[i] = 0; - } - SetMakeGeneralHistograms(kTRUE); } //________________________________________________________________________ AliJetResponseMaker::AliJetResponseMaker(const char *name) : AliAnalysisTaskEmcalJet(name, kTRUE), - fMCTracksName("MCParticles"), - fMCJetsName("MCJets"), - fMaxDistance(0.25), - fDoWeighting(kFALSE), - fEventWeightHist(kFALSE), - fMCFiducial(kFALSE), - fMCminEta(0), - fMCmaxEta(0), - fMCminPhi(0), - fMCmaxPhi(0), - fSelectPtHardBin(-999), - fDoMatching(kTRUE), - fPythiaHeader(0), - fEventWeight(0), - fPtHardBin(0), - fNTrials(0), - fMCTracks(0), - fMCJets(0), - fHistNTrials(0), - fHistEvents(0), - fHistMCJetsPhiEta(0), - fHistMCJetsPtArea(0), - fHistMCJetsPhiEtaFiducial(0), - fHistMCJetsPtAreaFiducial(0), - fHistMCJetsNEFvsPt(0), - fHistMCJetsZvsPt(0), - fHistJetsPhiEta(0), - fHistJetsPtArea(0), - fHistJetsCorrPtArea(0), - fHistJetsNEFvsPt(0), - fHistJetsZvsPt(0), - fHistMatchingLevelMCPt(0), - fHistClosestDeltaEtaPhiMCPt(0), - fHistClosestDeltaPtMCPt(0), - fHistClosestDeltaCorrPtMCPt(0), - fHistNonMatchedMCJetsPtArea(0), - fHistNonMatchedJetsPtArea(0), - fHistNonMatchedJetsCorrPtArea(0), - fHistPartvsDetecPt(0), - fHistPartvsDetecCorrPt(0), - fHistMissedMCJetsPtArea(0) + fMatching(kNoMatching), + fMatchingPar1(0), + fMatchingPar2(0), + fUseCellsToMatch(kFALSE), + fMinJetMCPt(1), + fHistoType(0), + fDeltaPtAxis(0), + fDeltaEtaDeltaPhiAxis(0), + fNEFAxis(0), + fZAxis(0), + fIsJet1Rho(kFALSE), + fIsJet2Rho(kFALSE), + fHistRejectionReason1(0), + fHistRejectionReason2(0), + fHistJets1(0), + fHistJets2(0), + fHistMatching(0), + fHistJets1PhiEta(0), + fHistJets1PtArea(0), + fHistJets1CorrPtArea(0), + fHistJets1NEFvsPt(0), + fHistJets1CEFvsCEFPt(0), + fHistJets1ZvsPt(0), + fHistJets2PhiEta(0), + fHistJets2PtArea(0), + fHistJets2CorrPtArea(0), + fHistJets2NEFvsPt(0), + fHistJets2CEFvsCEFPt(0), + fHistJets2ZvsPt(0), + fHistCommonEnergy1vsJet1Pt(0), + fHistCommonEnergy2vsJet2Pt(0), + fHistDistancevsJet1Pt(0), + fHistDistancevsJet2Pt(0), + fHistDistancevsCommonEnergy1(0), + fHistDistancevsCommonEnergy2(0), + fHistCommonEnergy1vsCommonEnergy2(0), + fHistDeltaEtaDeltaPhi(0), + fHistJet2PtOverJet1PtvsJet2Pt(0), + fHistJet1PtOverJet2PtvsJet1Pt(0), + fHistDeltaPtvsJet1Pt(0), + fHistDeltaPtvsJet2Pt(0), + fHistDeltaPtOverJet1PtvsJet1Pt(0), + fHistDeltaPtOverJet2PtvsJet2Pt(0), + fHistDeltaPtvsDistance(0), + fHistDeltaPtvsCommonEnergy1(0), + fHistDeltaPtvsCommonEnergy2(0), + fHistDeltaPtvsArea1(0), + fHistDeltaPtvsArea2(0), + fHistDeltaPtvsDeltaArea(0), + fHistJet1PtvsJet2Pt(0), + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0), + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0), + fHistDeltaCorrPtvsJet1CorrPt(0), + fHistDeltaCorrPtvsJet2CorrPt(0), + fHistDeltaCorrPtvsDistance(0), + fHistDeltaCorrPtvsCommonEnergy1(0), + fHistDeltaCorrPtvsCommonEnergy2(0), + fHistDeltaCorrPtvsArea1(0), + fHistDeltaCorrPtvsArea2(0), + fHistDeltaCorrPtvsDeltaArea(0), + fHistJet1CorrPtvsJet2CorrPt(0), + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0), + fHistDeltaMCPtOverJet2PtvsJet2Pt(0), + fHistDeltaMCPtvsJet1MCPt(0), + fHistDeltaMCPtvsJet2Pt(0), + fHistDeltaMCPtvsDistance(0), + fHistDeltaMCPtvsCommonEnergy1(0), + fHistDeltaMCPtvsCommonEnergy2(0), + fHistDeltaMCPtvsArea1(0), + fHistDeltaMCPtvsArea2(0), + fHistDeltaMCPtvsDeltaArea(0), + fHistJet1MCPtvsJet2Pt(0) { // Standard constructor. - for (Int_t i = 0; i < 11; i++) { - fHistEventWeight[i] = 0; - } - SetMakeGeneralHistograms(kTRUE); } @@ -135,263 +192,915 @@ AliJetResponseMaker::~AliJetResponseMaker() // Destructor } + //________________________________________________________________________ -void AliJetResponseMaker::UserCreateOutputObjects() +void AliJetResponseMaker::AllocateTH2() { - // Create user objects. + // Allocate TH2 histograms. - AliAnalysisTaskEmcalJet::UserCreateOutputObjects(); + fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2); + fHistJets1PhiEta->GetXaxis()->SetTitle("#eta"); + fHistJets1PhiEta->GetYaxis()->SetTitle("#phi"); + fOutput->Add(fHistJets1PhiEta); + + fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt); + fHistJets1PtArea->GetXaxis()->SetTitle("area"); + fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)"); + fHistJets1PtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets1PtArea); + + if (fIsJet1Rho) { + fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistJets1CorrPtArea->GetXaxis()->SetTitle("area"); + fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)"); + fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets1CorrPtArea); + } - const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234}; - const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000}; + fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistJets1ZvsPt->GetXaxis()->SetTitle("Z"); + fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)"); + fHistJets1ZvsPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets1ZvsPt); + + fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF"); + fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)"); + fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets1NEFvsPt); + + fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF"); + fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)"); + fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets1CEFvsCEFPt); + + // Jets 2 histograms + + fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2); + fHistJets2PhiEta->GetXaxis()->SetTitle("#eta"); + fHistJets2PhiEta->GetYaxis()->SetTitle("#phi"); + fOutput->Add(fHistJets2PhiEta); + + fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt); + fHistJets2PtArea->GetXaxis()->SetTitle("area"); + fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); + fHistJets2PtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets2PtArea); + + if (fIsJet2Rho) { + fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistJets2CorrPtArea->GetXaxis()->SetTitle("area"); + fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)"); + fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets2CorrPtArea); + } - fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11); - fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin"); - fHistNTrials->GetYaxis()->SetTitle("trials"); - fOutput->Add(fHistNTrials); + fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistJets2ZvsPt->GetXaxis()->SetTitle("Z"); + fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); + fHistJets2ZvsPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets2ZvsPt); + + fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF"); + fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); + fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets2NEFvsPt); + + fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF"); + fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)"); + fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJets2CEFvsCEFPt); + + // Matching histograms + + fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); + fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)"); + fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}"); + fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistCommonEnergy1vsJet1Pt); + + fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); + fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)"); + fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); + fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistCommonEnergy2vsJet2Pt); + + fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); + fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance"); + fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}"); + fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDistancevsJet1Pt); + + fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); + fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance"); + fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); + fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDistancevsJet2Pt); + + fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2); + fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance"); + fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)"); + fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDistancevsCommonEnergy1); + + fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2); + fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance"); + fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)"); + fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDistancevsCommonEnergy2); + + fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2); + fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)"); + fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)"); + fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistCommonEnergy1vsCommonEnergy2); + + fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2); + fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)"); + fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)"); + fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaEtaDeltaPhi); + + fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5); + fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); + fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}"); + fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt); + + fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5); + fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}"); + fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}"); + fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt); + + fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}"); + fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsJet1Pt); + + fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); + fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsJet2Pt); + + fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt", + fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); + fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}"); + fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}"); + fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt); + + fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt", + fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); + fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); + fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}"); + fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt); + + fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance"); + fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsDistance); + + fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)"); + fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsCommonEnergy1); + + fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)"); + fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsCommonEnergy2); + + fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}"); + fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsArea1); + + fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}"); + fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsArea2); + + fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea", + fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}"); + fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaPtvsDeltaArea); + + fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); + fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}"); + fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); + fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJet1PtvsJet2Pt); + + if (fIsJet1Rho || fIsJet2Rho) { + if (!fIsJet1Rho) + fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + else + fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", + 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + + fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}"); + fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt); + + if (!fIsJet2Rho) + fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + else + fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", + 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + + fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}"); + fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt); + + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", + 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5); + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}"); + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}"); + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt); + + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", + 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5); + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}"); + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}"); + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt); + + fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance"); + fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsDistance); + + fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)"); + fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1); + + fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)"); + fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2); + + fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}"); + fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsArea1); + + fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}"); + fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsArea2); + + fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea", + fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}"); + fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); + fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaCorrPtvsDeltaArea); + + if (!fIsJet1Rho) + fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + else if (!fIsJet2Rho) + fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", + 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); + else + fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", + 2*fNbins, -fMaxBinPt, fMaxBinPt, + 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}"); + fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}"); + fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJet1CorrPtvsJet2CorrPt); + } - fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11); - fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin"); - fHistEvents->GetYaxis()->SetTitle("total events"); - fOutput->Add(fHistEvents); + if (fIsEmbedded) { + fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}"); + fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsJet1MCPt); + + fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt", + fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); + fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsJet2Pt); + + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", + fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}"); + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}"); + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt); + + fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt", + fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); + fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); + fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}"); + fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt); + + fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance"); + fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsDistance); + + fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)"); + fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsCommonEnergy1); + + fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2", + fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)"); + fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsCommonEnergy2); + + fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}"); + fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsArea1); + + fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2", + fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}"); + fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsArea2); + + fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea", + fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt); + fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}"); + fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); + fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistDeltaMCPtvsDeltaArea); + + fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt", + fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); + fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}"); + fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); + fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJet1MCPtvsJet2Pt); + } +} - if (fEventWeightHist) { - for (Int_t i = 0; i < 11; i++) { - TString name(Form("fHistEventWeight_%d", i+1)); - fHistEventWeight[i] = new TH1F(name, name, 10, 0, 10); - fOutput->Add(fHistEventWeight[i]); - } +//________________________________________________________________________ +void AliJetResponseMaker::AllocateTHnSparse() +{ + // Allocate THnSparse histograms. + + TString title[25]= {""}; + Int_t nbins[25] = {0}; + Double_t min[25] = {0.}; + Double_t max[25] = {0.}; + Int_t dim = 0; + + title[dim] = "#phi"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1)); + dim++; + + title[dim] = "#eta"; + nbins[dim] = fNbins/4; + min[dim] = -1; + max[dim] = 1; + dim++; + + title[dim] = "p_{T}"; + nbins[dim] = fNbins; + min[dim] = 0; + max[dim] = 250; + dim++; + + title[dim] = "A_{jet}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.5; + dim++; + + if (fNEFAxis) { + title[dim] = "NEF"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; } - for (Int_t i = 1; i < 12; i++) { - fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1])); - fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1])); + if (fZAxis) { + title[dim] = "Z"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; } - fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 20, -2, 2, 32, 0, 6.4); - fHistJetsPhiEta->GetXaxis()->SetTitle("#eta"); - fHistJetsPhiEta->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistJetsPhiEta); - - fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); - fHistJetsPtArea->GetXaxis()->SetTitle("area"); - fHistJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); - fHistJetsPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistJetsPtArea); - - fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt); - fHistJetsCorrPtArea->GetXaxis()->SetTitle("area"); - fHistJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)"); - fHistJetsCorrPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistJetsCorrPtArea); - - fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); - fHistJetsZvsPt->GetXaxis()->SetTitle("Z"); - fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); - fOutput->Add(fHistJetsZvsPt); - - fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); - fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF"); - fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); - fOutput->Add(fHistJetsNEFvsPt); - - fHistMCJetsPhiEta = new TH2F("fHistMCJetsPhiEta", "fHistMCJetsPhiEta", 20, -2, 2, 32, 0, 6.4); - fHistMCJetsPhiEta->GetXaxis()->SetTitle("#eta"); - fHistMCJetsPhiEta->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistMCJetsPhiEta); - - fHistMCJetsPtArea = new TH2F("fHistMCJetsPtArea", "fHistMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); - fHistMCJetsPtArea->GetXaxis()->SetTitle("area"); - fHistMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); - fHistMCJetsPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistMCJetsPtArea); - - fHistMCJetsPhiEtaFiducial = new TH2F("fHistMCJetsPhiEtaFiducial", "fHistMCJetsPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4); - fHistMCJetsPhiEtaFiducial->GetXaxis()->SetTitle("#eta"); - fHistMCJetsPhiEtaFiducial->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistMCJetsPhiEtaFiducial); - - fHistMCJetsPtAreaFiducial = new TH2F("fHistMCJetsPtAreaFiducial", "fHistMCJetsPtAreaFiducial", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); - fHistMCJetsPtAreaFiducial->GetXaxis()->SetTitle("area"); - fHistMCJetsPtAreaFiducial->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); - fHistMCJetsPtAreaFiducial->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistMCJetsPtAreaFiducial); + title[dim] = "p_{T,particle}^{leading} (GeV/c)"; + nbins[dim] = 120; + min[dim] = 0; + max[dim] = 120; + dim++; + + Int_t dim1 = dim, dim2 = dim; + + if (fIsJet1Rho) { + title[dim1] = "p_{T}^{corr}"; + nbins[dim1] = fNbins*2; + min[dim1] = -250; + max[dim1] = 250; + dim1++; + } + + if (fIsEmbedded) { + title[dim1] = "p_{T}^{MC}"; + nbins[dim1] = fNbins; + min[dim1] = 0; + max[dim1] = 250; + dim1++; + } + + fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max); + for (Int_t i = 0; i < dim1; i++) + fHistJets1->GetAxis(i)->SetTitle(title[i]); + fOutput->Add(fHistJets1); + + if (fIsJet2Rho) { + title[dim2] = "p_{T}^{corr}"; + nbins[dim2] = fNbins*2; + min[dim2] = -250; + max[dim2] = 250; + dim2++; + } + + fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max); + for (Int_t i = 0; i < dim2; i++) + fHistJets2->GetAxis(i)->SetTitle(title[i]); + fOutput->Add(fHistJets2); - fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); - fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z"); - fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); - fOutput->Add(fHistMCJetsZvsPt); - - fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); - fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF"); - fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); - fOutput->Add(fHistMCJetsNEFvsPt); - - fHistMatchingLevelMCPt = new TH2F("fHistMatchingLevelMCPt", "fHistMatchingLevelMCPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); - fHistMatchingLevelMCPt->GetXaxis()->SetTitle("Matching level"); - fHistMatchingLevelMCPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); - fOutput->Add(fHistMatchingLevelMCPt); - - fHistClosestDeltaEtaPhiMCPt = new TH3F("fHistClosestDeltaEtaPhiMCPt", "fHistClosestDeltaEtaPhiMCPt", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2, 64, -1.6, 4.8, fNbins, fMinBinPt, fMaxBinPt); - fHistClosestDeltaEtaPhiMCPt->GetXaxis()->SetTitle("#Delta#eta"); - fHistClosestDeltaEtaPhiMCPt->GetYaxis()->SetTitle("#Delta#phi"); - fHistClosestDeltaEtaPhiMCPt->GetZaxis()->SetTitle("p_{T}^{gen}"); - fOutput->Add(fHistClosestDeltaEtaPhiMCPt); - - fHistClosestDeltaPtMCPt = new TH2F("fHistClosestDeltaPtMCPt", "fHistClosestDeltaPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt); - fHistClosestDeltaPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}"); - fHistClosestDeltaPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{rec} (GeV/c)"); - fHistClosestDeltaPtMCPt->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistClosestDeltaPtMCPt); - - fHistClosestDeltaCorrPtMCPt = new TH2F("fHistClosestDeltaCorrPtMCPt", "fHistClosestDeltaCorrPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt); - fHistClosestDeltaCorrPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}"); - fHistClosestDeltaCorrPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)"); - fHistClosestDeltaCorrPtMCPt->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistClosestDeltaCorrPtMCPt); - - fHistNonMatchedMCJetsPtArea = new TH2F("fHistNonMatchedMCJetsPtArea", "fHistNonMatchedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); - fHistNonMatchedMCJetsPtArea->GetXaxis()->SetTitle("area"); - fHistNonMatchedMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); - fHistNonMatchedMCJetsPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistNonMatchedMCJetsPtArea); - - fHistNonMatchedJetsPtArea = new TH2F("fHistNonMatchedJetsPtArea", "fHistNonMatchedJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); - fHistNonMatchedJetsPtArea->GetXaxis()->SetTitle("area"); - fHistNonMatchedJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); - fHistNonMatchedJetsPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistNonMatchedJetsPtArea); - - fHistNonMatchedJetsCorrPtArea = new TH2F("fHistNonMatchedJetsCorrPtArea", "fHistNonMatchedJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt); - fHistNonMatchedJetsCorrPtArea->GetXaxis()->SetTitle("area"); - fHistNonMatchedJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)"); - fHistNonMatchedJetsCorrPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistNonMatchedJetsCorrPtArea); - - fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); - fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}"); - fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}"); - fOutput->Add(fHistPartvsDetecPt); - - fHistPartvsDetecCorrPt = new TH2F("fHistPartvsDetecCorrPt", "fHistPartvsDetecCorrPt", (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); - fHistPartvsDetecCorrPt->GetXaxis()->SetTitle("p_{T}^{corr}"); - fHistPartvsDetecCorrPt->GetYaxis()->SetTitle("p_{T}^{gen}"); - fOutput->Add(fHistPartvsDetecCorrPt); - - fHistMissedMCJetsPtArea = new TH2F("fHistMissedMCJetsPtArea", "fHistMissedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); - fHistMissedMCJetsPtArea->GetXaxis()->SetTitle("area"); - fHistMissedMCJetsPtArea->GetYaxis()->SetTitle("p_{T} (GeV/c)"); - fHistMissedMCJetsPtArea->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistMissedMCJetsPtArea); + // Matching + + dim = 0; + + title[dim] = "p_{T,1}"; + nbins[dim] = fNbins; + min[dim] = 0; + max[dim] = 250; + dim++; + + title[dim] = "p_{T,2}"; + nbins[dim] = fNbins; + min[dim] = 0; + max[dim] = 250; + dim++; + + title[dim] = "A_{jet,1}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.5; + dim++; + + title[dim] = "A_{jet,2}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.5; + dim++; + + title[dim] = "distance"; + nbins[dim] = fNbins/2; + min[dim] = 0; + max[dim] = 1.2; + dim++; + + title[dim] = "CE1"; + nbins[dim] = fNbins/2; + min[dim] = 0; + max[dim] = 1.2; + dim++; + + title[dim] = "CE2"; + nbins[dim] = fNbins/2; + min[dim] = 0; + max[dim] = 1.2; + dim++; + + title[dim] = "p_{T,particle,1}^{leading} (GeV/c)"; + nbins[dim] = 120; + min[dim] = 0; + max[dim] = 120; + dim++; + + title[dim] = "p_{T,particle,2}^{leading} (GeV/c)"; + nbins[dim] = 120; + min[dim] = 0; + max[dim] = 120; + dim++; + + if (fDeltaPtAxis) { + title[dim] = "#deltaA_{jet}"; + nbins[dim] = fNbins/2; + min[dim] = -1; + max[dim] = 1; + dim++; + + title[dim] = "#deltap_{T}"; + nbins[dim] = fNbins*2; + min[dim] = -250; + max[dim] = 250; + dim++; + } + if (fIsJet1Rho) { + title[dim] = "p_{T,1}^{corr}"; + nbins[dim] = fNbins*2; + min[dim] = -250; + max[dim] = 250; + dim++; + } + if (fIsJet2Rho) { + title[dim] = "p_{T,2}^{corr}"; + nbins[dim] = fNbins*2; + min[dim] = -250; + max[dim] = 250; + dim++; + } + if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) { + title[dim] = "#deltap_{T}^{corr}"; + nbins[dim] = fNbins*2; + min[dim] = -250; + max[dim] = 250; + dim++; + } + if (fDeltaEtaDeltaPhiAxis) { + title[dim] = "#delta#eta"; + nbins[dim] = fNbins/2; + min[dim] = -1; + max[dim] = 1; + dim++; + + title[dim] = "#delta#phi"; + nbins[dim] = fNbins/2; + min[dim] = -TMath::Pi()/2; + max[dim] = TMath::Pi()*3/2; + dim++; + } + if (fIsEmbedded) { + title[dim] = "p_{T,1}^{MC}"; + nbins[dim] = fNbins; + min[dim] = 0; + max[dim] = 250; + dim++; + + if (fDeltaPtAxis) { + title[dim] = "#deltap_{T}^{MC}"; + nbins[dim] = fNbins*2; + min[dim] = -250; + max[dim] = 250; + dim++; + } + } - PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram + if (fNEFAxis) { + title[dim] = "NEF_{1}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; + + title[dim] = "NEF_{2}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; + } + + if (fZAxis) { + title[dim] = "Z_{1}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; + + title[dim] = "Z_{2}"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; + } + + fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max); + + for (Int_t i = 0; i < dim; i++) + fHistMatching->GetAxis(i)->SetTitle(title[i]); + + fOutput->Add(fHistMatching); } //________________________________________________________________________ -Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const -{ - // Return true if jet is accepted. +void AliJetResponseMaker::UserCreateOutputObjects() +{ + // Create user objects. - if (jet->Pt() <= fJetPtCut) - return kFALSE; - if (jet->Area() <= fJetAreaCut) - return kFALSE; + AliAnalysisTaskEmcalJet::UserCreateOutputObjects(); - return kTRUE; + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); + + if (!jets1 || !jets2) return; + + if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE; + else fIsJet1Rho = kTRUE; + if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE; + else fIsJet2Rho = kTRUE; + + fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250); + fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason"); + fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)"); + fHistRejectionReason1->GetZaxis()->SetTitle("counts"); + SetRejectionReasonLabels(fHistRejectionReason1->GetXaxis()); + fOutput->Add(fHistRejectionReason1); + + fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250); + fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason"); + fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)"); + fHistRejectionReason2->GetZaxis()->SetTitle("counts"); + SetRejectionReasonLabels(fHistRejectionReason2->GetXaxis()); + fOutput->Add(fHistRejectionReason2); + + if (fHistoType==0) + AllocateTH2(); + else + AllocateTHnSparse(); + + PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram } //________________________________________________________________________ -void AliJetResponseMaker::ExecOnce() +void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t LeadingPt, Double_t CorrPt, Double_t MCPt, Int_t Set) { - // Execute once. - - if (!fMCJetsName.IsNull() && !fMCJets) { - fMCJets = dynamic_cast(InputEvent()->FindListObject(fMCJetsName)); - if (!fMCJets) { - AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data())); - return; - } - else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) { - AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data())); - fMCJets = 0; + if (fHistoType==1) { + THnSparse *histo = 0; + if (Set==1) + histo = fHistJets1; + else if (Set==2) + histo = fHistJets2; + + if (!histo) return; - } - } - if (!fMCTracksName.IsNull() && !fMCTracks) { - fMCTracks = dynamic_cast(InputEvent()->FindListObject(fMCTracksName)); - if (!fMCTracks) { - AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data())); - return; + Double_t contents[20]={0}; + + for (Int_t i = 0; i < histo->GetNdimensions(); i++) { + TString title(histo->GetAxis(i)->GetTitle()); + if (title=="#phi") + contents[i] = Phi; + else if (title=="#eta") + contents[i] = Eta; + else if (title=="p_{T}") + contents[i] = Pt; + else if (title=="A_{jet}") + contents[i] = A; + else if (title=="NEF") + contents[i] = NEF; + else if (title=="Z") + contents[i] = LeadingPt/Pt; + else if (title=="p_{T}^{corr}") + contents[i] = CorrPt; + else if (title=="p_{T}^{MC}") + contents[i] = MCPt; + else if (title=="p_{T,particle}^{leading} (GeV/c)") + contents[i] = LeadingPt; + else + AliWarning(Form("Unable to fill dimension %s!",title.Data())); } - else { - TClass *cl = fMCTracks->GetClass(); - if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) { - AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data())); - fMCTracks = 0; - return; - } - } - } - AliAnalysisTaskEmcalJet::ExecOnce(); - - if (fMCFiducial) { - fMCminEta = fJetMinEta; - fMCmaxEta = fJetMaxEta; - fMCminPhi = fJetMinPhi; - fMCmaxPhi = fJetMaxPhi; + histo->Fill(contents); } else { - fMCminEta = fJetMinEta - fJetRadius; - fMCmaxEta = fJetMaxEta + fJetRadius; - fMCminPhi = fJetMinPhi - fJetRadius; - fMCmaxPhi = fJetMaxPhi + fJetRadius; + if (Set == 1) { + fHistJets1PtArea->Fill(A, Pt); + fHistJets1PhiEta->Fill(Eta, Phi); + fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt); + fHistJets1NEFvsPt->Fill(NEF, Pt); + fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt); + if (fIsJet1Rho) + fHistJets1CorrPtArea->Fill(A, CorrPt); + } + else if (Set == 2) { + fHistJets2PtArea->Fill(A, Pt); + fHistJets2PhiEta->Fill(Eta, Phi); + fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt); + fHistJets2NEFvsPt->Fill(NEF, Pt); + fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt); + if (fIsJet2Rho) + fHistJets2CorrPtArea->Fill(A, CorrPt); + } } } //________________________________________________________________________ -Bool_t AliJetResponseMaker::IsEventSelected() +void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, + Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, + Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2) { - // Check if event is selected + if (fHistoType==1) { + Double_t contents[20]={0}; + + for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) { + TString title(fHistMatching->GetAxis(i)->GetTitle()); + if (title=="p_{T,1}") + contents[i] = Pt1; + else if (title=="p_{T,2}") + contents[i] = Pt2; + else if (title=="A_{jet,1}") + contents[i] = A1; + else if (title=="A_{jet,2}") + contents[i] = A2; + else if (title=="distance") + contents[i] = d; + else if (title=="CE1") + contents[i] = CE1; + else if (title=="CE2") + contents[i] = CE2; + else if (title=="#deltaA_{jet}") + contents[i] = A1-A2; + else if (title=="#deltap_{T}") + contents[i] = Pt1-Pt2; + else if (title=="#delta#eta") + contents[i] = Eta1-Eta2; + else if (title=="#delta#phi") + contents[i] = Phi1-Phi2; + else if (title=="p_{T,1}^{corr}") + contents[i] = CorrPt1; + else if (title=="p_{T,2}^{corr}") + contents[i] = CorrPt2; + else if (title=="#deltap_{T}^{corr}") + contents[i] = CorrPt1-CorrPt2; + else if (title=="p_{T,1}^{MC}") + contents[i] = MCPt1; + else if (title=="#deltap_{T}^{MC}") + contents[i] = MCPt1-Pt2; + else if (title=="NEF_{1}") + contents[i] = NEF1; + else if (title=="NEF_{2}") + contents[i] = NEF2; + else if (title=="Z_{1}") + contents[i] = LeadingPt1/Pt1; + else if (title=="Z_{2}") + contents[i] = LeadingPt2/Pt2; + else if (title=="p_{T,particle,1}^{leading} (GeV/c)") + contents[i] = LeadingPt1; + else if (title=="p_{T,particle,2}^{leading} (GeV/c)") + contents[i] = LeadingPt2; + else + AliWarning(Form("Unable to fill dimension %s!",title.Data())); + } - if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) - return kFALSE; + fHistMatching->Fill(contents); + } + else { + fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1); + fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2); + + fHistDistancevsJet1Pt->Fill(d, Pt1); + fHistDistancevsJet2Pt->Fill(d, Pt2); + + fHistDistancevsCommonEnergy1->Fill(d, CE1); + fHistDistancevsCommonEnergy2->Fill(d, CE2); + fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2); + + fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2); + + fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1); + fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2); + + Double_t dpt = Pt1 - Pt2; + fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt); + fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt); + fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1); + fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2); + + fHistDeltaPtvsDistance->Fill(d, dpt); + fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt); + fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt); + + fHistDeltaPtvsArea1->Fill(A1, dpt); + fHistDeltaPtvsArea2->Fill(A2, dpt); + + Double_t darea = A1 - A2; + fHistDeltaPtvsDeltaArea->Fill(darea, dpt); + + fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2); + + if (fIsJet1Rho || fIsJet2Rho) { + Double_t dcorrpt = CorrPt1 - CorrPt2; + fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt); + fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt); + fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1); + fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2); + fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt); + fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt); + fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt); + fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt); + fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt); + fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt); + fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2); + } - return AliAnalysisTaskEmcalJet::IsEventSelected(); + if (fIsEmbedded) { + Double_t dmcpt = MCPt1 - Pt2; + fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt); + fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt); + fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1); + fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2); + fHistDeltaMCPtvsDistance->Fill(d, dmcpt); + fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt); + fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt); + fHistDeltaMCPtvsArea1->Fill(A1, dmcpt); + fHistDeltaMCPtvsArea2->Fill(A2, dmcpt); + fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt); + fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2); + } + } } //________________________________________________________________________ -Bool_t AliJetResponseMaker::RetrieveEventObjects() +void AliJetResponseMaker::ExecOnce() { - // Retrieve event objects. - - if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) - return kFALSE; - - fPythiaHeader = dynamic_cast(MCEvent()->GenEventHeader()); - - if (!fPythiaHeader) - return kFALSE; + // Execute once. - if (fDoWeighting) - fEventWeight = fPythiaHeader->EventWeight(); - else - fEventWeight = 1; + AliAnalysisTaskEmcalJet::ExecOnce(); - const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234}; - const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000}; + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); - Double_t pthard = fPythiaHeader->GetPtHard(); + if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; - for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) { - if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin]) - break; + if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) { + if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) { + AliWarning("Changing matching type from MC label to same collection..."); + fMatching = kSameCollections; + } + else { + AliWarning("Changing matching type from MC label to geometrical..."); + fMatching = kGeometrical; + } + } + else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) { + if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) { + AliWarning("Changing matching type from same collection to MC label..."); + fMatching = kMCLabel; + } + else { + AliWarning("Changing matching type from same collection to geometrical..."); + fMatching = kGeometrical; + } } - - fNTrials = fPythiaHeader->Trials(); - - return kTRUE; } //________________________________________________________________________ @@ -399,231 +1108,588 @@ Bool_t AliJetResponseMaker::Run() { // Find the closest jets - if (!fDoMatching) + if (fMatching == kNoMatching) return kTRUE; + else + return DoJetMatching(); +} - DoJetLoop(fJets, fMCJets, kFALSE); - DoJetLoop(fMCJets, fJets, kTRUE); - - const Int_t nMCJets = fMCJets->GetEntriesFast(); +//________________________________________________________________________ +Bool_t AliJetResponseMaker::DoJetMatching() +{ + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); - for (Int_t i = 0; i < nMCJets; i++) { + if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE; - AliEmcalJet* jet = static_cast(fMCJets->At(i)); + DoJetLoop(); - if (!jet) { - AliError(Form("Could not receive jet %d", i)); - continue; - } + AliEmcalJet* jet1 = 0; - if (!AcceptJet(jet)) - continue; + jets1->ResetCurrentID(); + while ((jet1 = jets1->GetNextJet())) { - if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi) - continue; + AliEmcalJet *jet2 = jet1->ClosestJet(); - if (jet->Pt() > fMaxBinPt) - continue; + if (!jet2) continue; + if (jet2->ClosestJet() != jet1) continue; + if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue; - if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && - jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found - jet->SetMatchedToClosest(); - jet->ClosestJet()->SetMatchedToClosest(); - } + // Matched jet found + jet1->SetMatchedToClosest(fMatching); + jet2->SetMatchedToClosest(fMatching); + AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f", + jet1->Pt(), jet1->Eta(), jet1->Phi(), + jet2->Pt(), jet2->Eta(), jet2->Phi())); } return kTRUE; } //________________________________________________________________________ -void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc) +void AliJetResponseMaker::DoJetLoop() { // Do the jet loop. - Int_t nJets1 = jets1->GetEntriesFast(); - Int_t nJets2 = jets2->GetEntriesFast(); + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); - for (Int_t i = 0; i < nJets1; i++) { + if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; - AliEmcalJet* jet1 = static_cast(jets1->At(i)); + AliEmcalJet* jet1 = 0; + AliEmcalJet* jet2 = 0; - if (!jet1) { - AliError(Form("Could not receive jet %d", i)); - continue; - } + jets2->ResetCurrentID(); + while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching(); + + jets1->ResetCurrentID(); + while ((jet1 = jets1->GetNextJet())) { + jet1->ResetMatching(); + + if (jet1->MCPt() < fMinJetMCPt) continue; - if (!AcceptJet(jet1)) - continue; + jets2->ResetCurrentID(); + while ((jet2 = jets2->GetNextJet())) { + SetMatchingLevel(jet1, jet2, fMatching); + } // jet2 loop + } // jet1 loop +} + +//________________________________________________________________________ +void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const +{ + Double_t deta = jet2->Eta() - jet1->Eta(); + Double_t dphi = jet2->Phi() - jet1->Phi(); + d = TMath::Sqrt(deta * deta + dphi * dphi); +} - if (!mc) { - if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi) +//________________________________________________________________________ +void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const +{ + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); + + if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; + + AliParticleContainer *tracks1 = jets1->GetParticleContainer(); + AliClusterContainer *clusters1 = jets1->GetClusterContainer(); + AliParticleContainer *tracks2 = jets2->GetParticleContainer(); + + // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated + d1 = jet1->Pt(); + d2 = jet2->Pt(); + Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background + + // remove completely tracks that are not MC particles (label == 0) + if (tracks1 && tracks1->GetArray()) { + for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) { + AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray()); + if (!track) { + AliWarning(Form("Could not find track %d!", iTrack)); continue; + } + + Int_t MClabel = TMath::Abs(track->GetLabel()); + MClabel -= fMCLabelShift; + if (MClabel != 0) continue; + + // this is not a MC particle; remove it completely + AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel)); + totalPt1 -= track->Pt(); + d1 -= track->Pt(); } - else { - if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi) + } + + // remove completely clusters that are not MC particles (label == 0) + if (fUseCellsToMatch && fCaloCells) { + for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { + AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray()); + if (!clus) { + AliWarning(Form("Could not find cluster %d!", iClus)); continue; + } + TLorentzVector part; + clus->GetMomentum(part, const_cast(fVertex)); + + for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) { + Int_t cellId = clus->GetCellAbsId(iCell); + Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell); + + Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId)); + MClabel -= fMCLabelShift; + if (MClabel != 0) continue; + + // this is not a MC particle; remove it completely + AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel)); + totalPt1 -= part.Pt() * cellFrac; + d1 -= part.Pt() * cellFrac; + } + } + } + else { + for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { + AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray()); + if (!clus) { + AliWarning(Form("Could not find cluster %d!", iClus)); + continue; + } + TLorentzVector part; + clus->GetMomentum(part, const_cast(fVertex)); + + Int_t MClabel = TMath::Abs(clus->GetLabel()); + MClabel -= fMCLabelShift; + if (MClabel != 0) continue; + + // this is not a MC particle; remove it completely + AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel)); + totalPt1 -= part.Pt(); + d1 -= part.Pt(); } + } - for (Int_t j = 0; j < nJets2; j++) { - - AliEmcalJet* jet2 = static_cast(jets2->At(j)); - - if (!jet2) { - AliError(Form("Could not receive jet %d", j)); + for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) { + Bool_t track2Found = kFALSE; + Int_t index2 = jet2->TrackAt(iTrack2); + + // now look for common particles in the track array + for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) { + AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray()); + if (!track) { + AliWarning(Form("Could not find track %d!", iTrack)); continue; - } - - if (!AcceptJet(jet2)) + } + Int_t MClabel = TMath::Abs(track->GetLabel()); + MClabel -= fMCLabelShift; + if (MClabel <= 0) continue; + + Int_t index = -1; + index = tracks2->GetIndexFromLabel(MClabel); + if (index < 0) { + AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel)); continue; + } + + if (index2 != index) continue; + + // found common particle + d1 -= track->Pt(); + + if (!track2Found) { + AliVParticle *MCpart = tracks2->GetParticle(index2); + AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!", + iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi())); + d2 -= MCpart->Pt(); + } + + track2Found = kTRUE; + } - if (mc) { - if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi) + // now look for common particles in the cluster array + if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle + for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { + AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters); + if (!clus) { + AliWarning(Form("Could not find cluster %d!", iClus)); continue; + } + TLorentzVector part; + clus->GetMomentum(part, const_cast(fVertex)); + + for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) { + Int_t cellId = clus->GetCellAbsId(iCell); + Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell); + + Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId)); + MClabel -= fMCLabelShift; + if (MClabel <= 0) continue; + + Int_t index1 = -1; + index1 = tracks2->GetIndexFromLabel(MClabel); + if (index1 < 0) { + AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel)); + continue; + } + + if (index2 != index1) continue; + + // found common particle + d1 -= part.Pt() * cellFrac; + + if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found) + AliVParticle *MCpart = tracks2->GetParticle(index2); + AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!", + iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi())); + d2 -= MCpart->Pt() * cellFrac; + } + + track2Found = kTRUE; + } } - else { - if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi) + } + else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it + for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { + AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters); + if (!clus) { + AliWarning(Form("Could not find cluster %d!", iClus)); continue; - } + } + TLorentzVector part; + clus->GetMomentum(part, const_cast(fVertex)); - Double_t deta = jet2->Eta() - jet1->Eta(); - Double_t dphi = jet2->Phi() - jet1->Phi(); - Double_t d = TMath::Sqrt(deta * deta + dphi * dphi); + Int_t MClabel = TMath::Abs(clus->GetLabel()); + MClabel -= fMCLabelShift; + if (MClabel <= 0) continue; - if (d < jet1->ClosestJetDistance()) { - jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance()); - jet1->SetClosestJet(jet2, d); - } - else if (d < jet1->SecondClosestJetDistance()) { - jet1->SetSecondClosestJet(jet2, d); + Int_t index = -1; + index = tracks2->GetIndexFromLabel(MClabel); + + if (index < 0) { + AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel)); + continue; + } + + if (index2 != index) continue; + + // found common particle + d1 -= part.Pt(); + + if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found) + AliVParticle *MCpart = tracks2->GetParticle(index2); + AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!", + iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi())); + + d2 -= MCpart->Pt(); + } + + track2Found = kTRUE; } } } -} -//________________________________________________________________________ -Bool_t AliJetResponseMaker::FillHistograms() -{ - // Fill histograms. + if (d1 < 0) + d1 = 0; - fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1); - fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials); - if (fEventWeightHist) - fHistEventWeight[fPtHardBin]->Fill(fPythiaHeader->EventWeight()); + if (d2 < 0) + d2 = 0; - const Int_t nMCJets = fMCJets->GetEntriesFast(); - - for (Int_t i = 0; i < nMCJets; i++) { - - AliEmcalJet* jet = static_cast(fMCJets->At(i)); + if (totalPt1 < 1) + d1 = -1; + else + d1 /= totalPt1; - if (!jet) { - AliError(Form("Could not receive jet %d", i)); - continue; - } + if (jet2->Pt() < 1) + d2 = -1; + else + d2 /= jet2->Pt(); +} - if (!AcceptJet(jet)) - continue; +//________________________________________________________________________ +void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const +{ + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); + + if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; + + AliParticleContainer *tracks1 = jets1->GetParticleContainer(); + AliClusterContainer *clusters1 = jets1->GetClusterContainer(); + AliParticleContainer *tracks2 = jets2->GetParticleContainer(); + AliClusterContainer *clusters2 = jets2->GetClusterContainer(); + + // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated + d1 = jet1->Pt(); + d2 = jet2->Pt(); + + if (tracks1 && tracks2) { + + for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) { + Int_t index2 = jet2->TrackAt(iTrack2); + for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) { + Int_t index1 = jet1->TrackAt(iTrack1); + if (index2 == index1) { // found common particle + AliVParticle *part1 = tracks1->GetParticle(index1); + if (!part1) { + AliWarning(Form("Could not find track %d!", index1)); + continue; + } + AliVParticle *part2 = tracks2->GetParticle(index2); + if (!part2) { + AliWarning(Form("Could not find track %d!", index2)); + continue; + } + + d1 -= part1->Pt(); + d2 -= part2->Pt(); + break; + } + } + } - if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi) - continue; + } - if (jet->Pt() > fMaxBinPt) - continue; + if (clusters1 && clusters2) { + + if (fUseCellsToMatch) { + const Int_t nClus1 = jet1->GetNumberOfClusters(); + + Int_t ncells1[nClus1]; + UShort_t *cellsId1[nClus1]; + Double_t *cellsFrac1[nClus1]; + Int_t *sortedIndexes1[nClus1]; + Double_t ptClus1[nClus1]; + for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) { + Int_t index1 = jet1->ClusterAt(iClus1); + AliVCluster *clus1 = clusters1->GetCluster(index1); + if (!clus1) { + AliWarning(Form("Could not find cluster %d!", index1)); + ncells1[iClus1] = 0; + cellsId1[iClus1] = 0; + cellsFrac1[iClus1] = 0; + sortedIndexes1[iClus1] = 0; + ptClus1[iClus1] = 0; + continue; + } + TLorentzVector part1; + clus1->GetMomentum(part1, const_cast(fVertex)); - if (jet->MatchedJet()) { + ncells1[iClus1] = clus1->GetNCells(); + cellsId1[iClus1] = clus1->GetCellsAbsId(); + cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction(); + sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]]; + ptClus1[iClus1] = part1.Pt(); - if (!AcceptBiasJet(jet->MatchedJet()) || - jet->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet->MatchedJet()->MaxClusterPt() > fMaxClusterPt || - jet->MatchedJet()->Pt() > fMaxBinPt) { - fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); + TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE); } - else { - fHistMatchingLevelMCPt->Fill(jet->ClosestJetDistance(), jet->Pt(), fEventWeight); - - Double_t deta = jet->MatchedJet()->Eta() - jet->Eta(); - Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi(); - fHistClosestDeltaEtaPhiMCPt->Fill(deta, dphi, jet->Pt(), fEventWeight); - - Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt(); - fHistClosestDeltaPtMCPt->Fill(jet->Pt(), dpt, fEventWeight); - fHistClosestDeltaCorrPtMCPt->Fill(jet->Pt(), dpt - fRhoVal * jet->MatchedJet()->Area(), fEventWeight); - - fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight); - fHistPartvsDetecCorrPt->Fill(jet->MatchedJet()->Pt() - fRhoVal * jet->MatchedJet()->Area(), jet->Pt(), fEventWeight); + + const Int_t nClus2 = jet2->GetNumberOfClusters(); + + const Int_t maxNcells2 = 11520; + Int_t sortedIndexes2[maxNcells2]; + for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) { + Int_t index2 = jet2->ClusterAt(iClus2); + AliVCluster *clus2 = clusters2->GetCluster(index2); + if (!clus2) { + AliWarning(Form("Could not find cluster %d!", index2)); + continue; + } + Int_t ncells2 = clus2->GetNCells(); + if (ncells2 >= maxNcells2) { + AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2)); + continue; + } + UShort_t *cellsId2 = clus2->GetCellsAbsId(); + Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction(); + + TLorentzVector part2; + clus2->GetMomentum(part2, const_cast(fVertex)); + Double_t ptClus2 = part2.Pt(); + + TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE); + + for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) { + if (sortedIndexes1[iClus1] == 0) + continue; + Int_t iCell1 = 0, iCell2 = 0; + while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) { + if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell + d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1]; + d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2; + iCell1++; + iCell2++; + } + else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) { + iCell2++; + } + else { + iCell1++; + } + } + } } } else { - fHistNonMatchedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); - fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); + for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) { + Int_t index2 = jet2->ClusterAt(iClus2); + for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) { + Int_t index1 = jet1->ClusterAt(iClus1); + if (index2 == index1) { // found common particle + AliVCluster *clus1 = clusters1->GetCluster(index1); + if (!clus1) { + AliWarning(Form("Could not find cluster %d!", index1)); + continue; + } + AliVCluster *clus2 = clusters2->GetCluster(index2); + if (!clus2) { + AliWarning(Form("Could not find cluster %d!", index2)); + continue; + } + TLorentzVector part1, part2; + clus1->GetMomentum(part1, const_cast(fVertex)); + clus2->GetMomentum(part2, const_cast(fVertex)); + + d1 -= part1.Pt(); + d2 -= part2.Pt(); + break; + } + } + } } + } - fHistMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); - fHistMCJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight); + if (d1 < 0) + d1 = 0; - fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight); - - for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) { - AliVParticle *track = jet->TrackAt(it, fMCTracks); - if (track) - fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight); - } + if (d2 < 0) + d2 = 0; - if (!AcceptBiasJet(jet)) - continue; - if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi) - continue; + if (jet1->Pt() > 0) + d1 /= jet1->Pt(); + else + d1 = -1; + + if (jet2->Pt() > 0) + d2 /= jet2->Pt(); + else + d2 = -1; +} + +//________________________________________________________________________ +void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) +{ + Double_t d1 = -1; + Double_t d2 = -1; + + switch (matching) { + case kGeometrical: + GetGeometricalMatchingLevel(jet1,jet2,d1); + d2 = d1; + break; + case kMCLabel: // jet1 = detector level and jet2 = particle level! + GetMCLabelMatchingLevel(jet1,jet2,d1,d2); + break; + case kSameCollections: + GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2); + break; + default: + ; + } + + if (d1 >= 0) { + + if (d1 < jet1->ClosestJetDistance()) { + jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance()); + jet1->SetClosestJet(jet2, d1); + } + else if (d1 < jet1->SecondClosestJetDistance()) { + jet1->SetSecondClosestJet(jet2, d1); + } + } + + if (d2 >= 0) { - fHistMCJetsPtAreaFiducial->Fill(jet->Area(), jet->Pt(), fEventWeight); - fHistMCJetsPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight); + if (d2 < jet2->ClosestJetDistance()) { + jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance()); + jet2->SetClosestJet(jet1, d2); + } + else if (d2 < jet2->SecondClosestJetDistance()) { + jet2->SetSecondClosestJet(jet1, d2); + } } +} - const Int_t nJets = fJets->GetEntriesFast(); +//________________________________________________________________________ +Bool_t AliJetResponseMaker::FillHistograms() +{ + // Fill histograms. - for (Int_t i = 0; i < nJets; i++) { + AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); + AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); - AliEmcalJet* jet = static_cast(fJets->At(i)); + if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE; - if (!jet) { - AliError(Form("Could not receive mc jet %d", i)); - continue; - } - - if (!AcceptJet(jet)) - continue; - if (!AcceptBiasJet(jet)) - continue; - if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt) - continue; - if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi) - continue; + AliEmcalJet* jet1 = 0; + AliEmcalJet* jet2 = 0; + + jets2->ResetCurrentID(); + while ((jet2 = jets2->GetNextJet())) { - if (!jet->MatchedJet()) { - fHistNonMatchedJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); - fHistNonMatchedJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight); - } + AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID())); + + if (jet2->Pt() < jets2->GetJetPtCut()) continue; - fHistJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); - fHistJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight); + Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2); + Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area(); - fHistJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight); + if (jets2->AcceptJet(jet2)) + FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2, + corrpt2, jet2->MCPt(), 2); + else + fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(), jet2->Pt()); - fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight); + jet1 = jet2->MatchedJet(); - for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) { - AliVParticle *track = jet->TrackAt(it, fTracks); - if (track) - fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight); + if (!jet1) continue; + if (!jets1->AcceptJet(jet1)) continue; + if (jet1->MCPt() < fMinJetMCPt) continue; + + Double_t d=-1, ce1=-1, ce2=-1; + if (jet2->GetMatchingType() == kGeometrical) { + if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported + GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2); + else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) + GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2); + + d = jet2->ClosestJetDistance(); } + else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) { + GetGeometricalMatchingLevel(jet1, jet2, d); - for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) { - AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters); - if (cluster) { - TLorentzVector nP; - cluster->GetMomentum(nP, fVertex); - fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight); - } + ce1 = jet1->ClosestJetDistance(); + ce2 = jet2->ClosestJetDistance(); } + + Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1); + Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area(); + + FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(), + jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(), + jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2); } + jets1->ResetCurrentID(); + while ((jet1 = jets1->GetNextJet())) { + if (!jets1->AcceptJet(jet1)) { + fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(), jet1->Pt()); + continue; + } + if (jet1->MCPt() < fMinJetMCPt) continue; + AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID())); + + Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1); + Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area(); + + FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1, + corrpt1, jet1->MCPt(), 1); + } return kTRUE; }