// $Id$ // // Emcal jet response matrix maker task. // // Author: S. Aiola #include "AliJetResponseMaker.h" #include #include #include #include #include "AliAnalysisManager.h" #include "AliVCluster.h" #include "AliVTrack.h" #include "AliEmcalJet.h" #include "AliLog.h" #include "AliRhoParameter.h" #include "AliNamedArrayI.h" #include "AliJetContainer.h" #include "AliParticleContainer.h" #include "AliClusterContainer.h" ClassImp(AliJetResponseMaker) //________________________________________________________________________ AliJetResponseMaker::AliJetResponseMaker() : AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE), fMatching(kNoMatching), fMatchingPar1(0), fMatchingPar2(0), fUseCellsToMatch(kFALSE), fMinJetMCPt(1), fHistoType(0), fDeltaPtAxis(0), fDeltaEtaDeltaPhiAxis(0), fNEFAxis(0), fZAxis(0), fIsJet1Rho(kFALSE), fIsJet2Rho(kFALSE), fHistRejectionReason1(0), fHistRejectionReason2(0), fHistJets1(0), fHistJets2(0), fHistMatching(0), fHistJets1PhiEta(0), fHistJets1PtArea(0), fHistJets1CorrPtArea(0), fHistJets1NEFvsPt(0), fHistJets1CEFvsCEFPt(0), fHistJets1ZvsPt(0), fHistJets2PhiEta(0), fHistJets2PtArea(0), fHistJets2CorrPtArea(0), fHistJets2NEFvsPt(0), fHistJets2CEFvsCEFPt(0), fHistJets2ZvsPt(0), fHistCommonEnergy1vsJet1Pt(0), fHistCommonEnergy2vsJet2Pt(0), fHistDistancevsJet1Pt(0), fHistDistancevsJet2Pt(0), fHistDistancevsCommonEnergy1(0), fHistDistancevsCommonEnergy2(0), fHistCommonEnergy1vsCommonEnergy2(0), fHistDeltaEtaDeltaPhi(0), fHistJet2PtOverJet1PtvsJet2Pt(0), fHistJet1PtOverJet2PtvsJet1Pt(0), fHistDeltaPtvsJet1Pt(0), fHistDeltaPtvsJet2Pt(0), fHistDeltaPtOverJet1PtvsJet1Pt(0), fHistDeltaPtOverJet2PtvsJet2Pt(0), fHistDeltaPtvsDistance(0), fHistDeltaPtvsCommonEnergy1(0), fHistDeltaPtvsCommonEnergy2(0), fHistDeltaPtvsArea1(0), fHistDeltaPtvsArea2(0), fHistDeltaPtvsDeltaArea(0), fHistJet1PtvsJet2Pt(0), fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0), fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0), fHistDeltaCorrPtvsJet1CorrPt(0), fHistDeltaCorrPtvsJet2CorrPt(0), fHistDeltaCorrPtvsDistance(0), fHistDeltaCorrPtvsCommonEnergy1(0), fHistDeltaCorrPtvsCommonEnergy2(0), fHistDeltaCorrPtvsArea1(0), fHistDeltaCorrPtvsArea2(0), fHistDeltaCorrPtvsDeltaArea(0), fHistJet1CorrPtvsJet2CorrPt(0), fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0), fHistDeltaMCPtOverJet2PtvsJet2Pt(0), fHistDeltaMCPtvsJet1MCPt(0), fHistDeltaMCPtvsJet2Pt(0), fHistDeltaMCPtvsDistance(0), fHistDeltaMCPtvsCommonEnergy1(0), fHistDeltaMCPtvsCommonEnergy2(0), fHistDeltaMCPtvsArea1(0), fHistDeltaMCPtvsArea2(0), fHistDeltaMCPtvsDeltaArea(0), fHistJet1MCPtvsJet2Pt(0) { // Default constructor. SetMakeGeneralHistograms(kTRUE); } //________________________________________________________________________ AliJetResponseMaker::AliJetResponseMaker(const char *name) : AliAnalysisTaskEmcalJet(name, kTRUE), fMatching(kNoMatching), fMatchingPar1(0), fMatchingPar2(0), fUseCellsToMatch(kFALSE), fMinJetMCPt(1), fHistoType(0), fDeltaPtAxis(0), fDeltaEtaDeltaPhiAxis(0), fNEFAxis(0), fZAxis(0), fIsJet1Rho(kFALSE), fIsJet2Rho(kFALSE), fHistRejectionReason1(0), fHistRejectionReason2(0), fHistJets1(0), fHistJets2(0), fHistMatching(0), fHistJets1PhiEta(0), fHistJets1PtArea(0), fHistJets1CorrPtArea(0), fHistJets1NEFvsPt(0), fHistJets1CEFvsCEFPt(0), fHistJets1ZvsPt(0), fHistJets2PhiEta(0), fHistJets2PtArea(0), fHistJets2CorrPtArea(0), fHistJets2NEFvsPt(0), fHistJets2CEFvsCEFPt(0), fHistJets2ZvsPt(0), fHistCommonEnergy1vsJet1Pt(0), fHistCommonEnergy2vsJet2Pt(0), fHistDistancevsJet1Pt(0), fHistDistancevsJet2Pt(0), fHistDistancevsCommonEnergy1(0), fHistDistancevsCommonEnergy2(0), fHistCommonEnergy1vsCommonEnergy2(0), fHistDeltaEtaDeltaPhi(0), fHistJet2PtOverJet1PtvsJet2Pt(0), fHistJet1PtOverJet2PtvsJet1Pt(0), fHistDeltaPtvsJet1Pt(0), fHistDeltaPtvsJet2Pt(0), fHistDeltaPtOverJet1PtvsJet1Pt(0), fHistDeltaPtOverJet2PtvsJet2Pt(0), fHistDeltaPtvsDistance(0), fHistDeltaPtvsCommonEnergy1(0), fHistDeltaPtvsCommonEnergy2(0), fHistDeltaPtvsArea1(0), fHistDeltaPtvsArea2(0), fHistDeltaPtvsDeltaArea(0), fHistJet1PtvsJet2Pt(0), fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0), fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0), fHistDeltaCorrPtvsJet1CorrPt(0), fHistDeltaCorrPtvsJet2CorrPt(0), fHistDeltaCorrPtvsDistance(0), fHistDeltaCorrPtvsCommonEnergy1(0), fHistDeltaCorrPtvsCommonEnergy2(0), fHistDeltaCorrPtvsArea1(0), fHistDeltaCorrPtvsArea2(0), fHistDeltaCorrPtvsDeltaArea(0), fHistJet1CorrPtvsJet2CorrPt(0), fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0), fHistDeltaMCPtOverJet2PtvsJet2Pt(0), fHistDeltaMCPtvsJet1MCPt(0), fHistDeltaMCPtvsJet2Pt(0), fHistDeltaMCPtvsDistance(0), fHistDeltaMCPtvsCommonEnergy1(0), fHistDeltaMCPtvsCommonEnergy2(0), fHistDeltaMCPtvsArea1(0), fHistDeltaMCPtvsArea2(0), fHistDeltaMCPtvsDeltaArea(0), fHistJet1MCPtvsJet2Pt(0) { // Standard constructor. SetMakeGeneralHistograms(kTRUE); } //________________________________________________________________________ AliJetResponseMaker::~AliJetResponseMaker() { // Destructor } //________________________________________________________________________ void AliJetResponseMaker::AllocateTH2() { // Allocate TH2 histograms. fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2); fHistJets1PhiEta->GetXaxis()->SetTitle("#eta"); fHistJets1PhiEta->GetYaxis()->SetTitle("#phi"); fOutput->Add(fHistJets1PhiEta); fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt); fHistJets1PtArea->GetXaxis()->SetTitle("area"); fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)"); fHistJets1PtArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets1PtArea); if (fIsJet1Rho) { fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistJets1CorrPtArea->GetXaxis()->SetTitle("area"); fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)"); fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets1CorrPtArea); } fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets1ZvsPt->GetXaxis()->SetTitle("Z"); fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)"); fHistJets1ZvsPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets1ZvsPt); fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF"); fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)"); fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets1NEFvsPt); fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF"); fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)"); fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets1CEFvsCEFPt); // Jets 2 histograms fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2); fHistJets2PhiEta->GetXaxis()->SetTitle("#eta"); fHistJets2PhiEta->GetYaxis()->SetTitle("#phi"); fOutput->Add(fHistJets2PhiEta); fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt); fHistJets2PtArea->GetXaxis()->SetTitle("area"); fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); fHistJets2PtArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets2PtArea); if (fIsJet2Rho) { fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistJets2CorrPtArea->GetXaxis()->SetTitle("area"); fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)"); fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets2CorrPtArea); } fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets2ZvsPt->GetXaxis()->SetTitle("Z"); fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); fHistJets2ZvsPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets2ZvsPt); fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF"); fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets2NEFvsPt); fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF"); fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)"); fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJets2CEFvsCEFPt); // Matching histograms fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)"); fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}"); fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistCommonEnergy1vsJet1Pt); fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)"); fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistCommonEnergy2vsJet2Pt); fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance"); fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}"); fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDistancevsJet1Pt); fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt); fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance"); fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDistancevsJet2Pt); fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2); fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance"); fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)"); fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDistancevsCommonEnergy1); fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2); fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance"); fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)"); fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDistancevsCommonEnergy2); fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2); fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)"); fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)"); fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistCommonEnergy1vsCommonEnergy2); fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2); fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)"); fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)"); fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaEtaDeltaPhi); fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5); fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}"); fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt); fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5); fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}"); fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}"); fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt); fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}"); fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsJet1Pt); fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsJet2Pt); fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}"); fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}"); fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt); fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}"); fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt); fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance"); fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsDistance); fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)"); fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsCommonEnergy1); fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)"); fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsCommonEnergy2); fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}"); fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsArea1); fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}"); fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsArea2); fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea", fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}"); fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaPtvsDeltaArea); fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}"); fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJet1PtvsJet2Pt); if (fIsJet1Rho || fIsJet2Rho) { if (!fIsJet1Rho) fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); else fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}"); fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt); if (!fIsJet2Rho) fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); else fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}"); fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt); fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5); fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}"); fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}"); fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt); fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5); fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}"); fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}"); fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt); fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance"); fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsDistance); fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)"); fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1); fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)"); fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2); fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}"); fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsArea1); fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}"); fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsArea2); fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea", fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}"); fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)"); fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaCorrPtvsDeltaArea); if (!fIsJet1Rho) fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); else if (!fIsJet2Rho) fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); else fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}"); fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}"); fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJet1CorrPtvsJet2CorrPt); } if (fIsEmbedded) { fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}"); fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsJet1MCPt); fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsJet2Pt); fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}"); fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}"); fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt); fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5); fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}"); fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}"); fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt); fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance"); fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsDistance); fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)"); fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsCommonEnergy1); fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)"); fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsCommonEnergy2); fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}"); fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsArea1); fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2", fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}"); fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsArea2); fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea", fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt); fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}"); fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)"); fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistDeltaMCPtvsDeltaArea); fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}"); fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}"); fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts"); fOutput->Add(fHistJet1MCPtvsJet2Pt); } } //________________________________________________________________________ void AliJetResponseMaker::AllocateTHnSparse() { // Allocate THnSparse histograms. TString title[25]= {""}; Int_t nbins[25] = {0}; Double_t min[25] = {0.}; Double_t max[25] = {0.}; Int_t dim = 0; title[dim] = "#phi"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1)); dim++; title[dim] = "#eta"; nbins[dim] = fNbins/4; min[dim] = -1; max[dim] = 1; dim++; title[dim] = "p_{T}"; nbins[dim] = fNbins; min[dim] = 0; max[dim] = 250; dim++; title[dim] = "A_{jet}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.5; dim++; if (fNEFAxis) { title[dim] = "NEF"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.2; dim++; } 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] = 120; dim++; Int_t dim1 = dim, dim2 = dim; if (fIsJet1Rho) { title[dim1] = "p_{T}^{corr}"; nbins[dim1] = fNbins*2; min[dim1] = -250; max[dim1] = 250; dim1++; } if (fIsEmbedded) { title[dim1] = "p_{T}^{MC}"; nbins[dim1] = fNbins; min[dim1] = 0; max[dim1] = 250; dim1++; } fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max); for (Int_t i = 0; i < dim1; i++) fHistJets1->GetAxis(i)->SetTitle(title[i]); fOutput->Add(fHistJets1); if (fIsJet2Rho) { title[dim2] = "p_{T}^{corr}"; nbins[dim2] = fNbins*2; min[dim2] = -250; max[dim2] = 250; dim2++; } fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max); for (Int_t i = 0; i < dim2; i++) fHistJets2->GetAxis(i)->SetTitle(title[i]); fOutput->Add(fHistJets2); // Matching dim = 0; title[dim] = "p_{T,1}"; nbins[dim] = fNbins; min[dim] = 0; max[dim] = 250; dim++; title[dim] = "p_{T,2}"; nbins[dim] = fNbins; min[dim] = 0; max[dim] = 250; dim++; title[dim] = "A_{jet,1}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.5; dim++; title[dim] = "A_{jet,2}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.5; dim++; title[dim] = "distance"; nbins[dim] = fNbins/2; min[dim] = 0; max[dim] = 1.2; dim++; title[dim] = "CE1"; nbins[dim] = fNbins/2; min[dim] = 0; max[dim] = 1.2; dim++; title[dim] = "CE2"; nbins[dim] = fNbins/2; min[dim] = 0; max[dim] = 1.2; dim++; title[dim] = "p_{T,particle,1}^{leading} (GeV/c)"; nbins[dim] = 120; min[dim] = 0; max[dim] = 120; dim++; title[dim] = "p_{T,particle,2}^{leading} (GeV/c)"; nbins[dim] = 120; min[dim] = 0; max[dim] = 120; dim++; if (fDeltaPtAxis) { title[dim] = "#deltaA_{jet}"; nbins[dim] = fNbins/2; min[dim] = -1; max[dim] = 1; dim++; title[dim] = "#deltap_{T}"; nbins[dim] = fNbins*2; min[dim] = -250; max[dim] = 250; dim++; } if (fIsJet1Rho) { title[dim] = "p_{T,1}^{corr}"; nbins[dim] = fNbins*2; min[dim] = -250; max[dim] = 250; dim++; } if (fIsJet2Rho) { title[dim] = "p_{T,2}^{corr}"; nbins[dim] = fNbins*2; min[dim] = -250; max[dim] = 250; dim++; } if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) { title[dim] = "#deltap_{T}^{corr}"; nbins[dim] = fNbins*2; min[dim] = -250; max[dim] = 250; dim++; } if (fDeltaEtaDeltaPhiAxis) { title[dim] = "#delta#eta"; nbins[dim] = fNbins/2; min[dim] = -1; max[dim] = 1; dim++; title[dim] = "#delta#phi"; nbins[dim] = fNbins/2; min[dim] = -TMath::Pi()/2; max[dim] = TMath::Pi()*3/2; dim++; } if (fIsEmbedded) { title[dim] = "p_{T,1}^{MC}"; nbins[dim] = fNbins; min[dim] = 0; max[dim] = 250; dim++; if (fDeltaPtAxis) { title[dim] = "#deltap_{T}^{MC}"; nbins[dim] = fNbins*2; min[dim] = -250; max[dim] = 250; dim++; } } if (fNEFAxis) { title[dim] = "NEF_{1}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.2; dim++; title[dim] = "NEF_{2}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.2; dim++; } if (fZAxis) { title[dim] = "Z_{1}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.2; dim++; title[dim] = "Z_{2}"; nbins[dim] = fNbins/4; min[dim] = 0; max[dim] = 1.2; dim++; } fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max); for (Int_t i = 0; i < dim; i++) fHistMatching->GetAxis(i)->SetTitle(title[i]); fOutput->Add(fHistMatching); } //________________________________________________________________________ void AliJetResponseMaker::UserCreateOutputObjects() { // Create user objects. AliAnalysisTaskEmcalJet::UserCreateOutputObjects(); AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets2) return; if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE; else fIsJet1Rho = kTRUE; if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE; else fIsJet2Rho = kTRUE; fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250); fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason"); fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)"); fHistRejectionReason1->GetZaxis()->SetTitle("counts"); SetRejectionReasonLabels(fHistRejectionReason1->GetXaxis()); fOutput->Add(fHistRejectionReason1); fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250); fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason"); fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)"); fHistRejectionReason2->GetZaxis()->SetTitle("counts"); SetRejectionReasonLabels(fHistRejectionReason2->GetXaxis()); fOutput->Add(fHistRejectionReason2); if (fHistoType==0) AllocateTH2(); else AllocateTHnSparse(); PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram } //________________________________________________________________________ void AliJetResponseMaker::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; if (Set==1) histo = fHistJets1; else if (Set==2) histo = fHistJets2; if (!histo) return; Double_t contents[20]={0}; for (Int_t i = 0; i < histo->GetNdimensions(); i++) { TString title(histo->GetAxis(i)->GetTitle()); if (title=="#phi") contents[i] = Phi; else if (title=="#eta") contents[i] = Eta; else if (title=="p_{T}") contents[i] = Pt; else if (title=="A_{jet}") contents[i] = A; else if (title=="NEF") contents[i] = NEF; else if (title=="Z") contents[i] = LeadingPt/Pt; else if (title=="p_{T}^{corr}") contents[i] = CorrPt; else if (title=="p_{T}^{MC}") contents[i] = MCPt; else if (title=="p_{T,particle}^{leading} (GeV/c)") contents[i] = LeadingPt; else AliWarning(Form("Unable to fill dimension %s!",title.Data())); } histo->Fill(contents); } else { if (Set == 1) { fHistJets1PtArea->Fill(A, Pt); fHistJets1PhiEta->Fill(Eta, Phi); fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt); fHistJets1NEFvsPt->Fill(NEF, Pt); fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt); if (fIsJet1Rho) fHistJets1CorrPtArea->Fill(A, CorrPt); } else if (Set == 2) { fHistJets2PtArea->Fill(A, Pt); fHistJets2PhiEta->Fill(Eta, Phi); fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt); fHistJets2NEFvsPt->Fill(NEF, Pt); fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt); if (fIsJet2Rho) fHistJets2CorrPtArea->Fill(A, CorrPt); } } } //________________________________________________________________________ void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2) { if (fHistoType==1) { Double_t contents[20]={0}; for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) { TString title(fHistMatching->GetAxis(i)->GetTitle()); if (title=="p_{T,1}") contents[i] = Pt1; else if (title=="p_{T,2}") contents[i] = Pt2; else if (title=="A_{jet,1}") contents[i] = A1; else if (title=="A_{jet,2}") contents[i] = A2; else if (title=="distance") contents[i] = d; else if (title=="CE1") contents[i] = CE1; else if (title=="CE2") contents[i] = CE2; else if (title=="#deltaA_{jet}") contents[i] = A1-A2; else if (title=="#deltap_{T}") contents[i] = Pt1-Pt2; else if (title=="#delta#eta") contents[i] = Eta1-Eta2; else if (title=="#delta#phi") contents[i] = Phi1-Phi2; else if (title=="p_{T,1}^{corr}") contents[i] = CorrPt1; else if (title=="p_{T,2}^{corr}") contents[i] = CorrPt2; else if (title=="#deltap_{T}^{corr}") contents[i] = CorrPt1-CorrPt2; else if (title=="p_{T,1}^{MC}") contents[i] = MCPt1; else if (title=="#deltap_{T}^{MC}") contents[i] = MCPt1-Pt2; else if (title=="NEF_{1}") contents[i] = NEF1; else if (title=="NEF_{2}") contents[i] = NEF2; else if (title=="Z_{1}") contents[i] = LeadingPt1/Pt1; else if (title=="Z_{2}") contents[i] = LeadingPt2/Pt2; else if (title=="p_{T,particle,1}^{leading} (GeV/c)") contents[i] = LeadingPt1; else if (title=="p_{T,particle,2}^{leading} (GeV/c)") contents[i] = LeadingPt2; else AliWarning(Form("Unable to fill dimension %s!",title.Data())); } fHistMatching->Fill(contents); } else { fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1); fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2); fHistDistancevsJet1Pt->Fill(d, Pt1); fHistDistancevsJet2Pt->Fill(d, Pt2); fHistDistancevsCommonEnergy1->Fill(d, CE1); fHistDistancevsCommonEnergy2->Fill(d, CE2); fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2); fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2); fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1); fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2); Double_t dpt = Pt1 - Pt2; fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt); fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt); fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1); fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2); fHistDeltaPtvsDistance->Fill(d, dpt); fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt); fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt); fHistDeltaPtvsArea1->Fill(A1, dpt); fHistDeltaPtvsArea2->Fill(A2, dpt); Double_t darea = A1 - A2; fHistDeltaPtvsDeltaArea->Fill(darea, dpt); fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2); if (fIsJet1Rho || fIsJet2Rho) { Double_t dcorrpt = CorrPt1 - CorrPt2; fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt); fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt); fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1); fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2); fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt); fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt); fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt); fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt); fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt); fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt); fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2); } if (fIsEmbedded) { Double_t dmcpt = MCPt1 - Pt2; fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt); fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt); fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1); fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2); fHistDeltaMCPtvsDistance->Fill(d, dmcpt); fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt); fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt); fHistDeltaMCPtvsArea1->Fill(A1, dmcpt); fHistDeltaMCPtvsArea2->Fill(A2, dmcpt); fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt); fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2); } } } //________________________________________________________________________ void AliJetResponseMaker::ExecOnce() { // Execute once. AliAnalysisTaskEmcalJet::ExecOnce(); AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) { if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) { AliWarning("Changing matching type from MC label to same collection..."); fMatching = kSameCollections; } else { AliWarning("Changing matching type from MC label to geometrical..."); fMatching = kGeometrical; } } else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) { if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) { AliWarning("Changing matching type from same collection to MC label..."); fMatching = kMCLabel; } else { AliWarning("Changing matching type from same collection to geometrical..."); fMatching = kGeometrical; } } } //________________________________________________________________________ Bool_t AliJetResponseMaker::Run() { // Find the closest jets if (fMatching == kNoMatching) return kTRUE; else return DoJetMatching(); } //________________________________________________________________________ Bool_t AliJetResponseMaker::DoJetMatching() { AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE; DoJetLoop(); AliEmcalJet* jet1 = 0; jets1->ResetCurrentID(); while ((jet1 = jets1->GetNextJet())) { AliEmcalJet *jet2 = jet1->ClosestJet(); if (!jet2) continue; if (jet2->ClosestJet() != jet1) continue; if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue; // 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() { // Do the jet loop. AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; AliEmcalJet* jet1 = 0; AliEmcalJet* jet2 = 0; jets2->ResetCurrentID(); while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching(); jets1->ResetCurrentID(); while ((jet1 = jets1->GetNextJet())) { jet1->ResetMatching(); if (jet1->MCPt() < fMinJetMCPt) continue; jets2->ResetCurrentID(); while ((jet2 = jets2->GetNextJet())) { SetMatchingLevel(jet1, jet2, fMatching); } // jet2 loop } // jet1 loop } //________________________________________________________________________ void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const { Double_t deta = jet2->Eta() - jet1->Eta(); Double_t dphi = jet2->Phi() - jet1->Phi(); d = TMath::Sqrt(deta * deta + dphi * dphi); } //________________________________________________________________________ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const { AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; AliParticleContainer *tracks1 = jets1->GetParticleContainer(); AliClusterContainer *clusters1 = jets1->GetClusterContainer(); AliParticleContainer *tracks2 = jets2->GetParticleContainer(); // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated d1 = jet1->Pt(); d2 = jet2->Pt(); Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background // remove completely tracks that are not MC particles (label == 0) if (tracks1 && tracks1->GetArray()) { for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) { AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray()); if (!track) { AliWarning(Form("Could not find track %d!", iTrack)); continue; } Int_t MClabel = TMath::Abs(track->GetLabel()); MClabel -= fMCLabelShift; if (MClabel != 0) continue; // this is not a MC particle; remove it completely AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel)); totalPt1 -= track->Pt(); d1 -= track->Pt(); } } // remove completely clusters that are not MC particles (label == 0) if (fUseCellsToMatch && fCaloCells) { for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray()); if (!clus) { AliWarning(Form("Could not find cluster %d!", iClus)); continue; } TLorentzVector part; clus->GetMomentum(part, const_cast(fVertex)); for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) { Int_t cellId = clus->GetCellAbsId(iCell); Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell); Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId)); MClabel -= fMCLabelShift; if (MClabel != 0) continue; // this is not a MC particle; remove it completely AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel)); totalPt1 -= part.Pt() * cellFrac; d1 -= part.Pt() * cellFrac; } } } else { for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray()); if (!clus) { AliWarning(Form("Could not find cluster %d!", iClus)); continue; } TLorentzVector part; clus->GetMomentum(part, const_cast(fVertex)); Int_t MClabel = TMath::Abs(clus->GetLabel()); MClabel -= fMCLabelShift; if (MClabel != 0) continue; // this is not a MC particle; remove it completely AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel)); totalPt1 -= part.Pt(); d1 -= part.Pt(); } } for (Int_t 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(fVertex)); for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) { Int_t cellId = clus->GetCellAbsId(iCell); Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell); Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId)); MClabel -= fMCLabelShift; if (MClabel <= 0) continue; Int_t index1 = -1; index1 = tracks2->GetIndexFromLabel(MClabel); if (index1 < 0) { AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel)); continue; } if (index2 != index1) continue; // found common particle d1 -= part.Pt() * cellFrac; if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found) AliVParticle *MCpart = tracks2->GetParticle(index2); AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!", iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi())); d2 -= MCpart->Pt() * cellFrac; } track2Found = kTRUE; } } } else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) { AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters); if (!clus) { AliWarning(Form("Could not find cluster %d!", iClus)); continue; } TLorentzVector part; clus->GetMomentum(part, const_cast(fVertex)); 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) 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; } } } if (d1 < 0) d1 = 0; if (d2 < 0) d2 = 0; if (totalPt1 < 1) d1 = -1; else d1 /= totalPt1; if (jet2->Pt() < 1) d2 = -1; else d2 /= jet2->Pt(); } //________________________________________________________________________ void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const { AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return; AliParticleContainer *tracks1 = jets1->GetParticleContainer(); AliClusterContainer *clusters1 = jets1->GetClusterContainer(); AliParticleContainer *tracks2 = jets2->GetParticleContainer(); AliClusterContainer *clusters2 = jets2->GetClusterContainer(); // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated d1 = jet1->Pt(); d2 = jet2->Pt(); if (tracks1 && tracks2) { for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) { Int_t index2 = jet2->TrackAt(iTrack2); for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) { Int_t index1 = jet1->TrackAt(iTrack1); if (index2 == index1) { // found common particle AliVParticle *part1 = tracks1->GetParticle(index1); if (!part1) { AliWarning(Form("Could not find track %d!", index1)); continue; } AliVParticle *part2 = tracks2->GetParticle(index2); if (!part2) { AliWarning(Form("Could not find track %d!", index2)); continue; } d1 -= part1->Pt(); d2 -= part2->Pt(); break; } } } } if (clusters1 && clusters2) { if (fUseCellsToMatch) { const Int_t nClus1 = jet1->GetNumberOfClusters(); Int_t ncells1[nClus1]; UShort_t *cellsId1[nClus1]; Double_t *cellsFrac1[nClus1]; Int_t *sortedIndexes1[nClus1]; Double_t ptClus1[nClus1]; for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) { Int_t index1 = jet1->ClusterAt(iClus1); AliVCluster *clus1 = clusters1->GetCluster(index1); if (!clus1) { AliWarning(Form("Could not find cluster %d!", index1)); ncells1[iClus1] = 0; cellsId1[iClus1] = 0; cellsFrac1[iClus1] = 0; sortedIndexes1[iClus1] = 0; ptClus1[iClus1] = 0; continue; } TLorentzVector part1; clus1->GetMomentum(part1, const_cast(fVertex)); ncells1[iClus1] = clus1->GetNCells(); cellsId1[iClus1] = clus1->GetCellsAbsId(); cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction(); sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]]; ptClus1[iClus1] = part1.Pt(); TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE); } const Int_t nClus2 = jet2->GetNumberOfClusters(); const Int_t maxNcells2 = 11520; Int_t sortedIndexes2[maxNcells2]; for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) { Int_t index2 = jet2->ClusterAt(iClus2); AliVCluster *clus2 = clusters2->GetCluster(index2); if (!clus2) { AliWarning(Form("Could not find cluster %d!", index2)); continue; } Int_t ncells2 = clus2->GetNCells(); if (ncells2 >= maxNcells2) { AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2)); continue; } UShort_t *cellsId2 = clus2->GetCellsAbsId(); Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction(); TLorentzVector part2; clus2->GetMomentum(part2, const_cast(fVertex)); Double_t ptClus2 = part2.Pt(); TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE); for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) { if (sortedIndexes1[iClus1] == 0) continue; Int_t iCell1 = 0, iCell2 = 0; while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) { if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1]; d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2; iCell1++; iCell2++; } else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) { iCell2++; } else { iCell1++; } } } } } else { for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) { Int_t index2 = jet2->ClusterAt(iClus2); for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) { Int_t index1 = jet1->ClusterAt(iClus1); if (index2 == index1) { // found common particle AliVCluster *clus1 = clusters1->GetCluster(index1); if (!clus1) { AliWarning(Form("Could not find cluster %d!", index1)); continue; } AliVCluster *clus2 = clusters2->GetCluster(index2); if (!clus2) { AliWarning(Form("Could not find cluster %d!", index2)); continue; } TLorentzVector part1, part2; clus1->GetMomentum(part1, const_cast(fVertex)); clus2->GetMomentum(part2, const_cast(fVertex)); d1 -= part1.Pt(); d2 -= part2.Pt(); break; } } } } } if (d1 < 0) d1 = 0; if (d2 < 0) d2 = 0; if (jet1->Pt() > 0) d1 /= jet1->Pt(); else d1 = -1; if (jet2->Pt() > 0) d2 /= jet2->Pt(); else d2 = -1; } //________________________________________________________________________ void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) { Double_t d1 = -1; Double_t d2 = -1; switch (matching) { case kGeometrical: GetGeometricalMatchingLevel(jet1,jet2,d1); d2 = d1; break; case kMCLabel: // jet1 = detector level and jet2 = particle level! GetMCLabelMatchingLevel(jet1,jet2,d1,d2); break; case kSameCollections: GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2); break; default: ; } if (d1 >= 0) { if (d1 < jet1->ClosestJetDistance()) { jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance()); jet1->SetClosestJet(jet2, d1); } else if (d1 < jet1->SecondClosestJetDistance()) { jet1->SetSecondClosestJet(jet2, d1); } } if (d2 >= 0) { if (d2 < jet2->ClosestJetDistance()) { jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance()); jet2->SetClosestJet(jet1, d2); } else if (d2 < jet2->SecondClosestJetDistance()) { jet2->SetSecondClosestJet(jet1, d2); } } } //________________________________________________________________________ Bool_t AliJetResponseMaker::FillHistograms() { // Fill histograms. AliJetContainer *jets1 = static_cast(fJetCollArray.At(0)); AliJetContainer *jets2 = static_cast(fJetCollArray.At(1)); if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE; AliEmcalJet* jet1 = 0; AliEmcalJet* jet2 = 0; jets2->ResetCurrentID(); while ((jet2 = jets2->GetNextJet())) { AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID())); if (jet2->Pt() < jets2->GetJetPtCut()) continue; Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2); Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area(); if (jets2->AcceptJet(jet2)) FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2, corrpt2, jet2->MCPt(), 2); else fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(), jet2->Pt()); jet1 = jet2->MatchedJet(); if (!jet1) continue; if (!jets1->AcceptJet(jet1)) continue; if (jet1->MCPt() < fMinJetMCPt) continue; Double_t d=-1, ce1=-1, ce2=-1; if (jet2->GetMatchingType() == kGeometrical) { if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2); else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2); d = jet2->ClosestJetDistance(); } else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) { GetGeometricalMatchingLevel(jet1, jet2, d); ce1 = jet1->ClosestJetDistance(); ce2 = jet2->ClosestJetDistance(); } Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1); Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area(); FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(), jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(), jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2); } jets1->ResetCurrentID(); while ((jet1 = jets1->GetNextJet())) { if (!jets1->AcceptJet(jet1)) { fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(), jet1->Pt()); continue; } if (jet1->MCPt() < fMinJetMCPt) continue; AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID())); Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1); Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area(); FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1, corrpt1, jet1->MCPt(), 1); } return kTRUE; }