PWG4/macros/underlyingevent/correct.C

   1 const char* objectName = "PWG4_LeadingTrackUE/histosLeadingTrackUE";
   2 Float_t gLeadingpTMin = 0.51;
   3 Int_t gUEHist = 0;
   4 Bool_t gCache = 0;
   5 void* gFirst = 0;
   6 void* gSecond = 0;
   7 Float_t gForceRange = -1;
   8 Int_t gEnergy = 900;
   9 Int_t gRegion = 0;
  10
  11 //const char* objectName = "PWG4_LeadingTrackUE/histosLeadingTrackUE_filterbit32";
  12
  13 void SetRanges(TAxis* axis)
  14 {
  15   if (strcmp(axis->GetTitle(), "leading p_{T} (GeV/c)") == 0)
  16     axis->SetRangeUser(0, (gEnergy == 900) ? 10 : 25);
  17
  18   if (strcmp(axis->GetTitle(), "multiplicity") == 0)
  19     axis->SetRangeUser(0, 10);
  20 }
  21
  22 void SetRanges(TH1* hist)
  23 {
  24   SetRanges(hist->GetXaxis());
  25   SetRanges(hist->GetYaxis());
  26   SetRanges(hist->GetZaxis());
  27 }
  28
  29 void Prepare1DPlot(TH1* hist)
  30 {
  31   hist->SetLineWidth(2);
  32   hist->SetStats(kFALSE);
  33
  34   hist->GetXaxis()->SetLabelOffset(0.02);
  35   hist->GetXaxis()->SetTitleOffset(1.3);
  36   hist->GetYaxis()->SetTitleOffset(1.3);
  37
  38   SetRanges(hist);
  39 }
  40
  41 TH1* GetSystematicUncertainty(TH1* corr, TH1* trackHist)
  42 {
  43   if (!trackHist)
  44   {
  45     systError = (TH1*) corr->Clone("systError");
  46     systError->Reset();
  47   }
  48   else  // for dphi evaluation
  49     systError = new TH1F("systError", "", 100, 0, 50);
  50
  51   Float_t constantUnc = 0;
  52
  53   // particle composition
  54   constantUnc += 0.8 ** 2;
  55
  56   // tpc efficiency
  57   if (gEnergy == 900 && gLeadingpTMin < 1.0)
  58     constantUnc += 1.0 ** 2;
  59   else if (gEnergy == 900 && gLeadingpTMin < 1.5)
  60     constantUnc += 0.5 ** 2;
  61   if (gEnergy == 7000 && gLeadingpTMin < 1.0)
  62     constantUnc += 1.0 ** 2;
  63   else if (gEnergy == 7000 && gLeadingpTMin < 1.5)
  64     constantUnc += 0.6 ** 2;
  65
  66   // track cuts
  67   if (gEnergy == 900 && gLeadingpTMin < 1.0)
  68     constantUnc += 2.5 ** 2;
  69   else if (gEnergy == 900 && gLeadingpTMin < 1.5)
  70     constantUnc += 2.0 ** 2;
  71   if (gEnergy == 7000)
  72     constantUnc += 3.0 ** 2;
  73
  74   // difference corrected with pythia and phojet
  75   if (gEnergy == 900 && gLeadingpTMin < 1.0)
  76     constantUnc += 0.6 ** 2;
  77   else if (gEnergy == 900 && gLeadingpTMin < 1.5)
  78     constantUnc += 0.8 ** 2;
  79
  80   if (gEnergy == 7000 && gLeadingpTMin < 1.0)
  81   {
  82     if (gUEHist == 0)
  83       constantUnc += 0.6 ** 2;
  84     if (gUEHist == 1)
  85       constantUnc += 0.8 ** 2;
  86     if (gUEHist == 2)
  87       constantUnc += 1.0 ** 2;
  88   }
  89   else if (gEnergy == 7000 && gLeadingpTMin < 1.5)
  90     constantUnc += 1.0 ** 2;
  91
  92   for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
  93     systError->SetBinContent(bin, constantUnc);
  94
  95   // mis id bias
  96   if (gUEHist == 0 || gUEHist == 2)
  97     systError->Fill(0.75, 4.0 ** 2);
  98   if (gUEHist == 1)
  99     systError->Fill(0.75, 5.0 ** 2);
 100
 101   if (gEnergy == 900)
 102   {
 103     if (gLeadingpTMin < 1.0)
 104       systError->Fill(1.25, 1.0 ** 2);
 105     else if (gLeadingpTMin < 1.5)
 106       systError->Fill(1.25, 2.0 ** 2);
 107   }
 108
 109   // non-closure in MC
 110   if (gEnergy == 900)
 111     for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
 112       systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 1.0 ** 2);
 113
 114   if (gEnergy == 7000)
 115   {
 116     if (gUEHist == 0 && gUEHist == 1)
 117       systError->Fill(0.75, 2.0 ** 2);
 118     if (gUEHist == 2)
 119       systError->Fill(0.75, 1.2 ** 2);
 120   }
 121
 122   // vertex efficiency
 123   systError->Fill(0.75, 1.0 ** 2);
 124
 125   // strangeness
 126   for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
 127   {
 128     if (gEnergy == 900)
 129       systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 0.5 ** 2);
 130     if (gEnergy == 7000 && systError->GetXaxis()->GetBinCenter(bin) < 1.5)
 131       systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 2.0 ** 2);
 132     else if (gEnergy == 7000)
 133       systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 1.0 ** 2);
 134   }
 135
 136   for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
 137     systError->SetBinContent(bin, TMath::Sqrt(systError->GetBinContent(bin)));
 138
 139   if (trackHist)
 140   {
 141     //new TCanvas; trackHist->Draw();
 142     //new TCanvas; systError->DrawCopy("");
 143
 144     Float_t uncFlat = 0;
 145     for (Int_t i=1; i<=trackHist->GetNbinsX(); i++)
 146       uncFlat += trackHist->GetBinContent(i) * systError->GetBinContent(systError->FindBin(trackHist->GetBinCenter(i)));
 147     uncFlat /= trackHist->Integral();
 148
 149     systError = (TH1F*) corr->Clone("systError");
 150     systError->Reset();
 151
 152     for (Int_t i=1; i<=systError->GetNbinsX(); i++)
 153       systError->SetBinContent(i, uncFlat);
 154
 155     //new TCanvas; systError->DrawCopy("");
 156   }
 157
 158   systError->SetFillColor(kGray);
 159   systError->SetFillStyle(1001);
 160   systError->SetMarkerStyle(0);
 161   systError->SetLineColor(0);
 162
 163   return systError;
 164 }
 165
 166 void DrawRatio(TH1* corr, TH1* mc, const char* name, TH1* syst = 0)
 167 {
 168   mc->SetLineColor(2);
 169
 170   TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
 171   canvas3->Range(0, 0, 1, 1);
 172
 173   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
 174   pad1->Draw();
 175
 176   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
 177   pad2->Draw();
 178
 179   pad1->SetRightMargin(0.01);
 180   pad2->SetRightMargin(0.01);
 181   pad1->SetTopMargin(0.05);
 182   pad1->SetLeftMargin(0.13);
 183   pad2->SetLeftMargin(0.13);
 184   pad2->SetBottomMargin(0.22);
 185
 186   // no border between them
 187   pad1->SetBottomMargin(0);
 188   pad2->SetTopMargin(0);
 189
 190   pad1->cd();
 191   pad1->SetGridx();
 192   pad1->SetGridy();
 193
 194   if (gUEHist != 2)
 195     TLegend* legend = new TLegend(0.15, 0.65, 0.55, 0.90);
 196   else
 197     TLegend* legend = new TLegend(0.55, 0.65, 0.95, 0.90);
 198
 199   legend->SetFillColor(0);
 200   legend->AddEntry(corr, "Corrected");
 201   legend->AddEntry(mc, "MC prediction");
 202   legend->SetTextSize(0.08);
 203
 204   Prepare1DPlot(corr);
 205
 206   TH2F* dummy = new TH2F("dummy", "", 100, corr->GetXaxis()->GetBinLowEdge(1), corr->GetXaxis()->GetBinUpEdge(corr->GetNbinsX()), 1000, 0, TMath::Max(corr->GetMaximum(), mc->GetMaximum()) * 1.1);
 207   dummy->SetStats(kFALSE);
 208   dummy->SetXTitle(corr->GetXaxis()->GetTitle());
 209   dummy->SetYTitle(corr->GetYaxis()->GetTitle());
 210   dummy->GetYaxis()->SetTitleOffset(1);
 211   Prepare1DPlot(dummy);
 212
 213   dummy->GetXaxis()->SetLabelSize(0.08);
 214   dummy->GetYaxis()->SetLabelSize(0.08);
 215   dummy->GetXaxis()->SetTitleSize(0.08);
 216   dummy->GetYaxis()->SetTitleSize(0.08);
 217   dummy->GetYaxis()->SetTitleOffset(0.8);
 218
 219   if (gForceRange > 0)
 220     dummy->GetYaxis()->SetRangeUser(0, gForceRange);
 221
 222   dummy->DrawCopy();
 223
 224   if (gUEHist != 2)
 225   {
 226     const char* regionStr[] = { "Toward Region", "Away Region", "Transverse Region" };
 227     latex = new TLatex(0.65, 0.1, regionStr[gRegion]);
 228     latex->SetTextSize(0.075);
 229     latex->SetNDC();
 230     latex->Draw();
 231   }
 232
 233   if (syst)
 234   {
 235     systError = (TH1*) syst->Clone("corrSystError");
 236     for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
 237     {
 238       systError->SetBinError(bin, corr->GetBinContent(bin) * syst->GetBinContent(bin) / 100);
 239       systError->SetBinContent(bin, corr->GetBinContent(bin));
 240     }
 241
 242     systError->Draw("E2 ][ SAME");
 243   }
 244
 245   corr->Draw("SAME");
 246   mc->Draw("SAME");
 247
 248   legend->Draw();
 249
 250   pad2->cd();
 251   //pad2->SetBottomMargin(0.15);
 252   //pad2->SetGridx();
 253   //pad2->SetGridy();
 254
 255   TH1* ratio = (TH1*) mc->Clone("ratio");
 256   ratio->Divide(corr);
 257
 258   Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
 259   Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
 260
 261   minR = 0.6;
 262   maxR = 1.4;
 263
 264   TH1F dummy3("dummy3", ";;Ratio: MC / corr", 100, corr->GetXaxis()->GetBinLowEdge(1), corr->GetXaxis()->GetBinUpEdge(corr->GetNbinsX()));
 265   dummy3.SetXTitle(corr->GetXaxis()->GetTitle());
 266   Prepare1DPlot(&dummy3);
 267   dummy3.SetStats(kFALSE);
 268   for (Int_t i=1; i<=100; ++i)
 269         dummy3.SetBinContent(i, 1);
 270   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 271   dummy3.SetLineWidth(2);
 272   dummy3.GetXaxis()->SetLabelSize(0.08);
 273   dummy3.GetYaxis()->SetLabelSize(0.08);
 274   dummy3.GetXaxis()->SetTitleSize(0.08);
 275   dummy3.GetYaxis()->SetTitleSize(0.08);
 276   dummy3.GetYaxis()->SetTitleOffset(0.8);
 277   dummy3.DrawCopy();
 278
 279   if (syst)
 280   {
 281     // for the ratio add in quadrature
 282     for (Int_t bin=1; bin<=syst->GetNbinsX(); bin++)
 283     {
 284       if (corr->GetBinError(bin) > 0)
 285         syst->SetBinError(bin, TMath::Sqrt(TMath::Power(syst->GetBinContent(bin) / 100, 2) + TMath::Power(corr->GetBinError(bin) / corr->GetBinContent(bin), 2)));
 286       else
 287         syst->SetBinError(bin, 0);
 288       syst->SetBinContent(bin, 1);
 289     }
 290
 291     syst->Draw("E2 ][ SAME");
 292     dummy3.DrawCopy("SAME");
 293   }
 294
 295   ratio->Draw("SAME");
 296
 297   ratio->Fit("pol0", "N");
 298
 299   canvas3->Modified();
 300
 301 }
 302
 303 void loadlibs()
 304 {
 305   gSystem->Load("libANALYSIS");
 306   gSystem->Load("libANALYSISalice");
 307   gSystem->Load("libCORRFW");
 308   gSystem->Load("libJETAN");
 309   gSystem->Load("libPWG4JetTasks");
 310 }
 311
 312 void CompareBias(const char* mcFile = "PWG4_JetTasksOutput.root", Int_t region, Int_t ueHist)
 313 {
 314   loadlibs();
 315
 316   TFile::Open(mcFile);
 317   list = (TList*) gFile->Get(objectName);
 318   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 319   SetupRanges(h);
 320
 321   const char* axis = "z";
 322   Float_t ptLeadMin = 0;
 323   Float_t ptLeadMax = -1;
 324
 325   if (ueHist == 2)
 326   {
 327     ptLeadMin = 0.51 + 0;
 328     ptLeadMax = 1.99 + 0;
 329     axis = "y";
 330   }
 331
 332   biasFromData = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepBiasStudy, AliUEHist::kCFStepReconstructed, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromData");
 333   biasFromData2 = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepBiasStudy2, AliUEHist::kCFStepReconstructed, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromData2");
 334   //biasFromData = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, axis)->Clone("biasFromData");
 335   biasFromMC   = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepTracked, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromMC");
 336   //biasFromMC   = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, axis)->Clone("biasFromMC");
 337
 338 /*  biasFromData->Draw();
 339   biasFromMC->SetLineColor(2);
 340   biasFromMC->Draw("SAME");
 341   return;*/
 342
 343   DrawRatio(biasFromData, biasFromMC, "bias: data vs MC");
 344   DrawRatio(biasFromData, biasFromData2, "bias: data vs data two step");
 345 }
 346
 347 void CompareBiasWithData(const char* mcFile = "PWG4_JetTasksOutput.root", const char* dataFile = "esd.root")
 348 {
 349   loadlibs();
 350
 351   file1 = TFile::Open(mcFile);
 352   list = (TList*) file1->Get(objectName);
 353   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 354   SetupRanges(h);
 355
 356   biasFromMC   = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepTracked, "z")->Clone("biasFromMC");
 357
 358   file2 = TFile::Open(dataFile);
 359   list = (TList*) file2->Get(objectName);
 360   AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 361   SetupRanges(h2);
 362
 363   biasFromData = (TH1*) h2->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepBiasStudy, AliUEHist::kCFStepReconstructed, "z")->Clone("biasFromData");
 364   biasFromData2 = (TH1*) h2->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepBiasStudy2, AliUEHist::kCFStepReconstructed, "z")->Clone("biasFromData2");
 365
 366   DrawRatio(biasFromData, biasFromMC, "bias: data vs MC");
 367   DrawRatio(biasFromData, biasFromData2, "bias: data vs data two step");
 368 }
 369
 370 void Compare(const char* fileName1, const char* fileName2, Int_t id, Int_t step1, Int_t step2, Int_t region, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1)
 371 {
 372   loadlibs();
 373
 374   if (!gCache || !gFirst)
 375   {
 376     file1 = TFile::Open(fileName1);
 377     list = (TList*) file1->Get(objectName);
 378     AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 379     SetupRanges(h);
 380
 381     file2 = TFile::Open(fileName2);
 382     list = (TList*) file2->Get(objectName);
 383     AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 384     SetupRanges(h2);
 385   }
 386
 387   if (gCache)
 388   {
 389     if (!gFirst)
 390     {
 391       gFirst = h;
 392       gSecond = h2;
 393     }
 394     else
 395     {
 396       AliUEHistograms* h = (AliUEHistograms*) gFirst;
 397       AliUEHistograms* h2 = (AliUEHistograms*) gSecond;
 398     }
 399   }
 400
 401   TH1* hist1 = h->GetUEHist(id)->GetUEHist(step1, region, ptLeadMin, ptLeadMax);
 402   TH1* hist2 = h2->GetUEHist(id)->GetUEHist(step2, region, ptLeadMin, ptLeadMax);
 403
 404   TH1* trackHist = 0;
 405   if (id == 2)
 406   {
 407     trackHist = h->GetUEHist(id)->GetTrackHist(region)->ShowProjection(2, 0);
 408     // only keep bins under consideration
 409     for (Int_t bin=1; bin<=trackHist->GetNbinsX(); bin++)
 410       if (bin < trackHist->FindBin(ptLeadMin) || bin > trackHist->FindBin(ptLeadMax))
 411         trackHist->SetBinContent(bin, 0);
 412   }
 413
 414   // systematic uncertainty
 415   TH1* syst = GetSystematicUncertainty(hist1, trackHist);
 416
 417   DrawRatio(hist1, hist2, Form("%s_%d_%d_%d_%.2f_%.2f", TString(fileName1).Tokenize(".")->First()->GetName(), id, step1, region, ptLeadMin, ptLeadMax), syst);
 418 }
 419
 420 void CompareEventHist(const char* fileName1, const char* fileName2, Int_t id, Int_t step, Int_t var)
 421 {
 422   loadlibs();
 423
 424   file1 = TFile::Open(fileName1);
 425   list = (TList*) file1->Get(objectName);
 426   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 427   SetupRanges(h);
 428
 429   file2 = TFile::Open(fileName2);
 430   list = (TList*) file2->Get(objectName);
 431   AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 432   SetupRanges(h2);
 433
 434   TH1* hist1 = h->GetUEHist(id)->GetEventHist()->ShowProjection(var, step);
 435   TH1* hist2 = h2->GetUEHist(id)->GetEventHist()->ShowProjection(var, step);
 436
 437   DrawRatio(hist1, hist2, "compare");
 438 }
 439
 440 void CompareStep(const char* fileName1, const char* fileName2, Int_t id, Int_t step, Int_t region, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1)
 441 {
 442   // fileName1 is labelled Corrected in the plot
 443
 444   loadlibs();
 445
 446   gUEHist = id;
 447   gRegion = region;
 448   Compare(fileName1, fileName2, id, step, step, region, ptLeadMin, ptLeadMax);
 449 }
 450
 451 void DrawStep(const char* fileName, Int_t id, Int_t step, Int_t region, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1)
 452 {
 453   loadlibs();
 454
 455   file1 = TFile::Open(fileName);
 456   list = (TList*) file1->Get(objectName);
 457   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 458   SetupRanges(h);
 459
 460   new TCanvas;
 461   h->GetUEHist(id)->GetUEHist(step, region, ptLeadMin, ptLeadMax)->Draw();
 462 }
 463
 464 void ProfileMultiplicity(const char* fileName = "PWG4_JetTasksOutput.root")
 465 {
 466   loadlibs();
 467
 468   TFile::Open(fileName);
 469   list = (TList*) gFile->Get(objectName);
 470   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 471
 472   new TCanvas;
 473   h->GetCorrelationMultiplicity()->Draw("colz");
 474   gPad->SetLogz();
 475
 476   new TCanvas;
 477   h->GetCorrelationMultiplicity()->ProfileX()->DrawCopy()->Fit("pol1", "", "", 1, 10);
 478 }
 479
 480 void SetupRanges(void* obj)
 481 {
 482   ((AliUEHistograms*) obj)->SetEtaRange(-0.79, 0.79);
 483   ((AliUEHistograms*) obj)->SetPtRange(gLeadingpTMin, 100);
 484   ((AliUEHistograms*) obj)->SetCombineMinMax(kTRUE);
 485 }
 486
 487 void DrawRatios(const char* name, void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1)
 488 {
 489   AliUEHist* corrected = (AliUEHist*) correctedVoid;
 490   AliUEHist* comparison = (AliUEHist*) comparisonVoid;
 491
 492   Int_t beginStep = AliUEHist::kCFStepAll;
 493   Int_t endStep = AliUEHist::kCFStepReconstructed;
 494
 495   if (compareStep != -1)
 496   {
 497     beginStep = compareStep;
 498     endStep = compareStep;
 499   }
 500
 501   Int_t beginRegion = 0;
 502   Int_t endRegion = 2;
 503
 504   if (compareRegion != -1)
 505   {
 506     beginRegion = compareRegion;
 507     endRegion = compareRegion;
 508   }
 509
 510   for (Int_t step=beginStep; step<=endStep; step++)
 511   {
 512     for (Int_t region=beginRegion; region <= endRegion; region++)
 513     {
 514       Printf("%f %f", ptLeadMin, ptLeadMax);
 515       TH1* corrHist = corrected->GetUEHist(step, region, ptLeadMin, ptLeadMax);
 516       TH1* mcHist   = comparison->GetUEHist(step, region, ptLeadMin, ptLeadMax);
 517
 518       DrawRatio(corrHist, mcHist, TString(Form("%s: step %d %s %s", name, step, corrected->GetStepTitle(step), corrected->GetRegionTitle(region))));
 519     }
 520   }
 521 }
 522
 523 void DrawRatios(void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 524 {
 525   AliUEHistograms* corrected = (AliUEHistograms*) correctedVoid;
 526   AliUEHistograms* comparison = (AliUEHistograms*) comparisonVoid;
 527
 528   if (1 && compareUEHist == 2)
 529   {
 530     for (Float_t ptLeadMin = 0.51; ptLeadMin < 10; ptLeadMin += 1.5)
 531       DrawRatios(TString(Form("UE %d pT %f", compareUEHist, ptLeadMin)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion, ptLeadMin, ptLeadMin + 0.48);
 532     return;
 533   }
 534
 535   if (compareUEHist == -1)
 536   {
 537     for (Int_t i=0; i<2; i++)
 538       DrawRatios(TString(Form("UE %d", i)), corrected->GetUEHist(i), comparison->GetUEHist(i), compareStep, compareRegion);
 539   }
 540   else
 541     DrawRatios(TString(Form("UE %d", compareUEHist)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion);
 542 }
 543
 544 void CompareEventsTracks(void* corrVoid, void* mcVoid, Int_t compareStep, Int_t compareRegion, Int_t compareUEHist = 0)
 545 {
 546   AliUEHist* corr = ((AliUEHistograms*) corrVoid)->GetUEHist(compareUEHist);
 547   AliUEHist* mc   = ((AliUEHistograms*) mcVoid)->GetUEHist(compareUEHist);
 548
 549   Float_t ptLeadMin = 0;
 550   Float_t ptLeadMax = -1;
 551   Int_t axis = 2;
 552
 553   if (compareUEHist == 2)
 554   {
 555     ptLeadMin = 1.01;
 556     ptLeadMax = 1.49;
 557     axis = 4;
 558   }
 559
 560   TH1* corrHist = corr->GetUEHist(compareStep, compareRegion, ptLeadMin, ptLeadMax);
 561   TH1* mcHist   = mc  ->GetUEHist(compareStep, compareRegion, ptLeadMin, ptLeadMax);
 562
 563   DrawRatio(corrHist, mcHist, Form("check"));
 564
 565   corr->SetBinLimits(corr->GetTrackHist(compareRegion)->GetGrid(compareStep));
 566   mc->SetBinLimits(corr->GetTrackHist(compareRegion)->GetGrid(compareStep));
 567
 568   corrHist =  corr->GetTrackHist(compareRegion)->GetGrid(compareStep)->Project(axis);
 569   mcHist   =  mc  ->GetTrackHist(compareRegion)->GetGrid(compareStep)->Project(axis);
 570   DrawRatio(corrHist, mcHist, Form("check2"));
 571
 572   corrHist =  corr->GetEventHist()->GetGrid(compareStep)->Project(0);
 573   mcHist   =  mc  ->GetEventHist()->GetGrid(compareStep)->Project(0);
 574   DrawRatio(corrHist, mcHist, Form("check3"));
 575 }
 576
 577 void correctMC(const char* fileNameCorrections, const char* fileNameESD = 0, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 578 {
 579   // corrects the reconstructed step in fileNameESD with fileNameCorr
 580   // if fileNameESD is 0 data from fileNameCorr is taken
 581   // afterwards the corrected distributions are compared with the MC stored in fileNameESD
 582
 583   loadlibs();
 584
 585   TFile::Open(fileNameCorrections);
 586   list = (TList*) gFile->Get(objectName);
 587   AliUEHistograms* corr = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 588   SetupRanges(corr);
 589   corr->ExtendTrackingEfficiency();
 590
 591   if (fileNameESD)
 592   {
 593     file2 = TFile::Open(fileNameESD);
 594     list = (TList*) file2->Get(objectName);
 595   }
 596   AliUEHistograms* testSample = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 597
 598   // copy to esd object
 599   AliUEHistograms* esd = new AliUEHistograms;
 600   esd->CopyReconstructedData(testSample);
 601
 602   SetupRanges(corr);
 603   SetupRanges(testSample);
 604   SetupRanges(esd);
 605
 606   esd->Correct(corr);
 607   esd->GetUEHist(2)->AdditionalDPhiCorrection(0);
 608
 609   DrawRatios(esd, testSample, compareStep, compareRegion, compareUEHist);
 610
 611   //CompareEventsTracks(esd, testSample, compareStep, compareRegion, compareUEHist);
 612 }
 613
 614 // function to compare only final step for all regions and distributions
 615
 616 void correctData(const char* fileNameCorrections, const char* fileNameESD, const char* contEnhancement, Float_t contEncUpTo = 1.0, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 617 {
 618   // corrects fileNameESD with fileNameCorrections and compares the two
 619
 620   loadlibs();
 621
 622   TFile::Open(fileNameCorrections);
 623   list = (TList*) gFile->Get(objectName);
 624   AliUEHistograms* corr = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 625
 626   file2 = TFile::Open(fileNameESD);
 627   list = (TList*) file2->Get(objectName);
 628   AliUEHistograms* esd = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 629
 630   SetupRanges(corr);
 631   SetupRanges(esd);
 632
 633   corr->ExtendTrackingEfficiency();
 634
 635   TFile::Open(contEnhancement);
 636   contEncHist = (TH1*) gFile->Get("histo");
 637   contEncHistFullRange = (TH1*) corr->GetUEHist(0)->GetTrackingEfficiency(1)->Clone("contEncHistFullRange");
 638
 639   contEncHistFullRange->Reset();
 640   for (Int_t i=1; i<=contEncHistFullRange->GetNbinsX(); i++)
 641   {
 642     contEncHistFullRange->SetBinContent(i, 1);
 643     if (i <= contEncHist->GetNbinsX() && contEncHist->GetXaxis()->GetBinCenter(i) < contEncUpTo && contEncHist->GetBinContent(i) > 0)
 644       contEncHistFullRange->SetBinContent(i, contEncHist->GetBinContent(i));
 645   }
 646   corr->SetContaminationEnhancement((TH1F*) contEncHistFullRange);
 647
 648   esd->Correct(corr);
 649   esd->GetUEHist(2)->AdditionalDPhiCorrection(0);
 650
 651   file3 = TFile::Open("corrected.root", "RECREATE");
 652   file3->mkdir("PWG4_LeadingTrackUE");
 653   file3->cd("PWG4_LeadingTrackUE");
 654   list->Write(0, TObject::kSingleKey);
 655   file3->Close();
 656
 657   DrawRatios(esd, corr, compareStep, compareRegion, compareUEHist);
 658 }
 659
 660 void ITSTPCEfficiency(const char* fileNameData, Int_t id, Int_t itsTPC = 0)
 661 {
 662   // its = 0; tpc = 1
 663
 664   // uncertainty from dN/dpT paper
 665   Double_t pTBins[] =  {0.0, 0.1, 0.15,  0.2,  0.25,  0.3,   0.35,  0.4,   0.45,  0.5,   0.6,   0.7,   0.8,   0.9,   1.0,   1.5,   2.0,   2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
 666   Float_t effITS[] =   {0.,  0.,  0.995, 0.98, 0.986, 0.996, 1,     1,     1,     1,     1,     1,     1,     1,     1,     1,     1,    };  // the last three are the same because i don't have entries
 667   Float_t effTPC[] =   {0.,  0,   1.042, 1.026,1.021, 1.018, 1.015, 1.015, 1.012, 1.012, 1.007, 1.0075,1.006, 1.006, 1.004, 1.004, 1.009 }; // the last bins put as if they were the same
 668
 669   TH1F* effHist = new TH1F("effHist", "effHist", 39, pTBins);
 670   for (Int_t i=0; i<39; i++)
 671   {
 672     Int_t bin = i;
 673     if (i > 16)
 674       bin = 16;
 675     effHist->SetBinContent(i+1, (itsTPC == 0) ? effITS[bin] : effTPC[bin]);
 676   }
 677
 678   new TCanvas; effHist->Draw();
 679
 680   EffectOfModifiedTrackingEfficiency(fileNameData, id, 2, effHist, 1, -1, (itsTPC == 0) ? "ITS" : "TPC");
 681 }
 682
 683
 684 void EffectOfModifiedTrackingEfficiency(const char* fileNameData, Int_t id, Int_t region, TH1* trackingEff, Int_t axis1, Int_t axis2 = -1, const char* name = "EffectOfModifiedTrackingEfficiency")
 685 {
 686   // trackingEff should contain the change in tracking efficiency, i.e. between before and after in the eta-pT plane
 687
 688   loadlibs();
 689
 690   TFile::Open(fileNameData);
 691   list = (TList*) gFile->Get(objectName);
 692   AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 693   SetupRanges(corrected);
 694
 695   AliUEHist* ueHist = corrected->GetUEHist(id);
 696
 697   Float_t ptLeadMin = -1;
 698   Float_t ptLeadMax = -1;
 699   if (id == 2)
 700   {
 701     // the uncertainty is flat in delta phi, so use this trick to get directly the uncertainty as function of leading pT
 702     //ptLeadMin = 1.01;
 703     //ptLeadMax = 1.99;
 704   }
 705
 706   // histogram before
 707   TH1* before = ueHist->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
 708
 709   // copy histogram
 710   // the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
 711   ueHist->CorrectTracks(AliUEHist::kCFStepAll, AliUEHist::kCFStepTriggered, (TH1*) 0, 0, -1);
 712
 713   // reapply tracking efficiency
 714   ueHist->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, trackingEff, axis1, axis2);
 715
 716   // histogram after
 717   TH1* after = ueHist->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
 718
 719   DrawRatio(before, after, name);
 720   gPad->GetCanvas()->SaveAs(Form("%s.png", name));
 721 }
 722
 723 void EffectOfTrackCuts(const char* fileNameData, Int_t id, const char* systFile)
 724 {
 725   loadlibs();
 726
 727   TFile::Open(fileNameData);
 728   list = (TList*) gFile->Get(objectName);
 729   AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 730   effHist = (TH2D*) corrected->GetUEHist(0)->GetTrackingEfficiency()->Clone("effHist");
 731
 732   file = TFile::Open(systFile);
 733
 734   Int_t maxSyst = 3;
 735   const char* systNames[] = { "NClusTPC", "Chi2TPC", "SigmaDCA" };
 736
 737   for (Int_t i=0; i<maxSyst; i++)
 738   {
 739     for (Int_t j=0; j<2; j++)
 740     {
 741       TString histName;
 742       histName.Form("%s_syst_%s", systNames[i], (j == 0) ? "up" : "down");
 743       systEffect = (TH2*) file->Get(histName);
 744
 745       // rebin
 746       effHist->Reset();
 747       for (Int_t x=1; x <= effHist->GetNbinsX(); x++)
 748         for (Int_t y=1; y <= effHist->GetNbinsY(); y++)
 749           effHist->SetBinContent(x, y, 1);
 750
 751       for (Int_t x=1; x <= systEffect->GetNbinsX(); x++)
 752         for (Int_t y=1; y <= systEffect->GetNbinsY(); y++)
 753           if (systEffect->GetBinContent(x, y) != 0)
 754             effHist->SetBinContent(effHist->GetXaxis()->FindBin(systEffect->GetYaxis()->GetBinCenter(y)), effHist->GetYaxis()->FindBin(systEffect->GetXaxis()->GetBinCenter(x)), systEffect->GetBinContent(x, y));
 755
 756
 757       //new TCanvas; systEffect->Draw("COLZ"); new TCanvas; effHist->Draw("COLZ");
 758
 759       EffectOfModifiedTrackingEfficiency(fileNameData, id, (id == 2) ? 0 : 2, effHist, 0, 1, histName);
 760
 761       //return;
 762     }
 763   }
 764 }
 765
 766 void ModifyComposition(const char* fileNameData, const char* fileNameCorrections, Int_t id, Bool_t verbose = kFALSE)
 767 {
 768   loadlibs();
 769
 770   TFile::Open(fileNameData);
 771   list = (TList*) gFile->Get(objectName);
 772   AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 773   SetupRanges(corrected);
 774
 775   file2 = TFile::Open(fileNameCorrections);
 776   list2 = (TList*) file2->Get(objectName);
 777   AliUEHistograms* corrections = (AliUEHistograms*) list2->FindObject("AliUEHistograms");
 778   SetupRanges(corrections);
 779
 780   ueHistData        = (AliUEHist*) corrected->GetUEHist(id);
 781   ueHistCorrections = (AliUEHist*) corrections->GetUEHist(id);
 782
 783   // copy histogram
 784   // the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
 785   ueHistData->CorrectTracks(AliUEHist::kCFStepAll, AliUEHist::kCFStepTriggered, (TH1*) 0, 0);
 786
 787   Int_t maxRegion = 3;
 788   Float_t ptLeadMin = -1;
 789   Float_t ptLeadMax = -1;
 790   if (id == 2)
 791   {
 792     maxRegion = 1;
 793     // the uncertainty is flat in delta phi, so use this trick to get directly the uncertainty as function of leading pT
 794     //ptLeadMin = 1.01;
 795     //ptLeadMax = 1.99;
 796   }
 797
 798   // histogram before
 799   TH1* before[3];
 800   for (Int_t region=0; region<maxRegion; region++)
 801     before[region] = ueHistData->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
 802
 803   //defaultEff = ueHistCorrections->GetTrackingEfficiency();
 804   defaultEff = ueHistCorrections->GetTrackingCorrection();
 805   //defaultEffpT = ueHistCorrections->GetTrackingEfficiency(1);
 806   defaultEffpT = ueHistCorrections->GetTrackingCorrection(1);
 807   defaultContainer = ueHistCorrections->GetTrackHistEfficiency();
 808
 809   c = new TCanvas;
 810   defaultEffpT->Draw("");
 811
 812   Float_t largestDeviation[3];
 813   for (Int_t i=0; i<maxRegion; i++)
 814     largestDeviation[i] = 0;
 815
 816   for (Int_t i=0; i<7; i++)
 817   {
 818     // case 0: // default
 819     // case 1: // + 30% kaons
 820     // case 2: // - 30% kaons
 821     // case 3: // + 30% protons
 822     // case 4: // - 30% protons
 823     // case 5: // + 30% others
 824     // case 6: // - 30% others
 825     Int_t correctionIndex = (i+1) / 2 + 1; // bin 1 == protons, bin 2 == kaons, ...
 826     Double_t scaleFactor = (i % 2 == 1) ? 1.3 : 0.7;
 827     if (i == 0)
 828       scaleFactor = 1;
 829
 830     newContainer = (AliCFContainer*) defaultContainer->Clone();
 831
 832     // modify, change all steps
 833     for (Int_t j=0; j<newContainer->GetNStep(); j++)
 834     {
 835       THnSparse* grid = newContainer->GetGrid(j)->GetGrid();
 836
 837       for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 838       {
 839         Int_t bins[5];
 840         Double_t value = grid->GetBinContent(binIdx, bins);
 841         Double_t error = grid->GetBinError(binIdx);
 842
 843         if (bins[2] != correctionIndex)
 844           continue;
 845
 846         value *= scaleFactor;
 847         error *= scaleFactor;
 848
 849         grid->SetBinContent(bins, value);
 850         grid->SetBinError(bins, error);
 851       }
 852     }
 853
 854     // put in corrections
 855     ueHistCorrections->SetTrackHistEfficiency(newContainer);
 856
 857     // ratio
 858     //modifiedEff = ueHistCorrections->GetTrackingEfficiency();
 859     modifiedEff = ueHistCorrections->GetTrackingCorrection();
 860     modifiedEff->Divide(modifiedEff, defaultEff);
 861     //modifiedEff->Draw("COLZ");
 862
 863     c->cd();
 864     //modifiedEffpT = ueHistCorrections->GetTrackingEfficiency(1);
 865     modifiedEffpT = ueHistCorrections->GetTrackingCorrection(1);
 866     modifiedEffpT->SetLineColor(i+1);
 867     modifiedEffpT->Draw("SAME");
 868
 869     // apply change in tracking efficiency
 870     ueHistData->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, modifiedEff, 0, 1);
 871
 872     for (Int_t region=0; region<maxRegion; region++)
 873     {
 874       // histogram after
 875       TH1* after = ueHistData->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
 876
 877       if (verbose)
 878         DrawRatio(before[region], (TH1*) after->Clone(), Form("Region %d Composition %d", region, i));
 879
 880       // ratio is flat, extract deviation
 881       after->Divide(before[region]);
 882       after->Fit("pol0", "0");
 883       Float_t deviation = 100.0 - 100.0 * after->GetFunction("pol0")->GetParameter(0);
 884       Printf("Deviation for region %d case %d is %.2f %%", region, i, deviation);
 885
 886       if (TMath::Abs(deviation) > largestDeviation[region])
 887         largestDeviation[region] = TMath::Abs(deviation);
 888     }
 889     //return;
 890   }
 891
 892   for (Int_t i=0; i<maxRegion; i++)
 893     Printf("Largest deviation in region %d is %f", i, largestDeviation[i]);
 894 }
 895
 896 void MergeList()
 897 {
 898   loadlibs();
 899
 900   ifstream in;
 901   in.open("list");
 902
 903   TFileMerger m;
 904
 905   TString line;
 906   while (in.good())
 907   {
 908     in >> line;
 909
 910     if (line.Length() == 0)
 911       continue;
 912
 913     TString fileName;
 914     fileName.Form("%s/%s/PWG4_JetTasksOutput.root", "maps", line.Data());
 915     Printf("%s", fileName.Data());
 916
 917     m.AddFile(fileName);
 918   }
 919
 920   m.SetFastMethod();
 921   m.OutputFile("merged.root");
 922   m.Merge();
 923 }
 924
 925 void MergeList2(const char* listFile, Bool_t onlyPrintEvents = kFALSE)
 926 {
 927   loadlibs();
 928
 929   ifstream in;
 930   in.open(listFile);
 931
 932   AliUEHistograms* final = 0;
 933   TList* finalList = 0;
 934
 935   TString line;
 936   while (in.good())
 937   {
 938     in >> line;
 939
 940     if (line.Length() == 0)
 941       continue;
 942
 943     TString fileName;
 944     fileName.Form("%s/%s/PWG4_JetTasksOutput.root", "maps", line.Data());
 945     Printf("%s", fileName.Data());
 946
 947     file = TFile::Open(fileName);
 948     if (!file)
 949       continue;
 950
 951     list = (TList*) file->Get(objectName);
 952
 953     if (!final)
 954     {
 955       final = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 956       //final->GetEventCount()->Draw(); return;
 957       Printf("Events: %d", (Int_t) final->GetEventCount()->ProjectionX()->GetBinContent(4));
 958       finalList = list;
 959     }
 960     else
 961     {
 962       additional = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 963       Printf("Events: %d", (Int_t) additional->GetEventCount()->ProjectionX()->GetBinContent(4));
 964
 965       if (!onlyPrintEvents)
 966       {
 967         TList list2;
 968         list2.Add(additional);
 969         final->Merge(&list2);
 970       }
 971       delete additional;
 972       file->Close();
 973     }
 974   }
 975
 976   if (onlyPrintEvents)
 977     return;
 978
 979   Printf("Total events (at step 0): %d", (Int_t) final->GetEventCount()->ProjectionX()->GetBinContent(4));
 980
 981   file3 = TFile::Open("merged.root", "RECREATE");
 982   file3->mkdir("PWG4_LeadingTrackUE");
 983   file3->cd("PWG4_LeadingTrackUE");
 984   finalList->Write(0, TObject::kSingleKey);
 985   file3->Close();
 986 }
 987
 988 void PlotAll(const char* correctedFile, const char* mcFile)
 989 {
 990   gCache = 1;
 991
 992   if (gEnergy == 900)
 993   {
 994     Float_t range[] = { 1.5, 2 };
 995   }
 996   else
 997   {
 998     Float_t range[] = { 3, 10 };
 999   }
1000
1001   for (Int_t id=0; id<3; id++)
1002   {
1003     if (id < 2)
1004       gForceRange = range[id];
1005     else
1006       gForceRange = -1;
1007
1008     if (id < 2)
1009     {
1010       for (Int_t region=0; region<3; region++)
1011       {
1012         CompareStep(correctedFile, mcFile, id, 0, region);
1013         gPad->GetCanvas()->SaveAs(Form("%s_%s_%d_%d.png", TString(correctedFile).Tokenize(".")->First()->GetName(), TString(mcFile).Tokenize(".")->First()->GetName(), id, region));
1014       }
1015     }
1016     else
1017     {
1018       Float_t leadingPtArr[] = { 0.50, 2.0, 4.0, 6.0, 10.0, 20.0, 50.0 };
1019       for (Int_t leadingPtID=0; leadingPtID<6; leadingPtID++)
1020       {
1021         CompareStep(correctedFile, mcFile, id, 0, 0, leadingPtArr[leadingPtID] + 0.01, leadingPtArr[leadingPtID+1] - 0.01);
1022         gPad->GetCanvas()->SaveAs(Form("%s_%s_%d_%.2f_%.2f.png", TString(correctedFile).Tokenize(".")->First()->GetName(), TString(mcFile).Tokenize(".")->First()->GetName(), id, leadingPtArr[leadingPtID], leadingPtArr[leadingPtID+1]));
1023       }
1024     }
1025   }
1026 }
1027