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 Int_t markers[] = {20,20,21,22,23,28,29};
  27 Int_t colors[]  = {1,2,3,4,6,8,102};
  28
  29 void loadlibs()
  30 {
  31   gSystem->Load("libTree");
  32   gSystem->Load("libVMC");
  33
  34   gSystem->Load("libSTEERBase");
  35   gSystem->Load("libANALYSIS");
  36   gSystem->Load("libPWG0base");
  37 }
  38
  39 void SetRanges(TAxis* axis)
  40 {
  41   if (strcmp(axis->GetTitle(), "#eta") == 0)
  42     axis->SetRangeUser(-1.7999, 1.7999);
  43   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
  44     axis->SetRangeUser(0, 9.9999);
  45   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
  46     axis->SetRangeUser(-15, 14.9999);
  47   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
  48     axis->SetRangeUser(0, 99.9999);
  49 }
  50
  51 void SetRanges(TH1* hist)
  52 {
  53   SetRanges(hist->GetXaxis());
  54   SetRanges(hist->GetYaxis());
  55   SetRanges(hist->GetZaxis());
  56 }
  57
  58 void Prepare3DPlot(TH3* hist)
  59 {
  60   hist->GetXaxis()->SetTitleOffset(1.5);
  61   hist->GetYaxis()->SetTitleOffset(1.5);
  62   hist->GetZaxis()->SetTitleOffset(1.5);
  63
  64   hist->SetStats(kFALSE);
  65 }
  66
  67 void Prepare2DPlot(TH2* hist)
  68 {
  69   hist->SetStats(kFALSE);
  70   hist->GetYaxis()->SetTitleOffset(1.4);
  71
  72   SetRanges(hist);
  73 }
  74
  75 void Prepare1DPlot(TH1* hist, Bool_t setRanges = kTRUE)
  76 {
  77   hist->SetLineWidth(2);
  78   hist->SetStats(kFALSE);
  79
  80   hist->GetXaxis()->SetTitleOffset(1.2);
  81   hist->GetYaxis()->SetTitleOffset(1.2);
  82
  83   if (setRanges)
  84     SetRanges(hist);
  85 }
  86
  87 void InitPad()
  88 {
  89   if (!gPad)
  90     return;
  91
  92   gPad->Range(0, 0, 1, 1);
  93   gPad->SetLeftMargin(0.15);
  94   //gPad->SetRightMargin(0.05);
  95   //gPad->SetTopMargin(0.13);
  96   //gPad->SetBottomMargin(0.1);
  97
  98   gPad->SetGridx();
  99   gPad->SetGridy();
 100 }
 101
 102 void InitPadCOLZ()
 103 {
 104   gPad->Range(0, 0, 1, 1);
 105   gPad->SetRightMargin(0.15);
 106   gPad->SetLeftMargin(0.12);
 107
 108   gPad->SetGridx();
 109   gPad->SetGridy();
 110 }
 111
 112 void Secondaries()
 113 {
 114   TFile* file = TFile::Open("systematics.root");
 115
 116   TH2F* secondaries = dynamic_cast<TH2F*> (file->Get("fSecondaries"));
 117   if (!secondaries)
 118   {
 119     printf("Could not read histogram\n");
 120     return;
 121   }
 122
 123   TCanvas* canvas = new TCanvas("Secondaries", "Secondaries", 1000, 1000);
 124   canvas->Divide(3, 3);
 125   for (Int_t i=1; i<=8; i++)
 126   {
 127     TH1D* hist = secondaries->ProjectionY(Form("proj_%d", i), i, i);
 128     hist->SetTitle(secondaries->GetXaxis()->GetBinLabel(i));
 129
 130     canvas->cd(i);
 131     hist->Draw();
 132   }
 133 }
 134
 135 void Track2Particle1DComposition(const char** fileNames, Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
 136 {
 137   gSystem->Load("libPWG0base");
 138
 139   TString canvasName;
 140   canvasName.Form("Track2Particle1DComposition");
 141   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
 142   canvas->Divide(3, 1);
 143
 144   TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
 145
 146   for (Int_t i=0; i<folderCount; ++i)
 147   {
 148     Correction1DCreatePlots(fileNames[i], folderNames[i], upperPtLimit);
 149
 150     TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("generated_x_div_measured_x"));
 151     TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject("generated_y_div_measured_y"));
 152     TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("generated_z_div_measured_z"));
 153
 154     Prepare1DPlot(corrX);
 155     Prepare1DPlot(corrY);
 156     Prepare1DPlot(corrZ);
 157
 158     const char* title = "Track2Particle Correction";
 159     corrX->SetTitle(title);
 160     corrY->SetTitle(title);
 161     corrZ->SetTitle(title);
 162
 163     corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
 164
 165     corrX->SetLineColor(i+1);
 166     corrY->SetLineColor(i+1);
 167     corrZ->SetLineColor(i+1);
 168
 169     canvas->cd(1);
 170     InitPad();
 171     corrX->DrawCopy(((i>0) ? "SAME" : ""));
 172
 173     canvas->cd(2);
 174     InitPad();
 175     corrY->DrawCopy(((i>0) ? "SAME" : ""));
 176
 177     canvas->cd(3);
 178     InitPad();
 179     corrZ->DrawCopy(((i>0) ? "SAME" : ""));
 180
 181     legend->AddEntry(corrZ, folderNames[i]);
 182   }
 183
 184   legend->Draw();
 185
 186   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
 187   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
 188 }
 189
 190 TH1** DrawRatios(const char* fileName = "systematics.root")
 191 {
 192   gSystem->Load("libPWG0base");
 193
 194   TFile* file = TFile::Open(fileName);
 195
 196   TString canvasName;
 197   canvasName.Form("DrawRatios");
 198   TCanvas* canvas = new TCanvas(canvasName, canvasName, 800, 400);
 199   canvas->Divide(2, 1);
 200   canvas->cd(1);
 201
 202   TH1** ptDists = new TH1*[4];
 203
 204   TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
 205
 206   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
 207   const char* particleNames[] = { "#pi", "K", "p", "other" };
 208   for (Int_t i=0; i<4; ++i)
 209   {
 210     AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(folderNames[i], folderNames[i]);
 211     dNdEtaCorrection->LoadHistograms(fileName, folderNames[i]);
 212
 213     TH3F* gene = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
 214
 215     gene->GetYaxis()->SetRangeUser(-0.8, 0.8);
 216     gene->GetXaxis()->SetRangeUser(-10, 10);
 217
 218     ptDists[i] = dynamic_cast<TH1*> (gene->Project3D("z")->Clone(Form("%s_z", folderNames[i])));
 219     ptDists[i]->SetTitle(";p_{T};Count");
 220     if (!ptDists[i])
 221     {
 222       printf("Problem getting distribution %d.\n", i);
 223       return 0;
 224     }
 225
 226     ptDists[i]->GetYaxis()->SetRangeUser(1, ptDists[i]->GetMaximum()*1.1);
 227     ptDists[i]->GetXaxis()->SetRangeUser(0, 9.9);
 228     ptDists[i]->SetLineColor(i+1);
 229     ptDists[i]->DrawCopy((i == 0) ? "" : "SAME");
 230     ptDists[i]->GetYaxis()->SetRange(0, 0);
 231
 232     legend->AddEntry(ptDists[i], particleNames[i]);
 233   }
 234   gPad->SetLogy();
 235
 236   TH1* total = dynamic_cast<TH1*> (ptDists[0]->Clone("total"));
 237
 238   for (Int_t i=1; i<4; ++i)
 239     total->Add(ptDists[i]);
 240
 241   canvas->cd(2);
 242   for (Int_t i=0; i<4; ++i)
 243   {
 244     ptDists[i]->Divide(total);
 245     ptDists[i]->SetStats(kFALSE);
 246     ptDists[i]->SetTitle(";p_{T};Fraction of total");
 247     ptDists[i]->GetYaxis()->SetRangeUser(0, 1);
 248     ptDists[i]->Draw((i == 0) ? "" : "SAME");
 249   }
 250   legend->SetFillColor(0);
 251   legend->Draw();
 252
 253   canvas->SaveAs("DrawRatios.gif");
 254
 255
 256   canvas = new TCanvas("PythiaRatios", "PythiaRatios", 500, 500);
 257   for (Int_t i=0; i<4; ++i)
 258   {
 259     TH1* hist = ptDists[i]->Clone();
 260     hist->GetXaxis()->SetRangeUser(0, 1.9);
 261     hist->Draw((i == 0) ? "" : "SAME");
 262   }
 263   legend->Draw();
 264
 265   canvas->SaveAs("pythiaratios.eps");
 266
 267   file->Close();
 268
 269   return ptDists;
 270 }
 271
 272 void DrawCompareToReal()
 273 {
 274   gROOT->ProcessLine(".L drawPlots.C");
 275
 276   const char* fileNames[] = { "correction_map.root", "new_compositions.root" };
 277   const char* folderNames[] = { "dndeta_correction", "PythiaRatios" };
 278
 279   Track2Particle1DComposition(fileNames, 2, folderNames);
 280 }
 281
 282 void DrawDifferentSpecies()
 283 {
 284   gROOT->ProcessLine(".L drawPlots.C");
 285
 286   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "systematics.root" };
 287   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
 288
 289   Track2Particle1DComposition(fileNames, 4, folderNames);
 290 }
 291
 292 void DrawpiKpAndCombined()
 293 {
 294   gROOT->ProcessLine(".L drawPlots.C");
 295
 296   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
 297   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
 298
 299   Track2Particle1DComposition(fileNames, 4, folderNames);
 300 }
 301
 302 void DrawSpeciesAndCombination()
 303 {
 304   gROOT->ProcessLine(".L drawPlots.C");
 305
 306   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "new_compositions.root" };
 307   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "PythiaRatios" };
 308
 309   Track2Particle1DComposition(fileNames, 4, folderNames);
 310 }
 311
 312 void DrawSimulatedVsCombined()
 313 {
 314   gROOT->ProcessLine(".L drawPlots.C");
 315
 316   const char* fileNames[] = { "new_compositions.root", "new_compositions.root" };
 317   const char* folderNames[] = { "Pythia", "PythiaRatios" };
 318
 319   Track2Particle1DComposition(fileNames, 2, folderNames);
 320 }
 321
 322 void DrawBoosts()
 323 {
 324   gROOT->ProcessLine(".L drawPlots.C");
 325
 326   const char* fileNames[] = { "new_compositions.root", "new_compositions.root", "new_compositions.root", "new_compositions.root" };
 327   const char* folderNames[] = { "PythiaRatios", "PiBoosted", "KBoosted", "pBoosted" };
 328
 329   Track2Particle1DComposition(fileNames, 4, folderNames);
 330 }
 331
 332 TH2F* HistogramDifferences(const char* filename, const char* folder1, const char* folder2)
 333 {
 334   gSystem->Load("libPWG0base");
 335
 336   AlidNdEtaCorrection* fdNdEtaCorrection[2];
 337
 338   TFile::Open(filename);
 339
 340   fdNdEtaCorrection[0] = new AlidNdEtaCorrection(folder1, folder1);
 341   fdNdEtaCorrection[0]->LoadHistograms(filename, folder1);
 342
 343   fdNdEtaCorrection[1] = new AlidNdEtaCorrection(folder2, folder2);
 344   fdNdEtaCorrection[1]->LoadHistograms(filename, folder2);
 345
 346   TH3F* hist1 = fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
 347   TH3F* hist2 = fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
 348
 349   //TH1F* difference = new TH1F("difference", Form(";#DeltaC_{pT, z, #eta} %s / %s;Count", folder2, folder1), 1000, 0.9, 1.1);
 350   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());
 351
 352   for (Int_t x=hist1->GetXaxis()->FindBin(-10); x<=hist1->GetXaxis()->FindBin(10); ++x)
 353     for (Int_t y=hist1->GetYaxis()->FindBin(-0.8); y<=hist1->GetYaxis()->FindBin(0.8); ++y)
 354       for (Int_t z=hist1->GetZaxis()->FindBin(0.3); z<=hist1->GetZaxis()->FindBin(9.9); ++z)
 355         if (hist1->GetBinContent(x, y, z) != 0)
 356           difference->Fill(hist2->GetBinContent(x, y, z) / hist1->GetBinContent(x, y, z), hist1->GetYaxis()->GetBinCenter(y));
 357
 358   difference->GetYaxis()->SetRangeUser(-0.8, 0.8);
 359
 360   printf("Over-/Underflow bins: %d %d\n", difference->GetBinContent(0), difference->GetBinContent(difference->GetNbinsX()+1));
 361
 362   return difference;
 363 }
 364
 365 void HistogramDifferences()
 366 {
 367   TH2F* KBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "KBoosted");
 368   TH2F* pBoosted = HistogramDifferences("new_compositions.root", "PythiaRatios", "pBoosted");
 369   TH2F* KReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "KReduced");
 370   TH2F* pReduced = HistogramDifferences("new_compositions.root", "PythiaRatios", "pReduced");
 371
 372   TCanvas* canvas = new TCanvas("HistogramDifferences", "HistogramDifferences", 1000, 1000);
 373   canvas->Divide(2, 2);
 374
 375   canvas->cd(1);
 376   KBoosted->GetXaxis()->SetRangeUser(-0.05, 0.05);
 377   KBoosted->Draw("COLZ");
 378
 379   canvas->cd(2);
 380   KReduced->GetXaxis()->SetRangeUser(-0.05, 0.05);
 381   KReduced->Draw("COLZ");
 382
 383   canvas->cd(3);
 384   pBoosted->GetXaxis()->SetRangeUser(-0.02, 0.02);
 385   pBoosted->Draw("COLZ");
 386
 387   canvas->cd(4);
 388   pReduced->GetXaxis()->SetRangeUser(-0.02, 0.02);
 389   pReduced->Draw("COLZ");
 390
 391   canvas->SaveAs("HistogramDifferences.gif");
 392
 393   canvas = new TCanvas("HistogramDifferencesProfile", "HistogramDifferencesProfile", 1000, 500);
 394   canvas->Divide(2, 1);
 395
 396   canvas->cd(1);
 397   gPad->SetBottomMargin(0.13);
 398   KBoosted->SetTitle("a) Ratio of correction maps");
 399   KBoosted->SetStats(kFALSE);
 400   KBoosted->GetXaxis()->SetTitleOffset(1.4);
 401   KBoosted->GetXaxis()->SetLabelOffset(0.02);
 402   KBoosted->Draw("COLZ");
 403
 404   canvas->cd(2);
 405   gPad->SetGridx();
 406   gPad->SetGridy();
 407
 408   TLegend* legend = new TLegend(0.73, 0.73, 0.98, 0.98);
 409
 410   for (Int_t i=0; i<4; ++i)
 411   {
 412     TH2F* hist = 0;
 413     TString name;
 414     switch (i)
 415     {
 416       case 0: hist = KBoosted; name = "K enhanced"; break;
 417       case 1: hist = KReduced; name = "K reduced"; break;
 418       case 2: hist = pBoosted; name = "p enhanced"; break;
 419       case 3: hist = pReduced; name = "p reduced"; break;
 420     }
 421
 422     TProfile* profile = hist->ProfileY();
 423     profile->SetTitle("b) Mean and RMS");
 424     profile->GetXaxis()->SetRange(hist->GetYaxis()->GetFirst(), hist->GetYaxis()->GetLast());
 425     profile->GetXaxis()->SetTitleOffset(1.2);
 426     profile->GetXaxis()->SetLabelOffset(0.02);
 427     profile->GetYaxis()->SetRangeUser(0.98, 1.02);
 428     profile->SetStats(kFALSE);
 429     profile->SetLineColor(i+1);
 430     profile->SetMarkerColor(i+1);
 431     profile->DrawCopy(((i > 0) ? "SAME" : ""));
 432
 433
 434     legend->AddEntry(profile, name);
 435   }
 436
 437   legend->Draw();
 438   canvas->SaveAs("particlecomposition_result.eps");
 439 }
 440
 441
 442 void ScalePtDependent(TH3F* hist, TF1* function)
 443 {
 444   // assumes that pt is the third dimension of hist
 445   // scales with function(pt)
 446
 447   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
 448   {
 449     Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z));
 450     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
 451
 452     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
 453       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
 454         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
 455   }
 456 }
 457
 458 void ScalePtDependent(TH3F* hist, TH1* function)
 459 {
 460   // assumes that pt is the third dimension of hist
 461   // scales with histogram(pt)
 462
 463   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
 464   {
 465     Double_t factor = function->GetBinContent(function->GetXaxis()->FindBin(hist->GetZaxis()->GetBinCenter(z)));
 466     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
 467
 468     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
 469       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
 470         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
 471   }
 472 }
 473
 474 const char* ChangeComposition(void** correctionPointer, Int_t index)
 475 {
 476   AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
 477
 478   switch (index)
 479   {
 480     case 0: // result from pp events
 481       {
 482         TFile::Open("pythiaratios.root");
 483
 484         for (Int_t i=0; i<4; ++i)
 485         {
 486           TString name;
 487           name.Form("correction_%d", i);
 488           fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
 489           fdNdEtaCorrection[i]->LoadHistograms("pythiaratios.root", name);
 490         }
 491       }
 492       return "Pythia";
 493       break;
 494
 495     case 1: // each species rated with pythia ratios
 496       /*TH1** ptDists = DrawRatios("pythiaratios.root");
 497
 498       for (Int_t i=0; i<3; ++i)
 499       {
 500         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDists[i]);
 501         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDists[i]);
 502       }*/
 503       return "PythiaRatios";
 504       break;
 505
 506       // one species enhanced / reduced
 507     case 2: // + 50% kaons
 508     case 3: // - 50% kaons
 509     case 4: // + 50% protons
 510     case 5: // - 50% protons
 511     case 6: // + 50% kaons + 50% protons
 512     case 7: // - 50% kaons - 50% protons
 513     case 8: // + 50% kaons - 50% protons
 514     case 9: // - 50% kaons + 50% protons
 515       TString* str = new TString;
 516       if (index < 6)
 517       {
 518         Int_t correctionIndex = index / 2;
 519         Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
 520
 521         fdNdEtaCorrection[correctionIndex]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
 522         str->Form("%s%s", (correctionIndex == 0) ? "Pi" : ((correctionIndex == 1) ? "K" : (correctionIndex == 2) ? "p" : "others"), (index % 2 == 0) ? "Boosted" : "Reduced");
 523       }
 524       else
 525       {
 526         Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
 527         fdNdEtaCorrection[1]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
 528         str->Form("%s%s", "K", (scaleFactor > 1) ? "Boosted" : "Reduced");
 529
 530         if (index >= 8)
 531           scaleFactor = (index % 2 == 0) ? 0.5 : 1.5;
 532         fdNdEtaCorrection[2]->GetTrack2ParticleCorrection()->GetTrackCorrection()->Scale(scaleFactor);
 533         *str += Form("%s%s", "p", (scaleFactor > 1) ? "Boosted" : "Reduced");
 534       }
 535
 536       return str->Data();
 537       break;
 538
 539       // each species rated with pythia ratios
 540     case 12: // + 50% pions
 541     case 13: // - 50% pions
 542     case 14: // + 50% kaons
 543     case 15: // - 50% kaons
 544     case 16: // + 50% protons
 545     case 17: // - 50% protons
 546       TH1** ptDists = DrawRatios("pythiaratios.root");
 547       Int_t functionIndex = (index - 2) / 2;
 548       Double_t scaleFactor = (index % 2 == 0) ? 1.5 : 0.5;
 549       ptDists[functionIndex]->Scale(scaleFactor);
 550
 551       for (Int_t i=0; i<3; ++i)
 552       {
 553         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDists[i]);
 554         ScalePtDependent(fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDists[i]);
 555       }
 556       TString* str = new TString;
 557       str->Form("%s%s", (functionIndex == 0) ? "Pi" : ((functionIndex == 1) ? "K" : "p"), (index % 2 == 0) ? "Boosted" : "Reduced");
 558       return str->Data();
 559       break;
 560
 561     case 999:
 562       TF1* ptDependence = new TF1("simple", "x", 0, 100);
 563       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram(), ptDependence);
 564       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram(), ptDependence);
 565       break;
 566
 567   }
 568
 569   return "noname";
 570 }
 571
 572 void Composition()
 573 {
 574   loadlibs();
 575
 576   gSystem->Unlink("new_compositions.root");
 577   gSystem->Unlink("new_compositions_analysis.root");
 578
 579   const char* names[] = { "pi", "K", "p", "other" };
 580   TH1* hRatios[20];
 581
 582   Int_t nCompositions = 10;
 583   Int_t counter = 0;
 584   for (Int_t comp = 1; comp < nCompositions; ++comp)
 585   {
 586     AlidNdEtaCorrection* fdNdEtaCorrection[4];
 587
 588     TFile::Open("correction_mapparticle-species.root");
 589
 590     for (Int_t i=0; i<4; ++i)
 591     {
 592       TString name;
 593       name.Form("dndeta_correction_%s", names[i]);
 594       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
 595       fdNdEtaCorrection[i]->LoadHistograms();
 596     }
 597
 598     const char* newName = ChangeComposition(fdNdEtaCorrection, comp);
 599
 600     Double_t geneCount[5];
 601     Double_t measCount[5];
 602     geneCount[4] = 0;
 603     measCount[4] = 0;
 604
 605     for (Int_t i=0; i<4; ++i)
 606     {
 607       geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram()->Integral();
 608       measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Integral();
 609
 610       geneCount[4] += geneCount[i];
 611       measCount[4] += measCount[i];
 612
 613       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]);
 614     }
 615
 616     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]);
 617
 618     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]);
 619
 620     TList* collection = new TList;
 621
 622     // skip "other" particle correction here
 623     // with them has to be dealt differently, maybe just increasing the neutral particles...
 624     for (Int_t i=1; i<3; ++i)
 625       collection->Add(fdNdEtaCorrection[i]);
 626
 627     fdNdEtaCorrection[0]->Merge(collection);
 628     fdNdEtaCorrection[0]->Finish();
 629
 630     delete collection;
 631
 632     // save everything
 633     TFile* file = TFile::Open("new_compositions.root", "UPDATE");
 634     fdNdEtaCorrection[0]->SetName(newName);
 635     fdNdEtaCorrection[0]->SaveHistograms();
 636     //file->Write();
 637     file->Close();
 638
 639     // correct dNdeta distribution with modified correction map
 640     TFile::Open("analysis_esd_raw.root");
 641
 642     dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
 643     fdNdEtaAnalysis->LoadHistograms();
 644
 645     fdNdEtaAnalysis->Finish(fdNdEtaCorrection[0], 0.2, 3, newName);
 646
 647     hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(newName);
 648     hRatios[counter]->SetTitle(newName);
 649     hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default composition}{modified composition}");
 650
 651     if (counter > 0)
 652       hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
 653
 654     file = TFile::Open("new_compositions_analysis.root", "UPDATE");
 655     hRatios[counter]->Write();
 656     file->Close();
 657
 658     delete fdNdEtaAnalysis;
 659
 660     counter++;
 661   }
 662
 663   /*
 664   gROOT->ProcessLine(".L drawPlots.C");
 665
 666   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" };
 667   const char* folderNames[] = { "PythiaRatios", "PiBoosted", "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced" };
 668
 669   Track2Particle1DComposition(fileNames, nCompositions, folderNames);
 670   */
 671 }
 672
 673
 674 void drawSystematics()
 675 {
 676   //Secondaries();
 677   //DrawDifferentSpecies();
 678   //Composition();
 679
 680   Sigma2VertexSimulation();
 681
 682 }
 683
 684 void DrawdNdEtaDifferences()
 685 {
 686   TH1* hists[5];
 687
 688   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
 689   legend->SetFillColor(0);
 690
 691   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
 692   canvas->Divide(2, 1);
 693
 694   canvas->cd(1);
 695
 696   for (Int_t i=0; i<5; ++i)
 697   {
 698     hists[i] = 0;
 699     TFile* file = 0;
 700     TString title;
 701
 702     switch(i)
 703     {
 704       case 0 : file = TFile::Open("systematics_dndeta_reference.root"); title = "standard composition"; break;
 705       case 1 : file = TFile::Open("systematics_dndeta_KBoosted.root"); title = "+ 50% kaons"; break;
 706       case 2 : file = TFile::Open("systematics_dndeta_KReduced.root"); title = "- 50% kaons"; break;
 707       case 3 : file = TFile::Open("systematics_dndeta_pBoosted.root"); title = "+ 50% protons"; break;
 708       case 4 : file = TFile::Open("systematics_dndeta_pReduced.root"); title = "- 50% protons"; break;
 709       default: return;
 710     }
 711
 712     if (file)
 713     {
 714       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
 715       hists[i]->SetTitle("a)");
 716
 717       Prepare1DPlot(hists[i], kFALSE);
 718       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
 719       hists[i]->GetYaxis()->SetRangeUser(6, 7);
 720       hists[i]->SetLineColor(colors[i]);
 721       hists[i]->SetMarkerColor(colors[i]);
 722       hists[i]->SetMarkerStyle(markers[i]);
 723       hists[i]->GetXaxis()->SetLabelOffset(0.015);
 724       hists[i]->GetYaxis()->SetTitleOffset(1.5);
 725       gPad->SetLeftMargin(0.12);
 726       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
 727
 728       legend->AddEntry(hists[i], title);
 729       hists[i]->SetTitle(title);
 730     }
 731   }
 732   legend->Draw();
 733
 734   canvas->cd(2);
 735   gPad->SetLeftMargin(0.14);
 736
 737   TLegend* legend2 = new TLegend(0.73, 0.73, 0.98, 0.98);
 738   legend2->SetFillColor(0);
 739
 740   for (Int_t i=1; i<5; ++i)
 741   {
 742     if (hists[i])
 743     {
 744       legend2->AddEntry(hists[i]);
 745
 746       hists[i]->Divide(hists[0]);
 747       hists[i]->SetTitle("b)");
 748       hists[i]->SetLineColor(colors[i-1]);
 749       hists[i]->SetMarkerColor(colors[i-1]);
 750       hists[i]->GetYaxis()->SetRangeUser(0.95, 1.05);
 751       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
 752       hists[i]->GetYaxis()->SetTitleOffset(1.8);
 753       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
 754     }
 755   }
 756
 757   legend2->Draw();
 758
 759   canvas->SaveAs("particlecomposition_result_detail.gif");
 760
 761   TCanvas* canvas2 = new TCanvas("DrawdNdEtaDifferences2", "DrawdNdEtaDifferences2", 700, 500);
 762
 763   for (Int_t i=1; i<5; ++i)
 764   {
 765     if (hists[i])
 766     {
 767       hists[i]->SetTitle("");
 768       hists[i]->GetYaxis()->SetTitleOffset(1.1);
 769       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
 770     }
 771   }
 772
 773   legend2->Draw();
 774
 775   canvas2->SaveAs("particlecomposition_result.gif");
 776   canvas2->SaveAs("particlecomposition_result.eps");
 777 }
 778
 779 void mergeCorrectionsWithDifferentCrosssections(Int_t correctionTarget = 3, Char_t* correctionFileName="correction_mapprocess-types.root", const char* analysisFileName = "analysis_esd_raw.root", const Char_t* outputFileName="systematics_vtxtrigger_compositions.root") {
 780   //
 781   // Function used to merge standard corrections with vertex
 782   // reconstruction corrections obtained by a certain mix of ND, DD
 783   // and SD events.
 784   //
 785   // the dn/deta spectrum is corrected and the ratios
 786   // (standard to changed x-section) of the different dN/deta
 787   // distributions are saved to a file.
 788   //
 789   // correctionTarget is of type AlidNdEtaCorrection::CorrectionType
 790   //    kINEL = 3
 791   //    kNSD = 4
 792
 793   loadlibs();
 794
 795   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
 796
 797   //Karel:
 798 //     fsd = 0.153 +- 0.031 (0.050 to take into account SD definition) --> change
 799 //     fdd = 0.080 +- 0.050 --> change
 800 //     fnd = 0.767 +- 0.059 --> keep (error small)
 801
 802 //  const Char_t* changes[]  = { "pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdlessddmore25", "sdmoreddless25"};
 803   //Float_t scalesND[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0,  1.0 };
 804   //Float_t scalesDD[] = {1.0, 1.5, 0.5, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
 805   //Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.5, 0.5, 1.5, 0.5, 0.5, 1.5, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
 806   //Float_t scalesDD[] = {1.0, 1.4, 0.6, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 1.25, 0.75, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75};
 807   //Float_t scalesSD[] = {1.0, 1.0, 1.0, 1.4, 0.6, 1.4, 0.6, 0.4, 1.6, 1.0,  1.0,  1.25, 0.75, 1.25, 0.75, 0.75, 1.25};
 808   Int_t nChanges = 9;
 809
 810   const Char_t* changes[]  = { "ua5","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
 811   Float_t scalesND[] = {0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767, 0.767};
 812   Float_t scalesDD[] = {0.080, 0.130, 0.030, 0.080, 0.080, 0.130, 0.030, 0.130, 0.030};
 813   Float_t scalesSD[] = {0.153, 0.153, 0.153, 0.203, 0.103, 0.203, 0.103, 0.103, 0.203};
 814
 815   for (Int_t i=0; i<9; i++)
 816     scalesND[i] = 1.0 - scalesDD[i] - scalesSD[i];
 817
 818   /*
 819   const Char_t* changes[]  = { "pythia", "qgsm", "phojet"};
 820   Float_t scalesND[] = {1.0, 1.10, 1.11};
 821   Float_t scalesSD[] = {1.0, 0.69, 0.86};
 822   Float_t scalesDD[] = {1.0, 0.98, 0.61};
 823   Int_t nChanges = 3;
 824   */
 825
 826   // cross section from Pythia
 827   // 14 TeV!
 828   Float_t sigmaND = 55.2;
 829   Float_t sigmaDD = 9.78;
 830   Float_t sigmaSD = 14.30;
 831
 832   // standard correction
 833   TFile::Open(correctionFileName);
 834   AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
 835   correctionStandard->LoadHistograms();
 836
 837   // dont take vertexreco from this one
 838   correctionStandard->GetVertexRecoCorrection()->Reset();
 839   // dont take triggerbias from this one
 840   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
 841   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
 842   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
 843
 844   AlidNdEtaCorrection* corrections[100];
 845   TH1F* hRatios[100];
 846
 847   Int_t counter = 0;
 848   for (Int_t j=2; j<3; j++) { // j = 0 (change vtx), j = 1 (change trg), j = 2 (change both)
 849
 850     for (Int_t i=0; i<nChanges; i++) {
 851       TFile::Open(correctionFileName);
 852
 853       TString name;
 854       name.Form("dndeta_correction_syst_%s_%s", typeName[j], changes[i]);
 855       AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
 856       current->LoadHistograms("dndeta_correction");
 857       current->Reset();
 858
 859       name.Form("dndeta_correction_ND");
 860       AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
 861       dNdEtaCorrectionND->LoadHistograms();
 862       name.Form("dndeta_correction_DD");
 863       AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
 864       dNdEtaCorrectionDD->LoadHistograms();
 865       name.Form("dndeta_correction_SD");
 866       AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
 867       dNdEtaCorrectionSD->LoadHistograms();
 868
 869       // calculating relative
 870       Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
 871       Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
 872       Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
 873       Float_t total = nd + dd + sd;
 874
 875       nd /= total;
 876       sd /= total;
 877       dd /= total;
 878
 879       Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
 880
 881       Float_t scaleND = scalesND[i] / nd;
 882       Float_t scaleDD = scalesDD[i] / dd;
 883       Float_t scaleSD = scalesSD[i] / sd;
 884
 885       Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);
 886
 887 /*      Float_t nd = 100 * sigmaND/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 888       Float_t dd = 100 * (scalesDD[i]*sigmaDD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 889       Float_t sd = 100 * (scalesSD[i]*sigmaSD)/(sigmaND + (scalesDD[i]*sigmaDD) + (scalesDD[i]*sigmaSD));
 890
 891       printf(Form("%s : ND=%.2f\%, DD=%.2f\%, SD=%.2f\% \n",changes[i],nd,dd,sd));*/
 892       current->SetTitle(Form("ND=%.2f\%,DD=%.2f\%,SD=%.2f\%",scaleND,scaleDD,scaleSD));
 893       current->SetTitle(name);
 894
 895       // scale
 896       if (j == 0 || j == 2)
 897       {
 898         dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
 899         dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
 900         dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
 901       }
 902       if (j == 1 || j == 2)
 903       {
 904         dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
 905         dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scalesDD[i]);
 906         dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scalesSD[i]);
 907
 908         dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
 909         dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
 910         dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);
 911
 912         dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
 913         dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
 914         dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
 915       }
 916
 917       //clear track in correction
 918       dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
 919       dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
 920       dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
 921
 922       TList collection;
 923       collection.Add(correctionStandard);
 924       collection.Add(dNdEtaCorrectionND);
 925       collection.Add(dNdEtaCorrectionDD);
 926       collection.Add(dNdEtaCorrectionSD);
 927
 928       current->Merge(&collection);
 929       current->Finish();
 930
 931       corrections[counter] = current;
 932
 933       // now correct dNdeta distribution with modified correction map
 934       TFile::Open(analysisFileName);
 935
 936       dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD");
 937       fdNdEtaAnalysis->LoadHistograms();
 938
 939       fdNdEtaAnalysis->Finish(current, 0.3, correctionTarget, Form("%d %d", j, i));
 940
 941       name = "ratio";
 942       if (j==0) name.Append("_vetexReco_");
 943       if (j==1) name.Append("_triggerBias_");
 944       if (j==2) name.Append("_vertexReco_triggerBias_");
 945       name.Append(changes[i]);
 946
 947       hRatios[counter] = (TH1F*) fdNdEtaAnalysis->GetdNdEtaHistogram()->Clone(name);
 948
 949       name.Form("ND #times %0.2f DD #times %0.2f, SD #times %0.2f",scaleND,scaleDD,scaleSD);
 950       hRatios[counter]->SetTitle(name.Data());
 951       hRatios[counter]->SetYTitle("dN_{ch}/d#eta ratio #frac{default cross-section}{modified cross-sections}");
 952
 953       if (counter > 0)
 954         hRatios[counter]->Divide(hRatios[0],hRatios[counter],1,1);
 955
 956       delete fdNdEtaAnalysis;
 957
 958       counter++;
 959     }
 960   }
 961
 962   TFile* fout = new TFile(outputFileName,"RECREATE");
 963
 964   // to make everything consistent
 965   hRatios[0]->Divide(hRatios[0],hRatios[0],1,1);
 966
 967   for (Int_t i=0; i<counter; i++)
 968   {
 969     corrections[i]->SaveHistograms();
 970     hRatios[i]->Write();
 971   }
 972
 973   fout->Write();
 974   fout->Close();
 975 }
 976
 977 void CreateCorrectionsWithUA5CrossSections(Int_t origin, const Char_t* correctionFileName="correction_mapprocess-types.root", const Char_t* outputFileName="correction_map2.root") {
 978   //
 979   // Function used to merge standard corrections with vertex
 980   // reconstruction corrections obtained by a certain mix of ND, DD
 981   // and SD events.
 982   //
 983   // origin:
 984   //   -1 = Pythia (test)
 985   //   0 = UA5
 986   //   1 = Data 1.8 TeV
 987   //   2 = Tel-Aviv
 988   //   3 = Durham
 989   //
 990
 991   loadlibs();
 992
 993   const Char_t* typeName[] = { "vertexreco", "trigger", "vtxtrigger" };
 994
 995   Float_t ref_SD = -1;
 996   Float_t ref_DD = -1;
 997   Float_t ref_ND = -1;
 998
 999   switch (origin)
1000   {
1001     case -1: // Pythia, as test
1002       ref_SD = 0.223788;
1003       ref_DD = 0.123315;
1004       ref_ND = 0.652897;
1005       break;
1006
1007     case 0: // UA5
1008       ref_SD = 0.153;
1009       ref_DD = 0.080;
1010       ref_ND = 0.767;
1011       break;
1012
1013     case 1: // data 1.8 TeV
1014       ref_SD = 0.152;
1015       ref_DD = 0.092;
1016       ref_ND = 1 - ref_SD - ref_DD;
1017       break;
1018
1019     case 2: // tel-aviv model
1020       ref_SD = 0.171;
1021       ref_DD = 0.094;
1022       ref_ND = 1 - ref_SD - ref_DD;
1023       break;
1024
1025     case 3: // durham model
1026       ref_SD = 0.190;
1027       ref_DD = 0.125;
1028       ref_ND = 1 - ref_SD - ref_DD;
1029       break;
1030
1031     default:
1032       return;
1033   }
1034
1035   //Karel (UA5):
1036 //     fsd = 0.153 +- 0.031
1037 //     fdd = 0.080 +- 0.050
1038 //     fnd = 0.767 +- 0.059
1039
1040 //       Karel (1.8 TeV):
1041 //
1042 //       Tel-Aviv model Sd/Inel = 0.171           Dd/Inel = 0.094
1043 //       Durham model   Sd/Inel = 0.190           Dd/Inel = 0.125
1044 //       Data           Sd/Inel = 0.152 +- 0.030  Dd/Inel = 0.092 +- 0.45
1045
1046
1047
1048   // standard correction
1049   TFile::Open(correctionFileName);
1050   AlidNdEtaCorrection* correctionStandard = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1051   correctionStandard->LoadHistograms();
1052
1053   // dont take vertexreco from this one
1054   correctionStandard->GetVertexRecoCorrection()->Reset();
1055   // dont take triggerbias from this one
1056   correctionStandard->GetTriggerBiasCorrectionINEL()->Reset();
1057   correctionStandard->GetTriggerBiasCorrectionNSD()->Reset();
1058   correctionStandard->GetTriggerBiasCorrectionND()->Reset();
1059
1060   AlidNdEtaCorrection* corrections[100];
1061   TH1F* hRatios[100];
1062
1063   Int_t counter = 0;
1064
1065   TFile::Open(correctionFileName);
1066
1067   AlidNdEtaCorrection* current = new AlidNdEtaCorrection("dndeta_correction_ua5", "dndeta_correction_ua5");
1068   current->LoadHistograms("dndeta_correction");
1069   current->Reset();
1070
1071   TString name;
1072   name.Form("dndeta_correction_ND");
1073   AlidNdEtaCorrection* dNdEtaCorrectionND = new AlidNdEtaCorrection(name,name);
1074   dNdEtaCorrectionND->LoadHistograms();
1075   name.Form("dndeta_correction_DD");
1076   AlidNdEtaCorrection* dNdEtaCorrectionDD = new AlidNdEtaCorrection(name,name);
1077   dNdEtaCorrectionDD->LoadHistograms();
1078   name.Form("dndeta_correction_SD");
1079   AlidNdEtaCorrection* dNdEtaCorrectionSD = new AlidNdEtaCorrection(name,name);
1080   dNdEtaCorrectionSD->LoadHistograms();
1081
1082   // calculating relative
1083   Float_t nd = dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1084   Float_t dd = dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1085   Float_t sd = dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral();
1086   Float_t total = nd + dd + sd;
1087
1088   nd /= total;
1089   sd /= total;
1090   dd /= total;
1091
1092   Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
1093
1094   Float_t scaleND = ref_ND / nd;
1095   Float_t scaleDD = ref_DD / dd;
1096   Float_t scaleSD = ref_SD / sd;
1097
1098   Printf("ND=%.2f, DD=%.2f, SD=%.2f",scaleND, scaleDD, scaleSD);
1099
1100   // scale
1101   dNdEtaCorrectionND->GetVertexRecoCorrection()->Scale(scaleND);
1102   dNdEtaCorrectionDD->GetVertexRecoCorrection()->Scale(scaleDD);
1103   dNdEtaCorrectionSD->GetVertexRecoCorrection()->Scale(scaleSD);
1104
1105   dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->Scale(scaleND);
1106   dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->Scale(scaleDD);
1107   dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->Scale(scaleSD);
1108
1109   dNdEtaCorrectionND->GetTriggerBiasCorrectionNSD()->Scale(scaleND);
1110   dNdEtaCorrectionDD->GetTriggerBiasCorrectionNSD()->Scale(scaleDD);
1111   dNdEtaCorrectionSD->GetTriggerBiasCorrectionNSD()->Scale(scaleSD);
1112
1113   dNdEtaCorrectionND->GetTriggerBiasCorrectionND()->Scale(scaleND);
1114   dNdEtaCorrectionDD->GetTriggerBiasCorrectionND()->Scale(scaleDD);
1115   dNdEtaCorrectionSD->GetTriggerBiasCorrectionND()->Scale(scaleSD);
1116
1117   //clear track in correction
1118   dNdEtaCorrectionND->GetTrack2ParticleCorrection()->Reset();
1119   dNdEtaCorrectionDD->GetTrack2ParticleCorrection()->Reset();
1120   dNdEtaCorrectionSD->GetTrack2ParticleCorrection()->Reset();
1121
1122   TList collection;
1123   collection.Add(correctionStandard);
1124   collection.Add(dNdEtaCorrectionND);
1125   collection.Add(dNdEtaCorrectionDD);
1126   collection.Add(dNdEtaCorrectionSD);
1127
1128   current->Merge(&collection);
1129   current->Finish();
1130
1131   TFile* fout = new TFile(outputFileName,"RECREATE");
1132   current->SaveHistograms();
1133
1134   fout->Write();
1135   fout->Close();
1136
1137   Printf("Trigger efficiencies:");
1138   Printf("ND: %.2f %%", 100.0 * dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1139   Printf("SD: %.2f %%", 100.0 * dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionSD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1140   Printf("DD: %.2f %%", 100.0 * dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1141   Printf("INEL: %.2f %%", 100.0 * current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() / current->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral());
1142   Printf("NSD: %.2f %%", 100.0 * (dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral() + dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetMeasuredHistogram()->Integral()) / (dNdEtaCorrectionND->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral() + dNdEtaCorrectionDD->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetGeneratedHistogram()->Integral()));
1143 }
1144
1145 DrawTriggerEfficiency(Char_t* fileName) {
1146
1147   gStyle->SetOptStat(0);
1148   gStyle->SetOptTitle(0);
1149   gStyle->SetOptFit(0);
1150
1151   gStyle->SetTextSize(0.04);
1152   gStyle->SetTitleSize(0.05,"xyz");
1153   //gStyle->SetTitleFont(133, "xyz");
1154   //gStyle->SetLabelFont(133, "xyz");
1155   //gStyle->SetLabelSize(17, "xyz");
1156   gStyle->SetLabelOffset(0.01, "xyz");
1157
1158   gStyle->SetTitleOffset(1.1, "y");
1159   gStyle->SetTitleOffset(1.1, "x");
1160   gStyle->SetEndErrorSize(0.0);
1161
1162   //##############################################
1163
1164   //making canvas and pads
1165   TCanvas *c = new TCanvas("trigger_eff", "",600,500);
1166
1167   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
1168
1169   p1->SetBottomMargin(0.15);
1170   p1->SetTopMargin(0.03);
1171   p1->SetLeftMargin(0.15);
1172   p1->SetRightMargin(0.03);
1173
1174   p1->SetGridx();
1175   p1->SetGridy();
1176
1177   p1->Draw();
1178   p1->cd();
1179
1180   gSystem->Load("libPWG0base");
1181
1182   AlidNdEtaCorrection* corrections[4];
1183   AliCorrectionMatrix2D* triggerBiasCorrection[4];
1184
1185   TH1F* hTriggerEffInv[4];
1186   TH1F* hTriggerEff[4];
1187
1188   Char_t* names[] = {"triggerBiasND", "triggerBiasDD", "triggerBiasSD", "triggerBiasINEL"};
1189
1190   for (Int_t i=0; i<4; i++) {
1191     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1192     corrections[i]->LoadHistograms(fileName, names[i]);
1193
1194     triggerBiasCorrection[i] = corrections[i]->GetTriggerBiasCorrectionINEL();
1195
1196
1197     hTriggerEffInv[i] = triggerBiasCorrection[i]->Get1DCorrection();
1198     hTriggerEff[i]    = (TH1F*)hTriggerEffInv[i]->Clone();
1199
1200     for (Int_t b=0; b<=hTriggerEff[i]->GetNbinsX(); b++) {
1201       hTriggerEff[i]->SetBinContent(b,1);
1202       hTriggerEff[i]->SetBinError(b,0);
1203     }
1204     hTriggerEff[i]->Divide(hTriggerEff[i],hTriggerEffInv[i]);
1205     hTriggerEff[i]->Scale(100);
1206   }
1207
1208   Int_t colors[] = {2,3,4,1};
1209   Float_t effs[4];
1210   for (Int_t i=0; i<4; i++) {
1211     hTriggerEff[i]->Fit("pol0","","",-20,20);
1212     effs[i] = ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->GetParameter(0);
1213     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineWidth(1);
1214     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineStyle(2);
1215     ((TF1*)hTriggerEff[i]->GetListOfFunctions()->At(0))->SetLineColor(colors[i]);
1216     cout << effs[i] << endl;
1217   }
1218
1219
1220   Char_t* text[] = {"ND", "DD", "SD", "INEL"};
1221   TLatex* latex[4];
1222
1223   TH2F* null = new TH2F("","",100,-25,35,100,0,110);
1224   null->SetXTitle("Vertex z [cm]");
1225   null->GetXaxis()->CenterTitle(kTRUE);
1226   null->SetYTitle("Trigger efficiency [%]");
1227   null->GetYaxis()->CenterTitle(kTRUE);
1228   null->Draw();
1229
1230
1231   for (Int_t i=0; i<4; i++) {
1232     hTriggerEff[i]->SetLineWidth(2);
1233     hTriggerEff[i]->SetLineColor(colors[i]);
1234
1235     hTriggerEff[i]->Draw("same");
1236
1237     latex[i] = new TLatex(22,effs[i]-1.5, Form("%s (%0.1f)",text[i],effs[i]));
1238     latex[i]->SetTextColor(colors[i]);
1239     latex[i]->Draw();
1240   }
1241
1242 }
1243
1244
1245 DrawSpectraPID(Char_t* fileName) {
1246
1247   gSystem->Load("libPWG0base");
1248
1249   Char_t* names[] = {"correction_0", "correction_1", "correction_2", "correction_3"};
1250   AlidNdEtaCorrection* corrections[4];
1251   AliCorrectionMatrix3D* trackToPartCorrection[4];
1252
1253   TH1F* measuredPt[4];
1254   TH1F* generatedPt[4];
1255   TH1F* ratioPt[4];
1256
1257   for (Int_t i=0; i<4; i++) {
1258     corrections[i] = new AlidNdEtaCorrection(names[i], names[i]);
1259     corrections[i]->LoadHistograms(fileName, names[i]);
1260
1261     trackToPartCorrection[i] = corrections[i]->GetTrack2ParticleCorrection();
1262
1263     Int_t binX1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(-10);
1264     Int_t binX2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetXaxis()->FindBin(10);
1265     Int_t binY1 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(-1);
1266     Int_t binY2 = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->GetYaxis()->FindBin(1);
1267
1268     measuredPt[i]  = (TH1F*)trackToPartCorrection[i]->GetMeasuredHistogram()->ProjectionZ(Form("m_%d",i),binX1,binX2,binY1,binY2);
1269     generatedPt[i] = (TH1F*)trackToPartCorrection[i]->GetGeneratedHistogram()->ProjectionZ(Form("g_%d",i),binX1,binX2,binY1,binY2);
1270     ratioPt[i] = (TH1F*)generatedPt[i]->Clone(Form("r_%d",i));
1271     ratioPt[i]->Divide(measuredPt[i], generatedPt[i], 1,1,"B");
1272   }
1273
1274   ratioPt[0]->Draw();
1275
1276   for (Int_t i=0; i<3; i++) {
1277     ratioPt[i]->SetLineColor(i+1);
1278     ratioPt[i]->SetLineWidth(2);
1279
1280     ratioPt[i]->Draw("same");
1281
1282   }
1283
1284   return;
1285   measuredPt[0]->SetLineColor(2);
1286   measuredPt[0]->SetLineWidth(5);
1287
1288   measuredPt[0]->Draw();
1289   generatedPt[0]->Draw("same");
1290 }
1291
1292 void changePtSpectrum(const char* fileName = "analysis_mc.root", Float_t ptCutOff = 0.2, const char* fileName2 = 0)
1293 {
1294   Float_t factor = 0.5;
1295
1296   TFile* file = TFile::Open(fileName);
1297   TH1F* hist = dynamic_cast<TH1F*> (file->Get("dndeta_check_pt")->Clone());
1298
1299   TH1* hist2 = 0;
1300   if (fileName2)
1301   {
1302     file2 = TFile::Open(fileName2);
1303     hist2 = dynamic_cast<TH1*> (file2->Get("dndeta_check_pt")->Clone());
1304     hist2->Scale(hist->GetBinContent(hist->FindBin(ptCutOff)) / hist2->GetBinContent(hist2->FindBin(ptCutOff)));
1305   }
1306
1307   //hist->Scale(1.0 / hist->Integral());
1308
1309   //hist->Rebin(3);
1310   //hist->Scale(1.0/3);
1311
1312   TH1F* clone1 = dynamic_cast<TH1F*> (hist->Clone("clone1"));
1313   TH1F* clone2 = dynamic_cast<TH1F*> (hist->Clone("clone2"));
1314
1315   TH1F* scale1 =  dynamic_cast<TH1F*> (hist->Clone("scale1"));
1316   TH1F* scale2 =  dynamic_cast<TH1F*> (hist->Clone("scale2"));
1317
1318   for (Int_t i=1; i <= hist->GetNbinsX(); ++i)
1319   {
1320     if (hist->GetBinCenter(i) > ptCutOff)
1321     {
1322       scale1->SetBinContent(i, 1);
1323       scale2->SetBinContent(i, 1);
1324     }
1325     else
1326     {
1327       // 90 % at pt = 0, 0% at pt = ptcutoff
1328       scale1->SetBinContent(i, 1 - (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
1329
1330       // 110% at pt = 0, ...
1331       scale2->SetBinContent(i, 1 + (ptCutOff - hist->GetBinCenter(i)) / ptCutOff * factor);
1332     }
1333     scale1->SetBinError(i, 0);
1334     scale2->SetBinError(i, 0);
1335   }
1336
1337   /*
1338   new TCanvas;
1339   scale1->Draw();
1340   scale2->SetLineColor(kRed);
1341   scale2->Draw("SAME");
1342   */
1343
1344   clone1->Multiply(scale1);
1345   clone2->Multiply(scale2);
1346
1347   Prepare1DPlot(hist);
1348   Prepare1DPlot(clone1);
1349   Prepare1DPlot(clone2);
1350
1351   /*hist->SetMarkerStyle(markers[0]);
1352   clone1->SetMarkerStyle(markers[0]);
1353   clone2->SetMarkerStyle(markers[0]);*/
1354
1355   hist->SetTitle(";p_{T} in GeV/c;dN_{ch}/dp_{T} in c/GeV");
1356   hist->GetXaxis()->SetRangeUser(0, 0.5);
1357   hist->GetYaxis()->SetRangeUser(0.01, clone2->GetMaximum() * 1.1);
1358   hist->GetYaxis()->SetTitleOffset(1);
1359
1360   TCanvas* canvas = new TCanvas("c", "c", 600, 600);
1361   gPad->SetGridx();
1362   gPad->SetGridy();
1363   gPad->SetTopMargin(0.05);
1364   gPad->SetRightMargin(0.05);
1365   gPad->SetBottomMargin(0.12);
1366   hist->Draw("H");
1367   clone1->SetLineColor(kRed);
1368   clone1->Draw("HSAME");
1369   clone2->SetLineColor(kBlue);
1370   clone2->Draw("HSAME");
1371   hist->Draw("HSAME");
1372
1373   if (hist2)
1374   {
1375     Prepare1DPlot(hist2);
1376     hist2->SetLineStyle(2);
1377     hist2->Draw("HSAME");
1378   }
1379
1380   Float_t fraction =  hist->Integral(hist->GetXaxis()->FindBin(ptCutOff), hist->GetNbinsX()) / hist->Integral(1, hist->GetNbinsX());
1381   Float_t fraction1 = clone1->Integral(clone1->GetXaxis()->FindBin(ptCutOff), clone1->GetNbinsX()) / clone1->Integral(1, clone1->GetNbinsX());
1382   Float_t fraction2 = clone2->Integral(clone2->GetXaxis()->FindBin(ptCutOff), clone2->GetNbinsX()) / clone2->Integral(1, clone2->GetNbinsX());
1383
1384   printf("%f %f %f\n", fraction, fraction1, fraction2);
1385   printf("Rel. %f %f\n", fraction1 / fraction, fraction2 / fraction);
1386
1387   //canvas->SaveAs("changePtSpectrum.gif");
1388   canvas->SaveAs("changePtSpectrum.eps");
1389 }
1390
1391 void FractionBelowPt()
1392 {
1393   gSystem->Load("libPWG0base");
1394
1395   AlidNdEtaCorrection* fdNdEtaCorrection[4];
1396
1397   TFile::Open("systematics.root");
1398
1399   for (Int_t i=0; i<4; ++i)
1400   {
1401     TString name;
1402     name.Form("correction_%d", i);
1403     fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
1404     fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
1405   }
1406
1407   Double_t geneCount[5];
1408   Double_t measCount[5];
1409   geneCount[4] = 0;
1410   measCount[4] = 0;
1411
1412   for (Int_t i=0; i<4; ++i)
1413   {
1414     TH3F* hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram();
1415     geneCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10),
1416                                   hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8),
1417                                   1, hist->GetZaxis()->FindBin(0.3));
1418
1419     hist = (TH3F*) fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram();
1420     measCount[i] = hist->Integral(hist->GetXaxis()->FindBin(-10), hist->GetXaxis()->FindBin(10), hist->GetYaxis()->FindBin(-0.8), hist->GetYaxis()->FindBin(0.8), 1, hist->GetZaxis()->FindBin(0.3));
1421
1422     geneCount[4] += geneCount[i];
1423     measCount[4] += measCount[i];
1424
1425     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]);
1426   }
1427
1428   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]);
1429
1430   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]);
1431 }
1432
1433
1434 mergeCorrectionsMisalignment(Char_t* alignedFile = "correction_map_aligned.root",
1435                                              Char_t* misalignedFile = "correction_map_misaligned.root",
1436                                              Char_t* outputFileName="correction_map_misaligned_single.root")
1437 {
1438   //
1439   // from the aligned and misaligned corrections, 3 new corrections are created
1440   // in these new corrections only one of the corrections (track2particle, vertex, trigger)
1441   // is taken from the misaligned input to allow study of the effect on the different
1442   // corrections
1443
1444   gSystem->Load("libPWG0base");
1445
1446   const Char_t* typeName[] = { "track2particle", "vertex", "trigger" };
1447
1448   AlidNdEtaCorrection* corrections[3];
1449   for (Int_t j=0; j<3; j++) { // j = 0 (track2particle), j = 1 (vertex), j = 2 (trigger)
1450     AlidNdEtaCorrection* alignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1451     alignedCorrection->LoadHistograms(alignedFile);
1452
1453     AlidNdEtaCorrection* misalignedCorrection = new AlidNdEtaCorrection("dndeta_correction","dndeta_correction");
1454     misalignedCorrection->LoadHistograms(misalignedFile);
1455
1456     TString name;
1457     name.Form("dndeta_correction_alignment_%s", typeName[j]);
1458     AlidNdEtaCorrection* current = new AlidNdEtaCorrection(name, name);
1459
1460     switch (j)
1461     {
1462       case 0:
1463         alignedCorrection->GetTrack2ParticleCorrection()->Reset();
1464         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1465         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1466         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1467         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1468         break;
1469
1470       case 1:
1471         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1472         misalignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1473         misalignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1474         misalignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1475         alignedCorrection->GetVertexRecoCorrection()->Reset();
1476         break;
1477
1478       case 2:
1479         misalignedCorrection->GetTrack2ParticleCorrection()->Reset();
1480         alignedCorrection->GetTriggerBiasCorrectionNSD()->Reset();
1481         alignedCorrection->GetTriggerBiasCorrectionND()->Reset();
1482         alignedCorrection->GetTriggerBiasCorrectionINEL()->Reset();
1483         misalignedCorrection->GetVertexRecoCorrection()->Reset();
1484         break;
1485
1486       default:
1487         return;
1488     }
1489
1490     TList collection;
1491     collection.Add(misalignedCorrection);
1492     collection.Add(alignedCorrection);
1493
1494     current->Merge(&collection);
1495     current->Finish();
1496
1497     corrections[j] = current;
1498   }
1499
1500   TFile* fout = new TFile(outputFileName, "RECREATE");
1501
1502   for (Int_t i=0; i<3; i++)
1503     corrections[i]->SaveHistograms();
1504
1505   fout->Write();
1506   fout->Close();
1507 }
1508
1509
1510 void DrawdNdEtaDifferencesAlignment()
1511 {
1512   TH1* hists[5];
1513
1514   TLegend* legend = new TLegend(0.3, 0.73, 0.70, 0.98);
1515   legend->SetFillColor(0);
1516
1517   TCanvas* canvas = new TCanvas("DrawdNdEtaDifferences", "DrawdNdEtaDifferences", 1000, 500);
1518   canvas->Divide(2, 1);
1519
1520   canvas->cd(1);
1521
1522   for (Int_t i=0; i<5; ++i)
1523   {
1524     hists[i] = 0;
1525     TFile* file = 0;
1526     TString title;
1527
1528     switch(i)
1529     {
1530       case 0 : file = TFile::Open("systematics_misalignment_aligned.root"); title = "aligned"; break;
1531       case 1 : file = TFile::Open("systematics_misalignment_misaligned.root"); title = "fully misaligned"; break;
1532       case 2 : file = TFile::Open("systematics_misalignment_track2particle.root"); title = "only track2particle"; break;
1533       case 3 : file = TFile::Open("systematics_misalignment_vertex.root"); title = "only vertex rec."; break;
1534       case 4 : file = TFile::Open("systematics_misalignment_trigger.root"); title = "only trigger bias"; break;
1535       default: return;
1536     }
1537
1538     if (file)
1539     {
1540       hists[i] = (TH1*) file->Get("dndeta/dndeta_dNdEta_corrected_2");
1541       hists[i]->SetTitle("");
1542
1543       Prepare1DPlot(hists[i], kFALSE);
1544       hists[i]->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
1545       hists[i]->GetYaxis()->SetRangeUser(6, 7);
1546       hists[i]->SetLineWidth(1);
1547       hists[i]->SetLineColor(colors[i]);
1548       hists[i]->SetMarkerColor(colors[i]);
1549       hists[i]->SetMarkerStyle(markers[i]);
1550       hists[i]->GetXaxis()->SetLabelOffset(0.015);
1551       hists[i]->GetYaxis()->SetTitleOffset(1.5);
1552       gPad->SetLeftMargin(0.12);
1553       hists[i]->DrawCopy(((i > 0) ? "SAME" : ""));
1554
1555       legend->AddEntry(hists[i], title);
1556       hists[i]->SetTitle(title);
1557     }
1558   }
1559   legend->Draw();
1560
1561   canvas->cd(2);
1562   gPad->SetLeftMargin(0.14);
1563
1564   TLegend* legend2 = new TLegend(0.63, 0.73, 0.98, 0.98);
1565   legend2->SetFillColor(0);
1566
1567   for (Int_t i=1; i<5; ++i)
1568   {
1569     if (hists[i])
1570     {
1571       legend2->AddEntry(hists[i]);
1572
1573       hists[i]->Divide(hists[0]);
1574       hists[i]->SetTitle("b)");
1575       hists[i]->GetYaxis()->SetRangeUser(0.9, 1.1);
1576       hists[i]->GetYaxis()->SetTitle("Ratio to standard composition");
1577       hists[i]->GetYaxis()->SetTitleOffset(1.8);
1578       hists[i]->DrawCopy(((i > 1) ? "SAME" : ""));
1579     }
1580   }
1581
1582   legend2->Draw();
1583
1584   canvas->SaveAs("misalignment_result.eps");
1585   canvas->SaveAs("misalignment_result.gif");
1586 }
1587
1588 void EvaluateMultiplicityEffect()
1589 {
1590   gSystem->Load("libPWG0base");
1591
1592   AlidNdEtaCorrection* fdNdEtaCorrectionLow[4];
1593   AlidNdEtaCorrection* fdNdEtaCorrectionHigh[4];
1594
1595   TFile::Open("systematics-low-multiplicity.root");
1596
1597   for (Int_t i=0; i<4; ++i)
1598   {
1599     TString name;
1600     name.Form("correction_%d", i);
1601     fdNdEtaCorrectionLow[i] = new AlidNdEtaCorrection(name, name);
1602     fdNdEtaCorrectionLow[i]->LoadHistograms("systematics-low-multiplicity.root", name);
1603   }
1604
1605   TList list;
1606   for (Int_t i=1; i<4; ++i)
1607     list.Add(fdNdEtaCorrectionLow[i]);
1608
1609   fdNdEtaCorrectionLow[0]->Merge(&list);
1610   fdNdEtaCorrectionLow[0]->Finish();
1611
1612   TFile::Open("systematics-high-multiplicity.root");
1613
1614   for (Int_t i=0; i<4; ++i)
1615   {
1616     TString name;
1617     name.Form("correction_%d", i);
1618     fdNdEtaCorrectionHigh[i] = new AlidNdEtaCorrection(name, name);
1619     fdNdEtaCorrectionHigh[i]->LoadHistograms("systematics-high-multiplicity.root", name);
1620   }
1621
1622   TList list2;
1623   for (Int_t i=1; i<4; ++i)
1624     list2.Add(fdNdEtaCorrectionHigh[i]);
1625
1626   fdNdEtaCorrectionHigh[0]->Merge(&list2);
1627   fdNdEtaCorrectionHigh[0]->Finish();
1628
1629   TH1F* outputLow = new TH1F("Track2ParticleLow", "Track2Particle at low multiplicity", 200, 0, 2);
1630   TH1F* outputHigh = new TH1F("Track2ParticleHigh", "Track2Particle at high multiplicity", 200, 0, 2);
1631
1632   TH3F* hist = fdNdEtaCorrectionLow[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1633   TH3F* hist2 = fdNdEtaCorrectionHigh[0]->GetTrack2ParticleCorrection()->GetCorrectionHistogram();
1634   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
1635     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
1636       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
1637       //for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
1638       {
1639         if (hist->GetBinContent(x, y, z) > 0)
1640           outputLow->Fill(hist->GetBinContent(x, y, z));
1641         //if (hist->GetBinContent(x, y, z) == 1)
1642         //  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));
1643
1644         if (hist2->GetBinContent(x, y, z) > 0)
1645           outputHigh->Fill(hist2->GetBinContent(x, y, z));
1646       }
1647
1648   TCanvas* canvas = new TCanvas("EvaluateMultiplicityEffect", "EvaluateMultiplicityEffect", 1000, 500);
1649   canvas->Divide(2, 1);
1650
1651   canvas->cd(1);
1652   outputLow->Draw();
1653   outputLow->Fit("gaus", "0");
1654   outputLow->GetFunction("gaus")->SetLineColor(2);
1655   outputLow->GetFunction("gaus")->DrawCopy("SAME");
1656
1657   canvas->cd(2);
1658   outputHigh->Draw();
1659   outputHigh->Fit("gaus", "0");
1660   outputHigh->GetFunction("gaus")->DrawCopy("SAME");
1661
1662   canvas->SaveAs(Form("%s.gif", canvas->GetName()));
1663 }
1664
1665 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")
1666 {
1667     gSystem->Load("libPWG0base");
1668
1669         AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(correctionMapFolder, correctionMapFolder);
1670     dNdEtaCorrection->LoadHistograms(correctionMapFile, correctionMapFolder);
1671
1672     dNdEtaCorrection->GetTrack2ParticleCorrection()->PlotBinErrors(xmin, xmax, ymin, ymax, zmin, zmax)->Draw();
1673 }
1674
1675
1676 void runStudy(const char* baseCorrectionMapFile = "correction_map.root", const char* baseCorrectionMapFolder = "dndeta_correction", const char* changedCorrectionMapFile = "correction_map.root", const char* changedCorrectionMapFolder = "dndeta_correction", const char* dataFile = "analysis_esd_raw.root", const char* output = "analysis_esd_syst.root")
1677 {
1678   gSystem->Load("libPWG0base");
1679
1680   TFile* file = TFile::Open(output, "RECREATE");
1681
1682   const Int_t max = 5;
1683   dNdEtaAnalysis* fdNdEtaAnalysis[5];
1684
1685   new TCanvas;
1686   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
1687   legend->SetFillColor(0);
1688
1689   for (Int_t i = 0; i < max; ++i)
1690   {
1691     TFile::Open(baseCorrectionMapFile);
1692     AlidNdEtaCorrection* baseCorrection = new AlidNdEtaCorrection(baseCorrectionMapFolder, baseCorrectionMapFolder);
1693     baseCorrection->LoadHistograms();
1694
1695     AlidNdEtaCorrection::CorrectionType correctionType = AlidNdEtaCorrection::kNone;
1696     const char* name = 0;
1697
1698     TFile::Open(changedCorrectionMapFile);
1699     switch (i)
1700     {
1701       case 0 :
1702         name = "default";
1703         break;
1704
1705       case 1 :
1706         baseCorrection->GetTrack2ParticleCorrection()->LoadHistograms(Form("%s/Track2Particle", changedCorrectionMapFolder));
1707         name = "Track2Particle";
1708         break;
1709
1710       case 2 :
1711         baseCorrection->GetVertexRecoCorrection()->LoadHistograms(Form("%s/VertexReconstruction", changedCorrectionMapFolder));
1712         name = "VertexReco";
1713         break;
1714
1715       case 3 :
1716         baseCorrection->GetTriggerBiasCorrectionINEL()->LoadHistograms(Form("%s/TriggerBias_MBToINEL", changedCorrectionMapFolder));
1717         name = "TriggerBias_MBToINEL";
1718         break;
1719
1720       case 4 :
1721         baseCorrection->LoadHistograms(changedCorrectionMapFolder);
1722         name = "all";
1723         break;
1724
1725       default: return;
1726     }
1727
1728     TFile::Open(dataFile);
1729     fdNdEtaAnalysis[i] = new dNdEtaAnalysis(name, name);
1730     fdNdEtaAnalysis[i]->LoadHistograms("dndeta");
1731
1732     fdNdEtaAnalysis[i]->Finish(baseCorrection, 0.3, AlidNdEtaCorrection::kINEL);
1733     file->cd();
1734     fdNdEtaAnalysis[i]->SaveHistograms();
1735
1736     TH1* hist = fdNdEtaAnalysis[i]->GetdNdEtaHistogram(0);
1737     hist->SetLineColor(colors[i]);
1738     hist->SetMarkerColor(colors[i]);
1739     hist->SetMarkerStyle(markers[i]+1);
1740     hist->DrawCopy((i == 0) ? "" : "SAME");
1741     legend->AddEntry(hist, name);
1742   }
1743
1744   legend->Draw();
1745 }
1746
1747 void ChangePtInCorrection(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
1748 {
1749   loadlibs();
1750   if (!TFile::Open(fileName))
1751     return;
1752
1753   AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection(dirName, dirName);
1754   if (!dNdEtaCorrection->LoadHistograms())
1755     return;
1756
1757   dNdEtaCorrection->Finish();
1758
1759   AliCorrection* corr = dNdEtaCorrection->GetCorrection(AlidNdEtaCorrection::kTrack2Particle);
1760
1761   Printf(">>>> Before");
1762   corr->PrintInfo(0);
1763
1764   Float_t factor = 0.5;
1765   Float_t ptCutOff = 0.2;
1766
1767   TH3* gene = corr->GetTrackCorrection()->GetGeneratedHistogram();
1768   TH3* meas = corr->GetTrackCorrection()->GetMeasuredHistogram();
1769
1770   for (Int_t z = 1; z <= gene->GetZaxis()->FindBin(ptCutOff - 0.01); z++)
1771   {
1772     Float_t localFactor = 1 - (ptCutOff - gene->GetZaxis()->GetBinCenter(z)) / ptCutOff * factor;
1773     Printf("%f %f", gene->GetZaxis()->GetBinCenter(z), localFactor);
1774     for (Int_t x = 1; x <= gene->GetNbinsX(); x++)
1775       for (Int_t y = 1; y <= gene->GetNbinsY(); y++)
1776       {
1777         gene->SetBinContent(x, y, z, gene->GetBinContent(x, y, z) * localFactor);
1778         meas->SetBinContent(x, y, z, meas->GetBinContent(x, y, z) * localFactor);
1779       }
1780   }
1781
1782   dNdEtaCorrection->Finish();
1783
1784   Printf(">>>> After");
1785   corr->PrintInfo(0);
1786 }