more plots
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 6 Nov 2006 08:30:53 +0000 (08:30 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 6 Nov 2006 08:30:53 +0000 (08:30 +0000)
PWG0/dNdEta/drawPlots.C
PWG0/dNdEta/drawSystematics.C
PWG0/dNdEta/makeCorrection2.C
PWG0/dNdEta/runMultiplicitySelector.C

index 6bd6047..4203705 100644 (file)
@@ -94,6 +94,114 @@ void InitPadCOLZ()
   gPad->SetGridy();
 }
 
+void ComparedNdEta(const char* ESDfolder = "dndeta", const char* MCfolder = "dndeta", const char* esdFile = "analysis_esd.root", const char* mcFile = "analysis_mc.root")
+{
+  gSystem->Load("libPWG0base");
+
+  TFile::Open(esdFile);
+  dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
+  fdNdEtaAnalysisESD->LoadHistograms();
+
+  TFile::Open(mcFile);
+  dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(MCfolder, MCfolder);
+  fdNdEtaAnalysisMC->LoadHistograms();
+  fdNdEtaAnalysisMC->Finish(0, 0.3);
+
+  for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
+    fdNdEtaAnalysisESD->GetdNdEtaHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaHistogram(i));
+
+  fdNdEtaAnalysisESD->DrawHistograms();
+}
+
+void CompareVertexDist(Int_t plot = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
+{
+  gSystem->Load("libPWG0base");
+
+  const char* ESDfolder = 0;
+
+  if (plot == 0) // all
+    ESDfolder = "dndeta";
+  else if (plot == 1) // mb
+    ESDfolder = "dndeta_mb";
+  else if (plot == 2) // mb vtx
+    ESDfolder = "dndeta_mbvtx";
+
+  TFile::Open("analysis_esd.root");
+  dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis(ESDfolder, ESDfolder);
+  fdNdEtaAnalysisESD->LoadHistograms();
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+  dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
+
+  TH2F* hist = 0;
+
+  if (plot == 0) // all
+    hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetGeneratedHistogram();
+  else if (plot == 1) // mb
+    hist = dNdEtaCorrection->GetTriggerBiasCorrection()->GetMeasuredHistogram();
+  else if (plot == 2) // mb vtx
+    hist = dNdEtaCorrection->GetVertexRecoCorrection()->GetMeasuredHistogram();
+
+  TH1* proj = hist->ProjectionX();
+
+  TH1* vertex = fdNdEtaAnalysisESD->GetVtxHistogram();
+  for (Int_t i=1; i<=vertex->GetNbinsX(); ++i)
+  {
+    Float_t value = proj->GetBinContent(proj->FindBin(vertex->GetBinCenter(i)));
+    if (value != 0)
+    {
+      printf("vtx = %f, esd = %f, corr = %f, ratio = %f\n", vertex->GetBinCenter(i), vertex->GetBinContent(i), value, vertex->GetBinContent(i) / value);
+      vertex->SetBinContent(i, vertex->GetBinContent(i) / value);
+    }
+  }
+
+  new TCanvas;
+  vertex->DrawCopy();
+}
+
+void CompareTrack2ParticleWithAnalysisData(const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction")
+{
+  gSystem->Load("libPWG0base");
+
+  TFile::Open("analysis_esd.root");
+  dNdEtaAnalysis* fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
+  fdNdEtaAnalysisESD->LoadHistograms();
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+  dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
+
+  //TH1* histESD = fdNdEtaAnalysisESD->GetUncorrectedHistogram();
+  //TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
+
+  TH1* histESD = fdNdEtaAnalysisESD->GetHistogram();
+  TH1* histCorr = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
+
+  TH1F* diff = new TH1F("diff", "diff", 100, 0.98, 1.02);
+
+  new TCanvas;
+  histESD->Draw();
+
+  new TCanvas;
+  histCorr->Draw();
+
+  for (Int_t x=1; x<=histESD->GetNbinsX(); ++x)
+    for (Int_t y=1; y<=histESD->GetNbinsY(); ++y)
+      for (Int_t z=1; z<=histESD->GetNbinsZ(); ++z)
+      {
+        Float_t value1 = histESD->GetBinContent(x, y, z);
+        Float_t value2 = histCorr->GetBinContent(histCorr->FindBin(histESD->GetXaxis()->GetBinCenter(x), histESD->GetYaxis()->GetBinCenter(y), histESD->GetZaxis()->GetBinCenter(z)));
+
+        if (value2 > 0 && value1 > 0)
+        {
+          printf("%f %f %f\n", value1, value2, value1 / value2);
+          diff->Fill(value1 / value2);
+        }
+      }
+
+  new TCanvas;
+  diff->Draw();
+}
+
 void dNdEta(Bool_t onlyESD = kFALSE)
 {
   TFile* file = TFile::Open("analysis_esd.root");
@@ -121,14 +229,14 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   histESDMB->SetMarkerStyle(21);
   histESDMBVtx->SetMarkerStyle(22);
 
-  TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 100, 0, histESDMBVtx->GetMaximum() * 1.1);
+  TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0.1, histESDMBVtx->GetMaximum() * 1.1);
   Prepare1DPlot(dummy);
   dummy->SetStats(kFALSE);
   dummy->SetXTitle("#eta");
   dummy->SetYTitle("dN_{ch}/d#eta");
   dummy->GetYaxis()->SetTitleOffset(1);
 
-  Float_t etaLimit = 1.1999;
+  Float_t etaLimit = 0.7999;
 
   histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
   histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
@@ -178,11 +286,38 @@ void dNdEta(Bool_t onlyESD = kFALSE)
   canvas2->SaveAs("dNdEta2.gif");
   canvas2->SaveAs("dNdEta2.eps");
 
-  TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 500);
-  //InitPad();
-  gPad->SetBottomMargin(0.12);
+  TH1* ratio = (TH1*) histMC->Clone("ratio");
+  TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
+
+  ratio->Divide(histESD);
+  ratioNoPt->Divide(histESDNoPt);
+
+  ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
+
+  ratio->SetLineColor(1);
+  ratioNoPt->SetLineColor(2);
+
+  TCanvas* canvas3 = new TCanvas("dNdEta", "dNdEta", 700, 700);
+  canvas3->Range(0, 0, 1, 1);
+  //canvas3->Divide(1, 2, 0, 0);
+
+  //canvas3->cd(1);
+  TPad* pad1 = new TPad("dNdEta_1", "", 0, 0.4, 0.98, 0.98);
+  pad1->Draw();
+
+  TPad* pad2 = new TPad("dNdEta_2", "", 0, 0.02, 0.98, 0.4);
+  pad2->Draw();
+
+  pad1->SetRightMargin(0.05);
+  pad2->SetRightMargin(0.05);
+
+  // no border between them
+  pad1->SetBottomMargin(0);
+  pad2->SetTopMargin(0);
 
-  TLegend* legend = new TLegend(0.35, 0.2, 0.6, 0.4);
+  pad1->cd();
+
+  TLegend* legend = new TLegend(0.4, 0.2, 0.65, 0.4);
   legend->SetFillColor(0);
   legend->AddEntry(histESDMBVtx, "triggered, vertex");
   legend->AddEntry(histESDMB, "triggered");
@@ -197,21 +332,31 @@ void dNdEta(Bool_t onlyESD = kFALSE)
 
   legend->Draw();
 
-  canvas3->SaveAs("dNdEta.gif");
-  canvas3->SaveAs("dNdEta.eps");
+  pad2->cd();
+  pad2->SetBottomMargin(0.15);
 
-  TH1* ratio = (TH1*) histMC->Clone("ratio");
-  TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
+  TH1F dummy3("dummy3", ";#eta;Ratio: MC / ESD", 1, -1.5, 1.5);
+  dummy3.SetStats(kFALSE);
+  dummy3.SetBinContent(1, 1);
+  dummy3.GetYaxis()->SetRangeUser(0.961, 1.049);
+  dummy3.SetLineWidth(2);
+  dummy3.GetXaxis()->SetLabelSize(0.06);
+  dummy3.GetYaxis()->SetLabelSize(0.06);
+  dummy3.GetXaxis()->SetTitleSize(0.06);
+  dummy3.GetYaxis()->SetTitleSize(0.06);
+  dummy3.GetYaxis()->SetTitleOffset(0.7);
+  dummy3.DrawCopy();
 
-  ratio->Divide(histESD);
-  ratioNoPt->Divide(histESDNoPt);
+  ratio->Draw("SAME");
 
-  TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
+  //pad2->Draw();
 
-  ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
+  canvas3->Modified();
 
-  ratio->SetLineColor(1);
-  ratioNoPt->SetLineColor(2);
+  canvas3->SaveAs("dNdEta.gif");
+  canvas3->SaveAs("dNdEta.eps");
+
+  TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
 
   ratio->Draw();
   ratioNoPt->Draw("SAME");
@@ -382,9 +527,13 @@ void TriggerBias1D(const char* fileName = "correction_map.root", const char* fol
   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderName, folderName);
   dNdEtaCorrection->LoadHistograms(fileName, folderName);
 
-  TH1* hist = dNdEtaCorrection->GetTriggerCorrection()->Get1DCorrection("x");
+  TH1* hist = dNdEtaCorrection->GetTriggerBiasCorrection()->Get1DCorrection("x");
+  TH1* hist2 = dNdEtaCorrection->GetTriggerBiasCorrection()->Get1DCorrection("y", -10, 10);
+
+  TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
+  canvas->Divide(2, 1);
 
-  TCanvas* canvas = new TCanvas("TriggerBias1D", "TriggerBias1D", 500, 500);
+  canvas->cd(1);
   InitPad();
 
   Prepare1DPlot(hist);
@@ -394,6 +543,22 @@ void TriggerBias1D(const char* fileName = "correction_map.root", const char* fol
   hist->GetYaxis()->SetTitleOffset(1.6);
   hist->Draw();
 
+  canvas->cd(2);
+  InitPad();
+
+  Prepare1DPlot(hist2);
+  hist2->SetTitle("");
+  hist2->GetYaxis()->SetTitle("correction factor");
+  hist2->GetXaxis()->SetRangeUser(0, 5);
+  hist2->GetYaxis()->SetTitleOffset(1.6);
+  hist2->GetXaxis()->SetTitle("multiplicity");
+  hist2->Draw();
+
+  TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
+  pave->SetFillColor(0);
+  pave->AddText("|z| < 10 cm");
+  pave->Draw();
+
   canvas->SaveAs("TriggerBias1D.eps");
 }
 
@@ -432,8 +597,12 @@ void VtxRecon1D(const char* fileName = "correction_map.root", const char* folder
   dNdEtaCorrection->LoadHistograms(fileName, folderName);
 
   TH1* hist = dNdEtaCorrection->GetVertexRecoCorrection()->Get1DCorrection("x");
+  TH1* hist2 = dNdEtaCorrection->GetVertexRecoCorrection()->Get1DCorrection("y", -10, 10);
+
+  TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 1000, 500);
+  canvas->Divide(2, 1);
 
-  TCanvas* canvas = new TCanvas("VtxRecon1D", "VtxRecon1D", 500, 500);
+  canvas->cd(1);
   InitPad();
 
   Prepare1DPlot(hist);
@@ -441,9 +610,34 @@ void VtxRecon1D(const char* fileName = "correction_map.root", const char* folder
   hist->GetYaxis()->SetTitle("correction factor");
   hist->GetYaxis()->SetRangeUser(1, 1.8);
   hist->GetYaxis()->SetTitleOffset(1.6);
-  hist->Draw();
+  hist->DrawCopy();
+
+  canvas->cd(2);
+  InitPad();
+
+  Prepare1DPlot(hist2);
+  hist2->SetTitle("");
+  hist2->GetYaxis()->SetTitle("correction factor");
+  hist2->GetXaxis()->SetRangeUser(0, 20);
+  hist2->GetYaxis()->SetTitleOffset(1.6);
+  hist2->GetXaxis()->SetTitle("multiplicity");
+  hist2->Draw();
+
+  TPaveText* pave = new TPaveText(0.6, 0.8, 0.8, 0.85, "NDC");
+  pave->SetFillColor(0);
+  pave->AddText("|z| < 10 cm");
+  pave->Draw();
 
   canvas->SaveAs("VtxRecon1D.eps");
+
+  for (Int_t i=1; i<=hist->GetNbinsX() / 2; ++i)
+    if (hist->GetBinContent(hist->GetNbinsX() + 1 - i) != 0)
+      hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinContent(hist->GetNbinsX() + 1 - i));
+
+  new TCanvas;
+  hist->GetXaxis()->SetRange(1, hist->GetNbinsX() / 2);
+  hist->GetYaxis()->SetRangeUser(0.8, 1.2);
+  hist->DrawCopy();
 }
 
 void Track2ParticleAsNumber(const char* fileName = "correction_map.root")
@@ -595,66 +789,78 @@ void CompareTrack2Particle1D(Float_t upperPtLimit = 9.9)
 {
   gSystem->Load("libPWG0base");
 
-  Track2Particle1DCreatePlots("correction_maponly-positive.root", "dndeta_correction", upperPtLimit);
+  // particle type
+  for (Int_t particle=0; particle<4; ++particle)
+  {
+    TString dirName;
+    dirName.Form("correction_%d", particle);
+    Track2Particle1DCreatePlots("systematics-detail-only-positive.root", dirName, upperPtLimit);
 
-  TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")->Clone("pos_x"));
-  TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")->Clone("pos_y"));
-  TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")->Clone("pos_z"));
+    TString tmpx, tmpy, tmpz;
+    tmpx.Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", dirName.Data(), dirName.Data());
+    tmpy.Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", dirName.Data(), dirName.Data());
+    tmpz.Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", dirName.Data(), dirName.Data());
 
-  Track2Particle1DCreatePlots("correction_maponly-negative.root", "dndeta_correction", upperPtLimit);
+    TH1* posX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("pos_x"));
+    TH1* posY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("pos_y"));
+    TH1* posZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("pos_z"));
 
-  TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x")->Clone("neg_x"));
-  TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y")->Clone("neg_y"));
-  TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z")->Clone("neg_z"));
+    Track2Particle1DCreatePlots("systematics-detail-only-negative.root", dirName, upperPtLimit);
 
-  //printf("%f %f %f %f\n", posX->GetBinContent(20), posX->GetBinError(20), negX->GetBinContent(20), negX->GetBinError(20));
+    TH1* negX = dynamic_cast<TH1*> (gROOT->FindObject(tmpx)->Clone("neg_x"));
+    TH1* negY = dynamic_cast<TH1*> (gROOT->FindObject(tmpy)->Clone("neg_y"));
+    TH1* negZ = dynamic_cast<TH1*> (gROOT->FindObject(tmpz)->Clone("neg_z"));
 
-  posX->Divide(negX);
-  posY->Divide(negY);
-  posZ->Divide(negZ);
+    posX->Divide(negX);
+    posY->Divide(negY);
+    posZ->Divide(negZ);
 
-  //printf("%f %f\n", posX->GetBinContent(20), posX->GetBinError(20));
+    Prepare1DPlot(posX);
+    Prepare1DPlot(posY);
+    Prepare1DPlot(posZ);
 
-  Prepare1DPlot(posX);
-  Prepare1DPlot(posY);
-  Prepare1DPlot(posZ);
+    Float_t min = 0.8;
+    Float_t max = 1.2;
 
-  Float_t min = 0.8;
-  Float_t max = 1.2;
+    posX->SetMinimum(min);
+    posX->SetMaximum(max);
+    posY->SetMinimum(min);
+    posY->SetMaximum(max);
+    posZ->SetMinimum(min);
+    posZ->SetMaximum(max);
 
-  posX->SetMinimum(min);
-  posX->SetMaximum(max);
-  posY->SetMinimum(min);
-  posY->SetMaximum(max);
-  posZ->SetMinimum(min);
-  posZ->SetMaximum(max);
+    posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 
-  posZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
+    posX->GetYaxis()->SetTitleOffset(1.7);
+    posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
+    posY->GetYaxis()->SetTitleOffset(1.7);
+    posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
+    posZ->GetYaxis()->SetTitleOffset(1.7);
+    posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
 
-  posX->GetYaxis()->SetTitleOffset(1.7);
-  posX->GetYaxis()->SetTitle("C_{+} / C_{-}");
-  posY->GetYaxis()->SetTitleOffset(1.7);
-  posY->GetYaxis()->SetTitle("C_{+} / C_{-}");
-  posZ->GetYaxis()->SetTitleOffset(1.7);
-  posZ->GetYaxis()->SetTitle("C_{+} / C_{-}");
+    posZ->GetXaxis()->SetRangeUser(0, 1);
 
-  TCanvas* canvas = new TCanvas("CompareTrack2Particle1D", "CompareTrack2Particle1D", 1200, 400);
-  canvas->Divide(3, 1);
+    TString canvasName;
+    canvasName.Form("PosNegRatios_%s_%f", ((particle == 0) ? "Pi" : ((particle == 1) ? "K" : ((particle == 2) ? "p" : "other"))), upperPtLimit);
 
-  canvas->cd(1);
-  InitPad();
-  posX->Draw();
+    TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
+    canvas->Divide(3, 1);
 
-  canvas->cd(2);
-  InitPad();
-  posY->Draw();
+    canvas->cd(1);
+    InitPad();
+    posX->DrawCopy();
 
-  canvas->cd(3);
-  InitPad();
-  posZ->Draw();
+    canvas->cd(2);
+    InitPad();
+    posY->DrawCopy();
+
+    canvas->cd(3);
+    InitPad();
+    posZ->DrawCopy();
 
-  canvas->SaveAs(Form("CompareTrack2Particle1D_%f.gif", upperPtLimit));
-  canvas->SaveAs(Form("CompareTrack2Particle1D_%f.eps", upperPtLimit));
+    canvas->SaveAs(Form("%s.gif", canvas->GetName()));
+    canvas->SaveAs(Form("%s.eps", canvas->GetName()));
+  }
 }
 
 void Track2Particle2DCreatePlots(const char* fileName = "correction_map.root")
index 6a8e4dd..7963746 100644 (file)
@@ -23,6 +23,8 @@ 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 SetRanges(TAxis* axis)
 {
@@ -257,6 +259,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 +333,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");
@@ -798,8 +864,9 @@ void DrawdNdEtaDifferences()
       Prepare1DPlot(hists[i], kFALSE);
       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
       hists[i]->GetYaxis()->SetRangeUser(6, 7);
-      hists[i]->SetLineColor(i+1);
-      hists[i]->SetMarkerColor(i+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);
@@ -825,6 +892,8 @@ void DrawdNdEtaDifferences()
 
       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);
@@ -834,11 +903,27 @@ void DrawdNdEtaDifferences()
 
   legend2->Draw();
 
-  canvas->SaveAs("particlecomposition_result.eps");
-  canvas->SaveAs("particlecomposition_result.gif");
+  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");
 }
 
-mergeCorrections4SystematicStudies(Char_t* standardCorrectionFileName="correction_map.root",
+mergeCorrectionsWithDifferentCrosssections(Char_t* standardCorrectionFileName="correction_map.root",
                                             Char_t* systematicCorrectionFileName="systematics.root",
                                             Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
   //
@@ -1076,7 +1161,7 @@ DrawVertexRecoSyst() {
   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");
 
@@ -1248,11 +1333,11 @@ DrawSpectraPID(Char_t* fileName) {
 
 void changePtSpectrum()
 {
-  TFile* file = TFile::Open("pt.root");
+  TFile* file = TFile::Open("analysis_mc.root");
   TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt"));
 
-  hist->Rebin(3);
-  hist->Scale(1.0/3);
+  //hist->Rebin(3);
+  //hist->Scale(1.0/3);
 
   TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
   TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
@@ -1277,9 +1362,12 @@ void changePtSpectrum()
       // 110% at pt = 0, ...
       scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
     }
+    scale1->SetBinError(i, 0);
+    scale2->SetBinError(i, 0);
   }
 
   new TCanvas;
+
   scale1->Draw();
   scale2->SetLineColor(kRed);
   scale2->Draw("SAME");
@@ -1291,17 +1379,23 @@ void changePtSpectrum()
   Prepare1DPlot(clone1);
   Prepare1DPlot(clone2);
 
+  /*hist->SetMarkerStyle(markers[0]);
+  clone1->SetMarkerStyle(markers[0]);
+  clone2->SetMarkerStyle(markers[0]);*/
+
   hist->SetTitle(";p_{T} in GeV/c;dN/dp_{T} in c/GeV");
   hist->GetXaxis()->SetRangeUser(0, 0.7);
-  hist->GetYaxis()->SetRangeUser(0, clone2->GetMaximum() * 1.1);
+  hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
   hist->GetYaxis()->SetTitleOffset(1);
 
   TCanvas* canvas = new TCanvas;
-  hist->Draw();
+  gPad->SetBottomMargin(0.12);
+  hist->Draw("H");
   clone1->SetLineColor(kRed);
-  clone1->Draw("SAME");
+  clone1->Draw("HSAME");
   clone2->SetLineColor(kBlue);
-  clone2->Draw("SAME");
+  clone2->Draw("HSAME");
+  hist->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());
@@ -1311,6 +1405,7 @@ void changePtSpectrum()
   printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
 
   canvas->SaveAs("changePtSpectrum.gif");
+  canvas->SaveAs("changePtSpectrum.eps");
 }
 
 void FractionBelowPt()
@@ -1354,3 +1449,245 @@ void FractionBelowPt()
 
   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();
+}
index c248439..21521ac 100644 (file)
 #include "../CreateESDChain.C"
 #include "../PWG0Helper.C"
 
-void makeCorrection2(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t debug = kFALSE, Bool_t aProof = kFALSE, const Char_t* option = "")
+void makeCorrection2(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t debug = kFALSE, Bool_t aProof = kFALSE, const Char_t* option = "", const Char_t* proofConnect = "proof01@lxb6046")
 {
   if (aProof)
-    connectProof("proof01@lxb6046");
+    connectProof(proofConnect);
 
   TString libraries("libEG;libGeom;libESD;libPWG0base;libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6");
   TString packages("PWG0base;PWG0dep");
@@ -41,3 +41,49 @@ void makeCorrection2(Char_t* dataDir, Int_t nRuns=20, Int_t offset = 0, Bool_t d
 
   Int_t result = executeQuery(chain, &inputList, selector, option);
 }
+
+void VerifyCorrection(Int_t draw = 0, const char* correctionMapFile = "correction_map.root", const char* correctionMapFolder = "dndeta_correction", const char* analysisDir = "dndetaESD", const char* checkDir = "dndetaMC")
+{
+  gSystem->Load("libPWG0base");
+
+  AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
+  dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
+  dNdEtaCorrection->SetName("dndeta_correction");
+  dNdEtaCorrection->SetTitle("dndeta_correction");
+
+  TFile::Open(correctionMapFile);
+
+  dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(analysisDir, analysisDir);
+  fdNdEtaAnalysis->LoadHistograms();
+
+  // correct with track2particle
+  TH3F* hist = fdNdEtaAnalysis->GetHistogram();
+  for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
+    for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
+      for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
+      {
+        Float_t correction = dNdEtaCorrection->GetTrack2ParticleCorrection(hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
+        Float_t value = hist->GetBinContent(x, y, z);
+
+        hist->SetBinContent(x, y, z, correction * value);
+      }
+
+  fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
+
+  dNdEtaAnalysis* fdNdEtaAnalysisMC = new dNdEtaAnalysis(checkDir, checkDir);
+  fdNdEtaAnalysisMC->LoadHistograms();
+
+  fdNdEtaAnalysisMC->Finish(0, 0.3);
+
+  if (draw == 0)
+    fdNdEtaAnalysis->DrawHistograms();
+  else if (draw == 1)
+    fdNdEtaAnalysisMC->DrawHistograms();
+  else
+  {
+    for (Int_t i=0; i<dNdEtaAnalysis::kVertexBinning; ++i)
+      fdNdEtaAnalysis->GetdNdEtaHistogram(i)->Divide(fdNdEtaAnalysisMC->GetdNdEtaHistogram(i));
+
+    fdNdEtaAnalysis->DrawHistograms();
+  }
+}
index 68b78ca..8f376d0 100644 (file)
@@ -7,7 +7,7 @@
 #include "../CreateESDChain.C"
 #include "../PWG0Helper.C"
 
-TChain* runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE)
+void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE)
 {
   if (aProof)
     connectProof("proof01@lxb6046");