#include "AliUEHist.h"
#include "AliCFContainer.h"
#include "THnSparse.h"
#include "TMath.h"
#include "TList.h"
#include "TCollection.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "AliLog.h"
#include "TCanvas.h"
#include "TF1.h"
#include "AliTHn.h"
#include "THn.h"
ClassImp(AliUEHist)
const Int_t AliUEHist::fgkCFSteps = 11;
AliUEHist::AliUEHist(const char* reqHist, const char* binning) :
TObject(),
fkRegions(4),
fEventHist(0),
fTrackHistEfficiency(0),
fFakePt(0),
fEtaMin(0),
fEtaMax(0),
fPtMin(0),
fPtMax(0),
fPartSpecies(-1),
fCentralityMin(0),
fCentralityMax(0),
fZVtxMin(0),
fZVtxMax(0),
fPt2Min(0),
fContaminationEnhancement(0),
fCombineMinMax(0),
fTrackEtaCut(0),
fWeightPerEvent(0),
fSkipScaleMixedEvent(kFALSE),
fCache(0),
fGetMultCacheOn(kFALSE),
fGetMultCache(0),
fHistogramType(reqHist)
{
for (UInt_t i=0; i<fkRegions; i++)
fTrackHist[i] = 0;
if (strlen(reqHist) == 0)
return;
Printf("Creating AliUEHist with %s (binning: %s)", reqHist, binning);
AliLog::SetClassDebugLevel("AliCFContainer", -1);
AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
const char* title = "";
Int_t nTrackVars = 4;
Int_t iTrackBin[7];
Double_t* trackBins[7];
const char* trackAxisTitle[7];
Int_t nEtaBins = -1;
Double_t* etaBins = GetBinning(binning, "eta", nEtaBins);
const char* etaTitle = "#eta";
iTrackBin[0] = nEtaBins;
trackBins[0] = etaBins;
trackAxisTitle[0] = etaTitle;
Int_t nDeltaEtaBins = -1;
Double_t* deltaEtaBins = GetBinning(binning, "delta_eta", nDeltaEtaBins);
trackBins[1] = GetBinning(binning, "p_t_assoc", iTrackBin[1]);
trackAxisTitle[1] = "p_{T} (GeV/c)";
Int_t npTBinsFine = -1;
Double_t* pTBinsFine = GetBinning(binning, "p_t_eff", npTBinsFine);
Int_t nLeadingpTBins = -1;
Double_t* leadingpTBins = GetBinning(binning, "p_t_leading", nLeadingpTBins);
Int_t nLeadingpTBins2 = -1;
Double_t* leadingpTBins2 = GetBinning(binning, "p_t_leading_course", nLeadingpTBins2);
Int_t nLeadingPhiBins = -1;
Double_t* leadingPhiBins = GetBinning(binning, "delta_phi", nLeadingPhiBins);
trackBins[3] = GetBinning(binning, "multiplicity", iTrackBin[3]);
trackAxisTitle[3] = "multiplicity";
const Int_t kNSpeciesBins = 4;
Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
const char* vertexTitle = "z-vtx (cm)";
Int_t nVertexBins = -1;
Double_t* vertexBins = GetBinning(binning, "vertex", nVertexBins);
Int_t nVertexBinsEff = -1;
Double_t* vertexBinsEff = GetBinning(binning, "vertex_eff", nVertexBinsEff);
Int_t useVtxAxis = 0;
Int_t axis = -1;
if (strcmp(reqHist, "NumberDensitypT") == 0)
{
axis = 0;
title = "d^{2}N_{ch}/d#varphid#eta";
}
else if (strcmp(reqHist, "NumberDensityPhi") == 0)
{
axis = 1;
title = "d^{2}N_{ch}/d#varphid#eta";
}
else if (TString(reqHist).BeginsWith("NumberDensityPhiCentrality"))
{
if (TString(reqHist).Contains("Vtx"))
useVtxAxis = 1;
reqHist = "NumberDensityPhiCentrality";
fHistogramType = reqHist;
axis = 2;
title = "d^{2}N_{ch}/d#varphid#eta";
}
else if (strcmp(reqHist, "SumpT") == 0)
{
axis = 0;
title = "d^{2}#Sigma p_{T}/d#varphid#eta";
}
else if (TString(reqHist).BeginsWith("TwoPlusOne"))
{
useVtxAxis = 1;
reqHist = "TwoPlusOne";
fHistogramType = reqHist;
axis = 3;
title = "d^{2}N_{ch}/d#varphid#eta";
}
else
AliFatal(Form("Invalid histogram requested: %s", reqHist));
UInt_t initRegions = fkRegions;
if (axis == 0)
{
trackBins[2] = leadingpTBins;
iTrackBin[2] = nLeadingpTBins;
trackAxisTitle[2] = "leading p_{T} (GeV/c)";
}
else if (axis == 1)
{
nTrackVars = 5;
initRegions = 1;
iTrackBin[2] = nLeadingpTBins2;
trackBins[2] = leadingpTBins2;
trackAxisTitle[2] = "leading p_{T} (GeV/c)";
iTrackBin[4] = nLeadingPhiBins;
trackBins[4] = leadingPhiBins;
trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
}
else if (axis == 2)
{
nTrackVars = 5;
initRegions = 1;
iTrackBin[0] = nDeltaEtaBins;
trackBins[0] = deltaEtaBins;
trackAxisTitle[0] = "#Delta#eta";
iTrackBin[2] = nLeadingpTBins2;
trackBins[2] = leadingpTBins2;
trackAxisTitle[2] = "leading p_{T} (GeV/c)";
trackAxisTitle[3] = "centrality";
iTrackBin[4] = nLeadingPhiBins;
trackBins[4] = leadingPhiBins;
trackAxisTitle[4] = "#Delta#varphi (rad)";
if (useVtxAxis > 0)
{
nTrackVars = 6;
iTrackBin[5] = nVertexBins;
trackBins[5] = vertexBins;
trackAxisTitle[5] = vertexTitle;
}
}
else if (axis == 3)
{
nTrackVars = 7;
initRegions = 1;
iTrackBin[0] = nDeltaEtaBins;
trackBins[0] = deltaEtaBins;
trackAxisTitle[0] = "#Delta#eta";
iTrackBin[2] = nLeadingpTBins;
trackBins[2] = leadingpTBins;
trackAxisTitle[2] = "Trigger 1 p_{T} (GeV/c)";
trackAxisTitle[3] = "centrality";
iTrackBin[4] = nLeadingPhiBins;
trackBins[4] = leadingPhiBins;
trackAxisTitle[4] = "#Delta#varphi (rad)";
iTrackBin[5] = nVertexBins;
trackBins[5] = vertexBins;
trackAxisTitle[5] = vertexTitle;
iTrackBin[6] = nLeadingpTBins2;
trackBins[6] = leadingpTBins2;
trackAxisTitle[6] = "Trigger 2 p_{T} (GeV/c)";
}
for (UInt_t i=0; i<initRegions; i++)
{
if (axis >= 2)
fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
else
fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
for (Int_t j=0; j<nTrackVars; j++)
{
fTrackHist[i]->SetBinLimits(j, trackBins[j]);
fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
}
SetStepNames(fTrackHist[i]);
}
Int_t nEventVars = 2;
Int_t iEventBin[4];
iEventBin[0] = iTrackBin[2];
iEventBin[1] = iTrackBin[3];
if (axis >= 2 && useVtxAxis)
{
nEventVars = 3;
iEventBin[2] = iTrackBin[5];
}
if (axis >= 3)
{
nEventVars = 4;
iEventBin[3] = iTrackBin[6];
}
fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
fEventHist->SetBinLimits(0, trackBins[2]);
fEventHist->SetVarTitle(0, trackAxisTitle[2]);
fEventHist->SetBinLimits(1, trackBins[3]);
fEventHist->SetVarTitle(1, trackAxisTitle[3]);
if (axis >= 2 && useVtxAxis)
{
fEventHist->SetBinLimits(2, trackBins[5]);
fEventHist->SetVarTitle(2, trackAxisTitle[5]);
}
if (axis >= 3)
{
fEventHist->SetBinLimits(3, trackBins[6]);
fEventHist->SetVarTitle(3, trackAxisTitle[6]);
}
SetStepNames(fEventHist);
iTrackBin[0] = nEtaBins;
iTrackBin[1] = npTBinsFine;
iTrackBin[2] = kNSpeciesBins;
iTrackBin[4] = nVertexBinsEff;
fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 6, 5, iTrackBin);
fTrackHistEfficiency->SetBinLimits(0, etaBins);
fTrackHistEfficiency->SetVarTitle(0, etaTitle);
fTrackHistEfficiency->SetBinLimits(1, pTBinsFine);
fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
fTrackHistEfficiency->SetBinLimits(2, speciesBins);
fTrackHistEfficiency->SetVarTitle(2, "particle species");
fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
fTrackHistEfficiency->SetBinLimits(4, vertexBinsEff);
fTrackHistEfficiency->SetVarTitle(4, vertexTitle);
fFakePt = new TH3F("fFakePt","fFakePt;p_{T,rec};p_{T};centrality", 200, 0, 20, 200, 0, 20, 20, 0, 100);
delete[] deltaEtaBins;
delete[] pTBinsFine;
delete[] leadingpTBins;
delete[] leadingpTBins2;
delete[] leadingPhiBins;
delete[] vertexBins;
delete[] vertexBinsEff;
}
TString AliUEHist::CombineBinning(TString defaultBinning, TString customBinning)
{
TString binningStr;
TObjArray* lines = defaultBinning.Tokenize("\n");
for (Int_t i=0; i<lines->GetEntriesFast(); i++)
{
TString line(lines->At(i)->GetName());
TString tag = line(0, line.Index(":")+1);
if (!customBinning.BeginsWith(tag) && !customBinning.Contains(TString("\n") + tag))
binningStr += line + "\n";
else
Printf("Using custom binning for %s", tag.Data());
}
delete lines;
binningStr += customBinning;
return binningStr;
}
Double_t* AliUEHist::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
{
TString config(configuration);
TObjArray* lines = config.Tokenize("\n");
for (Int_t i=0; i<lines->GetEntriesFast(); i++)
{
TString line(lines->At(i)->GetName());
if (line.BeginsWith(TString(tag) + ":"))
{
line.Remove(0, strlen(tag) + 1);
line.ReplaceAll(" ", "");
TObjArray* binning = line.Tokenize(",");
Double_t* bins = new Double_t[binning->GetEntriesFast()];
for (Int_t j=0; j<binning->GetEntriesFast(); j++)
bins[j] = TString(binning->At(j)->GetName()).Atof();
nBins = binning->GetEntriesFast() - 1;
delete binning;
delete lines;
return bins;
}
}
delete lines;
AliFatal(Form("Tag %s not found in %s", tag, configuration));
return 0;
}
AliUEHist::AliUEHist(const AliUEHist &c) :
TObject(),
fkRegions(4),
fEventHist(0),
fTrackHistEfficiency(0),
fFakePt(0),
fEtaMin(0),
fEtaMax(0),
fPtMin(0),
fPtMax(0),
fPartSpecies(-1),
fCentralityMin(0),
fCentralityMax(0),
fZVtxMin(0),
fZVtxMax(0),
fPt2Min(0),
fContaminationEnhancement(0),
fCombineMinMax(0),
fTrackEtaCut(0),
fWeightPerEvent(0),
fSkipScaleMixedEvent(kFALSE),
fCache(0),
fGetMultCacheOn(kFALSE),
fGetMultCache(0),
fHistogramType()
{
((AliUEHist &) c).Copy(*this);
}
void AliUEHist::SetStepNames(AliCFContainer* container)
{
for (Int_t i=0; i<fgkCFSteps; i++)
container->SetStepTitle(i, GetStepTitle((CFStep) i));
}
AliUEHist::~AliUEHist()
{
for (UInt_t i=0; i<fkRegions; i++)
{
if (fTrackHist[i])
{
delete fTrackHist[i];
fTrackHist[i] = 0;
}
}
if (fEventHist)
{
delete fEventHist;
fEventHist = 0;
}
if (fTrackHistEfficiency)
{
delete fTrackHistEfficiency;
fTrackHistEfficiency = 0;
}
if (fFakePt)
{
delete fFakePt;
fFakePt = 0;
}
if (fCache)
{
delete fCache;
fCache = 0;
}
}
AliUEHist &AliUEHist::operator=(const AliUEHist &c)
{
if (this != &c)
((AliUEHist &) c).Copy(*this);
return *this;
}
void AliUEHist::Copy(TObject& c) const
{
AliUEHist& target = (AliUEHist &) c;
for (UInt_t i=0; i<fkRegions; i++)
if (fTrackHist[i])
target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
if (fEventHist)
target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
if (fTrackHistEfficiency)
target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
if (fFakePt)
target.fFakePt = dynamic_cast<TH3F*> (fFakePt->Clone());
target.fEtaMin = fEtaMin;
target.fEtaMax = fEtaMax;
target.fPtMin = fPtMin;
target.fPtMax = fPtMax;
target.fPartSpecies = fPartSpecies;
target.fCentralityMin = fCentralityMin;
target.fCentralityMax = fCentralityMax;
target.fZVtxMin = fZVtxMin;
target.fZVtxMax = fZVtxMax;
target.fPt2Min = fPt2Min;
if (fContaminationEnhancement)
target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
target.fCombineMinMax = fCombineMinMax;
target.fTrackEtaCut = fTrackEtaCut;
target.fWeightPerEvent = fWeightPerEvent;
target.fSkipScaleMixedEvent = fSkipScaleMixedEvent;
target.fHistogramType = fHistogramType;
}
Long64_t AliUEHist::Merge(TCollection* list)
{
if (!list)
return 0;
if (list->IsEmpty())
return 1;
TIterator* iter = list->MakeIterator();
TObject* obj;
const UInt_t kMaxLists = fkRegions+3;
TList** lists = new TList*[kMaxLists];
for (UInt_t i=0; i<kMaxLists; i++)
lists[i] = new TList;
Int_t count = 0;
while ((obj = iter->Next())) {
AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
if (entry == 0)
continue;
for (UInt_t i=0; i<fkRegions; i++)
if (entry->fTrackHist[i])
lists[i]->Add(entry->fTrackHist[i]);
lists[fkRegions]->Add(entry->fEventHist);
lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
if (entry->fFakePt)
lists[fkRegions+2]->Add(entry->fFakePt);
count++;
}
for (UInt_t i=0; i<fkRegions; i++)
if (fTrackHist[i])
fTrackHist[i]->Merge(lists[i]);
fEventHist->Merge(lists[fkRegions]);
fTrackHistEfficiency->Merge(lists[fkRegions+1]);
if (fFakePt)
fFakePt->Merge(lists[fkRegions+2]);
for (UInt_t i=0; i<kMaxLists; i++)
delete lists[i];
delete[] lists;
return count+1;
}
void AliUEHist::SetBinLimits(THnBase* grid)
{
if (fEtaMax > fEtaMin)
grid->GetAxis(0)->SetRangeUser(fEtaMin, fEtaMax);
if (fPtMax > fPtMin)
grid->GetAxis(1)->SetRangeUser(fPtMin, fPtMax);
if (fPt2Min > 0)
grid->GetAxis(6)->SetRangeUser(fPt2Min, grid->GetAxis(6)->GetXmax() - 0.01);
}
void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
{
SetBinLimits(grid->GetGrid());
}
void AliUEHist::ResetBinLimits(THnBase* grid)
{
for (Int_t i=0; i<grid->GetNdimensions(); i++)
if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange))
grid->GetAxis(i)->SetRangeUser(0, -1);
}
void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
{
ResetBinLimits(grid->GetGrid());
}
void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
{
Int_t binBegin[4];
Int_t binEnd[4];
for (Int_t i=0; i<4; i++)
{
binBegin[i] = 1;
binEnd[i] = fTrackHist[0]->GetNBins(i);
}
if (fEtaMax > fEtaMin)
{
binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
binEnd[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
}
if (fPtMax > fPtMin)
{
binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
binEnd[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
}
if (ptLeadMax > ptLeadMin)
{
binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
binEnd[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
}
binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
for (UInt_t region=0; region<fkRegions; region++)
{
Int_t total = 0;
Int_t count = 0;
Int_t vars[4];
for (Int_t i=0; i<4; i++)
vars[i] = binBegin[i];
AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
while (1)
{
if (grid->GetElement(vars) == 0)
{
Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
grid->GetBinCenter(0, vars[0]),
grid->GetBinCenter(1, vars[1]),
grid->GetBinCenter(2, vars[2]),
grid->GetBinCenter(3, vars[3])
);
count++;
}
vars[3]++;
for (Int_t i=3; i>0; i--)
{
if (vars[i] == binEnd[i]+1)
{
vars[i] = binBegin[i];
vars[i-1]++;
}
}
if (vars[0] == binEnd[0]+1)
break;
total++;
}
Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
}
}
TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm, Long64_t* normEvents)
{
ResetBinLimits(fTrackHist[region]->GetGrid(step));
ResetBinLimits(fEventHist->GetGrid(step));
SetBinLimits(fTrackHist[region]->GetGrid(step));
if (fZVtxMax > fZVtxMin)
{
Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
}
TH1D* tracks = 0;
if (ptLeadMin < 0)
{
tracks = fTrackHist[region]->ShowProjection(2, step);
tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
if (fCombineMinMax && region == kMin)
{
ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
SetBinLimits(fTrackHist[kMax]->GetGrid(step));
TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
tracks->Add(tracks2);
ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
}
TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
Float_t phiRegion = TMath::TwoPi() / 3;
if (!fCombineMinMax && region == kMin)
phiRegion /= 2;
Float_t normalization = phiRegion;
if (etaNorm)
normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
tracks->Scale(1.0 / normalization);
TH1D* events = fEventHist->ShowProjection(0, step);
tracks->Divide(events);
}
else
{
if (multBinEnd >= multBinBegin)
{
Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
}
Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
Printf("Using leading pT range %d --> %d", firstBin, lastBin);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
if (twoD == 0)
tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
else
tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
Printf("Calculated histogram --> %f tracks", tracks->Integral());
fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
if (twoD == 10 || twoD == 11)
{
Float_t etaLimit = 1.0;
if (twoD == 10)
{
tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
tracks->Scale(1. / 0.75);
}
else if (twoD == 11)
{
TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
tracksTmp1->Add(tracksTmp2);
delete tracks;
tracks = tracksTmp1;
delete tracksTmp2;
tracks->Scale(1. / 0.25);
}
}
Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
if (etaNorm)
{
TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
if (strcmp(axis->GetTitle(), "#eta") == 0)
{
Printf("Normalizing using eta-axis range");
normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
}
else
{
Printf("Normalizing assuming accepted range of |eta| < 0.8");
normalization *= 0.8 * 2;
}
}
tracks->Scale(1.0 / normalization);
TH1D* events = fEventHist->ShowProjection(0, step);
Long64_t nEvents = (Long64_t) events->Integral(firstBin, lastBin);
Printf("Calculated histogram --> %lld events", nEvents);
if (normEvents)
*normEvents = nEvents;
if (nEvents > 0)
tracks->Scale(1.0 / nEvents);
delete events;
}
ResetBinLimits(fTrackHist[region]->GetGrid(step));
ResetBinLimits(fEventHist->GetGrid(step));
return tracks;
}
void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist)
{
if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
AliFatal("Histogram without vertex axis provided");
ResetBinLimits(fTrackHist[region]->GetGrid(step));
ResetBinLimits(fEventHist->GetGrid(step));
SetBinLimits(fTrackHist[region]->GetGrid(step));
if (multBinEnd >= multBinBegin)
{
Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
}
Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
Printf("Using leading pT range %d --> %d", firstBin, lastBin);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
*trackHist = (TH3*) fTrackHist[region]->GetGrid(step)->Project(4, 0, 5);
*eventHist = (TH1*) fEventHist->GetGrid(step)->Project(2);
ResetBinLimits(fTrackHist[region]->GetGrid(step));
ResetBinLimits(fEventHist->GetGrid(step));
}
void AliUEHist::GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist)
{
if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
AliFatal("Histogram without vertex axis provided");
THnBase* sparse = fTrackHist[region]->GetGrid(step)->GetGrid();
if (fGetMultCacheOn)
{
if (!fGetMultCache)
{
fGetMultCache = ChangeToThn(sparse);
}
sparse = fGetMultCache;
}
ResetBinLimits(sparse);
ResetBinLimits(fEventHist->GetGrid(step));
SetBinLimits(sparse);
Int_t firstBin = sparse->GetAxis(2)->FindBin(ptLeadMin);
Int_t lastBin = sparse->GetAxis(2)->FindBin(ptLeadMax);
Printf("Using leading pT range %d --> %d", firstBin, lastBin);
sparse->GetAxis(2)->SetRange(firstBin, lastBin);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
if (fPt2Min > 0)
{
Int_t firstBinPt2 = sparse->GetAxis(6)->FindBin(fPt2Min);
Int_t lastBinPt2 = sparse->GetAxis(6)->GetNbins();
fEventHist->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(firstBinPt2, lastBinPt2);
}
Int_t dimensions[] = { 4, 0, 5, 3 };
THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E");
*eventHist = (TH2*) fEventHist->GetGrid(step)->Project(2, 1);
ResetBinLimits(sparse);
ResetBinLimits(fEventHist->GetGrid(step));
*trackHist = ChangeToThn(tmpTrackHist);
delete tmpTrackHist;
}
TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t normalizePerTrigger, Int_t stepForMixed, Int_t* trigger)
{
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
TH2* totalTracks = 0;
THnBase* trackSameAll = 0;
THnBase* trackMixedAll = 0;
THnBase* trackMixedAllStep6 = 0;
TH2* eventSameAll = 0;
TH2* eventMixedAll = 0;
TH2* eventMixedAllStep6 = 0;
Long64_t totalEvents = 0;
Int_t nCorrelationFunctions = 0;
GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackSameAll, &eventSameAll);
mixed->GetHistsZVtxMult((stepForMixed == -1) ? step : (CFStep) stepForMixed, region, ptLeadMin, ptLeadMax, &trackMixedAll, &eventMixedAll);
if (stepForMixed == -1 && step == kCFStepBiasStudy && mixed->fEventHist->GetGrid(kCFStepReconstructed)->GetEntries() > 0 && !fSkipScaleMixedEvent)
{
Printf("Using mixed-event normalization factors from step %d", kCFStepReconstructed);
mixed->GetHistsZVtxMult(kCFStepReconstructed, region, ptLeadMin, ptLeadMax, &trackMixedAllStep6, &eventMixedAllStep6);
}
TAxis* multAxis = trackSameAll->GetAxis(3);
if (multBinEnd < multBinBegin)
{
multBinBegin = 1;
multBinEnd = multAxis->GetNbins();
}
for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++)
{
trackSameAll->GetAxis(3)->SetRange(multBin, multBin);
trackMixedAll->GetAxis(3)->SetRange(multBin, multBin);
if (trackMixedAllStep6)
trackMixedAllStep6->GetAxis(3)->SetRange(multBin, multBin);
Double_t mixedNorm = 1;
Double_t mixedNormError = 0;
if (!fSkipScaleMixedEvent)
{
TH2* tracksMixed = 0;
if (trackMixedAllStep6)
{
trackMixedAllStep6->GetAxis(2)->SetRange(0, -1);
tracksMixed = trackMixedAllStep6->Projection(1, 0, "E");
}
else
{
trackMixedAll->GetAxis(2)->SetRange(0, -1);
tracksMixed = trackMixedAll->Projection(1, 0, "E");
}
Float_t binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1);
if (stepForMixed == -1 && step == kCFStepBiasStudy && !trackMixedAllStep6)
{
Float_t phiExclude = 0.41;
mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetXaxis()->FindBin(-phiExclude)-1, tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
Double_t mixedNormError2 = 0;
Double_t mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(phiExclude)+1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2);
if (mixedNormError == 0 || mixedNormError2 == 0)
{
Printf("ERROR: Skipping multiplicity %d because mixed event is empty %f %f %f %f", multBin, mixedNorm, mixedNormError, mixedNorm2, mixedNormError2);
continue;
}
Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(-phiExclude) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
mixedNorm /= nBinsMixedNorm;
mixedNormError /= nBinsMixedNorm;
Int_t nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetXaxis()->FindBin(phiExclude) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
mixedNorm2 /= nBinsMixedNorm2;
mixedNormError2 /= nBinsMixedNorm2;
mixedNorm = mixedNorm / mixedNormError / mixedNormError + mixedNorm2 / mixedNormError2 / mixedNormError2;
mixedNormError = TMath::Sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2));
mixedNorm *= mixedNormError * mixedNormError;
}
else
{
mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
mixedNorm /= nBinsMixedNorm;
mixedNormError /= nBinsMixedNorm;
if (mixedNormError == 0)
{
Printf("ERROR: Skipping multiplicity %d because mixed event is empty %f %f", multBin, mixedNorm, mixedNormError);
continue;
}
}
if (fTrackEtaCut > 0)
{
Double_t finiteBinCorrection = -1.0 / (2*fTrackEtaCut) * binWidthEta / 2 + 1;
Printf("Finite bin correction: %f", finiteBinCorrection);
mixedNorm /= finiteBinCorrection;
mixedNormError /= finiteBinCorrection;
}
else
{
Printf("ERROR: fTrackEtaCut not set. Finite bin correction cannot be applied. Continuing anyway...");
}
Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin);
if (triggers <= 0)
{
Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
continue;
}
mixedNorm /= triggers;
mixedNormError /= triggers;
delete tracksMixed;
}
else
Printf("WARNING: Skipping mixed-event scaling! fSkipScaleMixedEvent IS set!");
if (mixedNorm <= 0)
{
Printf("ERROR: Skipping multiplicity %d because mixed event is empty at (0,0)", multBin);
continue;
}
TAxis* vertexAxis = trackSameAll->GetAxis(2);
Int_t vertexBinBegin = 1;
Int_t vertexBinEnd = vertexAxis->GetNbins();
if (fZVtxMax > fZVtxMin)
{
vertexBinBegin = vertexAxis->FindBin(fZVtxMin);
vertexBinEnd = vertexAxis->FindBin(fZVtxMax);
}
for (Int_t vertexBin = vertexBinBegin; vertexBin <= vertexBinEnd; vertexBin++)
{
trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
TH2* tracksSame = trackSameAll->Projection(1, 0, "E");
TH2* tracksMixed = trackMixedAll->Projection(1, 0, "E");
Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin);
if (triggers2 <= 0)
{
Printf("ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin);
}
else
{
if (!fSkipScaleMixedEvent)
tracksMixed->Scale(1.0 / triggers2 / mixedNorm);
else if (tracksMixed->Integral() > 0)
tracksMixed->Scale(1.0 / tracksMixed->Integral());
Double_t sums[] = { 0, 0, 0 };
Double_t errors[] = { 0, 0, 0 };
for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
{
sums[0] += tracksSame->GetBinContent(x, y);
errors[0] += tracksSame->GetBinError(x, y);
sums[1] += tracksMixed->GetBinContent(x, y);
errors[1] += tracksMixed->GetBinError(x, y);
}
tracksSame->Divide(tracksMixed);
for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
{
sums[2] += tracksSame->GetBinContent(x, y);
errors[2] += tracksSame->GetBinError(x, y);
}
for (Int_t x=0; x<3; x++)
if (sums[x] > 0)
errors[x] /= sums[x];
Printf("The correlation function %d %d has uncertainties %f %f %f (Ratio S/M %f)", multBin, vertexBin, errors[0], errors[1], errors[2], (errors[1] > 0) ? errors[0] / errors[1] : -1);
if (!totalTracks)
totalTracks = (TH2*) tracksSame->Clone("totalTracks");
else
totalTracks->Add(tracksSame);
totalEvents += eventSameAll->GetBinContent(vertexBin, multBin);
}
delete tracksSame;
delete tracksMixed;
nCorrelationFunctions++;
}
}
if (totalTracks) {
Double_t sums[] = { 0, 0, 0 };
Double_t errors[] = { 0, 0, 0 };
for (Int_t x=1; x<=totalTracks->GetNbinsX(); x++)
for (Int_t y=1; y<=totalTracks->GetNbinsY(); y++)
{
sums[0] += totalTracks->GetBinContent(x, y);
errors[0] += totalTracks->GetBinError(x, y);
}
if (sums[0] > 0)
errors[0] /= sums[0];
if (normalizePerTrigger)
{
Printf("Dividing %f tracks by %lld events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]);
if (totalEvents > 0)
totalTracks->Scale(1.0 / totalEvents);
}
if(trigger!=NULL)
*trigger = (Int_t)totalEvents;
Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1);
totalTracks->Scale(1.0 / normalization);
}
delete trackSameAll;
delete trackMixedAll;
delete trackMixedAllStep6;
delete eventSameAll;
delete eventMixedAll;
delete eventMixedAllStep6;
TH1::AddDirectory(oldStatus);
return totalTracks;
}
TH1* AliUEHist::GetTriggersAsFunctionOfMultiplicity(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
{
ResetBinLimits(fEventHist->GetGrid(step));
Int_t firstBin = fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(ptLeadMin);
Int_t lastBin = fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(ptLeadMax);
Printf("Using pT range %d --> %d", firstBin, lastBin);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
if (fZVtxMax > fZVtxMin)
{
fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
Printf("Restricting z-vtx: %f-%f", fZVtxMin, fZVtxMax);
}
TH1* eventHist = fEventHist->GetGrid(step)->Project(1);
ResetBinLimits(fEventHist->GetGrid(step));
return eventHist;
}
TH2* AliUEHist::GetSumOfRatios(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm, Bool_t useVertexBins)
{
TH2* totalTracks = 0;
Int_t totalEvents = 0;
Int_t vertexBin = 1;
TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
if (useVertexBins && !vertexAxis)
{
Printf("Vertex axis requested but not available");
return 0;
}
while (1)
{
if (useVertexBins)
{
SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
vertexBin++;
}
Int_t multIter = multBinBegin;
while (1)
{
Int_t multBinBeginLocal = multBinBegin;
Int_t multBinEndLocal = multBinEnd;
if (multBinEnd >= multBinBegin)
{
multBinBeginLocal = multIter;
multBinEndLocal = multIter;
multIter++;
}
Long64_t nEvents = 0;
TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
tracks->Scale(nEvents);
totalEvents += nEvents;
TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
tracks->Scale(mixedTwoD->Integral() / tracks->Integral());
tracks->Divide(mixedTwoD);
delete mixedTwoD;
if (!totalTracks)
totalTracks = tracks;
else
{
totalTracks->Add(tracks);
delete tracks;
}
if (multIter > multBinEnd)
break;
}
if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
break;
}
if (useVertexBins)
totalEvents = vertexAxis->GetNbins();
Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
if (totalEvents > 0)
totalTracks->Scale(1.0 / totalEvents);
return totalTracks;
}
TH1* AliUEHist::GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization)
{
ResetBinLimits(fTrackHist[region]->GetGrid(step));
ResetBinLimits(fEventHist->GetGrid(step));
TH1D* tracks = 0;
if (multBinEnd > multBinBegin)
Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
Printf("phi min = %.2f and max = %.2f requested; due to binning range will be from %.2f to %.2f, i.e. a %.0f%% larger window", phiMin, phiMax, phiMinReal, phiMaxReal, (phiMaxReal - phiMinReal) / (phiMax - phiMin) * 100 - 100);
Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
TH1D* events = fEventHist->ShowProjection(0, step);
for (Int_t bin=firstBin; bin<=lastBin; bin++)
{
fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
Float_t normalization = 1;
TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
if (strcmp(axis->GetTitle(), "#eta") == 0)
{
Printf("Normalizing using eta-axis range");
normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
}
else
{
Printf("Normalizating assuming accepted range of |eta| < 0.8");
normalization *= 0.8 * 2;
}
if (!skipPhiNormalization)
normalization *= phiMaxReal - phiMinReal;
tracksTmp->Scale(1.0 / normalization);
for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
{
tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
tracksTmp->SetBinError (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
}
Int_t nEvents = (Int_t) events->GetBinContent(bin);
if (nEvents > 0)
tracksTmp->Scale(1.0 / nEvents);
Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
if (!tracks)
tracks = tracksTmp;
else
{
tracks->Add(tracksTmp);
delete tracksTmp;
}
}
delete events;
ResetBinLimits(fTrackHist[region]->GetGrid(step));
ResetBinLimits(fEventHist->GetGrid(step));
return tracks;
}
void AliUEHist::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
{
if (grid != target)
target->Reset();
if (histogram != 0)
{
if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
}
if (grid->GetNdimensions() > 6)
AliFatal("Too many dimensions in THnSparse");
Int_t bins[6];
for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
{
Double_t value = grid->GetBinContent(binIdx, bins);
Double_t error = grid->GetBinError(binIdx);
if (histogram != 0)
{
if (var2 < 0)
{
value *= histogram->GetBinContent(bins[var1]);
error *= histogram->GetBinContent(bins[var1]);
}
else
{
value *= histogram->GetBinContent(bins[var1], bins[var2]);
error *= histogram->GetBinContent(bins[var1], bins[var2]);
}
}
target->SetBinContent(bins, value);
target->SetBinError(bins, error);
}
}
void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
{
for (UInt_t region=0; region<fkRegions; region++)
CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
}
void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
{
if (!fTrackHist[region])
return;
THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
Float_t entriesBefore = grid->GetEntries();
MultiplyHistograms(grid, target, trackCorrection, var1, var2);
Printf("AliUEHist::CorrectTracks: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0);
}
void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
{
AliCFGridSparse* grid = fEventHist->GetGrid(step1);
AliCFGridSparse* target = fEventHist->GetGrid(step2);
Float_t entriesBefore = grid->GetEntries();
MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);
Printf("AliUEHist::CorrectEvents: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0);
}
void AliUEHist::Correct(AliUEHist* corrections)
{
if (fHistogramType.Length() == 0)
{
Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
if (fTrackHist[kToward]->GetNVar() < 5)
{
if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#varphid#eta") == 0)
fHistogramType = "NumberDensitypT";
else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#eta") == 0)
fHistogramType = "SumpT";
}
else if (fTrackHist[kToward]->GetNVar() == 5)
{
if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
fHistogramType = "NumberDensityPhi";
else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
fHistogramType = "NumberDensityPhiCentrality";
}
if (fHistogramType.Length() == 0)
AliFatal("Cannot figure out histogram type. Try setting it manually...");
}
Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
{
Bool_t biasFromMC = kFALSE;
const char* projAxis = "z";
Int_t secondBin = -1;
if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
{
projAxis = "yz";
secondBin = 4;
}
for (UInt_t region = 0; region < fkRegions; region++)
{
if (!fTrackHist[region])
continue;
TH1* leadingBiasTracks = 0;
if (biasFromMC)
{
leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1);
Printf("WARNING: Using MC bias correction");
}
else
leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);
CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
if (region == kMin && fCombineMinMax)
{
CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
delete leadingBiasTracks;
break;
}
delete leadingBiasTracks;
}
TH1* leadingBiasEvents = 0;
if (biasFromMC)
leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2);
else
leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);
CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
delete leadingBiasEvents;
TH2D* contamination = corrections->GetTrackingContamination();
if (corrections->fContaminationEnhancement)
{
Printf("Applying contamination enhancement");
for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
}
CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
delete contamination;
TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
delete efficiencyCorrection;
efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
delete efficiencyCorrection;
CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
vertexCorrectionObs->Reset();
TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
func->SetParameters(0.1, 1, -0.7);
vertexCorrection->Fit(func, "0I", "", 0, 3);
for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
{
Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
if (xPos < 3)
vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
else
vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
}
#if 0
new TCanvas;
vertexCorrection->DrawCopy();
vertexCorrectionObs->SetLineColor(2);
vertexCorrectionObs->DrawCopy("same");
func->SetRange(0, 4);
func->DrawClone("same");
#endif
CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
delete vertexCorrectionObs;
delete vertexCorrection;
delete func;
CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
}
else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
{
if (fTrackHist[0]->GetNVar() <= 5)
{
CFStep step = kCFStepBiasStudy;
CorrectTracks(step, kCFStepTracked, 0, -1);
CorrectEvents(step, kCFStepTracked, 0, -1);
TH1* contamination = corrections->GetTrackingContamination(1);
if (0)
{
Printf("Applying contamination enhancement");
for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
{
printf("%f", contamination->GetBinContent(bin));
if (contamination->GetBinContent(bin) > 0)
contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
printf(" --> %f\n", contamination->GetBinContent(bin));
}
}
CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
delete contamination;
TH2* correlatedContamination = corrections->GetCorrelatedContamination();
if (0)
{
Printf("Applying contamination enhancement");
for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
{
printf("%f", correlatedContamination->GetBinContent(bin, bin2));
if (correlatedContamination->GetBinContent(bin, bin2) > 0)
correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
}
}
CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
delete correlatedContamination;
CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
delete efficiencyCorrection;
efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
delete efficiencyCorrection;
CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
}
else
{
Printf(">>>>>>>> Applying corrections in place to reduce memory consumption");
CFStep step = kCFStepBiasStudy;
TH1* contamination = corrections->GetTrackingContamination(1);
if (0)
{
Printf("Applying contamination enhancement");
for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
{
printf("%f", contamination->GetBinContent(bin));
if (contamination->GetBinContent(bin) > 0)
contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
printf(" --> %f\n", contamination->GetBinContent(bin));
}
}
CorrectTracks(step, step, contamination, 1);
delete contamination;
contamination = corrections->GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, -1, 2);
new TCanvas; contamination->DrawCopy();
CorrectEvents(step, step, contamination, 0);
CorrectTracks(step, step, contamination, 2);
delete contamination;
if (0)
{
TH2* correlatedContamination = corrections->GetCorrelatedContamination();
if (1)
{
Printf("Applying contamination enhancement");
for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
{
printf("%f", correlatedContamination->GetBinContent(bin, bin2));
if (correlatedContamination->GetBinContent(bin, bin2) > 0)
correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
}
}
CorrectCorrelatedContamination(step, 0, correlatedContamination);
delete correlatedContamination;
}
else
Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
CorrectTracks(step, step, efficiencyCorrection, 1, 3);
delete efficiencyCorrection;
efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
CorrectEvents(step, step, efficiencyCorrection, 0, 1);
CorrectTracks(step, step, efficiencyCorrection, 2, 3);
delete efficiencyCorrection;
}
}
else
AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
}
THnBase* AliUEHist::GetTrackEfficiencyND(CFStep step1, CFStep step2)
{
AliCFContainer* sourceContainer = fTrackHistEfficiency;
step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
ResetBinLimits(sourceContainer->GetGrid(step1));
ResetBinLimits(sourceContainer->GetGrid(step2));
if (fEtaMax > fEtaMin)
{
Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
}
Int_t dimensions[] = { 0, 1, 3, 4 };
THnBase* generated = sourceContainer->GetGrid(step1)->GetGrid()->ProjectionND(4, dimensions);
THnBase* measured = sourceContainer->GetGrid(step2)->GetGrid()->ProjectionND(4, dimensions);
ResetBinLimits(sourceContainer->GetGrid(step1));
ResetBinLimits(sourceContainer->GetGrid(step2));
measured->Divide(measured, generated, 1, 1, "B");
delete generated;
return measured;
}
TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
{
AliCFContainer* sourceContainer = 0;
if (source == 0)
{
if (!fCache)
{
fCache = (AliCFContainer*) fTrackHist[0]->Clone();
for (UInt_t i = 1; i < fkRegions; i++)
if (fTrackHist[i])
fCache->Add(fTrackHist[i]);
}
sourceContainer = fCache;
}
else if (source == 1 || source == 2)
{
sourceContainer = fTrackHistEfficiency;
step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
}
else
return 0;
ResetBinLimits(sourceContainer->GetGrid(step1));
ResetBinLimits(sourceContainer->GetGrid(step2));
if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
{
Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
}
if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
{
Printf("Restricted pt-range to %f %f", fPtMin, fPtMax);
sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
}
if (fPartSpecies != -1 && axis1 != 2 && axis2 != 2 && axis3 != 2)
{
Printf("Restricted to particle species %d", fPartSpecies);
sourceContainer->GetGrid(step1)->SetRangeUser(2, fPartSpecies, fPartSpecies);
sourceContainer->GetGrid(step2)->SetRangeUser(2, fPartSpecies, fPartSpecies);
}
if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
{
Printf("Restricted centrality range to %f %f", fCentralityMin, fCentralityMax);
sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
}
if (fZVtxMax > fZVtxMin && axis1 != 4 && axis2 != 4 && axis3 != 4)
{
Printf("Restricted z-vtx range to %f %f", fZVtxMin, fZVtxMax);
sourceContainer->GetGrid(step1)->SetRangeUser(4, fZVtxMin, fZVtxMax);
sourceContainer->GetGrid(step2)->SetRangeUser(4, fZVtxMin, fZVtxMax);
}
TH1* measured = 0;
TH1* generated = 0;
if (axis3 >= 0)
{
generated = sourceContainer->Project(step1, axis1, axis2, axis3);
measured = sourceContainer->Project(step2, axis1, axis2, axis3);
}
else if (axis2 >= 0)
{
generated = sourceContainer->Project(step1, axis1, axis2);
measured = sourceContainer->Project(step2, axis1, axis2);
}
else
{
generated = sourceContainer->Project(step1, axis1);
measured = sourceContainer->Project(step2, axis1);
}
Int_t binBegin[3];
Int_t binEnd[3];
binBegin[0] = 1;
binBegin[1] = 1;
binBegin[2] = 1;
binEnd[0] = generated->GetNbinsX();
binEnd[1] = generated->GetNbinsY();
binEnd[2] = generated->GetNbinsZ();
if (fEtaMax > fEtaMin)
{
if (axis1 == 0)
{
binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
binEnd[0] = generated->GetXaxis()->FindBin(fEtaMax);
}
if (axis2 == 0)
{
binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
binEnd[1] = generated->GetYaxis()->FindBin(fEtaMax);
}
if (axis3 == 0)
{
binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
binEnd[2] = generated->GetZaxis()->FindBin(fEtaMax);
}
}
if (fPtMax > fPtMin)
{
Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
if (axis1 == 1)
{
binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
binEnd[0] = generated->GetXaxis()->FindBin(ptMax);
}
if (axis2 == 1)
{
binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
binEnd[1] = generated->GetYaxis()->FindBin(ptMax);
}
if (axis3 == 1)
{
binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
binEnd[2] = generated->GetZaxis()->FindBin(ptMax);
}
}
Int_t total = 0;
Int_t count = 0;
Int_t vars[3];
for (Int_t i=0; i<3; i++)
vars[i] = binBegin[i];
const Int_t limit = 50;
while (1)
{
if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
{
Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
count++;
}
else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
{
Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
generated->GetBinContent(vars[0], vars[1]));
count++;
}
else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
{
Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
generated->GetBinContent(vars[0], vars[1], vars[2]));
count++;
}
vars[2]++;
if (vars[2] == binEnd[2]+1)
{
vars[2] = binBegin[2];
vars[1]++;
}
if (vars[1] == binEnd[1]+1)
{
vars[1] = binBegin[1];
vars[0]++;
}
if (vars[0] == binEnd[0]+1)
break;
total++;
}
Printf("Correction has %d empty bins (out of %d bins)", count, total);
if (source == 2)
{
TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
if (axis->GetNbins() < measured->GetNbinsX())
{
if (axis2 != -1)
{
TH1* tmp = measured;
measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
{
((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
measured->SetBinError(x, y, 0);
}
delete tmp;
tmp = generated;
generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
{
((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
generated->SetBinError(x, y, 0);
}
delete tmp;
}
else
{
TH1* tmp = measured;
measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
delete tmp;
tmp = generated;
generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
delete tmp;
}
}
else if (axis->GetNbins() > measured->GetNbinsX())
{
if (axis2 != -1)
AliFatal("Rebinning only works for 1d at present");
if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
Double_t newBins[] = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
TH1* tmp = measured;
measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
delete tmp;
tmp = measured;
measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
{
measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
}
delete tmp;
tmp = generated;
generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
delete tmp;
tmp = generated;
generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
{
generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
}
delete tmp;
}
}
measured->Divide(measured, generated, 1, 1, "B");
delete generated;
ResetBinLimits(sourceContainer->GetGrid(step1));
ResetBinLimits(sourceContainer->GetGrid(step2));
return measured;
}
TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
{
if (ptLeadMax > ptLeadMin)
{
fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
}
TH1* measured = 0;
TH1* generated = 0;
if (axis2 >= 0)
{
generated = fEventHist->Project(step1, axis1, axis2);
measured = fEventHist->Project(step2, axis1, axis2);
}
else
{
generated = fEventHist->Project(step1, axis1);
measured = fEventHist->Project(step2, axis1);
}
measured->Divide(measured, generated, 1, 1, "B");
delete generated;
if (ptLeadMax > ptLeadMin)
{
fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
}
return measured;
}
void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
{
if (hist1->GetNbinsZ() != hist2->GetNbinsX())
AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
{
for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
{
for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
{
if (hist2->GetBinContent(z) > 0)
{
hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
}
else
{
hist1->SetBinContent(x, y, z, 0);
hist1->SetBinError(x, y, z, 0);
}
}
}
}
}
TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
{
AliCFContainer* tmp = 0;
if (region == -1)
{
tmp = (AliCFContainer*) fTrackHist[0]->Clone();
for (UInt_t i = 1; i < fkRegions; i++)
if (fTrackHist[i])
tmp->Add(fTrackHist[i]);
}
else if (region == kMin && fCombineMinMax)
{
tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
tmp->Add(fTrackHist[kMax]);
}
else
tmp = fTrackHist[region];
ResetBinLimits(tmp->GetGrid(step1));
ResetBinLimits(fEventHist->GetGrid(step1));
SetBinLimits(tmp->GetGrid(step1));
ResetBinLimits(tmp->GetGrid(step2));
ResetBinLimits(fEventHist->GetGrid(step2));
SetBinLimits(tmp->GetGrid(step2));
TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
if (weighting == 0)
WeightHistogram(hist1, events1);
TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
if (weighting == 0)
WeightHistogram(hist2, events2);
TH1* generated = hist1;
TH1* measured = hist2;
if (weighting == 0 || weighting == 1)
{
if (axis)
{
if (leadPtMax > leadPtMin)
{
hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
}
if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
{
hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
}
generated = hist1->Project3D(axis);
measured = hist2->Project3D(axis);
delete hist1;
delete hist2;
}
delete events1;
delete events2;
}
else if (weighting == 2)
{
delete hist1;
delete hist2;
generated = events1;
measured = events2;
}
measured->Divide(generated);
delete generated;
ResetBinLimits(tmp->GetGrid(step1));
ResetBinLimits(tmp->GetGrid(step2));
if ((region == -1) || (region == kMin && fCombineMinMax))
delete tmp;
return measured;
}
void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
{
if (!fTrackHist[region])
return;
THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
Int_t var1 = 1;
Int_t var2 = 2;
if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
{
Int_t bins[6];
Double_t value = grid->GetBinContent(binIdx, bins);
Double_t error = grid->GetBinError(binIdx);
if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
continue;
if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
continue;
value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
grid->SetBinContent(bins, value);
grid->SetBinError(bins, error);
}
Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
}
TH2* AliUEHist::GetCorrelatedContamination()
{
Int_t step1 = kCFStepTrackedOnlyPrim;
Int_t step2 = kCFStepTracked;
fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01);
fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01);
TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01);
fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01);
TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
tracksStep1->Divide(tracksStep2);
TH1* triggersStep1 = fEventHist->Project(step1, 0);
TH1* triggersStep2 = fEventHist->Project(step2, 0);
TH1* singleParticle = GetTrackingContamination(1);
for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
if (singleParticle->GetBinContent(x) > 0 && triggersStep1->GetBinContent(y) > 0)
tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
else
tracksStep1->SetBinContent(x, y, 0);
delete singleParticle;
delete tracksStep2;
delete triggersStep1;
delete triggersStep2;
return tracksStep1;
}
TH2D* AliUEHist::GetTrackingEfficiency()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
}
TH2D* AliUEHist::GetFakeRate()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), 0, 1));
}
TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
}
TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
{
return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
}
TH1D* AliUEHist::GetFakeRate(Int_t axis)
{
return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), axis));
}
TH2D* AliUEHist::GetTrackingCorrection()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
}
TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
{
return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
}
TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
}
TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
}
TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
{
return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
}
TH2D* AliUEHist::GetTrackingContamination()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
}
TH2D* AliUEHist::GetTrackingContaminationCentrality()
{
return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
}
TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
{
return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
}
const char* AliUEHist::GetRegionTitle(Region region)
{
switch (region)
{
case kToward:
return "Towards";
case kAway:
return "Away";
case kMin:
return (fCombineMinMax) ? "Transverse" : "Min";
case kMax:
return "Max";
}
return 0;
}
const char* AliUEHist::GetStepTitle(CFStep step)
{
switch (step)
{
case kCFStepAll:
return "All events";
case kCFStepTriggered:
return "Triggered";
case kCFStepVertex:
return "Primary Vertex";
case kCFStepAnaTopology:
return "Required analysis topology";
case kCFStepTrackedOnlyPrim:
return "Tracked (matched MC, only primaries)";
case kCFStepTracked:
return "Tracked (matched MC, all)";
case kCFStepReconstructed:
return "Reconstructed";
case kCFStepRealLeading:
return "Correct leading particle identified";
case kCFStepBiasStudy:
return "Bias study applying tracking efficiency";
case kCFStepBiasStudy2:
return "Bias study applying tracking efficiency in two steps";
case kCFStepCorrected:
return "Corrected for efficiency on-the-fly";
}
return 0;
}
void AliUEHist::CopyReconstructedData(AliUEHist* from)
{
for (Int_t region=0; region<4; region++)
{
if (!fTrackHist[region])
continue;
fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
}
fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy, from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
}
void AliUEHist::DeepCopy(AliUEHist* from)
{
for (Int_t region=0; region<4; region++)
{
if (!fTrackHist[region] || !from->fTrackHist[region])
continue;
for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
{
Printf("Copying region %d step %d", region, step);
THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
THnSparse* source = from->fTrackHist[region]->GetGrid(step)->GetGrid();
target->Reset();
target->RebinnedAdd(source);
}
}
for (Int_t step=0; step<fEventHist->GetNStep(); step++)
{
Printf("Ev: Copying step %d", step);
THnSparse* target = fEventHist->GetGrid(step)->GetGrid();
THnSparse* source = from->fEventHist->GetGrid(step)->GetGrid();
target->Reset();
target->RebinnedAdd(source);
}
for (Int_t step=0; step<TMath::Min(fTrackHistEfficiency->GetNStep(), from->fTrackHistEfficiency->GetNStep()); step++)
{
if (!from->fTrackHistEfficiency->GetGrid(step))
continue;
Printf("Eff: Copying step %d", step);
THnSparse* target = fTrackHistEfficiency->GetGrid(step)->GetGrid();
THnSparse* source = from->fTrackHistEfficiency->GetGrid(step)->GetGrid();
target->Reset();
target->RebinnedAdd(source);
}
}
void AliUEHist::SymmetrizepTBins()
{
for (Int_t region=0; region<4; region++)
{
if (!fTrackHist[region])
continue;
for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
{
Printf("Copying region %d step %d", region, step);
THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
if (target->GetEntries() == 0)
continue;
THnSparse* source = (THnSparse*) target->Clone();
Int_t zVtxBins = 1;
if (target->GetNdimensions() > 5)
zVtxBins = target->GetAxis(5)->GetNbins();
for (Int_t i3 = 1; i3 <= target->GetAxis(3)->GetNbins(); i3++)
for (Int_t i5 = 1; i5 <= zVtxBins; i5++)
{
for (Int_t i1 = 1; i1 <= target->GetAxis(1)->GetNbins(); i1++)
for (Int_t i2 = 1; i2 <= target->GetAxis(2)->GetNbins(); i2++)
{
Int_t binA = target->GetAxis(1)->FindBin(target->GetAxis(2)->GetBinCenter(i2));
Int_t binT = target->GetAxis(2)->FindBin(target->GetAxis(1)->GetBinCenter(i1));
Printf("(%d %d) Copying from %d %d to %d %d", i3, i5, binA, binT, i1, i2);
for (Int_t i0 = 1; i0 <= target->GetAxis(0)->GetNbins(); i0++)
for (Int_t i4 = 1; i4 <= target->GetAxis(4)->GetNbins(); i4++)
{
Int_t binEta = target->GetAxis(0)->FindBin(-target->GetAxis(0)->GetBinCenter(i0));
Double_t phi = -target->GetAxis(4)->GetBinCenter(i4);
if (phi < -TMath::Pi()/2)
phi += TMath::TwoPi();
Int_t binPhi = target->GetAxis(4)->FindBin(phi);
Int_t binSource[] = { binEta, binA, binT, i3, binPhi, i5 };
Int_t binTarget[] = { i0, i1, i2, i3, i4, i5 };
Double_t value = source->GetBinContent(binSource);
Double_t error = source->GetBinError(binSource);
if (error == 0)
continue;
Double_t value2 = target->GetBinContent(binTarget);
Double_t error2 = target->GetBinError(binTarget);
Double_t sum = value;
Double_t err = error;
if (error2 > 0)
{
sum = value + value2;
err = TMath::Sqrt(error * error + error2 * error2);
}
target->SetBinContent(binTarget, sum);
target->SetBinError(binTarget, err);
}
}
}
delete source;
}
}
}
void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
{
Float_t fitRangeBegin = 5.01;
Float_t fitRangeEnd = 14.99;
Float_t extendRangeBegin = 10.01;
if (fTrackHistEfficiency->GetNVar() == 3)
{
TH1* obj = GetTrackingEfficiency(1);
if (verbose)
{
new TCanvas;
obj->Draw();
}
obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
obj = GetTrackingContamination(1);
if (verbose)
{
new TCanvas;
obj->Draw();
}
obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
Printf("AliUEHist::ExtendTrackingEfficiency: Fitted efficiency between %f and %f and got %f tracking efficiency and %f tracking contamination correction. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, trackingEff, trackingCont, extendRangeBegin, fEtaMin, fEtaMax);
for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++)
{
Int_t bins[3];
bins[0] = x;
bins[1] = y;
bins[2] = z;
fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
}
}
else if (fTrackHistEfficiency->GetNVar() == 4)
{
Float_t* trackingEff = 0;
Float_t* trackingCont = 0;
Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
Int_t nCentralityBins = 5;
Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
for (Int_t caseNo = 0; caseNo < 2; caseNo++)
{
Float_t* target = 0;
Int_t centralityBinsLocal = nCentralityBins;
if (caseNo == 0)
{
trackingEff = new Float_t[centralityBinsLocal];
target = trackingEff;
}
else
{
centralityBinsLocal = 1;
trackingCont = new Float_t[centralityBinsLocal];
target = trackingCont;
}
for (Int_t i=0; i<centralityBinsLocal; i++)
{
if (centralityBinsLocal == 1)
SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
else
SetCentralityRange(centralityBins[i] + 0.1, centralityBins[i+1] - 0.1);
TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
if (verbose)
{
new TCanvas;
proj->DrawCopy();
}
if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
target[i] = proj->GetFunction("pol0")->GetParameter(0);
else
target[i] = 0;
Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
}
}
for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++)
{
for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++)
{
Int_t bins[4];
bins[0] = x;
bins[1] = y;
bins[2] = z;
bins[3] = z2;
Int_t z2Bin = 0;
while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
z2Bin++;
fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
if (trackingCont[0] > 0)
fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
else
fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
}
}
delete[] trackingEff;
delete[] trackingCont;
}
SetCentralityRange(0, 0);
}
void AliUEHist::AdditionalDPhiCorrection(Int_t step)
{
Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
{
Int_t bins[6];
Double_t value = grid->GetBinContent(binIdx, bins);
Double_t error = grid->GetBinError(binIdx);
Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
if (TMath::Abs(binCenter) < 0.2)
{
value *= 0.985;
error *= 0.985;
}
else if (TMath::Abs(binCenter) < 0.3)
{
value *= 0.9925;
error *= 0.9925;
}
grid->SetBinContent(bins, value);
grid->SetBinError(bins, error);
}
}
void AliUEHist::Scale(Double_t factor)
{
for (Int_t i=0; i<4; i++)
if (fTrackHist[i])
fTrackHist[i]->Scale(factor);
fEventHist->Scale(factor);
fTrackHistEfficiency->Scale(factor);
}
void AliUEHist::Reset()
{
for (Int_t i=0; i<4; i++)
if (fTrackHist[i])
for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
for (Int_t step=0; step<fEventHist->GetNStep(); step++)
fEventHist->GetGrid(step)->GetGrid()->Reset();
for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
}
THnBase* AliUEHist::ChangeToThn(THnBase* sparse)
{
Int_t nBins[10];
for (Int_t i=0; i<sparse->GetNdimensions(); i++)
nBins[i] = sparse->GetAxis(i)->GetNbins();
THn* tmpTHn = new THnF(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse->GetNdimensions(), nBins, 0, 0);
for (Int_t i=0; i<sparse->GetNdimensions(); i++)
{
tmpTHn->SetBinEdges(i, sparse->GetAxis(i)->GetXbins()->GetArray());
tmpTHn->GetAxis(i)->SetTitle(sparse->GetAxis(i)->GetTitle());
}
tmpTHn->RebinnedAdd(sparse);
return tmpTHn;
}
void AliUEHist::CondenseBin(THnSparse* grid, THnSparse* target, Int_t axis, Float_t targetValue, Float_t from, Float_t to)
{
if (grid->GetNdimensions() > 6)
AliFatal("Too many dimensions in THnSparse");
Int_t targetBin = grid->GetAxis(axis)->FindBin(targetValue);
AliInfo(Form("Target bin on axis %d with value %f is %d", axis, targetValue, targetBin));
Int_t fromBin = 1;
Int_t toBin = grid->GetAxis(axis)->GetNbins();
if (to > from)
{
fromBin = grid->GetAxis(axis)->FindBin(from);
toBin = grid->GetAxis(axis)->FindBin(to);
AliInfo(Form("Only condensing between bin %d and %d", fromBin, toBin));
}
Int_t bins[6];
for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
{
Double_t value = grid->GetBinContent(binIdx, bins);
Double_t error = grid->GetBinError(binIdx);
if (bins[axis] >= fromBin && bins[axis] <= toBin)
bins[axis] = targetBin;
value += target->GetBinContent(bins);
error = TMath::Sqrt(error * error + target->GetBinError(bins) * target->GetBinError(bins));
target->SetBinContent(bins, value);
target->SetBinError(bins, error);
}
}
void AliUEHist::CondenseBin(CFStep step, Int_t trackAxis, Int_t eventAxis, Float_t targetValue, Float_t from, Float_t to, CFStep tmpStep)
{
fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
for (UInt_t i=0; i<fkRegions; i++)
if (fTrackHist[i])
fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
CorrectTracks(step, tmpStep, 0, -1);
CorrectEvents(step, tmpStep, 0, -1);
fEventHist->GetGrid(step)->GetGrid()->Reset();
for (UInt_t i=0; i<fkRegions; i++)
if (fTrackHist[i])
fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
for (UInt_t i=0; i<fkRegions; i++)
{
if (!fTrackHist[i])
continue;
THnSparse* grid = fTrackHist[i]->GetGrid(tmpStep)->GetGrid();
THnSparse* target = fTrackHist[i]->GetGrid(step)->GetGrid();
CondenseBin(grid, target, trackAxis, targetValue, from, to);
}
CondenseBin(fEventHist->GetGrid(tmpStep)->GetGrid(), fEventHist->GetGrid(step)->GetGrid(), eventAxis, targetValue, from, to);
fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
for (UInt_t i=0; i<fkRegions; i++)
if (fTrackHist[i])
fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
}