#include "AliJetResponseMaker.h"
-#include <TChain.h>
#include <TClonesArray.h>
-#include <TH1F.h>
#include <TH2F.h>
-#include <TList.h>
+#include <THnSparse.h>
#include <TLorentzVector.h>
#include "AliAnalysisManager.h"
-#include "AliCentrality.h"
-#include "AliFJWrapper.h"
#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
-#include "AliGenPythiaEventHeader.h"
-#include "AliMCEvent.h"
+#include "AliLog.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),
- fHistMCJetPhiEta(0),
- fHistMCJetsPt(0),
- fHistMCJetPhiEtaFiducial(0),
- fHistMCJetsPtFiducial(0),
- fHistMCJetsNEFvsPt(0),
- fHistMCJetsZvsPt(0),
- fHistJetPhiEta(0),
- fHistJetsPt(0),
- fHistJetsNEFvsPt(0),
- fHistJetsZvsPt(0),
- fHistClosestDistance(0),
- fHistClosestDeltaPhi(0),
- fHistClosestDeltaEta(0),
- fHistClosestDeltaPt(0),
- fHistNonMatchedMCJetPt(0),
- fHistNonMatchedJetPt(0),
- fHistPartvsDetecPt(0),
- fHistMissedMCJets(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),
+ 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),
- fHistMCJetPhiEta(0),
- fHistMCJetsPt(0),
- fHistMCJetPhiEtaFiducial(0),
- fHistMCJetsPtFiducial(0),
- fHistMCJetsNEFvsPt(0),
- fHistMCJetsZvsPt(0),
- fHistJetPhiEta(0),
- fHistJetsPt(0),
- fHistJetsNEFvsPt(0),
- fHistJetsZvsPt(0),
- fHistClosestDistance(0),
- fHistClosestDeltaPhi(0),
- fHistClosestDeltaEta(0),
- fHistClosestDeltaPt(0),
- fHistNonMatchedMCJetPt(0),
- fHistNonMatchedJetPt(0),
- fHistPartvsDetecPt(0),
- fHistMissedMCJets(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),
+ 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);
}
// 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++;
}
- fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
- fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
- fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistJetPhiEta);
-
- fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
- fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistJetsPt->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistJetsPt);
-
- fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
- fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
- fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [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} [GeV/c]");
- fOutput->Add(fHistJetsNEFvsPt);
-
- fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
- fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
- fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistMCJetPhiEta);
-
- fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
- fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistMCJetsPt->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistMCJetsPt);
-
- fHistMCJetPhiEtaFiducial = new TH2F("fHistMCJetPhiEtaFiducial", "fHistMCJetPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4);
- fHistMCJetPhiEtaFiducial->GetXaxis()->SetTitle("#eta");
- fHistMCJetPhiEtaFiducial->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistMCJetPhiEtaFiducial);
-
- fHistMCJetsPtFiducial = new TH1F("fHistMCJetsPtFiducial", "fHistMCJetsPtFiducial", fNbins, fMinBinPt, fMaxBinPt);
- fHistMCJetsPtFiducial->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistMCJetsPtFiducial->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistMCJetsPtFiducial);
+ 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} [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} [GeV/c]");
- fOutput->Add(fHistMCJetsNEFvsPt);
-
- fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
- fHistClosestDistance->GetXaxis()->SetTitle("d");
- fHistClosestDistance->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistClosestDistance);
-
- fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
- fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
- fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistClosestDeltaPhi);
-
- fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2);
- fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
- fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistClosestDeltaEta);
-
- fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
- fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
- fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistClosestDeltaPt);
-
- fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
- fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistNonMatchedMCJetPt);
-
- fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
- fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistNonMatchedJetPt);
-
- 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);
-
- fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
- fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistMissedMCJets->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistMissedMCJets);
+ // 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<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(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;
+
+ 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<TClonesArray*>(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<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
- if (!fMCTracks) {
- AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data()));
- return;
- }
- 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;
- }
+ 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()));
}
- }
-
- 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<AliGenPythiaEventHeader*>(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<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(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;
}
//________________________________________________________________________
{
// 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<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- for (Int_t i = 0; i < nMCJets; i++) {
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(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<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
+
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
- for (Int_t i = 0; i < nJets1; i++) {
+ AliEmcalJet* jet1 = 0;
+ AliEmcalJet* jet2 = 0;
+
+ jets2->ResetCurrentID();
+ while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
+
+ jets1->ResetCurrentID();
+ while ((jet1 = jets1->GetNextJet())) {
+ jet1->ResetMatching();
+
+ if (jet1->MCPt() < fMinJetMCPt) 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);
+}
- AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
+//________________________________________________________________________
+void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
+{
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(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;
+ }
- if (!jet1) {
- AliError(Form("Could not receive jet %d", i));
- continue;
- }
+ Int_t MClabel = TMath::Abs(track->GetLabel());
+ MClabel -= fMCLabelShift;
+ if (MClabel != 0) continue;
- if (!AcceptJet(jet1))
- 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();
+ }
+ }
- if (!mc) {
- if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
+ // 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<Double_t*>(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 {
- if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi)
+ }
+ 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<Double_t*>(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<AliEmcalJet*>(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();
+ }
- if (mc) {
- if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
+ track2Found = kTRUE;
+ }
+
+ // 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<Double_t*>(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<Double_t*>(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.
-
- 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());
-
- const Int_t nMCJets = fMCJets->GetEntriesFast();
- for (Int_t i = 0; i < nMCJets; i++) {
+ if (d1 < 0)
+ d1 = 0;
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
+ if (d2 < 0)
+ d2 = 0;
- if (!jet) {
- AliError(Form("Could not receive jet %d", i));
- continue;
- }
-
- if (!AcceptJet(jet))
- continue;
-
- if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi)
- continue;
-
- if (jet->Pt() > fMaxBinPt)
- continue;
+ if (totalPt1 < 1)
+ d1 = -1;
+ else
+ d1 /= totalPt1;
- if (jet->MatchedJet()) {
+ if (jet2->Pt() < 1)
+ d2 = -1;
+ else
+ d2 /= jet2->Pt();
+}
- if (!AcceptBiasJet(jet->MatchedJet()) ||
- jet->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
- jet->MatchedJet()->Pt() > fMaxBinPt) {
- fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
+//________________________________________________________________________
+void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
+{
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(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;
+ }
}
- else {
- fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
+ }
- Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
- fHistClosestDeltaEta->Fill(deta, fEventWeight);
+ }
- Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
- fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
+ 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<Double_t*>(fVertex));
- Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
- fHistClosestDeltaPt->Fill(dpt, fEventWeight);
+ ncells1[iClus1] = clus1->GetNCells();
+ cellsId1[iClus1] = clus1->GetCellsAbsId();
+ cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
+ sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
+ ptClus1[iClus1] = part1.Pt();
- fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
+ TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
+ }
+
+ 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<Double_t*>(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 {
- fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
- fHistMissedMCJets->Fill(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<Double_t*>(fVertex));
+ clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
+
+ d1 -= part1.Pt();
+ d2 -= part2.Pt();
+ break;
+ }
+ }
+ }
}
+ }
- fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
- fHistMCJetPhiEta->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;
-
- fHistMCJetsPtFiducial->Fill(jet->Pt(), fEventWeight);
- fHistMCJetPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight);
- }
+ if (jet1->Pt() > 0)
+ d1 /= jet1->Pt();
+ else
+ d1 = -1;
- const Int_t nJets = fJets->GetEntriesFast();
+ if (jet2->Pt() > 0)
+ d2 /= jet2->Pt();
+ else
+ d2 = -1;
+}
- for (Int_t i = 0; i < nJets; i++) {
+//________________________________________________________________________
+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:
+ ;
+ }
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
+ if (d1 >= 0) {
- if (!jet) {
- AliError(Form("Could not receive mc jet %d", i));
- continue;
- }
+ 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) {
- 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;
-
- if (!jet->MatchedJet()) {
- fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
+ if (d2 < jet2->ClosestJetDistance()) {
+ jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
+ jet2->SetClosestJet(jet1, d2);
}
+ else if (d2 < jet2->SecondClosestJetDistance()) {
+ jet2->SetSecondClosestJet(jet1, d2);
+ }
+ }
+}
- fHistJetsPt->Fill(jet->Pt(), fEventWeight);
+//________________________________________________________________________
+Bool_t AliJetResponseMaker::FillHistograms()
+{
+ // Fill histograms.
+
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
- fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
+ AliEmcalJet* jet1 = 0;
+ AliEmcalJet* jet2 = 0;
+
+ jets2->ResetCurrentID();
+ while ((jet2 = jets2->GetNextJet())) {
- 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);
+ AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
+
+ Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
+ Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
+
+ if (jets2->AcceptJet(jet2))
+ FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2,
+ corrpt2, jet2->MCPt(), 2);
+
+ jet1 = jet2->MatchedJet();
+
+ 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->GetNextAcceptJet())) {
+ 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;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff