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", 50, 0.05, 5.05);
 664   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
 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", 50, 0.05, 5.05);
 672   Double_t sigma3 = Sigma2VertexCount(tracks, 3);
 673   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
 674     ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(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;N#sigma;% included in 3 #sigma / % included in 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   ratio1->SetLineWidth(2);
 727   ratio2->SetLineWidth(2);
 728
 729   TLegend* legend = new TLegend(0.7, 0.8, 0.95, 0.95);
 730   legend->SetFillColor(0);
 731   legend->AddEntry(ratio1, "Gaussian");
 732   legend->AddEntry(ratio2, "Simulation");
 733
 734   ratio2->SetTitle("");
 735   ratio2->GetYaxis()->SetTitleOffset(1.5);
 736   ratio2->GetXaxis()->SetRangeUser(2, 4);
 737
 738   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
 739   InitPad();
 740
 741   ratio2->SetMarkerStyle(21);
 742   ratio1->SetMarkerStyle(22);
 743
 744   ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
 745   ratio2->SetLineColor(kRed);
 746   ratio2->SetMarkerColor(kRed);
 747   ratio2->Draw("PL");
 748   ratio1->Draw("SAMEPL");
 749
 750   legend->Draw();
 751
 752   canvas->SaveAs("Sigma2VertexCompare.eps");
 753 }
 754
 755 void drawSystematics()
 756 {
 757   //Secondaries();
 758   //DrawDifferentSpecies();
 759   //Composition();
 760
 761   Sigma2VertexSimulation();
 762
 763 }
 764
 765 void DrawdNdEtaDifferences()
 766 {
 767   TH1* hists[5];
 768
 769   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
 770   legend->SetFillColor(0);
 771
 772   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
 773   canvas->Divide(2, 1);
 774
 775   canvas->cd(1);
 776
 777   for (Int_t i=0; i<5; ++i)
 778   {
 779     hists[i] = 0;
 780     TFile* file = 0;
 781     TString title;
 782
 783     switch(i)
 784     {
 785       case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
 786       case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
 787       case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
 788       case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
 789       case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
 790       default: return;
 791     }
 792
 793     if (file)
 794     {
 795       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
 796       hists[i]->SetTitle("a)");
 797
 798       Prepare1DPlot(hists[i], kFALSE);
 799       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
 800       hists[i]->GetYaxis()->SetRangeUser(6, 7);
 801       hists[i]->SetLineColor(i+1);
 802       hists[i]->SetMarkerColor(i+1);
 803       hists[i]->GetXaxis()->SetLabelOffset(0.015);
 804       hists[i]->GetYaxis()->SetTitleOffset(1.5);
 805       gPad->SetLeftMargin(0.12);
 806       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
 807
 808       legend->AddEntry(hists[i], title);
 809       hists[i]->SetTitle(title);
 810     }
 811   }
 812   legend->Draw();
 813
 814   canvas->cd(2);
 815   gPad->SetLeftMargin(0.14);
 816
 817   TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
 818   legend2->SetFillColor(0);
 819
 820   for (Int_t i=1; i<5; ++i)
 821   {
 822     if (hists[i])
 823     {
 824       legend2->AddEntry(hists[i]);
 825
 826       hists[i]->Divide(hists[0]);
 827       hists[i]->SetTitle("b)");
 828       hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
 829       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
 830       hists[i]->GetYaxis()->SetTitleOffset(1.8);
 831       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
 832     }
 833   }
 834
 835   legend2->Draw();
 836
 837   canvas->SaveAs("particlecomposition_result.eps");
 838   canvas->SaveAs("particlecomposition_result.gif");
 839 }
 840
 841 mergeCorrections4SystematicStudies(Char_t* standardCorrectionFileName="correction_map.root",
 842                                              Char_t* systematicCorrectionFileName="systematics.root",
 843                                              Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
 844   //
 845   // Function used to merge standard corrections with vertex
 846   // reconstruction corrections obtained by a certain mix of ND, DD
 847   // and SD events.
 848   //
 849
 850   gSystem->Load("libPWG0base");
 851
 852   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
 853   const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
 854   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5};
 855   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5};
 856
 857   // cross section from Pythia
 858   Float_t sigmaND = 55.2;
 859   Float_t sigmaDD = 9.78;
 860   Float_t sigmaSD = 14.30;
 861
 862   AlidNdEtaCorrection* corrections[21];
 863   for (Int_t j=0; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
 864     AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
 865     correctionStandard->LoadHistograms(standardCorrectionFileName);
 866
 867     // dont take vertexreco from this one
 868     correctionStandard->GetVertexRecoCorrection()->Reset();
 869     // dont take triggerbias from this one
 870     correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
 871
 872     for (Int_t i=0; i<7; i++) {
 873       TString name;
 874       name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
 875       AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
 876
 877       name.Form("vertexRecoND");
 878       AlidNdEtaCorrection* dNdEtaCorrectionVtxND = new AlidNdEtaCorrection(name,name);
 879       dNdEtaCorrectionVtxND->LoadHistograms(systematicCorrectionFileName, name);
 880       name.Form("vertexRecoDD");
 881       AlidNdEtaCorrection* dNdEtaCorrectionVtxDD = new AlidNdEtaCorrection(name,name);
 882       dNdEtaCorrectionVtxDD->LoadHistograms(systematicCorrectionFileName, name);
 883       name.Form("vertexRecoSD");
 884       AlidNdEtaCorrection* dNdEtaCorrectionVtxSD = new AlidNdEtaCorrection(name,name);
 885       dNdEtaCorrectionVtxSD->LoadHistograms(systematicCorrectionFileName, name);
 886
 887       name.Form("triggerBiasND");
 888       AlidNdEtaCorrection* dNdEtaCorrectionTriggerND = new AlidNdEtaCorrection(name,name);
 889       dNdEtaCorrectionTriggerND->LoadHistograms(systematicCorrectionFileName, name);
 890       name.Form("triggerBiasDD");
 891       AlidNdEtaCorrection* dNdEtaCorrectionTriggerDD = new AlidNdEtaCorrection(name,name);
 892       dNdEtaCorrectionTriggerDD->LoadHistograms(systematicCorrectionFileName, name);
 893       name.Form("triggerBiasSD");
 894       AlidNdEtaCorrection* dNdEtaCorrectionTriggerSD = new AlidNdEtaCorrection(name,name);
 895       dNdEtaCorrectionTriggerSD->LoadHistograms(systematicCorrectionFileName, name);
 896
 897       // calculating relative
 898       Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 899       Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 900       Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 901
 902       printf(Form("%s : ND=%.1f\%, DD=%.1f\%, SD=%.1f\% \n",changes[i],nd,dd,sd));
 903       current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",nd,dd,sd));
 904       current->SetTitle(name);
 905
 906       // scale
 907       if (j == 0 || j == 2)
 908       {
 909         dNdEtaCorrectionVtxDD->GetVertexRecoCorrection()->GetMeasuredHistogram()->Scale(scalesDD[i]);
 910         dNdEtaCorrectionVtxSD->GetVertexRecoCorrection()->GetMeasuredHistogram()->Scale(scalesSD[i]);
 911         dNdEtaCorrectionVtxDD->GetVertexRecoCorrection()->GetGeneratedHistogram()->Scale(scalesDD[i]);
 912         dNdEtaCorrectionVtxSD->GetVertexRecoCorrection()->GetGeneratedHistogram()->Scale(scalesSD[i]);
 913       }
 914       if (j == 1 || j == 2)
 915       {
 916         dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionINEL()->GetMeasuredHistogram()->Scale(scalesDD[i]);
 917         dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionINEL()->GetMeasuredHistogram()->Scale(scalesSD[i]);
 918         dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionINEL()->GetGeneratedHistogram()->Scale(scalesDD[i]);
 919         dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionINEL()->GetGeneratedHistogram()->Scale(scalesSD[i]);
 920       }
 921
 922       //clear track, trigger in Vtx correction
 923       dNdEtaCorrectionVtxND->GetTrack2ParticleCorrection()->Reset();
 924       dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionNSD()->Reset();
 925       dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionND()->Reset();
 926       dNdEtaCorrectionVtxND->GetTriggerBiasCorrectionINEL()->Reset();
 927       dNdEtaCorrectionVtxSD->GetTrack2ParticleCorrection()->Reset();
 928       dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionNSD()->Reset();
 929       dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionND()->Reset();
 930       dNdEtaCorrectionVtxSD->GetTriggerBiasCorrectionINEL()->Reset();
 931       dNdEtaCorrectionVtxDD->GetTrack2ParticleCorrection()->Reset();
 932       dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionNSD()->Reset();
 933       dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionND()->Reset();
 934       dNdEtaCorrectionVtxDD->GetTriggerBiasCorrectionINEL()->Reset();
 935
 936       //clear track, vertexreco in trigger correction
 937       dNdEtaCorrectionTriggerND->GetTrack2ParticleCorrection()->Reset();
 938       dNdEtaCorrectionTriggerND->GetTriggerBiasCorrectionNSD()->Reset();
 939       dNdEtaCorrectionTriggerND->GetTriggerBiasCorrectionND()->Reset();
 940       dNdEtaCorrectionTriggerND->GetVertexRecoCorrection()->Reset();
 941       dNdEtaCorrectionTriggerSD->GetTrack2ParticleCorrection()->Reset();
 942       dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionNSD()->Reset();
 943       dNdEtaCorrectionTriggerSD->GetTriggerBiasCorrectionND()->Reset();
 944       dNdEtaCorrectionTriggerSD->GetVertexRecoCorrection()->Reset();
 945       dNdEtaCorrectionTriggerDD->GetTrack2ParticleCorrection()->Reset();
 946       dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionNSD()->Reset();
 947       dNdEtaCorrectionTriggerDD->GetTriggerBiasCorrectionND()->Reset();
 948       dNdEtaCorrectionTriggerDD->GetVertexRecoCorrection()->Reset();
 949
 950       TList collection;
 951       collection.Add(correctionStandard);
 952       collection.Add(dNdEtaCorrectionVtxND);
 953       collection.Add(dNdEtaCorrectionVtxDD);
 954       collection.Add(dNdEtaCorrectionVtxSD);
 955       collection.Add(dNdEtaCorrectionTriggerND);
 956       collection.Add(dNdEtaCorrectionTriggerDD);
 957       collection.Add(dNdEtaCorrectionTriggerSD);
 958
 959       current->Merge(&collection);
 960       current->Finish();
 961
 962       corrections[i+j*7] = current;
 963     }
 964   }
 965
 966   TFile* fout = new TFile(outputFileName,"RECREATE");
 967
 968   for (Int_t i=0; i<21; i++)
 969     corrections[i]->SaveHistograms();
 970
 971   fout->Write();
 972   fout->Close();
 973 }
 974
 975
 976 DrawVertexRecoSyst(const char* plot = "vtxreco") {
 977
 978   Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless"};
 979   Char_t* descr[]  =   {"",
 980                         "#sigma_{DD}' = 1.5#times#sigma_{DD}",
 981                         "#sigma_{DD}' = 0.5#times#sigma_{DD}",
 982                         "#sigma_{SD}' = 1.5#times#sigma_{SD}",
 983                         "#sigma_{SD}' = 0.5#times#sigma_{SD}",
 984                         "#sigma_{D}' = 1.5#times#sigma_{D}",
 985                         "#sigma_{D}' = 0.5#times#sigma_{D}"};
 986
 987   Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.5, 0.5, 1.5, 0.5};
 988   Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.5, 0.5};
 989
 990   Int_t markers[] = {22,22,22,22,22,22,22};
 991   Int_t colors[]  = {1,2,3,4,6,8,102};
 992
 993   // cross section from Pythia
 994   Float_t sigmaND = 55.2;
 995   Float_t sigmaDD = 9.78;
 996   Float_t sigmaSD = 14.30;
 997
 998
 999   TH1F* dNdEta[7];
1000
1001   TH1F* ratios[7];
1002
1003   TFile* fin;
1004
1005
1006   for (Int_t i=0; i<7; i++) {
1007     // calculating relative
1008     fin = TFile::Open(Form("systematics_%s_%s.root",plot,changes[i]));
1009
1010     dNdEta[i] = (TH1F*)(fin->Get("dndeta/dndeta_dNdEta_corrected_2"))->Clone();
1011
1012     for (Int_t b=0; b<dNdEta[i]->GetNbinsX(); b++) {
1013       if (TMath::Abs(dNdEta[i]->GetBinCenter(b))>0.9) {
1014         dNdEta[i]->SetBinContent(b,0);
1015         dNdEta[i]->SetBinError(b,0);
1016       }
1017     }
1018
1019     dNdEta[i]->Rebin(4);
1020
1021     dNdEta[i]->SetLineWidth(2);
1022     dNdEta[i]->SetLineColor(colors[i]);
1023     dNdEta[i]->SetMarkerStyle(markers[i]);
1024     dNdEta[i]->SetMarkerSize(0.8);
1025
1026
1027     ratios[i] = (TH1F*)dNdEta[i]->Clone("ratio");
1028     ratios[i]->Divide(ratios[i],dNdEta[0],1,1,"B");
1029
1030     ratios[i]->SetName(changes[i]);
1031     ratios[i]->SetTitle(changes[i]);
1032   }
1033
1034   //##########################################################
1035
1036   gStyle->SetOptStat(0);
1037   gStyle->SetOptTitle(0);
1038   gStyle->SetOptFit(0);
1039
1040   gStyle->SetTextSize(0.2);
1041   gStyle->SetTitleSize(0.05,"xyz");
1042   //gStyle->SetTitleFont(133, "xyz");
1043   //gStyle->SetLabelFont(133, "xyz");
1044   //gStyle->SetLabelSize(17, "xyz");
1045   gStyle->SetLabelOffset(0.01, "xyz");
1046
1047
1048   gStyle->SetTitleOffset(1.2, "y");
1049   gStyle->SetTitleOffset(1.2, "x");
1050   gStyle->SetEndErrorSize(0.0);
1051
1052   //##############################################
1053
1054   //making canvas and pads
1055   TCanvas *c = new TCanvas(Form("vertex_reco_syst_%s", plot), Form("vertex_reco_syst_%s", plot),600,500);
1056
1057   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1058
1059   p1->SetBottomMargin(0.15);
1060   p1->SetTopMargin(0.03);
1061   p1->SetLeftMargin(0.15);
1062   p1->SetRightMargin(0.03);
1063
1064   p1->SetGridx();
1065   p1->SetGridy();
1066
1067   p1->Draw();
1068   p1->cd();
1069
1070
1071
1072
1073   TH2F* null = new TH2F("","",100,-1.5,1.5,100,0.9601,1.099);
1074   null->SetXTitle("#eta");
1075   null->GetXaxis()->CenterTitle(kTRUE);
1076   null->SetYTitle("dN/d#eta / dN/d#eta_{pythia}");
1077   null->GetYaxis()->CenterTitle(kTRUE);
1078   null->Draw();
1079
1080   for (Int_t i=1; i<7; i++)
1081     ratios[i]->Draw("same");
1082   //     dNdEta[i]->Draw("same");
1083
1084   TLatex* text[7];
1085   for (Int_t i=1; i<7; i++) {
1086
1087     text[i] = new TLatex(0.75,0.95-0.05*i, descr[i]);
1088     text[i]->SetTextSize(0.04);
1089     text[i]->SetTextColor(colors[i]);
1090     text[i]->SetNDC(kTRUE);
1091     text[i]->Draw();
1092
1093
1094   }
1095   //text(0.2,0.88,"Effect of changing",0.045,1,kTRUE);
1096   //text(0.2,0.83,"relative cross-sections",0.045,1,kTRUE);
1097   //text(0.2,0.78,"(vertex reconstruction corr.)",0.043,13,kTRUE);
1098
1099   c->SaveAs(Form("%s.gif", c->GetName()));
1100
1101 }
1102
1103
1104
1105 DrawTriggerEfficiency(Char_t* fileName) {
1106
1107   gStyle->SetOptStat(0);
1108   gStyle->SetOptTitle(0);
1109   gStyle->SetOptFit(0);
1110
1111   gStyle->SetTextSize(0.04);
1112   gStyle->SetTitleSize(0.05,"xyz");
1113   //gStyle->SetTitleFont(133, "xyz");
1114   //gStyle->SetLabelFont(133, "xyz");
1115   //gStyle->SetLabelSize(17, "xyz");
1116   gStyle->SetLabelOffset(0.01, "xyz");
1117
1118   gStyle->SetTitleOffset(1.1, "y");
1119   gStyle->SetTitleOffset(1.1, "x");
1120   gStyle->SetEndErrorSize(0.0);
1121
1122   //##############################################
1123
1124   //making canvas and pads
1125   TCanvas *c = new TCanvas("trigger_eff", "",600,500);
1126
1127   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1128
1129   p1->SetBottomMargin(0.15);
1130   p1->SetTopMargin(0.03);
1131   p1->SetLeftMargin(0.15);
1132   p1->SetRightMargin(0.03);
1133
1134   p1->SetGridx();
1135   p1->SetGridy();
1136
1137   p1->Draw();
1138   p1->cd();
1139
1140   gSystem->Load("libPWG0base");
1141
1142   AlidNdEtaCorrection* corrections[4];
1143   AliCorrectionMatrix2D* triggerBiasCorrection[4];
1144
1145   TH1F* hTriggerEffInv[4];
1146   TH1F* hTriggerEff[4];
1147
1148   Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
1149
1150   for (Int_t i=0; i<4; i++) {
1151     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1152     corrections[i]->LoadHistograms(fileName, names[i]);
1153
1154     triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1155
1156
1157     hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1158     hTriggerEff[i]    = (TH1F*)hTriggerEffInv[i]->Clone();
1159
1160     for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1161       hTriggerEff[i]->SetBinContent(b,1);
1162       hTriggerEff[i]->SetBinError(b,0);
1163     }
1164     hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1165     hTriggerEff[i]->Scale(100);
1166   }
1167
1168   Int_t colors[] = {2,3,4,1};
1169   Float_t effs[4];
1170   for (Int_t i=0; i<4; i++) {
1171     hTriggerEff[i]->Fit("pol0","","",-20,20);
1172     effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1173     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1174     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1175     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1176     cout << effs[i] << endl;
1177   }
1178
1179
1180   Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1181   TLatex* latex[4];
1182
1183   TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1184   null->SetXTitle("Vertex z [cm]");
1185   null->GetXaxis()->CenterTitle(kTRUE);
1186   null->SetYTitle("Trigger efficiency [%]");
1187   null->GetYaxis()->CenterTitle(kTRUE);
1188   null->Draw();
1189
1190
1191   for (Int_t i=0; i<4; i++) {
1192     hTriggerEff[i]->SetLineWidth(2);
1193     hTriggerEff[i]->SetLineColor(colors[i]);
1194
1195     hTriggerEff[i]->Draw("same");
1196
1197     latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1198     latex[i]->SetTextColor(colors[i]);
1199     latex[i]->Draw();
1200   }
1201
1202 }
1203
1204
1205 DrawSpectraPID(Char_t* fileName) {
1206
1207   gSystem->Load("libPWG0base");
1208
1209   Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1210   AlidNdEtaCorrection* corrections[4];
1211   AliCorrectionMatrix3D* trackToPartCorrection[4];
1212
1213   TH1F* measuredPt[4];
1214   TH1F* generatedPt[4];
1215   TH1F* ratioPt[4];
1216
1217   for (Int_t i=0; i<4; i++) {
1218     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1219     corrections[i]->LoadHistograms(fileName, names[i]);
1220
1221     trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1222
1223     Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1224     Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1225     Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1226     Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1227
1228     measuredPt[i]  = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1229     generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1230     ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1231     ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1232   }
1233
1234   ratioPt[0]->Draw();
1235
1236   for (Int_t i=0; i<3; i++) {
1237     ratioPt[i]->SetLineColor(i+1);
1238     ratioPt[i]->SetLineWidth(2);
1239
1240     ratioPt[i]->Draw("same");
1241
1242   }
1243
1244   return;
1245   measuredPt[0]->SetLineColor(2);
1246   measuredPt[0]->SetLineWidth(5);
1247
1248   measuredPt[0]->Draw();
1249   generatedPt[0]->Draw("same");
1250 }
1251
1252 void changePtSpectrum()
1253 {
1254   TFile* file = TFile::Open("pt.root");
1255   TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt"));
1256
1257   hist->Rebin(3);
1258   hist->Scale(1.0/3);
1259
1260   TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1261   TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1262
1263   TH1F* scale1 =  dynamic_cast<TH1F*> (hist->Clone("scale1"));
1264   TH1F* scale2 =  dynamic_cast<TH1F*> (hist->Clone("scale2"));
1265
1266   Float_t ptCutOff = 0.3;
1267
1268   for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1269   {
1270     if (hist->GetBinCenter(i) > ptCutOff)
1271     {
1272       scale1->SetBinContent(i, 1);
1273       scale2->SetBinContent(i, 1);
1274     }
1275     else
1276     {
1277       // 90 % at pt = 0, 0% at pt = ptcutoff
1278       scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
1279
1280       // 110% at pt = 0, ...
1281       scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * 0.3);
1282     }
1283   }
1284
1285   new TCanvas;
1286   scale1->Draw();
1287   scale2->SetLineColor(kRed);
1288   scale2->Draw("SAME");
1289
1290   clone1->Multiply(scale1);
1291   clone2->Multiply(scale2);
1292
1293   Prepare1DPlot(hist);
1294   Prepare1DPlot(clone1);
1295   Prepare1DPlot(clone2);
1296
1297   hist->SetTitle(";p_{T} in GeV/c;dN/dp_{T} in c/GeV");
1298   hist->GetXaxis()->SetRangeUser(0, 0.7);
1299   hist->GetYaxis()->SetRangeUser(0, clone2->GetMaximum() * 1.1);
1300   hist->GetYaxis()->SetTitleOffset(1);
1301
1302   TCanvas* canvas = new TCanvas;
1303   hist->Draw();
1304   clone1->SetLineColor(kRed);
1305   clone1->Draw("SAME");
1306   clone2->SetLineColor(kBlue);
1307   clone2->Draw("SAME");
1308
1309   Float_t fraction =  hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1310   Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1311   Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1312
1313   printf("%f %f %f\n", fraction, fraction1, fraction2);
1314   printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1315
1316   canvas->SaveAs("changePtSpectrum.gif");
1317 }