Merge branch 'feature-movesplit'

[u/mrichter/AliRoot.git] / PWGUD / selectors / dNdEta / drawSystematics.C

/* $Id$ */

#if !defined(__CINT__) || defined(__MAKECINT__)

#include "AliPWG0Helper.h"
#include "dNdEtaAnalysis.h"
#include "AlidNdEtaCorrection.h"

#include <TCanvas.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TLine.h>
#include <TSystem.h>
#include <TLegend.h>
#include <TPad.h>
#include <TF1.h>

extern TPad* gPad;

void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10);

#endif

Int_t markers[] = {20,20,21,22,23,28,29};
Int_t colors[]  = {1,2,3,4,6,8,102};

void loadlibs()
{
  gSystem->Load("libTree");
  gSystem->Load("libVMC");

  gSystem->Load("libSTEERBase");
  gSystem->Load("libANALYSIS");
  gSystem->Load("libPWG0base");
}

void SetRanges(TAxis* axis)
{
  if (strcmp(axis->GetTitle(), "#eta") == 0)
    axis->SetRangeUser(-1.7999, 1.7999);
  if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
    axis->SetRangeUser(0, 9.9999);
  if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
    axis->SetRangeUser(-15, 14.9999);
  if (strcmp(axis->GetTitle(), "Ntracks") == 0)
    axis->SetRangeUser(0, 99.9999);
}

void SetRanges(TH1* hist)
{
  SetRanges(hist->GetXaxis());
  SetRanges(hist->GetYaxis());
  SetRanges(hist->GetZaxis());
}

void Prepare3DPlot(TH3* hist)
{
  hist->GetXaxis()->SetTitleOffset(1.5);
  hist->GetYaxis()->SetTitleOffset(1.5);
  hist->GetZaxis()->SetTitleOffset(1.5);

  hist->SetStats(kFALSE);
}

void Prepare2DPlot(TH2* hist)
{
  hist->SetStats(kFALSE);
  hist->GetYaxis()->SetTitleOffset(1.4);

  SetRanges(hist);
}

void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
{
  hist->SetLineWidth(2);
  hist->SetStats(kFALSE);

  hist->GetXaxis()->SetTitleOffset(1.2);
  hist->GetYaxis()->SetTitleOffset(1.2);

  if (setRanges)
    SetRanges(hist);
}

void InitPad()
{
  if (!gPad)
    return;

  gPad->Range(0, 0, 1, 1);
  gPad->SetLeftMargin(0.15);
  //gPad->SetRightMargin(0.05);
  //gPad->SetTopMargin(0.13);
  //gPad->SetBottomMargin(0.1);

  gPad->SetGridx();
  gPad->SetGridy();
}

void InitPadCOLZ()
{
  gPad->Range(0, 0, 1, 1);
  gPad->SetRightMargin(0.15);
  gPad->SetLeftMargin(0.12);

  gPad->SetGridx();
  gPad->SetGridy();
}

void Secondaries()
{
  TFile* file = TFile::Open("systematics.root");

  TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
  if (!secondaries)
  {
    printf("Could not read histogram\n");
    return;
  }

  TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
  canvas->Divide(3, 3);
  for (Int_t i=1; i<=8; i++)
  {
    TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
    hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));

    canvas->cd(i);
    hist->Draw();
  }
}

void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
{
  gSystem->Load("libPWG0base");

  TString canvasName;
  canvasName.Form("Track2Particle1DComposition");
  TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
  canvas->Divide(3, 1);

  TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);

  for (Int_t i=0; i<folderCount; ++i)
  {
    Correction1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);

    TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
    TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
    TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));

    Prepare1DPlot(corrX);
    Prepare1DPlot(corrY);
    Prepare1DPlot(corrZ);

    const char* title = "Track2Particle Correction";
    corrX->SetTitle(title);
    corrY->SetTitle(title);
    corrZ->SetTitle(title);

    corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);

    corrX->SetLineColor(i+1);
    corrY->SetLineColor(i+1);
    corrZ->SetLineColor(i+1);

    canvas->cd(1);
    InitPad();
    corrX->DrawCopy(((i>0) ? "SAME" : ""));

    canvas->cd(2);
    InitPad();
    corrY->DrawCopy(((i>0) ? "SAME" : ""));

    canvas->cd(3);
    InitPad();
    corrZ->DrawCopy(((i>0) ? "SAME" : ""));

    legend->AddEntry(corrZ, folderNames[i]);
  }

  legend->Draw();

  //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
  //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
}

TH1** DrawRatios(const char* fileName = "systematics.root")
{
  gSystem->Load("libPWG0base");

  TFile* file = TFile::Open(fileName);

  TString canvasName;
  canvasName.Form("DrawRatios");
  TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
  canvas->Divide(2, 1);
  canvas->cd(1);

  TH1** ptDists = new TH1*[4];

  TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);

  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
  const char* particleNames[] = { "#pi", "K", "p", "other" };
  for (Int_t i=0; i<4; ++i)
  {
    AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]);
    dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]);

    TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();

    gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
    gene->GetXaxis()->SetRangeUser(-10, 10);

    ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i])));
    ptDists[i]->SetTitle(";p_{T};Count");
    if (!ptDists[i])
    {
      printf("Problem getting distribution %d.\n", i);
      return 0;
    }

    ptDists[i]->GetYaxis()->SetRangeUser(1, ptDists[i]->GetMaximum()*1.1);
    ptDists[i]->GetXaxis()->SetRangeUser(0, 9.9);
    ptDists[i]->SetLineColor(i+1);
    ptDists[i]->DrawCopy((i == 0) ? "" : "SAME");
    ptDists[i]->GetYaxis()->SetRange(0, 0);

    legend->AddEntry(ptDists[i], particleNames[i]);
  }
  gPad->SetLogy();

  TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total"));

  for (Int_t i=1; i<4; ++i)
    total->Add(ptDists[i]);

  canvas->cd(2);
  for (Int_t i=0; i<4; ++i)
  {
    ptDists[i]->Divide(total);
    ptDists[i]->SetStats(kFALSE);
    ptDists[i]->SetTitle(";p_{T};Fraction of total");
    ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
    ptDists[i]->Draw((i == 0) ? "" : "SAME");
  }
  legend->SetFillColor(0);
  legend->Draw();

  canvas->SaveAs("DrawRatios.gif");


  canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500);
  for (Int_t i=0; i<4; ++i)
  {
    TH1* hist = ptDists[i]->Clone();
    hist->GetXaxis()->SetRangeUser(0, 1.9);
    hist->Draw((i == 0) ? "" : "SAME");
  }
  legend->Draw();

  canvas->SaveAs("pythiaratios.eps");

  file->Close();

  return ptDists;
}

void DrawCompareToReal()
{
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = { "correction_map.root", "new_compositions.root" };
  const char* folderNames[] = { "dndeta_correction", "PythiaRatios" };

  Track2Particle1DComposition(fileNames, 2, folderNames);
}

void DrawDifferentSpecies()
{
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };

  Track2Particle1DComposition(fileNames, 4, folderNames);
}

void DrawpiKpAndCombined()
{
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };

  Track2Particle1DComposition(fileNames, 4, folderNames);
}

void DrawSpeciesAndCombination()
{
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" };
  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" };

  Track2Particle1DComposition(fileNames, 4, folderNames);
}

void DrawSimulatedVsCombined()
{
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = { "new_compositions.root", "new_compositions.root" };
  const char* folderNames[] = { "Pythia", "PythiaRatios" };

  Track2Particle1DComposition(fileNames, 2, folderNames);
}

void DrawBoosts()
{
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
  const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" };

  Track2Particle1DComposition(fileNames, 4, folderNames);
}

TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
{
  gSystem->Load("libPWG0base");

  AlidNdEtaCorrection* fdNdEtaCorrection[2];

  TFile::Open(filename);

  fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1);
  fdNdEtaCorrection[0]->LoadHistograms(filename, folder1);

  fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2);
  fdNdEtaCorrection[1]->LoadHistograms(filename, folder2);

  TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
  TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();

  //TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s / %s;Count", folder2, folder1), 1000, 0.9, 1.1);
  TH2F* difference = new TH2F(Form("difference_%s_%s", folder1, folder2), Form(";#Sigma (C_{pT, z} %s / C_{pT, z} %s);#eta;Count", folder2, folder1), 100, 0.9, 1.1, hist1->GetYaxis()->GetNbins(), hist1->GetYaxis()->GetXmin(), hist1->GetYaxis()->GetXmax());

  for (Int_t x=hist1->GetXaxis()->FindBin(-10); x<=hist1->GetXaxis()->FindBin(10); ++x)
    for (Int_t y=hist1->GetYaxis()->FindBin(-0.8); y<=hist1->GetYaxis()->FindBin(0.8); ++y)
      for (Int_t z=hist1->GetZaxis()->FindBin(0.3); z<=hist1->GetZaxis()->FindBin(9.9); ++z)
        if (hist1->GetBinContent(x, y, z) != 0)
          difference->Fill(hist2->GetBinContent(x, y, z) / hist1->GetBinContent(x, y, z), hist1->GetYaxis()->GetBinCenter(y));

  difference->GetYaxis()->SetRangeUser(-0.8, 0.8);

  printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));

  return difference;
}

void HistogramDifferences()
{
  TH2F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted");
  TH2F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted");
  TH2F* KReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "KReduced");
  TH2F* pReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "pReduced");

  TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
  canvas->Divide(2, 2);

  canvas->cd(1);
  KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
  KBoosted->Draw("COLZ");

  canvas->cd(2);
  KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
  KReduced->Draw("COLZ");

  canvas->cd(3);
  pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02);
  pBoosted->Draw("COLZ");

  canvas->cd(4);
  pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02);
  pReduced->Draw("COLZ");

  canvas->SaveAs("HistogramDifferences.gif");

  canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500);
  canvas->Divide(2, 1);

  canvas->cd(1);
  gPad->SetBottomMargin(0.13);
  KBoosted->SetTitle("a) Ratio of correction maps");
  KBoosted->SetStats(kFALSE);
  KBoosted->GetXaxis()->SetTitleOffset(1.4);
  KBoosted->GetXaxis()->SetLabelOffset(0.02);
  KBoosted->Draw("COLZ");

  canvas->cd(2);
  gPad->SetGridx();
  gPad->SetGridy();

  TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);

  for (Int_t i=0; i<4; ++i)
  {
    TH2F* hist = 0;
    TString name;
    switch (i)
    {
      case 0: hist = KBoosted; name = "K enhanced"; break;
      case 1: hist = KReduced; name = "K reduced"; break;
      case 2: hist = pBoosted; name = "p enhanced"; break;
      case 3: hist = pReduced; name = "p reduced"; break;
    }

    TProfile* profile = hist->ProfileY();
    profile->SetTitle("b) Mean and RMS");
    profile->GetXaxis()->SetRange(hist->GetYaxis()->GetFirst(), hist->GetYaxis()->GetLast());
    profile->GetXaxis()->SetTitleOffset(1.2);
    profile->GetXaxis()->SetLabelOffset(0.02);
    profile->GetYaxis()->SetRangeUser(0.98, 1.02);
    profile->SetStats(kFALSE);
    profile->SetLineColor(i+1);
    profile->SetMarkerColor(i+1);
    profile->DrawCopy(((i > 0) ? "SAME" : ""));


    legend->AddEntry(profile, name);
  }

  legend->Draw();
  canvas->SaveAs("particlecomposition_result.eps");
}


void ScalePtDependent(TH3F* hist, TF1* function)
{
  // assumes that pt is the third dimension of hist
  // scales with function(pt)

  for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
  {
    Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z));
    printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);

    for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
      for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
        hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
  }
}

void ScalePtDependent(TH3F* hist, TH1* function)
{
  // assumes that pt is the third dimension of hist
  // scales with histogram(pt)

  for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
  {
    Double_t factor = function->GetBinContent(function->GetXaxis()->FindBin(hist->GetZaxis()->GetBinCenter(z)));
    printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);

    for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
      for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
        hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
  }
}

const char* ChangeComposition(void** correctionPointer, Int_t index)
{
  AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;

  switch (index)
  {
    case 0: // result from pp events
      {
        TFile::Open("pythiaratios.root");

        for (Int_t i=0; i<4; ++i)
        {
          TString name;
          name.Form("correction_%d", i);
          fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
          fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name);
        }
      }
      return "Pythia";
      break;

    case 1: // each species rated with pythia ratios
      /*TH1** ptDists = DrawRatios("pythiaratios.root");

      for (Int_t i=0; i<3; ++i)
      {
        ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
        ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
      }*/
      return "PythiaRatios";
      break;

      // one species enhanced / reduced
    case 2: // + 30% kaons
    case 3: // - 30% kaons
    case 4: // + 30% protons
    case 5: // - 30% protons
    case 6: // + 30% kaons + 30% protons
    case 7: // - 30% kaons - 30% protons
    case 8: // + 30% kaons - 30% protons
    case 9: // - 30% kaons + 30% protons
    case 10: // + 30% others
    case 11: // - 30% others
      TString* str = new TString;
      if (index < 6)
      {
        Int_t correctionIndex = index / 2;
        Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
  
        fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
        str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : (correctionIndex == 2) ? "p" : "others"), (index % 2 == 0) ? "Boosted" : "Reduced");
      }
      else if (index < 10)
      {
        Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
        fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
        str->Form("%s%s", "K", (scaleFactor > 1) ? "Boosted" : "Reduced");
        
        if (index >= 8)
          scaleFactor = (index % 2 == 0) ? 0.3 : 1.7;
        fdNdEtaCorrection[2]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
        *str += Form("%s%s", "p", (scaleFactor > 1) ? "Boosted" : "Reduced");
      }
      else
      {
        Double_t scaleFactor = (index % 2 == 0) ? 1.3 : 0.7;
        fdNdEtaCorrection[3]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
        str->Form("%s%s", "others", (scaleFactor > 1) ? "Boosted" : "Reduced");
      }

      return str->Data();
      break;

      // each species rated with pythia ratios
    case 12: // + 50% pions
    case 13: // - 50% pions
    case 14: // + 50% kaons
    case 15: // - 50% kaons
    case 16: // + 50% protons
    case 17: // - 50% protons
      TH1** ptDists = DrawRatios("pythiaratios.root");
      Int_t functionIndex = (index - 2) / 2;
      Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
      ptDists[functionIndex]->Scale(scaleFactor);

      for (Int_t i=0; i<3; ++i)
      {
        ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDists[i]);
        ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDists[i]);
      }
      TString* str = new TString;
      str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
      return str->Data();
      break;

    case 999:
      TF1* ptDependence = new TF1("simple", "x", 0, 100);
      ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDependence);
      ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDependence);
      break;

  }

  return "noname";
}

void Composition()
{
  loadlibs();

  gSystem->Unlink("new_compositions.root");
  gSystem->Unlink("new_compositions_analysis.root");
  
  const char* names[] = { "pi", "K", "p", "other" };
  TH1* hRatios[20];

  //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
  backgroundEvents = -1;    // use 0 bin from MC! for 2.36 TeV
  
  Printf("Subtracting %d background events!!!", backgroundEvents);
  gSystem->Sleep(1000);
  
  Int_t nCompositions = 12;
  Int_t counter = 0;
  for (Int_t comp = 1; comp < nCompositions; ++comp)
  {
    AlidNdEtaCorrection* fdNdEtaCorrection[4];

    TFile::Open("correction_mapparticle-species.root");

    for (Int_t i=0; i<4; ++i)
    {
      TString name;
      name.Form("dndeta_correction_%s", names[i]);
      fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
      fdNdEtaCorrection[i]->LoadHistograms();
    }

    const char* newName = ChangeComposition(fdNdEtaCorrection, comp);

    Double_t geneCount[5];
    Double_t measCount[5];
    geneCount[4] = 0;
    measCount[4] = 0;

    for (Int_t i=0; i<4; ++i)
    {
      geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Integral();
      measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Integral();

      geneCount[4] += geneCount[i];
      measCount[4] += measCount[i];

      printf("Particle %s: %d gene, %d meas\n", ((i == 0) ? "pi" : (i == 1) ? "K" : (i == 2) ? "p" : "others"), (Int_t) geneCount[i], (Int_t) measCount[i]);
    }

    printf("Generated ratios are:     %f pi, %f K, %f p, %f others\n", geneCount[0] / geneCount[4], geneCount[1] / geneCount[4], geneCount[2] / geneCount[4], geneCount[3] / geneCount[4]);

    printf("Reconstructed ratios are: %f pi, %f K, %f p, %f others\n", measCount[0] / measCount[4], measCount[1] / measCount[4], measCount[2] / measCount[4], measCount[3] / measCount[4]);

    TList* collection = new TList;

    // skip "other" particle correction here
    // with them has to be dealt differently, maybe just increasing the neutral particles...
    for (Int_t i=1; i<4; ++i)
      collection->Add(fdNdEtaCorrection[i]);

    fdNdEtaCorrection[0]->Merge(collection);
    fdNdEtaCorrection[0]->Finish();

    delete collection;

    // save everything
    TFile* file = TFile::Open("new_compositions.root", "UPDATE");
    fdNdEtaCorrection[0]->SetName(newName);
    fdNdEtaCorrection[0]->SaveHistograms();
    //file->Write();
    file->Close();
    
    // correct dNdeta distribution with modified correction map
    TFile::Open("analysis_esd_raw.root");

    dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
    fdNdEtaAnalysis->LoadHistograms();

    fdNdEtaAnalysis->Finish(fdNdEtaCorrection[0], 0.2, 3, newName);
    
    hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(newName);
    hRatios[counter]->SetTitle(newName);
    hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default composition}{modified composition}");

    if (counter > 0)
      hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);

    file = TFile::Open("new_compositions_analysis.root", "UPDATE");
    hRatios[counter]->Write();
    file->Close();
    
    delete fdNdEtaAnalysis;

    counter++;
  }

  /*
  gROOT->ProcessLine(".L drawPlots.C");

  const char* fileNames[] = {"new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
  const char* folderNames[] = { "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };

  Track2Particle1DComposition(fileNames, nCompositions, folderNames);
  */
}


void drawSystematics()
{
  //Secondaries();
  //DrawDifferentSpecies();
  //Composition();

  Sigma2VertexSimulation();

}

void DrawdNdEtaDifferences()
{
  TH1* hists[5];

  TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
  legend->SetFillColor(0);

  TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
  canvas->Divide(2, 1);

  canvas->cd(1);

  for (Int_t i=0; i<5; ++i)
  {
    hists[i] = 0;
    TFile* file = 0;
    TString title;

    switch(i)
    {
      case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
      case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
      case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
      case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
      case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
      default: return;
    }

    if (file)
    {
      hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
      hists[i]->SetTitle("a)");

      Prepare1DPlot(hists[i], kFALSE);
      hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
      hists[i]->GetYaxis()->SetRangeUser(6, 7);
      hists[i]->SetLineColor(colors[i]);
      hists[i]->SetMarkerColor(colors[i]);
      hists[i]->SetMarkerStyle(markers[i]);
      hists[i]->GetXaxis()->SetLabelOffset(0.015);
      hists[i]->GetYaxis()->SetTitleOffset(1.5);
      gPad->SetLeftMargin(0.12);
      hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));

      legend->AddEntry(hists[i], title);
      hists[i]->SetTitle(title);
    }
  }
  legend->Draw();

  canvas->cd(2);
  gPad->SetLeftMargin(0.14);

  TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
  legend2->SetFillColor(0);

  for (Int_t i=1; i<5; ++i)
  {
    if (hists[i])
    {
      legend2->AddEntry(hists[i]);

      hists[i]->Divide(hists[0]);
      hists[i]->SetTitle("b)");
      hists[i]->SetLineColor(colors[i-1]);
      hists[i]->SetMarkerColor(colors[i-1]);
      hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
      hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
      hists[i]->GetYaxis()->SetTitleOffset(1.8);
      hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
    }
  }

  legend2->Draw();

  canvas->SaveAs("particlecomposition_result_detail.gif");

  TCanvas* canvas2 = new TCanvas("DrawdNdEtaDifferences2", "DrawdNdEtaDifferences2", 700, 500);

  for (Int_t i=1; i<5; ++i)
  {
    if (hists[i])
    {
      hists[i]->SetTitle("");
      hists[i]->GetYaxis()->SetTitleOffset(1.1);
      hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
    }
  }

  legend2->Draw();

  canvas2->SaveAs("particlecomposition_result.gif");
  canvas2->SaveAs("particlecomposition_result.eps");
}

void mergeCorrectionsWithDifferentCrosssections(Int_t origin, Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName=0) {
  //
  // Function used to merge standard corrections with vertex
  // reconstruction corrections obtained by a certain mix of ND, DD
  // and SD events.
  //
  // the dn/deta spectrum is corrected and the ratios
  // (standard to changed x-section) of the different dN/deta
  // distributions are saved to a file.
  //
  // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
  //    kINEL = 3
  //    kNSD = 4
  //    kOnePart = 6

  if (outputFileName == 0)
  {
    if (correctionTarget == 3)
      outputFileName = "systematics_vtxtrigger_compositions_inel.root";
    if (correctionTarget == 4)
      outputFileName = "systematics_vtxtrigger_compositions_nsd.root";
    if (correctionTarget == 6)
      outputFileName = "systematics_vtxtrigger_compositions_onepart.root";
  }

  loadlibs();

  const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };

  //Karel:
//     fsd = 0.153 +- 0.031 (0.050 to take into account SD definition) --> change
//     fdd = 0.080 +- 0.050 --> change 
//     fnd = 0.767 +- 0.059 --> keep (error small)

//  const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdlessddmore25", "sdmoreddless25"};
  //Float_t scalesND[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0 };
  //Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
  //Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 0.5, 1.5, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
  //Float_t scalesDD[] = {1.0, 1.4, 0.6, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
  //Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 0.4, 1.6, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
/*  Int_t nChanges = 9;

  const Char_t* changes[]  = { "ua5","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
  Float_t scalesND[] = {0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767};
  Float_t scalesDD[] = {0.080, 0.130, 0.030, 0.080, 0.080, 0.130, 0.030, 0.130, 0.030};
  Float_t scalesSD[] = {0.153, 0.153, 0.153, 0.203, 0.103, 0.203, 0.103, 0.103, 0.203};*/
  
  Float_t ref_SD = -1;
  Float_t ref_DD = -1;
  Float_t ref_ND = -1;
  
  Float_t error_SD = -1;
  Float_t error_DD = -1;
  Float_t error_ND = -1;
  
  GetRelativeFractions(origin, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
  
  Printf("Using x-sections:\n SD: %f +- %f\n DD: %f +- %f\n ND: %f +- %f", ref_SD, error_SD, ref_DD, error_DD, ref_ND, error_ND);
  
  const Char_t* changes[]  = { "default","sdless","sdmore","ddless","ddmore", "dless", "dmore", "sdlessddmore", "sdmoreddless" };
  Int_t nChanges = 9;
  Float_t scalesSD[9];
  Float_t scalesDD[9];
  Float_t scalesND[9];
  
  if (1)
  {
    // sample 8 points on the error ellipse
    for (Int_t i=0; i<9; i++)
    {
      Float_t factorSD = 0;
      Float_t factorDD = 0;
      
      if (i > 0 && i < 3)
        factorSD = (i % 2 == 0) ? 1 : -1;
      else if (i >= 3 && i < 5)
        factorDD = (i % 2 == 0) ? 1 : -1;
      else if (i >= 5 && i < 9)
      {
        factorSD = ((i % 2 == 0) ? 1.0 : -1.0) / TMath::Sqrt(2);
        if (i == 5 || i == 6)
          factorDD = factorSD;
        else
          factorDD = -factorSD;
      }
      
      scalesSD[i] = ref_SD + factorSD * error_SD;
      scalesDD[i] = ref_DD + factorDD * error_DD;
      scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
      
      Printf("Case %d: SD: %f DD: %f ND: %f", i, scalesSD[i], scalesDD[i], scalesND[i]);
    }
  }
  else
  {
    Printf("WARNING: Special treatment for ratios active");
    gSystem->Sleep(1000);
    
    // constrained values by allowed changing of cross-sections
    Float_t pythiaScaling = 0.224 / 0.189;

    if (origin == 10)
    {
      // 900 GeV
      for (Int_t i=0; i<9; i++)
      {
        scalesSD[i] = 15.3;
        scalesDD[i] = 9.5;
      }

      scalesSD[1] = 15.7;
      scalesSD[2] = 17.6;
      scalesSD[3] = 13.5;
      scalesSD[4] = 17.6;

      scalesDD[5] = 15.5;
      scalesDD[6] = 8.8;
      scalesDD[7] = 13.8;
      scalesDD[8] = 7.6;
    }
    else if (origin == 20)
    {
      // 2.36 TeV
      pythiaScaling = 0.217 / 0.167;
      
      for (Int_t i=0; i<9; i++)
      {
        scalesSD[i] = 15.9;
        scalesDD[i] = 10.7;
      }

      scalesSD[1] = 13.5;
      scalesSD[2] = 15.2;
      scalesSD[3] = 13.5;
      scalesSD[4] = 17.6;

      scalesDD[5] = 13.8;
      scalesDD[6] = 7.6;
      scalesDD[7] = 13.8;
      scalesDD[8] = 7.6;
    }
    else
      AliFatal("Not supported");

    for (Int_t i=0; i<9; i++)
    {
      scalesSD[i] /= 100;
      scalesSD[i] *= pythiaScaling;
      scalesDD[i] /= 100;
      scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
      Printf("Case %d: SD: %f DD: %f ND: %f", i, scalesSD[i], scalesDD[i], scalesND[i]);
    }
  }
  
  Int_t backgroundEvents = 0;
  
  //backgroundEvents = 1162+434; // Michele for MB1, run 104892, 15.02.10
  //backgroundEvents = 6;          // Michele for V0AND, run 104892, 15.02.10
  
  //backgroundEvents = 4398+961;   // Michele for MB1, run 104824-52, 16.02.10
  //backgroundEvents = 19;         // Michele for V0AND, run 104824-52, 16.02.10
  
  backgroundEvents = -1;    // use 0 bin from MC! for 2.36 TeV
  
  Printf("Subtracting %d background events!!!", backgroundEvents);
  gSystem->Sleep(1000);
  
  /*
  const Char_t* changes[]  = { "pythia", "qgsm", "phojet"};
  Float_t scalesND[] = {1.0, 1.10, 1.11};
  Float_t scalesSD[] = {1.0, 0.69, 0.86};
  Float_t scalesDD[] = {1.0, 0.98, 0.61};
  Int_t nChanges = 3;
  */
  
  // cross section from Pythia
  // 14 TeV!
//   Float_t sigmaND = 55.2;
//   Float_t sigmaDD = 9.78;
//   Float_t sigmaSD = 14.30;

  // standard correction
  TFile::Open(correctionFileName);
  AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
  correctionStandard->LoadHistograms();

  // dont take vertexreco from this one
  correctionStandard->GetVertexRecoCorrection()->Reset();
  // dont take triggerbias from this one
  correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
  correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
  correctionStandard->GetTriggerBiasCorrectionND()->Reset();
  correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();

  AlidNdEtaCorrection* corrections[100];
  TH1F* hRatios[100];

  Int_t counter = 0;
  for (Int_t j=2; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)

    for (Int_t i=0; i<nChanges; i++) {
      TFile::Open(correctionFileName);

      TString name;
      name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
      AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
      current->LoadHistograms("dndeta_correction");
      current->Reset();

      name.Form("dndeta_correction_ND");
      AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
      dNdEtaCorrectionND->LoadHistograms();
      name.Form("dndeta_correction_DD");
      AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
      dNdEtaCorrectionDD->LoadHistograms();
      name.Form("dndeta_correction_SD");
      AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
      dNdEtaCorrectionSD->LoadHistograms();

      // calculating relative
      Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
      Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
      Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
      Float_t total = nd + dd + sd;
      
      nd /= total;
      sd /= total;
      dd /= total;
      
      Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
      
      Float_t scaleND = scalesND[i] / nd;
      Float_t scaleDD = scalesDD[i] / dd;
      Float_t scaleSD = scalesSD[i] / sd;
      
      Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);      
      
/*      Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
      Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
      Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));

      printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));*/
      current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",scaleND,scaleDD,scaleSD));
      current->SetTitle(name);

      // scale
      if (j == 0 || j == 2)
      {
        dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
        dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
        dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
      }
      if (j == 1 || j == 2)
      {
        dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
        dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
        dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);

        dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
        dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
        dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);

        dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
        dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
        dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);

        dNdEtaCorrectionND->GetTriggerBiasCorrectionOnePart()->Scale(scaleND);
        dNdEtaCorrectionDD->GetTriggerBiasCorrectionOnePart()->Scale(scaleDD);
        dNdEtaCorrectionSD->GetTriggerBiasCorrectionOnePart()->Scale(scaleSD);
      }

      //clear track in correction
      dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
      dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
      dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();

      TList collection;
      collection.Add(correctionStandard);
      collection.Add(dNdEtaCorrectionND);
      collection.Add(dNdEtaCorrectionDD);
      collection.Add(dNdEtaCorrectionSD);

      current->Merge(&collection);
      current->Finish();
      
      // print 0 bin efficiency
      TH1* hist2 = current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
      if (hist2->GetBinContent(1))
      {
        Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
        Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
      }

      corrections[counter] = current;

      // now correct dNdeta distribution with modified correction map
      TFile::Open(analysisFileName);

      dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
      fdNdEtaAnalysis->LoadHistograms();

      fdNdEtaAnalysis->Finish(current, 0.151, correctionTarget, Form("%d %d", j, i), backgroundEvents);

      name = "ratio";
      if (j==0) name.Append("_vetexReco_");
      if (j==1) name.Append("_triggerBias_");
      if (j==2) name.Append("_vertexReco_triggerBias_");
      name.Append(changes[i]);

      hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);

      name.Form("ND #times %0.2f DD #times %0.2f, SD #times %0.2f",scaleND,scaleDD,scaleSD);
      hRatios[counter]->SetTitle(name.Data());
      hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default cross-section}{modified cross-sections}");
      
      TF1* pol0 = new TF1("pol0", "[0]", -0.49, 0.49);
      pol0->SetParameter(0, 0);
      hRatios[counter]->Fit(pol0, "RN");
      Printf("Case %d: %f", i, pol0->GetParameter(0));
      
      if (counter > 0)
        hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);

      delete fdNdEtaAnalysis;

      counter++;
    }
  }

  TFile* fout = new TFile(outputFileName,"RECREATE");

  // to make everything consistent
  hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);

  for (Int_t i=0; i<counter; i++)
  {
    corrections[i]->SaveHistograms();
    hRatios[i]->Write();
  }

  fout->Write();
  fout->Close();
}

void GetRelativeFractions(Int_t origin, Float_t& ref_SD, Float_t& ref_DD, Float_t& ref_ND, Float_t& error_SD, Float_t& error_DD, Float_t& error_ND)
{
  // origin: 
  //   -1 = Pythia (test)
  //   0 = UA5
  //   1 = Data 1.8 TeV
  //   2 = Tel-Aviv
  //   3 = Durham
  //

  switch (origin)
  {
    case -10: // Pythia default at 7 GeV, 50% error
      Printf("PYTHIA x-sections");
      ref_SD = 0.192637; error_SD = ref_SD * 0.5;
      ref_DD = 0.129877; error_DD = ref_DD * 0.5;
      ref_ND = 0.677486; error_ND = 0;
      break;

    case -1: // Pythia default at 900 GeV, as test
      Printf("PYTHIA x-sections");
      ref_SD = 0.223788;
      ref_DD = 0.123315;
      ref_ND = 0.652897;
      break;
      
    case 0: // UA5
      Printf("UA5 x-sections a la first paper");
      ref_SD = 0.153; error_SD = 0.05;
      ref_DD = 0.080; error_DD = 0.05;
      ref_ND = 0.767; error_ND = 0;
      break;
      
    case 10: // UA5
      Printf("UA5 x-sections hadron level definition for Pythia"); 
      // Fractions in Pythia with UA5 cuts selection for SD
      // ND: 0.688662
      // SD: 0.188588 --> this should be 15.3
      // DD: 0.122750
      ref_SD = 0.224 * 0.153 / 0.189; error_SD = 0.023 * 0.224 / 0.189;
      ref_DD = 0.095;                 error_DD = 0.06; 
      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
      break;
    
    case 11: // UA5
      Printf("UA5 x-sections hadron level definition for Phojet"); 
      // Fractions in Phojet with UA5 cuts selection for SD
      // ND: 0.783573
      // SD: 0.151601 --> this should be 15.3
      // DD: 0.064827
      ref_SD = 0.191 * 0.153 / 0.152; error_SD = 0.023 * 0.191 / 0.152;
      ref_DD = 0.095;                 error_DD = 0.06; 
      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
      break;
      
    case 20: // E710, 1.8 TeV
      Printf("E710 x-sections hadron level definition for Pythia");
      // ND: 0.705709
      // SD: 0.166590 --> this should be 15.9
      // DD: 0.127701
      ref_SD = 0.217 * 0.159 / 0.167; error_SD = 0.024 * 0.217 / 0.167;
      ref_DD = 0.075 * 1.43;          error_DD = 0.02 * 1.43; 
      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
      break;
    
    case 21: // E710, 1.8 TeV
      Printf("E710 x-sections hadron level definition for Phojet"); 
      // ND: 0.817462
      // SD: 0.125506 --> this should be 15.9
      // DD: 0.057032
      ref_SD = 0.161 * 0.159 / 0.126; error_SD = 0.024 * 0.161 / 0.126;
      ref_DD = 0.075 * 1.43;         error_DD = 0.02 * 1.43;
      ref_ND = 1.0 - ref_SD - ref_DD; error_ND = 0;
      break;
    
    case 1: // data 1.8 TeV
      Printf("??? x-sections");
      ref_SD = 0.152;
      ref_DD = 0.092;
      ref_ND = 1 - ref_SD - ref_DD;
      break;
      
    case 2: // tel-aviv model
      Printf("Tel-aviv model x-sections");
      ref_SD = 0.171;
      ref_DD = 0.094;
      ref_ND = 1 - ref_SD - ref_DD;
      break;
    
    case 3: // durham model
      Printf("Durham model x-sections");
      ref_SD = 0.190;
      ref_DD = 0.125;
      ref_ND = 1 - ref_SD - ref_DD;
      break;
    
    default:
      AliFatal(Form("Unknown origin %d", origin));
  }
}

void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctionFileName="correction_mapprocess-types.root", const Char_t* outputFileName="correction_map2.root") {
  //
  // Function used to merge standard corrections with vertex
  // reconstruction corrections obtained by a certain mix of ND, DD
  // and SD events.
  //
  loadlibs();

  const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
  
  Float_t ref_SD = -1;
  Float_t ref_DD = -1;
  Float_t ref_ND = -1;
  
  Float_t error_SD = -1;
  Float_t error_DD = -1;
  Float_t error_ND = -1;
  
  GetRelativeFractions(origin, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
  
  Printf("Using x-sections:\n SD: %f +- %f\n DD: %f +- %f\n ND: %f +- %f", ref_SD, error_SD, ref_DD, error_DD, ref_ND, error_ND);
  
//Karel (UA5):
//     fsd = 0.153 +- 0.031
//     fdd = 0.080 +- 0.050
//     fnd = 0.767 +- 0.059

//       Karel (1.8 TeV):
//       
//       Tel-Aviv model Sd/Inel = 0.171           Dd/Inel = 0.094
//       Durham model   Sd/Inel = 0.190           Dd/Inel = 0.125
//       Data           Sd/Inel = 0.152 +- 0.030  Dd/Inel = 0.092 +- 0.45

  // standard correction
  TFile::Open(correctionFileName);
  AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
  correctionStandard->LoadHistograms();

  // dont take vertexreco from this one
  correctionStandard->GetVertexRecoCorrection()->Reset();
  // dont take triggerbias from this one
  correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
  correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
  correctionStandard->GetTriggerBiasCorrectionND()->Reset();
  correctionStandard->GetTriggerBiasCorrectionOnePart()->Reset();

  AlidNdEtaCorrection* corrections[100];
  TH1F* hRatios[100];

  Int_t counter = 0;
      
  TFile::Open(correctionFileName);

  AlidNdEtaCorrection* current = new AlidNdEtaCorrection("dndeta_correction_ua5", "dndeta_correction_ua5");
  current->LoadHistograms("dndeta_correction");
  current->Reset();

  TString name;
  name.Form("dndeta_correction_ND");
  AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
  dNdEtaCorrectionND->LoadHistograms();
  name.Form("dndeta_correction_DD");
  AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
  dNdEtaCorrectionDD->LoadHistograms();
  name.Form("dndeta_correction_SD");
  AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
  dNdEtaCorrectionSD->LoadHistograms();

  // calculating relative
  Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
  Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
  Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
  Float_t total = nd + dd + sd;
  
  nd /= total;
  sd /= total;
  dd /= total;
  
  Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
  
  Float_t scaleND = ref_ND / nd;
  Float_t scaleDD = ref_DD / dd;
  Float_t scaleSD = ref_SD / sd;
  
  Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);

  // scale
  dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
  dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
  dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
    
  dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
  dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
  dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);

  dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
  dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
  dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);

  dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
  dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
  dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);

  dNdEtaCorrectionND->GetTriggerBiasCorrectionOnePart()->Scale(scaleND);
  dNdEtaCorrectionDD->GetTriggerBiasCorrectionOnePart()->Scale(scaleDD);
  dNdEtaCorrectionSD->GetTriggerBiasCorrectionOnePart()->Scale(scaleSD);

  //clear track in correction
  dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
  dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
  dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();

  TList collection;
  collection.Add(correctionStandard);
  collection.Add(dNdEtaCorrectionND);
  collection.Add(dNdEtaCorrectionDD);
  collection.Add(dNdEtaCorrectionSD);

  current->Merge(&collection);
  current->Finish();

  // print 0 bin efficiency
  TH1* hist2 = current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->Get1DCorrection("y", -10, 10);
  if (hist2->GetBinContent(1) > 0)
  {
    Float_t triggerEff = 100.0 / hist2->GetBinContent(1);
    Printf("trigger eff in 0 bin is: %.2f %%", triggerEff);
  }
  
  TFile* fout = new TFile(outputFileName,"RECREATE");
  current->SaveHistograms();

  fout->Write();
  fout->Close();

  Printf("Trigger efficiencies:");
  Printf("ND: %.2f %%", 100.0 * dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
  Printf("SD: %.2f %%", 100.0 * dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
  Printf("DD: %.2f %%", 100.0 * dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
  Printf("INEL: %.2f %%", 100.0 * current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
  Printf("NSD: %.2f %%", 100.0 * (dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() + dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral()) / (dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral() + dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral()));
}

DrawTriggerEfficiency(Char_t* fileName) {

  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetOptFit(0);

  gStyle->SetTextSize(0.04);
  gStyle->SetTitleSize(0.05,"xyz");
  //gStyle->SetTitleFont(133, "xyz");
  //gStyle->SetLabelFont(133, "xyz");
  //gStyle->SetLabelSize(17, "xyz");
  gStyle->SetLabelOffset(0.01, "xyz");

  gStyle->SetTitleOffset(1.1, "y");
  gStyle->SetTitleOffset(1.1, "x");
  gStyle->SetEndErrorSize(0.0);

  //##############################################

  //making canvas and pads
  TCanvas *c = new TCanvas("trigger_eff", "",600,500);

  TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);

  p1->SetBottomMargin(0.15);
  p1->SetTopMargin(0.03);
  p1->SetLeftMargin(0.15);
  p1->SetRightMargin(0.03);
  
  p1->SetGridx();
  p1->SetGridy();

  p1->Draw();
  p1->cd();

  gSystem->Load("libPWG0base");

  AlidNdEtaCorrection* corrections[4];
  AliCorrectionMatrix2D* triggerBiasCorrection[4];

  TH1F* hTriggerEffInv[4];
  TH1F* hTriggerEff[4];

  Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
  
  for (Int_t i=0; i<4; i++) {
    corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
    corrections[i]->LoadHistograms(fileName, names[i]);    

    triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();

    
    hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
    hTriggerEff[i]    = (TH1F*)hTriggerEffInv[i]->Clone();
    
    for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
      hTriggerEff[i]->SetBinContent(b,1);
      hTriggerEff[i]->SetBinError(b,0);
    }
    hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
    hTriggerEff[i]->Scale(100);
  }

  Int_t colors[] = {2,3,4,1};
  Float_t effs[4];
  for (Int_t i=0; i<4; i++) {
    hTriggerEff[i]->Fit("pol0","","",-20,20);
    effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
    ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
    ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
    ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
    cout << effs[i] << endl;
  }


  Char_t* text[] = {"ND", "DD", "SD", "INEL"};
  TLatex* latex[4];

  TH2F* null = new TH2F("","",100,-25,35,100,0,110);
  null->SetXTitle("Vertex z [cm]");
  null->GetXaxis()->CenterTitle(kTRUE);
  null->SetYTitle("Trigger efficiency [%]");
  null->GetYaxis()->CenterTitle(kTRUE);
  null->Draw();


  for (Int_t i=0; i<4; i++) {
    hTriggerEff[i]->SetLineWidth(2);
    hTriggerEff[i]->SetLineColor(colors[i]);

    hTriggerEff[i]->Draw("same");

    latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
    latex[i]->SetTextColor(colors[i]);
    latex[i]->Draw();
  }
  
}


DrawSpectraPID(Char_t* fileName) {

  gSystem->Load("libPWG0base");

  Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
  AlidNdEtaCorrection* corrections[4];
  AliCorrectionMatrix3D* trackToPartCorrection[4];

  TH1F* measuredPt[4];
  TH1F* generatedPt[4];
  TH1F* ratioPt[4];

  for (Int_t i=0; i<4; i++) {
    corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
    corrections[i]->LoadHistograms(fileName, names[i]);    

    trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();

    Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
    Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
    Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
    Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);

    measuredPt[i]  = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
    generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
    ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
    ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
  }
  
  ratioPt[0]->Draw();
  
  for (Int_t i=0; i<3; i++) {
    ratioPt[i]->SetLineColor(i+1);
    ratioPt[i]->SetLineWidth(2);
    
    ratioPt[i]->Draw("same");
    
  }

  return;
  measuredPt[0]->SetLineColor(2);
  measuredPt[0]->SetLineWidth(5);

  measuredPt[0]->Draw();
  generatedPt[0]->Draw("same");
}

void changePtSpectrum(const char* fileName = "analysis_mc.root", Float_t ptCutOff = 0.2, const char* fileName2 = 0)
{
  Float_t factor = 0.5;

  TFile* file = TFile::Open(fileName);
  TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt")->Clone());
  
  TH1* hist2 = 0;
  if (fileName2)
  {
    file2 = TFile::Open(fileName2);
    hist2 = dynamic_cast<TH1*> (file2->Get("dndeta_check_pt")->Clone());
    hist2->Scale(hist->GetBinContent(hist->FindBin(ptCutOff)) / hist2->GetBinContent(hist2->FindBin(ptCutOff)));
  }
  
  //hist->Scale(1.0 / hist->Integral());

  //hist->Rebin(3);
  //hist->Scale(1.0/3);

  TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
  TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));

  TH1F* scale1 =  dynamic_cast<TH1F*> (hist->Clone("scale1"));
  TH1F* scale2 =  dynamic_cast<TH1F*> (hist->Clone("scale2"));

  for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
  {
    if (hist->GetBinCenter(i) > ptCutOff)
    {
      scale1->SetBinContent(i, 1);
      scale2->SetBinContent(i, 1);
    }
    else
    {
      // 90 % at pt = 0, 0% at pt = ptcutoff
      scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);

      // 110% at pt = 0, ...
      scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
    }
    scale1->SetBinError(i, 0);
    scale2->SetBinError(i, 0);
  }

  /*
  new TCanvas;
  scale1->Draw();
  scale2->SetLineColor(kRed);
  scale2->Draw("SAME");
  */

  clone1->Multiply(scale1);
  clone2->Multiply(scale2);

  Prepare1DPlot(hist);
  Prepare1DPlot(clone1);
  Prepare1DPlot(clone2);

  /*hist->SetMarkerStyle(markers[0]);
  clone1->SetMarkerStyle(markers[0]);
  clone2->SetMarkerStyle(markers[0]);*/

  hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
  hist->GetXaxis()->SetRangeUser(0, 0.5);
  hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
  hist->GetYaxis()->SetTitleOffset(1);

  TCanvas* canvas = new TCanvas("c", "c", 600, 600);
  gPad->SetGridx();
  gPad->SetGridy();
  gPad->SetTopMargin(0.05);
  gPad->SetRightMargin(0.05);
  gPad->SetBottomMargin(0.12);
  hist->Draw("H");
  clone1->SetLineColor(kRed);
  clone1->Draw("HSAME");
  clone2->SetLineColor(kBlue);
  clone2->Draw("HSAME");
  hist->Draw("HSAME");
  
  if (hist2)
  {
    Prepare1DPlot(hist2);
    hist2->SetLineStyle(2);
    hist2->Draw("HSAME");
  }

  Float_t fraction =  hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
  Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
  Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());

  printf("%f %f %f\n", fraction, fraction1, fraction2);
  printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);

  //canvas->SaveAs("changePtSpectrum.gif");
  canvas->SaveAs("changePtSpectrum.eps");
}

void FractionBelowPt()
{
  gSystem->Load("libPWG0base");

  AlidNdEtaCorrection* fdNdEtaCorrection[4];

  TFile::Open("systematics.root");

  for (Int_t i=0; i<4; ++i)
  {
    TString name;
    name.Form("correction_%d", i);
    fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
    fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
  }

  Double_t geneCount[5];
  Double_t measCount[5];
  geneCount[4] = 0;
  measCount[4] = 0;

  for (Int_t i=0; i<4; ++i)
  {
    TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
    geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
                                  hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
                                  1, hist->GetZaxis()->FindBin(0.3));

    hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
    measCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10), hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8), 1, hist->GetZaxis()->FindBin(0.3));

    geneCount[4] += geneCount[i];
    measCount[4] += measCount[i];

    printf("Particle %s: %d gene, %d meas\n", ((i == 0) ? "pi" : (i == 1) ? "K" : (i == 2) ? "p" : "others"), (Int_t) geneCount[i], (Int_t) measCount[i]);
  }

  printf("Generated ratios are:     %f pi, %f K, %f p, %f others\n", geneCount[0] / geneCount[4], geneCount[1] / geneCount[4], geneCount[2] / geneCount[4], geneCount[3] / geneCount[4]);

  printf("Reconstructed ratios are: %f pi, %f K, %f p, %f others\n", measCount[0] / measCount[4], measCount[1] / measCount[4], measCount[2] / measCount[4], measCount[3] / measCount[4]);
}


mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
					     Char_t* misalignedFile = "correction_map_misaligned.root",
					     Char_t* outputFileName="correction_map_misaligned_single.root")
{
  //
  // from the aligned and misaligned corrections, 3 new corrections are created
  // in these new corrections only one of the corrections (track2particle, vertex, trigger)
  // is taken from the misaligned input to allow study of the effect on the different
  // corrections

  gSystem->Load("libPWG0base");

  const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };

  AlidNdEtaCorrection* corrections[3];
  for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
    AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
    alignedCorrection->LoadHistograms(alignedFile);

    AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
    misalignedCorrection->LoadHistograms(misalignedFile);

    TString name;
    name.Form("dndeta_correction_alignment_%s", typeName[j]);
    AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);

    switch (j)
    {
      case 0:
        alignedCorrection->GetTrack2ParticleCorrection()->Reset();
        misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
        misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
        misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
        misalignedCorrection->GetVertexRecoCorrection()->Reset();
        break;

      case 1:
        misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
        misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
        misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
        misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
        alignedCorrection->GetVertexRecoCorrection()->Reset();
        break;

      case 2:
        misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
        alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
        alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
        alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
        misalignedCorrection->GetVertexRecoCorrection()->Reset();
        break;

      default:
        return;
    }

    TList collection;
    collection.Add(misalignedCorrection);
    collection.Add(alignedCorrection);

    current->Merge(&collection);
    current->Finish();

    corrections[j] = current;
  }

  TFile* fout = new TFile(outputFileName, "RECREATE");

  for (Int_t i=0; i<3; i++)
    corrections[i]->SaveHistograms();

  fout->Write();
  fout->Close();
}


void DrawdNdEtaDifferencesAlignment()
{
  TH1* hists[5];

  TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
  legend->SetFillColor(0);

  TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
  canvas->Divide(2, 1);

  canvas->cd(1);

  for (Int_t i=0; i<5; ++i)
  {
    hists[i] = 0;
    TFile* file = 0;
    TString title;

    switch(i)
    {
      case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
      case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
      case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
      case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
      case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
      default: return;
    }

    if (file)
    {
      hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
      hists[i]->SetTitle("");

      Prepare1DPlot(hists[i], kFALSE);
      hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
      hists[i]->GetYaxis()->SetRangeUser(6, 7);
      hists[i]->SetLineWidth(1);
      hists[i]->SetLineColor(colors[i]);
      hists[i]->SetMarkerColor(colors[i]);
      hists[i]->SetMarkerStyle(markers[i]);
      hists[i]->GetXaxis()->SetLabelOffset(0.015);
      hists[i]->GetYaxis()->SetTitleOffset(1.5);
      gPad->SetLeftMargin(0.12);
      hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));

      legend->AddEntry(hists[i], title);
      hists[i]->SetTitle(title);
    }
  }
  legend->Draw();

  canvas->cd(2);
  gPad->SetLeftMargin(0.14);

  TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
  legend2->SetFillColor(0);

  for (Int_t i=1; i<5; ++i)
  {
    if (hists[i])
    {
      legend2->AddEntry(hists[i]);

      hists[i]->Divide(hists[0]);
      hists[i]->SetTitle("b)");
      hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
      hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
      hists[i]->GetYaxis()->SetTitleOffset(1.8);
      hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
    }
  }

  legend2->Draw();

  canvas->SaveAs("misalignment_result.eps");
  canvas->SaveAs("misalignment_result.gif");
}

void EvaluateMultiplicityEffect()
{
  gSystem->Load("libPWG0base");

  AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
  AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];

  TFile::Open("systematics-low-multiplicity.root");

  for (Int_t i=0; i<4; ++i)
  {
    TString name;
    name.Form("correction_%d", i);
    fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
    fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
  }

  TList list;
  for (Int_t i=1; i<4; ++i)
    list.Add(fdNdEtaCorrectionLow[i]);

  fdNdEtaCorrectionLow[0]->Merge(&list);
  fdNdEtaCorrectionLow[0]->Finish();

  TFile::Open("systematics-high-multiplicity.root");

  for (Int_t i=0; i<4; ++i)
  {
    TString name;
    name.Form("correction_%d", i);
    fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
    fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
  }

  TList list2;
  for (Int_t i=1; i<4; ++i)
    list2.Add(fdNdEtaCorrectionHigh[i]);

  fdNdEtaCorrectionHigh[0]->Merge(&list2);
  fdNdEtaCorrectionHigh[0]->Finish();

  TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
  TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);

  TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
  TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
  for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
    for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
      for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
      //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
      {
        if (hist->GetBinContent(x, y, z) > 0)
          outputLow->Fill(hist->GetBinContent(x, y, z));
        //if (hist->GetBinContent(x, y, z) == 1)
        //  printf("z = %f, eta = %f, pt = %f: %f %f %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z), hist->GetBinContent(x, y, z), fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->GetBinContent(x, y, z), fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->GetBinContent(x, y, z));

        if (hist2->GetBinContent(x, y, z) > 0)
          outputHigh->Fill(hist2->GetBinContent(x, y, z));
      }

  TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
  canvas->Divide(2, 1);

  canvas->cd(1);
  outputLow->Draw();
  outputLow->Fit("gaus", "0");
  outputLow->GetFunction("gaus")->SetLineColor(2);
  outputLow->GetFunction("gaus")->DrawCopy("SAME");

  canvas->cd(2);
  outputHigh->Draw();
  outputHigh->Fit("gaus", "0");
  outputHigh->GetFunction("gaus")->DrawCopy("SAME");

  canvas->SaveAs(Form("%s.gif", canvas->GetName()));
}

void PlotErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
{
    gSystem->Load("libPWG0base");
	
	AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
    dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
    
    dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
}


void runStudy(const char* baseCorrectionMapFile = "correction_map.root", const char* baseCorrectionMapFolder = "dndeta_correction", const char* changedCorrectionMapFile = "correction_map.root", const char* changedCorrectionMapFolder = "dndeta_correction", const char* dataFile = "analysis_esd_raw.root", const char* output = "analysis_esd_syst.root")
{
  gSystem->Load("libPWG0base");

  TFile* file = TFile::Open(output, "RECREATE");

  const Int_t max = 5;
  dNdEtaAnalysis* fdNdEtaAnalysis[5];

  new TCanvas;
  TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
  legend->SetFillColor(0);

  for (Int_t i = 0; i < max; ++i)
  {
    TFile::Open(baseCorrectionMapFile);
    AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
    baseCorrection->LoadHistograms();

    AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
    const char* name = 0;

    TFile::Open(changedCorrectionMapFile);
    switch (i)
    {
      case 0 :
        name = "default";
        break;

      case 1 :
        baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
        name = "Track2Particle";
        break;

      case 2 :
        baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
        name = "VertexReco";
        break;

      case 3 :
        baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
        name = "TriggerBias_MBToINEL";
        break;

      case 4 :
        baseCorrection->LoadHistograms(changedCorrectionMapFolder);
        name = "all";
        break;

      default: return;
    }

    TFile::Open(dataFile);
    fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
    fdNdEtaAnalysis[i]->LoadHistograms("dndeta");

    fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
    file->cd();
    fdNdEtaAnalysis[i]->SaveHistograms();

    TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
    hist->SetLineColor(colors[i]);
    hist->SetMarkerColor(colors[i]);
    hist->SetMarkerStyle(markers[i]+1);
    hist->DrawCopy((i == 0) ? "" : "SAME");
    legend->AddEntry(hist, name);
  }

  legend->Draw();
}

void ChangePtInCorrection(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;

  dNdEtaCorrection->Finish();

  AliCorrection* corr = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kTrack2Particle);
 
  Printf(">>>> Before");
  corr->PrintInfo(0);

  Float_t factor = 0.5;
  Float_t ptCutOff = 0.2;
  
  TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
  TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
  
  for (Int_t z = 1; z <= gene->GetZaxis()->FindBin(ptCutOff - 0.01); z++)
  {
    Float_t localFactor = 1 - (ptCutOff - gene->GetZaxis()->GetBinCenter(z)) / ptCutOff * factor;
    Printf("%f %f", gene->GetZaxis()->GetBinCenter(z), localFactor);
    for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
      for (Int_t y = 1; y <= gene->GetNbinsY(); y++)
      {
        gene->SetBinContent(x, y, z, gene->GetBinContent(x, y, z) * localFactor);
        meas->SetBinContent(x, y, z, meas->GetBinContent(x, y, z) * localFactor);
      }
  }
  
  dNdEtaCorrection->Finish();
  
  Printf(">>>> After");
  corr->PrintInfo(0);
}

