#include "AliJetResponseMaker.h"
#include <TClonesArray.h>
-#include <TH1F.h>
#include <TH2F.h>
-#include <TH3F.h>
+#include <THnSparse.h>
#include <TLorentzVector.h>
-#include <TSystem.h>
-#include <TFile.h>
-#include <TKey.h>
-#include <TProfile.h>
#include "AliAnalysisManager.h"
#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
-#include "AliGenPythiaEventHeader.h"
-#include "AliMCEvent.h"
#include "AliLog.h"
#include "AliRhoParameter.h"
+#include "AliNamedArrayI.h"
+#include "AliJetContainer.h"
+#include "AliParticleContainer.h"
+#include "AliClusterContainer.h"
ClassImp(AliJetResponseMaker)
//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker() :
AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
- fTracks2Name(""),
- fCalo2Name(""),
- fJets2Name(""),
- fRho2Name(""),
- fPtBiasJet2Track(0),
- fPtBiasJet2Clus(0),
- fAreCollections1MC(kFALSE),
- fAreCollections2MC(kTRUE),
fMatching(kNoMatching),
fMatchingPar1(0),
fMatchingPar2(0),
- fJet2MinEta(-999),
- fJet2MaxEta(-999),
- fJet2MinPhi(-999),
- fJet2MaxPhi(-999),
- fSelectPtHardBin(-999),
- fPythiaHeader(0),
- fPtHardBin(0),
- fNTrials(0),
- fTracks2(0),
- fCaloClusters2(0),
- fJets2(0),
- fRho2(0),
- fRho2Val(0),
- fTracks2Map(0),
- fHistTrialsAfterSel(0),
- fHistEventsAfterSel(0),
- fHistTrials(0),
- fHistXsection(0),
- fHistEvents(0),
+ fUseCellsToMatch(kFALSE),
+ fMinJetMCPt(1),
+ fHistoType(0),
+ fDeltaPtAxis(0),
+ fDeltaEtaDeltaPhiAxis(0),
+ fNEFAxis(0),
+ fZAxis(0),
+ fIsJet1Rho(kFALSE),
+ fIsJet2Rho(kFALSE),
+ fHistJets1(0),
+ fHistJets2(0),
+ fHistMatching(0),
fHistJets1PhiEta(0),
fHistJets1PtArea(0),
fHistJets1CorrPtArea(0),
- fHistLeadingJets1PtArea(0),
- fHistLeadingJets1CorrPtArea(0),
+ fHistJets1NEFvsPt(0),
+ fHistJets1CEFvsCEFPt(0),
+ fHistJets1ZvsPt(0),
fHistJets2PhiEta(0),
fHistJets2PtArea(0),
fHistJets2CorrPtArea(0),
- fHistLeadingJets2PtArea(0),
- fHistLeadingJets2CorrPtArea(0),
- fHistJets2PhiEtaAcceptance(0),
- fHistJets2PtAreaAcceptance(0),
- fHistJets2CorrPtAreaAcceptance(0),
- fHistLeadingJets2PtAreaAcceptance(0),
- fHistLeadingJets2CorrPtAreaAcceptance(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),
- fHistDeltaEtaPhi(0),
fHistDeltaPtvsJet1Pt(0),
fHistDeltaPtvsJet2Pt(0),
- fHistDeltaPtvsMatchingLevel(0),
- fHistDeltaCorrPtvsJet1Pt(0),
- fHistDeltaCorrPtvsJet2Pt(0),
- fHistDeltaCorrPtvsMatchingLevel(0),
- fHistNonMatchedJets1PtArea(0),
- fHistNonMatchedJets2PtArea(0),
- fHistNonMatchedJets1CorrPtArea(0),
- fHistNonMatchedJets2CorrPtArea(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),
- fHistMissedJets2PtArea(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.
//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker(const char *name) :
AliAnalysisTaskEmcalJet(name, kTRUE),
- fTracks2Name("MCParticles"),
- fCalo2Name(""),
- fJets2Name("MCJets"),
- fRho2Name(""),
- fPtBiasJet2Track(0),
- fPtBiasJet2Clus(0),
- fAreCollections1MC(kFALSE),
- fAreCollections2MC(kTRUE),
fMatching(kNoMatching),
fMatchingPar1(0),
fMatchingPar2(0),
- fJet2MinEta(-999),
- fJet2MaxEta(-999),
- fJet2MinPhi(-999),
- fJet2MaxPhi(-999),
- fSelectPtHardBin(-999),
- fPythiaHeader(0),
- fPtHardBin(0),
- fNTrials(0),
- fTracks2(0),
- fCaloClusters2(0),
- fJets2(0),
- fRho2(0),
- fRho2Val(0),
- fTracks2Map(0),
- fHistTrialsAfterSel(0),
- fHistEventsAfterSel(0),
- fHistTrials(0),
- fHistXsection(0),
- fHistEvents(0),
+ fUseCellsToMatch(kFALSE),
+ fMinJetMCPt(1),
+ fHistoType(0),
+ fDeltaPtAxis(0),
+ fDeltaEtaDeltaPhiAxis(0),
+ fNEFAxis(0),
+ fZAxis(0),
+ fIsJet1Rho(kFALSE),
+ fIsJet2Rho(kFALSE),
+ fHistJets1(0),
+ fHistJets2(0),
+ fHistMatching(0),
fHistJets1PhiEta(0),
fHistJets1PtArea(0),
fHistJets1CorrPtArea(0),
- fHistLeadingJets1PtArea(0),
- fHistLeadingJets1CorrPtArea(0),
+ fHistJets1NEFvsPt(0),
+ fHistJets1CEFvsCEFPt(0),
+ fHistJets1ZvsPt(0),
fHistJets2PhiEta(0),
fHistJets2PtArea(0),
fHistJets2CorrPtArea(0),
- fHistLeadingJets2PtArea(0),
- fHistLeadingJets2CorrPtArea(0),
- fHistJets2PhiEtaAcceptance(0),
- fHistJets2PtAreaAcceptance(0),
- fHistJets2CorrPtAreaAcceptance(0),
- fHistLeadingJets2PtAreaAcceptance(0),
- fHistLeadingJets2CorrPtAreaAcceptance(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),
- fHistDeltaEtaPhi(0),
fHistDeltaPtvsJet1Pt(0),
fHistDeltaPtvsJet2Pt(0),
- fHistDeltaPtvsMatchingLevel(0),
- fHistDeltaCorrPtvsJet1Pt(0),
- fHistDeltaCorrPtvsJet2Pt(0),
- fHistDeltaCorrPtvsMatchingLevel(0),
- fHistNonMatchedJets1PtArea(0),
- fHistNonMatchedJets2PtArea(0),
- fHistNonMatchedJets1CorrPtArea(0),
- fHistNonMatchedJets2CorrPtArea(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),
- fHistMissedJets2PtArea(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.
// Destructor
}
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
-{
- //
- // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
- // Get the pt hard bin from the file path
- // This is to called in Notify and should provide the path to the AOD/ESD file
- // (Partially copied from AliAnalysisHelperJetTasks)
-
- TString file(currFile);
- fXsec = 0;
- fTrials = 1;
-
- if(file.Contains("root_archive.zip#")){
- Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
- Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
- Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
- file.Replace(pos+1,pos2-pos1,"");
- }
- else {
- // not an archive take the basename....
- file.ReplaceAll(gSystem->BaseName(file.Data()),"");
- }
- Printf("%s",file.Data());
-
- // Get the pt hard bin
- TString strPthard(file);
- strPthard.Remove(strPthard.Last('/'));
- strPthard.Remove(strPthard.Last('/'));
- strPthard.Remove(0,strPthard.Last('/')+1);
- if (strPthard.IsDec())
- pthard = strPthard.Atoi();
- else
- AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
-
- TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
- if(!fxsec){
- // next trial fetch the histgram file
- fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
- if(!fxsec){
- // not a severe condition but inciate that we have no information
- return kFALSE;
- }
- else{
- // find the tlist we want to be independtent of the name so use the Tkey
- TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
- if(!key){
- fxsec->Close();
- return kFALSE;
- }
- TList *list = dynamic_cast<TList*>(key->ReadObj());
- if(!list){
- fxsec->Close();
- return kFALSE;
- }
- fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
- fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
- fxsec->Close();
- }
- } // no tree pyxsec.root
- else {
- TTree *xtree = (TTree*)fxsec->Get("Xsection");
- if(!xtree){
- fxsec->Close();
- return kFALSE;
- }
- UInt_t ntrials = 0;
- Double_t xsection = 0;
- xtree->SetBranchAddress("xsection",&xsection);
- xtree->SetBranchAddress("ntrials",&ntrials);
- xtree->GetEntry(0);
- fTrials = ntrials;
- fXsec = xsection;
- fxsec->Close();
- }
- return kTRUE;
-}
-
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::UserNotify()
-{
- TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
- if (!tree) {
- AliError(Form("%s - UserNotify: No current tree!",GetName()));
- return kFALSE;
- }
-
- Float_t xsection = 0;
- Float_t trials = 0;
- Int_t pthard = 0;
-
- TFile *curfile = tree->GetCurrentFile();
- if (!curfile) {
- AliError(Form("%s - UserNotify: No current file!",GetName()));
- return kFALSE;
- }
-
- PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
-
- fHistTrials->SetBinContent(pthard + 1, fHistTrials->GetBinContent(pthard + 1) + trials);
- fHistXsection->SetBinContent(pthard + 1, fHistXsection->GetBinContent(pthard + 1) + xsection);
- fHistEvents->SetBinContent(pthard + 1, fHistEvents->GetBinContent(pthard + 1) + 1);
-
- return kTRUE;
-}
//________________________________________________________________________
-void AliJetResponseMaker::UserCreateOutputObjects()
+void AliJetResponseMaker::AllocateTH2()
{
- // Create user objects.
-
- AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
-
- fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
- fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
- fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
- fOutput->Add(fHistTrialsAfterSel);
-
- fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
- fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
- fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
- fOutput->Add(fHistEventsAfterSel);
-
- fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
- fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
- fHistTrials->GetYaxis()->SetTitle("trials");
- fOutput->Add(fHistTrials);
-
- fHistXsection = new TH1F("fHistXsection", "fHistXsection", 11, 0, 11);
- fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
- fHistXsection->GetYaxis()->SetTitle("xsection");
- fOutput->Add(fHistXsection);
-
- fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
- fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
- fHistEvents->GetYaxis()->SetTitle("total events");
- fOutput->Add(fHistEvents);
-
- const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
- const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
-
- for (Int_t i = 1; i < 12; i++) {
- fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
- fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
-
- fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
- fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
- fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
- }
+ // 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", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+ 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);
-
- fHistLeadingJets1PtArea = new TH2F("fHistLeadingJets1PtArea", "fHistLeadingJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 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()) {
- fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+
+ 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);
-
- fHistLeadingJets1CorrPtArea = new TH2F("fHistLeadingJets1CorrPtArea", "fHistLeadingJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistLeadingJets1CorrPtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
- fHistLeadingJets1CorrPtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets1CorrPtArea);
}
+ 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);
- fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
- fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
- fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistJets2PhiEtaAcceptance);
-
- fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+ 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);
- fHistLeadingJets2PtArea = new TH2F("fHistLeadingJets2PtArea", "fHistLeadingJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
- fHistLeadingJets2PtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistLeadingJets2PtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2PtArea);
-
- fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
- fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistJets2PtAreaAcceptance);
-
- fHistLeadingJets2PtAreaAcceptance = new TH2F("fHistLeadingJets2PtAreaAcceptance", "fHistLeadingJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
- fHistLeadingJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistLeadingJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistLeadingJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2PtAreaAcceptance);
-
- if (!fRho2Name.IsNull()) {
- fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+ 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);
-
- fHistLeadingJets2CorrPtArea = new TH2F("fHistLeadingJets2CorrPtArea", "fHistLeadingJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistLeadingJets2CorrPtArea->GetXaxis()->SetTitle("area");
- fHistLeadingJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
- fHistLeadingJets2CorrPtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2CorrPtArea);
-
- fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
- fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistJets2CorrPtAreaAcceptance);
-
- fHistLeadingJets2CorrPtAreaAcceptance = new TH2F("fHistLeadingJets2CorrPtAreaAcceptance", "fHistLeadingJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistLeadingJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
- fHistLeadingJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
- fHistLeadingJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistLeadingJets2CorrPtAreaAcceptance);
}
+ 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}");
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}");
fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
- fHistDeltaEtaPhi = new TH2F("fHistDeltaEtaPhi", "fHistDeltaEtaPhi", 200, -1, 1, 250, -1.6, 4.8);
- fHistDeltaEtaPhi->GetXaxis()->SetTitle("#Delta#eta");
- fHistDeltaEtaPhi->GetYaxis()->SetTitle("#Delta#phi");
- fHistDeltaEtaPhi->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaEtaPhi);
-
- fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+ 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 = 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);
- fHistDeltaPtvsMatchingLevel = new TH2F("fHistDeltaPtvsMatchingLevel", "fHistDeltaPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistDeltaPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
- fHistDeltaPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T} (GeV/c)");
- fHistDeltaPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaPtvsMatchingLevel);
-
- if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
- fHistDeltaCorrPtvsJet1Pt = new TH2F("fHistDeltaCorrPtvsJet1Pt", "fHistDeltaCorrPtvsJet1Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistDeltaCorrPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
- fHistDeltaCorrPtvsJet1Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
- fHistDeltaCorrPtvsJet1Pt->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaCorrPtvsJet1Pt);
-
- fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
- fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
- fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
-
- fHistDeltaCorrPtvsMatchingLevel = new TH2F("fHistDeltaCorrPtvsMatchingLevel", "fHistDeltaCorrPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
- fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
- fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaCorrPtvsMatchingLevel);
- }
-
- fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
- fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
- fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
- fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistNonMatchedJets1PtArea);
-
- fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
- fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
- fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistNonMatchedJets2PtArea);
-
- if (!fRhoName.IsNull()) {
- fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
- fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
- fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistNonMatchedJets1CorrPtArea);
- }
-
- if (!fRho2Name.IsNull()) {
- fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
- fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistNonMatchedJets2CorrPtArea);
- }
+ 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->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistJet1PtvsJet2Pt);
- if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
- if (fRhoName.IsNull())
- fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
- else if (fRho2Name.IsNull())
- fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+ 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
- fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+ 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);
}
- fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
- fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
- fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
- fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
- fOutput->Add(fHistMissedJets2PtArea);
-
- PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
+ 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);
+ }
}
//________________________________________________________________________
-Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
-{
- // Return true if jet is accepted.
+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 (jet->Pt() <= fJetPtCut)
- return kFALSE;
- if (jet->Area() <= fJetAreaCut)
- return kFALSE;
+ if (fZAxis) {
+ title[dim] = "Z";
+ nbins[dim] = fNbins/4;
+ min[dim] = 0;
+ max[dim] = 1.2;
+ dim++;
+ }
- return kTRUE;
-}
+ title[dim] = "p_{T,particle}^{leading} (GeV/c)";
+ nbins[dim] = 120;
+ min[dim] = 0;
+ max[dim] = 120;
+ dim++;
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
-{
- // Accept jet with a bias.
+ Int_t dim1 = dim, dim2 = dim;
- if (fLeadingHadronType == 0) {
- if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
+ if (fIsJet1Rho) {
+ title[dim1] = "p_{T}^{corr}";
+ nbins[dim1] = fNbins*2;
+ min[dim1] = -250;
+ max[dim1] = 250;
+ dim1++;
}
- else if (fLeadingHadronType == 1) {
- if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
+
+ if (fIsEmbedded) {
+ title[dim1] = "p_{T}^{MC}";
+ nbins[dim1] = fNbins;
+ min[dim1] = 0;
+ max[dim1] = 250;
+ dim1++;
}
- else {
- if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
+
+ 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++;
}
- return kTRUE;
+ 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::ExecOnce()
+void AliJetResponseMaker::UserCreateOutputObjects()
{
- // Execute once.
+ // Create user objects.
- 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;
- }
- }
+ AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
- 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;
- }
- }
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- if (fAreCollections2MC) {
- fTracks2Map = dynamic_cast<TH1*>(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 TH1I("tracksMap","tracksMap",9999,0,1);
- for (Int_t i = 0; i < 9999; i++) {
- fTracks2Map->SetBinContent(i,i);
- }
- }
- }
- }
+ if (!jets1 || !jets2) return;
+
+ if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
+ else fIsJet1Rho = kTRUE;
+ if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
+ else fIsJet2Rho = kTRUE;
+
+ if (fHistoType==0)
+ AllocateTH2();
+ else
+ AllocateTHnSparse();
- if (!fCalo2Name.IsNull() && !fCaloClusters2) {
- fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
- if (!fCaloClusters2) {
- AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
+ 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;
- } 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;
- }
+
+ 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()));
}
- }
- 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;
+ 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);
}
}
-
- if (fJet2MinEta == -999)
- fJet2MinEta = fJetMinEta - fJetRadius;
- if (fJet2MaxEta == -999)
- fJet2MaxEta = fJetMaxEta + fJetRadius;
- if (fJet2MinPhi == -999)
- fJet2MinPhi = fJetMinPhi - fJetRadius;
- if (fJet2MaxPhi == -999)
- fJet2MaxPhi = fJetMaxPhi + fJetRadius;
-
- AliAnalysisTaskEmcalJet::ExecOnce();
}
//________________________________________________________________________
-Bool_t AliJetResponseMaker::IsEventSelected()
+void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2,
+ Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2,
+ Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
{
- // Check if event is selected
+ if (fHistoType==1) {
+ Double_t contents[20]={0};
+
+ for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
+ TString title(fHistMatching->GetAxis(i)->GetTitle());
+ if (title=="p_{T,1}")
+ contents[i] = Pt1;
+ else if (title=="p_{T,2}")
+ contents[i] = Pt2;
+ else if (title=="A_{jet,1}")
+ contents[i] = A1;
+ else if (title=="A_{jet,2}")
+ contents[i] = A2;
+ else if (title=="distance")
+ contents[i] = d;
+ else if (title=="CE1")
+ contents[i] = CE1;
+ else if (title=="CE2")
+ contents[i] = CE2;
+ else if (title=="#deltaA_{jet}")
+ contents[i] = A1-A2;
+ else if (title=="#deltap_{T}")
+ contents[i] = Pt1-Pt2;
+ else if (title=="#delta#eta")
+ contents[i] = Eta1-Eta2;
+ else if (title=="#delta#phi")
+ contents[i] = Phi1-Phi2;
+ else if (title=="p_{T,1}^{corr}")
+ contents[i] = CorrPt1;
+ else if (title=="p_{T,2}^{corr}")
+ contents[i] = CorrPt2;
+ else if (title=="#deltap_{T}^{corr}")
+ contents[i] = CorrPt1-CorrPt2;
+ else if (title=="p_{T,1}^{MC}")
+ contents[i] = MCPt1;
+ else if (title=="#deltap_{T}^{MC}")
+ contents[i] = MCPt1-Pt2;
+ else if (title=="NEF_{1}")
+ contents[i] = NEF1;
+ else if (title=="NEF_{2}")
+ contents[i] = NEF2;
+ else if (title=="Z_{1}")
+ contents[i] = LeadingPt1/Pt1;
+ else if (title=="Z_{2}")
+ contents[i] = LeadingPt2/Pt2;
+ else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
+ contents[i] = LeadingPt1;
+ else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
+ contents[i] = LeadingPt2;
+ else
+ AliWarning(Form("Unable to fill dimension %s!",title.Data()));
+ }
- if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
- return kFALSE;
+ fHistMatching->Fill(contents);
+ }
+ else {
+ fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
+ fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
+
+ fHistDistancevsJet1Pt->Fill(d, Pt1);
+ fHistDistancevsJet2Pt->Fill(d, Pt2);
+
+ fHistDistancevsCommonEnergy1->Fill(d, CE1);
+ fHistDistancevsCommonEnergy2->Fill(d, CE2);
+ fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
+
+ fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
+
+ fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
+ fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
+
+ Double_t dpt = Pt1 - Pt2;
+ fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
+ fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
+ fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
+ fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
+
+ fHistDeltaPtvsDistance->Fill(d, dpt);
+ fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
+ fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
+
+ fHistDeltaPtvsArea1->Fill(A1, dpt);
+ fHistDeltaPtvsArea2->Fill(A2, dpt);
+
+ Double_t darea = A1 - A2;
+ fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
+
+ fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
+
+ if (fIsJet1Rho || fIsJet2Rho) {
+ Double_t dcorrpt = CorrPt1 - CorrPt2;
+ fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
+ fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
+ fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
+ fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
+ fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
+ fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
+ fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
+ fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
+ fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
+ fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
+ fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
+ }
- return AliAnalysisTaskEmcalJet::IsEventSelected();
+ if (fIsEmbedded) {
+ Double_t dmcpt = MCPt1 - Pt2;
+ fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
+ fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
+ fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
+ fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
+ fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
+ fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
+ fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
+ fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
+ fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
+ fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
+ fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
+ }
+ }
}
//________________________________________________________________________
-Bool_t AliJetResponseMaker::RetrieveEventObjects()
+void AliJetResponseMaker::ExecOnce()
{
- // Retrieve event objects.
+ // Execute once.
- if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
- return kFALSE;
+ AliAnalysisTaskEmcalJet::ExecOnce();
- if (fRho2)
- fRho2Val = fRho2->GetVal();
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
- const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
-
- if (MCEvent())
- fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
- if (fPythiaHeader) {
- Double_t pthard = fPythiaHeader->GetPtHard();
-
- for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
- if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
- break;
+ if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
+ if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
+ AliWarning("Changing matching type from MC label to same collection...");
+ fMatching = kSameCollections;
+ }
+ else {
+ AliWarning("Changing matching type from MC label to geometrical...");
+ fMatching = kGeometrical;
+ }
+ }
+ else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
+ if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
+ AliWarning("Changing matching type from same collection to MC label...");
+ fMatching = kMCLabel;
+ }
+ else {
+ AliWarning("Changing matching type from same collection to geometrical...");
+ fMatching = kGeometrical;
}
-
- fNTrials = fPythiaHeader->Trials();
}
-
- return kTRUE;
}
//________________________________________________________________________
//________________________________________________________________________
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);
- }
+ // 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;
-
- if (order) {
- jets1 = fJets2;
- jets2 = fJets;
- }
- else {
- jets1 = fJets;
- jets2 = fJets2;
- }
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- Int_t nJets1 = jets1->GetEntriesFast();
- Int_t nJets2 = jets2->GetEntriesFast();
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
- for (Int_t i = 0; i < nJets1; i++) {
+ AliEmcalJet* jet1 = 0;
+ AliEmcalJet* jet2 = 0;
- AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
-
- if (!jet1) {
- AliError(Form("Could not receive jet %d", i));
- 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 (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
- continue;
- }
+ 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::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());
- 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->GetNbinsX()-2) {
- index = fTracks2Map->GetBinContent(MClabel);
+ TLorentzVector part;
+ clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+
+ for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
+ Int_t cellId = clus->GetCellAbsId(iCell);
+ Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
+
+ Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
+ MClabel -= fMCLabelShift;
+ if (MClabel != 0) continue;
+
+ // this is not a MC particle; remove it completely
+ AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
+ totalPt1 -= part.Pt() * cellFrac;
+ d1 -= part.Pt() * cellFrac;
}
+ }
+ }
+ else {
+ for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+ AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
+ if (!clus) {
+ AliWarning(Form("Could not find cluster %d!", iClus));
+ continue;
+ }
+ TLorentzVector part;
+ clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+ Int_t MClabel = TMath::Abs(clus->GetLabel());
+ MClabel -= fMCLabelShift;
+ if (MClabel != 0) continue;
+
+ // this is not a MC particle; remove it completely
+ AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
+ totalPt1 -= part.Pt();
+ d1 -= part.Pt();
+ }
+ }
+
+ for (Int_t 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 (index2 == index) { // found common particle
- track2Found = kTRUE;
- d1 -= track->Pt();
- AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
+ 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();
- break;
}
+
+ track2Found = kTRUE;
}
- 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));
-
- if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
+
+ // 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));
- 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->GetNbinsX()-2) {
- index = fTracks2Map->GetBinContent(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());
- 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->GetNbinsX()-2) {
- index = fTracks2Map->GetBinContent(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;
}
}
}
- if (d1 <= 0 || totalPt1 < 1)
+
+ if (d1 < 0)
d1 = 0;
+
+ if (d2 < 0)
+ d2 = 0;
+
+ if (totalPt1 < 1)
+ d1 = -1;
else
d1 /= totalPt1;
- if (jet2->Pt() > 0 && d2 > 0)
- d2 /= jet2->Pt();
+ if (jet2->Pt() < 1)
+ d2 = -1;
else
- d2 = 0;
+ d2 /= jet2->Pt();
}
//________________________________________________________________________
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) {
-
- 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);
- AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(index));
- if (!clus) {
- AliWarning(Form("Could not find cluster %d!", index));
+ 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;
}
- AliVCluster *clus2 = static_cast<AliVCluster*>(fCaloClusters2->At(index2));
+ TLorentzVector part1;
+ clus1->GetMomentum(part1, const_cast<Double_t*>(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;
}
- TLorentzVector part, part2;
- clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+ Int_t ncells2 = clus2->GetNCells();
+ if (ncells2 >= maxNcells2) {
+ AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
+ continue;
+ }
+ UShort_t *cellsId2 = clus2->GetCellsAbsId();
+ Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
+
+ TLorentzVector part2;
clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
+ Double_t ptClus2 = part2.Pt();
- d1 -= part.Pt();
- d2 -= part2.Pt();
- break;
+ 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<Double_t*>(fVertex));
+ clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
+
+ d1 -= part1.Pt();
+ d2 -= part2.Pt();
+ break;
+ }
+ }
}
}
-
}
- if (jet1->Pt() > 0 && d1 > 0)
+ if (d1 < 0)
+ d1 = 0;
+
+ if (d2 < 0)
+ d2 = 0;
+
+ if (jet1->Pt() > 0)
d1 /= jet1->Pt();
else
- d1 = 0;
+ d1 = -1;
- if (jet2->Pt() > 0 && d2 > 0)
+ if (jet2->Pt() > 0)
d2 /= jet2->Pt();
else
- d2 = 0;
+ d2 = -1;
}
//________________________________________________________________________
;
}
- if (d1 > 0) {
+ if (d1 >= 0) {
if (d1 < jet1->ClosestJetDistance()) {
jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
}
}
- if (d2 > 0) {
+ if (d2 >= 0) {
if (d2 < jet2->ClosestJetDistance()) {
jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
{
// Fill histograms.
- static Int_t indexes[9999] = {-1};
+ AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+ AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
- fHistEventsAfterSel->SetBinContent(fPtHardBin + 1, fHistEventsAfterSel->GetBinContent(fPtHardBin + 1) + 1);
- fHistTrialsAfterSel->SetBinContent(fPtHardBin + 1, fHistTrialsAfterSel->GetBinContent(fPtHardBin + 1) + fNTrials);
+ if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
- GetSortedArray(indexes, fJets2, fRho2Val);
+ AliEmcalJet* jet1 = 0;
+ AliEmcalJet* jet2 = 0;
+
+ jets2->ResetCurrentID();
+ while ((jet2 = jets2->GetNextJet())) {
- const Int_t nJets2 = fJets2->GetEntriesFast();
+ AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
- Int_t naccJets2 = 0;
- Int_t naccJets2Acceptance = 0;
+ Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
+ Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
- for (Int_t i = 0; i < nJets2; i++) {
+ if (jets2->AcceptJet(jet2))
+ FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2,
+ corrpt2, jet2->MCPt(), 2);
- AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(indexes[i]));
+ jet1 = jet2->MatchedJet();
- if (!jet2) {
- AliError(Form("Could not receive jet2 %d", i));
- continue;
- }
+ if (!jet1) continue;
+ if (!jets1->AcceptJet(jet1)) continue;
+ if (jet1->MCPt() < fMinJetMCPt) continue;
- if (!AcceptJet(jet2))
- continue;
-
- if (AcceptBiasJet(jet2) &&
- (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
-
- fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
- fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
+ 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);
- if (naccJets2Acceptance < fNLeadingJets)
- fHistLeadingJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
-
- if (!fRho2Name.IsNull()) {
- fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
- if (naccJets2Acceptance < fNLeadingJets)
- fHistLeadingJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
- }
-
- naccJets2Acceptance++;
+ d = jet2->ClosestJetDistance();
}
+ else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
+ GetGeometricalMatchingLevel(jet1, jet2, d);
- if (!AcceptBiasJet2(jet2))
- continue;
-
- if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
- continue;
-
- fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
- fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
-
- if (naccJets2 < fNLeadingJets)
- fHistLeadingJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-
- if (!fRho2Name.IsNull()) {
- fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
- if (naccJets2 < fNLeadingJets)
- fHistLeadingJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
+ ce1 = jet1->ClosestJetDistance();
+ ce2 = jet2->ClosestJetDistance();
}
- naccJets2++;
-
- if (jet2->MatchedJet()) {
-
- if (!AcceptBiasJet(jet2->MatchedJet()) ||
- jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt) {
- fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
- }
- else {
- if (jet2->MatchedJet()->Pt() > fMaxBinPt)
- fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-
- Double_t d1=-1, d2=-1;
- if (jet2->GetMatchingType() == kGeometrical) {
-
- if (fAreCollections2MC && !fAreCollections1MC)
- GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
- else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
- GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
-
- fHistDistancevsCommonEnergy1->Fill(jet2->ClosestJetDistance(), d1);
- fHistDistancevsCommonEnergy2->Fill(jet2->ClosestJetDistance(), d2);
-
- fHistDistancevsJet1Pt->Fill(jet2->ClosestJetDistance(), jet2->MatchedJet()->Pt());
- fHistDistancevsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
-
- fHistCommonEnergy1vsJet1Pt->Fill(d1, jet2->MatchedJet()->Pt());
- fHistCommonEnergy2vsJet2Pt->Fill(d2, jet2->Pt());
- }
- else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
- GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d1);
-
- fHistDistancevsCommonEnergy1->Fill(d1, jet2->MatchedJet()->ClosestJetDistance());
- fHistDistancevsCommonEnergy2->Fill(d1, jet2->ClosestJetDistance());
-
- fHistDistancevsJet1Pt->Fill(d1, jet2->MatchedJet()->Pt());
- fHistDistancevsJet2Pt->Fill(d1, jet2->Pt());
-
- fHistCommonEnergy1vsJet1Pt->Fill(jet2->MatchedJet()->ClosestJetDistance(), jet2->MatchedJet()->Pt());
- fHistCommonEnergy2vsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
- }
-
- Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
- Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
- fHistDeltaEtaPhi->Fill(deta, dphi, jet2->Pt());
-
- Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
- fHistDeltaPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt);
- fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
- fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
-
- fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(), jet2->Pt() / jet2->MatchedJet()->Pt());
- fHistJet1PtOverJet2PtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), jet2->MatchedJet()->Pt() / jet2->Pt());
+ Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
+ Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
- fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
-
- if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
- dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
- fHistDeltaCorrPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt);
- fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
- fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
- fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
- }
- }
- }
- else {
- fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
- fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
-
- if (!fRho2Name.IsNull())
- fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->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);
}
- GetSortedArray(indexes, fJets, fRhoVal);
-
- const Int_t nJets1 = fJets->GetEntriesFast();
- Int_t naccJets1 = 0;
-
- for (Int_t i = 0; i < nJets1; i++) {
-
- AliDebug(2,Form("Processing jet %d", i));
- AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(indexes[i]));
-
- if (!jet1) {
- AliError(Form("Could not receive jet %d", i));
- continue;
- }
-
- if (!AcceptJet(jet1))
- continue;
-
- if (!AcceptBiasJet(jet1))
- continue;
-
- if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
- continue;
-
- if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
- continue;
-
- if (!jet1->MatchedJet()) {
- fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
- if (!fRhoName.IsNull())
- fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
- }
+ jets1->ResetCurrentID();
+ while ((jet1 = jets1->GetNextAcceptJet())) {
+ if (jet1->MCPt() < fMinJetMCPt) continue;
+ AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
- fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
- fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
+ Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
+ Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
- if (naccJets1 < fNLeadingJets)
- fHistLeadingJets1PtArea->Fill(jet1->Area(), jet1->Pt());
-
- if (!fRhoName.IsNull()) {
- fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
-
- if (naccJets1 < fNLeadingJets)
- fHistLeadingJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
- }
-
- naccJets1++;
+ 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