#include "AliExternalTrackParam.h"
#include "AliTrackerBase.h"
#include "AliLog.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALGeoParams.h"
+#include "AliPicoTrack.h"
#include "AliAnalysisTaskSAQA.h"
AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
AliAnalysisTaskEmcalJet("AliAnalysisTaskSAQA", kTRUE),
fCellEnergyCut(0.1),
- fDoTrigger(kFALSE),
- fRepropagateTracks(kFALSE),
- fTrgClusName("ClustersL1GAMMAFEE"),
- fTrgClusters(0),
- fHistCentrality(0),
- fHistZVertex(0),
+ fParticleLevel(kFALSE),
+ fNclusters(0),
+ fNtracks(0),
+ fNjets(0),
fHistTracksCent(0),
fHistClusCent(0),
+ fHistJetsCent(0),
fHistClusTracks(0),
+ fHistJetsParts(0),
fHistCellsCent(0),
fHistCellsTracks(0),
- fHistMaxL1FastORCent(0),
- fHistMaxL1ClusCent(0),
- fHistMaxL1ThrCent(0),
- fHistTracksPt(0),
- fHistTrPhiEta(0),
fHistTrEmcPhiEta(0),
fHistTrPhiEtaNonProp(0),
fHistDeltaEtaPt(0),
fHistDeltaPhiPt(0),
- fHistDeltaEtaNewProp(0),
- fHistDeltaPhiNewProp(0),
- fHistClusPhiEtaEnergy(0),
fHistNCellsEnergy(0),
+ fHistFcrossEnergy(0),
fHistClusTimeEnergy(0),
- fHistCellsEnergy(0),
+ fHistCellsAbsIdEnergy(0),
fHistChVSneCells(0),
fHistChVSneClus(0),
fHistChVSneCorrCells(0)
{
// Default constructor.
- for (Int_t i = 0; i < 5; i++) {
- fHistTrackPhi[i] = 0;
- fHistTrackEta[i] = 0;
- }
-
- for (Int_t i = 0; i < 6; i++) {
- fHistTrackPhiPt[i] = 0;
- fHistTrackEtaPt[i] = 0;
- }
-
for (Int_t i = 0; i < 4; i++) {
- fHistJetsPhiEta[i] = 0;
- fHistJetsPtNonBias[i] = 0;
- fHistJetsPtTrack[i] = 0;
- fHistJetsPtClus[i] = 0;
- fHistJetsPt[i] = 0;
+ for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+ fHistTrPhiEtaPtNegLab[i] = 0;
+ fHistClusPhiEtaEnergy[i] = 0;
+ fHistJetsPhiEtaPt[i] = 0;
fHistJetsPtArea[i] = 0;
- fHistJetsPtAreaNonBias[i] = 0;
}
+
+ SetMakeGeneralHistograms(kTRUE);
}
//________________________________________________________________________
AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
AliAnalysisTaskEmcalJet(name, kTRUE),
fCellEnergyCut(0.1),
- fDoTrigger(kFALSE),
- fRepropagateTracks(kFALSE),
- fTrgClusName("ClustersL1GAMMAFEE"),
- fTrgClusters(0),
- fHistCentrality(0),
- fHistZVertex(0),
+ fParticleLevel(kFALSE),
+ fNclusters(0),
+ fNtracks(0),
+ fNjets(0),
fHistTracksCent(0),
fHistClusCent(0),
+ fHistJetsCent(0),
fHistClusTracks(0),
+ fHistJetsParts(0),
fHistCellsCent(0),
fHistCellsTracks(0),
- fHistMaxL1FastORCent(0),
- fHistMaxL1ClusCent(0),
- fHistMaxL1ThrCent(0),
- fHistTracksPt(0),
- fHistTrPhiEta(0),
fHistTrEmcPhiEta(0),
fHistTrPhiEtaNonProp(0),
fHistDeltaEtaPt(0),
fHistDeltaPhiPt(0),
- fHistDeltaEtaNewProp(0),
- fHistDeltaPhiNewProp(0),
- fHistClusPhiEtaEnergy(0),
fHistNCellsEnergy(0),
+ fHistFcrossEnergy(0),
fHistClusTimeEnergy(0),
- fHistCellsEnergy(0),
+ fHistCellsAbsIdEnergy(0),
fHistChVSneCells(0),
fHistChVSneClus(0),
fHistChVSneCorrCells(0)
{
// Standard constructor.
- for (Int_t i = 0; i < 5; i++) {
- fHistTrackPhi[i] = 0;
- fHistTrackEta[i] = 0;
- }
-
- for (Int_t i = 0; i < 6; i++) {
- fHistTrackPhiPt[i] = 0;
- fHistTrackEtaPt[i] = 0;
- }
-
for (Int_t i = 0; i < 4; i++) {
- fHistJetsPhiEta[i] = 0;
- fHistJetsPtNonBias[i] = 0;
- fHistJetsPtTrack[i] = 0;
- fHistJetsPtClus[i] = 0;
- fHistJetsPt[i] = 0;
+ for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+ fHistTrPhiEtaPtNegLab[i] = 0;
+ fHistClusPhiEtaEnergy[i] = 0;
+ fHistJetsPhiEtaPt[i] = 0;
fHistJetsPtArea[i] = 0;
- fHistJetsPtAreaNonBias[i] = 0;
}
+
+ SetMakeGeneralHistograms(kTRUE);
}
//________________________________________________________________________
{
// Create histograms
- OpenFile(1);
- fOutput = new TList();
- fOutput->SetOwner();
+ AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
- fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
- fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
- fHistCentrality->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistCentrality);
+ if (!fTracksName.IsNull()) {
+ fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
+ fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
+ fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
+ fOutput->Add(fHistTracksCent);
- fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
- fHistZVertex->GetXaxis()->SetTitle("z");
- fHistZVertex->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistZVertex);
+ TString histname;
- fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
- fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
- fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
- fOutput->Add(fHistTracksCent);
+ Int_t nlabels = 4;
+ if (fParticleLevel)
+ nlabels = 1;
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ for (Int_t j = 0; j < nlabels; j++) {
+ histname = Form("fHistTrPhiEtaPt_%d_%d",i,j);
+ fHistTrPhiEtaPt[i][j] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+ fHistTrPhiEtaPt[i][j]->GetXaxis()->SetTitle("#eta");
+ fHistTrPhiEtaPt[i][j]->GetYaxis()->SetTitle("#phi");
+ fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+ fOutput->Add(fHistTrPhiEtaPt[i][j]);
+ }
+ if (!fParticleLevel) {
+ histname = Form("fHistTrPhiEtaPtNegLab_%d",i);
+ fHistTrPhiEtaPtNegLab[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+ fHistTrPhiEtaPtNegLab[i]->GetXaxis()->SetTitle("#eta");
+ fHistTrPhiEtaPtNegLab[i]->GetYaxis()->SetTitle("#phi");
+ fHistTrPhiEtaPtNegLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+ fOutput->Add(fHistTrPhiEtaPtNegLab[i]);
+ }
+ }
+
+ fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+ fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
+ fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
+ fOutput->Add(fHistTrEmcPhiEta);
+
+ fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
+ fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
+ fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
+ fOutput->Add(fHistTrPhiEtaNonProp);
+
+ fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
+ fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
+ fOutput->Add(fHistDeltaEtaPt);
+
+ fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
+ fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
+ fOutput->Add(fHistDeltaPhiPt);
+ }
- if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
+ if (!fCaloName.IsNull()) {
fHistClusCent = new TH2F("fHistClusCent","Clusters vs. centrality", 100, 0, 100, fNbins, 0, 2000);
fHistClusCent->GetXaxis()->SetTitle("Centrality (%)");
fHistClusCent->GetYaxis()->SetTitle("No. of clusters");
fHistCellsTracks->GetYaxis()->SetTitle("No. of EMCal cells");
fOutput->Add(fHistCellsTracks);
- fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins, -1e-6, 1e-6);
- fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
- fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
- fOutput->Add(fHistClusTimeEnergy);
+ TString histname;
- if (fDoTrigger) {
- fHistMaxL1FastORCent = new TH2F("fHistMaxL1FastORCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
- fHistMaxL1FastORCent->GetXaxis()->SetTitle("Centrality [%]");
- fHistMaxL1FastORCent->GetYaxis()->SetTitle("Maximum L1 FastOR");
- fOutput->Add(fHistMaxL1FastORCent);
-
- fHistMaxL1ClusCent = new TH2F("fHistMaxL1ClusCent","fHistMaxL1ClusCent", 100, 0, 100, 250, 0, 250);
- fHistMaxL1ClusCent->GetXaxis()->SetTitle("Centrality [%]");
- fHistMaxL1ClusCent->GetYaxis()->SetTitle("Maximum L1 trigger cluster");
- fOutput->Add(fHistMaxL1ClusCent);
-
- fHistMaxL1ThrCent = new TH2F("fHistMaxL1ThrCent","fHistMaxL1ThrCent", 100, 0, 100, 250, 0, 250);
- fHistMaxL1ThrCent->GetXaxis()->SetTitle("Centrality [%]");
- fHistMaxL1ThrCent->GetYaxis()->SetTitle("Maximum L1 threshold");
- fOutput->Add(fHistMaxL1ThrCent);
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ histname = "fHistClusPhiEtaEnergy_";
+ histname += i;
+ fHistClusPhiEtaEnergy[i] = new TH3F(histname, histname, 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+ fHistClusPhiEtaEnergy[i]->GetXaxis()->SetTitle("#eta");
+ fHistClusPhiEtaEnergy[i]->GetYaxis()->SetTitle("#phi");
+ fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
+ fOutput->Add(fHistClusPhiEtaEnergy[i]);
}
- }
-
- fHistTracksPt = new TH1F("fHistTracksPt","p_{T} spectrum of reconstructed tracks", fNbins, fMinBinPt, fMaxBinPt);
- fHistTracksPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistTracksPt->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistTracksPt);
-
- fHistTrPhiEta = new TH2F("fHistTrPhiEta","Phi-Eta distribution of tracks", 80, -2, 2, 128, 0, 6.4);
- fHistTrPhiEta->GetXaxis()->SetTitle("#eta");
- fHistTrPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistTrPhiEta);
-
- fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 80, -2, 2, 128, 0, 6.4);
- fHistTrEmcPhiEta->GetXaxis()->SetTitle("#eta");
- fHistTrEmcPhiEta->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistTrEmcPhiEta);
-
- fHistTrPhiEtaNonProp = new TH2F("fHistTrPhiEtaNonProp","fHistTrPhiEtaNonProp", 80, -2, 2, 128, 0, 6.4);
- fHistTrPhiEtaNonProp->GetXaxis()->SetTitle("#eta");
- fHistTrPhiEtaNonProp->GetYaxis()->SetTitle("#phi");
- fOutput->Add(fHistTrPhiEtaNonProp);
-
- fHistDeltaEtaPt = new TH2F("fHistDeltaEtaPt","fHistDeltaEtaPt", fNbins, fMinBinPt, fMaxBinPt, 80, -0.5, 0.5);
- fHistDeltaEtaPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistDeltaEtaPt->GetYaxis()->SetTitle("#delta#eta");
- fOutput->Add(fHistDeltaEtaPt);
-
- fHistDeltaPhiPt = new TH2F("fHistDeltaPhiPt","fHistDeltaPhiPt", fNbins, fMinBinPt, fMaxBinPt, 256, -1.6, 4.8);
- fHistDeltaPhiPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistDeltaPhiPt->GetYaxis()->SetTitle("#delta#phi");
- fOutput->Add(fHistDeltaPhiPt);
-
- if (fRepropagateTracks) {
- fHistDeltaEtaNewProp = new TH1F("fHistDeltaEtaNewProp","fHistDeltaEtaNewProp", 800, -0.4, 0.4);
- fHistDeltaEtaNewProp->GetXaxis()->SetTitle("#delta#eta");
- fHistDeltaEtaNewProp->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaEtaNewProp);
-
- fHistDeltaPhiNewProp = new TH1F("fHistDeltaPhiNewProp","fHistDeltaPhiNewProp", 2560, -1.6, 3.2);
- fHistDeltaPhiNewProp->GetXaxis()->SetTitle("#delta#phi");
- fHistDeltaPhiNewProp->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistDeltaPhiNewProp);
- }
- if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
- fHistClusPhiEtaEnergy = new TH3F("fHistClusPhiEtaEnergy","Phi-Eta-Energy distribution of clusters", fNbins, fMinBinPt, fMaxBinPt, 80, -2, 2, 128, 0, 6.4);
- fHistClusPhiEtaEnergy->GetXaxis()->SetTitle("E [GeV]");
- fHistClusPhiEtaEnergy->GetYaxis()->SetTitle("#eta");
- fHistClusPhiEtaEnergy->GetZaxis()->SetTitle("#phi");
- fOutput->Add(fHistClusPhiEtaEnergy);
+ fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins, -1e-6, 1e-6);
+ fHistClusTimeEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
+ fHistClusTimeEnergy->GetYaxis()->SetTitle("Time");
+ fOutput->Add(fHistClusTimeEnergy);
fHistNCellsEnergy = new TH2F("fHistNCellsEnergy","Number of cells vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, 30, 0, 30);
- fHistNCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
+ fHistNCellsEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
- fOutput->Add(fHistNCellsEnergy);
- }
-
- if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
-
- fHistCellsEnergy = new TH1F("fHistCellsEnergy","Energy spectrum of cells", fNbins, fMinBinPt, fMaxBinPt);
- fHistCellsEnergy->GetXaxis()->SetTitle("E [GeV]");
- fHistCellsEnergy->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistCellsEnergy);
+ fOutput->Add(fHistNCellsEnergy);
+
+ fHistFcrossEnergy = new TH2F("fHistFcrossEnergy","fHistFcrossEnergy", fNbins, fMinBinPt, fMaxBinPt, 200, -3.5, 1.5);
+ fHistFcrossEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
+ fHistFcrossEnergy->GetYaxis()->SetTitle("F_{cross}");
+ fOutput->Add(fHistFcrossEnergy);
+
+ fHistCellsAbsIdEnergy = new TH2F("fHistCellsAbsIdEnergy","fHistCellsAbsIdEnergy", 11600,0,11599,(Int_t)(fNbins / 2), fMinBinPt, fMaxBinPt / 2);
+ fHistCellsAbsIdEnergy->GetXaxis()->SetTitle("cell abs. Id");
+ fHistCellsAbsIdEnergy->GetYaxis()->SetTitle("E_{cluster} (GeV)");
+ fHistCellsAbsIdEnergy->GetZaxis()->SetTitle("counts");
+ fOutput->Add(fHistCellsAbsIdEnergy);
fHistChVSneCells = new TH2F("fHistChVSneCells","Charged energy vs. neutral (cells) energy",
(Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
- fHistChVSneCells->GetXaxis()->SetTitle("E [GeV]");
- fHistChVSneCells->GetYaxis()->SetTitle("P [GeV/c]");
+ fHistChVSneCells->GetXaxis()->SetTitle("Energy (GeV)");
+ fHistChVSneCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
fOutput->Add(fHistChVSneCells);
fHistChVSneClus = new TH2F("fHistChVSneClus","Charged energy vs. neutral (clusters) energy",
(Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
- fHistChVSneClus->GetXaxis()->SetTitle("E [GeV]");
- fHistChVSneClus->GetYaxis()->SetTitle("P [GeV/c]");
+ fHistChVSneClus->GetXaxis()->SetTitle("Energy (GeV)");
+ fHistChVSneClus->GetYaxis()->SetTitle("Momentum (GeV/c)");
fOutput->Add(fHistChVSneClus);
fHistChVSneCorrCells = new TH2F("fHistChVSneCorrCells","Charged energy vs. neutral (corrected cells) energy",
(Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, (Int_t)(fNbins * 2.5), fMinBinPt , fMaxBinPt * 2.5);
- fHistChVSneCorrCells->GetXaxis()->SetTitle("E [GeV]");
- fHistChVSneCorrCells->GetYaxis()->SetTitle("P [GeV/c]");
+ fHistChVSneCorrCells->GetXaxis()->SetTitle("Energy (GeV)");
+ fHistChVSneCorrCells->GetYaxis()->SetTitle("Momentum (GeV/c)");
fOutput->Add(fHistChVSneCorrCells);
}
-
- for (Int_t i = 0; i < 5; i++) {
- TString histnamephi("fHistTrackPhi_");
- histnamephi += i;
- fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 128, 0, 6.4);
- fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
- fOutput->Add(fHistTrackPhi[i]);
-
- TString histnameeta("fHistTrackEta_");
- histnameeta += i;
- fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -2, 2);
- fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
- fOutput->Add(fHistTrackEta[i]);
- }
-
- for (Int_t i = 0; i < 6; i++) {
- TString histnamephipt("fHistTrackPhiPt_");
- histnamephipt += i;
- fHistTrackPhiPt[i] = new TH1F(histnamephipt.Data(),histnamephipt.Data(), 128, 0, 6.4);
- fHistTrackPhiPt[i]->GetXaxis()->SetTitle("Phi");
- fOutput->Add(fHistTrackPhiPt[i]);
-
- TString histnameetapt("fHistTrackEtaPt_");
- histnameetapt += i;
- fHistTrackEtaPt[i] = new TH1F(histnameetapt.Data(),histnameetapt.Data(), 100, -2, 2);
- fHistTrackEtaPt[i]->GetXaxis()->SetTitle("Eta");
- fOutput->Add(fHistTrackEtaPt[i]);
- }
-
- fHistTrackPhi[0]->SetLineColor(kRed);
- fHistTrackEta[0]->SetLineColor(kRed);
- fHistTrackPhi[1]->SetLineColor(kBlue);
- fHistTrackEta[1]->SetLineColor(kBlue);
- fHistTrackPhi[2]->SetLineColor(kGreen);
- fHistTrackEta[2]->SetLineColor(kGreen);
- fHistTrackPhi[3]->SetLineColor(kOrange);
- fHistTrackEta[3]->SetLineColor(kOrange);
- fHistTrackPhi[4]->SetLineColor(kBlack);
- fHistTrackEta[4]->SetLineColor(kBlack);
-
+
if (!fJetsName.IsNull()) {
+ fHistJetsCent = new TH2F("fHistJetsCent","Jets vs. centrality", 100, 0, 100, 150, -0.5, 149.5);
+ fHistJetsCent->GetXaxis()->SetTitle("Centrality (%)");
+ fHistJetsCent->GetYaxis()->SetTitle("No. of jets");
+ fOutput->Add(fHistJetsCent);
+
+ fHistJetsParts = new TH2F("fHistJetsParts","Jets vs. centrality", fNbins, 0, 6000, 150, -0.5, 149.5);
+ fHistJetsParts->GetXaxis()->SetTitle("No. of particles");
+ fHistJetsParts->GetYaxis()->SetTitle("No. of jets");
+ fOutput->Add(fHistJetsParts);
TString histname;
- for (Int_t i = 0; i < 4; i++) {
- histname = "fHistJetsPhiEta_";
- histname += i;
- fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 80, -2, 2, 128, 0, 6.4);
- fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
- fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
- fHistJetsPhiEta[i]->GetZaxis()->SetTitle("p_{T} [GeV/c]");
- fOutput->Add(fHistJetsPhiEta[i]);
-
- histname = "fHistJetsPtNonBias_";
- histname += i;
- fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
- fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistJetsPtNonBias[i]);
-
- histname = "fHistJetsPtTrack_";
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ histname = "fHistJetsPhiEtaPt_";
histname += i;
- fHistJetsPtTrack[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
- fHistJetsPtTrack[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistJetsPtTrack[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistJetsPtTrack[i]);
-
- if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
- histname = "fHistJetsPtClus_";
- histname += i;
- fHistJetsPtClus[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
- fHistJetsPtClus[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistJetsPtClus[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistJetsPtClus[i]);
- }
-
- histname = "fHistJetsPt_";
- histname += i;
- fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
- fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
- fOutput->Add(fHistJetsPt[i]);
-
- histname = "fHistJetsPtAreaNonBias_";
- histname += i;
- fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
- fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
- fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
- fOutput->Add(fHistJetsPtAreaNonBias[i]);
+ fHistJetsPhiEtaPt[i] = new TH3F(histname.Data(), histname.Data(), 100, -1.2, 1.2, 201, 0, TMath::Pi() * 2.01, (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5);
+ fHistJetsPhiEtaPt[i]->GetXaxis()->SetTitle("#eta");
+ fHistJetsPhiEtaPt[i]->GetYaxis()->SetTitle("#phi");
+ fHistJetsPhiEtaPt[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+ fOutput->Add(fHistJetsPhiEtaPt[i]);
histname = "fHistJetsPtArea_";
histname += i;
- fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
- fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+ fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), (Int_t)(fNbins * 2.5), fMinBinPt, fMaxBinPt * 2.5, 30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
+ fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
fOutput->Add(fHistJetsPtArea[i]);
}
if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
return kFALSE;
- if (!fTrgClusName.IsNull() && fDoTrigger && !fTrgClusters) {
- fTrgClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
- if (!fTrgClusters) {
- AliError(Form("%s: Could not retrieve trigger clusters %s!", GetName(), fTrgClusName.Data()));
- return kFALSE;
- }
- else {
- TClass *cl = fTrgClusters->GetClass();
- if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
- AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fTrgClusName.Data()));
- fTrgClusters = 0;
- return kFALSE;
- }
- }
- }
+ fNclusters = 0;
+ fNtracks = 0;
+ fNjets = 0;
return kTRUE;
}
+
//________________________________________________________________________
Bool_t AliAnalysisTaskSAQA::FillHistograms()
{
// Fill histograms.
- fHistCentrality->Fill(fCent);
- if (fTracks)
- fHistTracksCent->Fill(fCent, fTracks->GetEntriesFast());
- if (fCaloClusters)
- fHistClusCent->Fill(fCent, fCaloClusters->GetEntriesFast());
-
- fHistZVertex->Fill(fVertex[2]);
-
- Float_t trackSum = DoTrackLoop();
+ Float_t clusSum = 0;
+ Float_t trackSum = 0;
- DoJetLoop();
+ if (fTracks) {
+ trackSum = DoTrackLoop();
- if (fAnaType == kEMCAL || fAnaType == kEMCALOnly) {
+ fHistTracksCent->Fill(fCent, fNtracks);
+ }
- if (fTracks && fCaloClusters)
- fHistClusTracks->Fill(fTracks->GetEntriesFast(), fCaloClusters->GetEntriesFast());
+ if (fCaloClusters) {
+ clusSum = DoClusterLoop();
- Float_t clusSum = DoClusterLoop();
+ fHistClusCent->Fill(fCent, fNclusters);
+ fHistClusTracks->Fill(fNtracks, fNclusters);
Float_t cellSum = 0, cellCutSum = 0;
Int_t ncells = DoCellLoop(cellSum, cellCutSum);
if (fTracks)
- fHistCellsTracks->Fill(fTracks->GetEntriesFast(), ncells);
+ fHistCellsTracks->Fill(fNtracks, ncells);
fHistCellsCent->Fill(fCent, ncells);
fHistChVSneCells->Fill(cellSum, trackSum);
fHistChVSneClus->Fill(clusSum, trackSum);
fHistChVSneCorrCells->Fill(cellCutSum, trackSum);
+ }
- if (fDoTrigger) {
- Float_t maxTrgClus = DoTriggerClusLoop();
- fHistMaxL1ClusCent->Fill(fCent, maxTrgClus);
+ if (fJets) {
+ DoJetLoop();
- Int_t maxL1amp = -1;
- Int_t maxL1thr = -1;
-
- DoTriggerPrimitives(maxL1amp, maxL1thr);
-
- if (maxL1amp > -1)
- fHistMaxL1FastORCent->Fill(fCent, maxL1amp);
-
- if (maxL1thr > -1)
- fHistMaxL1ThrCent->Fill(fCent, maxL1thr);
- }
+ fHistJetsCent->Fill(fCent, fNjets);
+ fHistJetsParts->Fill(fNtracks + fNclusters, fNjets);
}
return kTRUE;
const Int_t ncells = cells->GetNumberOfCells();
for (Int_t pos = 0; pos < ncells; pos++) {
- Float_t amp = cells->GetAmplitude(pos);
- fHistCellsEnergy->Fill(amp);
+ Float_t amp = cells->GetAmplitude(pos);
+ Int_t absId = cells->GetCellNumber(pos);
+ fHistCellsAbsIdEnergy->Fill(absId,amp);
sum += amp;
if (amp < fCellEnergyCut)
continue;
return ncells;
}
+//________________________________________________________________________
+Double_t AliAnalysisTaskSAQA::GetFcross(AliVCluster *cluster, AliVCaloCells *cells)
+{
+ Int_t AbsIdseed = -1;
+ Double_t Eseed = 0;
+ for (Int_t i = 0; i < cluster->GetNCells(); i++) {
+ if (cells->GetCellAmplitude(cluster->GetCellAbsId(i)) > AbsIdseed) {
+ Eseed = cells->GetCellAmplitude(cluster->GetCellAbsId(i));
+ AbsIdseed = cluster->GetCellAbsId(i);
+ }
+ }
+
+ if (Eseed < 1e-9)
+ return 100;
+
+ Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
+ fGeom->GetCellIndex(AbsIdseed,imod,iTower,iIphi,iIeta);
+ fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi,iIeta,iphi,ieta);
+
+ //Get close cells index and energy, not in corners
+
+ Int_t absID1 = -1;
+ Int_t absID2 = -1;
+
+ if (iphi < AliEMCALGeoParams::fgkEMCALRows-1) absID1 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
+ if (iphi > 0) absID2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
+
+ // In case of cell in eta = 0 border, depending on SM shift the cross cell index
+
+ Int_t absID3 = -1;
+ Int_t absID4 = -1;
+
+ if (ieta == AliEMCALGeoParams::fgkEMCALCols-1 && !(imod%2)) {
+ absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod+1, iphi, 0);
+ absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
+ }
+ else if (ieta == 0 && imod%2) {
+ absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
+ absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod-1, iphi, AliEMCALGeoParams::fgkEMCALCols-1);
+ }
+ else {
+ if (ieta < AliEMCALGeoParams::fgkEMCALCols-1)
+ absID3 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1);
+ if (ieta > 0)
+ absID4 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1);
+ }
+
+ Double_t ecell1 = cells->GetCellAmplitude(absID1);
+ Double_t ecell2 = cells->GetCellAmplitude(absID2);
+ Double_t ecell3 = cells->GetCellAmplitude(absID3);
+ Double_t ecell4 = cells->GetCellAmplitude(absID4);
+
+ Double_t Ecross = ecell1 + ecell2 + ecell3 + ecell4;
+
+ Double_t Fcross = 1 - Ecross/Eseed;
+
+ return Fcross;
+}
+
//________________________________________________________________________
Float_t AliAnalysisTaskSAQA::DoClusterLoop()
{
if (!fCaloClusters)
return 0;
+ AliVCaloCells *cells = InputEvent()->GetEMCALCells();
+
Float_t sum = 0;
// Cluster loop
continue;
}
- if (!AcceptCluster(cluster, kTRUE))
+ if (!AcceptCluster(cluster))
continue;
sum += cluster->E();
TLorentzVector nPart;
cluster->GetMomentum(nPart, fVertex);
- fHistClusPhiEtaEnergy->Fill(cluster->E(), nPart.Eta(), nPart.Phi());
+ fHistClusPhiEtaEnergy[fCentBin]->Fill(nPart.Eta(), nPart.Phi(), cluster->E());
fHistNCellsEnergy->Fill(cluster->E(), cluster->GetNCells());
fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
+
+ if (cells)
+ fHistFcrossEnergy->Fill(cluster->E(), GetFcross(cluster, cells));
+
+ fNclusters++;
}
return sum;
AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
- if (vtrack && !AcceptTrack(vtrack, kTRUE))
+ if (vtrack && !AcceptTrack(vtrack))
continue;
-
- fHistTracksPt->Fill(track->Pt());
-
- sum += track->P();
- Int_t label = track->GetLabel();
-
- fHistTrPhiEta->Fill(track->Eta(), track->Phi());
-
- fHistTrackEta[4]->Fill(track->Eta());
- fHistTrackPhi[4]->Fill(track->Phi());
+ fNtracks++;
- if (label >= 0 && label < 4) {
- fHistTrackEta[label]->Fill(track->Eta());
- fHistTrackPhi[label]->Fill(track->Phi());
- }
+ sum += track->P();
- if (track->Pt() < 0.5) {
- fHistTrackPhiPt[0]->Fill(track->Phi());
- fHistTrackEtaPt[0]->Fill(track->Eta());
- }
- else if (track->Pt() < 1) {
- fHistTrackPhiPt[1]->Fill(track->Phi());
- fHistTrackEtaPt[1]->Fill(track->Eta());
- }
- else if (track->Pt() < 2) {
- fHistTrackPhiPt[2]->Fill(track->Phi());
- fHistTrackEtaPt[2]->Fill(track->Eta());
- }
- else if (track->Pt() < 3) {
- fHistTrackPhiPt[3]->Fill(track->Phi());
- fHistTrackEtaPt[3]->Fill(track->Eta());
- }
- else if (track->Pt() < 5) {
- fHistTrackPhiPt[4]->Fill(track->Phi());
- fHistTrackEtaPt[4]->Fill(track->Eta());
+ if (fParticleLevel) {
+ fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
}
else {
- fHistTrackPhiPt[5]->Fill(track->Phi());
- fHistTrackEtaPt[5]->Fill(track->Eta());
+ fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
+ if (fHistTrPhiEtaPtNegLab[fCentBin] && track->GetLabel() < 0)
+ fHistTrPhiEtaPtNegLab[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
+
+ Int_t type = 0;
+
+ AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
+ if (ptrack)
+ type = ptrack->GetTrackType();
+
+ if (type >= 0 && type < 3)
+ fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
+ else
+ AliWarning(Form("%s: track type %d not recognized!", GetName(), type));
}
if (!vtrack)
fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
-
- if (fRepropagateTracks && vtrack->GetTrackEtaOnEMCal() > -2) {
- Float_t propeta = -999, propphi = -999;
- PropagateTrack(vtrack, propeta, propphi);
- fHistDeltaEtaNewProp->Fill(propeta - vtrack->GetTrackEtaOnEMCal());
- fHistDeltaPhiNewProp->Fill(propphi - vtrack->GetTrackPhiOnEMCal());
- }
}
return sum;
}
-//____________________________________________________________________________
-void AliAnalysisTaskSAQA::PropagateTrack(AliVTrack *track, Float_t &eta, Float_t &phi)
-{
- eta = -999;
- phi = -999;
-
- if (!track)
- return;
-
- if (track->Pt() == 0)
- return;
-
- // init the magnetic field if not already on
- if(!TGeoGlobalMagField::Instance()->GetField()) {
- AliInfo("Init the magnetic field\n");
- AliAODEvent* aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
- if (aodevent) {
- Double_t curSol = 30000*aodevent->GetMagneticField()/5.00668;
- Double_t curDip = 6000 *aodevent->GetMuonMagFieldScale();
- AliMagF *field = AliMagF::CreateFieldMap(curSol,curDip);
- TGeoGlobalMagField::Instance()->SetField(field);
- }
- }
-
- Double_t cv[21];
- for (Int_t i = 0; i < 21; i++) cv[i] = 0;
-
- Double_t pos[3], mom[3];
- track->GetXYZ(pos);
- track->GetPxPyPz(mom);
- AliExternalTrackParam *trackParam = new AliExternalTrackParam(pos, mom, cv, track->Charge());
-
- if(!AliTrackerBase::PropagateTrackToBxByBz(trackParam, 430., 0.139, 20, kTRUE, 0.8, -1)) return;
- Double_t trkPos[3] = {0., 0., 0.};
- if(!trackParam->GetXYZ(trkPos)) return;
- TVector3 trkPosVec(trkPos[0], trkPos[1], trkPos[2]);
- eta = trkPosVec.Eta();
- phi = trkPosVec.Phi();
- if(phi < 0)
- phi += 2 * TMath::Pi();
-}
-
//________________________________________________________________________
void AliAnalysisTaskSAQA::DoJetLoop()
{
continue;
}
- if (!AcceptJet(jet, kFALSE))
+ if (!AcceptJet(jet))
continue;
- fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
- fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());
-
- if (jet->MaxTrackPt() > fPtBiasJetTrack)
- fHistJetsPtTrack[fCentBin]->Fill(jet->Pt());
-
- if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
- fHistJetsPtClus[fCentBin]->Fill(jet->Pt());
-
- if (!AcceptBiasJet(jet))
- continue;
+ fNjets++;
- fHistJetsPt[fCentBin]->Fill(jet->Pt());
+ fHistJetsPhiEtaPt[fCentBin]->Fill(jet->Eta(), jet->Phi(), jet->Pt());
fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
-
- fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
}
}
-//________________________________________________________________________
-Float_t AliAnalysisTaskSAQA::DoTriggerClusLoop()
-{
- // Do trigger cluster loop.
-
- if (!fTrgClusters)
- return 0;
-
- Int_t ntrgclusters = fTrgClusters->GetEntriesFast();
- Float_t maxTrgClus = 0;
-
- for (Int_t iClusters = 0; iClusters < ntrgclusters; iClusters++) {
- AliVCluster* cluster = static_cast<AliVCluster*>(fTrgClusters->At(iClusters));
- if (!cluster) {
- AliError(Form("Could not receive cluster %d", iClusters));
- continue;
- }
-
- if (!(cluster->IsEMCAL())) continue;
-
- if (cluster->E() > maxTrgClus)
- maxTrgClus = cluster->E();
-
- }
- return maxTrgClus;
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSAQA::DoTriggerPrimitives(Int_t &maxL1amp, Int_t &maxL1thr)
-{
- // Do trigger primitives loop.
-
- AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
-
- if (!triggers || triggers->GetEntries() == 0)
- return;
-
- triggers->Reset();
- Int_t L1amp = 0;
- Int_t L1thr = 0;
- maxL1amp = -1;
- maxL1thr = -1;
-
- while (triggers->Next()) {
- triggers->GetL1TimeSum(L1amp);
- if (maxL1amp < L1amp)
- maxL1amp = L1amp;
-
- triggers->GetL1Threshold(L1thr);
- if (maxL1thr < L1thr)
- maxL1thr = L1thr;
- }
-}
//________________________________________________________________________
void AliAnalysisTaskSAQA::Terminate(Option_t *)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff