PWG4/macros/underlyingevent/correct.C

   1 const char* objectName = "PWG4_LeadingTrackUE/histosLeadingTrackUE";
   2 Float_t gLeadingpTMin = 0.51;
   3 //const char* objectName = "PWG4_LeadingTrackUE/histosLeadingTrackUE_filterbit32";
   4
   5 void SetRanges(TAxis* axis)
   6 {
   7   if (strcmp(axis->GetTitle(), "leading p_{T} (GeV/c)") == 0)
   8     axis->SetRangeUser(0, 10);
   9
  10   if (strcmp(axis->GetTitle(), "multiplicity") == 0)
  11     axis->SetRangeUser(0, 10);
  12 }
  13
  14 void SetRanges(TH1* hist)
  15 {
  16   SetRanges(hist->GetXaxis());
  17   SetRanges(hist->GetYaxis());
  18   SetRanges(hist->GetZaxis());
  19 }
  20
  21 void Prepare1DPlot(TH1* hist)
  22 {
  23   hist->SetLineWidth(2);
  24   hist->SetStats(kFALSE);
  25
  26   hist->GetXaxis()->SetLabelOffset(0.02);
  27   hist->GetXaxis()->SetTitleOffset(1.3);
  28   hist->GetYaxis()->SetTitleOffset(1.3);
  29
  30   SetRanges(hist);
  31 }
  32
  33 void DrawRatio(TH1* corr, TH1* mc, const char* name)
  34 {
  35   mc->SetLineColor(2);
  36
  37   TCanvas* canvas3 = new TCanvas(name, name, 600, 600);
  38   canvas3->Range(0, 0, 1, 1);
  39
  40   TPad* pad1 = new TPad(Form("%s_1", name), "", 0, 0.5, 0.98, 0.98);
  41   pad1->Draw();
  42
  43   TPad* pad2 = new TPad(Form("%s_2", name), "", 0, 0.02, 0.98, 0.5);
  44   pad2->Draw();
  45
  46   pad1->SetRightMargin(0.01);
  47   pad2->SetRightMargin(0.01);
  48   pad1->SetTopMargin(0.05);
  49   pad1->SetLeftMargin(0.13);
  50   pad2->SetLeftMargin(0.13);
  51   pad2->SetBottomMargin(0.22);
  52
  53   // no border between them
  54   pad1->SetBottomMargin(0);
  55   pad2->SetTopMargin(0);
  56
  57   pad1->cd();
  58   pad1->SetGridx();
  59   pad1->SetGridy();
  60
  61   TLegend* legend = new TLegend(0.15, 0.65, 0.55, 0.90);
  62   legend->SetFillColor(0);
  63   legend->AddEntry(corr, "Corrected");
  64   legend->AddEntry(mc, "MC prediction");
  65   legend->SetTextSize(0.08);
  66
  67   Prepare1DPlot(corr);
  68
  69   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);
  70   dummy->SetStats(kFALSE);
  71   dummy->SetXTitle(corr->GetXaxis()->GetTitle());
  72   dummy->SetYTitle(corr->GetYaxis()->GetTitle());
  73   dummy->GetYaxis()->SetTitleOffset(1);
  74   Prepare1DPlot(dummy);
  75
  76   dummy->GetXaxis()->SetLabelSize(0.08);
  77   dummy->GetYaxis()->SetLabelSize(0.08);
  78   dummy->GetXaxis()->SetTitleSize(0.08);
  79   dummy->GetYaxis()->SetTitleSize(0.08);
  80   dummy->GetYaxis()->SetTitleOffset(0.8);
  81   dummy->DrawCopy();
  82
  83   corr->Draw("SAME");
  84   mc->Draw("SAME");
  85
  86   legend->Draw();
  87
  88   pad2->cd();
  89   //pad2->SetBottomMargin(0.15);
  90   //pad2->SetGridx();
  91   //pad2->SetGridy();
  92
  93   TH1* ratio = (TH1*) mc->Clone("ratio");
  94   ratio->Divide(corr);
  95
  96   Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
  97   Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
  98
  99   minR = 0.8;
 100   maxR = 1.2;
 101
 102   TH1F dummy3("dummy3", ";;Ratio: MC / corr", 100, corr->GetXaxis()->GetBinLowEdge(1), corr->GetXaxis()->GetBinUpEdge(corr->GetNbinsX()));
 103   dummy3.SetXTitle(corr->GetXaxis()->GetTitle());
 104   Prepare1DPlot(&dummy3);
 105   dummy3.SetStats(kFALSE);
 106   for (Int_t i=1; i<=100; ++i)
 107         dummy3.SetBinContent(i, 1);
 108   dummy3.GetYaxis()->SetRangeUser(minR, maxR);
 109   dummy3.SetLineWidth(2);
 110   dummy3.GetXaxis()->SetLabelSize(0.08);
 111   dummy3.GetYaxis()->SetLabelSize(0.08);
 112   dummy3.GetXaxis()->SetTitleSize(0.08);
 113   dummy3.GetYaxis()->SetTitleSize(0.08);
 114   dummy3.GetYaxis()->SetTitleOffset(0.8);
 115   dummy3.DrawCopy();
 116
 117   ratio->Draw("SAME");
 118
 119   canvas3->Modified();
 120 }
 121
 122 void loadlibs()
 123 {
 124   gSystem->Load("libANALYSIS");
 125   gSystem->Load("libANALYSISalice");
 126   gSystem->Load("libCORRFW");
 127   gSystem->Load("libJETAN");
 128   gSystem->Load("libPWG4JetTasks");
 129 }
 130
 131 void CompareBias(const char* mcFile = "PWG4_JetTasksOutput.root", Int_t region, Int_t ueHist)
 132 {
 133   loadlibs();
 134
 135   TFile::Open(mcFile);
 136   list = (TList*) gFile->Get(objectName);
 137   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 138   SetupRanges(h);
 139
 140   const char* axis = "z";
 141   Float_t ptLeadMin = 0;
 142   Float_t ptLeadMax = -1;
 143
 144   if (ueHist == 2)
 145   {
 146     ptLeadMin = 0.51 + 0;
 147     ptLeadMax = 1.99 + 0;
 148     axis = "y";
 149   }
 150
 151   biasFromData = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepBiasStudy, AliUEHist::kCFStepReconstructed, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromData");
 152   biasFromData2 = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepBiasStudy2, AliUEHist::kCFStepReconstructed, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromData2");
 153   //biasFromData = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, axis)->Clone("biasFromData");
 154   biasFromMC   = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepTracked, region, axis, ptLeadMin, ptLeadMax)->Clone("biasFromMC");
 155   //biasFromMC   = (TH1*) h->GetUEHist(ueHist)->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepRealLeading, region, axis)->Clone("biasFromMC");
 156
 157 /*  biasFromData->Draw();
 158   biasFromMC->SetLineColor(2);
 159   biasFromMC->Draw("SAME");
 160   return;*/
 161
 162   DrawRatio(biasFromData, biasFromMC, "bias: data vs MC");
 163   DrawRatio(biasFromData, biasFromData2, "bias: data vs data two step");
 164 }
 165
 166 void CompareBiasWithData(const char* mcFile = "PWG4_JetTasksOutput.root", const char* dataFile = "esd.root")
 167 {
 168   loadlibs();
 169
 170   file1 = TFile::Open(mcFile);
 171   list = (TList*) file1->Get(objectName);
 172   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 173   SetupRanges(h);
 174
 175   biasFromMC   = (TH1*) h->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepReconstructed, AliUEHist::kCFStepTracked, "z")->Clone("biasFromMC");
 176
 177   file2 = TFile::Open(dataFile);
 178   list = (TList*) file2->Get(objectName);
 179   AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 180   SetupRanges(h2);
 181
 182   biasFromData = (TH1*) h2->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepBiasStudy, AliUEHist::kCFStepReconstructed, "z")->Clone("biasFromData");
 183   biasFromData2 = (TH1*) h2->GetNumberDensitypT()->GetBias(AliUEHist::kCFStepBiasStudy2, AliUEHist::kCFStepReconstructed, "z")->Clone("biasFromData2");
 184
 185   DrawRatio(biasFromData, biasFromMC, "bias: data vs MC");
 186   DrawRatio(biasFromData, biasFromData2, "bias: data vs data two step");
 187 }
 188
 189 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)
 190 {
 191   loadlibs();
 192
 193   file1 = TFile::Open(fileName1);
 194   list = (TList*) file1->Get(objectName);
 195   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 196   SetupRanges(h);
 197
 198   file2 = TFile::Open(fileName2);
 199   list = (TList*) file2->Get(objectName);
 200   AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 201   SetupRanges(h2);
 202
 203   TH1* hist1 = h->GetUEHist(id)->GetUEHist(step1, region, ptLeadMin, ptLeadMax);
 204   TH1* hist2 = h2->GetUEHist(id)->GetUEHist(step2, region, ptLeadMin, ptLeadMax);
 205
 206   DrawRatio(hist1, hist2, "compare");
 207 }
 208
 209 void CompareEventHist(const char* fileName1, const char* fileName2, Int_t id, Int_t step, Int_t var)
 210 {
 211   loadlibs();
 212
 213   file1 = TFile::Open(fileName1);
 214   list = (TList*) file1->Get(objectName);
 215   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 216   SetupRanges(h);
 217
 218   file2 = TFile::Open(fileName2);
 219   list = (TList*) file2->Get(objectName);
 220   AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 221   SetupRanges(h2);
 222
 223   TH1* hist1 = h->GetUEHist(id)->GetEventHist()->ShowProjection(var, step);
 224   TH1* hist2 = h2->GetUEHist(id)->GetEventHist()->ShowProjection(var, step);
 225
 226   DrawRatio(hist1, hist2, "compare");
 227 }
 228
 229 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)
 230 {
 231   loadlibs();
 232
 233   Compare(fileName1, fileName2, id, step, step, region, ptLeadMin, ptLeadMax);
 234 }
 235
 236 void DrawStep(const char* fileName, Int_t id, Int_t step, Int_t region, Float_t ptLeadMin = -1, Float_t ptLeadMax = -1)
 237 {
 238   loadlibs();
 239
 240   file1 = TFile::Open(fileName);
 241   list = (TList*) file1->Get(objectName);
 242   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 243   SetupRanges(h);
 244
 245   new TCanvas;
 246   h->GetUEHist(id)->GetUEHist(step, region, ptLeadMin, ptLeadMax)->Draw();
 247 }
 248
 249 void ProfileMultiplicity(const char* fileName = "PWG4_JetTasksOutput.root")
 250 {
 251   loadlibs();
 252
 253   TFile::Open(fileName);
 254   list = (TList*) gFile->Get(objectName);
 255   AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 256
 257   new TCanvas;
 258   h->GetCorrelationMultiplicity()->Draw("colz");
 259   gPad->SetLogz();
 260
 261   new TCanvas;
 262   h->GetCorrelationMultiplicity()->ProfileX()->DrawCopy()->Fit("pol1", "", "", 1, 10);
 263 }
 264
 265 void SetupRanges(void* obj)
 266 {
 267   ((AliUEHistograms*) obj)->SetEtaRange(-0.79, 0.79);
 268   ((AliUEHistograms*) obj)->SetPtRange(gLeadingpTMin, 100);
 269   ((AliUEHistograms*) obj)->SetCombineMinMax(kTRUE);
 270 }
 271
 272 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)
 273 {
 274   AliUEHist* corrected = (AliUEHist*) correctedVoid;
 275   AliUEHist* comparison = (AliUEHist*) comparisonVoid;
 276
 277   Int_t beginStep = AliUEHist::kCFStepAll;
 278   Int_t endStep = AliUEHist::kCFStepReconstructed;
 279
 280   if (compareStep != -1)
 281   {
 282     beginStep = compareStep;
 283     endStep = compareStep;
 284   }
 285
 286   Int_t beginRegion = 0;
 287   Int_t endRegion = 2;
 288
 289   if (compareRegion != -1)
 290   {
 291     beginRegion = compareRegion;
 292     endRegion = compareRegion;
 293   }
 294
 295   for (Int_t step=beginStep; step<=endStep; step++)
 296   {
 297     for (Int_t region=beginRegion; region <= endRegion; region++)
 298     {
 299       Printf("%f %f", ptLeadMin, ptLeadMax);
 300       TH1* corrHist = corrected->GetUEHist(step, region, ptLeadMin, ptLeadMax);
 301       TH1* mcHist   = comparison->GetUEHist(step, region, ptLeadMin, ptLeadMax);
 302
 303       DrawRatio(corrHist, mcHist, TString(Form("%s: step %d %s %s", name, step, corrected->GetStepTitle(step), corrected->GetRegionTitle(region))));
 304     }
 305   }
 306 }
 307
 308 void DrawRatios(void* correctedVoid, void* comparisonVoid, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 309 {
 310   AliUEHistograms* corrected = (AliUEHistograms*) correctedVoid;
 311   AliUEHistograms* comparison = (AliUEHistograms*) comparisonVoid;
 312
 313   if (compareUEHist == 2)
 314   {
 315     for (Float_t ptLeadMin = 1.01; ptLeadMin < 4; ptLeadMin += 2)
 316       DrawRatios(TString(Form("UE %d pT %f", compareUEHist, ptLeadMin)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion, ptLeadMin, ptLeadMin + 0.48);
 317     return;
 318   }
 319
 320   if (compareUEHist == -1)
 321   {
 322     for (Int_t i=0; i<2; i++)
 323       DrawRatios(TString(Form("UE %d", i)), corrected->GetUEHist(i), comparison->GetUEHist(i), compareStep, compareRegion);
 324   }
 325   else
 326     DrawRatios(TString(Form("UE %d", compareUEHist)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion);
 327 }
 328
 329 void CompareEventsTracks(void* corrVoid, void* mcVoid, Int_t compareStep, Int_t compareRegion, Int_t compareUEHist = 0)
 330 {
 331   AliUEHist* corr = ((AliUEHistograms*) corrVoid)->GetUEHist(compareUEHist);
 332   AliUEHist* mc   = ((AliUEHistograms*) mcVoid)->GetUEHist(compareUEHist);
 333
 334   Float_t ptLeadMin = 0;
 335   Float_t ptLeadMax = -1;
 336   Int_t axis = 2;
 337
 338   if (compareUEHist == 2)
 339   {
 340     ptLeadMin = 1.01;
 341     ptLeadMax = 1.49;
 342     axis = 4;
 343   }
 344
 345   TH1* corrHist = corr->GetUEHist(compareStep, compareRegion, ptLeadMin, ptLeadMax);
 346   TH1* mcHist   = mc  ->GetUEHist(compareStep, compareRegion, ptLeadMin, ptLeadMax);
 347
 348   DrawRatio(corrHist, mcHist, Form("check"));
 349
 350   corr->SetBinLimits(corr->GetTrackHist(compareRegion)->GetGrid(compareStep));
 351   mc->SetBinLimits(corr->GetTrackHist(compareRegion)->GetGrid(compareStep));
 352
 353   corrHist =  corr->GetTrackHist(compareRegion)->GetGrid(compareStep)->Project(axis);
 354   mcHist   =  mc  ->GetTrackHist(compareRegion)->GetGrid(compareStep)->Project(axis);
 355   DrawRatio(corrHist, mcHist, Form("check2"));
 356
 357   corrHist =  corr->GetEventHist()->GetGrid(compareStep)->Project(0);
 358   mcHist   =  mc  ->GetEventHist()->GetGrid(compareStep)->Project(0);
 359   DrawRatio(corrHist, mcHist, Form("check3"));
 360 }
 361
 362 void correctMC(const char* fileNameCorrections, const char* fileNameESD = 0, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 363 {
 364   // corrects the reconstructed step in fileNameESD with fileNameCorr
 365   // if fileNameESD is 0 data from fileNameCorr is taken
 366   // afterwards the corrected distributions are compared with the MC stored in fileNameESD
 367
 368   loadlibs();
 369
 370   TFile::Open(fileNameCorrections);
 371   list = (TList*) gFile->Get(objectName);
 372   AliUEHistograms* corr = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 373   SetupRanges(corr);
 374   corr->ExtendTrackingEfficiency();
 375
 376   if (fileNameESD)
 377   {
 378     file2 = TFile::Open(fileNameESD);
 379     list = (TList*) file2->Get(objectName);
 380   }
 381   AliUEHistograms* testSample = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 382
 383   // copy to esd object
 384   AliUEHistograms* esd = new AliUEHistograms;
 385   esd->CopyReconstructedData(testSample);
 386
 387   SetupRanges(corr);
 388   SetupRanges(testSample);
 389   SetupRanges(esd);
 390
 391   esd->Correct(corr);
 392
 393   DrawRatios(esd, testSample, compareStep, compareRegion, compareUEHist);
 394
 395   //CompareEventsTracks(esd, testSample, compareStep, compareRegion, compareUEHist);
 396 }
 397
 398 // function to compare only final step for all regions and distributions
 399
 400 void correctData(const char* fileNameCorrections, const char* fileNameESD, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
 401 {
 402   // corrects fileNameESD with fileNameCorrections and compares the two
 403
 404   loadlibs();
 405
 406   TFile::Open(fileNameCorrections);
 407   list = (TList*) gFile->Get(objectName);
 408   AliUEHistograms* corr = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 409
 410   file2 = TFile::Open(fileNameESD);
 411   list = (TList*) file2->Get(objectName);
 412   AliUEHistograms* esd = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 413
 414   SetupRanges(corr);
 415   SetupRanges(esd);
 416
 417   corr->ExtendTrackingEfficiency();
 418
 419   esd->Correct(corr);
 420
 421   file3 = TFile::Open("corrected.root", "RECREATE");
 422   file3->mkdir("PWG4_LeadingTrackUE");
 423   file3->cd("PWG4_LeadingTrackUE");
 424   list->Write(0, TObject::kSingleKey);
 425   file3->Close();
 426
 427   DrawRatios(esd, corr, compareStep, compareRegion, compareUEHist);
 428 }
 429
 430 void ITSTPCEfficiency(const char* fileNameData, Int_t id, Int_t itsTPC = 0)
 431 {
 432   // its = 0; tpc = 1
 433
 434   // uncertainty from dN/dpT paper
 435   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};
 436   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
 437   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
 438
 439   TH1F* effHist = new TH1F("effHist", "effHist", 39, pTBins);
 440   for (Int_t i=0; i<39; i++)
 441   {
 442     Int_t bin = i;
 443     if (i > 16)
 444       bin = 16;
 445     effHist->SetBinContent(i+1, (itsTPC == 0) ? effITS[bin] : effTPC[bin]);
 446   }
 447
 448   new TCanvas; effHist->Draw();
 449
 450   EffectOfModifiedTrackingEfficiency(fileNameData, id, effHist, 1, -1, (itsTPC == 0) ? "ITS" : "TPC");
 451 }
 452
 453
 454 void EffectOfModifiedTrackingEfficiency(const char* fileNameData, Int_t id, TH1* trackingEff, Int_t axis1, Int_t axis2 = -1, const char* name = "EffectOfModifiedTrackingEfficiency")
 455 {
 456   // trackingEff should contain the change in tracking efficiency, i.e. between before and after in the eta-pT plane
 457
 458   loadlibs();
 459
 460   TFile::Open(fileNameData);
 461   list = (TList*) gFile->Get(objectName);
 462   AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 463   SetupRanges(corrected);
 464
 465   AliUEHist* ueHist = corrected->GetUEHist(id);
 466   Int_t region = 2;
 467
 468   // histogram before
 469   TH1* before = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
 470
 471   // copy histogram
 472   // the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
 473   ueHist->CorrectTracks(AliUEHist::kCFStepAll, AliUEHist::kCFStepTriggered, (TH1*) 0, 0, -1);
 474
 475   // reapply tracking efficiency
 476   ueHist->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, trackingEff, axis1, axis2);
 477
 478   // histogram after
 479   TH1* after = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
 480
 481   DrawRatio(before, after, name);
 482   gPad->GetCanvas()->SaveAs(Form("%s.png", name));
 483 }
 484
 485 void EffectOfTrackCuts(const char* fileNameData, Int_t id, const char* systFile)
 486 {
 487   loadlibs();
 488
 489   TFile::Open(fileNameData);
 490   list = (TList*) gFile->Get(objectName);
 491   AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 492   effHist = (TH2D*) corrected->GetUEHist(0)->GetTrackingEfficiency()->Clone("effHist");
 493
 494   file = TFile::Open(systFile);
 495
 496   Int_t maxSyst = 3;
 497   const char* systNames[] = { "NClusTPC", "Chi2TPC", "SigmaDCA" };
 498
 499   for (Int_t i=0; i<maxSyst; i++)
 500   {
 501     for (Int_t j=0; j<2; j++)
 502     {
 503       TString histName;
 504       histName.Form("%s_syst_%s", systNames[i], (j == 0) ? "up" : "down");
 505       systEffect = (TH2*) file->Get(histName);
 506
 507       // rebin
 508       effHist->Reset();
 509       for (Int_t x=1; x <= effHist->GetNbinsX(); x++)
 510         for (Int_t y=1; y <= effHist->GetNbinsY(); y++)
 511           effHist->SetBinContent(x, y, 1);
 512
 513       for (Int_t x=1; x <= systEffect->GetNbinsX(); x++)
 514         for (Int_t y=1; y <= systEffect->GetNbinsY(); y++)
 515           if (systEffect->GetBinContent(x, y) != 0)
 516             effHist->SetBinContent(effHist->GetXaxis()->FindBin(systEffect->GetYaxis()->GetBinCenter(y)), effHist->GetYaxis()->FindBin(systEffect->GetXaxis()->GetBinCenter(x)), systEffect->GetBinContent(x, y));
 517
 518
 519       //new TCanvas; systEffect->Draw("COLZ"); new TCanvas; effHist->Draw("COLZ");
 520
 521       EffectOfModifiedTrackingEfficiency(fileNameData, id, effHist, 0, 1, histName);
 522
 523       //return;
 524     }
 525   }
 526 }
 527
 528 void ModifyComposition(const char* fileNameData, const char* fileNameCorrections, Int_t id, Bool_t verbose = kFALSE)
 529 {
 530   loadlibs();
 531
 532   TFile::Open(fileNameData);
 533   list = (TList*) gFile->Get(objectName);
 534   AliUEHistograms* corrected = (AliUEHistograms*) list->FindObject("AliUEHistograms");
 535   SetupRanges(corrected);
 536
 537   file2 = TFile::Open(fileNameCorrections);
 538   list2 = (TList*) file2->Get(objectName);
 539   AliUEHistograms* corrections = (AliUEHistograms*) list2->FindObject("AliUEHistograms");
 540   SetupRanges(corrections);
 541
 542   ueHistData        = (AliUEHist*) corrected->GetUEHist(id);
 543   ueHistCorrections = (AliUEHist*) corrections->GetUEHist(id);
 544
 545   // copy histogram
 546   // the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
 547   ueHistData->CorrectTracks(AliUEHist::kCFStepAll, AliUEHist::kCFStepTriggered, (TH1*) 0, 0);
 548
 549   Int_t maxRegion = 3;
 550   Float_t ptLeadMin = -1;
 551   Float_t ptLeadMax = -1;
 552   if (id == 2)
 553   {
 554     maxRegion = 1;
 555     // the uncertainty is flat in delta phi, so use this trick to get directly the uncertainty as function of leading pT
 556     //ptLeadMin = 1.01;
 557     //ptLeadMax = 1.99;
 558   }
 559
 560   // histogram before
 561   TH1* before[3];
 562   for (Int_t region=0; region<maxRegion; region++)
 563     before[region] = ueHistData->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
 564
 565   defaultEff = ueHistCorrections->GetTrackingEfficiency();
 566   defaultEffpT = ueHistCorrections->GetTrackingEfficiency(1);
 567   defaultContainer = ueHistCorrections->GetTrackHistEfficiency();
 568
 569   c = new TCanvas;
 570   defaultEffpT->Draw("");
 571
 572   Float_t largestDeviation[3];
 573   for (Int_t i=0; i<maxRegion; i++)
 574     largestDeviation[i] = 0;
 575
 576   for (Int_t i=0; i<7; i++)
 577   {
 578     // case 0: // default
 579     // case 1: // + 30% kaons
 580     // case 2: // - 30% kaons
 581     // case 3: // + 30% protons
 582     // case 4: // - 30% protons
 583     // case 5: // + 30% others
 584     // case 6: // - 30% others
 585     Int_t correctionIndex = (i+1) / 2 + 1; // bin 1 == protons, bin 2 == kaons, ...
 586     Double_t scaleFactor = (i % 2 == 1) ? 1.3 : 0.7;
 587     if (i == 0)
 588       scaleFactor = 1;
 589
 590     newContainer = (AliCFContainer*) defaultContainer->Clone();
 591
 592     // modify, change all steps
 593     for (Int_t j=0; j<newContainer->GetNStep(); j++)
 594     {
 595       THnSparse* grid = newContainer->GetGrid(j)->GetGrid();
 596
 597       for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 598       {
 599         Int_t bins[5];
 600         Double_t value = grid->GetBinContent(binIdx, bins);
 601         Double_t error = grid->GetBinError(binIdx);
 602
 603         if (bins[2] != correctionIndex)
 604           continue;
 605
 606         value *= scaleFactor;
 607         error *= scaleFactor;
 608
 609         grid->SetBinContent(bins, value);
 610         grid->SetBinError(bins, error);
 611       }
 612     }
 613
 614     // put in corrections
 615     ueHistCorrections->SetTrackHistEfficiency(newContainer);
 616
 617     // ratio
 618     modifiedEff = ueHistCorrections->GetTrackingEfficiency();
 619     modifiedEff->Divide(modifiedEff, defaultEff);
 620     //modifiedEff->Draw("COLZ");
 621
 622     c->cd();
 623     modifiedEffpT = ueHistCorrections->GetTrackingEfficiency(1);
 624     modifiedEffpT->SetLineColor(i+1);
 625     modifiedEffpT->Draw("SAME");
 626
 627     // apply change in tracking efficiency
 628     ueHistData->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, modifiedEff, 0, 1);
 629
 630     for (Int_t region=0; region<maxRegion; region++)
 631     {
 632       // histogram after
 633       TH1* after = ueHistData->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
 634
 635       if (verbose)
 636         DrawRatio(before[region], (TH1*) after->Clone(), Form("Region %d Composition %d", region, i));
 637
 638       // ratio is flat, extract deviation
 639       after->Divide(before[region]);
 640       after->Fit("pol0", "0");
 641       Float_t deviation = 100.0 - 100.0 * after->GetFunction("pol0")->GetParameter(0);
 642       Printf("Deviation for region %d case %d is %.2f %%", region, i, deviation);
 643
 644       if (TMath::Abs(deviation) > largestDeviation[region])
 645         largestDeviation[region] = TMath::Abs(deviation);
 646     }
 647     //return;
 648   }
 649
 650   for (Int_t i=0; i<maxRegion; i++)
 651     Printf("Largest deviation in region %d is %f", i, largestDeviation[i]);
 652 }