#include <TCollection.h>
#include <TIterator.h>
#include <TList.h>
+#include <TLegend.h>
+
+#include "dNdEtaCorrection.h"
//____________________________________________________________________
ClassImp(dNdEtaAnalysis)
hEtaVsVtx->SetYTitle("#eta");
hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
+ hEtaVsVtxCheck = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_check", name)));
hVtx = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
for (Int_t i=0; i<kVertexBinning; ++i)
{
//____________________________________________________________________
void
-dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) {
- hEtaVsVtx->Fill(vtx, eta, weight);
+dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t c) {
hEtaVsVtxUncorrected->Fill(vtx,eta);
+ hEtaVsVtxCheck->Fill(vtx, eta, c);
}
//____________________________________________________________________
}
//____________________________________________________________________
-void
-dNdEtaAnalysis::Finish() {
+void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction)
+{
+ // correct with correction values if available
+ // TODO what do we do with the error?
+ if (!correction)
+ printf("INFO: No correction applied\n");
+
+ // this can be replaced by TH2F::Divide if we agree that the binning will be always the same
+ for (Int_t iVtx=0; iVtx<=hEtaVsVtxUncorrected->GetNbinsX(); iVtx++)
+ {
+ for (Int_t iEta=0; iEta<=hEtaVsVtxUncorrected->GetNbinsY(); iEta++)
+ {
+ Float_t correctionValue = 1;
+ if (correction)
+ correctionValue = correction->GetCorrection(hEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), hEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta));
+
+ Float_t value = hEtaVsVtxUncorrected->GetBinContent(iVtx, iEta);
+ Float_t error = hEtaVsVtxUncorrected->GetBinError(iVtx, iEta);
+
+ Float_t correctedValue = value * correctionValue;
+ Float_t correctedError = error * correctionValue;
+
+ if (correctedValue != 0)
+ {
+ hEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue);
+ hEtaVsVtx->SetBinError(iVtx, iEta, correctedError);
+ }
+ }
+ }
+
+ // normalize with n events (per vtx)
+ for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) {
+ Float_t nEvents = hVtx->GetBinContent(iVtx);
+ Float_t nEventsError = hVtx->GetBinError(iVtx);
- // first normalize with n events (per vtx)
- for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) {
- Float_t nEvents = hVtx->GetBinContent(i_vtx);
- Float_t nEventsError = hVtx->GetBinError(i_vtx);
-
if (nEvents==0) continue;
-
- for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
- Float_t value = hEtaVsVtx->GetBinContent(i_vtx, i_eta)/nEvents;
+
+ for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
+ Float_t value = hEtaVsVtx->GetBinContent(iVtx, iEta) / nEvents;
+ if (value==0) continue;
+ Float_t error = hEtaVsVtx->GetBinError(iVtx, iEta)/nEvents;
+ error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta)/
+ hEtaVsVtx->GetBinContent(iVtx, iEta),2) +
+ TMath::Power(nEventsError/nEvents,2));
+ hEtaVsVtx->SetBinContent(iVtx, iEta, value);
+ hEtaVsVtx->SetBinError(iVtx, iEta, error);
+ }
+
+ //debug
+ for (Int_t iEta=0; iEta<=hEtaVsVtxCheck->GetNbinsY(); iEta++) {
+ Float_t value = hEtaVsVtxCheck->GetBinContent(iVtx, iEta) / nEvents;
if (value==0) continue;
- Float_t error = hEtaVsVtx->GetBinError(i_vtx, i_eta)/nEvents;
- error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta)/
- hEtaVsVtx->GetBinContent(i_vtx, i_eta),2) +
+ Float_t error = hEtaVsVtxCheck->GetBinError(iVtx, iEta)/nEvents;
+ error = TMath::Sqrt(TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta)/
+ hEtaVsVtxCheck->GetBinContent(iVtx, iEta),2) +
TMath::Power(nEventsError/nEvents,2));
- hEtaVsVtx->SetBinContent(i_vtx, i_eta, value);
- hEtaVsVtx->SetBinError(i_vtx, i_eta, error);
+ hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value);
+ hEtaVsVtxCheck->SetBinError(iVtx, iEta, error);
}
}
// then take the wieghted average for each eta
// is this the right way to do it???
- for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
+ for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
Float_t sum = 0;
Float_t sumError2 = 0;
Int_t nMeasurements = 0;
- Float_t sum_xw = 0;
- Float_t sum_w = 0;
+ Float_t sumXw = 0;
+ Float_t sumW = 0;
// do we have several histograms for different vertex positions?
- Int_t vertexBinWidth = hVtx->GetNbinsX() / kVertexBinning;
+ Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
{
- Int_t vertexBinBegin = vertexBinWidth * vertexPos;
- Int_t vertexBinEnd = vertexBinBegin + vertexBinWidth;
+ Int_t vertexBinBegin = 0;
+ Int_t vertexBinEnd = hVtx->GetNbinsX();
+
+ // the first histogram is always for the whole vertex range
+ if (vertexPos > 0)
+ {
+ vertexBinBegin = vertexBinWidth * (vertexPos-1);
+ vertexBinEnd = vertexBinBegin + vertexBinWidth;
+ }
- for (Int_t i_vtx=vertexBinBegin; i_vtx<=vertexBinEnd; i_vtx++) {
- if (hVtx->GetBinContent(i_vtx)==0) continue;
- if (hEtaVsVtx->GetBinContent(i_vtx, i_eta)==0) continue;
+ for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) {
+ if (hVtx->GetBinContent(iVtx)==0) continue;
+ if (hEtaVsVtx->GetBinContent(iVtx, iEta)==0) continue;
- Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
- sum_xw = sum_xw + hEtaVsVtx->GetBinContent(i_vtx, i_eta)*w;
- sum_w = sum_w + w;
+ Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
+ sumXw = sumXw + hEtaVsVtx->GetBinContent(iVtx, iEta)*w;
+ sumW = sumW + w;
- sum = sum + hEtaVsVtx->GetBinContent(i_vtx, i_eta);
- sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
+ sum = sum + hEtaVsVtx->GetBinContent(iVtx, iEta);
+ sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
nMeasurements++;
}
Float_t dndeta = 0;
dndeta = sum/Float_t(nMeasurements);
error = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
- dndeta = sum_xw/sum_w;
- error = 1/TMath::Sqrt(sum_w);
+ dndeta = sumXw/sumW;
+ error = 1/TMath::Sqrt(sumW);
- dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(i_eta);
- error = error/hdNdEta[vertexPos]->GetBinWidth(i_eta);
+ dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
+ error = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
- hdNdEta[vertexPos]->SetBinContent(i_eta, dndeta);
- hdNdEta[vertexPos]->SetBinError(i_eta, error);
+ hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
+ hdNdEta[vertexPos]->SetBinError(iEta, error);
}
}
-
}
}
gDirectory->cd(GetName());
hEtaVsVtx ->Write();
+ hEtaVsVtxUncorrected->Write();
hVtx ->Write();
for (Int_t i=0; i<kVertexBinning; ++i)
hdNdEta[i] ->Write();
hEtaVsVtx->Draw("COLZ");
canvas->cd(2);
+ if (hEtaVsVtxCheck)
+ hEtaVsVtxCheck->Draw("COLZ");
+
+ canvas->cd(3);
if (hEtaVsVtxUncorrected)
hEtaVsVtxUncorrected->Draw("COLZ");
- canvas->cd(3);
+ /*canvas->cd(3);
if (hVtx)
- hVtx->Draw();
+ hVtx->Draw();*/
canvas->cd(4);
if (hdNdEta[0])
hdNdEta[0]->Draw();
+
+ // histograms for different vertices?
+ if (kVertexBinning > 0)
+ {
+ // doesnt work, but i dont get it, giving up...
+ /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
+ //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
+ //printf("%d\n", yPads);
+ //canvas2->Divide(2, yPads);
+
+ TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
+
+ for (Int_t i=1; i<kVertexBinning; ++i)
+ {
+ //canvas2->cd(i-1);
+ //printf("%d\n", i);
+ if (hdNdEta[i])
+ {
+ hdNdEta[i]->SetLineColor(i);
+ hdNdEta[i]->Draw((i == 1) ? "" : "SAME");
+ legend->AddEntry(hdNdEta[i], Form("Vtx Bin %d", i-1));
+ }
+ }
+
+ legend->Draw();
+ }
}
Long64_t dNdEtaAnalysis::Merge(TCollection* list)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff