]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/dNdEta/drawSystematics.C
added correction for events with vertex but 0 tracks
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematics.C
index 2285ab3a6a7387fb976adefcbe62d18b8ab3db42..97dd60ef43c0922d008da11201783343d8a6b5a4 100644 (file)
@@ -23,6 +23,18 @@ void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", c
 
 #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)
 {
@@ -83,8 +95,8 @@ void InitPad()
   //gPad->SetTopMargin(0.13);
   //gPad->SetBottomMargin(0.1);
 
-  //gPad->SetGridx();
-  //gPad->SetGridy();
+  gPad->SetGridx();
+  gPad->SetGridy();
 }
 
 void InitPadCOLZ()
@@ -257,6 +269,60 @@ TH1** DrawRatios(const char* fileName = "systematics.root")
   return ptDists;
 }
 
+void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
+{
+  gROOT->ProcessLine(".L drawPlots.C");
+  gSystem->Load("libPWG0base");
+
+  const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
+  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
+  const char* legendNames[] = { "#pi", "K", "p", "standard" };
+  Int_t folderCount = 3;
+
+  /*const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
+  const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3", "dndeta_correction" };
+  const char* legendNames[] = { "#pi", "K", "p", "others", "standard" };
+  Int_t folderCount = 5;*/
+
+  TString canvasName;
+  canvasName.Form("Track2Particle1DComposition");
+  TCanvas* canvas = new TCanvas(canvasName, canvasName, 700, 500);
+  canvas->SetGridx();
+  canvas->SetGridy();
+  canvas->SetBottomMargin(0.12);
+  //InitPad();
+
+  TLegend* legend = new TLegend(0.8, 0.7, 0.95, 0.95);
+  legend->SetFillColor(0);
+
+  Int_t mycolors[] = {1, 2, 4};
+
+  for (Int_t i=0; i<folderCount; ++i)
+  {
+    Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
+
+    TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folderNames[i], folderNames[i])));
+
+    Prepare1DPlot(corrZ);
+
+    corrZ->SetTitle("");
+    corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
+    corrZ->GetYaxis()->SetRangeUser(0.51, 6);
+    corrZ->SetMarkerColor(mycolors[i]);
+    corrZ->SetLineColor(mycolors[i]);
+    corrZ->SetMarkerStyle(markers[i+1]);
+    corrZ->GetYaxis()->SetTitle("correction factor");
+
+    corrZ->DrawCopy(((i>0) ? "SAMEP" : "P"));
+
+    legend->AddEntry(corrZ, legendNames[i]);
+  }
+
+  legend->Draw();
+
+  canvas->SaveAs("ptcutoff_species.eps");
+}
+
 void DrawCompareToReal()
 {
   gROOT->ProcessLine(".L drawPlots.C");
@@ -277,6 +343,16 @@ void DrawDifferentSpecies()
   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");
@@ -600,167 +676,735 @@ void Composition()
   Track2Particle1DComposition(fileNames, nCompositions, folderNames);
 }
 
-Double_t ConvSigma1To2D(Double_t sigma)
-{
-  return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
-}
 
-Double_t ConvDistance1DTo2D(Double_t distance)
+void drawSystematics()
 {
-  return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
+  //Secondaries();
+  //DrawDifferentSpecies();
+  //Composition();
+
+  Sigma2VertexSimulation();
+
 }
 
-Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
+void DrawdNdEtaDifferences()
 {
-  Double_t count = 0;
+  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);
 
-  //nSigma = ConvSigma1To2D(nSigma);
+  canvas->cd(1);
 
-  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
-    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+  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)
     {
-      Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
-      Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
+      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();
 
-      Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
+  canvas->cd(2);
+  gPad->SetLeftMargin(0.14);
 
-      d = ConvDistance1DTo2D(d);
+  TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
+  legend2->SetFillColor(0);
 
-      if (d < nSigma)
-        count += tracks->GetBinContent(x, y);
+  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");
 
-  return count;
+  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");
 }
 
-TH2F* Sigma2VertexGaussianTracksHist()
-{
-  TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
+mergeCorrectionsWithDifferentCrosssections(Char_t* correctionFileName="correction_map.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
+  //
+  // 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.
+  //
+
+  loadlibs();
+
+  const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
+  const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25"};
+  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};
+
+  // cross section from Pythia
+  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();
+
+  AlidNdEtaCorrection* corrections[100];
+  TH1F* hRatios[100];
+
+  Int_t counter = 0;
+  for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
+
+    for (Int_t i=0; i<17; 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 = 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\%",nd,dd,sd));
+      current->SetTitle(name);
+
+      // scale
+      if (j == 0 || j == 2)
+      {
+        dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scalesSD[i]);
+      }
+      if (j == 1 || j == 2)
+      {
+        dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
+
+        dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scalesSD[i]);
+
+        dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scalesDD[i]);
+        dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scalesSD[i]);
+      }
+
+      //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();
+
+      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.3, AlidNdEtaCorrection::kINEL);
+      //fdNdEtaAnalysis->Finish(current, 0.3, AlidNdEtaCorrection::kNSD);
+
+      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.Append(Form(" (DD #times %0.2f, SD #times %0.2f)",scalesDD[i],scalesSD[i]));
+      hRatios[counter]->SetTitle(name.Data());
+      hRatios[counter]->SetYTitle("ratio (Pythia x-sections)/(changed x-sections)");
+
+      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);
 
-  TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
-  gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
+  for (Int_t i=0; i<51; i++)
+  {
+    corrections[i]->SaveHistograms();
+    hRatios[i]->Write();
+  }
+
+  fout->Write();
+  fout->Close();
+}
+
+
+DrawVertexRecoSyst() {
+  // Draws the ratio of the dN/dEta obtained with changed SD and DD
+  // cross-sections vertex reco corrections to the dN/dEta obtained
+  // from the standard pythia cross-sections 
+  //
+  // The files with the vertex reco corrections for different
+  // processes (and with the other standard corrections) are expected
+  // to have the names "analysis_esd_X.root", where the Xs are defined
+  // in the array changes below.
+
+  Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
+  Char_t* descr[]  =   {"",
+                       "#sigma_{DD}' = 1.5#times#sigma_{DD}",
+                       "#sigma_{DD}' = 0.5#times#sigma_{DD}",
+                       "#sigma_{SD}' = 1.5#times#sigma_{SD}",
+                       "#sigma_{SD}' = 0.5#times#sigma_{SD}",
+                       "#sigma_{D}' = 1.5#times#sigma_{D}",
+                       "#sigma_{D}' = 0.5#times#sigma_{D}"};
+
+  Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.5, 0.5, 1.5, 0.5};
+  Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.5, 0.5};
+
+  Int_t markers[] = {20,20,21,22,23,28,29};
+  Int_t colors[]  = {1,2,3,4,6,8,102};
+
+  // cross section from Pythia
+  Float_t sigmaND = 55.2;
+  Float_t sigmaDD = 9.78;
+  Float_t sigmaSD = 14.30;
+
+  TH1F* dNdEta[7];
+  TH1F* ratios[7];
+
+  TFile* fin;
+
+  for (Int_t i=0; i<7; i++) {
+    // calculating relative
+    fin = TFile::Open(Form("analysis_esd_%s.root",changes[i]));
+
+    dNdEta[i] = (TH1F*)(fin->Get("dndeta/dndeta_dNdEta_corrected_2"))->Clone();
+
+    for (Int_t b=0; b<dNdEta[i]->GetNbinsX(); b++) {
+      if (TMath::Abs(dNdEta[i]->GetBinCenter(b))>0.9) {
+       dNdEta[i]->SetBinContent(b,0);
+       dNdEta[i]->SetBinError(b,0);
+      }
+    }
+
+    dNdEta[i]->Rebin(4);
+
+    dNdEta[i]->SetLineWidth(2);
+    dNdEta[i]->SetLineColor(colors[i]);
+    dNdEta[i]->SetMarkerStyle(markers[i]);
+    dNdEta[i]->SetMarkerSize(0.9);
+    dNdEta[i]->SetMarkerColor(colors[i]);
+
+    ratios[i] = (TH1F*)dNdEta[i]->Clone("ratio");
+    ratios[i]->Divide(ratios[i],dNdEta[0],1,1,"B");
+    
+    ratios[i]->SetName(changes[i]);
+    ratios[i]->SetTitle(changes[i]);
+  }
+  
+  //##########################################################
+  
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptFit(0);
 
-  for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
-    for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
-      tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
+  gStyle->SetTextSize(0.2);
+  gStyle->SetTitleSize(0.05,"xyz");
+  gStyle->SetLabelOffset(0.01, "xyz");
 
-  //normalize
-  tracks->Scale(1.0 / tracks->Integral());
 
-  return tracks;
+  gStyle->SetTitleOffset(1.2, "y");
+  gStyle->SetTitleOffset(1.2, "x");
+  gStyle->SetEndErrorSize(0.0);
+
+  //##############################################
+
+  //making canvas and pads
+  TCanvas *c = new TCanvas(Form("vertex_reco_syst_%s", plot), Form("vertex_reco_syst_%s", plot),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();
+  
+  
+  TH2F* null = new TH2F("","",100,-1.5,1.5,100,0.9601,1.099);
+  null->SetXTitle("#eta");
+  null->GetXaxis()->CenterTitle(kTRUE);
+  null->SetYTitle("dN/d#eta / dN/d#eta_{pythia}");
+  null->GetYaxis()->CenterTitle(kTRUE);
+  null->Draw();
+
+  for (Int_t i=1; i<7; i++) 
+    ratios[i]->Draw("same");
+
+  TLegend* legend = new TLegend(0.6, 0.6, 0.95, 0.95);
+  legend->SetFillColor(0);
+
+  TLatex* text[7];
+  for (Int_t i=1; i<7; i++) {
+    legend->AddEntry(dNdEta[i], descr[i]);
+  }
+
+  legend->Draw();
+  //text(0.2,0.88,"Effect of changing",0.045,1,kTRUE);
+  //text(0.2,0.83,"relative cross-sections",0.045,1,kTRUE);
+  //text(0.2,0.78,"(vertex reconstruction corr.)",0.043,13,kTRUE);
+
+  c->SaveAs(Form("%s.gif", c->GetName()));
+  c->SaveAs(Form("%s.eps", c->GetName()));
 }
 
-TH1F* Sigma2VertexGaussian()
-{
-  TH2F* tracks = Sigma2VertexGaussianTracksHist();
 
-  TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
-  canvas->Divide(2, 2);
 
-  canvas->cd(1);
-  tracks->Draw("COLZ");
+DrawTriggerEfficiency(Char_t* fileName) {
 
-  TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 10, 0.25, 5.25);
-  for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
-    ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
-  ratio->SetMarkerStyle(21);
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  gStyle->SetOptFit(0);
 
-  canvas->cd(2);
-  ratio->DrawCopy("P");
+  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");
 
-  TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 10, 0.25, 5.25);
-  Double_t sigma3 = Sigma2VertexCount(tracks, 3);
-  for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
-    ratio2->Fill(nSigma, sigma3 / Sigma2VertexCount(tracks, nSigma));
-  ratio2->SetMarkerStyle(21);
+  gStyle->SetTitleOffset(1.1, "y");
+  gStyle->SetTitleOffset(1.1, "x");
+  gStyle->SetEndErrorSize(0.0);
 
-  canvas->cd(3);
-  ratio2->DrawCopy("P");
+  //##############################################
+
+  //making canvas and pads
+  TCanvas *c = new TCanvas("trigger_eff", "",600,500);
 
-  canvas->SaveAs("Sigma2Vertex.eps");
+  TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
 
-  return ratio2;
+  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");
 }
 
-TH1F* Sigma2VertexSimulation()
+void changePtSpectrum(const char* fileName = "analysis_mc.root")
 {
-  TFile* file = TFile::Open("systematics.root");
+  TFile* file = TFile::Open(fileName);
+  TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt"));
 
-  TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex"));
-  if (!sigmavertex)
+  //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"));
+
+  Float_t ptCutOff = 0.3;
+
+  for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
   {
-    printf("Could not read histogram\n");
-    return;
+    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 * 0.3);
+
+      // 110% at pt = 0, ...
+      scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
+    }
+    scale1->SetBinError(i, 0);
+    scale2->SetBinError(i, 0);
   }
 
-  TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included #sigma / % included n #sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
+  new TCanvas;
 
-  for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
-    ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i));
+  scale1->Draw();
+  scale2->SetLineColor(kRed);
+  scale2->Draw("SAME");
 
-  TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
-  canvas->Divide(2, 1);
+  clone1->Multiply(scale1);
+  clone2->Multiply(scale2);
 
-  canvas->cd(1);
-  sigmavertex->SetMarkerStyle(21);
-  sigmavertex->Draw("P");
+  Prepare1DPlot(hist);
+  Prepare1DPlot(clone1);
+  Prepare1DPlot(clone2);
 
-  canvas->cd(2);
-  ratio->SetMarkerStyle(21);
-  ratio->DrawCopy("P");
+  /*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.7);
+  hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
+  hist->GetYaxis()->SetTitleOffset(1);
+
+  TCanvas* canvas = new TCanvas;
+  gPad->SetBottomMargin(0.12);
+  hist->Draw("H");
+  clone1->SetLineColor(kRed);
+  clone1->Draw("HSAME");
+  clone2->SetLineColor(kBlue);
+  clone2->Draw("HSAME");
+  hist->Draw("HSAME");
 
-  return ratio;
+  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 Sigma2VertexCompare()
+void FractionBelowPt()
 {
-  TH1F* ratio1 = Sigma2VertexGaussian();
-  TH1F* ratio2 = Sigma2VertexSimulation();
+  gSystem->Load("libPWG0base");
 
-  ratio1->SetStats(kFALSE);
-  ratio2->SetStats(kFALSE);
+  AlidNdEtaCorrection* fdNdEtaCorrection[4];
 
-  ratio1->SetMarkerStyle(0);
-  ratio2->SetMarkerStyle(0);
+  TFile::Open("systematics.root");
 
-  ratio1->SetLineWidth(2);
-  ratio2->SetLineWidth(2);
+  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);
+  }
 
-  TLegend* legend = new TLegend(0.7, 0.8, 0.95, 0.95);
-  legend->SetFillColor(0);
-  legend->AddEntry(ratio1, "Gaussian");
-  legend->AddEntry(ratio2, "Simulation");
+  Double_t geneCount[5];
+  Double_t measCount[5];
+  geneCount[4] = 0;
+  measCount[4] = 0;
 
-  ratio2->SetTitle("");
-  ratio2->GetYaxis()->SetTitleOffset(1.5);
-  ratio2->GetXaxis()->SetRangeUser(2, 4);
+  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));
 
-  TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
-  InitPad();
+    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));
 
-  ratio2->SetLineColor(kRed);
-  ratio2->Draw();
-  ratio1->Draw("SAME");
+    geneCount[4] += geneCount[i];
+    measCount[4] += measCount[i];
 
-  legend->Draw();
+    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]);
 
-  canvas->SaveAs("Sigma2VertexCompare.eps");
+  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]);
 }
 
-void drawSystematics()
+
+mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
+                                            Char_t* misalignedFile = "correction_map_misaligned.root",
+                                            Char_t* outputFileName="correction_map_misaligned_single.root")
 {
-  //Secondaries();
-  //DrawDifferentSpecies();
-  //Composition();
+  //
+  // 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
 
-  Sigma2VertexSimulation();
+  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 DrawdNdEtaDifferences()
+
+void DrawdNdEtaDifferencesAlignment()
 {
   TH1* hists[5];
 
-  TLegend* legend = new TLegend(0.6, 0.73, 0.98, 0.98);
+  TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
   legend->SetFillColor(0);
 
   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
@@ -770,55 +1414,225 @@ void DrawdNdEtaDifferences()
 
   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;
+      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;
     }
 
-    hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
-    hists[i]->SetTitle("a)");
-
-    hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
-    hists[i]->SetLineColor(i+1);
-    hists[i]->SetMarkerColor(i+1);
-    hists[i]->GetXaxis()->SetLabelOffset(0.015);
-    Prepare1DPlot(hists[i], kFALSE);
-    hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
-
-    legend->AddEntry(hists[i], title);
-    hists[i]->SetTitle(title);
+    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.73, 0.73, 0.98, 0.98);
+  TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
   legend2->SetFillColor(0);
 
   for (Int_t i=1; i<5; ++i)
   {
-    legend2->AddEntry(hists[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" : ""));
+    }
+  }
 
-    hists[i]->Divide(hists[0]);
-    hists[i]->SetTitle("b)");
-    hists[i]->GetYaxis()->SetRangeUser(0.98, 1.02);
-    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);
   }
 
-  legend2->Draw();
+  TList list;
+  for (Int_t i=1; i<4; ++i)
+    list.Add(fdNdEtaCorrectionLow[i]);
 
-  canvas->SaveAs("particlecomposition_result.eps");
-  canvas->SaveAs("particlecomposition_result.gif");
+  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();
 }