// $Id$ // // Jet deltaPt task. // // Author: S.Aiola #include #include #include #include #include #include #include #include "AliVCluster.h" #include "AliVParticle.h" #include "AliVTrack.h" #include "AliEmcalJet.h" #include "AliRhoParameter.h" #include "AliLog.h" #include "AliJetContainer.h" #include "AliParticleContainer.h" #include "AliClusterContainer.h" #include "AliAnalysisTaskDeltaPt.h" ClassImp(AliAnalysisTaskDeltaPt) //________________________________________________________________________ AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() : AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE), fMCJetPtThreshold(1), fMinRC2LJ(-1), fRCperEvent(-1), fConeRadius(0.2), fConeMinEta(-0.9), fConeMaxEta(0.9), fConeMinPhi(0), fConeMaxPhi(TMath::Pi()*2), fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fEmbJetsCont(0), fEmbTracksCont(0), fEmbCaloClustersCont(0), fRandTracksCont(0), fRandCaloClustersCont(0), fHistRCPhiEta(0), fHistRCPt(0), fHistRCPtExLJ(0), fHistRCPtExPartialLJ(0), fHistRCPtRand(0), fHistRhoVSRCPt(0), fHistDeltaPtRCvsEP(0), fHistDeltaPtRCExLJ(0), fHistDeltaPtRCExPartialLJ(0), fHistDeltaPtRCRand(0), fHistEmbJetsPtArea(0), fHistEmbJetsCorrPtArea(0), fHistEmbPartPtvsJetPt(0), fHistEmbPartPtvsJetCorrPt(0), fHistJetPtvsJetCorrPt(0), fHistDistLeadPart2JetAxis(0), fHistEmbBkgArea(0), fHistRhoVSEmbBkg(0), fHistDeltaPtEmbArea(0), fHistDeltaPtEmbvsEP(0), fHistRCPtExLJVSDPhiLJ(0), fHistRCPtExPartialLJVSDPhiLJ(0), fHistEmbJetsPhiEta(0), fHistLeadPartPhiEta(0) { // Default constructor. fHistRCPt = 0; fHistRCPtExLJ = 0; fHistRCPtExPartialLJ = 0; fHistRCPtRand = 0; fHistRhoVSRCPt = 0; fHistDeltaPtRCvsEP = 0; fHistDeltaPtRCExLJ = 0; fHistDeltaPtRCExPartialLJ = 0; fHistDeltaPtRCRand = 0; fHistEmbJetsPtArea = 0; fHistEmbJetsCorrPtArea = 0; fHistEmbPartPtvsJetPt = 0; fHistEmbPartPtvsJetCorrPt = 0; fHistJetPtvsJetCorrPt = 0; fHistDistLeadPart2JetAxis = 0; fHistEmbBkgArea = 0; fHistRhoVSEmbBkg = 0; fHistDeltaPtEmbArea = 0; fHistDeltaPtEmbvsEP = 0; SetMakeGeneralHistograms(kTRUE); } //________________________________________________________________________ AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) : AliAnalysisTaskEmcalJet(name, kTRUE), fMCJetPtThreshold(1), fMinRC2LJ(-1), fRCperEvent(-1), fConeRadius(0.2), fConeMinEta(-0.9), fConeMaxEta(0.9), fConeMinPhi(0), fConeMaxPhi(TMath::Pi()*2), fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fEmbJetsCont(0), fEmbTracksCont(0), fEmbCaloClustersCont(0), fRandTracksCont(0), fRandCaloClustersCont(0), fHistRCPhiEta(0), fHistRCPt(0), fHistRCPtExLJ(0), fHistRCPtExPartialLJ(0), fHistRCPtRand(0), fHistRhoVSRCPt(0), fHistDeltaPtRCvsEP(0), fHistDeltaPtRCExLJ(0), fHistDeltaPtRCExPartialLJ(0), fHistDeltaPtRCRand(0), fHistEmbJetsPtArea(0), fHistEmbJetsCorrPtArea(0), fHistEmbPartPtvsJetPt(0), fHistEmbPartPtvsJetCorrPt(0), fHistJetPtvsJetCorrPt(0), fHistDistLeadPart2JetAxis(0), fHistEmbBkgArea(0), fHistRhoVSEmbBkg(0), fHistDeltaPtEmbArea(0), fHistDeltaPtEmbvsEP(0), fHistRCPtExLJVSDPhiLJ(0), fHistRCPtExPartialLJVSDPhiLJ(0), fHistEmbJetsPhiEta(0), fHistLeadPartPhiEta(0) { // Standard constructor. fHistRCPt = 0; fHistRCPtExLJ = 0; fHistRCPtExPartialLJ = 0; fHistRCPtRand = 0; fHistRhoVSRCPt = 0; fHistDeltaPtRCvsEP = 0; fHistDeltaPtRCExLJ = 0; fHistDeltaPtRCExPartialLJ = 0; fHistDeltaPtRCRand = 0; fHistEmbJetsPtArea = 0; fHistEmbJetsCorrPtArea = 0; fHistEmbPartPtvsJetPt = 0; fHistEmbPartPtvsJetCorrPt = 0; fHistJetPtvsJetCorrPt = 0; fHistDistLeadPart2JetAxis = 0; fHistEmbBkgArea = 0; fHistRhoVSEmbBkg = 0; fHistDeltaPtEmbArea = 0; fHistDeltaPtEmbvsEP = 0; SetMakeGeneralHistograms(kTRUE); } //________________________________________________________________________ void AliAnalysisTaskDeltaPt::AllocateHistogramArrays() { fHistRCPt = new TH1*[fNcentBins]; fHistRCPtExLJ = new TH1*[fNcentBins]; fHistRCPtExPartialLJ = new TH1*[fNcentBins]; fHistRCPtRand = new TH1*[fNcentBins]; fHistRhoVSRCPt = new TH2*[fNcentBins]; fHistDeltaPtRCvsEP = new TH2*[fNcentBins]; fHistDeltaPtRCExLJ = new TH1*[fNcentBins]; fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins]; fHistDeltaPtRCRand = new TH1*[fNcentBins]; fHistEmbJetsPtArea = new TH3*[fNcentBins]; fHistEmbJetsCorrPtArea = new TH3*[fNcentBins]; fHistEmbPartPtvsJetPt = new TH2*[fNcentBins]; fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins]; fHistJetPtvsJetCorrPt = new TH2*[fNcentBins]; fHistDistLeadPart2JetAxis = new TH1*[fNcentBins]; fHistEmbBkgArea = new TH2*[fNcentBins]; fHistRhoVSEmbBkg = new TH2*[fNcentBins]; fHistDeltaPtEmbArea = new TH2*[fNcentBins]; fHistDeltaPtEmbvsEP = new TH2*[fNcentBins]; for (Int_t i = 0; i < fNcentBins; i++) { fHistRCPt[i] = 0; fHistRCPtExLJ[i] = 0; fHistRCPtExPartialLJ[i] = 0; fHistRCPtRand[i] = 0; fHistRhoVSRCPt[i] = 0; fHistDeltaPtRCvsEP[i] = 0; fHistDeltaPtRCExLJ[i] = 0; fHistDeltaPtRCExPartialLJ[i] = 0; fHistDeltaPtRCRand[i] = 0; fHistEmbJetsPtArea[i] = 0; fHistEmbJetsCorrPtArea[i] = 0; fHistEmbPartPtvsJetPt[i] = 0; fHistEmbPartPtvsJetCorrPt[i] = 0; fHistJetPtvsJetCorrPt[i] = 0; fHistDistLeadPart2JetAxis[i] = 0; fHistEmbBkgArea[i] = 0; fHistRhoVSEmbBkg[i] = 0; fHistDeltaPtEmbArea[i] = 0; fHistDeltaPtEmbvsEP[i] = 0; } } //________________________________________________________________________ void AliAnalysisTaskDeltaPt::UserCreateOutputObjects() { // Create user output. AliAnalysisTaskEmcalJet::UserCreateOutputObjects(); AllocateHistogramArrays(); fJetsCont = GetJetContainer("Jets"); fTracksCont = GetParticleContainer("Tracks"); fCaloClustersCont = GetClusterContainer("CaloClusters"); fEmbJetsCont = GetJetContainer("EmbJets"); fEmbTracksCont = GetParticleContainer("EmbTracks"); fEmbCaloClustersCont = GetClusterContainer("EmbCaloClusters"); fRandTracksCont = GetParticleContainer("RandTracks"); fRandCaloClustersCont = GetClusterContainer("RandCaloClusters"); if (fTracksCont || fCaloClustersCont) { 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 (fJetsCont) { 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); fHistRCPtExPartialLJVSDPhiLJ = new TH2F("fHistRCPtExPartialLJVSDPhiLJ","fHistRCPtExPartialLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8); fHistRCPtExPartialLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})"); fHistRCPtExPartialLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi"); fOutput->Add(fHistRCPtExPartialLJVSDPhiLJ); } } if (fEmbJetsCont) { 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); 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; const Int_t nbinsZ = 12; Double_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000}; Double_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt); Double_t *binsCorrPt = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt); Double_t *binsArea = GenerateFixedBinArray(50, 0, 2); for (Int_t i = 0; i < fNcentBins; i++) { if (fTracksCont || fCaloClustersCont) { 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 = "fHistDeltaPtRCvsEP_"; histname += i; fHistDeltaPtRCvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt); fHistDeltaPtRCvsEP[i]->GetXaxis()->SetTitle("#phi_{RC} - #psi_{RP}"); fHistDeltaPtRCvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})"); fHistDeltaPtRCvsEP[i]->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtRCvsEP[i]); if (fJetsCont) { 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]); histname = "fHistRCPtExPartialLJ_"; histname += i; fHistRCPtExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2); fHistRCPtExPartialLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})"); fHistRCPtExPartialLJ[i]->GetYaxis()->SetTitle("counts"); fOutput->Add(fHistRCPtExPartialLJ[i]); histname = "fHistDeltaPtRCExPartialLJ_"; histname += i; fHistDeltaPtRCExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt); fHistDeltaPtRCExPartialLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})"); fHistDeltaPtRCExPartialLJ[i]->GetYaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtRCExPartialLJ[i]); } } if (fRandTracksCont || fRandCaloClustersCont) { 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 (fEmbJetsCont) { histname = "fHistEmbJetsPtArea_"; histname += i; fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ); fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area"); fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})"); fOutput->Add(fHistEmbJetsPtArea[i]); histname = "fHistEmbJetsCorrPtArea_"; histname += i; fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ); fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area"); fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})"); fOutput->Add(fHistEmbJetsCorrPtArea[i]); histname = "fHistEmbPartPtvsJetPt_"; histname += 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 = "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; 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 = "fHistEmbBkgArea_"; histname += i; fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt); fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area"); 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("#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(), 50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt); fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area"); fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})"); fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtEmbArea[i]); histname = "fHistDeltaPtEmbvsEP_"; histname += i; fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt); fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}"); fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})"); fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtEmbvsEP[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 } //________________________________________________________________________ Bool_t AliAnalysisTaskDeltaPt::FillHistograms() { // Fill histograms. // ************ // Random cones // _________________________________ const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi(); Float_t RCpt = 0; Float_t RCeta = 0; Float_t RCphi = 0; if (fTracksCont || fCaloClustersCont) { for (Int_t i = 0; i < fRCperEvent; i++) { // Simple random cones RCpt = 0; RCeta = 0; RCphi = 0; GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0); if (RCpt > 0) { fHistRCPhiEta->Fill(RCeta, RCphi); fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt); fHistRCPt[fCentBin]->Fill(RCpt); Double_t ep = RCphi - fEPV0; while (ep < 0) ep += TMath::Pi(); while (ep >= TMath::Pi()) ep -= TMath::Pi(); fHistDeltaPtRCvsEP[fCentBin]->Fill(ep, RCpt - rcArea * fRhoVal); } if (fJetsCont) { // Random cones far from leading jet AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho"); RCpt = 0; RCeta = 0; RCphi = 0; GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 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); } //partial exclusion if(fBeamType == kpA) { RCpt = 0; RCeta = 0; RCphi = 0; GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE); 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; fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi); } fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt); fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal); } } } } } // Random cones with randomized particles if (fRandTracksCont || fRandCaloClustersCont) { RCpt = 0; RCeta = 0; RCphi = 0; GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0); if (RCpt > 0) { fHistRCPtRand[fCentBin]->Fill(RCpt); fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal); } } // ************ // Embedding // _________________________________ if (fEmbJetsCont) { AliEmcalJet *embJet = NextEmbeddedJet(kTRUE); while (embJet != 0) { TLorentzVector mom; fEmbJetsCont->GetLeadingHadronMomentum(mom,embJet); Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi())); fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt()); fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal); fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi()); fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet); fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt()); fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), mom.Pt()); 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()); Double_t ep = embJet->Phi() - fEPV0; while (ep < 0) ep += TMath::Pi(); while (ep >= TMath::Pi()) ep -= TMath::Pi(); fHistDeltaPtEmbvsEP[fCentBin]->Fill(ep, embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt()); embJet = NextEmbeddedJet(); } } return kTRUE; } //________________________________________________________________________ AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Bool_t reset) { // Get the next accepted embedded jet. Int_t i = reset ? 0 : -1; AliEmcalJet* jet = fEmbJetsCont->GetNextAcceptJet(i); while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet(); return jet; } //________________________________________________________________________ void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet, Bool_t bPartialExclusion) const { // Get rigid cone. eta = -999; phi = -999; pt = 0; if (!tracks && !clusters) return; Float_t LJeta = 999; Float_t LJphi = 999; if (jet) { LJeta = jet->Eta(); LJphi = jet->Phi(); } Float_t maxEta = fConeMaxEta; Float_t minEta = fConeMinEta; Float_t maxPhi = fConeMaxPhi; Float_t minPhi = fConeMinPhi; if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2; if (minPhi < 0) minPhi = 0; Float_t dLJ = 0; Int_t repeats = 0; Bool_t reject = kTRUE; do { eta = gRandom->Rndm() * (maxEta - minEta) + minEta; phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi; dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi)); if(bPartialExclusion) { reject = kFALSE; TRandom3 rnd; rnd.SetSeed(0); Double_t ncoll = GetNColl(); Double_t prob = 0.; if(ncoll>0) prob = 1./ncoll; if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone } repeats++; } while (dLJ < fMinRC2LJ && repeats < 999 && reject); if (repeats == 999) { AliWarning(Form("%s: Could not get random cone!", GetName())); return; } if (clusters) { AliVCluster* cluster = clusters->GetNextAcceptCluster(0); while (cluster) { TLorentzVector nPart; cluster->GetMomentum(nPart, const_cast(fVertex)); Float_t cluseta = nPart.Eta(); Float_t clusphi = nPart.Phi(); if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi())) clusphi += 2 * TMath::Pi(); if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi())) clusphi -= 2 * TMath::Pi(); Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi)); if (d <= fConeRadius) pt += nPart.Pt(); cluster = clusters->GetNextAcceptCluster(); } } if (tracks) { AliVParticle* track = tracks->GetNextAcceptParticle(0); while(track) { Float_t tracketa = track->Eta(); Float_t trackphi = track->Phi(); if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi())) trackphi += 2 * TMath::Pi(); if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi())) trackphi -= 2 * TMath::Pi(); Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi)); if (d <= fConeRadius) pt += track->Pt(); track = tracks->GetNextAcceptParticle(); } } } //________________________________________________________________________ void AliAnalysisTaskDeltaPt::SetConeEtaPhiEMCAL() { // Set default cuts for full cones SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius); SetConePhiLimits(1.4+fConeRadius,TMath::Pi()-fConeRadius); } //________________________________________________________________________ void AliAnalysisTaskDeltaPt::SetConeEtaPhiTPC() { // Set default cuts for charged cones SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius); SetConePhiLimits(-10, 10); } //________________________________________________________________________ void AliAnalysisTaskDeltaPt::ExecOnce() { // Initialize the analysis. AliAnalysisTaskEmcalJet::ExecOnce(); if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0; if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0; if (fEmbTracksCont && fEmbTracksCont->GetArray() == 0) fEmbTracksCont = 0; if (fEmbCaloClustersCont && fEmbCaloClustersCont->GetArray() == 0) fEmbCaloClustersCont = 0; if (fRandTracksCont && fRandTracksCont->GetArray() == 0) fRandTracksCont = 0; if (fRandCaloClustersCont && fRandCaloClustersCont->GetArray() == 0) fRandCaloClustersCont = 0; if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0; if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0; if (fRCperEvent < 0) { Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi); Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius; fRCperEvent = TMath::FloorNint(area / rcArea - 0.5); if (fRCperEvent == 0) fRCperEvent = 1; } if (fMinRC2LJ < 0) fMinRC2LJ = fConeRadius * 1.5; const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2; if (fMinRC2LJ > maxDist) { AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. " "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist)); fMinRC2LJ = maxDist; } } //________________________________________________________________________ Double_t AliAnalysisTaskDeltaPt::GetNColl() const { // Get NColl - returns value of corresponding bin // only works for pA // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from if(fBeamType == kpA) { const Int_t nNCollBins = 7; Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.}; Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.}; Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52}; for(Int_t i = 0; i=centMin[i] && fCent