From 51a8c8fd92f151dd589ce85868cca2cd7da6ca70 Mon Sep 17 00:00:00 2001 From: jklein Date: Mon, 11 Aug 2014 16:59:56 +0200 Subject: [PATCH] - adding jet v2 task using EMCAL jet framework Jason Mueller (CERN summer student 2014) & Alice Ohlson --- PWGJE/CMakelibPWGJEEMCALJetTasks.pkg | 1 + .../UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx | 736 ++++++++++++++++++ .../UserTasks/AliAnalysisTaskEmcalJetv2QA.h | 107 +++ .../macros/AddTaskEmcalJetv2QA.C | 93 +++ PWGJE/PWGJEEMCALJetTasksLinkDef.h | 1 + 5 files changed, 938 insertions(+) create mode 100644 PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx create mode 100644 PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.h create mode 100644 PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetv2QA.C diff --git a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg index ff2849d5b3f..001d5fef758 100644 --- a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg +++ b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg @@ -76,6 +76,7 @@ set ( SRCS EMCALJetTasks/UserTasks/AliAnalysisTaskJetMassResponseDet.cxx EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx + EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassBase.cxx EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx new file mode 100644 index 00000000000..689bd876b0e --- /dev/null +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.cxx @@ -0,0 +1,736 @@ +// Jet v2 task using QA method, based on jet sample task (S.Aiola). +// +// Authors: Jason Mueller (CERN summer student 2014) & Alice Ohlson + + +#include +#include +#include +#include +#include +#include +#include + +#include "AliVCluster.h" +#include "AliAODCaloCluster.h" +#include "AliESDCaloCluster.h" +#include "AliVTrack.h" +#include "AliEmcalJet.h" +#include "AliRhoParameter.h" +#include "AliLog.h" +#include "AliJetContainer.h" +#include "AliParticleContainer.h" +#include "AliClusterContainer.h" +#include "AliPicoTrack.h" + +#include "AliAnalysisTaskEmcalJetv2QA.h" + +ClassImp(AliAnalysisTaskEmcalJetv2QA) + +//________________________________________________________________________ +AliAnalysisTaskEmcalJetv2QA::AliAnalysisTaskEmcalJetv2QA() : +AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetv2QA", kTRUE), + nCentBins(0), + nJetPtBins(0), + centBins(0x0), + jetPtBins(0x0), + fJetv2(0x0), + doPtWeight(kFALSE), + fHistTracksPt(0), + fHistClustersPt(0), + fHistLeadingJetPt(0), + fHistLeadingJetPtCorr(0), + fHistJetsPhiEta(0), + fHistJetsPtArea(0), + fHistJetsPtLeadHad(0), + fHistJetsCorrPtArea(0), + fHistPtDEtaDPhiTrackClus(0), + fHistPtDEtaDPhiClusTrack(0), + fDPhiJet(0), + fDPhiJetPythia(0), + fDPhiEP(0), + hGx(0), + hGy2(0), + hGxGy2(0), + hGy4(0), + hGx2(0), + hGx2Gy2(0), + hGxGy4(0), + hGy6(0), + hGx2Gy4(0), + hGxGy6(0), + hGy8(0), + hGy(0), + hN(0), + htv2std(0), + htjv2std(0), + htj2v2std(0), + hV0jv2std(0), + htdPsi(0), + htjdPsi(0), + htj2dPsi(0), + hV0jdPsi(0), + hAx(0), + hAxDijet(0), + hQx(0), + hQy(0), + hEventData(0), + hNTracks(0), + hNTracksCent(0), + hGxTracks(0), + hGyTracks(0), + hGy2Tracks(0), + hGxGy2Tracks(0), + hGy4Tracks(0), + htEPRes(0), + htjEPRes(0), + htj2EPRes(0), + fJetsCont(0), + fTracksCont(0), + fCaloClustersCont(0) +{ + // Default constructor. + + cout << endl; + cout << "*******************************************************" << endl; + cout << "*** AliAnalysisTaskEmcalJetv2QA constructor! ***" << endl; + cout << "*******************************************************" << endl; + cout << endl; + + + + SetMakeGeneralHistograms(kTRUE); +} + +//________________________________________________________________________ +AliAnalysisTaskEmcalJetv2QA::AliAnalysisTaskEmcalJetv2QA(const char *name) : + AliAnalysisTaskEmcalJet(name, kTRUE), + nCentBins(0), + nJetPtBins(0), + centBins(0x0), + jetPtBins(0x0), + fJetv2(0x0), + doPtWeight(kFALSE), + fHistTracksPt(0), + fHistClustersPt(0), + fHistLeadingJetPt(0), + fHistLeadingJetPtCorr(0), + fHistJetsPhiEta(0), + fHistJetsPtArea(0), + fHistJetsPtLeadHad(0), + fHistJetsCorrPtArea(0), + fHistPtDEtaDPhiTrackClus(0), + fHistPtDEtaDPhiClusTrack(0), + fDPhiJet(0), + fDPhiJetPythia(0), + fDPhiEP(0), + hGx(0), + hGy2(0), + hGxGy2(0), + hGy4(0), + hGx2(0), + hGx2Gy2(0), + hGxGy4(0), + hGy6(0), + hGx2Gy4(0), + hGxGy6(0), + hGy8(0), + hGy(0), + hN(0), + htv2std(0), + htjv2std(0), + htj2v2std(0), + hV0jv2std(0), + htdPsi(0), + htjdPsi(0), + htj2dPsi(0), + hV0jdPsi(0), + hAx(0), + hAxDijet(0), + hQx(0), + hQy(0), + hEventData(0), + hNTracks(0), + hNTracksCent(0), + hGxTracks(0), + hGyTracks(0), + hGy2Tracks(0), + hGxGy2Tracks(0), + hGy4Tracks(0), + htEPRes(0), + htjEPRes(0), + htj2EPRes(0), + fJetsCont(0), + fTracksCont(0), + fCaloClustersCont(0) +{ + // Standard constructor. + + + SetMakeGeneralHistograms(kTRUE); +} + +//________________________________________________________________________ +AliAnalysisTaskEmcalJetv2QA::~AliAnalysisTaskEmcalJetv2QA() +{ + // Destructor. +} + +//________________________________________________________________________ +void AliAnalysisTaskEmcalJetv2QA::UserCreateOutputObjects() +{ + // Create user output. + + AliAnalysisTaskEmcalJet::UserCreateOutputObjects(); + + fJetsCont = GetJetContainer(0); + if(fJetsCont) { //get particles and clusters connected to jets + fTracksCont = fJetsCont->GetParticleContainer(); + fCaloClustersCont = fJetsCont->GetClusterContainer(); + } else { //no jets, just analysis tracks and clusters + fTracksCont = GetParticleContainer(0); + fCaloClustersCont = GetClusterContainer(0); + } + + nCentBins = 6; // set number of centrality bins + nJetPtBins = 4; // set number of jetPtBins + centBins = new Double_t[nCentBins+1]; + jetPtBins = new Double_t[nJetPtBins+1]; + centBins[0] = 0.; centBins[1] = 5.; centBins[2] = 10.; centBins[3] = 20.; centBins[4] = 30.; centBins[5] = 40.; centBins[6] = 50.; // set edges of bins + jetPtBins[0] = 50.; jetPtBins[1] = 70.; jetPtBins[2] = 90.; jetPtBins[3] = 120.; jetPtBins[4] = 200.; // set edges of bins + + fTracksCont->SetClassName("AliVTrack"); + fCaloClustersCont->SetClassName("AliVCluster"); + + fHistTracksPt = new TH1F("fHistTracksPt","fHistTracksPt", fNbins / 2, fMinBinPt, fMaxBinPt / 2); + fHistTracksPt->GetXaxis()->SetTitle("p_{T,track} (GeV/c)"); + fHistTracksPt->GetYaxis()->SetTitle("counts"); + fOutput->Add(fHistTracksPt); + + fHistClustersPt = new TH1F("fHistClustersPt", "fHistClustersPt", fNbins / 2, fMinBinPt, fMaxBinPt / 2); + fHistClustersPt->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)"); + fHistClustersPt->GetYaxis()->SetTitle("counts"); + fOutput->Add(fHistClustersPt); + + fHistLeadingJetPt = new TH1F("fHistLeadingJetPt", "fHistLeadingJetPt", fNbins, fMinBinPt, fMaxBinPt); + fHistLeadingJetPt->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)"); + fHistLeadingJetPt->GetYaxis()->SetTitle("counts"); + fOutput->Add(fHistLeadingJetPt); + + fHistLeadingJetPtCorr = new TH1F("fHistLeadingJetPtCorr", "fHistLeadingJetPtCorr", fNbins, fMinBinPt, fMaxBinPt); + fHistLeadingJetPtCorr->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)"); + fHistLeadingJetPtCorr->GetYaxis()->SetTitle("counts"); + fOutput->Add(fHistLeadingJetPtCorr); + + fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200); + fHistJetsPhiEta->GetXaxis()->SetTitle("#eta"); + fHistJetsPhiEta->GetYaxis()->SetTitle("#phi"); + fOutput->Add(fHistJetsPhiEta); + + fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3); + fHistJetsPtArea->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)"); + fHistJetsPtArea->GetYaxis()->SetTitle("area"); + fOutput->Add(fHistJetsPtArea); + + fHistJetsPtLeadHad = new TH2F("fHistJetsPtLeadHad", "fHistJetsPtLeadHad", fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2); + fHistJetsPtLeadHad->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)"); + fHistJetsPtLeadHad->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)"); + fHistJetsPtLeadHad->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJetsPtLeadHad); + + fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3); + fHistJetsCorrPtArea->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]"); + fHistJetsCorrPtArea->GetYaxis()->SetTitle("area"); + fOutput->Add(fHistJetsCorrPtArea); + + fHistPtDEtaDPhiTrackClus = new TH3F("fHistPtDEtaDPhiTrackClus","fHistPtDEtaDPhiTrackClus;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",100,0.,100.,100,-0.1,0.1,100,-0.1,0.1); + fOutput->Add(fHistPtDEtaDPhiTrackClus); + + fHistPtDEtaDPhiClusTrack = new TH3F("fHistPtDEtaDPhiClusTrack","fHistPtDEtaDPhiClusTrack;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",100,0.,100.,100,-0.1,0.1,100,-0.1,0.1); + fOutput->Add(fHistPtDEtaDPhiClusTrack); + + fDPhiJet = new TH1F("fDPhiJet","fDPhiJet",90, -TMath::Pi()/3, 5*TMath::Pi()/3); + fDPhiJet->GetXaxis()->SetTitle("#Delta#varphi"); + fDPhiJet->GetYaxis()->SetTitle("counts"); + fOutput->Add(fDPhiJet); + + fDPhiJetPythia = new TH1F("fDPhiJetPythia","fDPhiJetPythia",90, -TMath::Pi()/3, 5*TMath::Pi()/3); + fDPhiJetPythia->GetXaxis()->SetTitle("#Delta#varphi"); + fDPhiJetPythia->GetYaxis()->SetTitle("counts"); + fOutput->Add(fDPhiJetPythia); + + fDPhiEP = new TH1F("fDPhiEP","fDPhiEP",90, 0, 2*TMath::Pi()); + fDPhiEP->GetXaxis()->SetTitle("#Delta#varphi"); + fDPhiEP->GetYaxis()->SetTitle("counts"); + fOutput->Add(fDPhiEP); + + hGx = new TH2D("hGx", "Gx v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGx->GetXaxis()->SetTitle("Centrality (%)"); + hGx->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGx); + + hGy2 = new TH2D("hGy2", "Gy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGy2->GetXaxis()->SetTitle("Centrality (%)"); + hGy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGy2); + + hGxGy2 = new TH2D("hGxGy2", "GxGy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGxGy2->GetXaxis()->SetTitle("Centrality (%)"); + hGxGy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGxGy2); + + hGy4 = new TH2D("hGy4", "Gy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGy4->GetXaxis()->SetTitle("Centrality (%)"); + hGy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGy4); + + hGx2 = new TH2D("hGx2", "Gx2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGx2->GetXaxis()->SetTitle("Centrality (%)"); + hGx2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGx2); + + hGx2Gy2 = new TH2D("hGx2Gy2", "Gx2Gy2 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGx2Gy2->GetXaxis()->SetTitle("Centrality (%)"); + hGx2Gy2->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGx2Gy2); + + hGxGy4 = new TH2D("hGxGy4", "GxGy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGxGy4->GetXaxis()->SetTitle("Centrality (%)"); + hGxGy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGxGy4); + + hGy6 = new TH2D("hGy6", "Gy6 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGy6->GetXaxis()->SetTitle("Centrality (%)"); + hGy6->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGy6); + + hGx2Gy4 = new TH2D("hGx2Gy4", "Gx2Gy4 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGx2Gy4->GetXaxis()->SetTitle("Centrality (%)"); + hGx2Gy4->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGx2Gy4); + + hGxGy6 = new TH2D("hGxGy6", "GxGy6 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGxGy6->GetXaxis()->SetTitle("Centrality (%)"); + hGxGy6->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGxGy6); + + hGy8 = new TH2D("hGy8", "Gy8 v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGy8->GetXaxis()->SetTitle("Centrality (%)"); + hGy8->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGy8); + + hGy = new TH2D("hGy", "Gy v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hGy->GetXaxis()->SetTitle("Centrality (%)"); + hGy->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hGy); + + hN = new TH2D("hN", "N v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hN->GetXaxis()->SetTitle("Centrality (%)"); + hN->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hN); + + htv2std = new TH2D("htv2std", "v2std v Centrality v JetPt w/o jets", nCentBins, centBins, nJetPtBins, jetPtBins); + htv2std->GetXaxis()->SetTitle("Centrality (%)"); + htv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(htv2std); + + htjv2std = new TH2D("htjv2std", "v2std v Centrality v JetPt w/ jets", nCentBins, centBins, nJetPtBins, jetPtBins); + htjv2std->GetXaxis()->SetTitle("Centrality (%)"); + htjv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(htjv2std); + + htj2v2std = new TH2D("htj2v2std", "v2std v Centrality v JetPt w/ trackPt < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins); + htj2v2std->GetXaxis()->SetTitle("Centrality (%)"); + htj2v2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(htj2v2std); + + hV0jv2std = new TH2D("hV0jv2std", "v2std v Centrality v JetPt", nCentBins, centBins, nJetPtBins, jetPtBins); + hV0jv2std->GetXaxis()->SetTitle("Centrality (%)"); + hV0jv2std->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hV0jv2std); + + Int_t ndpsibins = 100; + Double_t dpsibins[101]; + for(Int_t t = 0; t < 101; t++) dpsibins[t] = TMath::Pi()*t/50.; + + htdPsi = new TH3D("htdPsi", "JetAxis - EventPlane w/o jets", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins); + htdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}"); + fOutput->Add(htdPsi); + + htjdPsi = new TH3D("htjdPsi", "JetAxis - EventPlane w/ jets", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins); + htjdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}"); + fOutput->Add(htjdPsi); + + htj2dPsi = new TH3D("htj2dPsi", "JetAxis - EventPlane w/ trackPt < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins); + htj2dPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}"); + fOutput->Add(htj2dPsi); + + hV0jdPsi = new TH3D("hV0jdPsi", "JetAxis - EventPlane", nCentBins, centBins, nJetPtBins, jetPtBins, ndpsibins, dpsibins); + hV0jdPsi->GetZaxis()->SetTitle("#Psi_{jet} - #Psi_{EP}"); + fOutput->Add(hV0jdPsi); + + hQx = new TH2D("hQx", "Qx Distribution in EP frame", 100, -0.3, 0.3, nCentBins, centBins); + hQx->GetXaxis()->SetTitle("Qx"); + hQx->GetYaxis()->SetTitle("Centrality (%)"); + fOutput->Add(hQx); + + hQy = new TH2D("hQy", "Qy Distribution in EP frame", 100, -0.3, 0.3, nCentBins, centBins); + hQy->GetXaxis()->SetTitle("Qy"); + hQy->GetYaxis()->SetTitle("Centrality (%)"); + fOutput->Add(hQy); + + hAx = new TH2D("hAx", "Ax Distribution in EP frame w/o Dijets", 100, -35, 70, nJetPtBins, jetPtBins); + hAx->GetXaxis()->SetTitle("Ax"); + hAx->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hAx); + + hAxDijet = new TH2D("hAxDijet", "Ax Distribution in EP frame w/ Dijets", 100, -35, 70, nJetPtBins, jetPtBins); + hAxDijet->GetXaxis()->SetTitle("Ax"); + hAxDijet->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(hAxDijet); + + hEventData = new TH1F("hEventData","Events Kept and Discarded", 9, 0, 9); + hEventData->GetYaxis()->SetTitle("counts"); + fOutput->Add(hEventData); + + hNTracks = new TH2F("hNTracks","Number of Tracks Per Event", 100, 0, 3000, 3, 0, 3); + hNTracks->GetXaxis()->SetTitle("# tracks"); + fOutput->Add(hNTracks); + + hNTracksCent = new TH2D("hNTracksCent", "NTracks by centrality", 100, 0, 3000, (Int_t)centBins[nCentBins], centBins[0], centBins[nCentBins]); + hNTracksCent->GetXaxis()->SetTitle("# tracks"); + hNTracksCent->GetYaxis()->SetTitle("Centrality (%)"); + fOutput->Add(hNTracksCent); + + hGxTracks = new TH2D("hGxTracks", "Gx by NTracks", 200, -200, 200, 100, 0, 3000); + hGxTracks->GetXaxis()->SetTitle("Gx"); + hGxTracks->GetYaxis()->SetTitle("# tracks"); + fOutput->Add(hGxTracks); + + hGyTracks = new TH2D("hGyTracks", "Gy by NTracks", 200, -200, 200, 100, 0, 3000); + hGyTracks->GetXaxis()->SetTitle("Gy"); + hGyTracks->GetYaxis()->SetTitle("# tracks"); + fOutput->Add(hGyTracks); + + hGy2Tracks = new TH2D("hGy2Tracks", "Gy2 by NTracks", 100, 0, 20000, 100, 0, 3000); + hGy2Tracks->GetXaxis()->SetTitle("Gy2"); + hGy2Tracks->GetYaxis()->SetTitle("# tracks"); + fOutput->Add(hGy2Tracks); + + hGxGy2Tracks = new TH2D("hGxGy2Tracks", "GxGy2 by NTracks", 100, -100000, 100000, 100, 0, 3000); + hGxGy2Tracks->GetXaxis()->SetTitle("GxGy2"); + hGxGy2Tracks->GetYaxis()->SetTitle("# tracks"); + fOutput->Add(hGxGy2Tracks); + + hGy4Tracks = new TH2D("hGy4Tracks", "Gy4 by NTracks", 100, 0, 100000000, 100, 0, 3000); + hGy4Tracks->GetXaxis()->SetTitle("Gy4"); + hGy4Tracks->GetYaxis()->SetTitle("# tracks"); + fOutput->Add(hGy4Tracks); + + htEPRes = new TH2D("htEPRes", "EP Resolution w/o Jets", nCentBins, centBins, nJetPtBins, jetPtBins); + htEPRes->GetXaxis()->SetTitle("Centrality (%)"); + htEPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(htEPRes); + + htjEPRes = new TH2D("htjEPRes", "EP Resolution w/ Jets", nCentBins, centBins, nJetPtBins, jetPtBins); + htjEPRes->GetXaxis()->SetTitle("Centrality (%)"); + htjEPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(htjEPRes); + + htj2EPRes = new TH2D("htj2EPRes", "EP Resolution w/ trackPT < 2 GeV", nCentBins, centBins, nJetPtBins, jetPtBins); + htj2EPRes->GetXaxis()->SetTitle("Centrality (%)"); + htj2EPRes->GetYaxis()->SetTitle("Leading Jet pT (GeV/c)"); + fOutput->Add(htj2EPRes); + + PostData(1, fOutput); // Post data for ALL output slots > 0 here. +} + +//________________________________________________________________________ +Bool_t AliAnalysisTaskEmcalJetv2QA::FillHistograms() +{ + // Fill histograms. + + return kTRUE; +} + + +//________________________________________________________________________ +void AliAnalysisTaskEmcalJetv2QA::ExecOnce() { + + AliAnalysisTaskEmcalJet::ExecOnce(); + + if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0; + if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0; + if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0; + +} + +//________________________________________________________________________ +Bool_t AliAnalysisTaskEmcalJetv2QA::Run() // this part loops over each event +{ + // Run analysis code here, if needed. It will be executed before FillHistograms(). + + Double_t jetPhi = -999; // angle of leading jet + Double_t jetPt = -999; // pt of leading jet + Double_t jetArea = -999; + Double_t trackPt = -999; + Double_t phi = -999; // track phi + Double_t dPhi = -999; // track phi - jet axis + Double_t dPhiQA = -999; // track phi - EP + Double_t tSin = 0; // used for std EP calc + Double_t tCos = 0; // used for std EP calc + Double_t jSin = 0; // used for std EP calc + Double_t jCos = 0; // used for std EP calc + Double_t tSin2 = 0; // used for std EP calc with trackPt < 2 GeV + Double_t tCos2 = 0; // used for std EP calc with trackPt < 2 GeV + Double_t qx = 0; // used for Qx distribution + Double_t qy = 0; // used for Qy distribution + Double_t ax = 0; // used for Ax distribution + Double_t tEP = -999; // EP w/o jets + Double_t tjEP = -999; // EP w/ jets + Double_t tjEP2 = -999; // EP w/ jets + Double_t dPsi = -999; // jet axis - EP + Double_t gx = 0; // used for G moment calc + Double_t gy = 0; // used for G moment calc + Int_t isDijet = 0; // if 0, no dijet. if 1, dijet. + Int_t nTracksBkgnd = 0; // used to keep track of number of tracks in background + Int_t nTracksJet =0; // used to keep track of number of tracks in jets + + if(fJetsCont) + { + + AliEmcalJet *jettest = fJetsCont->GetNextAcceptJet(0); + while(jettest) { + + fHistJetsPtArea->Fill(jettest->Pt(), jettest->Area()); + fHistJetsPhiEta->Fill(jettest->Eta(), jettest->Phi()); + + Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jettest); + fHistJetsPtLeadHad->Fill(jettest->Pt(), ptLeading); + + if (fHistJetsCorrPtArea) { + Float_t corrPt = jettest->Pt() - fJetsCont->GetRhoVal() * jettest->Area(); + fHistJetsCorrPtArea->Fill(corrPt, jettest->Area()); + } + jettest = fJetsCont->GetNextAcceptJet(); + } + + + AliEmcalJet *jet = fJetsCont->GetLeadingJet(); + if(jet) + { + jetPhi = jet->Phi(); // get leading jet phi (jet axis) + jetPt = jet->Pt(); // get leading jet pT to filter out low pT jet events + jetArea = jet->Area(); + } + } + + // event cuts + if(!fTracksCont) + { + hEventData->Fill("!fTracksCont",1); + return kFALSE; + } + if(!fJetsCont) + { + hEventData->Fill("!fJetsCont",1); + return kFALSE; + } + if(jetPt == -999) + { + hEventData->Fill("jetPt=-999",1); + return kFALSE; + } + if(jetPt < jetPtBins[0]) + { + hEventData->Fill("leadingJetPt jetPtBins[nJetPtBins]) + { + hEventData->Fill("leadingJetPt>jetPtMax",1); + return kFALSE; + } + if(fCent < centBins[0]) + { + hEventData->Fill("cent centBins[nCentBins]) + { + hEventData->Fill("cent>centMax",1); + return kFALSE; + } + hEventData->Fill("good event",1); + + AliEmcalJet *dijet = fJetsCont->GetNextAcceptJet(0); // check for dijet events + while(dijet) + { + if(dijet->Pt() > 50 && fabs(jetPhi-dijet->Phi()-TMath::Pi()) < 0.4) // loop over jets with pT>50 and exclude leading jet and check that angular separation is < 0.4 + isDijet = 1; + dijet = fJetsCont->GetNextAcceptJet(); + } + + if (fCaloClustersCont) + { + AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); + while(cluster) { + TLorentzVector nPart; + cluster->GetMomentum(nPart, fVertex); + fHistClustersPt->Fill(nPart.Pt()); + cluster = fCaloClustersCont->GetNextAcceptCluster(); + } + } + + fHistLeadingJetPt->Fill(jetPt); + fHistLeadingJetPtCorr->Fill(jetPt-fJetsCont->GetRhoVal()*jetArea); + + AliVTrack *track = static_cast(fTracksCont->GetNextAcceptParticle(0)); + while(track) + { // loop over all particles (including jet tracks) + trackPt = track->Pt(); + fHistTracksPt->Fill(trackPt); + phi = track->Phi(); // get track phi + + dPhi = phi-jetPhi; // get track phi - jet axis + if(dPhi < 0) dPhi += TMath::TwoPi(); + if(dPhi > TMath::TwoPi()) dPhi -= TMath::TwoPi(); + + dPhiQA = phi-fEPV0; // get track phi - EP + if(dPhiQA < 0) dPhiQA += TMath::TwoPi(); + if(dPhiQA > TMath::TwoPi()) dPhiQA -= TMath::TwoPi(); + fDPhiEP->Fill(dPhiQA); + + // fill jet-hadron correlation just to check if track labels make sense... + if(track->Pt()>1.) + { + Double_t dphiJet = dPhi; + if(dphiJet > 5*TMath::Pi()/3) dphiJet -= 2*TMath::Pi(); + if(dphiJet < -TMath::Pi()/3) dphiJet += 2*TMath::Pi(); + if(track->GetLabel() == 0) + fDPhiJet->Fill(dphiJet); + else + fDPhiJetPythia->Fill(dphiJet); + } + + Double_t weight = 1.; + if(doPtWeight) weight = trackPt; + + gx += weight*cos(2*dPhi); + gy += weight*sin(2*dPhi); + + if(track->GetLabel() == 0) + { // sum for std EP method + tSin += weight*sin(2*phi); // bkgnd has label = 0 + tCos += weight*cos(2*phi); + qx += weight*cos(2*dPhiQA); + qy += weight*sin(2*dPhiQA); + nTracksBkgnd += 1; + } + else + { + jSin += weight*sin(2*phi); // jets have label =/= 0 + jCos += weight*cos(2*phi); + ax += weight*cos(2*dPhi); + nTracksJet += 1; + } + + if(track->Pt() < 2.) + { // sum for std EP method w/ trackPt < 2 GeV + tSin2 += weight*sin(2*phi); + tCos2 += weight*cos(2*phi); + } + + + track = static_cast(fTracksCont->GetNextAcceptParticle()); // increment to next track + } // close loop over particles + + hNTracks->Fill(nTracksBkgnd,"Bkgnd Tracks",1); + hNTracks->Fill(nTracksJet,"Jet Tracks",1); + hNTracks->Fill(nTracksBkgnd+nTracksJet,"Total Tracks",1); + hNTracksCent->Fill(nTracksBkgnd+nTracksJet,fCent,1); + + if(nTracksBkgnd == 0) + { + hEventData->Fill("no tracks",1); + hEventData->Fill("good event",-1); + return kFALSE; + } + + Double_t v2weight = 1+2*fJetv2*cos(2*(jetPhi-fEPV0)); // set v2 weight for event + + hGx->Fill(fCent, jetPt, v2weight*gx); // fill histograms for QA method + hGy2->Fill(fCent, jetPt, v2weight*gy*gy); + hGxGy2->Fill(fCent, jetPt, v2weight*gx*gy*gy); + hGy4->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy); + hGx2->Fill(fCent, jetPt, v2weight*gx*gx); + hGx2Gy2->Fill(fCent, jetPt, v2weight*gx*gx*gy*gy); + hGxGy4->Fill(fCent, jetPt, v2weight*gx*gy*gy*gy*gy); + hGy6->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy*gy*gy); + hGx2Gy4->Fill(fCent, jetPt, v2weight*gx*gx*gy*gy*gy*gy); + hGxGy6->Fill(fCent, jetPt, v2weight*gx*gy*gy*gy*gy*gy*gy); + hGy8->Fill(fCent, jetPt, v2weight*gy*gy*gy*gy*gy*gy*gy*gy); + hGy->Fill(fCent, jetPt, v2weight*gy); + hN->Fill(fCent, jetPt, v2weight); + + hGxTracks->Fill(gx,nTracksBkgnd+nTracksJet); + hGyTracks->Fill(gy,nTracksBkgnd+nTracksJet); + hGy2Tracks->Fill(gy*gy,nTracksBkgnd+nTracksJet); + hGxGy2Tracks->Fill(gx*gy*gy,nTracksBkgnd+nTracksJet); + hGy4Tracks->Fill(gy*gy*gy*gy,nTracksBkgnd+nTracksJet); + hQx->Fill(qx/nTracksBkgnd,fCent); + hQy->Fill(qy/nTracksBkgnd,fCent); + + if(isDijet == 0) + hAx->Fill(ax,jetPt); + if(isDijet == 1) + hAxDijet->Fill(ax,jetPt); + + tEP = 0.5*atan2(tSin,tCos); // calculate EP w/o jets + tjEP = 0.5*atan2((tSin+jSin),(tCos+jCos)); // calculate EP w/ jets + tjEP2 = 0.5*atan2(tSin2,tCos2); // calculate EP w/ trackPt < 2 GeV + + htEPRes->Fill(fCent, jetPt, v2weight*cos(2*(tEP-fEPV0))); + htjEPRes->Fill(fCent, jetPt, v2weight*cos(2*(tjEP-fEPV0))); + htj2EPRes->Fill(fCent, jetPt, v2weight*cos(2*(tjEP2-fEPV0))); + + + dPsi = jetPhi-tEP; + if(dPsi < 0) dPsi += TMath::TwoPi(); + if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi(); + + htv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/o jets + htdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/o jets + + dPsi = jetPhi-tjEP; + if(dPsi < 0) dPsi += TMath::TwoPi(); + if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi(); + + htjv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/ jets + htjdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/ jets + + dPsi = jetPhi-tjEP2; + if(dPsi < 0) dPsi += TMath::TwoPi(); + if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi(); + + htj2v2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data w/ trackPt < 2 GeV + htj2dPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EP w/ trackPt < 2 GeV + + dPsi = jetPhi-fEPV0; + if(dPsi < 0) dPsi += TMath::TwoPi(); + if(dPsi > TMath::TwoPi()) dPsi -= TMath::TwoPi(); + + hV0jv2std->Fill(fCent, jetPt, v2weight*cos(2*dPsi)); // fill histogram with v2 data + hV0jdPsi->Fill(fCent,jetPt,dPsi); // fill histogram with jet axis - EPV0 + + + return kTRUE; // If return kFALSE FillHistogram() will NOT be executed. +} + +//________________________________________________________________________ +void AliAnalysisTaskEmcalJetv2QA::Terminate(Option_t *) +{ + // Called once at the end of the analysis. + if(centBins) delete [] centBins; + if(jetPtBins) delete [] jetPtBins; +} diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.h new file mode 100644 index 00000000000..923f2e70376 --- /dev/null +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetv2QA.h @@ -0,0 +1,107 @@ +// Jet v2 task using QA method, based on jet sample task (S.Aiola). +// +// Authors: Jason Mueller (CERN summer student 2014) & Alice Ohlson + +#ifndef ALIANALYSISTASKEMCALJETVNQA_H +#define ALIANALYSISTASKEMCALJETVNQA_H + +// $Id$ + +class TH1; +class TH2; +class TH3; +class AliJetContainer; +class AliParticleContainer; +class AliClusterContainer; + +#include "AliAnalysisTaskEmcalJet.h" + +class AliAnalysisTaskEmcalJetv2QA : public AliAnalysisTaskEmcalJet { + public: + + AliAnalysisTaskEmcalJetv2QA(); + AliAnalysisTaskEmcalJetv2QA(const char *name); + virtual ~AliAnalysisTaskEmcalJetv2QA(); + + void UserCreateOutputObjects(); + void Terminate(Option_t *option); + Double_t GetJetv2(){return fJetv2;}; + void SetJetv2(Double_t jetv2){fJetv2 = jetv2;}; + Bool_t GetDoPtWeight(){return doPtWeight;}; + void SetDoPtWeight(Bool_t ptWeight){doPtWeight = ptWeight;}; + + protected: + Int_t nCentBins; + Int_t nJetPtBins; + Double_t *centBins; + Double_t *jetPtBins; + void ExecOnce(); + Bool_t FillHistograms(); + Bool_t Run(); + Double_t fJetv2; + Bool_t doPtWeight; + + // General histograms + TH1F *fHistTracksPt; //!Track pt spectrum + TH1F *fHistClustersPt; //!Cluster pt spectrum + TH1F *fHistLeadingJetPt; //!Leading jet pt spectrum + TH1F *fHistLeadingJetPtCorr; //!Leading jet pt spectrum, background subtracted + TH2F *fHistJetsPhiEta; //!Phi-Eta distribution of jets + TH2F *fHistJetsPtArea; //!Jet pt vs. area + TH2F *fHistJetsPtLeadHad; //!Jet pt vs. leading hadron + TH2F *fHistJetsCorrPtArea; //!Jet pt - bkg vs. area + TH3F *fHistPtDEtaDPhiTrackClus; //!track pt, delta eta, delta phi to matched cluster + TH3F *fHistPtDEtaDPhiClusTrack; //!cluster pt, delta eta, delta phi to matched track + TH1F *fDPhiJet; //! + TH1F *fDPhiJetPythia; //! + TH1F *fDPhiEP; //! + TH2D *hGx; //! + TH2D *hGy2; //! + TH2D *hGxGy2; //! + TH2D *hGy4; //! + TH2D *hGx2; //! + TH2D *hGx2Gy2; //! + TH2D *hGxGy4; //! + TH2D *hGy6; //! + TH2D *hGx2Gy4; //! + TH2D *hGxGy6; //! + TH2D *hGy8; //! + TH2D *hGy; //! + TH2D *hN; //! + TH2D *htv2std; //! + TH2D *htjv2std; //! + TH2D *htj2v2std; //! + TH2D *hV0jv2std; //! + TH3D *htdPsi; //! + TH3D *htjdPsi; //! + TH3D *htj2dPsi; //! + TH3D *hV0jdPsi; //! + TH2D *hAx; //! + TH2D *hAxDijet; //! + TH2D *hQx; //! + TH2D *hQy; //! + TH1F *hEventData; //! + TH2F *hNTracks; //! + TH2D *hNTracksCent; //! + TH2D *hGxTracks; //! + TH2D *hGyTracks; //! + TH2D *hGy2Tracks; //! + TH2D *hGxGy2Tracks; //! + TH2D *hGy4Tracks; //! + TH2D *htEPRes; //! + TH2D *htjEPRes; //! + TH2D *htj2EPRes; //! + + AliJetContainer *fJetsCont; //!Jets + AliParticleContainer *fTracksCont; //!Tracks + AliClusterContainer *fCaloClustersCont; //!Clusters + + + + private: + AliAnalysisTaskEmcalJetv2QA(const AliAnalysisTaskEmcalJetv2QA&); // not implemented + AliAnalysisTaskEmcalJetv2QA &operator=(const AliAnalysisTaskEmcalJetv2QA&); // not implemented + + ClassDef(AliAnalysisTaskEmcalJetv2QA, 3) // jet v2QA analysis task +}; +#endif diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetv2QA.C b/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetv2QA.C new file mode 100644 index 00000000000..8ff9d226174 --- /dev/null +++ b/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJetv2QA.C @@ -0,0 +1,93 @@ +// Jet v2 task using QA method, based on jet sample task (S.Aiola). +// +// Authors: Jason Mueller (CERN summer student 2014) & Alice Ohlson + +AliAnalysisTaskEmcalJetv2QA* AddTaskEmcalJetv2QA( + const char *ntracks = "Tracks", + const char *nclusters = "CaloClusters", + const char *njets = "Jets", + const char *nrho = "Rho", + Int_t nCentBins = 1, + Double_t jetradius = 0.2, + Double_t jetptcut = 1, + Double_t jetareacut = 0.6, + const char *type = "TPC", + Int_t leadhadtype = 0, + const char *taskname = "AliAnalysisTaskEmcalJetv2QA" +) +{ + // Get the pointer to the existing analysis manager via the static access method. + //============================================================================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) + { + ::Error("AddTaskEmcalJetv2QA", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) + { + ::Error("AddTaskEmcalJetv2QA", "This task requires an input event handler"); + return NULL; + } + + //------------------------------------------------------- + // Init the task and do settings + //------------------------------------------------------- + + TString name(taskname); + if (strcmp(njets,"")) { + name += "_"; + name += njets; + } + if (strcmp(nrho,"")) { + name += "_"; + name += nrho; + } + if (strcmp(type,"")) { + name += "_"; + name += type; + } + + Printf("name: %s",name.Data()); + + AliAnalysisTaskEmcalJetv2QA* jetTask = new AliAnalysisTaskEmcalJetv2QA(name); + jetTask->SetCentRange(0.,100.); + jetTask->SetNCentBins(nCentBins); + + AliParticleContainer *trackCont = jetTask->AddParticleContainer(ntracks); + trackCont->SetClassName("AliVTrack"); + AliClusterContainer *clusterCont = jetTask->AddClusterContainer(nclusters); + + TString strType(type); + AliJetContainer *jetCont = jetTask->AddJetContainer(njets,strType,jetradius); + if(jetCont) { + jetCont->SetRhoName(nrho); + jetCont->ConnectParticleContainer(trackCont); + jetCont->ConnectClusterContainer(clusterCont); + jetCont->SetZLeadingCut(0.98,0.98); + jetCont->SetPercAreaCut(0.6); + jetCont->SetJetPtCut(jetptcut); + jetCont->SetLeadingHadronType(leadhadtype); + } + + //------------------------------------------------------- + // Final settings, pass to manager and set the containers + //------------------------------------------------------- + + mgr->AddTask(jetTask); + + // Create containers for input/output + AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ; + TString contname(name); + contname += "_histos"; + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), + TList::Class(),AliAnalysisManager::kOutputContainer, + Form("%s", AliAnalysisManager::GetCommonFileName())); + mgr->ConnectInput (jetTask, 0, cinput1 ); + mgr->ConnectOutput (jetTask, 1, coutput1 ); + + return jetTask; +} diff --git a/PWGJE/PWGJEEMCALJetTasksLinkDef.h b/PWGJE/PWGJEEMCALJetTasksLinkDef.h index 4abe11c965e..0397dea9715 100644 --- a/PWGJE/PWGJEEMCALJetTasksLinkDef.h +++ b/PWGJE/PWGJEEMCALJetTasksLinkDef.h @@ -58,6 +58,7 @@ #pragma link C++ class AliAnalysisTaskJetMassResponseDet+; #pragma link C++ class AliAnalysisTaskJetMatching+; #pragma link C++ class AliAnalysisTaskJetV2+; +#pragma link C++ class AliAnalysisTaskEmcalJetv2QA+; #pragma link C++ class AliAnalysisTaskRhoMass+; #pragma link C++ class AliAnalysisTaskRhoMassBase+; #pragma link C++ class AliAnalysisTaskSAJF+; -- 2.43.0