+TH1* GetSystematicUncertainty(TH1* corr, TH1* trackHist)
+{
+ if (!trackHist)
+ {
+ systError = (TH1*) corr->Clone("systError");
+ systError->Reset();
+ }
+ else // for dphi evaluation
+ systError = new TH1F("systError", "", 100, 0, 50);
+
+ Float_t constantUnc = 0;
+
+ // particle composition
+ constantUnc += 0.8 ** 2;
+
+ // tpc efficiency
+ if (gEnergy == 900 && gLeadingpTMin < 1.0)
+ constantUnc += 1.0 ** 2;
+ else if (gEnergy == 900 && gLeadingpTMin < 1.5)
+ constantUnc += 0.5 ** 2;
+ if (gEnergy == 7000 && gLeadingpTMin < 1.0)
+ constantUnc += 1.0 ** 2;
+ else if (gEnergy == 7000 && gLeadingpTMin < 1.5)
+ constantUnc += 0.6 ** 2;
+
+ // track cuts
+ if (gEnergy == 900 && gLeadingpTMin < 1.0)
+ constantUnc += 2.5 ** 2;
+ else if (gEnergy == 900 && gLeadingpTMin < 1.5)
+ constantUnc += 2.0 ** 2;
+ if (gEnergy == 7000)
+ constantUnc += 3.0 ** 2;
+
+ // difference corrected with pythia and phojet
+ if (gEnergy == 900 && gLeadingpTMin < 1.0)
+ constantUnc += 0.6 ** 2;
+ else if (gEnergy == 900 && gLeadingpTMin < 1.5)
+ constantUnc += 0.8 ** 2;
+
+ if (gEnergy == 7000 && gLeadingpTMin < 1.0)
+ {
+ if (gUEHist == 0)
+ constantUnc += 0.6 ** 2;
+ if (gUEHist == 1)
+ constantUnc += 0.8 ** 2;
+ if (gUEHist == 2)
+ constantUnc += 1.0 ** 2;
+ }
+ else if (gEnergy == 7000 && gLeadingpTMin < 1.5)
+ constantUnc += 1.0 ** 2;
+
+ for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
+ systError->SetBinContent(bin, constantUnc);
+
+ // mis id bias
+ if (gUEHist == 0 || gUEHist == 2)
+ systError->Fill(0.75, 4.0 ** 2);
+ if (gUEHist == 1)
+ systError->Fill(0.75, 5.0 ** 2);
+
+ if (gEnergy == 900)
+ {
+ if (gLeadingpTMin < 1.0)
+ systError->Fill(1.25, 1.0 ** 2);
+ else if (gLeadingpTMin < 1.5)
+ systError->Fill(1.25, 2.0 ** 2);
+ }
+
+ // non-closure in MC
+ if (gEnergy == 900)
+ for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
+ systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 1.0 ** 2);
+
+ if (gEnergy == 7000)
+ {
+ if (gUEHist == 0 && gUEHist == 1)
+ systError->Fill(0.75, 2.0 ** 2);
+ if (gUEHist == 2)
+ systError->Fill(0.75, 1.2 ** 2);
+ }
+
+ // vertex efficiency
+ systError->Fill(0.75, 1.0 ** 2);
+
+ // strangeness
+ for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
+ {
+ if (gEnergy == 900)
+ systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 0.5 ** 2);
+ if (gEnergy == 7000 && systError->GetXaxis()->GetBinCenter(bin) < 1.5)
+ systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 2.0 ** 2);
+ else if (gEnergy == 7000)
+ systError->Fill(systError->GetXaxis()->GetBinCenter(bin), 1.0 ** 2);
+ }
+
+ for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
+ systError->SetBinContent(bin, TMath::Sqrt(systError->GetBinContent(bin)));
+
+ if (trackHist)
+ {
+ //new TCanvas; trackHist->Draw();
+ //new TCanvas; systError->DrawCopy("");
+
+ Float_t uncFlat = 0;
+ for (Int_t i=1; i<=trackHist->GetNbinsX(); i++)
+ uncFlat += trackHist->GetBinContent(i) * systError->GetBinContent(systError->FindBin(trackHist->GetBinCenter(i)));
+ uncFlat /= trackHist->Integral();
+
+ systError = (TH1F*) corr->Clone("systError");
+ systError->Reset();
+
+ for (Int_t i=1; i<=systError->GetNbinsX(); i++)
+ systError->SetBinContent(i, uncFlat);
+
+ //new TCanvas; systError->DrawCopy("");
+ }
+
+ systError->SetFillColor(kGray);
+ systError->SetFillStyle(1001);
+ systError->SetMarkerStyle(0);
+ systError->SetLineColor(0);
+
+ return systError;
+}
+
+void DrawRatio(TH1* corr, TH1* mc, const char* name, TH1* syst = 0)