Float_t FitAverage(TH1* hist, Int_t n, Float_t* begin, Float_t *end, Int_t color, Int_t& totalBins)
{
  Float_t average = 0;
  totalBins = 0;
  
  for (Int_t i=0; i<n; i++)
  {
    func = new TF1("func", "[0]", hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->FindBin(begin[i])), hist->GetXaxis()->GetBinUpEdge(hist->GetXaxis()->FindBin(end[i])));
    Int_t bins = hist->GetXaxis()->FindBin(end[i]) - hist->GetXaxis()->FindBin(begin[i]) + 1;
    func->SetParameter(0, 1);
    func->SetLineColor(color);

    hist->Fit(func, "RNQ");
    func->Draw("SAME");
    
    average += func->GetParameter(0) * bins;
    totalBins += bins;
  }
  
  return average / totalBins;
}

void SPDIntegrateGaps(Bool_t all, const char* mcFile = "../../../LHC10b2/v0or/spd/analysis_esd_raw.root")
{
  Float_t eta = 1.29;
  Int_t binBegin = ((TH2*) gFile->Get("fEtaPhi"))->GetXaxis()->FindBin(-eta);
  Int_t binEnd   = ((TH2*) gFile->Get("fEtaPhi"))->GetXaxis()->FindBin(eta);
  
  Printf("eta range: %f bins: %d %d", eta, binBegin, binEnd);
  
  if (!all)
    Printf("Eta smaller than 0 side");
  
  c = new TCanvas;
  TFile::Open("analysis_esd_raw.root");
  hist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("hist", binBegin, (all) ? binEnd : 40);
  hist->Rebin(2);
  hist->SetStats(0);
  hist->Sumw2();
  hist->Draw("HIST");
  gPad->SetGridx();
  gPad->SetGridy();
  
  TFile::Open(mcFile);  
  mcHist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("mcHist", binBegin, (all) ? binEnd : 40);
  mcHist->Rebin(2);
  mcHist->SetLineColor(2);
  mcHist->Scale(hist->Integral() / mcHist->Integral());
  mcHist->Draw("SAME");
  
  Float_t add = 0;
  Int_t bins;
  
  Float_t okRangeBegin[] = { 0.04, 0.67, 1.34 };
  Float_t okRangeEnd[] =   { 0.55, 1.24, 1.63 };
  Float_t gapRangeBegin[] = { 0.6, 1.27  };
  Float_t gapRangeEnd[] =   { 0.65, 1.32 };
  Float_t averageOK  = FitAverage(hist, 3, okRangeBegin, okRangeEnd, 1, bins);
  Float_t averageGap = FitAverage(hist, 2, gapRangeBegin, gapRangeEnd, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin2[] = { 2.4, 2.65 };
  Float_t okRangeEnd2[] =   { 2.55, 3.2 };
  Float_t gapRangeBegin2[] = { 2.59, 3.3 };
  Float_t gapRangeEnd2[] =   { 2.61, 3.3 };
  averageOK  = FitAverage(hist, 2, okRangeBegin2, okRangeEnd2, 1, bins);
  averageGap = FitAverage(hist, 2, gapRangeBegin2, gapRangeEnd2, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin3[] = { 3.55, 3.9 };
  Float_t okRangeEnd3[] =   { 3.8, 4.15 };
  Float_t gapRangeBegin3[] = { 3.83  };
  Float_t gapRangeEnd3[] =   { 3.86 };
  averageOK  = FitAverage(hist, 2, okRangeBegin3, okRangeEnd3, 1, bins);
  averageGap = FitAverage(hist, 1, gapRangeBegin3, gapRangeEnd3, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin4[] = { 4.2, 4.5 };
  Float_t okRangeEnd4[] =   { 4.4, 4.7 };
  Float_t gapRangeBegin4[] = { 4.45  };
  Float_t gapRangeEnd4[] =   { 4.45 };
  averageOK  = FitAverage(hist, 2, okRangeBegin4, okRangeEnd4, 1, bins);
  averageGap = FitAverage(hist, 1, gapRangeBegin4, gapRangeEnd4, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin5[] = { 5.4, 5.7 };
  Float_t okRangeEnd5[] =   { 5.6, 5.8 };
  Float_t gapRangeBegin5[] = { 5.63  };
  Float_t gapRangeEnd5[] =   { 5.67 };
  averageOK  = FitAverage(hist, 2, okRangeBegin5, okRangeEnd5, 1, bins);
  averageGap = FitAverage(hist, 1, gapRangeBegin5, gapRangeEnd5, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Printf("This adds %.2f %% to the total number of tracklets (%f)", 100.0 * add / hist->Integral(), hist->Integral());
  c->SaveAs("gap1.png");
  
  add1 = add;
  total1 = hist->Integral();

  if (all)
    return;

  Printf("\nEta larger than 0 side");
  
  c = new TCanvas;
  TFile::Open("analysis_esd_raw.root");
  hist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("hist2", 41, binEnd);
  hist->Rebin(2);
  hist->SetStats(0);
  hist->Sumw2();
  hist->Draw("HIST");
  gPad->SetGridx();
  gPad->SetGridy();
  
  TFile::Open(mcFile);  
  mcHist = ((TH2*) gFile->Get("fEtaPhi"))->ProjectionY("mcHist", 41, binEnd);
  mcHist->Rebin(2);
  mcHist->SetLineColor(2);
  mcHist->Scale(hist->Integral() / mcHist->Integral());
  mcHist->Draw("SAME");
  
  add = 0;
  
  Float_t okRangeBegin[] = { 0.04, 0.67, 1.34 };
  Float_t okRangeEnd[] =   { 0.55, 1.24, 1.63 };
  Float_t gapRangeBegin[] = { 0.6, 1.27  };
  Float_t gapRangeEnd[] =   { 0.65, 1.32 };
  Float_t averageOK  = FitAverage(hist, 3, okRangeBegin, okRangeEnd, 1, bins);
  Float_t averageGap = FitAverage(hist, 2, gapRangeBegin, gapRangeEnd, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin2[] = { 2.32, 2.65 };
  Float_t okRangeEnd2[] =   { 2.55, 3.2 };
  Float_t gapRangeBegin2[] = { 2.59, 3.3 };
  Float_t gapRangeEnd2[] =   { 2.61, 3.3 };
  averageOK  = FitAverage(hist, 2, okRangeBegin2, okRangeEnd2, 1, bins);
  averageGap = FitAverage(hist, 2, gapRangeBegin2, gapRangeEnd2, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin3[] = { 3.55, 3.9 };
  Float_t okRangeEnd3[] =   { 3.8, 4.15 };
  Float_t gapRangeBegin3[] = { 3.83  };
  Float_t gapRangeEnd3[] =   { 3.86 };
  averageOK  = FitAverage(hist, 2, okRangeBegin3, okRangeEnd3, 1, bins);
  averageGap = FitAverage(hist, 1, gapRangeBegin3, gapRangeEnd3, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);
  
  Float_t okRangeBegin4[] = { 4.2, 4.5 };
  Float_t okRangeEnd4[] =   { 4.4, 4.7 };
  Float_t gapRangeBegin4[] = { 4.45  };
  Float_t gapRangeEnd4[] =   { 4.45 };
  averageOK  = FitAverage(hist, 2, okRangeBegin4, okRangeEnd4, 1, bins);
  averageGap = FitAverage(hist, 1, gapRangeBegin4, gapRangeEnd4, 2, bins);
  add += bins * (averageOK - averageGap);
  Printf("Average OK: %f %f %d: %f", averageOK, averageGap, bins, add);

  Printf("This adds %.2f %% to the total number of tracklets (%f)", 100.0 * add / hist->Integral(), hist->Integral());
  c->SaveAs("gap2.png");
  
  Printf("In total we add %.2f %%.", 100.0 * (add1 + add) / (total1 + hist->Integral()));
}