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