DrawdNdeta.C can now draw centrality specific results and compare to other data (MC truth, other results, etc.)
Reduced the total number of events to normalise to to a simpler expresion.
*
*/
AliAnalysisTask*
-AddTaskCentraldNdeta(const char* trig="INEL",
- Double_t vzMin=-10,
- Double_t vzMax=10,
- Float_t centLow=0,
- Float_t centHigh=100)
+AddTaskCentraldNdeta(const char* trig = "INEL",
+ Double_t vzMin = -10,
+ Double_t vzMax = +10,
+ Bool_t useCent = false,
+ Bool_t cutEdges = false)
{
// analysis manager
AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
AliCentraldNdetaTask* task = new AliCentraldNdetaTask("Central");
task->SetVertexRange(vzMin, vzMax);
task->SetTriggerMask(trig);
- // task->SetUseShapeCorrection(false);
- task->AddCentralityBin( 0, 0); // All bin - integrate over centrality
- task->AddCentralityBin( 0, 5);
- task->AddCentralityBin( 5, 10);
- task->AddCentralityBin(10, 20);
- task->AddCentralityBin(20, 30);
- task->AddCentralityBin(30, 40);
- task->AddCentralityBin(40, 50);
- task->AddCentralityBin(50, 60);
- task->AddCentralityBin(60,100);
+ task->SetCutEdges(cutEdges);
+ task->SetUseShapeCorrection(false);
+ if (useCent) {
+ Short_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+ task->SetCentralityAxis(11, bins);
+ }
mgr->AddTask(task);
// create containers for input/output
*
*/
AliAnalysisTask*
-AddTaskForwarddNdeta(const char* trig="INEL",
- Double_t vzMin=-10,
- Double_t vzMax=10,
- Float_t centlow = 0,
- Float_t centhigh = 100,
- Bool_t cutEdges = false)
+AddTaskForwarddNdeta(const char* trig = "INEL",
+ Double_t vzMin = -10,
+ Double_t vzMax = +10,
+ Bool_t useCent = false,
+ Bool_t cutEdges = false)
{
// analysis manager
AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
task->SetVertexRange(vzMin, vzMax);
task->SetTriggerMask(trig);
task->SetCutEdges(cutEdges);
- // task->SetUseShapeCorrection(false);
- task->AddCentralityBin( 0, 0); // All bin - integrate over centrality
- task->AddCentralityBin( 0, 5);
- task->AddCentralityBin( 5, 10);
- task->AddCentralityBin(10, 20);
- task->AddCentralityBin(20, 30);
- task->AddCentralityBin(30, 40);
- task->AddCentralityBin(40, 50);
- task->AddCentralityBin(50, 60);
- task->AddCentralityBin(60,100);
+ task->SetUseShapeCorrection(false);
+ if (useCent) {
+ Short_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+ task->SetCentralityAxis(11, bins);
+ }
mgr->AddTask(task);
// create containers for input/output
fUseShapeCorr(true),
fSNNString(0),
fSysString(0),
- fCent(0)
+ fCent(0),
+ fCentAxis(0)
{
//
// Constructor
fUseShapeCorr(true),
fSNNString(0),
fSysString(0),
- fCent(0)
+ fCent(0),
+ fCentAxis(0)
{
//
// Constructor
//
- fListOfCentralities = new TList;
+ fListOfCentralities = new TObjArray(1);
// Output slot #1 writes into a TH1 container
DefineOutput(1, TList::Class());
fUseShapeCorr(o.fUseShapeCorr),
fSNNString(o.fSNNString),
fSysString(o.fSysString),
- fCent(o.fCent)
+ fCent(o.fCent),
+ fCentAxis(o.fCentAxis)
{}
//____________________________________________________________________
//________________________________________________________________________
void
-AliBasedNdetaTask::AddCentralityBin(Short_t low, Short_t high)
+AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
+{
+ if (!fCentAxis) {
+ fCentAxis = new TAxis();
+ fCentAxis->SetName("centAxis");
+ fCentAxis->SetTitle("Centrality [%]");
+ }
+ TArrayD dbins(n+1);
+ for (UShort_t i = 0; i <= n; i++)
+ dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
+ fCentAxis->Set(n, dbins.GetArray());
+}
+
+//________________________________________________________________________
+void
+AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
{
//
// Add a centrality bin
// low Low cut
// high High cut
//
- fListOfCentralities->Add(MakeCentralityBin(GetName(), low, high));
+ CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
+ AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
+ fListOfCentralities->AddAtAndExpand(bin, at);
}
//________________________________________________________________________
fSums->SetOwner();
// Automatically add 'all' centrality bin if nothing has been defined.
- if (fListOfCentralities->GetEntries() <= 0) AddCentralityBin(0,0);
+ AddCentralityBin(0, 0, 0);
+ if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
+ const TArrayD* bins = fCentAxis->GetXbins();
+ Int_t nbin = fCentAxis->GetNbins();
+ for (Int_t i = 0; i < nbin; i++)
+ AddCentralityBin(i+1, Short_t((*bins)[i]), Short_t((*bins)[i+1]));
+ }
+ fSums->Add(fCentAxis);
+
// Centrality histogram
fCent = new TH1D("cent", "Centrality", 100, 0, 100);
return;
}
AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
-
+
// Fill centrality histogram
- fCent->Fill(forward->GetCentrality());
+ Float_t cent = forward->GetCentrality();
+ fCent->Fill(cent);
// Get the histogram(s)
TH2D* data = GetHistogram(aod, false);
TH2D* dataMC = GetHistogram(aod, true);
// Loop over centrality bins
- TIter next(fListOfCentralities);
- CentralityBin* bin = 0;
- while ((bin = static_cast<CentralityBin*>(next())))
- bin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
+ CentralityBin* allBin = static_cast<CentralityBin*>(fListOfCentralities->At(0));
+ allBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
+
+ // Find this centrality bin
+ if (fCentAxis && fCentAxis->GetNbins() > 0) {
+ Int_t icent = fCentAxis->FindBin(cent);
+ CentralityBin* thisBin = 0;
+ if (icent >= 1 && icent <= fCentAxis->GetNbins())
+ thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
+ if (thisBin)
+ thisBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
+ }
// Here, we get the update
if (!fSNNString) {
fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
+ fCentAxis = static_cast<TAxis*>(fSums->FindObject("centAxis"));
if(fSysString && fSNNString &&
fSysString->GetUniqueID() == AliForwardUtil::kPP)
// Output collision system string
if (fSysString) fOutput->Add(fSysString->Clone());
+ // Output centrality axis
+ if (fCentAxis) fOutput->Add(fCentAxis);
+
// Output trigger string
TNamed* trigString =
new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
if (!forward) return false;
// Check the centrality class unless this is the 'all' bin
- if (!IsAllBin()) {
- Double_t centrality = forward->GetCentrality();
- if (centrality < fLow || centrality >= fHigh) return false;
- }
+ // if (!IsAllBin()) {
+ // Double_t centrality = forward->GetCentrality();
+ // if (centrality < fLow || centrality >= fHigh) return false;
+ // }
fTriggers->AddBinContent(kAll);
if (forward->IsTriggerBits(AliAODForwardMult::kB))
return;
}
if (!fSum) {
- AliError("Couldn't find histogram 'base' in list");
+ AliError(Form("Couldn't find histogram '%s' in list", GetName()));
return;
}
Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
Int_t nAccepted = Int_t(fTriggers->GetBinContent(kAccepted));
Int_t nGood = nB - nA - nC + 2 * nE;
+
- Double_t alpha = Double_t(nAccepted) / nWithVertex;
- Double_t vNorm = nAccepted + alpha*(nTriggered - nWithVertex);
- Double_t vtxEff = Double_t(nAccepted) / vNorm;
- vtxEff = vtxEff * trigEff;
-
- if (nTriggered <= 0) {
+ if (nTriggered <= 0.1) {
AliError("Number of triggered events <= 0");
return;
}
- if (nGood <= 0) {
+ if (nWithVertex <= 0.1) {
+ AliError("Number of events with vertex <= 0");
+ return;
+ }
+ if (nGood <= 0.1) {
AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
nGood, nB, nA, nC, nE));
nGood = nMB;
}
+
+
+ // Scaling
+ // N_A + N_A/N_V (N_T - N_V) = N_A + N_A/N_V*N_T - N_A
+ // = N_A/N_V * N_T
+ //
+ // where
+ // N_A = nAccepted
+ // N_V = nWithVertex
+ // N_T = nTriggered
+ // so that the normalisation is simply
+ // N_A / E_V
+ // where
+ // E_V=N_V/N_T is the vertex efficiency from data
+ //
+ // Double_t alpha = Double_t(nAccepted) / nWithVertex;
+ // Double_t vNorm = nAccepted + alpha*(nTriggered - nWithVertex);
+ AliInfo(Form("Calculating event normalisation as "
+ "1 / trigEff * nAccepted * nTriggered / nWithVertex "
+ "= 1/%f * %d * %d / %d = %f",
+ trigEff, nAccepted, nTriggered, nWithVertex,
+ 1. / trigEff * nAccepted * double(nTriggered) / nWithVertex));
+ Double_t vNorm = double(nWithVertex) / double(nTriggered);
+ Double_t vtxEff = vNorm;
+ Double_t ntotal = vtxEff * trigEff;
+
+
AliInfo(Form("Total of %9d events for %s\n"
- " of these %9d are minimum bias\n"
+ " of these %9d are triggered minimum bias\n"
" of these %9d has a %s trigger\n"
" of these %9d has a vertex\n"
" of these %9d were in [%+4.1f,%+4.1f]cm\n"
" A|C = %9d (%9d+%-9d)\n"
" E = %9d\n"
" Implies %9d good triggers\n"
- " Vertex efficiency: %f (%f)",
+ " Vertex efficiency: %f\n"
+ " Trigger efficiency: %f\n"
+ " Total number of events: %f\n",
nAll, GetTitle(), nMB, nTriggered,
AliAODForwardMult::GetTriggerString(triggerMask),
nWithVertex, nAccepted,
vzMin, vzMax,
- nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
+ nB, nA+nC, nA, nC, nE, nGood, vtxEff, trigEff, ntotal));
const char* name = GetName();
dndeta->Divide(norm);
// Scale by the vertex efficiency
- dndeta->Scale(vtxEff, "width");
+ dndeta->Scale(ntotal, "width");
SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name));
if (symmetrice)
fOutput->Add(Symmetrice(Rebin(dndetaMC, rebin, cutEdges)));
}
+ if (!IsAllBin()) return;
+
+ // Temporary stuff
+ TFile* forward = TFile::Open("forward.root", "READ");
+ if (!forward) {
+ AliWarning(Form("No forward.root file found"));
+ return;
+ }
+
+ TH1D* shapeCorrProj = 0;
+ if (shapeCorr) {
+ shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
+ shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
+ shapeCorrProj->SetDirectory(0);
+ fOutput->Add(shapeCorrProj);
+ }
+
+ TList* official = static_cast<TList*>(forward->Get("official"));
+ if (official) {
+ TH1F* histEta = static_cast<TH1F*>(official->FindObject("fHistEta"));
+ if (histEta) {
+ TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
+ histEta->GetXaxis()->GetXmin(),
+ histEta->GetXaxis()->GetXmax());
+ for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
+ oEta->SetBinContent(i, histEta->GetBinContent(i));
+ oEta->SetBinError(i, histEta->GetBinError(i));
+ }
+ if (shapeCorrProj) oEta->Divide(shapeCorrProj);
+ oEta->Scale(ntotal/nAccepted, "width");
+ oEta->SetDirectory(0);
+ oEta->SetMarkerStyle(21);
+ oEta->SetMarkerColor(kCyan+3);
+ fOutput->Add(oEta);
+ fOutput->Add(Rebin(oEta, rebin, false));
+ }
+ else
+ AliWarning(Form("Couldn't find histogram fHistEta in list %s",
+ official->GetName()));
+ }
+ else
+ AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
+
+ TList* tracks = static_cast<TList*>(forward->Get("tracks"));
+ if (tracks) {
+ TH1F* histEta = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
+ if (histEta) {
+ TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
+ histEta->GetXaxis()->GetXmin(),
+ histEta->GetXaxis()->GetXmax());
+ for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
+ oEta->SetBinContent(i, histEta->GetBinContent(i));
+ oEta->SetBinError(i, histEta->GetBinError(i));
+ }
+ if (shapeCorrProj) oEta->Divide(shapeCorrProj);
+ oEta->Scale(ntotal/nAccepted, "width");
+ oEta->SetDirectory(0);
+ oEta->SetMarkerStyle(22);
+ oEta->SetMarkerColor(kMagenta+2);
+ fOutput->Add(oEta);
+ fOutput->Add(Rebin(oEta, rebin, false));
+ }
+ else
+ AliWarning(Form("Couldn't find histogram fHistEta in list %s",
+ tracks->GetName()));
+ }
+ else
+ AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
+
+ forward->Close();
}
//
#ifndef ALIBASEDNDETATASK_H
#define ALIBASEDNDETATASK_H
#include <AliAnalysisTaskSE.h>
+class TAxis;
class TList;
class TH2D;
class TH1D;
class AliAODEvent;
class AliAODForwardMult;
+class TObjArray;
/**
* Task to determine the
* @name Task configuration
*/
/**
- * Add a centrality bin
- *
- * @param low Low cut
- * @param high High cut
- */
- void AddCentralityBin(Short_t low, Short_t high);
- /**
* Set the vertex range to use
*
* @param min Minimum (in centermeter)
*/
void SetTriggerMask(const char* mask);
/**
+ * Set the centrality bins to use.
+ *
+ * @code
+ * UShort_t bins[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+ * task->SetCentralityBins(11, bins);
+ * @endcode
+ *
+ * @param n Number of bins (elements in @a bins minus 1)
+ * @param bins Bin limits
+ */
+ void SetCentralityAxis(UShort_t n, Short_t* bins);
+ /**
* Trigger efficiency for selected trigger(s)
*
* @param e Trigger efficiency
*/
virtual TH2D* GetHistogram(const AliAODEvent* aod, Bool_t mc=false) = 0;
/**
+ * Add a centrality bin
+ *
+ * @param low Low cut
+ * @param high High cut
+ */
+ void AddCentralityBin(UShort_t at, Short_t low, Short_t high);
+ /**
* Make a centrality bin
*
* @param name Name used for histograms
Bool_t fCorrEmpty; // Correct for empty bins
Double_t fTriggerEff; // Trigger efficiency for selected trigger(s)
TH1* fShapeCorr; // Shape correction
- TList* fListOfCentralities; // Centrality bins
+ TObjArray* fListOfCentralities; // Centrality bins
Bool_t fUseShapeCorr; // Whether to use shape correction
TNamed* fSNNString; // sqrt(s_NN) string
TNamed* fSysString; // Collision system string
TH1D* fCent; // Centrality distribution
+ TAxis* fCentAxis; // Centrality axis
ClassDef(AliBasedNdetaTask,2); // Determine multiplicity in base area
};
// name Name of object
//
fRingHistos.SetName(GetName());
+ fRingHistos.SetOwner();
fRingHistos.Add(new RingHistos(1, 'I'));
fRingHistos.Add(new RingHistos(2, 'I'));
fRingHistos.Add(new RingHistos(2, 'O'));
case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
}
- if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
+ if (idx < 0 || idx >= fRingHistos.GetEntries()) {
+ AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
+ return 0;
+ }
return static_cast<RingHistos*>(fRingHistos.At(idx));
}
UShort_t nt= (q == 0 ? 512 : 256);
TH2D* h = hists.Get(d,r);
RingHistos* rh= GetRingHistos(d,r);
- rh->fEmptyStrips->Reset();
- rh->fTotalStrips->Reset();
- rh->fBasicHits->Reset();
+ if (!rh) {
+ AliError(Form("No ring histogram found for FMD%d%c", d, r));
+ fRingHistos.ls();
+ return false;
+ }
+ rh->ResetPoissonHistos();
for (UShort_t s=0; s<ns; s++) {
for (UShort_t t=0; t<nt; t++) {
fEmptyVsTotal->SetDirectory(0);
fEmptyVsTotal->SetXTitle("Total # strips");
fEmptyVsTotal->SetYTitle("Empty # strips");
- fEmptyVsTotal->SetZTitle("Correlation");
-
- //Inserted by HHD
-
- //Float_t nbinlimitsFMD3I[8] = {-3.4,-3.2,-3,-2.8,-2.6,-2.4,-2.2,-2};
- Int_t nbins = 40;
- /*if(fName.Contains("FMD1I")) nbins = 7;
- if(fName.Contains("FMD2I")) nbins = 8;
- if(fName.Contains("FMD2O")) nbins = 4;
- if(fName.Contains("FMD3I")) nbins = 8;
- if(fName.Contains("FMD3O")) nbins = 4;*/
- Float_t lowlimit = -4, highlimit = 6;
- /*if(fName.Contains("FMD1I")) {lowlimit = 3.6; highlimit = 5;}
- if(fName.Contains("FMD2I")) {lowlimit = 2.2; highlimit = 3.8;}
- if(fName.Contains("FMD2O")) {lowlimit = 1.6; highlimit = 2.4;}
- if(fName.Contains("FMD3I")) {lowlimit = -2.4; highlimit = -1.6;}
- if(fName.Contains("FMD3O")) {lowlimit = -3.5; highlimit = -2.1;}
-
- std::cout<<nbins<<" "<<lowlimit<<" "<<highlimit<<std::endl;
- */
- fTotalStrips = new TH2D(Form("total%s", fName.Data()),
- Form("Total number of strips in %s", fName.Data()),
- nbins,
- lowlimit,
- highlimit,
- (fRing == 'I' || fRing == 'i' ? 5 : 10),
- 0, 2*TMath::Pi());
- fEmptyStrips = new TH2D(Form("empty%s", fName.Data()),
- Form("Empty number of strips in %s", fName.Data()),
- nbins,
- lowlimit,
- highlimit,
- (fRing == 'I' || fRing == 'i' ? 5 : 10),
- 0, 2*TMath::Pi());
- fBasicHits = new TH2D(Form("basichits%s", fName.Data()),
- Form("Basic number of hits in %s", fName.Data()),
- 200,
- -4,
- 6,
- (fRing == 'I' || fRing == 'i' ? 20 : 40),
- 0, 2*TMath::Pi());
-
- fTotalStrips->SetDirectory(0);
- fEmptyStrips->SetDirectory(0);
- fBasicHits->SetDirectory(0);
- fTotalStrips->SetXTitle("#eta");
- fEmptyStrips->SetXTitle("#eta");
- fBasicHits->SetXTitle("#eta");
- fTotalStrips->SetYTitle("#varphi [radians]");
- fEmptyStrips->SetYTitle("#varphi [radians]");
- fBasicHits->SetYTitle("#varphi [radians]");
- fTotalStrips->Sumw2();
- fEmptyStrips->Sumw2();
- fBasicHits->Sumw2();
-
+ fEmptyVsTotal->SetZTitle("Correlation");
}
//____________________________________________________________________
AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
if (fEmptyStrips) delete fEmptyStrips;
}
- //____________________________________________________________________
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::RingHistos::ResetPoissonHistos()
+{
+ if (fTotalStrips) {
+ fTotalStrips->Reset();
+ fEmptyStrips->Reset();
+ fBasicHits->Reset();
+ return;
+ }
+ //Inserted by HHD
+
+ //Float_t nbinlimitsFMD3I[8] = {-3.4,-3.2,-3,-2.8,-2.6,-2.4,-2.2,-2};
+ Int_t nbins = 40;
+ /*if(fName.Contains("FMD1I")) nbins = 7;
+ if(fName.Contains("FMD2I")) nbins = 8;
+ if(fName.Contains("FMD2O")) nbins = 4;
+ if(fName.Contains("FMD3I")) nbins = 8;
+ if(fName.Contains("FMD3O")) nbins = 4;*/
+ Float_t lowlimit = -4, highlimit = 6;
+ /*if(fName.Contains("FMD1I")) {lowlimit = 3.6; highlimit = 5;}
+ if(fName.Contains("FMD2I")) {lowlimit = 2.2; highlimit = 3.8;}
+ if(fName.Contains("FMD2O")) {lowlimit = 1.6; highlimit = 2.4;}
+ if(fName.Contains("FMD3I")) {lowlimit = -2.4; highlimit = -1.6;}
+ if(fName.Contains("FMD3O")) {lowlimit = -3.5; highlimit = -2.1;}
+
+ std::cout<<nbins<<" "<<lowlimit<<" "<<highlimit<<std::endl;
+ */
+ fTotalStrips = new TH2D(Form("total%s", fName.Data()),
+ Form("Total number of strips in %s", fName.Data()),
+ nbins,
+ lowlimit,
+ highlimit,
+ (fRing == 'I' || fRing == 'i' ? 5 : 10),
+ 0, 2*TMath::Pi());
+ fEmptyStrips = new TH2D(Form("empty%s", fName.Data()),
+ Form("Empty number of strips in %s", fName.Data()),
+ nbins,
+ lowlimit,
+ highlimit,
+ (fRing == 'I' || fRing == 'i' ? 5 : 10),
+ 0, 2*TMath::Pi());
+ fBasicHits = new TH2D(Form("basichits%s", fName.Data()),
+ Form("Basic number of hits in %s", fName.Data()),
+ 200,
+ -4,
+ 6,
+ (fRing == 'I' || fRing == 'i' ? 20 : 40),
+ 0, 2*TMath::Pi());
+
+ fTotalStrips->SetDirectory(0);
+ fEmptyStrips->SetDirectory(0);
+ fBasicHits->SetDirectory(0);
+ fTotalStrips->SetXTitle("#eta");
+ fEmptyStrips->SetXTitle("#eta");
+ fBasicHits->SetXTitle("#eta");
+ fTotalStrips->SetYTitle("#varphi [radians]");
+ fEmptyStrips->SetYTitle("#varphi [radians]");
+ fBasicHits->SetYTitle("#varphi [radians]");
+ fTotalStrips->Sumw2();
+ fEmptyStrips->Sumw2();
+ fBasicHits->Sumw2();
+}
+
+//____________________________________________________________________
void
AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
{
* @param nEvents Number of events
*/
void ScaleHistograms(TList* dir, Int_t nEvents);
+ /**
+ * Create Poisson histograms
+ */
+ void ResetPoissonHistos();
TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
TProfile* fEtaVsN; // Average uncorrected Nch vs eta
TProfile* fCorr; // Average correction vs eta
TH2D* fDensity; // Distribution inclusive Nch
TH2D* fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
- TH2D* fTotalStrips; //! Total number of strips in a region
- TH2D* fEmptyStrips; //! Total number of strips in a region
- TH2D* fBasicHits ; //! Total number basic hits in a region
+ TH2D* fTotalStrips; // Total number of strips in a region
+ TH2D* fEmptyStrips; // Total number of strips in a region
+ TH2D* fBasicHits ; // Total number basic hits in a region
TH2D* fEmptyVsTotal; // # of empty strips vs total number of strips
- ClassDef(RingHistos,1);
+ ClassDef(RingHistos,3);
};
/**
* Get the ring histogram container
fHType(0),
fHWords(0),
fHCent(0),
+ fHCentVsQual(0),
fLowFluxCut(1000),
fMaxVzErr(0.2),
fList(0),
fEnergy(0),
fField(999),
fCollisionSystem(kUnknown),
- fDebug(0)
+ fDebug(0),
+ fCentAxis(0)
{
//
// Constructor
fHType(0),
fHWords(0),
fHCent(0),
+ fHCentVsQual(0),
fLowFluxCut(1000),
fMaxVzErr(0.2),
fList(0),
fEnergy(0),
fField(999),
fCollisionSystem(kUnknown),
- fDebug(0)
+ fDebug(0),
+ fCentAxis(0)
{
//
// Constructor
fHType(o.fHType),
fHWords(o.fHWords),
fHCent(o.fHCent),
+ fHCentVsQual(o.fHCentVsQual),
fLowFluxCut(o.fLowFluxCut),
fMaxVzErr(o.fMaxVzErr),
fList(o.fList),
fEnergy(o.fEnergy),
fField(o.fField),
fCollisionSystem(o.fCollisionSystem),
- fDebug(0)
+ fDebug(0),
+ fCentAxis(0)
{
//
// Copy constructor
//
// Destructor
//
- if (fHEventsTr) delete fHEventsTr;
- if (fHEventsTrVtx) delete fHEventsTrVtx;
- if (fHTriggers) delete fHTriggers;
- if (fHType) delete fHType;
- if (fHWords) delete fHWords;
- if (fHCent) delete fHCent;
if (fList) delete fList;
}
//____________________________________________________________________
fHType = o.fHType;
fHWords = o.fHWords;
fHCent = o.fHCent;
+ fHCentVsQual = o.fHCentVsQual;
fLowFluxCut = o.fLowFluxCut;
fMaxVzErr = o.fMaxVzErr;
fDebug = o.fDebug;
if (fHType) fList->Add(fHType);
if (fHWords) fList->Add(fHWords);
if (fHCent) fList->Add(fHCent);
+ if (fHCentVsQual) fList->Add(fHCentVsQual);
}
return *this;
}
// Parameters:
// vtxAxis Vertex axis in use
//
+
+ // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
+ // ----- 92 number --------- ---- 1 ---
+ TArrayD limits(93);
+ for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
+
+ fCentAxis = new TAxis(limits.GetSize()-1, limits.GetArray());
fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
vtxAxis.GetNbins(),
vtxAxis.GetXmin(),
fHWords->SetBit(TH1::kCanRebin);
fList->Add(fHWords);
- fHCent = new TH1F("cent", "Centrality", 101, -1.5, 100.5);
+ fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
fHCent->SetFillColor(kBlue+1);
fHCent->SetFillStyle(3001);
fHCent->SetStats(0);
fHCent->SetXTitle("Centrality [%]");
fHCent->SetYTitle("Events");
fList->Add(fHCent);
+
+ fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality",
+ 5, 0, 5, limits.GetSize()-1, limits.GetArray());
+ fHCentVsQual->SetXTitle("Quality");
+ fHCentVsQual->SetYTitle("Centrality [%]");
+ fHCentVsQual->SetZTitle("Events");
+ fHCentVsQual->GetXaxis()->SetBinLabel(1, "OK");
+ fHCentVsQual->GetXaxis()->SetBinLabel(2, "Outside v_{z} cut");
+ fHCentVsQual->GetXaxis()->SetBinLabel(3, "V0 vs SPD outlier");
+ fHCentVsQual->GetXaxis()->SetBinLabel(4, "V0 vs TPC outlier");
+ fHCentVsQual->GetXaxis()->SetBinLabel(5, "V0 vs ZDC outlier");
+ fList->Add(fHCentVsQual);
}
//____________________________________________________________________
fHType->Fill(lowFlux ? 0 : 1);
// --- Read centrality information
- cent = -10;
- if (!ReadCentrality(event, cent)) {
+ cent = -10;
+ UShort_t qual = 0;
+ if (!ReadCentrality(event, cent, qual)) {
if (fDebug > 3)
AliWarning("Failed to get centrality");
}
+ fHCent->Fill(cent);
+ if (qual == 0) fHCentVsQual->Fill(0., cent);
+ else {
+ for (UShort_t i = 0; i < 4; i++)
+ if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
+ }
// --- Get the vertex information ----------------------------------
vz = 0;
//____________________________________________________________________
Bool_t
-AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
+AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd,
+ Double_t& cent,
+ UShort_t& qual) const
{
//
// Read centrality from event
// Return:
// False on error, true otherwise
//
+ cent = -1;
+ qual = 0;
AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
- if (centObj) {
- // AliInfo(Form("Got centrality object %p with quality %d",
- // centObj, centObj->GetQuality()));
- // centObj->Print();
- cent = centObj->GetCentralityPercentile("V0M");
- }
- // AliInfo(Form("Centrality is %f", cent));
- fHCent->Fill(cent);
+ if (!centObj) return true;
+
+ // AliInfo(Form("Got centrality object %p with quality %d",
+ // centObj, centObj->GetQuality()));
+ // centObj->Print();
+ cent = centObj->GetCentralityPercentile("V0M");
+ qual = centObj->GetQuality();
return true;
}
class TH1D;
class TH1I;
class TH1F;
+class TH2F;
class TAxis;
class TList;
*
* @return False on error, true otherwise
*/
- virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent);
+ virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent,
+ UShort_t& qual) const;
TH1I* fHEventsTr; //! Histogram of events w/trigger
TH1I* fHEventsTrVtx; //! Events w/trigger and vertex
TH1I* fHType; //! Type (low/high flux) of event
TH1I* fHWords; //! Trigger words
TH1F* fHCent; //! Centrality
+ TH2F* fHCentVsQual; //! Centrality vs quality
Int_t fLowFluxCut; // Low flux cut
Double_t fMaxVzErr; // Maximum error on v_z
TList* fList; //! Histogram container
Short_t fField; // L3 magnetic field [kG]
UShort_t fCollisionSystem; // Collision system
Int_t fDebug; // Debug level
- ClassDef(AliFMDEventInspector,2); // Inspect the event
+ TAxis* fCentAxis; // Centrality axis used in histograms
+ ClassDef(AliFMDEventInspector,3); // Inspect the event
};
#endif
#include "AliGenGeVSimEventHeader.h"
#include "AliHeader.h"
#include <TList.h>
+#include <TH2F.h>
+
//====================================================================
AliFMDMCEventInspector::AliFMDMCEventInspector()
: AliFMDEventInspector(),
fHVertex(0),
fHPhiR(0),
- fHB(0)
+ fHB(0),
+ fHBvsPart(0),
+ fHBvsBin(0),
+ fHBvsCent(0),
+ fHVzComp(0),
+ fHCentVsPart(0),
+ fHCentVsBin(0)
{
//
// Constructor
: AliFMDEventInspector("fmdEventInspector"),
fHVertex(0),
fHPhiR(0),
- fHB(0)
+ fHB(0),
+ fHBvsPart(0),
+ fHBvsBin(0),
+ fHBvsCent(0),
+ fHVzComp(0),
+ fHCentVsPart(0),
+ fHCentVsBin(0)
{
//
// Constructor
: AliFMDEventInspector(o),
fHVertex(0),
fHPhiR(0),
- fHB(0)
+ fHB(0),
+ fHBvsPart(0),
+ fHBvsBin(0),
+ fHBvsCent(0),
+ fHVzComp(0),
+ fHCentVsPart(0),
+ fHCentVsBin(0)
{
//
// Copy constructor
//
AliFMDEventInspector::Init(vtxAxis);
- fHVertex = new TH1F("vertex", "True vertex distribution",
- vtxAxis.GetNbins(),
- vtxAxis.GetXmin(),
- vtxAxis.GetXmax());
+ Int_t maxPart = 450;
+ Int_t maxBin = 225;
+ Int_t maxB = 25;
+ Int_t nVtx = vtxAxis.GetNbins();
+ Double_t lVtx = vtxAxis.GetXmin();
+ Double_t hVtx = vtxAxis.GetXmax();
+ fHVertex = new TH1F("vertex", "True vertex distribution", nVtx, lVtx, hVtx);
fHVertex->SetXTitle("v_{z} [cm]");
fHVertex->SetYTitle("# of events");
fHVertex->SetFillColor(kGreen+1);
fHPhiR->SetDirectory(0);
fList->Add(fHPhiR);
- fHB = new TH1F("b", "Impact parameter", 125, 0, 25);
+ fHB = new TH1F("b", "Impact parameter", 5*maxB, 0, maxB);
fHB->SetXTitle("b [fm]");
fHB->SetYTitle("# of events");
fHB->SetFillColor(kGreen+1);
fHB->SetFillStyle(3001);
fHB->SetDirectory(0);
fList->Add(fHB);
+
+ fHBvsPart = new TH2F("bVsParticipants", "Impact parameter vs Participants",
+ 5*maxB, 0, maxB, maxPart, -.5, maxPart-.5);
+ fHBvsPart->SetXTitle("b [fm]");
+ fHBvsPart->SetYTitle("# of participants");
+ fHBvsPart->SetZTitle("Events");
+ fHBvsPart->SetDirectory(0);
+ fList->Add(fHBvsPart);
+
+ fHBvsBin = new TH2F("bVsBinary", "Impact parameter vs Binary Collisions",
+ 5*maxB, 0, maxB, maxBin, -.5, maxBin-.5);
+ fHBvsBin->SetXTitle("b [fm]");
+ fHBvsBin->SetYTitle("# of binary collisions");
+ fHBvsBin->SetZTitle("Events");
+ fHBvsBin->SetDirectory(0);
+ fList->Add(fHBvsBin);
+ fHBvsCent = new TH2F("bVsCentrality", "Impact parameter vs Centrality",
+ 5*maxB, 0, maxB, fCentAxis->GetNbins(),
+ fCentAxis->GetXbins()->GetArray());
+ fHBvsCent->SetXTitle("b [fm]");
+ fHBvsCent->SetYTitle("Centrality [%]");
+ fHBvsCent->SetZTitle("Event");
+ fHBvsCent->SetDirectory(0);
+ fList->Add(fHBvsCent);
+
+
+ fHVzComp = new TH2F("vzComparison", "v_{z} truth vs reconstructed",
+ 10*nVtx, lVtx, hVtx, 10*nVtx, lVtx, hVtx);
+ fHVzComp->SetXTitle("True v_{z} [cm]");
+ fHVzComp->SetYTitle("Reconstructed v_{z} [cm]");
+ fHVzComp->SetZTitle("Events");
+ fHVzComp->SetDirectory(0);
+ fList->Add(fHVzComp);
+
+ fHCentVsPart = new TH2F("centralityVsParticipans",
+ "# of participants vs Centrality",
+ maxPart, -.5, maxPart-.5, fCentAxis->GetNbins(),
+ fCentAxis->GetXbins()->GetArray());
+ fHCentVsPart->SetXTitle("Participants");
+ fHCentVsPart->SetYTitle("Centrality [%]");
+ fHCentVsPart->SetZTitle("Event");
+ fHCentVsPart->SetDirectory(0);
+ fList->Add(fHCentVsPart);
+
+ fHCentVsBin = new TH2F("centralityVsBinary",
+ "# of binary collisions vs Centrality",
+ maxBin, -.5, maxBin-.5, fCentAxis->GetNbins(),
+ fCentAxis->GetXbins()->GetArray());
+ fHCentVsBin->SetXTitle("Binary collisions");
+ fHCentVsBin->SetYTitle("Centrality [%]");
+ fHCentVsBin->SetZTitle("Event");
+ fHCentVsBin->SetDirectory(0);
+ fList->Add(fHCentVsBin);
}
//____________________________________________________________________
UShort_t& ivz,
Double_t& vz,
Double_t& b,
+ Int_t& npart,
+ Int_t& nbin,
Double_t& phiR)
{
//
//Assign MC only triggers : True NSD etc.
AliHeader* header = event->Header();
AliGenEventHeader* genHeader = header->GenEventHeader();
+ AliCollisionGeometry* colGeometry =
+ dynamic_cast<AliCollisionGeometry*>(genHeader);
AliGenPythiaEventHeader* pythiaHeader =
dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
AliGenDPMjetEventHeader* dpmHeader =
dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
AliGenGeVSimEventHeader* gevHeader =
dynamic_cast<AliGenGeVSimEventHeader*>(genHeader);
- AliGenHijingEventHeader* hijingHeader =
- dynamic_cast<AliGenHijingEventHeader*>(genHeader);
AliGenHerwigEventHeader* herwigHeader =
dynamic_cast<AliGenHerwigEventHeader*>(genHeader);
+ // AliGenHijingEventHeader* hijingHeader =
+ // dynamic_cast<AliGenHijingEventHeader*>(genHeader);
// AliGenHydjetEventHeader* hydjetHeader =
// dynamic_cast<AliGenHydjetEventHeader*>(genHeader);
- AliGenEposEventHeader* eposHeader =
- dynamic_cast<AliGenEposEventHeader*>(genHeader);
+ // AliGenEposEventHeader* eposHeader =
+ // dynamic_cast<AliGenEposEventHeader*>(genHeader);
// Check if this is a single diffractive event
- Bool_t sd = kFALSE;
- Double_t phi = -1111;
- b = -1;
+ Bool_t sd = kFALSE;
+ Double_t phi = -1111;
+ npart = 0;
+ nbin = 0;
+ b = -1;
+ if (colGeometry) {
+ b = colGeometry->ImpactParameter();
+ phi = colGeometry->ReactionPlaneAngle();
+ npart = (colGeometry->ProjectileParticipants() +
+ colGeometry->TargetParticipants());
+ nbin = colGeometry->NN();
+ }
if(pythiaHeader) {
Int_t pythiaType = pythiaHeader->ProcessType();
if (pythiaType==92 || pythiaType==93) sd = kTRUE;
- b = pythiaHeader->GetImpactParameter();
+ b = pythiaHeader->GetImpactParameter();
+ npart = 2; // Always 2 protons
+ nbin = 1; // Always 1 binary collision
}
- if(dpmHeader) {
+ if(dpmHeader) { // Also an AliCollisionGeometry
Int_t processType = dpmHeader->ProcessType();
if (processType == 5 || processType == 6) sd = kTRUE;
- b = dpmHeader->ImpactParameter();
- phi = dpmHeader->ReactionPlaneAngle();
-
}
if (gevHeader) {
phi = gevHeader->GetEventPlane();
}
- if (hijingHeader) {
- b = hijingHeader->ImpactParameter();
- phi = hijingHeader->ReactionPlaneAngle();
- }
if (herwigHeader) {
Int_t processType = herwigHeader->ProcessType();
// This is a guess
if (processType == 5 || processType == 6) sd = kTRUE;
+ npart = 2; // Always 2 protons
+ nbin = 1; // Always 1 binary collision
}
+ // if (hijingHeader) {
+ // b = hijingHeader->ImpactParameter();
+ // phi = hijingHeader->ReactionPlaneAngle();
+ // }
// if (hydjetHeader) {
// b = hydjetHeader->ImpactParameter();
// phi = hydjetHeader->ReactionPlaneAngle();
// }
- if (eposHeader) {
- b = eposHeader->ImpactParameter();
- phi = eposHeader->ReactionPlaneAngle();
- }
+ // if (eposHeader) {
+ // b = eposHeader->ImpactParameter();
+ // phi = eposHeader->ReactionPlaneAngle();
+ // }
// Normalize event plane angle to [0,2pi]
if (phi <= -1111) phiR = -1;
fHVertex->Fill(vz);
fHPhiR->Fill(phiR);
fHB->Fill(b);
+ fHBvsPart->Fill(b, npart);
+ fHBvsBin->Fill(b, nbin);
// Check for the vertex bin
ivz = fHEventsTr->GetXaxis()->FindBin(vz);
}
//____________________________________________________________________
Bool_t
-AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd, Double_t& cent)
+AliFMDMCEventInspector::ReadCentrality(const AliESDEvent* esd,
+ Double_t& cent,
+ UShort_t& qual) const
{
//
// Read centrality from event
// Return:
// False on error, true otherwise
//
+ cent = -1;
+ qual = 0;
AliCentrality* centObj = const_cast<AliESDEvent*>(esd)->GetCentrality();
- if (centObj) {
- // AliInfo(Form("Got centrality object %p with quality %d",
- // centObj, centObj->GetQuality()));
- // centObj->Print();
- if (centObj->GetQuality() == 0x8)
- cent = centObj->GetCentralityPercentileUnchecked("V0M");
- else
- cent = centObj->GetCentralityPercentile("V0M");
- }
- // AliInfo(Form("Centrality is %f", cent));
- fHCent->Fill(cent);
+ if (!centObj) return true;
+
+ qual = centObj->GetQuality();
+ if (qual == 0x8) // Ignore ZDC outliers
+ cent = centObj->GetCentralityPercentileUnchecked("V0M");
+ else
+ cent = centObj->GetCentralityPercentile("V0M");
return true;
}
-
+//____________________________________________________________________
+Bool_t
+AliFMDMCEventInspector::CompareResults(Double_t vz, Double_t trueVz,
+ Double_t cent, Double_t b,
+ Int_t npart, Int_t nbin)
+{
+ fHVzComp->Fill(trueVz, vz);
+ fHBvsCent->Fill(b, cent);
+ fHCentVsPart->Fill(npart, cent);
+ fHCentVsBin->Fill(nbin, cent);
+
+ return true;
+}
//
// EOF
//
#define ALIFMDMCEVENTINSPECTOR_H
#include "AliFMDEventInspector.h"
class AliMCEvent;
+class TH2F;
/**
* This class inspects the event
UShort_t& ivz,
Double_t& vz,
Double_t& b,
+ Int_t& npart,
+ Int_t& nbin,
Double_t& phiR);
+ /**
+ * Compare the result of analysing the ESD for
+ * the inclusive charged particle density to analysing
+ * MC truth
+ *
+ * @param esd
+ * @param mc
+ *
+ * @return
+ */
+ virtual Bool_t CompareResults(Double_t vz, Double_t trueVz,
+ Double_t cent, Double_t b,
+ Int_t npart, Int_t nbin);
protected:
/**
* Read centrality from event
*
* @return False on error, true otherwise
*/
- virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent);
+ virtual Bool_t ReadCentrality(const AliESDEvent* esd, Double_t& cent,
+ UShort_t& qual) const;
- TH1F* fHVertex; // Histogram of vertex
- TH1F* fHPhiR; // Histogram of event plane
- TH1F* fHB; // Histogram of impact parameter
- ClassDef(AliFMDMCEventInspector,1); // Inspect the event
+ TH1F* fHVertex; // Histogram of vertex
+ TH1F* fHPhiR; // Histogram of event plane
+ TH1F* fHB; // Histogram of impact parameter
+ TH2F* fHBvsPart; // Impact parameter vs # participants
+ TH2F* fHBvsBin; // Impact parameter vs # participants
+ TH2F* fHBvsCent; // Impact parameter vs centrality
+ TH2F* fHVzComp; // True vs reconstructed vz
+ TH2F* fHCentVsPart; // Centrality versus # participants
+ TH2F* fHCentVsBin; // Centrality versus # binary collisions
+ ClassDef(AliFMDMCEventInspector,2); // Inspect the event
};
#endif
Double_t vzMC = 0;
Double_t phiR = 0;
Double_t b = 0;
+ Int_t npart = 0;
+ Int_t nbin = 0;
// UInt_t foundMC =
- fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, phiR);
-
+ fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b,
+ npart, nbin, phiR);
+ fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
//Store all events
MarkEventForStore();
//MarkEventForStore();
// Set trigger bits, and mark this event for storage
fAODFMD.SetTriggerBits(triggers);
- fMCAODFMD.SetTriggerBits(triggers);
+ fAODFMD.SetSNN(fEventInspector.GetEnergy());
+ fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
fAODFMD.SetCentrality(cent);
+
+ fMCAODFMD.SetTriggerBits(triggers);
+ fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
+ fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
+ fMCAODFMD.SetCentrality(cent);
//All events should be stored - HHD
//if (isAccepted) MarkEventForStore();
*
* @param filename File containing the data
*/
- void Run(const char* filename="forward_dndeta.C")
+ void Run(const char* filename="forward_dndeta.root")
{
- Double_t max = 0;
+ Double_t max = 0, rmax=0, amax=0;
- if (!Open(filename, max)) return;
- Info("Run", "Got data from %s with maximum %f", filename, max);
+ gStyle->SetPalette(1);
- // Create our stack of results
- THStack* results = fResults; // StackResults(max);
+ // --- Open input file -------------------------------------------
+ TFile* file = TFile::Open(filename, "READ");
+ if (!file) {
+ Error("Open", "Cannot open %s", filename);
+ return;
+ }
+ // --- Get forward list ------------------------------------------
+ TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
+ if (!forward) {
+ Error("Open", "Couldn't find list ForwardResults");
+ return;
+ }
+ // --- Get information on the run --------------------------------
+ FetchInformation(forward);
+ // --- Set the macro pathand load other data script --------------
+ gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2",
+ gROOT->GetMacroPath()));
+ gROOT->LoadMacro("OtherData.C");
+
+ // --- Get the central results -----------------------------------
+ TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
+ if (!clusters) Warning("Open", "Couldn't find list CentralResults");
- // Create our stack of other results
- TMultiGraph* other = 0;
- if (fShowOthers || fShowAlice) other = StackOther(max);
+ // --- Make our containtes ---------------------------------------
+ fResults = new THStack("results", "Results");
+ fRatios = new THStack("ratios", "Ratios");
+ fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
+ fOthers = new TMultiGraph();
- Double_t smax = 0;
- THStack* ratios = 0;
- if (fShowRatios) ratios = StackRatios(other, smax);
+ // --- Loop over input data --------------------------------------
+ FetchResults(forward, "Forward", max, rmax, amax);
+ FetchResults(clusters, "Central", max, rmax, amax);
- Double_t amax = 0;
- THStack* leftright = 0;
- if (fShowLeftRight) leftright = StackLeftRight(amax);
+ // --- Get trigger information -----------------------------------
+ TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
+ if (sums) {
+ TList* all = static_cast<TList*>(sums->FindObject("all"));
+ if (all) {
+ fTriggers = FetchResult(all, "triggers");
+ if (!fTriggers) all->ls();
+ }
+ else {
+ Warning("Run", "List all not found in ForwardSums");
+ sums->ls();
+ }
+ }
+ else {
+ Warning("Run", "No ForwardSums directory found in %s", file->GetName());
+ file->ls();
+ }
+
+ // --- Check our stacks ------------------------------------------
+ if (!fResults->GetHists() ||
+ fResults->GetHists()->GetEntries() <= 0) {
+ Error("Run", "No histograms in result stack!");
+ return;
+ }
+ if (!fOthers->GetListOfGraphs() ||
+ fOthers->GetListOfGraphs()->GetEntries() <= 0) {
+ Warning("Run", "No other data found - disabling that");
+ fShowOthers = false;
+ }
+ if (!fRatios->GetHists() ||
+ fRatios->GetHists()->GetEntries() <= 0) {
+ Warning("Run", "No ratio data found - disabling that");
+ // fRatios->ls();
+ fShowRatios = false;
+ }
+ if (!fLeftRight->GetHists() ||
+ fLeftRight->GetHists()->GetEntries() <= 0) {
+ Warning("Run", "No left/right data found - disabling that");
+ // fLeftRight->ls();
+ fShowLeftRight = false;
+ }
+
+ // --- Close the input file --------------------------------------
+ file->Close();
+
+
- Plot(results, other, max, ratios, smax, leftright, amax);
+ // --- Plot results ----------------------------------------------
+ Plot(max, rmax, amax);
}
//__________________________________________________________________
fSysString = static_cast<TNamed*>(results->FindObject("sys"));
if (!fVtxAxis)
fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
+ if (!fCentAxis)
+ fCentAxis = static_cast<TAxis*>(results->FindObject("centAxis"));
if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
if (!fSNNString) fSNNString = new TNamed("sNN", "unknown");
if (!fSysString) fSysString = new TNamed("sys", "unknown");
fVtxAxis->SetTitle("v_{z} range unspecified");
}
+ TString centTxt;
+ if (fCentAxis) {
+ Int_t nCent = fCentAxis->GetNbins();
+ centTxt = Form("\n Centrality: %d bins", nCent);
+ for (Int_t i = 0; i <= nCent; i++)
+ centTxt.Append(Form("%c%d", i == 0 ? ' ' : ',',
+ fCentAxis->GetXbins()->At(i)));
+ }
Info("FetchInformation",
"Initialized for\n"
" Trigger: %s (%d)\n"
" sqrt(sNN): %s (%dGeV)\n"
" System: %s (%d)\n"
- " Vz range: %s (%f,%f)",
+ " Vz range: %s (%f,%f)%s",
fTrigString->GetTitle(), fTrigString->GetUniqueID(),
fSNNString->GetTitle(), fSNNString->GetUniqueID(),
fSysString->GetTitle(), fSysString->GetUniqueID(),
- fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
+ fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
+ centTxt.Data());
}
+ //__________________________________________________________________
+ TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
+ {
+ TMultiGraph* thisOther = 0;
+ if (!fShowOthers && !fShowAlice) return 0;
+
+ Bool_t onlya = (fShowOthers ? false : true);
+ UShort_t sys = (fSysString ? fSysString->GetUniqueID() : 0);
+ UShort_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
+ UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
+ Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d,%d,%d,%d);",
+ sys,snn,trg,
+ centLow,centHigh,onlya));
+ if (!ret) {
+ Warning("FetchResults", "No other data found for sys=%d, sNN=%d, "
+ "trigger=%d %d%%-%d%% central %s",
+ sys, snn, trg, centLow, centHigh,
+ onlya ? " (ALICE results)" : "all");
+ return 0;
+ }
+ thisOther = reinterpret_cast<TMultiGraph*>(ret);
+ return thisOther;
+ }
//__________________________________________________________________
/**
- * Open input file, and find data
+ * Get the results from the top-level list
*
- * @param filename File name
- *
- * @return true on success
+ * @param list List
+ * @param name name
+ * @param max On return, maximum of data
+ * @param rmax On return, maximum of ratios
+ * @param amax On return, maximum of left-right comparisons
*/
- Bool_t Open(const char* filename, Double_t& max)
+ void FetchResults(const TList* list,
+ const char* name,
+ Double_t& max,
+ Double_t& rmax,
+ Double_t& amax)
{
- TFile* file = TFile::Open(filename, "READ");
- if (!file) {
- Error("Open", "Cannot open %s", filename);
- return false;
+ UShort_t n = fCentAxis ? fCentAxis->GetNbins() : 0;
+ if (n == 0) {
+ TList* all = static_cast<TList*>(list->FindObject("all"));
+ if (!all)
+ Error("FetchResults", "Couldn't find list 'all' in %s",
+ list->GetName());
+ else
+ FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax);
+ return;
}
- TList* results = static_cast<TList*>(file->Get("ForwardResults"));
- if (!results) {
- Error("Open", "Couldn't find list ForwardResults");
- return false;
+ Int_t nCol = gStyle->GetNumberOfColors();
+ for (UShort_t i = 0; i < n; i++) {
+ UShort_t centLow = fCentAxis->GetXbins()->At(i);
+ UShort_t centHigh = fCentAxis->GetXbins()->At(i+1);
+ TString lname = Form("cent%03d_%03d", centLow, centHigh);
+ TList* thisCent = static_cast<TList*>(list->FindObject(lname));
+
+ Float_t fc = (centLow+double(centHigh-centLow)/2) / 100;
+ Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
+ Int_t col = gStyle->GetColorPalette(icol);
+ Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d",
+ centLow, centHigh, icol, nCol, fc, col);
+
+ TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh);
+ if (!thisCent)
+ Error("FetchResults", "Couldn't find list '%s' in %s",
+ lname.Data(), list->GetName());
+ else
+ FetchResults(thisCent, name, FetchOthers(centLow, centHigh),
+ col, centTxt.Data(), max, rmax, amax);
}
-
- // Get information on the run
- FetchInformation(results);
-
- TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
- if (!clusters) Warning("Open", "Couldn't find list CentralResults");
-
- fResults = new THStack("results", "Results");
- FetchResults(fResults, results, "Forward", max);
- FetchResults(fResults, clusters, "Central", max);
-
- TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
- if (sums) fTriggers = FetchResult(sums, "triggers");
-
- file->Close();
-
- return true;
+ }
+ //__________________________________________________________________
+ void SetAttributes(TH1* h, Int_t color)
+ {
+ if (!h) return;
+ if (color < 0) return;
+ h->SetLineColor(color);
+ h->SetMarkerColor(color);
+ h->SetFillColor(color);
+ }
+ //__________________________________________________________________
+ void SetAttributes(TGraph* g, Int_t color)
+ {
+ if (!g) return;
+ if (color < 0) return;
+ g->SetLineColor(color);
+ g->SetMarkerColor(color);
+ g->SetFillColor(color);
}
//__________________________________________________________________
- void FetchResults(THStack* stack, const TList* list, const char* name,
- Double_t& max)
+ void ModifyTitle(TNamed* h, const char* centTxt)
+ {
+ if (!centTxt || !h) return;
+ h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
+ }
+
+ //__________________________________________________________________
+ /**
+ * Fetch results for a particular centrality bin
+ *
+ * @param list List
+ * @param name Name
+ * @param centLow Low cut on centrality
+ * @param centHigh high cut on centrality
+ * @param max On return, data maximum
+ * @param rmax On return, ratio maximum
+ * @param amax On return, left-right maximum
+ */
+ void FetchResults(const TList* list,
+ const char* name,
+ TMultiGraph* thisOther,
+ Int_t color,
+ const char* centTxt,
+ Double_t& max,
+ Double_t& rmax,
+ Double_t& amax)
{
TH1* dndeta = FetchResult(list, Form("dndeta%s", name));
TH1* dndetaMC = FetchResult(list, Form("dndeta%sMC", name));
TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
TH1* dndetaSym = 0;
TH1* dndetaMCSym = 0;
+ TH1* tracks = FetchResult(list, "tracks");
+ SetAttributes(dndeta, color);
+ SetAttributes(dndetaMC, color+2);
+ SetAttributes(dndetaTruth,color);
+ SetAttributes(dndetaSym, color);
+ SetAttributes(dndetaMCSym,color+2);
+ SetAttributes(tracks, color+3);
+ ModifyTitle(dndeta, centTxt);
+ ModifyTitle(dndetaMC, centTxt);
+ ModifyTitle(dndetaTruth,centTxt);
+ ModifyTitle(dndetaSym, centTxt);
+ ModifyTitle(dndetaMCSym,centTxt);
+ ModifyTitle(tracks, centTxt);
+
- max = TMath::Max(max, AddHistogram(stack, dndetaTruth, "e5 p"));
- max = TMath::Max(max, AddHistogram(stack, dndetaMC, "", dndetaMCSym));
- max = TMath::Max(max, AddHistogram(stack, dndeta, "", dndetaSym));
+ max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5 p"));
+ max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym));
+ max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym));
+ max = TMath::Max(max, AddHistogram(fResults, tracks));
- Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
- dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
-
- if (dndetaMCSym) fNumerators.Add(dndetaMCSym);
- if (dndetaMC) fNumerators.Add(dndetaMC);
- if (dndetaSym) fNumerators.Add(dndetaSym);
- if (dndeta) fNumerators.Add(dndeta);
- if (dndetaTruth) fDenominators.Add(dndetaTruth);
- }
- //__________________________________________________________________
- /**
- * Make a histogram stack of results
- *
- * @param max On return, the maximum value in the stack
- *
- * @return Newly allocated stack
- */
- TMultiGraph* StackOther(Double_t& max)
- {
- gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
- Int_t error = 0;
- Bool_t onlya = (fShowOthers ? false : true);
- Int_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
- UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
- Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d);",
- snn,trg,onlya));
- if (error) {
- Error("StackOther", "Failed to execute GetData(%d,%d,%d)",
- snn, trg, onlya);
- return 0;
- }
- if (!ret) {
- Warning("StackOther", "No other data found for sNN=%d, trigger=%d",
- snn, trg);
- return 0;
- }
- TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret);
+ // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
+ // dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
- TGraphAsymmErrors* o = 0;
- TIter next(other->GetListOfGraphs());
- while ((o = static_cast<TGraphAsymmErrors*>(next()))) {
- max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
- fDenominators.Add(o);
+ if (fShowLeftRight) {
+ fLeftRight->Add(Asymmetry(dndeta, amax));
+ fLeftRight->Add(Asymmetry(dndetaMC, amax));
+ fLeftRight->Add(Asymmetry(tracks, amax));
}
- return other;
- }
- //__________________________________________________________________
- /**
- * Make a histogram stack of ratios of results to other data
- *
- * @param max On return, the maximum value in the stack
- *
- * @return Newly allocated stack
- */
- THStack* StackRatios(TMultiGraph* others, Double_t& max)
- {
- THStack* ratios = new THStack("ratios", "Ratios");
-
- if (others) {
- TObject* ua5_1 = 0;
- TObject* ua5_2 = 0;
- TObject* alice = 0;
- TObject* cms = 0;
- TObject* denom = 0;
- TIter nextD(&fDenominators);
- while ((denom = nextD())) {
- TIter nextN(&fNumerators);
- TObject* numer = 0;
- while ((numer = nextN())) {
- ratios->Add(Ratio(numer, denom, max));
- }
- TString oName(denom->GetName());
- oName.ToLower();
- if (oName.Contains("ua5")) {
- if (ua5_1) ua5_2 = denom;
- else ua5_1 = denom;
- }
- if (oName.Contains("alice")) alice = denom;
- if (oName.Contains("cms")) cms = denom;
+ if (thisOther) {
+ TIter next(thisOther->GetListOfGraphs());
+ TGraph* g = 0;
+ while ((g = static_cast<TGraph*>(next()))) {
+ fRatios->Add(Ratio(dndeta, g, rmax));
+ fRatios->Add(Ratio(dndetaSym, g, rmax));
+ SetAttributes(g, color);
+ ModifyTitle(g, centTxt);
+ if (!fOthers->GetListOfGraphs() ||
+ !fOthers->GetListOfGraphs()->FindObject(g->GetName()))
+ fOthers->Add(g);
}
- if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max));
- if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max));
- if (cms && alice) ratios->Add(Ratio(alice, cms, max));
+ // fOthers->Add(thisOther);
}
- // Check if we have ratios
- if (!ratios->GetHists() ||
- (ratios->GetHists()->GetEntries() <= 0)) {
- delete ratios;
- ratios = 0;
+ if (tracks) {
+ if (!fRatios->GetHists()->FindObject(tracks->GetName()))
+ fRatios->Add(Ratio(dndeta, tracks, rmax));
}
- return ratios;
- }
- //__________________________________________________________________
- /**
- * Make a histogram stack of the left-right asymmetry
- *
- * @param max On return, the maximum value in the stack
- *
- * @return Newly allocated stack
- */
- THStack* StackLeftRight(Double_t& max)
- {
- THStack* ret = new THStack("leftright", "Left-right asymmetry");
- // ret->Add(Asymmetry(fForward, max));
- // ret->Add(Asymmetry(fForwardMC, max));
- if (!ret->GetHists() ||
- (ret->GetHists()->GetEntries() <= 0)) {
- delete ret;
- ret = 0;
+ if (dndetaTruth) {
+ fRatios->Add(Ratio(dndeta, dndetaTruth, rmax));
+ fRatios->Add(Ratio(dndetaSym, dndetaTruth, rmax));
+ fRatios->Add(Ratio(dndetaMC, dndetaTruth, rmax));
+ fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax));
}
- return ret;
}
//__________________________________________________________________
/**
* Plot the results
- *
- * @param results Results
- * @param others Other data
* @param max Max value
- * @param ratios Stack of ratios (optional)
* @param rmax Maximum diviation from 1 of ratios
- * @param leftright Stack of left-right asymmetry (optional)
* @param amax Maximum diviation from 1 of asymmetries
*/
- void Plot(THStack* results,
- TMultiGraph* others,
- Double_t max,
- THStack* ratios,
+ void Plot(Double_t max,
Double_t rmax,
- THStack* leftright,
Double_t amax)
{
gStyle->SetOptTitle(0);
gStyle->SetTitleFont(132, "xyz");
gStyle->SetLabelFont(132, "xyz");
- Int_t h = 800;
- Int_t w = 800; // h / TMath::Sqrt(2);
- if (!ratios) w *= 1.4;
- if (!leftright) w *= 1.4;
+ Int_t h = 800;
+ Int_t w = 800; // h / TMath::Sqrt(2);
+ Double_t y1 = 0;
+ Double_t y2 = 0;
+ Double_t y3 = 0;
+ if (!fShowRatios) w *= 1.4;
+ else y1 = 0.3;
+ if (!fShowLeftRight) w *= 1.4;
+ else {
+ Double_t y11 = y1;
+ y1 = (y11 > 0.0001 ? 0.4 : 0.2);
+ y2 = (y11 > 0.0001 ? y11 : 0.2);
+ }
TCanvas* c = new TCanvas("Results", "Results", w, h);
c->SetFillColor(0);
c->SetBorderSize(0);
c->SetBorderMode(0);
- Double_t y1 = 0;
- Double_t y2 = 0;
- Double_t y3 = 0;
- if (ratios) y1 = 0.3;
- if (leftright) {
- if (y1 > 0.0001) {
- y2 = 0.2;
- y1 = 0.4;
- }
- else {
- y1 = 0.2;
- y2 = y1;
- }
- }
- PlotResults(results, others, max, y1);
+#if 0
+ Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s",
+ y1, y2, y2, (fShowRatios ? "ratios" : ""),
+ (fShowLeftRight ? "right/left" : ""));
+#endif
+ PlotResults(max, y1);
c->cd();
- PlotRatios(ratios, rmax, y2, y1);
+ PlotRatios(rmax, y2, y1);
c->cd( );
- PlotLeftRight(leftright, amax, y3, y2);
+ PlotLeftRight(amax, y3, y2);
c->cd();
Int_t vMin = fVtxAxis->GetXmin();
Int_t vMax = fVtxAxis->GetXmax();
TString trg(fTrigString->GetTitle());
- Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
+ Int_t nev = 0;
+ if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
trg = trg.Strip(TString::kBoth);
TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
fSysString->GetTitle(),
}
//__________________________________________________________________
/**
+ * Build main legend
+ *
+ * @param x1
+ * @param y1
+ * @param x2
+ * @param y2
+ */
+ void BuildLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
+ {
+ TLegend* l = new TLegend(x1,y1,x2,y2);
+ l->SetNColumns(2);
+ l->SetFillColor(0);
+ l->SetFillStyle(0);
+ l->SetBorderSize(0);
+ l->SetTextFont(132);
+
+ TIter next(fResults->GetHists());
+ TObject* hist = 0;
+ while ((hist = next())) {
+ TString n(hist->GetTitle());
+ if (n.Contains("mirrored")) continue;
+ l->AddEntry(hist, hist->GetTitle(), "pl");
+ }
+ TIter nexto(fOthers->GetListOfGraphs());
+ while ((hist = nexto())) {
+ TString n(hist->GetTitle());
+ if (n.Contains("mirrored")) continue;
+ l->AddEntry(hist, hist->GetTitle(), "pl");
+ }
+ TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
+ d1->SetLineColor(kBlack);
+ d1->SetMarkerColor(kBlack);
+ d1->SetMarkerStyle(20);
+ TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
+ d2->SetLineColor(kBlack);
+ d2->SetMarkerColor(kBlack);
+ d2->SetMarkerStyle(24);
+
+ l->Draw();
+ }
+ //__________________________________________________________________
+ /**
* Plot the results
*
- * @param results Results
- * @param others Other data
* @param max Maximum
* @param yd Bottom position of pad
*/
- void PlotResults(THStack* results, TMultiGraph* others,
- Double_t max, Double_t yd)
+ void PlotResults(Double_t max, Double_t yd)
{
// Make a sub-pad for the result itself
TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
p1->Draw();
p1->cd();
- Info("PlotResults", "Plotting results with max=%f", max);
- results->ls();
- results->SetMaximum(1.15*max);
- results->SetMinimum(yd > 0.00001 ? -0.1 : 0);
+ // Info("PlotResults", "Plotting results with max=%f", max);
+ fResults->SetMaximum(1.15*max);
+ fResults->SetMinimum(yd > 0.00001 ? -0.1 : 0);
- FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+ FixAxis(fResults, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
p1->Clear();
- results->DrawClone("nostack e1");
+ fResults->DrawClone("nostack e1");
- fRangeParam->fSlave1Axis = results->GetXaxis();
+ fRangeParam->fSlave1Axis = fResults->GetXaxis();
fRangeParam->fSlave1Pad = p1;
// Draw other data
- if (others) {
+ if (fShowOthers || fShowAlice) {
TGraphAsymmErrors* o = 0;
- TIter nextG(others->GetListOfGraphs());
+ TIter nextG(fOthers->GetListOfGraphs());
while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
o->DrawClone("same p");
}
// Make a legend in the main result pad
+ BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
+#if 0
TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
l->SetNColumns(2);
l->SetFillColor(0);
l->SetFillStyle(0);
l->SetBorderSize(0);
l->SetTextFont(132);
+#endif
// Put a title on top
TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
tt->Draw();
// Put number of accepted events on the plot
- Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
+ Int_t nev = 0;
+ if (fTriggers) nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
et->SetNDC();
et->SetTextFont(132);
}
// results->Draw("nostack e1 same");
- fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName());
+ fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
fRangeParam->fSlave1Pad = p1;
/**
* Plot the ratios
*
- * @param ratios Ratios to plot (if any)
* @param max Maximum diviation from 1
* @param y1 Lower y coordinate of pad
* @param y2 Upper y coordinate of pad
*/
- void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2)
+ void PlotRatios(Double_t max, Double_t y1, Double_t y2)
{
- if (!ratios) return;
+ if (!fShowRatios) return;
+
bool isBottom = (y1 < 0.0001);
Double_t yd = y2 - y1;
// Make a sub-pad for the result itself
p2->cd();
// Fix up axis
- FixAxis(ratios, 1/yd/1.7, "Ratios", 7);
+ FixAxis(fRatios, 1/yd/1.7, "Ratios", 7);
- ratios->SetMaximum(1+TMath::Max(.22,1.05*max));
- ratios->SetMinimum(1-TMath::Max(.32,1.05*max));
+ fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
+ fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
p2->Clear();
- ratios->DrawClone("nostack e1");
+ fRatios->DrawClone("nostack e1");
// Make a legend
band->DrawClone("X L same");
// Replot the ratios on top
- ratios->DrawClone("nostack e1 same");
+ fRatios->DrawClone("nostack e1 same");
if (isBottom) {
- fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName());
+ fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
fRangeParam));
}
else {
- fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName());
+ fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
fRangeParam->fSlave2Pad = p2;
}
}
/**
* Plot the asymmetries
*
- * @param ratios Asymmetries to plot (if any)
* @param max Maximum diviation from 1
* @param y1 Lower y coordinate of pad
* @param y2 Upper y coordinate of pad
*/
- void PlotLeftRight(THStack* leftright, Double_t max,
- Double_t y1, Double_t y2)
+ void PlotLeftRight(Double_t max, Double_t y1, Double_t y2)
{
- if (!leftright) return;
+ if (!fShowLeftRight) return;
+
bool isBottom = (y1 < 0.0001);
Double_t yd = y2 - y1;
// Make a sub-pad for the result itself
p3->cd();
TH1* dummy = 0;
- if (leftright->GetHists()->GetEntries() == 1) {
+ if (fLeftRight->GetHists()->GetEntries() == 1) {
// Add dummy histogram
dummy = new TH1F("dummy","", 10, -6, 6);
dummy->SetLineColor(0);
dummy->SetFillColor(0);
dummy->SetMarkerColor(0);
- leftright->Add(dummy);
+ fLeftRight->Add(dummy);
}
// Fix up axis
- FixAxis(leftright, 1/yd/1.7, "Right/Left", 4);
+ FixAxis(fLeftRight, 1/yd/1.7, "Right/Left", 4);
- leftright->SetMaximum(1+TMath::Max(.12,1.05*max));
- leftright->SetMinimum(1-TMath::Max(.15,1.05*max));
+ fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
+ fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
p3->Clear();
- leftright->DrawClone("nostack e1");
+ fLeftRight->DrawClone("nostack e1");
// Make a legend
band->Draw("3 same");
band->DrawClone("X L same");
- leftright->DrawClone("nostack e1 same");
+ fLeftRight->DrawClone("nostack e1 same");
if (isBottom) {
- fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName());
+ fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
fRangeParam));
}
else {
- fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName());
+ fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
fRangeParam->fSlave2Pad = p3;
}
}
{
if (!list) return 0;
- TList* all = static_cast<TList*>(list->FindObject("all"));
- if (!all) {
- Warning("GetResult", "No 'all' list find in %s", list->GetName());
- // list->ls();
- return 0;
- }
-
- TH1* ret = static_cast<TH1*>(all->FindObject(name));
+ TH1* ret = static_cast<TH1*>(list->FindObject(name));
if (!ret) {
// all->ls();
Warning("GetResult", "Histogram %s not found", name);
*
* @return Maximum of histogram
*/
- Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const
+ Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const
{
// Check if we have input
if (!hist) return 0;
// Rebin if needed
Rebin(hist);
+ // Info("AddHistogram", "Adding %s to %s",
+ // hist->GetName(), stack->GetName());
stack->Add(hist, option);
+ // stack->ls();
return hist->GetMaximum();
}
//__________________________________________________________________
*
* @return Maximum of histogram
*/
- Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option,
- TH1*& sym) const
+ Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
+ Option_t* option="") const
{
// Check if we have input
if (!hist) return 0;
sym = Symmetrice(hist);
stack->Add(sym, option);
+ // Info("AddHistogram", "Adding %s and %s to %s",
+ // hist->GetName(), sym->GetName(), stack->GetName());
return hist->GetMaximum();
}
// Double_t dBin = (high - low) / oBins;
// Int_t tBins = Int_t(2*high/dBin+.5);
// ret->SetBins(tBins, -high, high);
+ ret->SetDirectory(0);
ret->Reset();
ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
ret->SetYTitle("Right/Left");
ret->SetMarkerStyle(g->GetMarkerStyle());
ret->SetMarkerColor(h->GetMarkerColor());
ret->SetMarkerSize(0.9*g->GetMarkerSize());
+ ret->SetDirectory(0);
Double_t xlow = g->GetX()[0];
Double_t xhigh = g->GetX()[g->GetN()-1];
if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
t1->Divide(h2);
t1->SetMarkerColor(h1->GetMarkerColor());
t1->SetMarkerStyle(h2->GetMarkerStyle());
+ t1->SetDirectory(0);
max = TMath::Max(RatioMax(t1), max);
return t1;
}
h->SetMarkerStyle(g2->GetMarkerStyle());
h->SetMarkerColor(g1->GetMarkerColor());
h->SetMarkerSize(0.9*g2->GetMarkerSize());
+ h->SetDirectory(0);
Double_t low = g2->GetX()[0];
Double_t high = g2->GetX()[g2->GetN()-1];
void FixAxis(THStack* stack, Double_t s, const char* ytitle,
Int_t ynDiv=210, Bool_t force=true)
{
+ if (!stack) {
+ Warning("FixAxis", "No stack passed for %s!", ytitle);
+ return;
+ }
if (force) stack->Draw("nostack e1");
TH1* h = stack->GetHistogram();
- if (!h) return;
-
+ if (!h) {
+ Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
+ return;
+ }
+
h->SetXTitle("#eta");
h->SetYTitle(ytitle);
TAxis* xa = h->GetXaxis();
xa->SetTitleSize(s*xa->GetTitleSize());
xa->SetLabelSize(s*xa->GetLabelSize());
xa->SetTickLength(s*xa->GetTickLength());
+
+ if (stack != fResults) {
+ TAxis* rxa = fResults->GetXaxis();
+ xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
+ }
}
if (ya) {
ya->SetTitle(ytitle);
//__________________________________________________________________
- Bool_t fShowOthers; // Show other data
- Bool_t fShowAlice; // Show ALICE published data
- Bool_t fShowRatios; // Show ratios
- Bool_t fShowLeftRight;// Show asymmetry
- UShort_t fRebin; // Rebinning factor
- Bool_t fCutEdges; // Whether to cut edges
- TString fTitle; // Title on plot
- TNamed* fTrigString; // Trigger string (read, or set)
- TNamed* fSNNString; // Energy string (read, or set)
- TNamed* fSysString; // Collision system string (read or set)
- TAxis* fVtxAxis; // Vertex cuts (read or set)
- TList fNumerators; // Numerators in ratios
- TList fDenominators; // Denominators in ratios
- THStack* fResults; // Stack of results
- THStack* fRatios; // Stack of ratios
- // TH1* fForward; // Results
- // TH1* fForwardMC; // MC results
- // TH1* fForwardHHD; // Old results
- // TH1* fTruth; // MC truth
- // TH1* fCentral; // Central region data
- // TH1* fForwardSym; // Symmetric extension
- // TH1* fForwardMCSym; // Symmetric extension
- // TH1* fForwardHHDSym;// Symmetric extension
- TH1* fTriggers; // Number of triggers
- RangeParam* fRangeParam; // Parameter object for range zoom
+ Bool_t fShowOthers; // Show other data
+ Bool_t fShowAlice; // Show ALICE published data
+ Bool_t fShowRatios; // Show ratios
+ Bool_t fShowLeftRight;// Show asymmetry
+ UShort_t fRebin; // Rebinning factor
+ Bool_t fCutEdges; // Whether to cut edges
+ TString fTitle; // Title on plot
+ TNamed* fTrigString; // Trigger string (read, or set)
+ TNamed* fSNNString; // Energy string (read, or set)
+ TNamed* fSysString; // Collision system string (read or set)
+ TAxis* fVtxAxis; // Vertex cuts (read or set)
+ TAxis* fCentAxis; // Centrality axis
+ THStack* fResults; // Stack of results
+ THStack* fRatios; // Stack of ratios
+ THStack* fLeftRight; // Left-right asymmetry
+ TMultiGraph* fOthers; // Older data
+ TH1* fTriggers; // Number of triggers
+ RangeParam* fRangeParam; // Parameter object for range zoom
};
// --- AOD output handler ------------------------------------------
AliAODHandler* aodHandler = new AliAODHandler();
mgr->SetOutputEventHandler(aodHandler);
+ aodHandler->SetNeedsHeaderReplication();
aodHandler->SetOutputFileName("AliAOD.root");
// --- Add tasks ---------------------------------------------------
mgr->SetSkipTerminate(false);
// Some informative output
mgr->PrintStatus();
- // mgr->SetDebugLevel(3);
+ if (proof) mgr->SetDebugLevel(3);
if (mgr->GetDebugLevel() < 1 && !proof)
mgr->SetUseProgressBar(kTRUE,100);
*
* @ingroup pwg2_forward_scripts
*/
-void MakedNdeta(const char* aoddir=".",
- Int_t nEvents=-1,
- const char* trig="INEL",
- Float_t centLow = 0,
- Float_t centHigh = 100,
- Double_t vzMin=-10,
- Double_t vzMax=10,
- Int_t proof=0)
+void MakedNdeta(const char* aoddir = ".",
+ Int_t nEvents = -1,
+ const char* trig = "INEL",
+ Bool_t useCent = false,
+ Double_t vzMin = -10,
+ Double_t vzMax = +10,
+ Int_t proof = 0)
{
// --- Libraries to load -------------------------------------------
gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
// --- Add tasks ---------------------------------------------------
// Forward
gROOT->LoadMacro("AddTaskForwarddNdeta.C");
- AddTaskForwarddNdeta(trig, vzMin, vzMax, centLow, centHigh, true);
+ AddTaskForwarddNdeta(trig, vzMin, vzMax, useCent, true);
// Central
gROOT->LoadMacro("AddTaskCentraldNdeta.C");
- AddTaskCentraldNdeta(trig, vzMin, vzMax, centLow, centHigh);
+ AddTaskCentraldNdeta(trig, vzMin, vzMax, useCent);
// --- Run the analysis --------------------------------------------
/**
* Get a multi graph of data for a given energy and trigger type
*
+ * @param sys Collision system (1: pp, 2: PbPb)
* @param energy Energy in GeV (900, 2360, 7000)
* @param type Bit pattern of trigger type
* - 0x1 INEL
* @ingroup pwg2_forward_otherdata
*/
TMultiGraph*
-GetData(Int_t energy, Int_t type, bool aliceOnly=false)
+GetData(UShort_t sys,
+ UShort_t energy,
+ UShort_t type=0x1,
+ UShort_t centLow=0,
+ UShort_t centHigh=0,
+ bool aliceOnly=false)
{
- TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d", energy, type),"");
+ TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d_%03d_%03d",
+ energy, type, centLow, centHigh),"");
TString tn;
TString en;
- if (TMath::Abs(energy-900) < 10) {
- en.Append(" #sqrt{s}=900GeV");
- if (type & 0x1) {
- tn.Append(" INEL");
- if (!aliceOnly) mp->Add(UA5Inel(false));
- if (!aliceOnly) mp->Add(UA5Inel(true));
- mp->Add(AliceCentralInel900());
- }
- if (type & 0x4) {
- tn.Append(" NSD");
- if (!aliceOnly) mp->Add(UA5Nsd(false));
- if (!aliceOnly) mp->Add(UA5Nsd(true));
- mp->Add(AliceCentralNsd900());
- if (!aliceOnly) mp->Add(CMSNsd900());
+ TString sn;
+ TString cn;
+ if (sys == 1) {
+ sn = ", pp(p#bar{p})";
+ if (energy < 1000)
+ en = Form(", #sqrt{s}=%dGeV", energy);
+ else
+ en = Form(", #sqrt{s}=%f4.2TeV", float(energy)/1000);
+ if (!(type & 0x7))
+ Warning("GetData", "Unknown trigger mask 0x%x", type);
+
+ if (TMath::Abs(energy-900) < 10) {
+ if (type & 0x1) {
+ tn.Append(" INEL");
+ if (!aliceOnly) mp->Add(UA5Inel(false));
+ if (!aliceOnly) mp->Add(UA5Inel(true));
+ mp->Add(AliceCentralInel900());
+ }
+ if (type & 0x4) {
+ tn.Append(" NSD");
+ if (!aliceOnly) mp->Add(UA5Nsd(false));
+ if (!aliceOnly) mp->Add(UA5Nsd(true));
+ mp->Add(AliceCentralNsd900());
+ if (!aliceOnly) mp->Add(CMSNsd900());
+ }
+ if (type & 0x2) {
+ tn.Append(" INEL>0");
+ mp->Add(AliceCentralInelGt900());
+ }
}
- if (type & 0x2) {
- tn.Append(" INEL>0");
- mp->Add(AliceCentralInelGt900());
+ else if (TMath::Abs(energy-2360) < 10) {
+ if (type & 0x1) {
+ tn.Append(" INEL");
+ mp->Add(AliceCentralInel2360());
+ }
+ if (type & 0x4) {
+ tn.Append(" NSD");
+ mp->Add(AliceCentralNsd2360());
+ if (!aliceOnly) mp->Add(CMSNsd2360());
+ }
+ if (type & 0x2) {
+ tn.Append(" INEL>0");
+ mp->Add(AliceCentralInelGt2360());
+ }
}
- }
- if (TMath::Abs(energy-2360) < 10) {
- en.Append(" #sqrt{s}=2.36TeV");
- if (type & 0x1) {
- tn.Append(" INEL");
- mp->Add(AliceCentralInel2360());
- }
- if (type & 0x4) {
- tn.Append(" NSD");
- mp->Add(AliceCentralNsd2360());
- if (!aliceOnly) mp->Add(CMSNsd2360());
- }
- if (type & 0x1) {
- tn.Append(" INEL>0");
- mp->Add(AliceCentralInelGt2360());
+ else if (TMath::Abs(energy-7000) < 10) {
+ if (type & 0x1) {
+ tn.Append(" INEL");
+ }
+ if (type & 0x4) {
+ tn.Append(" NSD");
+ if (!aliceOnly) mp->Add(CMSNsd7000());
+ }
+ if (type & 0x2) {
+ tn.Append(" INEL>0");
+ mp->Add(AliceCentralInelGt7000());
+ }
}
+ else
+ Warning("GetData", "No other results for sys=%d, energy=%d",
+ sys, energy);
}
- if (TMath::Abs(energy-7000) < 10) {
- en.Append(" #sqrt{s}=7TeV");
- if (type & 0x1) {
- tn.Append(" INEL");
- }
- if (type & 0x4) {
- tn.Append(" NSD");
- if (!aliceOnly) mp->Add(CMSNsd7000());
- }
- if (type & 0x1) {
- tn.Append(" INEL>0");
- mp->Add(AliceCentralInelGt7000());
- }
+ else if (sys == 2) {
+ // Nothing for PbPb so far
+ cn = Form(", %d%%-%d%% central", centLow, centHigh);
+ sn = ", PbPb";
+ if (energy < 1000)
+ en = Form(", #sqrt{s_{NN}}=%dGeV", energy);
+ else
+ en = Form(", #sqrt{s_{NN}}=%f4.2TeV", float(energy)/1000);
+ Warning("GetData", "No other data for PbP b yet");
}
- mp->SetTitle(Form("1/N dN_{ch}/d#eta, pp(p#bar{p}), %s, %s",
- en.Data(), tn.Data()));
+ else
+ Warning("GetData", "Unknown system %d", sys);
+ TString tit(Form("1/N dN_{ch}/d#eta%s%s%s%s",
+ sn.Data(), en.Data(), tn.Data(), cn.Data()));
+ mp->SetTitle(tit.Data());
if (!mp->GetListOfGraphs() || mp->GetListOfGraphs()->GetEntries() <= 0) {
delete mp;
mp = 0;
* @ingroup pwg2_forward_otherdata
*/
void
-OtherData(Int_t energy=900, Int_t type=0x1, bool aliceOnly=false)
+OtherData(UShort_t sys=1,
+ UShort_t energy=900,
+ UShort_t type=0x1,
+ UShort_t centLow=0,
+ UShort_t centHigh=5,
+ bool aliceOnly=false)
{
- TMultiGraph* mp = GetData(energy, type, aliceOnly);
+ TMultiGraph* mp = GetData(sys, energy, type, centLow, centHigh, aliceOnly);
if (!mp) return;
gStyle->SetTitleX(0.1);
rebin=1
vzmin=-10
vzmax=10
-centlow=0
-centhigh=100
batch=0
gdb=0
proof=0
if test $dopass2 -gt 0 ; then
rotate ${output2}
- args=(\(\".\",$nev,\"$type\",$centlow,$centhigh,$vzmin,$vzmax,$proof\))
+ args=(\(\".\",$nev,\"$type\",$cent,$vzmin,$vzmax,$proof\))
if test "x$pass1" = "xMakeELossFits.C" ; then
args=(\(\"${output1}\"\))
fi
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff