#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
+#include <TH3F.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TRandom3.h>
//________________________________________________________________________
AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
- fMCAna(kFALSE),
+ fMCJetPtThreshold(1),
fMinRC2LJ(-1),
fEmbJetsName(""),
fEmbTracksName(""),
fEmbCaloClusters(0),
fRandTracks(0),
fRandCaloClusters(0),
- fEmbeddedClusterId(-1),
- fEmbeddedTrackId(-1),
+ fEmbeddedClusterNIds(0),
+ fEmbeddedTrackNIds(0),
fHistRCPhiEta(0),
fHistRCPtExLJVSDPhiLJ(0),
- fHistEmbJetPhiEta(0),
- fHistEmbPartPhiEta(0)
+ fHistEmbJetsPhiEta(0),
+ fHistLeadPartPhiEta(0)
{
// Default constructor.
fHistDeltaPtRC[i] = 0;
fHistDeltaPtRCExLJ[i] = 0;
fHistDeltaPtRCRand[i] = 0;
+ fHistEmbRejectedJetsPhiEta[i] = 0;
+ fHistEmbRejectedJetsPtArea[i] = 0;
+ fHistEmbNotFoundPt[i] = 0;
fHistEmbNotFoundPhiEta[i] = 0;
fHistEmbJetsPtArea[i] = 0;
fHistEmbJetsCorrPtArea[i] = 0;
- fHistEmbPartPt[i] = 0;
- fHistDistEmbPartJetAxis[i] = 0;
+ fHistEmbPartPtvsJetPt[i] = 0;
+ fHistEmbPartPtvsJetCorrPt[i] = 0;
+ fHistJetPtvsJetCorrPt[i] = 0;
+ fHistDistLeadPart2JetAxis[i] = 0;
fHistEmbBkgArea[i] = 0;
fHistRhoVSEmbBkg[i] = 0;
fHistDeltaPtEmbArea[i] = 0;
}
+ memset(fEmbeddedClusterIds, -1, 999*sizeof(Int_t));
+ memset(fEmbeddedTrackIds, -1, 999*sizeof(Int_t));
+
SetMakeGeneralHistograms(kTRUE);
}
//________________________________________________________________________
AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
AliAnalysisTaskEmcalJet(name, kTRUE),
- fMCAna(kFALSE),
+ fMCJetPtThreshold(1),
fMinRC2LJ(-1),
fEmbJetsName(""),
fEmbTracksName(""),
fEmbCaloClusters(0),
fRandTracks(0),
fRandCaloClusters(0),
- fEmbeddedClusterId(-1),
- fEmbeddedTrackId(-1),
+ fEmbeddedClusterNIds(0),
+ fEmbeddedTrackNIds(0),
fHistRCPhiEta(0),
fHistRCPtExLJVSDPhiLJ(0),
- fHistEmbJetPhiEta(0),
- fHistEmbPartPhiEta(0)
+ fHistEmbJetsPhiEta(0),
+ fHistLeadPartPhiEta(0)
{
// Standard constructor.
fHistDeltaPtRC[i] = 0;
fHistDeltaPtRCExLJ[i] = 0;
fHistDeltaPtRCRand[i] = 0;
+ fHistEmbRejectedJetsPhiEta[i] = 0;
+ fHistEmbRejectedJetsPtArea[i] = 0;
+ fHistEmbNotFoundPt[i] = 0;
fHistEmbNotFoundPhiEta[i] = 0;
fHistEmbJetsPtArea[i] = 0;
fHistEmbJetsCorrPtArea[i] = 0;
- fHistEmbPartPt[i] = 0;
- fHistDistEmbPartJetAxis[i] = 0;
+ fHistEmbPartPtvsJetPt[i] = 0;
+ fHistEmbPartPtvsJetCorrPt[i] = 0;
+ fHistJetPtvsJetCorrPt[i] = 0;
+ fHistDistLeadPart2JetAxis[i] = 0;
fHistEmbBkgArea[i] = 0;
fHistRhoVSEmbBkg[i] = 0;
fHistDeltaPtEmbArea[i] = 0;
}
+ memset(fEmbeddedClusterIds, -1, 999*sizeof(Int_t));
+ memset(fEmbeddedTrackIds, -1, 999*sizeof(Int_t));
+
SetMakeGeneralHistograms(kTRUE);
}
AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
- fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
- fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
- fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistRCPhiEta);
-
- fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
- fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone #it{p}_{T} (GeV/#it{c})");
- fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
- fOutput->Add(fHistRCPtExLJVSDPhiLJ);
-
- if (!fJetsName.IsNull()) {
- fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
- fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
- fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistEmbJetPhiEta);
+ if (!fTracksName.IsNull() || !fCaloName.IsNull()) {
+ fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+ fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
+ fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
+ fOutput->Add(fHistRCPhiEta);
+
+ if (!fJetsName.IsNull()) {
+ fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
+ fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
+ fOutput->Add(fHistRCPtExLJVSDPhiLJ);
+ }
+ }
+
+ if (!fEmbJetsName.IsNull()) {
+ fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+ fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
+ fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
+ fOutput->Add(fHistEmbJetsPhiEta);
- fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
- fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
- fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistEmbPartPhiEta);
+ fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+ fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
+ fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
+ fOutput->Add(fHistLeadPartPhiEta);
}
TString histname;
- for (Int_t i = 0; i < 4; i++) {
- histname = "fHistRCPt_";
- histname += i;
- fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
- fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone #it{p}_{T} (GeV/#it{c})");
- fHistRCPt[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistRCPt[i]);
-
- histname = "fHistRCPtExLJ_";
- histname += i;
- fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
- fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone #it{p}_{T}^{RC} (GeV/#it{c})");
- fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistRCPtExLJ[i]);
-
- histname = "fHistRCPtRand_";
- histname += i;
- fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
- fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone #it{p}_{T}^{RC} (GeV/#it{c})");
- fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistRCPtRand[i]);
-
- histname = "fHistRhoVSRCPt_";
- histname += i;
- fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
- fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
- fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("rigid cone #it{p}_{T} (GeV/#it{c})");
- fOutput->Add(fHistRhoVSRCPt[i]);
-
- histname = "fHistDeltaPtRC_";
- histname += i;
- fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
- fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
- fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaPtRC[i]);
-
- histname = "fHistDeltaPtRCExLJ_";
- histname += i;
- fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
- fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
- fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaPtRCExLJ[i]);
-
- histname = "fHistDeltaPtRCRand_";
- histname += i;
- fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
- fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
- fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaPtRCRand[i]);
+ const Int_t nbinsZ = 12;
+ Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
+
+ Float_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
+ Float_t *binsCorrPt = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
+ Float_t *binsArea = GenerateFixedBinArray(50, 0, 2);
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ if (!fTracksName.IsNull() || !fCaloName.IsNull()) {
+ histname = "fHistRCPt_";
+ histname += i;
+ fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+ fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ fHistRCPt[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistRCPt[i]);
+
+ histname = "fHistRhoVSRCPt_";
+ histname += i;
+ fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+ fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
+ fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ fOutput->Add(fHistRhoVSRCPt[i]);
+
+ histname = "fHistDeltaPtRC_";
+ histname += i;
+ fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+ fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+ fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistDeltaPtRC[i]);
+
+ if (!fJetsName.IsNull()) {
+ histname = "fHistRCPtExLJ_";
+ histname += i;
+ fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+ fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
+ fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistRCPtExLJ[i]);
+
+ histname = "fHistDeltaPtRCExLJ_";
+ histname += i;
+ fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+ fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+ fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistDeltaPtRCExLJ[i]);
+ }
+ }
+
+ if (!fRandTracksName.IsNull() || !fRandCaloName.IsNull()) {
+ histname = "fHistRCPtRand_";
+ histname += i;
+ fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
+ fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
+ fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistRCPtRand[i]);
+
+ histname = "fHistDeltaPtRCRand_";
+ histname += i;
+ fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
+ fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
+ fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistDeltaPtRCRand[i]);
+ }
if (!fEmbJetsName.IsNull()) {
histname = "fHistEmbJetsPtArea_";
histname += i;
- fHistEmbJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+ fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
- fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("embedded jet #it{p}_{T}^{raw} (GeV/#it{c})");
+ fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
fOutput->Add(fHistEmbJetsPtArea[i]);
histname = "fHistEmbJetsCorrPtArea_";
histname += i;
- fHistEmbJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins * 2, -fMaxBinPt, fMaxBinPt);
+ fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
- fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("embedded jet #it{p}_{T}^{corr} (GeV/#it{c})");
+ fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
fOutput->Add(fHistEmbJetsCorrPtArea[i]);
- histname = "fHistEmbPartPt_";
+ histname = "fHistEmbPartPtvsJetPt_";
histname += i;
- fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
- fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle #it{p}_{T}^{emb} (GeV/#it{c})");
- fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistEmbPartPt[i]);
+ fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+ fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
+ fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
+ fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
+ fOutput->Add(fHistEmbPartPtvsJetPt[i]);
- histname = "fHistDistEmbPartJetAxis_";
+ histname = "fHistEmbPartPtvsJetCorrPt_";
+ histname += i;
+ fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
+ fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
+ fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
+ fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
+ fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
+ fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
+
+ histname = "fHistJetPtvsJetCorrPt_";
histname += i;
- fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
- fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
- fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistDistEmbPartJetAxis[i]);
+ fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
+ fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
+ fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
+ fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
+ fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
+ fOutput->Add(fHistJetPtvsJetCorrPt[i]);
+
+ histname = "fHistDistLeadPart2JetAxis_";
+ histname += i;
+ fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
+ fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
+ fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistDistLeadPart2JetAxis[i]);
histname = "fHistEmbNotFoundPhiEta_";
histname += i;
- fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
+ fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
fOutput->Add(fHistEmbNotFoundPhiEta[i]);
+ histname = "fHistEmbNotFoundPt_";
+ histname += i;
+ fHistEmbNotFoundPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
+ fHistEmbNotFoundPt[i]->GetXaxis()->SetTitle("#it{p}_{T,const}^{emb} (GeV/#it{c})");
+ fHistEmbNotFoundPt[i]->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistEmbNotFoundPt[i]);
+
+ histname = "fHistEmbRejectedJetsPhiEta_";
+ histname += i;
+ fHistEmbRejectedJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+ fHistEmbRejectedJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
+ fHistEmbRejectedJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
+ fOutput->Add(fHistEmbRejectedJetsPhiEta[i]);
+
+ histname = "fHistEmbRejectedJetsPtArea_";
+ histname += i;
+ fHistEmbRejectedJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
+ fHistEmbRejectedJetsPtArea[i]->GetXaxis()->SetTitle("area");
+ fHistEmbRejectedJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
+ fOutput->Add(fHistEmbRejectedJetsPtArea[i]);
+
histname = "fHistEmbBkgArea_";
histname += i;
- fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+ fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
- fHistEmbBkgArea[i]->GetYaxis()->SetTitle("background of embedded track (GeV/#it{c})");
+ fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
fOutput->Add(fHistEmbBkgArea[i]);
histname = "fHistRhoVSEmbBkg_";
histname += i;
fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
- fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("background of embedded track (GeV/#it{c})");
+ fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
fOutput->Add(fHistRhoVSEmbBkg[i]);
histname = "fHistDeltaPtEmbArea_";
histname += i;
- fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins * 2, -fMaxBinPt, fMaxBinPt);
+ fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
+ 50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
fOutput->Add(fHistDeltaPtEmbArea[i]);
}
}
+ delete[] binsPt;
+ delete[] binsCorrPt;
+ delete[] binsArea;
+
PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
}
{
// Fill histograms.
- Int_t *sortedJets = GetSortedArray(fJets);
-
- AliEmcalJet* jet = 0;
-
- if (sortedJets && sortedJets[0] > 0)
- jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
-
// ************
// Random cones
// _________________________________
const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
-
- for (Int_t i = 0; i < fRCperEvent; i++) {
- // Simple random cones
- Float_t RCpt = 0;
- Float_t RCeta = 0;
- Float_t RCphi = 0;
- GetRandomCone(RCpt, RCeta, RCphi, 0);
- if (RCpt > 0) {
- fHistRCPhiEta->Fill(RCeta, RCphi);
- fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
-
- fHistRCPt[fCentBin]->Fill(RCpt);
- fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
- }
+ Float_t RCpt = 0;
+ Float_t RCeta = 0;
+ Float_t RCphi = 0;
- // Random cones far from leading jet
- RCpt = 0;
- RCeta = 0;
- RCphi = 0;
- GetRandomCone(RCpt, RCeta, RCphi, jet);
- if (RCpt > 0) {
- if (jet) {
- Float_t dphi = RCphi - jet->Phi();
- if (dphi > 4.8) dphi -= TMath::Pi() * 2;
- if (dphi < -1.6) dphi += TMath::Pi() * 2;
- fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
+ if (fTracks || fCaloClusters) {
+
+ for (Int_t i = 0; i < fRCperEvent; i++) {
+ // Simple random cones
+ RCpt = 0;
+ RCeta = 0;
+ RCphi = 0;
+ GetRandomCone(RCpt, RCeta, RCphi, 0);
+ if (RCpt > 0) {
+ fHistRCPhiEta->Fill(RCeta, RCphi);
+ fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
+
+ fHistRCPt[fCentBin]->Fill(RCpt);
+ fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+ }
+
+ if (fJets) {
+
+ // Random cones far from leading jet
+ static Int_t sortedJets[9999] = {-1};
+ GetSortedArray(sortedJets, fJets);
+
+ AliEmcalJet* jet = 0;
+
+ if (sortedJets[0] >= 0)
+ jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
+
+ RCpt = 0;
+ RCeta = 0;
+ RCphi = 0;
+ GetRandomCone(RCpt, RCeta, RCphi, jet);
+ if (RCpt > 0) {
+ if (jet) {
+ Float_t dphi = RCphi - jet->Phi();
+ if (dphi > 4.8) dphi -= TMath::Pi() * 2;
+ if (dphi < -1.6) dphi += TMath::Pi() * 2;
+ fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
+ }
+ fHistRCPtExLJ[fCentBin]->Fill(RCpt);
+ fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
+ }
}
- fHistRCPtExLJ[fCentBin]->Fill(RCpt);
- fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
}
-
- // Random cones with randomized particles
+ }
+
+ // Random cones with randomized particles
+ if (fRandTracks || fRandCaloClusters) {
RCpt = 0;
RCeta = 0;
RCphi = 0;
fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
}
}
-
+
// ************
// Embedding
// _________________________________
- if (!fEmbJets)
- return kTRUE;
-
- AliEmcalJet *embJet = 0;
- TObject *embPart = 0;
-
- DoEmbJetLoop(embJet, embPart);
-
- if (embJet) {
- Double_t probePt = 0;
- Double_t probeEta = 0;
- Double_t probePhi = 0;
-
- AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
- if (cluster) {
- TLorentzVector nPart;
- cluster->GetMomentum(nPart, fVertex);
+ if (fEmbJets) {
+
+ AliEmcalJet *embJet = NextEmbeddedJet(0);
+
+ Int_t countEmbJets = 0;
+
+ while (embJet != 0) {
+ AliDebug(2,Form("Elaborating embedded jet n. %d", countEmbJets));
+ countEmbJets++;
+
+ if (!AcceptJet(embJet)) {
+ AliDebug(2,"Embedded jet not accepted, skipping...");
+ fHistEmbRejectedJetsPhiEta[fCentBin]->Fill(embJet->Eta(), embJet->Phi());
+ fHistEmbRejectedJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt());
+
+ embJet = NextEmbeddedJet();
+ continue;
+ }
+
+ Double_t maxClusterPt = 0;
+ Double_t maxClusterEta = 0;
+ Double_t maxClusterPhi = 0;
- probeEta = nPart.Eta();
- probePhi = nPart.Phi();
- probePt = nPart.Pt();
- }
- else {
- AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
- if (track) {
- probeEta = track->Eta();
- probePhi = track->Phi();
- probePt = track->Pt();
+ Double_t maxTrackPt = 0;
+ Double_t maxTrackEta = 0;
+ Double_t maxTrackPhi = 0;
+
+ Double_t maxPartPt = 0;
+ Double_t maxPartEta = 0;
+ Double_t maxPartPhi = 0;
+
+ if (fLeadingHadronType == 1 || fLeadingHadronType == 2) {
+ AliVCluster *cluster = embJet->GetLeadingCluster(fEmbCaloClusters);
+ if (cluster) {
+ TLorentzVector nPart;
+ cluster->GetMomentum(nPart, fVertex);
+
+ maxClusterEta = nPart.Eta();
+ maxClusterPhi = nPart.Phi();
+ maxClusterPt = nPart.Pt();
+ }
+ }
+
+ if (fLeadingHadronType == 0 || fLeadingHadronType == 2) {
+ AliVParticle *track = embJet->GetLeadingTrack(fEmbTracks);
+ if (track) {
+ maxTrackEta = track->Eta();
+ maxTrackPhi = track->Phi();
+ maxTrackPt = track->Pt();
+ }
+ }
+
+ if (maxTrackPt > maxClusterPt) {
+ maxPartPt = maxTrackPt;
+ maxPartEta = maxTrackEta;
+ maxPartPhi = maxTrackPhi;
}
else {
- AliWarning(Form("%s - Embedded jet found but embedded particle not found (wrong type?)!", GetName()));
- return kTRUE;
+ maxPartPt = maxClusterPt;
+ maxPartEta = maxClusterEta;
+ maxPartPhi = maxClusterPhi;
}
- }
- Double_t distProbeJet = TMath::Sqrt((embJet->Eta() - probeEta) * (embJet->Eta() - probeEta) + (embJet->Phi() - probePhi) * (embJet->Phi() - probePhi));
-
- fHistEmbPartPt[fCentBin]->Fill(probePt);
- fHistEmbPartPhiEta->Fill(probeEta, probePhi);
-
- fHistDistEmbPartJetAxis[fCentBin]->Fill(distProbeJet);
-
- fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt());
- fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area());
- fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
-
- fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - probePt);
- fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - probePt);
- fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - probePt);
+
+ Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - maxPartEta) * (embJet->Eta() - maxPartEta) + (embJet->Phi() - maxPartPhi) * (embJet->Phi() - maxPartPhi));
+
+ fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
+ fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
+ fHistLeadPartPhiEta->Fill(maxPartEta, maxPartPhi);
+ fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
+
+ fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), maxPartPt);
+ fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), maxPartPt);
+ fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
+ fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
+
+ fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
+ fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
+ fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
- }
- else {
- if (fEmbTracks)
- DoEmbTrackLoop();
- if (fEmbCaloClusters)
- DoEmbClusterLoop();
- if (fEmbTracks && fEmbeddedTrackId >= 0) {
- AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
- fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
- }
- else if (fEmbCaloClusters && fEmbeddedClusterId >= 0) {
- AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
- TLorentzVector nPart;
- cluster2->GetMomentum(nPart, fVertex);
- fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
+ embJet = NextEmbeddedJet();
}
- else {
- AliWarning(Form("%s - Embedded particle not found!", GetName()));
+
+ if (countEmbJets==0) {
+ AliDebug(1,"No embedded jets found!");
+ if (fEmbTracks) {
+ DoEmbTrackLoop();
+ for (Int_t i = 0; i < fEmbeddedTrackNIds; i++) {
+ AliDebug(2,Form("Embedded track %d found!",i));
+ AliVParticle *track2 = static_cast<AliVParticle*>(fEmbTracks->At(fEmbeddedTrackIds[i]));
+ if (!track2) continue;
+ fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
+ fHistEmbNotFoundPt[fCentBin]->Fill(track2->Pt());
+ }
+ }
+
+ if (fEmbCaloClusters) {
+ DoEmbClusterLoop();
+ for (Int_t i = 0; i < fEmbeddedClusterNIds; i++) {
+ AliDebug(2,Form("Embedded cluster %d found!",i));
+ AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterIds[i]));
+ TLorentzVector nPart;
+ cluster2->GetMomentum(nPart, fVertex);
+ fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
+ fHistEmbNotFoundPt[fCentBin]->Fill(nPart.Pt());
+ }
+ }
}
}
{
// Do track loop.
+ fEmbeddedTrackNIds = 0;
+
if (!fEmbTracks)
return;
- fEmbeddedTrackId = -1;
-
Int_t ntracks = fEmbTracks->GetEntriesFast();
for (Int_t i = 0; i < ntracks; i++) {
AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
- if (vtrack && !AcceptTrack(vtrack, kTRUE))
+ if (vtrack && !AcceptTrack(vtrack))
continue;
- if (track->GetLabel() == 100) {
- fEmbeddedTrackId = i;
- continue;
+ if (TMath::Abs(track->GetLabel()) > fMinMCLabel) {
+ if (fEmbeddedTrackNIds >= 999) {
+ AliWarning("The number of embedded tracks exceds 999!");
+ break;
+ }
+ fEmbeddedTrackIds[fEmbeddedTrackNIds] = i;
+ fEmbeddedTrackNIds++;
}
}
}
{
// Do cluster loop.
+ fEmbeddedClusterNIds = 0;
+
if (!fEmbCaloClusters)
return;
- fEmbeddedClusterId = -1;
-
Int_t nclusters = fEmbCaloClusters->GetEntriesFast();
for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
continue;
}
- if (!AcceptCluster(cluster, kTRUE))
+ if (!AcceptCluster(cluster))
continue;
- if (cluster->Chi2() == 100) {
- fEmbeddedClusterId = iClusters;
- continue;
+ if (cluster->GetLabel() > fMinMCLabel) {
+ if (fEmbeddedClusterNIds >= 999) {
+ AliWarning("The number of embedded clusters exceds 999!");
+ break;
+ }
+ fEmbeddedClusterIds[fEmbeddedClusterNIds] = iClusters;
+ fEmbeddedClusterNIds++;
}
}
}
//________________________________________________________________________
-void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
+AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Int_t i)
{
// Do the embedded jet loop.
+ static Int_t iJet = 0;
+
+ if (i >= 0)
+ iJet = i;
+ else
+ iJet++;
+
if (!fEmbJets)
- return;
+ return 0;
TLorentzVector maxClusVect;
- Int_t nembjets = fEmbJets->GetEntriesFast();
+ const Int_t nembjets = fEmbJets->GetEntriesFast();
- for (Int_t ij = 0; ij < nembjets; ij++) {
+ for (; iJet < nembjets; iJet++) {
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
+ AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(iJet));
if (!jet) {
- AliError(Form("Could not receive jet %d", ij));
+ AliError(Form("Could not receive jet %d", iJet));
continue;
}
-
- if (!AcceptJet(jet))
- continue;
- if (!jet->IsMC())
+ if (jet->MCPt() < fMCJetPtThreshold)
continue;
-
- AliVParticle *maxTrack = 0;
-
- for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
- AliVParticle *track = jet->TrackAt(it, fEmbTracks);
-
- if (!track)
- continue;
-
- if (track->GetLabel() != 100)
- continue;
-
- if (!maxTrack || track->Pt() > maxTrack->Pt())
- maxTrack = track;
- }
-
- AliVCluster *maxClus = 0;
-
- for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
- AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
-
- if (!cluster)
- continue;
-
- if (cluster->Chi2() != 100)
- continue;
-
- TLorentzVector nPart;
- cluster->GetMomentum(nPart, fVertex);
-
- if (!maxClus || nPart.Et() > maxClusVect.Et()) {
- new (&maxClusVect) TLorentzVector(nPart);
- maxClus = cluster;
- }
- }
-
- if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
- embPart = maxClus;
- embJet = jet;
- }
- else if (maxTrack) {
- embPart = maxTrack;
- embJet = jet;
- }
-
- return; // MC jet found, exit
+
+ return jet;
}
+
+ return 0;
}
//________________________________________________________________________
continue;
}
- if (!AcceptCluster(cluster, fMCAna))
+ if (!AcceptCluster(cluster))
continue;
TLorentzVector nPart;
AliError(Form("Could not retrieve track %d",iTracks));
continue;
}
-
- if (!AcceptTrack(track, fMCAna))
+
+ if (!AcceptTrack(track))
continue;
Float_t tracketa = track->Eta();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff