PWG0/dNdEta/drawSystematics.C

   1 /* $Id$ */
   2
   3 #if !defined(__CINT__) || defined(__MAKECINT__)
   4
   5 #include "AliPWG0Helper.h"
   6 #include "dNdEtaAnalysis.h"
   7 #include "AlidNdEtaCorrection.h"
   8
   9 #include <TCanvas.h>
  10 #include <TFile.h>
  11 #include <TH1.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TLine.h>
  15 #include <TSystem.h>
  16 #include <TLegend.h>
  17 #include <TPad.h>
  18 #include <TF1.h>
  19
  20 extern TPad* gPad;
  21
  22 void Track2Particle1DCreatePlots(const char* fileName = "correction_map.root", const char* folderName = "dndeta_correction", Float_t upperPtLimit = 10);
  23
  24 #endif
  25
  26
  27 void SetRanges(TAxis* axis)
  28 {
  29   if (strcmp(axis->GetTitle(), "#eta") == 0)
  30     axis->SetRangeUser(-1.7999, 1.7999);
  31   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
  32     axis->SetRangeUser(0, 9.9999);
  33   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
  34     axis->SetRangeUser(-15, 14.9999);
  35   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
  36     axis->SetRangeUser(0, 99.9999);
  37 }
  38
  39 void SetRanges(TH1* hist)
  40 {
  41   SetRanges(hist->GetXaxis());
  42   SetRanges(hist->GetYaxis());
  43   SetRanges(hist->GetZaxis());
  44 }
  45
  46 void Prepare3DPlot(TH3* hist)
  47 {
  48   hist->GetXaxis()->SetTitleOffset(1.5);
  49   hist->GetYaxis()->SetTitleOffset(1.5);
  50   hist->GetZaxis()->SetTitleOffset(1.5);
  51
  52   hist->SetStats(kFALSE);
  53 }
  54
  55 void Prepare2DPlot(TH2* hist)
  56 {
  57   hist->SetStats(kFALSE);
  58   hist->GetYaxis()->SetTitleOffset(1.4);
  59
  60   SetRanges(hist);
  61 }
  62
  63 void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
  64 {
  65   hist->SetLineWidth(2);
  66   hist->SetStats(kFALSE);
  67
  68   hist->GetXaxis()->SetTitleOffset(1.2);
  69   hist->GetYaxis()->SetTitleOffset(1.2);
  70
  71   if (setRanges)
  72     SetRanges(hist);
  73 }
  74
  75 void InitPad()
  76 {
  77   if (!gPad)
  78     return;
  79
  80   gPad->Range(0, 0, 1, 1);
  81   gPad->SetLeftMargin(0.15);
  82   //gPad->SetRightMargin(0.05);
  83   //gPad->SetTopMargin(0.13);
  84   //gPad->SetBottomMargin(0.1);
  85
  86   //gPad->SetGridx();
  87   //gPad->SetGridy();
  88 }
  89
  90 void InitPadCOLZ()
  91 {
  92   gPad->Range(0, 0, 1, 1);
  93   gPad->SetRightMargin(0.15);
  94   gPad->SetLeftMargin(0.12);
  95
  96   gPad->SetGridx();
  97   gPad->SetGridy();
  98 }
  99
 100 void Secondaries()
 101 {
 102   TFile* file = TFile::Open("systematics.root");
 103
 104   TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
 105   if (!secondaries)
 106   {
 107     printf("Could not read histogram\n");
 108     return;
 109   }
 110
 111   TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
 112   canvas->Divide(3, 3);
 113   for (Int_t i=1; i<=8; i++)
 114   {
 115     TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
 116     hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
 117
 118     canvas->cd(i);
 119     hist->Draw();
 120   }
 121 }
 122
 123 void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
 124 {
 125   gSystem->Load("libPWG0base");
 126
 127   TString canvasName;
 128   canvasName.Form("Track2Particle1DComposition");
 129   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 130   canvas->Divide(3, 1);
 131
 132   TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
 133
 134   for (Int_t i=0; i<folderCount; ++i)
 135   {
 136     Track2Particle1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
 137
 138     TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_x_div_meas_%s_nTrackToNPart_x", folderNames[i], folderNames[i])));
 139     TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_y_div_meas_%s_nTrackToNPart_y", folderNames[i], folderNames[i])));
 140     TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject(Form("gene_%s_nTrackToNPart_z_div_meas_%s_nTrackToNPart_z", folderNames[i], folderNames[i])));
 141
 142     Prepare1DPlot(corrX);
 143     Prepare1DPlot(corrY);
 144     Prepare1DPlot(corrZ);
 145
 146     const char* title = "Track2Particle Correction";
 147     corrX->SetTitle(title);
 148     corrY->SetTitle(title);
 149     corrZ->SetTitle(title);
 150
 151     corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 152
 153     corrX->SetLineColor(i+1);
 154     corrY->SetLineColor(i+1);
 155     corrZ->SetLineColor(i+1);
 156
 157     canvas->cd(1);
 158     InitPad();
 159     corrX->DrawCopy(((i>0) ? "SAME" : ""));
 160
 161     canvas->cd(2);
 162     InitPad();
 163     corrY->DrawCopy(((i>0) ? "SAME" : ""));
 164
 165     canvas->cd(3);
 166     InitPad();
 167     corrZ->DrawCopy(((i>0) ? "SAME" : ""));
 168
 169     legend->AddEntry(corrZ, folderNames[i]);
 170   }
 171
 172   legend->Draw();
 173
 174   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
 175   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
 176 }
 177
 178 TH1** DrawRatios(const char* fileName = "systematics.root")
 179 {
 180   gSystem->Load("libPWG0base");
 181
 182   TFile* file = TFile::Open(fileName);
 183
 184   TString canvasName;
 185   canvasName.Form("DrawRatios");
 186   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
 187   canvas->Divide(2, 1);
 188   canvas->cd(1);
 189
 190   TH1** ptDists = new TH1*[4];
 191
 192   TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
 193
 194   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
 195   const char* particleNames[] = { "#pi", "K", "p", "other" };
 196   for (Int_t i=0; i<4; ++i)
 197   {
 198     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]);
 199     dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]);
 200
 201     TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 202
 203     gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 204     gene->GetXaxis()->SetRangeUser(-10, 10);
 205
 206     ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i])));
 207     ptDists[i]->SetTitle(";p_{T};Count");
 208     if (!ptDists[i])
 209     {
 210       printf("Problem getting distribution %d.\n", i);
 211       return 0;
 212     }
 213
 214     ptDists[i]->GetYaxis()->SetRangeUser(1, ptDists[i]->GetMaximum()*1.1);
 215     ptDists[i]->GetXaxis()->SetRangeUser(0, 9.9);
 216     ptDists[i]->SetLineColor(i+1);
 217     ptDists[i]->DrawCopy((i == 0) ? "" : "SAME");
 218     ptDists[i]->GetYaxis()->SetRange(0, 0);
 219
 220     legend->AddEntry(ptDists[i], particleNames[i]);
 221   }
 222   gPad->SetLogy();
 223
 224   TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total"));
 225
 226   for (Int_t i=1; i<4; ++i)
 227     total->Add(ptDists[i]);
 228
 229   canvas->cd(2);
 230   for (Int_t i=0; i<4; ++i)
 231   {
 232     ptDists[i]->Divide(total);
 233     ptDists[i]->SetStats(kFALSE);
 234     ptDists[i]->SetTitle(";p_{T};Fraction of total");
 235     ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
 236     ptDists[i]->Draw((i == 0) ? "" : "SAME");
 237   }
 238   legend->SetFillColor(0);
 239   legend->Draw();
 240
 241   canvas->SaveAs("DrawRatios.gif");
 242
 243
 244   canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500);
 245   for (Int_t i=0; i<4; ++i)
 246   {
 247     TH1* hist = ptDists[i]->Clone();
 248     hist->GetXaxis()->SetRangeUser(0, 1.9);
 249     hist->Draw((i == 0) ? "" : "SAME");
 250   }
 251   legend->Draw();
 252
 253   canvas->SaveAs("pythiaratios.eps");
 254
 255   file->Close();
 256
 257   return ptDists;
 258 }
 259
 260 void DrawCompareToReal()
 261 {
 262   gROOT->ProcessLine(".L drawPlots.C");
 263
 264   const char* fileNames[] = { "correction_map.root", "new_compositions.root" };
 265   const char* folderNames[] = { "dndeta_correction", "PythiaRatios" };
 266
 267   Track2Particle1DComposition(fileNames, 2, folderNames);
 268 }
 269
 270 void DrawDifferentSpecies()
 271 {
 272   gROOT->ProcessLine(".L drawPlots.C");
 273
 274   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
 275   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
 276
 277   Track2Particle1DComposition(fileNames, 4, folderNames);
 278 }
 279
 280 void DrawSpeciesAndCombination()
 281 {
 282   gROOT->ProcessLine(".L drawPlots.C");
 283
 284   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" };
 285   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" };
 286
 287   Track2Particle1DComposition(fileNames, 4, folderNames);
 288 }
 289
 290 void DrawSimulatedVsCombined()
 291 {
 292   gROOT->ProcessLine(".L drawPlots.C");
 293
 294   const char* fileNames[] = { "new_compositions.root", "new_compositions.root" };
 295   const char* folderNames[] = { "Pythia", "PythiaRatios" };
 296
 297   Track2Particle1DComposition(fileNames, 2, folderNames);
 298 }
 299
 300 void DrawBoosts()
 301 {
 302   gROOT->ProcessLine(".L drawPlots.C");
 303
 304   const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
 305   const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" };
 306
 307   Track2Particle1DComposition(fileNames, 4, folderNames);
 308 }
 309
 310 TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
 311 {
 312   gSystem->Load("libPWG0base");
 313
 314   AlidNdEtaCorrection* fdNdEtaCorrection[2];
 315
 316   TFile::Open(filename);
 317
 318   fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1);
 319   fdNdEtaCorrection[0]->LoadHistograms(filename, folder1);
 320
 321   fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2);
 322   fdNdEtaCorrection[1]->LoadHistograms(filename, folder2);
 323
 324   TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
 325   TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
 326
 327   //TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s / %s;Count", folder2, folder1), 1000, 0.9, 1.1);
 328   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());
 329
 330   for (Int_t x=hist1->GetXaxis()->FindBin(-10); x<=hist1->GetXaxis()->FindBin(10); ++x)
 331     for (Int_t y=hist1->GetYaxis()->FindBin(-0.8); y<=hist1->GetYaxis()->FindBin(0.8); ++y)
 332       for (Int_t z=hist1->GetZaxis()->FindBin(0.3); z<=hist1->GetZaxis()->FindBin(9.9); ++z)
 333         if (hist1->GetBinContent(x, y, z) != 0)
 334           difference->Fill(hist2->GetBinContent(x, y, z) / hist1->GetBinContent(x, y, z), hist1->GetYaxis()->GetBinCenter(y));
 335
 336   difference->GetYaxis()->SetRangeUser(-0.8, 0.8);
 337
 338   printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
 339
 340   return difference;
 341 }
 342
 343 void HistogramDifferences()
 344 {
 345   TH2F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted");
 346   TH2F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted");
 347   TH2F* KReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "KReduced");
 348   TH2F* pReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "pReduced");
 349
 350   TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
 351   canvas->Divide(2, 2);
 352
 353   canvas->cd(1);
 354   KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
 355   KBoosted->Draw("COLZ");
 356
 357   canvas->cd(2);
 358   KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
 359   KReduced->Draw("COLZ");
 360
 361   canvas->cd(3);
 362   pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02);
 363   pBoosted->Draw("COLZ");
 364
 365   canvas->cd(4);
 366   pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02);
 367   pReduced->Draw("COLZ");
 368
 369   canvas->SaveAs("HistogramDifferences.gif");
 370
 371   canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500);
 372   canvas->Divide(2, 1);
 373
 374   canvas->cd(1);
 375   gPad->SetBottomMargin(0.13);
 376   KBoosted->SetTitle("a) Ratio of correction maps");
 377   KBoosted->SetStats(kFALSE);
 378   KBoosted->GetXaxis()->SetTitleOffset(1.4);
 379   KBoosted->GetXaxis()->SetLabelOffset(0.02);
 380   KBoosted->Draw("COLZ");
 381
 382   canvas->cd(2);
 383   gPad->SetGridx();
 384   gPad->SetGridy();
 385
 386   TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
 387
 388   for (Int_t i=0; i<4; ++i)
 389   {
 390     TH2F* hist = 0;
 391     TString name;
 392     switch (i)
 393     {
 394       case 0: hist = KBoosted; name = "K enhanced"; break;
 395       case 1: hist = KReduced; name = "K reduced"; break;
 396       case 2: hist = pBoosted; name = "p enhanced"; break;
 397       case 3: hist = pReduced; name = "p reduced"; break;
 398     }
 399
 400     TProfile* profile = hist->ProfileY();
 401     profile->SetTitle("b) Mean and RMS");
 402     profile->GetXaxis()->SetRange(hist->GetYaxis()->GetFirst(), hist->GetYaxis()->GetLast());
 403     profile->GetXaxis()->SetTitleOffset(1.2);
 404     profile->GetXaxis()->SetLabelOffset(0.02);
 405     profile->GetYaxis()->SetRangeUser(0.98, 1.02);
 406     profile->SetStats(kFALSE);
 407     profile->SetLineColor(i+1);
 408     profile->SetMarkerColor(i+1);
 409     profile->DrawCopy(((i > 0) ? "SAME" : ""));
 410
 411
 412     legend->AddEntry(profile, name);
 413   }
 414
 415   legend->Draw();
 416   canvas->SaveAs("particlecomposition_result.eps");
 417 }
 418
 419
 420 void ScalePtDependent(TH3F* hist, TF1* function)
 421 {
 422   // assumes that pt is the third dimension of hist
 423   // scales with function(pt)
 424
 425   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
 426   {
 427     Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z));
 428     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
 429
 430     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
 431       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
 432         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
 433   }
 434 }
 435
 436 void ScalePtDependent(TH3F* hist, TH1* function)
 437 {
 438   // assumes that pt is the third dimension of hist
 439   // scales with histogram(pt)
 440
 441   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
 442   {
 443     Double_t factor = function->GetBinContent(function->GetXaxis()->FindBin(hist->GetZaxis()->GetBinCenter(z)));
 444     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
 445
 446     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
 447       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
 448         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
 449   }
 450 }
 451
 452 const char* ChangeComposition(void** correctionPointer, Int_t index)
 453 {
 454   AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
 455
 456   switch (index)
 457   {
 458     case 0: // result from pp events
 459       {
 460         TFile::Open("pythiaratios.root");
 461
 462         for (Int_t i=0; i<4; ++i)
 463         {
 464           TString name;
 465           name.Form("correction_%d", i);
 466           fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
 467           fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name);
 468         }
 469       }
 470       return "Pythia";
 471       break;
 472
 473     case 1: // each species rated with pythia ratios
 474       /*TH1** ptDists = DrawRatios("pythiaratios.root");
 475
 476       for (Int_t i=0; i<3; ++i)
 477       {
 478         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
 479         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
 480       }*/
 481       return "PythiaRatios";
 482       break;
 483
 484       // one species enhanced / reduced
 485     case 2: // + 50% pions
 486     case 3: // - 50% pions
 487     case 4: // + 50% kaons
 488     case 5: // - 50% kaons
 489     case 6: // + 50% protons
 490     case 7: // - 50% protons
 491       Int_t correctionIndex = (index - 2) / 2;
 492       Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
 493
 494       fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(scaleFactor);
 495       fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(scaleFactor);
 496
 497       TString* str = new TString;
 498       str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
 499       return str->Data();
 500       break;
 501
 502       // each species rated with pythia ratios
 503     case 12: // + 50% pions
 504     case 13: // - 50% pions
 505     case 14: // + 50% kaons
 506     case 15: // - 50% kaons
 507     case 16: // + 50% protons
 508     case 17: // - 50% protons
 509       TH1** ptDists = DrawRatios("pythiaratios.root");
 510       Int_t functionIndex = (index - 2) / 2;
 511       Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
 512       ptDists[functionIndex]->Scale(scaleFactor);
 513
 514       for (Int_t i=0; i<3; ++i)
 515       {
 516         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
 517         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
 518       }
 519       TString* str = new TString;
 520       str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
 521       return str->Data();
 522       break;
 523
 524     case 999:
 525       TF1* ptDependence = new TF1("simple", "x", 0, 100);
 526       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence);
 527       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence);
 528       break;
 529
 530   }
 531
 532   return "noname";
 533 }
 534
 535 void Composition()
 536 {
 537   gSystem->Load("libPWG0base");
 538
 539   gSystem->Unlink("new_compositions.root");
 540
 541   Int_t nCompositions = 8;
 542   for (Int_t comp = 0; comp < nCompositions; ++comp)
 543   {
 544     AlidNdEtaCorrection* fdNdEtaCorrection[4];
 545
 546     TFile::Open("systematics.root");
 547
 548     for (Int_t i=0; i<4; ++i)
 549     {
 550       TString name;
 551       name.Form("correction_%d", i);
 552       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
 553       fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
 554     }
 555
 556     const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
 557
 558     Double_t geneCount[5];
 559     Double_t measCount[5];
 560     geneCount[4] = 0;
 561     measCount[4] = 0;
 562
 563     for (Int_t i=0; i<4; ++i)
 564     {
 565       geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral();
 566       measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral();
 567
 568       geneCount[4] += geneCount[i];
 569       measCount[4] += measCount[i];
 570
 571       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]);
 572     }
 573
 574     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]);
 575
 576     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]);
 577
 578     TList* collection = new TList;
 579
 580     for (Int_t i=0; i<4; ++i)
 581       collection->Add(fdNdEtaCorrection[i]);
 582
 583     AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(newName, newName);
 584     newComposition->Merge(collection);
 585     newComposition->Finish();
 586
 587     delete collection;
 588
 589     TFile* file = TFile::Open("new_compositions.root", "UPDATE");
 590     newComposition->SaveHistograms();
 591     //file->Write();
 592     file->Close();
 593   }
 594
 595   gROOT->ProcessLine(".L drawPlots.C");
 596
 597   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" };
 598   const char* folderNames[] = { "Pythia", "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };
 599
 600   Track2Particle1DComposition(fileNames, nCompositions, folderNames);
 601 }
 602
 603 Double_t ConvSigma1To2D(Double_t sigma)
 604 {
 605   return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
 606 }
 607
 608 Double_t ConvDistance1DTo2D(Double_t distance)
 609 {
 610   return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
 611 }
 612
 613 Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
 614 {
 615   Double_t count = 0;
 616
 617   //nSigma = ConvSigma1To2D(nSigma);
 618
 619   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
 620     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
 621     {
 622       Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
 623       Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
 624
 625       Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
 626
 627       d = ConvDistance1DTo2D(d);
 628
 629       if (d < nSigma)
 630         count += tracks->GetBinContent(x, y);
 631     }
 632
 633   return count;
 634 }
 635
 636 TH2F* Sigma2VertexGaussianTracksHist()
 637 {
 638   TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
 639
 640   TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
 641   gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
 642
 643   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
 644     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
 645       tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
 646
 647   //normalize
 648   tracks->Scale(1.0 / tracks->Integral());
 649
 650   return tracks;
 651 }
 652
 653 TH1F* Sigma2VertexGaussian()
 654 {
 655   TH2F* tracks = Sigma2VertexGaussianTracksHist();
 656
 657   TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
 658   canvas->Divide(2, 2);
 659
 660   canvas->cd(1);
 661   tracks->Draw("COLZ");
 662
 663   TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 10, 0.25, 5.25);
 664   for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
 665     ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
 666   ratio->SetMarkerStyle(21);
 667
 668   canvas->cd(2);
 669   ratio->DrawCopy("P");
 670
 671   TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 10, 0.25, 5.25);
 672   Double_t sigma3 = Sigma2VertexCount(tracks, 3);
 673   for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
 674     ratio2->Fill(nSigma, sigma3 / Sigma2VertexCount(tracks, nSigma));
 675   ratio2->SetMarkerStyle(21);
 676
 677   canvas->cd(3);
 678   ratio2->DrawCopy("P");
 679
 680   canvas->SaveAs("Sigma2Vertex.eps");
 681
 682   return ratio2;
 683 }
 684
 685 TH1F* Sigma2VertexSimulation()
 686 {
 687   TFile* file = TFile::Open("systematics.root");
 688
 689   TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex"));
 690   if (!sigmavertex)
 691   {
 692     printf("Could not read histogram\n");
 693     return;
 694   }
 695
 696   TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;Nsigma;% included 3 sigma / % included n sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
 697
 698   for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
 699     ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i));
 700
 701   TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
 702   canvas->Divide(2, 1);
 703
 704   canvas->cd(1);
 705   sigmavertex->SetMarkerStyle(21);
 706   sigmavertex->Draw("P");
 707
 708   canvas->cd(2);
 709   ratio->SetMarkerStyle(21);
 710   ratio->DrawCopy("P");
 711
 712   return ratio;
 713 }
 714
 715 void Sigma2VertexCompare()
 716 {
 717   TH1F* ratio1 = Sigma2VertexGaussian();
 718   TH1F* ratio2 = Sigma2VertexSimulation();
 719
 720   ratio1->SetStats(kFALSE);
 721   ratio2->SetStats(kFALSE);
 722
 723   ratio1->SetMarkerStyle(0);
 724   ratio2->SetMarkerStyle(0);
 725
 726   TLegend* legend = new TLegend(0.647177,0.775424,0.961694,0.966102);
 727   legend->AddEntry(ratio1, "Gaussian");
 728   legend->AddEntry(ratio2, "Simulation");
 729
 730   ratio1->GetXaxis()->SetTitleOffset(1.5);
 731
 732   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
 733   InitPad();
 734
 735   ratio1->Draw();
 736   ratio2->SetLineColor(kRed);
 737   ratio2->Draw("SAME");
 738
 739   legend->Draw();
 740 }
 741
 742 void drawSystematics()
 743 {
 744   //Secondaries();
 745   //DrawDifferentSpecies();
 746   //Composition();
 747
 748   Sigma2VertexSimulation();
 749 }
 750
 751 void DrawdNdEtaDifferences()
 752 {
 753   TH1* hists[5];
 754
 755   TLegend* legend = new TLegend(0.6, 0.73, 0.98, 0.98);
 756   legend->SetFillColor(0);
 757
 758   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
 759   canvas->Divide(2, 1);
 760
 761   canvas->cd(1);
 762
 763   for (Int_t i=0; i<5; ++i)
 764   {
 765     TFile* file = 0;
 766     TString title;
 767
 768     switch(i)
 769     {
 770       case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
 771       case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
 772       case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
 773       case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
 774       case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
 775       default: return;
 776     }
 777
 778     hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
 779     hists[i]->SetTitle("a)");
 780
 781     hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
 782     hists[i]->SetLineColor(i+1);
 783     hists[i]->SetMarkerColor(i+1);
 784     hists[i]->GetXaxis()->SetLabelOffset(0.015);
 785     Prepare1DPlot(hists[i], kFALSE);
 786     hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
 787
 788     legend->AddEntry(hists[i], title);
 789     hists[i]->SetTitle(title);
 790   }
 791   legend->Draw();
 792
 793   canvas->cd(2);
 794   gPad->SetLeftMargin(0.14);
 795
 796   TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
 797   legend2->SetFillColor(0);
 798
 799   for (Int_t i=1; i<5; ++i)
 800   {
 801     legend2->AddEntry(hists[i]);
 802
 803     hists[i]->Divide(hists[0]);
 804     hists[i]->SetTitle("b)");
 805     hists[i]->GetYaxis()->SetRangeUser(0.98, 1.02);
 806     hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
 807     hists[i]->GetYaxis()->SetTitleOffset(1.8);
 808     hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
 809
 810   }
 811
 812   legend2->Draw();
 813
 814   canvas->SaveAs("particlecomposition_result.eps");
 815   canvas->SaveAs("particlecomposition_result.gif");
 816 }