+void GetAverageCorrectionFactor(Float_t etaRange = 1.5, Float_t vertexRange = 9.9, const char* rawFile = "analysis_esd_raw.root", const char* mcFile = "analysis_mc.root")
+{
+ loadlibs();
+
+ TFile::Open(rawFile);
+ dNdEtaAnalysis* raw = new dNdEtaAnalysis("dndeta", "dndeta");
+ raw->LoadHistograms("fdNdEtaAnalysisESD");
+ raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->GetXaxis()->SetRangeUser(-vertexRange, vertexRange);
+ tracks = raw->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Project3D("y");
+ events = raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("events", 0, raw->GetData()->GetEventCorrection()->GetMeasuredHistogram()->GetNbinsY() + 1);
+ Float_t nEvents = events->Integral(events->FindBin(-vertexRange), events->FindBin(vertexRange));
+ tracks->Scale(1.0 / nEvents / tracks->GetBinWidth(1));
+
+ TFile::Open(mcFile);
+ dNdEtaAnalysis* mc = new dNdEtaAnalysis("dndeta", "dndeta");
+ mc->LoadHistograms("dndetaTrVtx");
+ mcH = mc->GetdNdEtaPtCutOffCorrectedHistogram(0);
+
+ new TCanvas;
+ mcH->SetLineColor(2);
+ mcH->DrawCopy();
+ tracks->DrawCopy("SAME");
+
+ new TCanvas;
+ mcH->GetYaxis()->SetRangeUser(0, 5);
+ mcH->Divide(tracks);
+ mcH->DrawCopy();
+ mcH->Fit("pol0", "", "", -etaRange, etaRange);
+}
+
+void TrackCuts_Comparison_MC(char* histName, Int_t plotWhich = 0, const char* fileName = "correction_map.root", Bool_t mirror = kFALSE)
+{
+ // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
+ // --> manually disable it in the run.C
+ //
+ // plotWhich: 0 = only before
+ // 1 = both
+ // 2 = only after
+ //
+ // mirror: kTRUE --> project negative values on the positive side
+
+
+ file = TFile::Open(fileName);
+
+ Int_t count = 0;
+ Int_t colors[] = { 1, 2, 3, 4, 5, 6 };
+
+ TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
+ TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
+ TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
+
+ TCanvas* c1 = new TCanvas(histName, histName, 800, 1200);
+ c1->Divide(1, 2);
+ //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
+ //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
+
+ const char* folders2[] = { "before_cuts", "after_cuts" };
+ Bool_t first = kTRUE;
+ for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
+ {
+ const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
+ const char* names[] = { "all", "primaries", "secondaries" };
+ TH1* base = 0;
+ TH1* prim = 0;
+ TH1* sec = 0;
+ for (Int_t i = 0; i < 3; i++)
+ {
+ TString folder;
+ folder.Form("%s/%s/%s", folders1[i], folders2[j], histName);
+ TH1* hist = (TH1*) file->Get(folder);
+
+ if (mirror)
+ {
+ for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
+ {
+ Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
+ if (bin != newBin)
+ {
+ hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
+ hist->SetBinContent(bin, 0);
+ }
+ }
+ }
+
+ legend->AddEntry(hist, Form("%s %s", names[i], folders2[j]));
+
+ c1->cd(1);
+ hist->SetLineColor(colors[count]);
+ hist->DrawCopy((count == 0) ? "" : "SAME");
+
+ switch (i)
+ {
+ case 0: base = hist; break;
+ case 1: prim = hist; break;
+ case 2: sec = hist; break;
+ }
+
+ count++;
+ }
+
+ TH1* eff = (TH1*) prim->Clone("eff"); eff->Reset();
+ TH1* purity = (TH1*) prim->Clone("purity"); purity->Reset();
+
+ for (Int_t bin = 1; bin <= prim->GetNbinsX(); bin++)
+ {
+ eff->SetBinContent(bin, prim->Integral(1, bin) / prim->Integral(1, prim->GetNbinsX() + 1));
+ if (prim->Integral(1, bin) + sec->Integral(1, bin) > 0)
+ purity->SetBinContent(bin, sec->Integral(1, bin) / (prim->Integral(1, bin) + sec->Integral(1, bin)));
+ }
+
+ eff->GetYaxis()->SetRangeUser(0, 1);
+ eff->SetLineColor(colors[0+j*2]);
+ eff->SetStats(kFALSE);
+ purity->SetLineColor(colors[1+j*2]);
+
+ legend3->AddEntry(eff, Form("%s: efficiency", folders2[j]));
+ legend3->AddEntry(purity, Form("%s: contamination", folders2[j]));
+
+ c1->cd(2);
+ eff->DrawCopy((first) ? "" : "SAME");
+ first = kFALSE;
+ purity->DrawCopy("SAME");
+ }
+
+ c1->cd(1)->SetLogy();
+ c1->cd(1)->SetGridx();
+ c1->cd(1)->SetGridy();
+ legend->Draw();
+
+ //c2->cd();
+ // c2->SetGridx();
+ // c2->SetGridy();
+ //legend2->Draw();
+
+ c1->cd(2)->SetGridx();
+ c1->cd(2)->SetGridy();
+ legend3->Draw();
+
+ //c1->SaveAs(Form("%s.png", histName));
+}
+
+void TrackCuts_Comparison_Data(char* histName, Int_t plotWhich, const char* fileName1, const char* fileName2, Bool_t mirror = kFALSE, const char* label1 = "file1", const char* label2 = "file2")
+{
+ // for the nsigmaplot it is needed to run with all cuts except the nsigmatovertex
+ // --> manually disable it in the run.C
+ //
+ // plotWhich: 0 = only before
+ // 1 = both
+ // 2 = only after
+ //
+ // mirror: kTRUE --> project negative values on the positive side
+
+
+ Int_t count = 0;
+ Int_t colors[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
+
+ TLegend* legend = new TLegend(0.5, 0.7, 1, 1);
+ legend->SetFillColor(0);
+ TLegend* legend2 = new TLegend(0.4, 0.6, 1, 1);
+ TLegend* legend3 = new TLegend(0.6, 0.5, 1, 0.7);
+
+ TCanvas* c1 = new TCanvas(histName, histName, 600, 600);
+ //TCanvas* c2 = new TCanvas("c2", "c2", 800, 600);
+ //TCanvas* c3 = new TCanvas("c3", "c3", 800, 600);
+
+ const char* folders2[] = { "before_cuts", "after_cuts" };
+ Bool_t first = kTRUE;
+ for (Int_t j = ((plotWhich < 2) ? 0 : 1); j < ((plotWhich > 0) ? 2 : 1); j++)
+ {
+ const char* folders1[] = { "esd_track_cuts", "esd_track_cuts_primaries", "esd_track_cuts_secondaries" };
+ const char* names[] = { "all", "primaries", "secondaries" };
+
+ Float_t normalize[3];
+
+ for (Int_t i = 0; i < 2; i++)
+ {
+ file = TFile::Open((i == 0) ? fileName1 : fileName2);
+
+ for (Int_t k = 1; k < 3; k++)
+ {
+ TString folder;
+ folder.Form("%s/%s/%s", folders1[k], folders2[j], histName);
+ Printf("%s", folder.Data());
+ TH1* hist = (TH1*) file->Get(folder);
+
+ if (mirror)
+ {
+ for (Int_t bin=1; bin<=hist->GetXaxis()->FindBin(-0.0001); bin++)
+ {
+ Int_t newBin = hist->GetXaxis()->FindBin(-hist->GetXaxis()->GetBinCenter(bin));
+ if (bin != newBin)
+ {
+ hist->Fill(-hist->GetXaxis()->GetBinCenter(bin), hist->GetBinContent(bin));
+ hist->SetBinContent(bin, 0);
+ }
+ }
+ }
+
+ if (i == 0)
+ {
+ normalize[k] = hist->Integral();
+ }
+ else
+ hist->Scale(normalize[k] / hist->Integral());
+
+ legend->AddEntry(hist, Form("%s %s %s", (i == 0) ? label1 : label2, (k == 1) ? "primaries" : "secondaries", folders2[j]));
+
+ c1->cd();
+ hist->SetStats(0);
+ hist->SetLineColor(colors[count]);
+ hist->DrawCopy((count == 0) ? "" : "SAME");
+
+ count++;
+ }
+ }
+
+ }
+
+ //c1->SetLogy();
+ c1->SetGridx();
+ c1->SetGridy();
+ legend->Draw();
+}
+
+void TrackCuts_DCA()
+{
+ file = TFile::Open("correction_map.root");
+ hist = (TH2*) file->Get("esd_track_cuts/before_cuts/dXYvsDZ");
+
+ TCanvas* c1 = new TCanvas("c1", "c1", 600, 600);
+ c1->SetLogz();
+ c1->SetRightMargin(0.12);
+ c1->SetBottomMargin(0.12);
+
+ hist->SetStats(kFALSE);
+ hist->Draw("COLZ");
+
+ ellipse = new TEllipse(0, 0, 4);
+ ellipse->SetLineWidth(2);
+ ellipse->SetLineStyle(2);
+ ellipse->SetFillStyle(0);
+ ellipse->Draw();
+
+ c1->SaveAs("trackcuts_dca_2d.eps");
+}
+
+void FindNSigma(TH2* hist, Int_t nSigma = 1)
+{
+ TH1* proj = hist->ProjectionY();
+ proj->Reset();
+
+ for (Int_t bin=1; bin<=proj->GetNbinsX(); bin++)
+ {
+ if (hist->Integral(1, hist->GetNbinsX(), bin, bin) == 0)
+ continue;
+
+ Int_t limit = -1;
+ for (limit = 1; limit<=hist->GetNbinsX(); limit++)
+ }
+}
+
+void ShowOnlyAccepted(TH2* input, const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
+{
+ TH2* proj = GetCorrection(fileName, dirName, ptmin);
+
+ for (Int_t y=1; y<=proj->GetNbinsY(); ++y)
+ for (Int_t x=1; x<=proj->GetNbinsX(); ++x)
+ if (proj->GetBinContent(x, y) > 5 || proj->GetBinContent(x, y) == 0)
+ {
+ proj->SetBinContent(x, y, 0);
+ }
+ else
+ proj->SetBinContent(x, y, 1);
+
+
+ input->Multiply(proj);
+}
+
+void MakeGaussianProfile(const char* histName = "fVertexCorrelation", Bool_t subtractMean = kFALSE)
+{
+ TFile::Open("correction_map.root");
+
+ TH2* hist2d = (TH2*) gFile->Get(histName);
+ hist2d->Sumw2();
+
+ TH1* result = hist2d->ProjectionX("result");
+ result->GetYaxis()->SetTitle(hist2d->GetYaxis()->GetTitle());
+ result->Reset();
+
+ for (Int_t x=1; x<hist2d->GetNbinsX(); ++x)
+ {
+ hist = hist2d->ProjectionY(Form("temp_%d", x), x, x);
+ if (hist->GetEntries() == 0)
+ continue;
+ if (hist->Fit("gaus") == 0)
+ {
+ func = hist->GetFunction("gaus");
+ mean = func->GetParameter(1);
+ error = func->GetParError(1);
+
+ if (subtractMean)
+ mean = hist2d->GetXaxis()->GetBinCenter(x) - mean;
+
+ result->SetBinContent(x, mean);
+ result->SetBinError(x, error);
+
+ if (x % 10 == 0)
+ {
+ new TCanvas;
+ ((TH1*) hist->Clone())->DrawCopy();
+ }
+ }
+ //break;
+ }
+
+ new TCanvas;
+ result->GetYaxis()->SetRangeUser(-0.2, 0.2);
+ result->Draw();
+}
+
+TH2* GetAcceptance(void* corr2d_void)
+{
+ corr2d = (AliCorrectionMatrix2D*) corr2d_void;
+ corr_xy = (TH2*) corr2d->GetCorrectionHistogram()->Clone("acceptance");
+
+ // fold in acceptance
+ for (Int_t x=1; x<=corr_xy->GetNbinsX(); ++x)
+ for (Int_t y=1; y<=corr_xy->GetNbinsY(); ++y)
+ {
+ if (corr_xy->GetBinContent(x, y) > 1.5)
+ corr_xy->SetBinContent(x, y, 0);
+
+ if (corr_xy->GetBinContent(x, y) > 0)
+ corr_xy->SetBinContent(x, y, 1);
+
+ corr_xy->SetBinError(x, y, 0);
+ }
+
+ return corr_xy;
+}
+
+void ZeroOutsideAcceptance(TH2* acc, TH3* hist)
+{
+ for (Int_t x=0; x<=acc->GetNbinsX()+1; ++x)
+ for (Int_t y=0; y<=acc->GetNbinsY()+1; ++y)
+ {
+ if (acc->GetBinContent(x, y) > 2 || acc->GetBinContent(x, y) == 0)
+ {
+ for (Int_t z=0; z<=hist->GetNbinsZ()+1; ++z)
+ {
+ hist->SetBinContent(x, y, z, 0);
+ hist->SetBinError(x, y, z, 0);
+ }
+ }
+ }
+}
+
+void DrawPhi()
+{
+ loadlibs();
+
+ TFile::Open("correction_map.root");
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+ if (!dNdEtaCorrection->LoadHistograms())
+ return 0;
+
+ TFile::Open("analysis_esd.root");
+ dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
+ fdNdEtaAnalysis->LoadHistograms();
+
+ // acc. map!
+ //acc = GetAcceptance(dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000));
+ acc = dNdEtaCorrection->GetCorrection(1)->GetTrackCorrection()->Get2DCorrection("yx", 0, 1000)->GetCorrectionHistogram();
+ //new TCanvas; acc->Draw("COLZ");
+
+ histG = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
+ ZeroOutsideAcceptance(acc, histG);
+ //new TCanvas; histG->Project3D("yx")->Draw("COLZ");
+ //histG->GetYaxis()->SetRangeUser(-0.9, 0.9);
+
+ histM = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
+ ZeroOutsideAcceptance(acc, histM);
+ //histM->GetYaxis()->SetRangeUser(-0.9, 0.9);
+
+ TFile::Open("analysis_mc.root");
+ dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndetaTrVtxMC", "dndetaTrVtxMC");
+ fdNdEtaAnalysis2->LoadHistograms("dndetaTrVtx");
+
+ histMC = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
+ ZeroOutsideAcceptance(acc, histMC);
+ //new TCanvas; histMC->Project3D("yx2")->Draw("COLZ");
+
+ //histG->GetZaxis()->SetRangeUser(1,2); histMC->GetZaxis()->SetRangeUser(1,2);
+ new TCanvas; a = histG->Project3D("yx3"); a->Add(histMC->Project3D("yx4"), -1); a->Draw("COLZ");
+
+ //histMC->GetYaxis()->SetRangeUser(-0.9, 0.9);
+
+ c = new TCanvas;
+
+ histG->GetXaxis()->SetRangeUser(-9.9, 9.9);
+ histG->Project3D("z")->DrawCopy();
+
+ histM->GetXaxis()->SetRangeUser(-9.9, 9.9);
+ proj = histM->Project3D("z2");
+ proj->SetLineColor(2);
+ proj->DrawCopy("SAME");
+
+ histMC->GetXaxis()->SetRangeUser(-9.9, 9.9);
+ projMC = histMC->Project3D("z3");
+ projMC->SetLineColor(4);
+ projMC->DrawCopy("SAME");
+}
+
+void PrintEventStats(Int_t corrID = 3, const char* fileName = "correction_map.root", const char* dir = "dndeta_correction")
+{
+ loadlibs();
+
+ /*
+ TFile::Open("analysis_mc.root");
+ fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD");
+ fdNdEtaAnalysis->LoadHistograms();
+ trackHist = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetGeneratedHistogram();
+ eventHist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetGeneratedHistogram();
+ */
+
+ TFile::Open(fileName);
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dir, dir);
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+ trackHist = dNdEtaCorrection->GetCorrection(corrID)->GetTrackCorrection()->GetGeneratedHistogram();
+ eventHist = dNdEtaCorrection->GetCorrection(corrID)->GetEventCorrection()->GetGeneratedHistogram();
+
+ trackHist->GetXaxis()->SetRange(0, trackHist->GetNbinsX()+1);
+ trackHist->GetZaxis()->SetRange(0, trackHist->GetNbinsZ()+1);
+ eta = trackHist->Project3D("y");
+
+ events = eventHist->Integral(0, eventHist->GetNbinsX()+1, 0, eventHist->GetNbinsY()+1);
+
+ eta->Scale(1.0 / events);
+
+ Float_t avgN = eta->Integral(0, eta->GetNbinsX()+1);
+ Printf("<N> = %f", avgN);
+
+ eta->Scale(1.0 / eta->GetXaxis()->GetBinWidth(1));
+
+ Printf("dndeta | eta = 0 is %f", (eta->GetBinContent(eta->FindBin(0.01)) + eta->GetBinContent(eta->FindBin(-0.01))) / 2);
+ Printf("dndeta in |eta| < 0.5 is %f", eta->Integral(eta->FindBin(-0.49), eta->FindBin(0.49)) / (eta->FindBin(0.49) - eta->FindBin(-0.49) + 1));
+ Printf("dndeta in |eta| < 1 is %f", eta->Integral(eta->FindBin(-0.99), eta->FindBin(0.99)) / (eta->FindBin(0.99) - eta->FindBin(-0.99) + 1));
+ Printf("dndeta in |eta| < 1.5 is %f", eta->Integral(eta->FindBin(-1.49), eta->FindBin(1.49)) / (eta->FindBin(1.49) - eta->FindBin(-1.49) + 1));
+
+ stats = (TH2*) gFile->Get("fEventStats");
+ proj = stats->ProjectionX();
+ gROOT->ProcessLine(".L PrintHist.C");
+ PrintHist2D(stats);
+ PrintHist(proj);
+
+ Float_t ua5_SD = 0.153;
+ Float_t ua5_DD = 0.080;
+ Float_t ua5_ND = 0.767;
+
+ Printf("+++ FRACTIONS +++");
+
+ Printf("ND: %f", proj->GetBinContent(3) / proj->GetBinContent(1));
+ Printf("SD: %f", proj->GetBinContent(4) / proj->GetBinContent(1));
+ Printf("DD: %f", proj->GetBinContent(5) / proj->GetBinContent(1));
+
+ Printf("+++ TRIGGER EFFICIENCIES +++");
+
+ Printf("INEL = %.1f", 100. * (proj->GetBinContent(1) - stats->GetBinContent(1, 1) - stats->GetBinContent(1, 3)) / proj->GetBinContent(1));
+ Printf("NSD = %.1f", 100. * (proj->GetBinContent(2) - stats->GetBinContent(2, 1) - stats->GetBinContent(2, 3)) / proj->GetBinContent(2));
+
+
+ Float_t trigND = 100. * (proj->GetBinContent(3) - stats->GetBinContent(3, 1) - stats->GetBinContent(3, 3)) / proj->GetBinContent(3);
+ Float_t trigSD = 100. * (proj->GetBinContent(4) - stats->GetBinContent(4, 1) - stats->GetBinContent(4, 3)) / proj->GetBinContent(4);
+ Float_t trigDD = 100. * (proj->GetBinContent(5) - stats->GetBinContent(5, 1) - stats->GetBinContent(5, 3)) / proj->GetBinContent(5);
+
+ Printf("ND = %.1f", trigND);
+ Printf("SD = %.1f", trigSD);
+ Printf("DD = %.1f", trigDD);
+
+ Float_t trigINELUA5 = ua5_SD * trigSD + ua5_DD * trigDD + ua5_ND * trigND;
+ Float_t trigNSDUA5 = (ua5_DD * trigDD + ua5_ND * trigND) / (ua5_DD + ua5_ND);
+ Printf("INEL (UA5) = %.1f", trigINELUA5);
+ Printf("NSD (UA5) = %.1f", trigNSDUA5);
+
+ Printf("+++ VERTEX EFFICIENCIES +++");
+
+ Printf("INEL = %.1f", 100. * (stats->GetBinContent(1, 3) + stats->GetBinContent(1, 4)) / proj->GetBinContent(1));
+ Printf("NSD = %.1f", 100. * (stats->GetBinContent(2, 3) + stats->GetBinContent(2, 4)) / proj->GetBinContent(2));
+
+ Float_t vtxND = 100. * (stats->GetBinContent(3, 3) + stats->GetBinContent(3, 4)) / proj->GetBinContent(3);
+ Float_t vtxSD = 100. * (stats->GetBinContent(4, 3) + stats->GetBinContent(4, 4)) / proj->GetBinContent(4);
+ Float_t vtxDD = 100. * (stats->GetBinContent(5, 3) + stats->GetBinContent(5, 4)) / proj->GetBinContent(5);
+ Printf("ND = %.1f", vtxND);
+ Printf("SD = %.1f", vtxSD);
+ Printf("DD = %.1f", vtxDD);
+
+ Float_t vtxINELUA5 = ua5_SD * vtxSD + ua5_DD * vtxDD + ua5_ND * vtxND;
+ Float_t vtxNSDUA5 = (ua5_DD * vtxDD + ua5_ND * vtxND) / (ua5_DD + ua5_ND);
+ Printf("INEL (UA5) = %.1f", vtxINELUA5);
+ Printf("NSD (UA5) = %.1f", vtxNSDUA5);
+
+ Printf("+++ TRIGGER + VERTEX EFFICIENCIES +++");
+
+ Printf("INEL = %.1f", 100. * stats->GetBinContent(1, 4) / proj->GetBinContent(1));
+ Printf("NSD = %.1f", 100. * stats->GetBinContent(2, 4) / proj->GetBinContent(2));
+ Printf("ND = %.1f", 100. * stats->GetBinContent(3, 4) / proj->GetBinContent(3));
+ Printf("SD = %.1f", 100. * stats->GetBinContent(4, 4) / proj->GetBinContent(4));
+ Printf("DD = %.1f", 100. * stats->GetBinContent(5, 4) / proj->GetBinContent(5));
+
+
+
+ for (Int_t i=7; i<=proj->GetNbinsX(); i++)
+ if (proj->GetBinContent(i) > 0)
+ Printf("bin %d (process type %d) = %.2f", i, (Int_t) proj->GetXaxis()->GetBinCenter(i), 100.0 * (proj->GetBinContent(i) - stats->GetBinContent(i, 1)) / proj->GetBinContent(i));
+
+ //eta->Draw();
+}
+
+void TestAsymmetry()
+{
+ loadlibs();
+
+ TFile* file2 = TFile::Open("analysis_mc.root");
+
+ dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+ fdNdEtaAnalysis->LoadHistograms();
+ fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
+
+ fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy();
+
+ hist = (TH1*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram();
+ hist2 = (TH1*) hist->Clone("hist2");
+ for (Int_t x=1; x<=hist->GetNbinsX(); x++)
+ for (Int_t y=1; y<=hist->GetNbinsY(); y++)
+ for (Int_t z=1; z<=hist->GetNbinsZ(); z++)
+ {
+ Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
+ hist->SetBinContent(x, y, z, hist2->GetBinContent(hist->GetNbinsX() / 2, TMath::Min(y, hist->GetNbinsY() + 1 - y), z));
+ }
+
+ hist = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
+ for (Int_t x=1; x<=hist->GetNbinsX(); x++)
+ for (Int_t y=1; y<=hist->GetNbinsY(); y++)
+ {
+ //Printf("%d %d %d %d", x, y, z, hist->GetNbinsY() + 1 - y);
+ hist->SetBinContent(x, y, hist->GetBinContent(hist->GetNbinsX() / 2, y));
+ }
+
+ fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone, "...");
+ fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerColor(2);
+ fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetLineColor(2);
+ fdNdEtaAnalysis->GetdNdEtaHistogram(0)->SetMarkerStyle(5);
+ fdNdEtaAnalysis->GetdNdEtaHistogram(0)->DrawCopy("SAMEP");
+}
+
+void DeltaPhiFromPt(Float_t smearing = 0.005)
+{
+ loadlibs();
+
+ TFile::Open("analysis_mc.root");
+ hist = (TH1*) gFile->Get("dndeta_check_pt");
+
+ dPhiHist = new TH1F("dPhiHist", ";#Delta phi", 400, -0.1, 0.1);
+ dPhiHist2 = new TH1F("dPhiHist2", ";#Delta phi", 400, -0.1, 0.1);
+
+ for (Int_t i=1; i<=hist->GetNbinsX(); i++)
+ {
+ Float_t pt = hist->GetBinCenter(i);
+ Float_t deltaPhi = (0.076 - 0.039) / 2 / (pt / 0.15);
+
+ if (smearing > 0)
+ {
+ gaus = new TF1("mygaus", "gaus(0)", -0.1, 0.1);
+ gaus->SetParameters(1, -deltaPhi, smearing);
+
+ dPhiHist->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
+
+ dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
+ gaus->SetParameters(1, deltaPhi, smearing);
+ dPhiHist2->FillRandom("mygaus", hist->GetBinContent(i) / 2 * 1000);
+ }
+ else
+{
+dPhiHist->Fill(deltaPhi, hist->GetBinContent(i) / 2);
+dPhiHist2->Fill(deltaPhi, hist->GetBinContent(i) / 2);
+dPhiHist2->Fill(-deltaPhi, hist->GetBinContent(i) / 2);
+}
+ }
+
+ new TCanvas;
+ dPhiHist->Draw();
+ dPhiHist2->SetLineColor(2);
+ dPhiHist2->Draw("SAME");
+ gPad->SetLogy();
+
+ TFile::Open("trackletsDePhi.root");
+ //TFile::Open("tmp/correction_maponly-positive.root");
+ //TFile::Open("tmp/correction_map.root");
+ //tracklets = (TH1*) gFile->Get(Form("fDeltaPhi_%d", 0));
+ tracklets = (TH1*) gFile->Get("DePhiPPTracklets");
+ tracklets->Scale(1.0 / tracklets->GetMaximum() * dPhiHist->GetMaximum());
+ tracklets->SetLineColor(4);
+ tracklets->Draw("SAME");
+}
+
+void VertexDistributions()
+{
+ loadlibs();
+
+ TFile::Open("correction_map.root");
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+
+ all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all");
+ trigger = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("trigger");
+ vtx = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram()->ProjectionX("vtx");
+
+ nottriggered = (TH1*) all->Clone("nottriggered");
+ nottriggered->Add(trigger, -1);
+
+ novertex = (TH1*) trigger->Clone("novertex");
+ novertex->Add(vtx, -1);
+
+ temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
+ highmult = temphist->ProjectionX("highmult", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
+ //all = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram()->ProjectionX("all", temphist->GetYaxis()->FindBin(10), temphist->GetNbinsY());
+
+ for (Int_t i=1; i<=trigger->GetNbinsX(); i++)
+ {
+ all->SetBinContent(i, all->GetBinContent(i) / all->GetBinWidth(i));
+ trigger->SetBinContent(i, trigger->GetBinContent(i) / trigger->GetBinWidth(i));
+ vtx->SetBinContent(i, vtx->GetBinContent(i) / vtx->GetBinWidth(i));
+ nottriggered->SetBinContent(i, nottriggered->GetBinContent(i) / nottriggered->GetBinWidth(i));
+ novertex->SetBinContent(i, novertex->GetBinContent(i) / novertex->GetBinWidth(i));
+ highmult->SetBinContent(i, highmult->GetBinContent(i) / highmult->GetBinWidth(i));
+ }
+
+ new TCanvas;
+ vtx->SetTitle("");
+ vtx->SetStats(0);
+ vtx->DrawCopy("HIST");
+
+ all->Scale(1.0 / all->Integral());
+ nottriggered->Scale(1.0 / nottriggered->Integral());
+ novertex->Scale(1.0 / novertex->Integral());
+ highmult->Scale(1.0 / highmult->Integral());
+
+ new TCanvas;
+ all->Draw("HIST");
+ novertex->SetLineColor(2);
+ novertex->Draw("HISTSAME");
+ highmult->SetLineColor(3);
+ highmult->Draw("HISTSAME");
+
+ legend = new TLegend(0.5, 0.5, 0.8, 0.8);
+ legend->SetFillColor(0);
+ legend->AddEntry(all, "all");
+ legend->AddEntry(novertex, "no vertex");
+ legend->AddEntry(highmult, "mult > 10");
+ legend->Draw();
+
+ new TCanvas;
+ trigger->Scale(1.0 / trigger->Integral());
+ vtx->Scale(1.0 / vtx->Integral());
+
+ trigger->Divide(vtx);
+
+ trigger->Draw();
+ //vtx->SetLineColor(2);
+ //vtx->Draw("SAME");
+
+ //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kVertexReco)->GetEventCorrection()->GetMeasuredHistogram();
+ temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetGeneratedHistogram();
+ //temphist = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kINEL)->GetEventCorrection()->GetMeasuredHistogram();
+
+ temphist = (TH2*) gFile->Get("fTemp1");
+
+ new TCanvas;
+ legend = new TLegend(0.7, 0.7, 0.99, 0.99);
+ legend->SetFillColor(0);
+
+ Bool_t first = kTRUE;
+ for (Int_t i=0; i<20; i+=5)
+ {
+ highmult = temphist->ProjectionX("highmult", i+1, i+1+4);
+ highmult->Rebin(10);
+ Printf("%f", highmult->Integral());
+ if (highmult->Integral() <= 0)
+ continue;
+
+ for (Int_t j=1; j<=trigger->GetNbinsX(); j++)
+ highmult->SetBinContent(j, highmult->GetBinContent(j) / highmult->GetBinWidth(j));
+
+ highmult->Scale(1.0 / highmult->Integral());
+ highmult->SetLineColor((i/5)+1);
+ highmult->GetYaxis()->SetRangeUser(0, 0.15);
+ if (first)
+ {
+ highmult->DrawCopy();
+ first = kFALSE;
+ }
+ else
+ highmult->DrawCopy("SAME");
+ legend->AddEntry(highmult->Clone(), Form("%d <= N <= %d", i, i+4));
+ }
+ legend->Draw();
+
+}
+
+void PlotPt1DCorrection()
+{
+ const char* files[] = { "field.root", "field_onlyprim.root", "nofield.root", "nofield_onlyprim.root" };
+ const char* names[] = { "B: all", "B: primaries", "No B: all", "No B: primaries" };
+ Int_t colors[] = { 1, 2, 3, 4 };
+
+ loadlibs();
+
+ dummy = new TH2F("dummy", ";p_{T};correction", 100, 0, 1.4, 100, 0.5, 3);
+ dummy->SetStats(0);
+ //dummy->GetYaxis()->SetTitleOffset(1.3);
+ dummy->Draw();
+
+ legend = new TLegend(0.48, 0.57, 0.88, 0.88);
+ legend->SetFillColor(0);
+
+ for (Int_t i=0; i<4; i++)
+ {
+ TFile::Open(files[i]);
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction");
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+
+ hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z", -9.9, 9.9, -0.79, 0.79);
+ hist->SetLineColor(colors[i]);
+ hist->SetLineWidth(2);
+ hist->SetMarkerColor(colors[i]);
+ hist->Draw("SAME");
+
+ legend->AddEntry(hist, names[i], "L");
+ }
+
+ legend->Draw();
+}
+
+void FitDiamond()
+{
+ TFile::Open("analysis_esd_raw.root");
+
+ hist = (TH3*) gFile->Get("vertex_check");
+
+ gStyle->SetOptFit(1);
+
+ TH1* proj[3];
+ proj[0] = hist->ProjectionX();
+ proj[1] = hist->ProjectionY();
+ proj[2] = hist->ProjectionZ();
+
+ for (Int_t i=0; i<3; i++)
+ {
+ c = new TCanvas;
+ proj[i]->Draw();
+ proj[i]->Fit("gaus");
+
+ c->SaveAs(Form("FitDiamond_%d.png", i));
+ }
+}
+
+void CompareDiamond(const char* mc, const char* data)
+{
+ TFile::Open(mc);
+
+ hist = (TH3*) gFile->Get("vertex_check");
+
+ gStyle->SetOptFit(1);
+
+ TH1* proj[3];
+ proj[0] = hist->ProjectionX("vertex_check_px");
+ proj[1] = hist->ProjectionY("vertex_check_py");
+ proj[2] = hist->ProjectionZ("vertex_check_pz");
+
+ TFile::Open(data);
+
+ hist = (TH3*) gFile->Get("vertex_check");
+
+ TH1* proj2[3];
+ proj2[0] = hist->ProjectionX("vertex_check_px2");
+ proj2[1] = hist->ProjectionY("vertex_check_py2");
+ proj2[2] = hist->ProjectionZ("vertex_check_pz2");
+
+ for (Int_t i=0; i<3; i++)
+ {
+ CompareQualityHists(proj[i], proj2[i], 1, 1);
+ }
+}
+
+void FitDiamondVsMult()
+{
+ TFile::Open("analysis_esd_raw.root");
+
+ fVertexVsMult = (TH3*) gFile->Get("fVertexVsMult");
+ fVertexVsMult->GetZaxis()->SetTitle("multiplicity");
+
+ TH2* proj[2];
+ proj[0] = (TH2*) fVertexVsMult->Project3D("xz");
+ proj[1] = (TH2*) fVertexVsMult->Project3D("yz");
+
+ gStyle->SetPadGridX(kTRUE);
+ gStyle->SetPadGridY(kTRUE);
+
+ Int_t max = 40;
+
+ for (Int_t i=0; i<2; i++)
+ {
+ proj[i]->Rebin2D(4, 1);
+ proj[i]->FitSlicesY();
+
+ c = new TCanvas(Form("c_%d", i), Form("c_%d", i), 800, 400);
+ c->Divide(2, 1);
+
+ c->cd(1);
+ hist = (TH1*) gROOT->FindObject(Form("fVertexVsMult_%sz_1", (i == 0) ? "x" : "y"));
+ hist->GetXaxis()->SetRangeUser(0, max);
+ hist->GetYaxis()->SetRangeUser(-0.4, 0.4);
+ hist->Draw();
+
+ c->cd(2);
+ hist = (TH1*) gROOT->FindObject(Form("fVertexVsMult_%sz_2", (i == 0) ? "x" : "y"));
+ hist->GetXaxis()->SetRangeUser(0, max);
+ hist->GetYaxis()->SetRangeUser(0, 0.2);
+ hist->Draw();
+
+ c->SaveAs(Form("FitDiamondVsMult_%d.png", i));
+ }
+}
+
+void CompareQualityHists(const char* fileName1, const char* fileName2, const char* plotName, Int_t rebin1 = 1, Int_t rebin2 = 1, const char* exec = 0)
+{
+ file1 = TFile::Open(fileName1);
+ hist1 = (TH1*) file1->Get(plotName);
+
+ file2 = TFile::Open(fileName2);
+ hist2 = (TH1*) file2->Get(plotName);
+
+ hist1->SetStats(0);
+
+ Printf("Entries in histograms: %f %f", hist1->Integral(), hist2->Integral());
+
+ if (exec)
+ {
+ hist1 = (TH1*) gROOT->ProcessLine(Form(exec, hist1, "hist1a"));
+ hist2 = (TH1*) gROOT->ProcessLine(Form(exec, hist2, "hist2a"));
+ hist1->Sumw2();
+ hist2->Sumw2();
+ Printf("Entries in histograms: %f %f", hist1->Integral(), hist2->Integral());
+ }
+
+ CompareQualityHists(hist1, hist2, rebin1, rebin2);
+}
+
+void CompareQualityHists(TH1* hist1, TH1* hist2, Int_t rebin1 = 1, Int_t rebin2 = 1)
+{
+ hist1->SetLineColor(1);
+ hist2->SetLineColor(2);
+
+ if (rebin1 != 0 && rebin2 != 0)
+ {
+ hist1->Rebin(TMath::Abs(rebin1));
+ hist2->Rebin(TMath::Abs(rebin2));
+ }
+
+ //hist2 = hist2->Rebin(hist1->GetNbinsX(), Form("%s_rebinned", hist2->GetName()), hist1->GetXaxis()->GetXbins()->GetArray());
+
+ //hist1->Scale(1.0 / 0.83);
+
+//hist1->GetXaxis()->SetRangeUser(0, 50);
+/* hist1->GetYaxis()->SetRangeUser(0.9, 1.2);
+ hist1->Scale(1.0 / 0.808751);*/
+
+ //hist1->Scale(1.0 / 1.24632);
+ //hist1->Scale(1.0 / 1.23821);
+ //hist1->Scale(1.0 / 1.26213);
+
+ if (rebin1 > 0 && rebin2 > 0)
+ {
+ hist1->Scale(hist2->Integral() / hist1->Integral() * hist2->GetXaxis()->GetBinWidth(1) / hist1->GetXaxis()->GetBinWidth(1) / rebin1 * rebin2);
+
+ //hist1->Scale(0.5);
+ //hist2->Scale(0.5);
+ }
+
+ c = new TCanvas;
+ if (strcmp(hist1->GetName(), "fDeltaTheta") == 0 || strcmp(hist1->GetName(), "fDeltaPhi") == 0 || strcmp(hist1->GetName(), "fMultVtx") == 0 || TString(hist1->GetName()).BeginsWith("vertex_check"))
+ c->SetLogy();
+
+ if (TString(hist1->GetName()).BeginsWith("fMultiplicityESD"))
+ {
+ c->SetLogy();
+ loadlibs();
+ AliPWG0Helper::NormalizeToBinWidth(hist1);
+ AliPWG0Helper::NormalizeToBinWidth(hist2);
+ }
+
+ Printf("Means: %f %f %e", hist1->GetMean(), hist2->GetMean(), 1.0 - hist2->GetMean() / hist1->GetMean());
+
+ //hist1->GetYaxis()->SetRangeUser(0.01, hist1->GetMaximum() * 1.3);
+ hist1->DrawCopy("HISTE");
+ hist2->DrawCopy("HISTE SAME");
+ gPad->SetGridx();
+ gPad->SetGridy();
+ //gPad->SetLogy();
+ c->SaveAs(Form("%s_1.png", hist1->GetName()));
+
+ if (rebin1 == rebin2)
+ {
+ for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
+ if (hist1->GetBinContent(i) == 0 && hist2->GetBinContent(i) > 0 || hist1->GetBinContent(i) > 0 && hist2->GetBinContent(i) == 0)
+ Printf("Inconsistent bin %d: %f %f", i, hist1->GetBinContent(i), hist2->GetBinContent(i));
+
+ c2 = new TCanvas;
+ hist1->GetYaxis()->SetRangeUser(0.5, 1.5);
+ hist1->Divide(hist2);
+ hist1->DrawCopy("HIST");
+ gPad->SetGridx();
+ gPad->SetGridy();
+ c2->SaveAs(Form("%s_2.png", hist1->GetName()));
+
+ /*
+ for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
+ if (hist1->GetBinContent(i) > 0.9 && hist1->GetBinContent(i) < 1.1)
+ hist1->SetBinContent(i, 0);
+
+ new TCanvas;
+ hist1->SetMarkerStyle(20);
+ hist1->DrawCopy("P");
+ */
+ }
+}
+
+void DrawClustersVsTracklets()
+{
+ TFile::Open("analysis_esd_raw.root");
+
+ hist = (TH2*) gFile->Get("fTrackletsVsClusters");
+
+ c = new TCanvas("c", "c", 600, 600);
+ c->SetRightMargin(0.05);
+ c->SetTopMargin(0.05);
+
+ hist->SetStats(0);
+ hist->GetYaxis()->SetRangeUser(0, 400);
+ hist->GetYaxis()->SetTitleOffset(1.3);
+ hist->GetXaxis()->SetRangeUser(0, 30);
+ hist->Draw("BOX");
+
+ func = new TF1("func", "80 + x * 11", 0, 30);
+ func->Draw("SAME");
+
+ c->SaveAs("clusters_vs_tracklets.eps");
+}
+
+void VertexPlotBackgroundNote()
+{
+ TFile::Open("all.root");
+
+ hist = (TH3*) gFile->Get("vertex_check");
+ proj = (TH1*) hist->ProjectionZ()->Clone("all");
+ proj->Rebin(2);
+
+ proj->Draw();
+
+ TFile::Open("analysis_esd_raw.root");
+ hist = (TH3*) gFile->Get("vertex_check");
+ proj = (TH1*) hist->ProjectionZ()->Clone("afterbg");
+ proj->Rebin(2);
+
+ proj->SetLineColor(2);
+ proj->Draw("SAME");
+}
+
+void BackgroundAnalysis(const char* signal, const char* background)
+{
+ TFile::Open(signal);
+ signalHist = (TH2*) gFile->Get("fTrackletsVsClusters");
+
+ TFile::Open(background);
+ backgroundHist = (TH2*) gFile->Get("fTrackletsVsClusters");
+
+ Printf("For events with >= 1 tracklet:");
+
+ func = new TF1("func", "[0] + x * 11", 0, 30);
+ for (Int_t a = 50; a <= 100; a += 10)
+ {
+ func->SetParameter(0, a);
+
+ Float_t signalCount = 0;
+ Float_t backgroundCount = 0;
+ for (Int_t x = 2; x <= signalHist->GetNbinsX(); x++)
+ {
+ signalCount += signalHist->Integral(x, x, signalHist->GetYaxis()->FindBin(func->Eval(signalHist->GetXaxis()->GetBinCenter(x))), signalHist->GetNbinsY());
+ backgroundCount += backgroundHist->Integral(x, x, signalHist->GetYaxis()->FindBin(func->Eval(signalHist->GetXaxis()->GetBinCenter(x))), signalHist->GetNbinsY());
+ }
+
+ Float_t signalFraction = 100.0 * signalCount / signalHist->Integral(2, signalHist->GetNbinsX(), 1, signalHist->GetNbinsY());
+ Float_t backgroundFraction = 100.0 * backgroundCount / backgroundHist->Integral(2, signalHist->GetNbinsX(), 1, signalHist->GetNbinsY());
+
+ Printf("Cut at a = %d; Removed %.2f %% of the background (%.0f events); Removed %.2f %% of the signal", a, backgroundFraction, backgroundCount, signalFraction);
+ }
+}
+
+void ZPhiPlots()
+{
+ TFile::Open("analysis_esd_raw.root");
+
+ for (Int_t i=0; i<2; i++)
+ {
+ hist = (TH2*) gFile->Get(Form("fZPhi_%d", i));
+
+ c = new TCanvas;
+ hist->SetStats(0);
+ hist->Draw("COLZ");
+ c->SaveAs(Form("ZPhi_%d.png", i));
+ }
+}
+
+void DrawStats(Bool_t all = kFALSE)
+{
+ if (all)
+ {
+ Int_t count = 4;
+ const char* list[] = { "CINT1B-ABCE-NOPF-ALL/spd", "CINT1A-ABCE-NOPF-ALL/spd", "CINT1C-ABCE-NOPF-ALL/spd", "CINT1-E-NOPF-ALL/spd" };
+ }
+ else
+ {
+ Int_t count = 1;
+ const char* list[] = { "." };
+ }
+
+ for (Int_t i=0; i<count; i++)
+ {
+ TFile::Open(Form("%s/analysis_esd_raw.root", list[i]));
+
+ hist = (TH2*) gFile->Get("fStats2");
+
+ c = new TCanvas(list[i], list[i], 800, 600);
+ gPad->SetBottomMargin(0.2);
+ gPad->SetLeftMargin(0.2);
+ gPad->SetRightMargin(0.2);
+ hist->Draw("TEXT");
+ hist->SetMarkerSize(2);
+ //hist->GetYaxis()->SetRangeUser(0, 0);
+
+ gROOT->Macro("increaseFonts.C");
+
+ c->SaveAs(Form("%s/stats.png", list[i]));
+ }
+}
+
+void CompareMCDataTrigger(const char* mcFile, const char* dataFile)
+{
+ TH1* stat[2];
+
+ TFile::Open(mcFile);
+ mc = (TH1*) gFile->Get("trigger_histograms_/fHistFiredBitsSPD");
+ stat[0] = (TH1*) gFile->Get("fHistStatistics");
+
+ TFile::Open(dataFile);
+ data = (TH1*) gFile->Get("trigger_histograms_+CINT1B-ABCE-NOPF-ALL/fHistFiredBitsSPD");
+ if (!data)
+ data = (TH1*) gFile->Get("trigger_histograms_+CSMBB-ABCE-NOPF-ALL/fHistFiredBitsSPD");
+
+ stat[1] = (TH1*) gFile->Get("fHistStatistics");
+
+ CompareQualityHists(mc, data);
+
+ for (Int_t i=0; i<2; i++)
+ {
+ Float_t total = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("Trigger class"), 1);
+ Float_t spd = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("FO >= 2"), 1);
+ Float_t v0A = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("V0A"), 1);
+ Float_t v0C = stat[i]->GetBinContent(stat[i]->GetXaxis()->FindBin("V0C"), 1);
+
+ Printf("%s:\nSPD / V0A: %.3f\nSPD / V0C: %.3f\nV0A / V0C: %.3f", (i == 0) ? "MC " : "Data", spd / v0A, spd / v0C, v0A / v0C);
+ Printf("SPD / Total: %.3f\nV0A / Total: %.3f\nV0C / Total: %.3f\n", spd / total, v0A / total, v0C / total);
+ }
+}
+
+void CompareMCDatadNdEta(const char* mcFile, const char* dataFile)
+{
+ //CompareQualityHists(mcFile, dataFile, "fEtaPhi", 4, 4, "((TH2*)%p)->ProjectionY(\"%s\", 1, 40)");
+ //CompareQualityHists(mcFile, dataFile, "fEtaPhi", 4, 4, "((TH2*)%p)->ProjectionY(\"%s\", 41, 80)");
+
+ CompareQualityHists(mcFile, dataFile, "fEtaPhi", 1, 1, "((TH2*)%p)->ProjectionX(\"%s\", 271, 360)");
+}
+
+void TrigVsTrigVtx(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
+{
+ loadlibs();
+ if (!TFile::Open(fileName))
+ return;
+
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+
+ TH2* eTrig = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
+ TH2* eTrigVtx = dNdEtaCorrection->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
+
+ eTrig_proj = eTrig->ProjectionY("y1", eTrig->GetYaxis()->FindBin(-9.9), eTrig->GetYaxis()->FindBin(9.9));
+ eTrigVtx_proj = eTrigVtx->ProjectionY("y2", eTrig->GetYaxis()->FindBin(-9.9), eTrig->GetYaxis()->FindBin(9.9));
+
+ new TCanvas;
+ eTrig_proj->Draw();
+ eTrig_proj->GetXaxis()->SetRangeUser(0, 20);
+ eTrigVtx_proj->SetLineColor(2);
+ eTrigVtx_proj->Draw("SAME");
+
+ gPad->SetLogy();
+}
+
+void PrintAverageNSDCorrectionFactors()
+{
+ // factors estimated from MC, can be slighly different with data b/c correction is applies as function of measured multiplicity
+
+ loadlibs();
+
+ if (!TFile::Open("correction_mapprocess-types.root"))
+ return;
+
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND");
+ if (!dNdEtaCorrection->LoadHistograms())
+ return;
+
+ AlidNdEtaCorrection* dNdEtaCorrection2 = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD");
+ if (!dNdEtaCorrection2->LoadHistograms())
+ return;
+
+ // for scaling factors see drawSystematics.C; GetRelativeFractions()
+ // 900 GeV
+ //dNdEtaCorrection->Scale(1.06);
+ //dNdEtaCorrection->Add(dNdEtaCorrection2, 9.5 / 12.3);
+ // 2.36 TeV
+ dNdEtaCorrection->Scale(1.036);
+ dNdEtaCorrection->Add(dNdEtaCorrection2, 0.075 * 1.43 / 0.127);
+
+ Printf("event adding: %f", dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral());
+
+ Printf("track adding: %f", dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetGeneratedHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral());
+
+ AlidNdEtaCorrection* dNdEtaCorrection3 = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD");
+ if (!dNdEtaCorrection3->LoadHistograms())
+ return;
+
+ // 900 GeV
+ //dNdEtaCorrection3->Scale(0.153 / 0.189);
+ // 2.36 TeV
+ dNdEtaCorrection3->Scale(0.159 / 0.166);
+ dNdEtaCorrection->Add(dNdEtaCorrection3);
+
+ Printf("event subtraction: %f", dNdEtaCorrection3->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral());
+
+ Printf("track subtraction: %f", dNdEtaCorrection3->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrection->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetMeasuredHistogram()->Integral());
+
+ dNdEtaCorrection->GetTriggerBiasCorrectionNSD()->PrintInfo(0.0);
+}
+