#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),
- fTracks2Name(""),
- fCalo2Name(""),
- fJets2Name(""),
- fRho2Name(""),
- fJet2Radius(0.4),
- fJet2AreaCut(-1),
- fPtBiasJet2Track(0),
- fPtBiasJet2Clus(0),
- fJet2MinEta(-999),
- fJet2MaxEta(-999),
- fJet2MinPhi(-999),
- fJet2MaxPhi(-999),
- fMaxClusterPt2(1000),
- fMaxTrackPt2(1000),
- fAreCollections1MC(kFALSE),
- fAreCollections2MC(kTRUE),
fMatching(kNoMatching),
fMatchingPar1(0),
fMatchingPar2(0),
fDeltaEtaDeltaPhiAxis(0),
fNEFAxis(0),
fZAxis(0),
- fDoJet2Histogram(0),
- fTracks2(0),
- fCaloClusters2(0),
- fJets2(0),
- fRho2(0),
- fRho2Val(0),
- fTracks2Map(0),
- fHistLeadingJets1PtArea(0),
- fHistLeadingJets1CorrPtArea(0),
- fHistLeadingJets2PtArea(0),
- fHistLeadingJets2CorrPtArea(0),
- fHistLeadingJets2PtAreaAcceptance(0),
- fHistLeadingJets2CorrPtAreaAcceptance(0),
+ fIsJet1Rho(kFALSE),
+ fIsJet2Rho(kFALSE),
+ fHistRejectionReason1(0),
+ fHistRejectionReason2(0),
fHistJets1(0),
fHistJets2(0),
- fHistJets2Acceptance(0),
fHistMatching(0),
fHistJets1PhiEta(0),
fHistJets1PtArea(0),
fHistJets2PhiEta(0),
fHistJets2PtArea(0),
fHistJets2CorrPtArea(0),
- fHistJets2PhiEtaAcceptance(0),
- fHistJets2PtAreaAcceptance(0),
- fHistJets2CorrPtAreaAcceptance(0),
fHistJets2NEFvsPt(0),
fHistJets2CEFvsCEFPt(0),
fHistJets2ZvsPt(0),
//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker(const char *name) :
AliAnalysisTaskEmcalJet(name, kTRUE),
- fTracks2Name("MCParticles"),
- fCalo2Name(""),
- fJets2Name("MCJets"),
- fRho2Name(""),
- fJet2Radius(0.4),
- fJet2AreaCut(-1),
- fPtBiasJet2Track(0),
- fPtBiasJet2Clus(0),
- fJet2MinEta(-999),
- fJet2MaxEta(-999),
- fJet2MinPhi(-999),
- fJet2MaxPhi(-999),
- fMaxClusterPt2(1000),
- fMaxTrackPt2(1000),
- fAreCollections1MC(kFALSE),
- fAreCollections2MC(kTRUE),
fMatching(kNoMatching),
fMatchingPar1(0),
fMatchingPar2(0),
fDeltaEtaDeltaPhiAxis(0),
fNEFAxis(0),
fZAxis(0),
- fDoJet2Histogram(0),
- fTracks2(0),
- fCaloClusters2(0),
- fJets2(0),
- fRho2(0),
- fRho2Val(0),
- fTracks2Map(0),
- fHistLeadingJets1PtArea(0),
- fHistLeadingJets1CorrPtArea(0),
- fHistLeadingJets2PtArea(0),
- fHistLeadingJets2CorrPtArea(0),
- fHistLeadingJets2PtAreaAcceptance(0),
- fHistLeadingJets2CorrPtAreaAcceptance(0),
+ fIsJet1Rho(kFALSE),
+ fIsJet2Rho(kFALSE),
+ fHistRejectionReason1(0),
+ fHistRejectionReason2(0),
fHistJets1(0),
fHistJets2(0),
- fHistJets2Acceptance(0),
fHistMatching(0),
fHistJets1PhiEta(0),
fHistJets1PtArea(0),
fHistJets2PhiEta(0),
fHistJets2PtArea(0),
fHistJets2CorrPtArea(0),
- fHistJets2PhiEtaAcceptance(0),
- fHistJets2PtAreaAcceptance(0),
- fHistJets2CorrPtAreaAcceptance(0),
fHistJets2NEFvsPt(0),
fHistJets2CEFvsCEFPt(0),
fHistJets2ZvsPt(0),
fHistJets1PtArea->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistJets1PtArea);
- if (!fRhoName.IsNull()) {
+ if (fIsJet1Rho) {
fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea",
fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
fHistJets2PtArea->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistJets2PtArea);
- if (!fRho2Name.IsNull()) {
+ if (fIsJet2Rho) {
fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea",
fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
fOutput->Add(fHistJets2CorrPtArea);
}
- fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
- fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
- fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistJets2PhiEtaAcceptance);
-
- fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance",
- fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
- fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistJets2PtAreaAcceptance);
-
- if (!fRho2Name.IsNull()) {
- fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance",
- fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
- fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistJets2CorrPtAreaAcceptance);
- }
-
fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistJet1PtvsJet2Pt);
- if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
- if (fRhoName.IsNull())
+ if (fIsJet1Rho || fIsJet2Rho) {
+ if (!fIsJet1Rho)
fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
else
fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt);
- if (fRho2Name.IsNull())
+ if (!fIsJet2Rho)
fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
else
fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
- if (fRhoName.IsNull())
+ if (!fIsJet1Rho)
fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- else if (fRho2Name.IsNull())
+ else if (!fIsJet2Rho)
fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
else
{
// Allocate THnSparse histograms.
- TString title[20]= {""};
- Int_t nbins[20] = {0};
- Double_t min[20] = {0.};
- Double_t max[20] = {0.};
+ 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";
max[dim] = 1.5;
dim++;
- title[dim] = "NEF";
- nbins[dim] = fNbins/4;
- min[dim] = 0;
- max[dim] = 1.2;
- dim++;
+ if (fNEFAxis) {
+ title[dim] = "NEF";
+ nbins[dim] = fNbins/4;
+ min[dim] = 0;
+ max[dim] = 1.2;
+ dim++;
+ }
- title[dim] = "Z";
- nbins[dim] = fNbins/4;
+ if (fZAxis) {
+ title[dim] = "Z";
+ nbins[dim] = fNbins/4;
+ min[dim] = 0;
+ max[dim] = 1.2;
+ dim++;
+ }
+
+ title[dim] = "p_{T,particle}^{leading} (GeV/c)";
+ nbins[dim] = 120;
min[dim] = 0;
- max[dim] = 1.2;
+ max[dim] = 120;
dim++;
Int_t dim1 = dim, dim2 = dim;
- if (!fRhoName.IsNull()) {
+ if (fIsJet1Rho) {
title[dim1] = "p_{T}^{corr}";
nbins[dim1] = fNbins*2;
min[dim1] = -250;
fHistJets1->GetAxis(i)->SetTitle(title[i]);
fOutput->Add(fHistJets1);
- if (!fRho2Name.IsNull()) {
+ if (fIsJet2Rho) {
title[dim2] = "p_{T}^{corr}";
nbins[dim2] = fNbins*2;
min[dim2] = -250;
dim2++;
}
- if (fDoJet2Histogram) {
- 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);
- }
-
- fHistJets2Acceptance = new THnSparseD("fHistJets2Acceptance","fHistJets2Acceptance",dim2,nbins,min,max);
+ fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
for (Int_t i = 0; i < dim2; i++)
- fHistJets2Acceptance->GetAxis(i)->SetTitle(title[i]);
- fOutput->Add(fHistJets2Acceptance);
+ fHistJets2->GetAxis(i)->SetTitle(title[i]);
+ fOutput->Add(fHistJets2);
// Matching
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;
max[dim] = 250;
dim++;
}
- if (!fRhoName.IsNull()) {
+ if (fIsJet1Rho) {
title[dim] = "p_{T,1}^{corr}";
nbins[dim] = fNbins*2;
min[dim] = -250;
max[dim] = 250;
dim++;
}
- if (!fRho2Name.IsNull()) {
+ if (fIsJet2Rho) {
title[dim] = "p_{T,2}^{corr}";
nbins[dim] = fNbins*2;
min[dim] = -250;
max[dim] = 250;
dim++;
}
- if (fDeltaPtAxis && (!fRhoName.IsNull() || !fRho2Name.IsNull())) {
+ if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
title[dim] = "#deltap_{T}^{corr}";
nbins[dim] = fNbins*2;
min[dim] = -250;
AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
- // Jets 1 histograms
- fHistLeadingJets1PtArea = new TH2F("fHistLeadingJets1PtArea", "fHistLeadingJets1PtArea",
- fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
- fHistLeadingJets1PtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
- fHistLeadingJets1PtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets1PtArea);
-
- if (!fRhoName.IsNull()) {
- fHistLeadingJets1CorrPtArea = new TH2F("fHistLeadingJets1CorrPtArea", "fHistLeadingJets1CorrPtArea",
- fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistLeadingJets1CorrPtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
- fHistLeadingJets1CorrPtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets1CorrPtArea);
- }
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- fHistLeadingJets2PtAreaAcceptance = new TH2F("fHistLeadingJets2PtAreaAcceptance", "fHistLeadingJets2PtAreaAcceptance",
- fNbins/2, 1, 2, fNbins, fMinBinPt, fMaxBinPt);
- fHistLeadingJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistLeadingJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistLeadingJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2PtAreaAcceptance);
-
- fHistLeadingJets2PtArea = new TH2F("fHistLeadingJets2PtArea", "fHistLeadingJets2PtArea",
- fNbins/2, 0, 1, fNbins, fMinBinPt, fMaxBinPt);
- fHistLeadingJets2PtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistLeadingJets2PtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2PtArea);
-
- if (!fRho2Name.IsNull()) {
- fHistLeadingJets2CorrPtAreaAcceptance = new TH2F("fHistLeadingJets2CorrPtAreaAcceptance", "fHistLeadingJets2CorrPtAreaAcceptance",
- fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistLeadingJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistLeadingJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
- fHistLeadingJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2CorrPtAreaAcceptance);
-
- fHistLeadingJets2CorrPtArea = new TH2F("fHistLeadingJets2CorrPtArea", "fHistLeadingJets2CorrPtArea",
- fNbins/2, 0, 1, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistLeadingJets2CorrPtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
- fHistLeadingJets2CorrPtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2CorrPtArea);
- }
+ if (!jets1 || !jets2) return;
+
+ if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
+ else fIsJet1Rho = kTRUE;
+ if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
+ else fIsJet2Rho = kTRUE;
+
+ fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250);
+ fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason");
+ fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+ fHistRejectionReason1->GetZaxis()->SetTitle("counts");
+ SetRejectionReasonLabels(fHistRejectionReason1->GetXaxis());
+ fOutput->Add(fHistRejectionReason1);
+
+ fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250);
+ fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason");
+ fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
+ fHistRejectionReason2->GetZaxis()->SetTitle("counts");
+ SetRejectionReasonLabels(fHistRejectionReason2->GetXaxis());
+ fOutput->Add(fHistRejectionReason2);
if (fHistoType==0)
AllocateTH2();
}
//________________________________________________________________________
-void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t Z, Double_t CorrPt, Double_t MCPt, Int_t Set)
+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)
{
if (fHistoType==1) {
THnSparse *histo = 0;
histo = fHistJets1;
else if (Set==2)
histo = fHistJets2;
- else if (Set==3)
- histo = fHistJets2Acceptance;
if (!histo)
return;
else if (title=="NEF")
contents[i] = NEF;
else if (title=="Z")
- contents[i] = 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()));
}
if (Set == 1) {
fHistJets1PtArea->Fill(A, Pt);
fHistJets1PhiEta->Fill(Eta, Phi);
- fHistJets1ZvsPt->Fill(Z, Pt);
+ fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt);
fHistJets1NEFvsPt->Fill(NEF, Pt);
fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
- if (!fRhoName.IsNull())
+ if (fIsJet1Rho)
fHistJets1CorrPtArea->Fill(A, CorrPt);
}
else if (Set == 2) {
fHistJets2PtArea->Fill(A, Pt);
fHistJets2PhiEta->Fill(Eta, Phi);
- if (!fRho2Name.IsNull())
- fHistJets2CorrPtArea->Fill(A, CorrPt);
- }
- else if (Set == 3) {
- fHistJets2PtAreaAcceptance->Fill(A, Pt);
- fHistJets2PhiEtaAcceptance->Fill(Eta, Phi);
- fHistJets2ZvsPt->Fill(Z, Pt);
+ fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt);
fHistJets2NEFvsPt->Fill(NEF, Pt);
fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
- if (!fRho2Name.IsNull())
- fHistJets2CorrPtAreaAcceptance->Fill(A, CorrPt);
+ if (fIsJet2Rho)
+ fHistJets2CorrPtArea->Fill(A, CorrPt);
}
}
}
//________________________________________________________________________
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 Z1, Double_t Z2)
+ Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
{
if (fHistoType==1) {
Double_t contents[20]={0};
else if (title=="NEF_{2}")
contents[i] = NEF2;
else if (title=="Z_{1}")
- contents[i] = Z1;
+ contents[i] = LeadingPt1/Pt1;
else if (title=="Z_{2}")
- contents[i] = Z2;
+ 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()));
}
fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
- if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
+ if (fIsJet1Rho || fIsJet2Rho) {
Double_t dcorrpt = CorrPt1 - CorrPt2;
fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
}
}
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
-{
- // Return true if jet is accepted.
-
- if (jet->Pt() < fJetPtCut)
- return kFALSE;
- if (jet->Area() < fJetAreaCut)
- return kFALSE;
- if (jet->MCPt() < fMinJetMCPt)
- return kFALSE;
-
- return kTRUE;
-}
-
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
-{
- // Accept jet with a bias.
-
- if (fLeadingHadronType == 0) {
- if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
- }
- else if (fLeadingHadronType == 1) {
- if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
- }
- else {
- if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
- }
-
- return kTRUE;
-}
-
//________________________________________________________________________
void AliJetResponseMaker::ExecOnce()
{
// Execute once.
- if (!fJets2Name.IsNull() && !fJets2) {
- fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
- if (!fJets2) {
- AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
- return;
- }
- else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
- AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
- fJets2 = 0;
- return;
- }
- }
-
- if (!fTracks2Name.IsNull() && !fTracks2) {
- fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
- if (!fTracks2) {
- AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
- return;
- }
- else {
- TClass *cl = fTracks2->GetClass();
- if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
- AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
- fTracks2 = 0;
- return;
- }
- }
-
- if (fAreCollections2MC) {
- fTracks2Map = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
- // this is needed to map the MC labels with the indexes of the MC particle collection
- // if teh map is not given, the MC labels are assumed to be consistent with the indexes (which is not the case if AliEmcalMCTrackSelector is used)
- if (!fTracks2Map) {
- AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data()));
- fTracks2Map = new AliNamedArrayI("tracksMap",9999);
- for (Int_t i = 0; i < 9999; i++) {
- fTracks2Map->AddAt(i,i);
- }
- }
- }
- }
+ AliAnalysisTaskEmcalJet::ExecOnce();
- if (!fCalo2Name.IsNull() && !fCaloClusters2) {
- fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
- if (!fCaloClusters2) {
- AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
- return;
- } else {
- TClass *cl = fCaloClusters2->GetClass();
- if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
- AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
- fCaloClusters2 = 0;
- return;
- }
- }
- }
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- if (!fRho2Name.IsNull() && !fRho2) {
- fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
- if (!fRho2) {
- AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
- fInitialized = kFALSE;
- return;
- }
- }
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
- if (fPercAreaCut >= 0) {
- if (fJet2AreaCut >= 0)
- AliInfo(Form("%s: jet 2 area cut will be calculated as a percentage of the average area, given value will be overwritten", GetName()));
- fJet2AreaCut = fPercAreaCut * fJet2Radius * fJet2Radius * TMath::Pi();
- }
- if (fJet2AreaCut < 0)
- fJet2AreaCut = 0;
-
- if (fJet2MinEta == -999)
- fJet2MinEta = fJetMinEta - fJetRadius;
- if (fJet2MaxEta == -999)
- fJet2MaxEta = fJetMaxEta + fJetRadius;
- if (fJet2MinPhi == -999)
- fJet2MinPhi = fJetMinPhi - fJetRadius;
- if (fJet2MaxPhi == -999)
- fJet2MaxPhi = fJetMaxPhi + fJetRadius;
-
- if (fMatching == kMCLabel && (!fAreCollections2MC || fAreCollections1MC)) {
- if (fAreCollections2MC == fAreCollections1MC) {
+ if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
+ if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
AliWarning("Changing matching type from MC label to same collection...");
fMatching = kSameCollections;
}
fMatching = kGeometrical;
}
}
- else if (fMatching == kSameCollections && fAreCollections2MC != fAreCollections1MC) {
- if (fAreCollections2MC && !fAreCollections1MC) {
+ 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;
}
fMatching = kGeometrical;
}
}
-
- AliAnalysisTaskEmcalJet::ExecOnce();
-}
-
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::RetrieveEventObjects()
-{
- // Retrieve event objects.
-
- if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
- return kFALSE;
-
- if (fRho2)
- fRho2Val = fRho2->GetVal();
-
- return kTRUE;
}
//________________________________________________________________________
//________________________________________________________________________
Bool_t AliJetResponseMaker::DoJetMatching()
{
- DoJetLoop(kFALSE);
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- const Int_t nJets = fJets->GetEntriesFast();
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
- for (Int_t i = 0; i < nJets; i++) {
+ DoJetLoop();
- AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
+ AliEmcalJet* jet1 = 0;
- if (!jet1) {
- AliError(Form("Could not receive jet %d", i));
- continue;
- }
+ jets1->ResetCurrentID();
+ while ((jet1 = jets1->GetNextJet())) {
- if (!AcceptJet(jet1))
- continue;
+ AliEmcalJet *jet2 = jet1->ClosestJet();
- if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
- continue;
+ if (!jet2) continue;
+ if (jet2->ClosestJet() != jet1) continue;
+ if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
- if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 &&
- jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) { // Matched jet found
- jet1->SetMatchedToClosest(fMatching);
- jet1->ClosestJet()->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(),
- jet1->MatchedJet()->Pt(), jet1->MatchedJet()->Eta(), jet1->MatchedJet()->Phi()));
- }
+ // 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(Bool_t order)
+void AliJetResponseMaker::DoJetLoop()
{
// Do the jet loop.
- TClonesArray *jets1 = 0;
- TClonesArray *jets2 = 0;
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- if (order) {
- jets1 = fJets2;
- jets2 = fJets;
- }
- else {
- jets1 = fJets;
- jets2 = fJets2;
- }
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
- Int_t nJets1 = jets1->GetEntriesFast();
- Int_t nJets2 = jets2->GetEntriesFast();
+ AliEmcalJet* jet1 = 0;
+ AliEmcalJet* jet2 = 0;
- 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));
- continue;
- }
-
- jet2->ResetMatching();
-
-
- if (!AcceptJet(jet2))
- continue;
-
- if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
- continue;
- }
+ jets2->ResetCurrentID();
+ while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
- for (Int_t i = 0; i < nJets1; i++) {
-
- AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
-
- if (!jet1) {
- AliError(Form("Could not receive jet %d", i));
- continue;
- }
-
+ jets1->ResetCurrentID();
+ while ((jet1 = jets1->GetNextJet())) {
jet1->ResetMatching();
+
+ if (jet1->MCPt() < fMinJetMCPt) continue;
- if (!AcceptJet(jet1))
- continue;
-
- if (order) {
- if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
- continue;
- }
- else {
- if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
- continue;
- }
-
- 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));
- continue;
- }
- if (!AcceptJet(jet2))
- continue;
-
- if (order) {
- if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
- continue;
- }
- else {
- if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
- continue;
- }
-
+ jets2->ResetCurrentID();
+ while ((jet2 = jets2->GetNextJet())) {
SetMatchingLevel(jet1, jet2, fMatching);
} // jet2 loop
-
} // jet1 loop
}
//________________________________________________________________________
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
- for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
- Bool_t track2Found = kFALSE;
- Int_t index2 = jet2->TrackAt(iTrack2);
+ // 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,fTracks);
+ AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
if (!track) {
AliWarning(Form("Could not find track %d!", iTrack));
continue;
}
+
Int_t MClabel = TMath::Abs(track->GetLabel());
- if (MClabel > fMCLabelShift)
- MClabel -= fMCLabelShift;
- Int_t index = -1;
-
- if (MClabel == 0) {// 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();
+ MClabel -= fMCLabelShift;
+ if (MClabel != 0) continue;
+
+ // this is not a MC particle; remove it completely
+ AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
+ totalPt1 -= track->Pt();
+ d1 -= track->Pt();
+ }
+ }
+
+ // 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;
}
- else if (MClabel < fTracks2Map->GetSize()) {
- index = fTracks2Map->At(MClabel);
- }
+ TLorentzVector part;
+ clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
- if (index < 0) {
- AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
- continue;
- }
+ for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
+ Int_t cellId = clus->GetCellAbsId(iCell);
+ Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
- if (index2 == index) { // found common particle
- track2Found = kTRUE;
- d1 -= track->Pt();
- AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(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();
- break;
+ Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
+ MClabel -= fMCLabelShift;
+ if (MClabel != 0) continue;
+
+ // this is not a MC particle; remove it completely
+ AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
+ totalPt1 -= part.Pt() * cellFrac;
+ d1 -= part.Pt() * cellFrac;
}
}
+ }
+ else {
for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
- AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
+ 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));
- if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
+ 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 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;
+ }
+ Int_t MClabel = TMath::Abs(track->GetLabel());
+ MClabel -= fMCLabelShift;
+ if (MClabel <= 0) continue;
+
+ Int_t index = -1;
+ index = tracks2->GetIndexFromLabel(MClabel);
+ if (index < 0) {
+ AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
+ continue;
+ }
+
+ if (index2 != index) continue;
+
+ // found common particle
+ d1 -= track->Pt();
+
+ if (!track2Found) {
+ AliVParticle *MCpart = tracks2->GetParticle(index2);
+ AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+ iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
+ d2 -= MCpart->Pt();
+ }
+
+ track2Found = kTRUE;
+ }
+
+ // 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));
- if (MClabel > fMCLabelShift)
- MClabel -= fMCLabelShift;
- Int_t index = -1;
+ MClabel -= fMCLabelShift;
+ if (MClabel <= 0) continue;
- if (MClabel == 0) {// 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;
- continue;
- }
- else if (MClabel < fTracks2Map->GetSize()) {
- index = fTracks2Map->At(MClabel);
- }
-
- if (index < 0) {
+ 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 == index) { // 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 = static_cast<AliVParticle*>(fTracks2->At(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;
- }
- break;
+
+ 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 { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
- Int_t MClabel = TMath::Abs(clus->GetLabel());
- if (MClabel > fMCLabelShift)
- MClabel -= fMCLabelShift;
- Int_t index = -1;
-
- if (MClabel == 0) {// 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();
+ }
+ 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;
}
- else if (MClabel < fTracks2Map->GetSize()) {
- index = fTracks2Map->At(MClabel);
- }
-
+ TLorentzVector part;
+ clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+
+ Int_t MClabel = TMath::Abs(clus->GetLabel());
+ MClabel -= fMCLabelShift;
+ if (MClabel <= 0) continue;
+
+ 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) { // 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 = static_cast<AliVParticle*>(fTracks2->At(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();
- }
- break;
+
+ 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;
}
}
}
//________________________________________________________________________
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 (fTracks && fTracks2) {
+ if (tracks1 && tracks2) {
for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
Int_t index2 = jet2->TrackAt(iTrack2);
- for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
- Int_t index = jet1->TrackAt(iTrack);
- if (index2 == index) { // found common particle
- AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
- if (!part) {
- AliWarning(Form("Could not find track %d!", index));
+ 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 = static_cast<AliVParticle*>(fTracks2->At(index2));
+ AliVParticle *part2 = tracks2->GetParticle(index2);
if (!part2) {
AliWarning(Form("Could not find track %d!", index2));
continue;
}
- d1 -= part->Pt();
+ d1 -= part1->Pt();
d2 -= part2->Pt();
break;
}
}
- if (fCaloClusters && fCaloClusters2) {
+ if (clusters1 && clusters2) {
if (fUseCellsToMatch) {
const Int_t nClus1 = jet1->GetNumberOfClusters();
Double_t ptClus1[nClus1];
for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
Int_t index1 = jet1->ClusterAt(iClus1);
- AliVCluster *clus1 = static_cast<AliVCluster*>(fCaloClusters->At(index1));
+ AliVCluster *clus1 = clusters1->GetCluster(index1);
if (!clus1) {
AliWarning(Form("Could not find cluster %d!", index1));
ncells1[iClus1] = 0;
Int_t sortedIndexes2[maxNcells2];
for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
Int_t index2 = jet2->ClusterAt(iClus2);
- AliVCluster *clus2 = static_cast<AliVCluster*>(fCaloClusters2->At(index2));
+ AliVCluster *clus2 = clusters2->GetCluster(index2);
if (!clus2) {
AliWarning(Form("Could not find cluster %d!", index2));
continue;
else {
for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
Int_t index2 = jet2->ClusterAt(iClus2);
- for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
- Int_t index = jet1->ClusterAt(iClus);
- if (index2 == index) { // found common particle
- AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(index));
- if (!clus) {
- AliWarning(Form("Could not find cluster %d!", index));
+ 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 = static_cast<AliVCluster*>(fCaloClusters2->At(index2));
+ AliVCluster *clus2 = clusters2->GetCluster(index2);
if (!clus2) {
AliWarning(Form("Could not find cluster %d!", index2));
continue;
}
- TLorentzVector part, part2;
- clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+ TLorentzVector part1, part2;
+ clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
- d1 -= part.Pt();
+ d1 -= part1.Pt();
d2 -= part2.Pt();
break;
}
{
// Fill histograms.
- static Int_t indexes[9999] = {-1};
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- const Int_t nJets2 = GetSortedArray(indexes, fJets2, fRho2Val);
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
- Int_t naccJets2 = 0;
- Int_t naccJets2Acceptance = 0;
-
- for (Int_t i = 0; i < nJets2; i++) {
- AliDebug(2,Form("Processing jet (2) %d", indexes[i]));
+ AliEmcalJet* jet1 = 0;
+ AliEmcalJet* jet2 = 0;
+
+ jets2->ResetCurrentID();
+ while ((jet2 = jets2->GetNextJet())) {
- AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(indexes[i]));
+ AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
- if (!jet2) {
- AliError(Form("Could not receive jet2 %d", i));
- continue;
- }
+ if (jet2->Pt() < jets2->GetJetPtCut()) continue;
- // Jet 2 cuts
- if (!AcceptJet(jet2))
- continue;
- if (!AcceptBiasJet2(jet2))
- continue;
- if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
- continue;
- if (jet2->MaxTrackPt() > fMaxTrackPt2 || jet2->MaxClusterPt() > fMaxClusterPt2)
- continue;
+ Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
+ Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
- if (naccJets2 < fNLeadingJets) {
- fHistLeadingJets2PtArea->Fill(jet2->Area(), jet2->Pt());
- if (!fRho2Name.IsNull())
- fHistLeadingJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
- }
+ if (jets2->AcceptJet(jet2))
+ FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2,
+ corrpt2, jet2->MCPt(), 2);
+ else
+ fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(), jet2->Pt());
- if (fDoJet2Histogram)
- FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), jet2->Pt() - fRho2Val * jet2->Area(), jet2->MCPt(), 2);
+ jet1 = jet2->MatchedJet();
- naccJets2++;
+ if (!jet1) continue;
+ if (!jets1->AcceptJet(jet1)) continue;
+ if (jet1->MCPt() < fMinJetMCPt) continue;
- // Verify also jet cuts 1 on jet 2
- if (AcceptBiasJet(jet2) &&
- (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
-
- if (naccJets2Acceptance < fNLeadingJets) {
- fHistLeadingJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
- if (!fRho2Name.IsNull())
- fHistLeadingJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
- }
+ 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);
- FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), jet2->Pt() - fRho2Val * jet2->Area(), jet2->MCPt(), 3);
- naccJets2Acceptance++;
+ d = jet2->ClosestJetDistance();
}
+ else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
+ GetGeometricalMatchingLevel(jet1, jet2, d);
- if (jet2->MatchedJet()) {
-
- if (AcceptJet(jet2->MatchedJet()) &&
- jet2->MatchedJet()->Eta() > fJetMinEta && jet2->MatchedJet()->Eta() < fJetMaxEta &&
- jet2->MatchedJet()->Phi() > fJetMinPhi && jet2->MatchedJet()->Phi() < fJetMaxPhi &&
- AcceptBiasJet(jet2->MatchedJet()) &&
- jet2->MatchedJet()->MaxTrackPt() < fMaxTrackPt && jet2->MatchedJet()->MaxClusterPt() < fMaxClusterPt) {
-
- Double_t d=-1, ce1=-1, ce2=-1;
- if (jet2->GetMatchingType() == kGeometrical) {
- if (fAreCollections2MC && !fAreCollections1MC) // the other way around is not supported
- GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, ce1, ce2);
- else if (fAreCollections1MC == fAreCollections2MC)
- GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, ce1, ce2);
-
- d = jet2->ClosestJetDistance();
- }
- else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
- GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d);
-
- ce1 = jet2->MatchedJet()->ClosestJetDistance();
- ce2 = jet2->ClosestJetDistance();
- }
-
- Double_t corrpt1 = jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area();
- Double_t corrpt2 = jet2->Pt() - fRho2Val * jet2->Area();
-
- FillMatchingHistos(jet2->MatchedJet()->Pt(), jet2->Pt(), jet2->MatchedJet()->Eta(), jet2->Eta(), jet2->MatchedJet()->Phi(), jet2->Phi(),
- jet2->MatchedJet()->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet2->MatchedJet()->MCPt(),
- jet2->MatchedJet()->NEF(), jet2->NEF(), jet2->MatchedJet()->MaxPartPt()/jet2->MatchedJet()->Pt(), jet2->MaxPartPt()/jet2->Pt());
- }
+ ce1 = jet1->ClosestJetDistance();
+ ce2 = jet2->ClosestJetDistance();
}
- }
-
- const Int_t nJets1 = GetSortedArray(indexes, fJets, fRhoVal);
-
- Int_t naccJets1 = 0;
- for (Int_t i = 0; i < nJets1; i++) {
- AliDebug(2,Form("Processing jet (1) %d", indexes[i]));
- AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(indexes[i]));
+ Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
+ Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
- if (!jet1) {
- AliError(Form("Could not receive jet %d", i));
- continue;
- }
+ 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);
+ }
- if (!AcceptJet(jet1))
- continue;
- if (!AcceptBiasJet(jet1))
- continue;
- if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
+ jets1->ResetCurrentID();
+ while ((jet1 = jets1->GetNextJet())) {
+ if (!jets1->AcceptJet(jet1)) {
+ fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(), jet1->Pt());
continue;
- if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
- continue;
-
- if (naccJets1 < fNLeadingJets){
- fHistLeadingJets1PtArea->Fill(jet1->Area(), jet1->Pt());
- if (!fRhoName.IsNull())
- fHistLeadingJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
}
+ if (jet1->MCPt() < fMinJetMCPt) continue;
+ AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
- FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), jet1->MaxPartPt()/jet1->Pt(), jet1->Pt() - fRhoVal * jet1->Area(), jet1->MCPt(), 1);
- naccJets1++;
- }
+ 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