//
// Author: S.Aiola
-#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#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),
+ fParticleLevel(kFALSE),
+ fIsMC(kFALSE),
fNclusters(0),
fNtracks(0),
fNjets(0),
fHistJetsParts(0),
fHistCellsCent(0),
fHistCellsTracks(0),
+ fHistTrNegativeLabels(0),
+ fHistTrZeroLabels(0),
fHistTrEmcPhiEta(0),
fHistTrPhiEtaNonProp(0),
fHistDeltaEtaPt(0),
fHistDeltaPhiPt(0),
fHistNCellsEnergy(0),
+ fHistFcrossEnergy(0),
fHistClusTimeEnergy(0),
fHistCellsAbsIdEnergy(0),
fHistChVSneCells(0),
{
// Default constructor.
- for (Int_t i = 0; i < 5; i++) {
- fHistTrackPhi[i] = 0;
- fHistTrackEta[i] = 0;
- }
-
for (Int_t i = 0; i < 4; i++) {
- fHistTrPhiEtaPt[i] = 0;
+ for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+ fHistTrPhiEtaPtZeroLab[i] = 0;
fHistClusPhiEtaEnergy[i] = 0;
fHistJetsPhiEtaPt[i] = 0;
fHistJetsPtArea[i] = 0;
AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
AliAnalysisTaskEmcalJet(name, kTRUE),
fCellEnergyCut(0.1),
+ fParticleLevel(kFALSE),
+ fIsMC(kFALSE),
fNclusters(0),
fNtracks(0),
fNjets(0),
fHistJetsParts(0),
fHistCellsCent(0),
fHistCellsTracks(0),
+ fHistTrNegativeLabels(0),
+ fHistTrZeroLabels(0),
fHistTrEmcPhiEta(0),
fHistTrPhiEtaNonProp(0),
fHistDeltaEtaPt(0),
fHistDeltaPhiPt(0),
fHistNCellsEnergy(0),
+ fHistFcrossEnergy(0),
fHistClusTimeEnergy(0),
fHistCellsAbsIdEnergy(0),
fHistChVSneCells(0),
{
// Standard constructor.
- for (Int_t i = 0; i < 5; i++) {
- fHistTrackPhi[i] = 0;
- fHistTrackEta[i] = 0;
- }
-
for (Int_t i = 0; i < 4; i++) {
- fHistTrPhiEtaPt[i] = 0;
+ for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+ fHistTrPhiEtaPtZeroLab[i] = 0;
fHistClusPhiEtaEnergy[i] = 0;
fHistJetsPhiEtaPt[i] = 0;
fHistJetsPtArea[i] = 0;
AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
if (!fTracksName.IsNull()) {
+ if (!fParticleLevel && fIsMC) {
+ fHistTrNegativeLabels = new TH1F("fHistTrNegativeLabels","fHistTrNegativeLabels", 500, 0, 1);
+ fHistTrNegativeLabels->GetXaxis()->SetTitle("% of negative labels");
+ fHistTrNegativeLabels->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistTrNegativeLabels);
+
+ fHistTrZeroLabels = new TH1F("fHistTrZeroLabels","fHistTrZeroLabels", 500, 0, 1);
+ fHistTrZeroLabels->GetXaxis()->SetTitle("% of negative labels");
+ fHistTrZeroLabels->GetYaxis()->SetTitle("counts");
+ fOutput->Add(fHistTrZeroLabels);
+ }
+
fHistTracksCent = new TH2F("fHistTracksCent","Tracks vs. centrality", 100, 0, 100, fNbins, 0, 4000);
fHistTracksCent->GetXaxis()->SetTitle("Centrality (%)");
fHistTracksCent->GetYaxis()->SetTitle("No. of tracks");
TString histname;
- for (Int_t i = 0; i < 4; i++) {
- histname = "fHistTrPhiEtaPt_";
- histname += i;
- fHistTrPhiEtaPt[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
- fHistTrPhiEtaPt[i] ->GetXaxis()->SetTitle("#eta");
- fHistTrPhiEtaPt[i] ->GetYaxis()->SetTitle("#phi");
- fHistTrPhiEtaPt[i] ->GetZaxis()->SetTitle("p_{T} (GeV/c)");
- fOutput->Add(fHistTrPhiEtaPt[i]);
+ 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 && fIsMC) {
+ histname = Form("fHistTrPhiEtaPtZeroLab_%d",i);
+ fHistTrPhiEtaPtZeroLab[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+ fHistTrPhiEtaPtZeroLab[i]->GetXaxis()->SetTitle("#eta");
+ fHistTrPhiEtaPtZeroLab[i]->GetYaxis()->SetTitle("#phi");
+ fHistTrPhiEtaPtZeroLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+ fOutput->Add(fHistTrPhiEtaPtZeroLab[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);
-
- for (Int_t i = 0; i < 5; i++) {
- TString histnamephi("fHistTrackPhi_");
- histnamephi += i;
- fHistTrackPhi[i] = new TH1F(histnamephi.Data(),histnamephi.Data(), 201, 0, TMath::Pi() * 2.01);
- fHistTrackPhi[i]->GetXaxis()->SetTitle("Phi");
- fOutput->Add(fHistTrackPhi[i]);
+ if (!fParticleLevel) {
+ 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);
- TString histnameeta("fHistTrackEta_");
- histnameeta += i;
- fHistTrackEta[i] = new TH1F(histnameeta.Data(),histnameeta.Data(), 100, -1, 1);
- fHistTrackEta[i]->GetXaxis()->SetTitle("Eta");
- fOutput->Add(fHistTrackEta[i]);
+ 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);
}
-
- 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 (!fCaloName.IsNull()) {
TString histname;
- for (Int_t i = 0; i < 4; i++) {
+ 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("Energy (GeV)");
+ fHistClusPhiEtaEnergy[i]->GetZaxis()->SetTitle("E_{cluster} (GeV)");
fOutput->Add(fHistClusPhiEtaEnergy[i]);
}
fHistClusTimeEnergy = new TH2F("fHistClusTimeEnergy","Time vs. energy of clusters", fNbins, fMinBinPt, fMaxBinPt, fNbins, -1e-6, 1e-6);
- fHistClusTimeEnergy->GetXaxis()->SetTitle("Energy (GeV)");
+ 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("Energy (GeV)");
+ fHistNCellsEnergy->GetXaxis()->SetTitle("E_{cluster} (GeV)");
fHistNCellsEnergy->GetYaxis()->SetTitle("N_{cells}");
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("Energy (GeV)");
+ fHistCellsAbsIdEnergy->GetYaxis()->SetTitle("E_{cluster} (GeV)");
fHistCellsAbsIdEnergy->GetZaxis()->SetTitle("counts");
fOutput->Add(fHistCellsAbsIdEnergy);
TString histname;
- for (Int_t i = 0; i < 4; i++) {
+ for (Int_t i = 0; i < fNcentBins; i++) {
histname = "fHistJetsPhiEtaPt_";
histname += 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);
if (fTracks) {
trackSum = DoTrackLoop();
-
+ AliDebug(2,Form("%d tracks found in the event", fTracks->GetEntriesFast()));
fHistTracksCent->Fill(fCent, fNtracks);
}
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();
fHistClusTimeEnergy->Fill(cluster->E(), cluster->GetTOF());
+ if (cells)
+ fHistFcrossEnergy->Fill(cluster->E(), GetFcross(cluster, cells));
+
fNclusters++;
}
Float_t sum = 0;
Int_t ntracks = fTracks->GetEntriesFast();
- Int_t nclusters = 0;
- if (fCaloClusters)
- nclusters = fCaloClusters->GetEntriesFast();
+ Int_t neg = 0;
+ Int_t zero = 0;
for (Int_t i = 0; i < ntracks; i++) {
continue;
}
- AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
-
- if (vtrack && !AcceptTrack(vtrack, kTRUE))
+ if (!AcceptTrack(track))
continue;
fNtracks++;
sum += track->P();
-
- fHistTrPhiEtaPt[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
-
- Int_t label = track->GetLabel();
- if (label >= 0 && label < 4) {
- fHistTrackEta[label]->Fill(track->Eta());
- fHistTrackPhi[label]->Fill(track->Phi());
+ if (fParticleLevel) {
+ fHistTrPhiEtaPt[fCentBin][0]->Fill(track->Eta(), track->Phi(), track->Pt());
+ }
+ else {
+ fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
+ if (track->GetLabel() == 0) {
+ zero++;
+ if (fHistTrPhiEtaPtZeroLab[fCentBin])
+ fHistTrPhiEtaPtZeroLab[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
+ }
+
+ if (track->GetLabel() < 0)
+ neg++;
+
+ 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
+ AliDebug(2,Form("%s: track type %d not recognized!", GetName(), type));
}
- fHistTrackEta[4]->Fill(track->Eta());
- fHistTrackPhi[4]->Fill(track->Phi());
+ AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
if (!vtrack)
continue;
- if (vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999)
+ if ((vtrack->GetTrackEtaOnEMCal() == -999 || vtrack->GetTrackPhiOnEMCal() == -999) && fHistTrPhiEtaNonProp)
fHistTrPhiEtaNonProp->Fill(vtrack->Eta(), vtrack->Phi());
- fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
- fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
- fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
+ if (fHistTrEmcPhiEta)
+ fHistTrEmcPhiEta->Fill(vtrack->GetTrackEtaOnEMCal(), vtrack->GetTrackPhiOnEMCal());
+ if (fHistDeltaEtaPt)
+ fHistDeltaEtaPt->Fill(vtrack->Pt(), vtrack->Eta() - vtrack->GetTrackEtaOnEMCal());
+ if (fHistDeltaPhiPt)
+ fHistDeltaPhiPt->Fill(vtrack->Pt(), vtrack->Phi() - vtrack->GetTrackPhiOnEMCal());
}
-
+
+ if (fHistTrNegativeLabels)
+ fHistTrNegativeLabels->Fill(1. * neg / ntracks);
+
+ if (fHistTrZeroLabels)
+ fHistTrZeroLabels->Fill(1. * zero / ntracks);
+
return sum;
}