PWGUD/selectors/dNdEta/drawSystematicsNew.C

   1
   2 Int_t markers[] = {20,21,22,23,28,29};
   3 Int_t colors[]  = {1,2,3,4,6,8,102};
   4
   5 void loadlibs()
   6 {
   7   gSystem->Load("libTree");
   8   gSystem->Load("libVMC");
   9
  10   gSystem->Load("libSTEERBase");
  11   gSystem->Load("libANALYSIS");
  12   gSystem->Load("libPWG0base");
  13 }
  14
  15 void InitPad()
  16 {
  17   if (!gPad)
  18     return;
  19
  20   gPad->Range(0, 0, 1, 1);
  21   gPad->SetLeftMargin(0.15);
  22   //gPad->SetRightMargin(0.05);
  23   //gPad->SetTopMargin(0.13);
  24   //gPad->SetBottomMargin(0.1);
  25
  26   gPad->SetGridx();
  27   gPad->SetGridy();
  28 }
  29
  30 void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
  31 {
  32   //gROOT->ProcessLine(".L drawPlots.C");
  33   gSystem->Load("libPWG0base");
  34
  35   const char* fileNames[] = { "systematics.root", "systematics.root", "systematics.root", "correction_map.root" };
  36   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "dndeta_correction" };
  37   const char* legendNames[] = { "#pi", "K", "p", "standard" };
  38   Int_t folderCount = 3;
  39
  40   TH2F* h2DCorrections[4];
  41   TH1F* h1DCorrections[4];
  42   for (Int_t i=0; i<4; i++) {
  43     TFile::Open(fileNames[i]);
  44     AlidNdEtaCorrection* correctionTmp = new AlidNdEtaCorrection(folderNames[i],folderNames[i]);
  45     correctionTmp->LoadHistograms();
  46
  47     //    h2DCorrections[i] = correctionTmp->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get2DCorrectionHistogram("yz",-1,1);
  48
  49     h1DCorrections[i] = correctionTmp->GetTrack2ParticleCorrection()->GetTrackCorrection()->Get1DCorrectionHistogram("z",-10,10,0.8,0.8);
  50   }
  51
  52   TH2F* null = new TH2F("","",100,0.1,10,100,0.5,9.99);
  53   null->SetXTitle("p_{T} (GeV/c)");
  54   null->SetXTitle("Correction");
  55
  56   null->Draw();
  57
  58   //h1DCorrections[0]->SetMaximum(5);
  59   //h1DCorrections[0]->Draw();
  60   //  h2DCorrections[0]->Draw("colz");
  61   for (Int_t i=1; i<4; i++) {
  62     h1DCorrections[i]->Draw("same");
  63   }
  64
  65 }
  66
  67 TPad* DrawChange(Bool_t spd, const char* basename, const char** changes, Int_t nChanges, Int_t nDraw, Int_t* colors, const char** names = 0, Float_t scale = 0.10)
  68 {
  69   Float_t etaMax = 1.05;
  70   if (spd)
  71     etaMax = 1.49;
  72
  73   TH1F* hRatios[100];
  74   for(Int_t i=0; i<nChanges; i++) {
  75     Printf("%d", i);
  76     hRatios[i] = (TH1F*)gFile->Get(Form("%s%s",basename,changes[i]));
  77     hRatios[i]->SetLineWidth(1);
  78     hRatios[i]->SetMarkerStyle(22);
  79     hRatios[i]->SetMarkerSize(0.8);
  80
  81     Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
  82     Printf("%s %s: %.2f %%" , changes[i], hRatios[i]->GetTitle(), (average - 1) * 100);
  83   }
  84
  85   TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
  86   p->SetRightMargin(0.2);
  87   p->SetLeftMargin(0.13);
  88
  89   TH2F* null = new TH2F("","",100,-etaMax, etaMax,100,1. - scale,1. + scale);
  90   null->GetXaxis()->SetTitle("#eta");
  91   null->GetYaxis()->SetTitle(hRatios[0]->GetYaxis()->GetTitle());
  92   null->Draw();
  93
  94   line = new TLine(-etaMax, 1, etaMax, 1);
  95   line->Draw();
  96
  97   TLatex* text[100];
  98
  99   for(Int_t i=1; i<nDraw; i++) {
 100     hRatios[i]->SetLineColor(colors[i]);
 101     hRatios[i]->SetMarkerColor(colors[i]);
 102     hRatios[i]->Draw("HISTPL SAME");
 103
 104     TString str(hRatios[i]->GetTitle());
 105
 106     if (names)
 107       str = names[i];
 108
 109     text[i] = new TLatex(etaMax + 0.03,hRatios[i]->GetBinContent(hRatios[i]->FindBin(etaMax-0.1))-0.002,str.Data());
 110     text[i]->SetTextAlign(11);
 111     text[i]->SetTextColor(colors[i]);
 112     text[i]->Draw();
 113   }
 114
 115   return p;
 116 }
 117
 118 TPad* DrawChangeRatios(Bool_t spd, TFile* file1, TFile* file2, const char* basename, const char** changes, Int_t nChanges, Int_t nDraw, Int_t* colors, const char** names = 0, Float_t scale = 0.10)
 119 {
 120   Float_t etaMax = 1.05;
 121   if (spd)
 122     etaMax = 1.79;
 123
 124   TH1F* hRatios[100];
 125   for(Int_t i=0; i<nChanges; i++) {
 126     hRatios[i] = (TH1F*)file1->Get(Form("%s%s",basename,changes[i]));
 127     hRatios[i]->SetLineWidth(1);
 128     hRatios[i]->SetMarkerStyle(22);
 129     hRatios[i]->SetMarkerSize(0.8);
 130
 131     hRatio2 = (TH1F*)file2->Get(Form("%s%s",basename,changes[i]));
 132     hRatios[i]->Divide(hRatio2);
 133
 134     Float_t average = hRatios[i]->Integral(hRatios[i]->FindBin(-1), hRatios[i]->FindBin(1)) / (hRatios[i]->FindBin(1) - hRatios[i]->FindBin(-1) + 1);
 135     Printf("%s: %.2f %%" , hRatios[i]->GetTitle(), (average - 1) * 100);
 136   }
 137
 138   TPad* p = DrawCanvasAndPad("syst_changeInXsection",700,400);
 139   p->SetRightMargin(0.2);
 140   p->SetLeftMargin(0.13);
 141
 142   TH2F* null = new TH2F("","",100,-etaMax, etaMax,100,1. - scale,1. + scale);
 143   null->GetXaxis()->SetTitle("#eta");
 144   null->GetYaxis()->SetTitle(hRatios[0]->GetYaxis()->GetTitle());
 145   null->Draw();
 146
 147   line = new TLine(-etaMax, 1, etaMax, 1);
 148   line->Draw();
 149
 150   TLatex* text[100];
 151
 152   for(Int_t i=1; i<nDraw; i++) {
 153     hRatios[i]->SetLineColor(colors[i]);
 154     hRatios[i]->SetMarkerColor(colors[i]);
 155     hRatios[i]->Draw("HISTPL SAME");
 156
 157     TString str(hRatios[i]->GetTitle());
 158
 159     if (names)
 160       str = names[i];
 161
 162     text[i] = new TLatex(etaMax + 0.03,hRatios[i]->GetBinContent(hRatios[i]->FindBin(etaMax-0.1))-0.002,str.Data());
 163     text[i]->SetTextAlign(11);
 164     text[i]->SetTextColor(colors[i]);
 165     text[i]->Draw();
 166   }
 167
 168   return p;
 169 }
 170
 171 void DrawEffectOfChangeInCrossSection(Bool_t spd = kFALSE, const char* fileName = "systematics_vtxtrigger_compositions.root")
 172 {
 173   TFile::Open(fileName);
 174
 175 //  const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
 176   const Char_t* changes[]  = { "default","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdlessddmore", "sdmoreddless" };
 177   //const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
 178   //const Int_t nChanges = 3;
 179   Int_t colors[] = {1,1,4,1,2,2,4,2,1};
 180
 181   c = DrawChange(spd, "ratio_vertexReco_triggerBias_", changes, 9, 9, colors, 0);
 182   c->SaveAs("cross_sections.eps");
 183 }
 184
 185 void DrawEffectOfChangeInCrossSection2Files(Bool_t spd = kFALSE, const char* fileName1 = "systematics_vtxtrigger_compositions.root", const char* fileName2)
 186 {
 187   file1 = TFile::Open(fileName1);
 188   file2 = TFile::Open(fileName2);
 189
 190   const Char_t* changes[]  = {"pythia","ddmore","ddless","sdmore","sdless", "dmore", "dless", "sdmoreddless", "sdlessddmore", "ddmore25","ddless25","sdmore25","sdless25", "dmore25", "dless25", "sdmoreddless25", "sdlessddmore25" };
 191   //const Char_t* changes[]  = { "pythia", "qgsm", "phojet" };
 192   //const Int_t nChanges = 3;
 193   Int_t colors[] = {1,1,4,1,2,2,4,2,1};
 194
 195   c = DrawChangeRatios(spd, file1, file2, "ratio_vertexReco_triggerBias_", changes, 17, 9, colors, 0);
 196   c->SaveAs("cross_sections.eps");
 197 }
 198
 199 void DrawEffectOfChangeInComposition(Bool_t spd = kFALSE, const char* fileName = "new_compositions_analysis.root")
 200 {
 201   TFile::Open(fileName);
 202
 203   const Char_t* changes[]  = { "KBoosted", "KReduced", "pBoosted", "pReduced", "KBoostedpBoosted", "KReducedpReduced", "KBoostedpReduced", "KReducedpBoosted", "othersBoosted", "othersReduced" };
 204   const char*   names[]    = { "K #times 1.3", "K #times 0.7", "p #times 1.3", "p #times 0.7", "K #times 1.3, p #times 1.3", "K #times 0.7, p #times 0.7", "K #times 1.3, p #times 0.7", "K #times 0.7, p #times 1.3", "O #times 1.3", "O #times 0.7" };
 205   //const Char_t* changes[]  = { "PythiaRatios", "PiBoosted",      "PiReduced", "KBoosted", "KReduced", "pBoosted", "pReduced", "othersBoosted", "othersReduced" };
 206   //const char*   names[]    = { "",             "#pi #times 1.5", "#pi #times 0.5", "K #times 1.5", "K #times 0.5", "p #times 1.5", "p #times 0.5",  "others #times 1.5", "others #times 0.5" };
 207   Int_t colors[] = {1,1,2,2,1,2,1,4,4,1,1};
 208
 209   c = DrawChange(spd, "", changes, 10, 10, colors, names, 0.03);
 210   //c->SaveAs("compositions.eps");
 211 }
 212
 213 TPad* DrawCanvasAndPad(const Char_t* name, Int_t sizeX=600, Int_t sizeY=500) {
 214
 215   gStyle->SetOptStat(0);
 216   gStyle->SetOptTitle(0);
 217   gStyle->SetOptFit(0);
 218
 219   gStyle->SetTextSize(0.04);
 220   gStyle->SetTitleSize(0.05,"xyz");
 221   //gStyle->SetTitleFont(133, "xyz");
 222   //gStyle->SetLabelFont(133, "xyz");
 223   //gStyle->SetLabelSize(17, "xyz");
 224   gStyle->SetLabelOffset(0.01, "xyz");
 225
 226   gStyle->SetTitleOffset(1.1, "y");
 227   gStyle->SetTitleOffset(1.1, "x");
 228   gStyle->SetEndErrorSize(0.0);
 229
 230   //##############################################
 231
 232   //making canvas and pads
 233   TCanvas *c = new TCanvas(name,name,sizeX,sizeY);
 234
 235   TPad* p1 = new TPad("pad1","", 0, 0.0, 1.0, 1.0, 0, 0, 0);
 236
 237   p1->SetBottomMargin(0.15);
 238   p1->SetTopMargin(0.03);
 239   p1->SetLeftMargin(0.15);
 240   p1->SetRightMargin(0.03);
 241
 242   p1->SetGridx();
 243   p1->SetGridy();
 244
 245   p1->Draw();
 246   p1->cd();
 247
 248   return p1;
 249 }
 250
 251 void MisalignmentShowRawTrackPlots(const char* dirName = "fdNdEtaAnalysisESD")
 252 {
 253   loadlibs();
 254
 255   TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
 256   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
 257   fdNdEtaAnalysis->LoadHistograms();
 258
 259   TFile* file2 = TFile::Open("fullA-sim/analysis_esd_raw.root");
 260   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
 261   fdNdEtaAnalysis2->LoadHistograms();
 262
 263   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
 264   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
 265
 266   // normalize to number of events;
 267   TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
 268   TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
 269   Int_t event1Count = event1->Integral();
 270   Int_t event2Count = event2->Integral();
 271   track1->Scale(1.0 / event1Count);
 272   track2->Scale(1.0 / event2Count);
 273
 274   const Float_t innerLimit = 0.49;
 275   const Float_t outerLimit = 0.99;
 276
 277   track1->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
 278   track2->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
 279   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze1");
 280   TH1* fullRange = gROOT->FindObject("track1_ze1_div_track2_ze1");
 281
 282   track1->GetYaxis()->SetRangeUser(-innerLimit, innerLimit);
 283   track2->GetYaxis()->SetRangeUser(-innerLimit, innerLimit);
 284   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze2");
 285   TH1* central = gROOT->FindObject("track1_ze2_div_track2_ze2");
 286   central->SetLineColor(1);
 287   central->SetMarkerStyle(21);
 288
 289   for (Int_t x=1; x<track1->GetXaxis()->GetNbins(); ++x)
 290     for (Int_t y=track1->GetYaxis()->FindBin(-innerLimit); y<track1->GetYaxis()->FindBin(innerLimit); ++y)
 291       for (Int_t z=1; z<track1->GetZaxis()->GetNbins(); ++z)
 292       {
 293         track1->SetBinContent(x, y, z, 0);
 294         track1->SetBinError(x, y, z, 0);
 295         track2->SetBinContent(x, y, z, 0);
 296         track2->SetBinError(x, y, z, 0);
 297       }
 298
 299   track1->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
 300   track2->GetYaxis()->SetRangeUser(-outerLimit, outerLimit);
 301   AliPWG0Helper::CreateDividedProjections(track1, track2, "ze3");
 302   TH1* peripheral = gROOT->FindObject("track1_ze3_div_track2_ze3");
 303   peripheral->SetLineColor(2);
 304   peripheral->SetMarkerStyle(22);
 305   peripheral->SetMarkerColor(2);
 306
 307   TH2* tmp = new TH2F("tmp", ";p_{T} [GeV/c]      ;#frac{tracks full misalignment during rec.}{tracks ideal misalignment during rec.}", 1, 0.1, 10, 1, 0.9, 1.3);
 308
 309   tmp->SetStats(kFALSE);
 310   //tmp->GetXaxis()->SetNoExponent();
 311
 312   Float_t ptStart = 0.1;
 313
 314   fullRange->GetXaxis()->SetRangeUser(ptStart, 9.9);
 315   central->GetXaxis()->SetRangeUser(ptStart, 9.9);
 316   peripheral->GetXaxis()->SetRangeUser(ptStart, 9.9);
 317
 318   TCanvas* canvas = new TCanvas("MisalignmentShowRawTrackPlots", "MisalignmentShowRawTrackPlots", 700, 400);
 319   gPad->SetLogx();
 320   gPad->SetGridx();
 321   gPad->SetGridy();
 322
 323   TLegend* legend = new TLegend(0.2, 0.7, 0.4, 0.8);
 324
 325   legend->AddEntry(central, "|#eta| < 0.5");
 326   legend->AddEntry(peripheral, "0.5 < |#eta| < 1.0");
 327
 328   legend->SetFillColor(0);
 329
 330   tmp->Draw();
 331   //fullRange->Draw("SAME");
 332   central->Draw("SAME");
 333   peripheral->Draw("SAME");
 334
 335   legend->Draw();
 336
 337   canvas->SaveAs("syst_mis_ntracks.eps");
 338 }
 339
 340
 341 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files, const char** dirs, const char** names, Int_t* histID)
 342 {
 343   loadlibs();
 344   gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C");
 345
 346   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
 347   canvas->Divide(2, 1);
 348   canvas->cd(2)->SetGridx();
 349   canvas->cd(2)->SetGridy();
 350
 351   TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
 352   legend->SetFillColor(0);
 353
 354   TH1* base = 0;
 355
 356   for (Int_t i = 0; i < n; ++i)
 357   {
 358     TFile::Open(files[i]);
 359
 360     dNdEtaAnalysis* tmp = new dNdEtaAnalysis(dirs[i], dirs[i]);
 361     tmp->LoadHistograms();
 362
 363     TH1* hist = tmp->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
 364
 365     if (i == 0)
 366       base = hist;
 367
 368     legend->AddEntry(hist, names[i]);
 369
 370     hist->SetMarkerColor(colors[i]);
 371     hist->SetMarkerStyle(markers[i]);
 372
 373     canvas->cd(1);
 374     hist->DrawCopy((i == 0) ? "" : "SAME");
 375
 376     if (i != 0)
 377     {
 378       PrintIntegratedDeviation(hist, base, names[i]);
 379
 380       canvas->cd(2);
 381       hist->Divide(hist, base, 1, 1);
 382       hist->GetYaxis()->SetRangeUser(0.9, 1.1);
 383       hist->DrawCopy((i == 1) ? "" : "SAME");
 384     }
 385   }
 386
 387   legend->Draw();
 388
 389   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 390 }
 391
 392 void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files1, const char** files2, const char** dirs1, const char** dirs2, const char** names, Int_t* histID)
 393 {
 394   gSystem->Load("libPWG0base");
 395
 396   TCanvas* canvas = new TCanvas(canvasName, canvasName, 700, 400);
 397   canvas->SetLeftMargin(0.12);
 398
 399   TLegend* legend = new TLegend(0.35, 0.7, 0.65, 0.85);
 400   legend->SetFillColor(0);
 401
 402   for (Int_t i = 0; i < n; ++i)
 403   {
 404     TFile::Open(files1[i]);
 405
 406     dNdEtaAnalysis* tmp1 = new dNdEtaAnalysis(dirs1[i], dirs1[i]);
 407     tmp1->LoadHistograms();
 408
 409     TFile::Open(files2[i]);
 410
 411     dNdEtaAnalysis* tmp2 = new dNdEtaAnalysis(dirs2[i], dirs2[i]);
 412     tmp2->LoadHistograms();
 413
 414     TH1* hist1 = tmp1->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
 415     TH1* hist2 = tmp2->GetdNdEtaPtCutOffCorrectedHistogram(histID[i]);
 416
 417     TH1* division = hist1->Clone();
 418
 419     division->Divide(hist1, hist2, 1, 1, "B");
 420
 421     division->SetMarkerColor(colors[i]);
 422     division->SetMarkerStyle(markers[i]);
 423
 424     legend->AddEntry(division, names[i]);
 425
 426     division->SetTitle("");
 427     division->GetYaxis()->SetTitle("#frac{dN_{ch}/d#eta using MC vtx}{dN_{ch}/d#eta using ESD vtx}");
 428     division->SetStats(kFALSE);
 429     division->GetYaxis()->SetTitleOffset(1.3);
 430     division->GetXaxis()->SetRangeUser(-0.99, 0.99);
 431     division->GetYaxis()->SetRangeUser(0.981, 1.02);
 432     division->DrawCopy((i == 0) ? "" : "SAME");
 433   }
 434
 435   gPad->SetGridx();
 436   gPad->SetGridy();
 437
 438   legend->Draw();
 439
 440   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 441 }
 442
 443 void MaterialBudgetChange()
 444 {
 445   const char* files[] =
 446     { "Material-normal-mcvtx/analysis_esd.root",
 447       "Material-increased-mcvtx/analysis_esd.root",
 448       "Material-decreased-mcvtx/analysis_esd.root" };
 449
 450   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
 451   const char* names[] = { "no change", "+ 10 % material", "- 10 % material" };
 452   Int_t hist[] = { 0, 0, 0 };
 453
 454   drawdNdEtaRatios("MaterialBudgetChange", 3, files, dirs, names, hist);
 455 }
 456
 457 void MisalignmentChange()
 458 {
 459   const char* files[] =
 460     { "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-sim-mcvtx/analysis_esd.root",
 461       "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-simrec-mcvtx/analysis_esd.root" };
 462
 463   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
 464   const char* names[] = { "no change", "increased in sim/rec", "increased in sim" };
 465   Int_t hist[] = { 0, 0, 0 };
 466
 467   drawdNdEtaRatios("MisalignmentChange", 2, files, dirs, names, hist);
 468 }
 469
 470 void dNdEtaVertexRanges()
 471 {
 472   const char* files[] =
 473     { "analysis_esd.root",
 474       "analysis_esd.root",
 475       "analysis_esd.root" };
 476
 477   const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
 478   const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
 479   Int_t hist[] = { 0, 1, 2 };
 480
 481   drawdNdEtaRatios("dNdEtaVertexRanges", 3, files, dirs, names, hist);
 482 }
 483
 484 void vertexShiftStudy(Int_t histID)
 485 {
 486   const char* files[] = { "maps/idealA/mc-vertex/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.05/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.1/analysis_esd.root", "maps/idealA/mc-vertex-shift-dep/analysis_esd.root" };
 487   const char* dirs[] = { "dndeta", "dndeta", "dndeta", "dndeta", "dndeta" };
 488   const char* names[] = { "mc vtx", "esd vtx", "+ 0.05 cm", "+ 0.1 cm", "old vtx shift" };
 489   Int_t hist[] = { histID, histID, histID, histID, histID };
 490
 491   drawdNdEtaRatios("syst_vertex_shift1", 5, files, dirs, names, hist);
 492
 493   const char* files1[] = { "maps/idealA/mc-vertex/analysis_esd.root", "maps/idealA/mc-vertex/analysis_esd.root", "maps/idealA/mc-vertex/analysis_esd.root" };
 494   const char* files2[] = { "results/idealC-idealA/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "results/idealC-idealA/analysis_esd.root" };
 495   const char* dirs1[] = { "dndeta", "dndeta", "dndeta"};
 496   const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
 497   Int_t hist[] = { 0, 1, 2 };
 498
 499   drawdNdEtaRatios("syst_vertex_shift2", 3, files1, files2, dirs, dirs, names, hist);
 500 }
 501
 502 void vertexShift()
 503 {
 504   TFile::Open("vertex.root");
 505
 506   TH2* hist = gFile->Get("fVertexCorr");
 507   TProfile* prof = hist->ProfileX();
 508
 509   prof->SetStats(kFALSE);
 510   prof->SetTitle(";MC vtx-z [cm];mean (ESD vtx-z - MC vtx-z) [cm]");
 511   prof->GetYaxis()->SetTitleOffset(1.2);
 512
 513   prof->SetLineWidth(2);
 514
 515   TCanvas* canvas = new TCanvas("syst_vertex_shift", "syst_vertex_shift", 700, 400);
 516
 517   gPad->SetGridx();
 518   gPad->SetGridy();
 519
 520   prof->Draw();
 521
 522   canvas->SaveAs(Form("%s.eps", canvas->GetName()));
 523 }
 524
 525 void CompareRawTrackPlots(const char* fileName1, const char* fileName2, Float_t ptCut = 0.0, Int_t multCut = 1)
 526 {
 527   loadlibs();
 528
 529   const char* dirName = "fdNdEtaAnalysisESD";
 530
 531   TFile* file = TFile::Open(fileName1);
 532   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
 533   fdNdEtaAnalysis->LoadHistograms();
 534
 535   TFile* file2 = TFile::Open(fileName2);
 536   dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
 537   fdNdEtaAnalysis2->LoadHistograms();
 538
 539   TH3* track1 = fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track1");
 540   TH3* track2 = fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("track2");
 541
 542   TH2* event1 = fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram();
 543   TH2* event2 = fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram();
 544   Int_t event1Count = event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), multCut, event1->GetNbinsY() + 1);
 545   Int_t event2Count = event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), multCut, event1->GetNbinsY() + 1);
 546
 547   Float_t nTrack1 = track1->Integral(track1->GetXaxis()->FindBin(-9.9), track1->GetXaxis()->FindBin(9.9), track1->GetYaxis()->FindBin(-0.79), track1->GetYaxis()->FindBin(0.79), track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
 548   Float_t nTrack2 = track2->Integral(track2->GetXaxis()->FindBin(-9.9), track2->GetXaxis()->FindBin(9.9), track2->GetYaxis()->FindBin(-0.79), track2->GetYaxis()->FindBin(0.79), track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
 549
 550   Printf("%d tracks in %d events in first sample; %d tracks in %d events in second sample", (Int_t) nTrack1, event1Count, (Int_t) nTrack2, event2Count);
 551
 552   // normalize to number of events;
 553   nTrack1 /= event1Count;
 554   nTrack2 /= event2Count;
 555
 556   Printf("There are %.2f tracks/event in the first sample and %.2f tracks/event in the second sample. %.2f %% difference (with pt cut at %.2f GeV/c)", nTrack1, nTrack2, 100.0 * (nTrack1 - nTrack2) / nTrack1, ptCut);
 557
 558   gROOT->cd();
 559
 560   AliPWG0Helper::CreateDividedProjections(track1, track2);
 561
 562   new TCanvas; gROOT->FindObject("track1_yx_div_track2_yx")->Draw("COLZ");
 563   new TCanvas; gROOT->FindObject("track1_zx_div_track2_zx")->Draw("COLZ");
 564   new TCanvas; gROOT->FindObject("track1_zy_div_track2_zy")->Draw("COLZ");
 565
 566   for (Int_t i=0; i<3; i++)
 567   {
 568     char c = 'x' + (char) i;
 569
 570     /*proj1 = track1->Project3D(Form("%ce2", c));
 571     proj2 = track2->Project3D(Form("%ce2", c));
 572     AliPWG0Helper::NormalizeToBinWidth(proj1);
 573     AliPWG0Helper::NormalizeToBinWidth(proj2);
 574
 575     new TCanvas;
 576     proj1->DrawCopy();
 577     proj2->SetLineColor(2);
 578     proj2->SetMarkerColor(2);
 579     proj2->DrawCopy("SAME");*/
 580
 581     AliPWG0Helper::CreateDividedProjections(track1, track2, Form("%ce2", c));
 582     TH1* pt = gROOT->FindObject(Form("track1_%ce2_div_track2_%ce2", c, c));
 583     new TCanvas; pt->DrawCopy();
 584     gPad->SetGridx(); gPad->SetGridy();
 585   }
 586
 587   event1_x = event1->ProjectionX("event1_x");
 588   event1_y = event1->ProjectionY("event1_y");
 589   event2_x = event2->ProjectionX("event2_x");
 590   event2_y = event2->ProjectionY("event2_y");
 591
 592   new TCanvas; event1_x->DrawCopy(); event2_x->SetLineColor(2); event2_x->DrawCopy("SAME");
 593
 594   event1_x->Divide(event2_x);
 595   event1_y->Divide(event2_y);
 596
 597   new TCanvas; event1_x->DrawCopy();
 598   new TCanvas; event1_y->DrawCopy();
 599
 600   event1->Divide(event2);
 601   new TCanvas;
 602   event1->Draw("COLZ");
 603   event1->SetMinimum(0.5);
 604   event1->SetMaximum(2);
 605
 606
 607 }
 608
 609 void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta", Float_t ptCut = 0.3)
 610 {
 611   loadlibs();
 612
 613   TFile* file = TFile::Open(fileName);
 614   dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
 615   fdNdEtaAnalysis->LoadHistograms();
 616   fdNdEtaAnalysis->GetData()->PrintInfo(ptCut);
 617 }
 618
 619 Double_t ConvSigma1To2D(Double_t sigma)
 620 {
 621   return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
 622 }
 623
 624 Double_t ConvDistance1DTo2D(Double_t distance)
 625 {
 626   return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
 627 }
 628
 629 Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
 630 {
 631   Double_t count = 0;
 632
 633   //nSigma = ConvSigma1To2D(nSigma);
 634
 635   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
 636     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
 637     {
 638       Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
 639       Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
 640
 641       Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
 642
 643       d = ConvDistance1DTo2D(d);
 644
 645       if (d < nSigma)
 646         count += tracks->GetBinContent(x, y);
 647     }
 648
 649   return count;
 650 }
 651
 652 TH2F* Sigma2VertexGaussianTracksHist()
 653 {
 654   TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
 655
 656   TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
 657   gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
 658
 659   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
 660     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
 661       tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
 662
 663   //normalize
 664   tracks->Scale(1.0 / tracks->Integral());
 665
 666   return tracks;
 667 }
 668
 669 TH1F* Sigma2VertexGaussian()
 670 {
 671   TH2F* tracks = Sigma2VertexGaussianTracksHist();
 672
 673   TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
 674   canvas->Divide(2, 2);
 675
 676   canvas->cd(1);
 677   tracks->Draw("COLZ");
 678
 679   TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
 680   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
 681     ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
 682   ratio->SetMarkerStyle(21);
 683
 684   canvas->cd(2);
 685   ratio->DrawCopy("P");
 686
 687   TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
 688   Double_t sigma3 = Sigma2VertexCount(tracks, 4);
 689   for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
 690     ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
 691   ratio2->SetMarkerStyle(21);
 692
 693   canvas->cd(3);
 694   ratio2->DrawCopy("P");
 695
 696   canvas->SaveAs("Sigma2Vertex.eps");
 697
 698   return ratio2;
 699 }
 700
 701 TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
 702 {
 703   TFile* file = TFile::Open(fileName);
 704
 705   TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
 706   TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
 707   if (!sigmavertex || !sigmavertexPrim)
 708   {
 709     printf("Could not read histogram(s)\n");
 710     return;
 711   }
 712
 713   // calculate ratio
 714   TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
 715
 716   // calculate contamination
 717   TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
 718   contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
 719
 720   for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
 721   {
 722     ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
 723     contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
 724   }
 725
 726   // print stats
 727   for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
 728   {
 729     Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
 730     Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
 731     Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
 732     Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
 733   }
 734
 735   TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
 736   canvas->Divide(2, 1);
 737
 738   canvas->cd(1);
 739   sigmavertex->SetMarkerStyle(21);
 740   sigmavertex->Draw("P");
 741
 742   canvas->cd(2);
 743   ratio->SetMarkerStyle(21);
 744   ratio->DrawCopy("P");
 745
 746   contamination->DrawCopy("SAME");
 747
 748   TH1F** returnContainer = new TH1F*[2];
 749   returnContainer[0] = ratio;
 750   returnContainer[1] = contamination;
 751
 752   return returnContainer;
 753 }
 754
 755 void Sigma2VertexCompare(const char* fileName = "systematics.root")
 756 {
 757   TH1F* ratio1 = Sigma2VertexGaussian();
 758
 759   TH1F** hists = Sigma2VertexSimulation(fileName);
 760   TH1F* ratio2 = hists[0];
 761   TH1F* contamination = hists[1];
 762
 763   ratio1->SetStats(kFALSE);
 764   ratio2->SetStats(kFALSE);
 765
 766   ratio1->SetMarkerStyle(0);
 767   ratio2->SetMarkerStyle(0);
 768
 769   ratio1->SetLineWidth(2);
 770   ratio2->SetLineWidth(2);
 771
 772   TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
 773   legend->SetFillColor(0);
 774   legend->AddEntry(ratio1, "Gaussian");
 775   legend->AddEntry(ratio2, "Simulation");
 776   legend->AddEntry(contamination, "1 + Contamination");
 777
 778   ratio2->SetTitle("");
 779   ratio2->GetYaxis()->SetTitleOffset(1.5);
 780   ratio2->GetXaxis()->SetRangeUser(2, 5);
 781
 782   TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
 783   InitPad();
 784
 785   ratio2->SetMarkerStyle(21);
 786   ratio1->SetMarkerStyle(22);
 787
 788   ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
 789   ratio2->SetLineColor(kRed);
 790   ratio2->SetMarkerColor(kRed);
 791   ratio2->Draw("PL");
 792   ratio1->Draw("SAMEPL");
 793
 794   contamination->Draw("SAME");
 795
 796   legend->Draw();
 797
 798   canvas->SaveAs("Sigma2VertexCompare.eps");
 799 }