const char* objectName = "PWG4_LeadingTrackUE/histosLeadingTrackUE";
Float_t gLeadingpTMin = 0.51;
+Int_t gUEHist = 0;
+Bool_t gCache = 0;
+void* gFirst = 0;
+void* gSecond = 0;
+Float_t gForceRange = -1;
+Int_t gEnergy = 900;
+Int_t gRegion = 0;
+
//const char* objectName = "PWG4_LeadingTrackUE/histosLeadingTrackUE_filterbit32";
void SetRanges(TAxis* axis)
{
if (strcmp(axis->GetTitle(), "leading p_{T} (GeV/c)") == 0)
- axis->SetRangeUser(0, 10);
+ axis->SetRangeUser(0, (gEnergy == 900) ? 10 : 25);
if (strcmp(axis->GetTitle(), "multiplicity") == 0)
axis->SetRangeUser(0, 10);
SetRanges(hist);
}
-void DrawRatio(TH1* corr, TH1* mc, const char* name)
+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)
{
mc->SetLineColor(2);
pad1->SetGridx();
pad1->SetGridy();
- TLegend* legend = new TLegend(0.15, 0.65, 0.55, 0.90);
+ if (gUEHist != 2)
+ TLegend* legend = new TLegend(0.15, 0.65, 0.55, 0.90);
+ else
+ TLegend* legend = new TLegend(0.55, 0.65, 0.95, 0.90);
+
legend->SetFillColor(0);
legend->AddEntry(corr, "Corrected");
legend->AddEntry(mc, "MC prediction");
dummy->GetXaxis()->SetTitleSize(0.08);
dummy->GetYaxis()->SetTitleSize(0.08);
dummy->GetYaxis()->SetTitleOffset(0.8);
+
+ if (gForceRange > 0)
+ dummy->GetYaxis()->SetRangeUser(0, gForceRange);
+
dummy->DrawCopy();
-
+
+ if (gUEHist != 2)
+ {
+ const char* regionStr[] = { "Toward Region", "Away Region", "Transverse Region" };
+ latex = new TLatex(0.65, 0.1, regionStr[gRegion]);
+ latex->SetTextSize(0.075);
+ latex->SetNDC();
+ latex->Draw();
+ }
+
+ if (syst)
+ {
+ systError = (TH1*) syst->Clone("corrSystError");
+ for (Int_t bin=1; bin<=systError->GetNbinsX(); bin++)
+ {
+ systError->SetBinError(bin, corr->GetBinContent(bin) * syst->GetBinContent(bin) / 100);
+ systError->SetBinContent(bin, corr->GetBinContent(bin));
+ }
+
+ systError->Draw("E2 ][ SAME");
+ }
+
corr->Draw("SAME");
mc->Draw("SAME");
Float_t minR = TMath::Min(0.91, ratio->GetMinimum() * 0.95);
Float_t maxR = TMath::Max(1.09, ratio->GetMaximum() * 1.05);
- minR = 0.8;
- maxR = 1.2;
+ minR = 0.6;
+ maxR = 1.4;
TH1F dummy3("dummy3", ";;Ratio: MC / corr", 100, corr->GetXaxis()->GetBinLowEdge(1), corr->GetXaxis()->GetBinUpEdge(corr->GetNbinsX()));
dummy3.SetXTitle(corr->GetXaxis()->GetTitle());
dummy3.GetYaxis()->SetTitleSize(0.08);
dummy3.GetYaxis()->SetTitleOffset(0.8);
dummy3.DrawCopy();
+
+ if (syst)
+ {
+ // for the ratio add in quadrature
+ for (Int_t bin=1; bin<=syst->GetNbinsX(); bin++)
+ {
+ if (corr->GetBinError(bin) > 0)
+ syst->SetBinError(bin, TMath::Sqrt(TMath::Power(syst->GetBinContent(bin) / 100, 2) + TMath::Power(corr->GetBinError(bin) / corr->GetBinContent(bin), 2)));
+ else
+ syst->SetBinError(bin, 0);
+ syst->SetBinContent(bin, 1);
+ }
+
+ syst->Draw("E2 ][ SAME");
+ dummy3.DrawCopy("SAME");
+ }
ratio->Draw("SAME");
+
+ ratio->Fit("pol0", "N");
canvas3->Modified();
+
}
void loadlibs()
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)
{
loadlibs();
-
- file1 = TFile::Open(fileName1);
- list = (TList*) file1->Get(objectName);
- AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
- SetupRanges(h);
-
- file2 = TFile::Open(fileName2);
- list = (TList*) file2->Get(objectName);
- AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
- SetupRanges(h2);
-
+
+ if (!gCache || !gFirst)
+ {
+ file1 = TFile::Open(fileName1);
+ list = (TList*) file1->Get(objectName);
+ AliUEHistograms* h = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+ SetupRanges(h);
+
+ file2 = TFile::Open(fileName2);
+ list = (TList*) file2->Get(objectName);
+ AliUEHistograms* h2 = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+ SetupRanges(h2);
+ }
+
+ if (gCache)
+ {
+ if (!gFirst)
+ {
+ gFirst = h;
+ gSecond = h2;
+ }
+ else
+ {
+ AliUEHistograms* h = (AliUEHistograms*) gFirst;
+ AliUEHistograms* h2 = (AliUEHistograms*) gSecond;
+ }
+ }
+
TH1* hist1 = h->GetUEHist(id)->GetUEHist(step1, region, ptLeadMin, ptLeadMax);
TH1* hist2 = h2->GetUEHist(id)->GetUEHist(step2, region, ptLeadMin, ptLeadMax);
- DrawRatio(hist1, hist2, "compare");
+ TH1* trackHist = 0;
+ if (id == 2)
+ {
+ trackHist = h->GetUEHist(id)->GetTrackHist(region)->ShowProjection(2, 0);
+ // only keep bins under consideration
+ for (Int_t bin=1; bin<=trackHist->GetNbinsX(); bin++)
+ if (bin < trackHist->FindBin(ptLeadMin) || bin > trackHist->FindBin(ptLeadMax))
+ trackHist->SetBinContent(bin, 0);
+ }
+
+ // systematic uncertainty
+ TH1* syst = GetSystematicUncertainty(hist1, trackHist);
+
+ DrawRatio(hist1, hist2, Form("%s_%d_%d_%d_%.2f_%.2f", TString(fileName1).Tokenize(".")->First()->GetName(), id, step1, region, ptLeadMin, ptLeadMax), syst);
}
void CompareEventHist(const char* fileName1, const char* fileName2, Int_t id, Int_t step, Int_t var)
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)
{
+ // fileName1 is labelled Corrected in the plot
+
loadlibs();
+ gUEHist = id;
+ gRegion = region;
Compare(fileName1, fileName2, id, step, step, region, ptLeadMin, ptLeadMax);
}
AliUEHistograms* corrected = (AliUEHistograms*) correctedVoid;
AliUEHistograms* comparison = (AliUEHistograms*) comparisonVoid;
- if (compareUEHist == 2)
+ if (1 && compareUEHist == 2)
{
- for (Float_t ptLeadMin = 1.01; ptLeadMin < 4; ptLeadMin += 2)
+ for (Float_t ptLeadMin = 0.51; ptLeadMin < 10; ptLeadMin += 1.5)
DrawRatios(TString(Form("UE %d pT %f", compareUEHist, ptLeadMin)), corrected->GetUEHist(compareUEHist), comparison->GetUEHist(compareUEHist), compareStep, compareRegion, ptLeadMin, ptLeadMin + 0.48);
return;
}
SetupRanges(esd);
esd->Correct(corr);
+ esd->GetUEHist(2)->AdditionalDPhiCorrection(0);
DrawRatios(esd, testSample, compareStep, compareRegion, compareUEHist);
// function to compare only final step for all regions and distributions
-void correctData(const char* fileNameCorrections, const char* fileNameESD, Int_t compareStep = -1, Int_t compareRegion = 2, Int_t compareUEHist = 0)
+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)
{
// corrects fileNameESD with fileNameCorrections and compares the two
corr->ExtendTrackingEfficiency();
+ TFile::Open(contEnhancement);
+ contEncHist = (TH1*) gFile->Get("histo");
+ contEncHistFullRange = (TH1*) corr->GetUEHist(0)->GetTrackingEfficiency(1)->Clone("contEncHistFullRange");
+
+ contEncHistFullRange->Reset();
+ for (Int_t i=1; i<=contEncHistFullRange->GetNbinsX(); i++)
+ {
+ contEncHistFullRange->SetBinContent(i, 1);
+ if (i <= contEncHist->GetNbinsX() && contEncHist->GetXaxis()->GetBinCenter(i) < contEncUpTo && contEncHist->GetBinContent(i) > 0)
+ contEncHistFullRange->SetBinContent(i, contEncHist->GetBinContent(i));
+ }
+ corr->SetContaminationEnhancement((TH1F*) contEncHistFullRange);
+
esd->Correct(corr);
+ esd->GetUEHist(2)->AdditionalDPhiCorrection(0);
file3 = TFile::Open("corrected.root", "RECREATE");
file3->mkdir("PWG4_LeadingTrackUE");
new TCanvas; effHist->Draw();
- EffectOfModifiedTrackingEfficiency(fileNameData, id, effHist, 1, -1, (itsTPC == 0) ? "ITS" : "TPC");
+ EffectOfModifiedTrackingEfficiency(fileNameData, id, 2, effHist, 1, -1, (itsTPC == 0) ? "ITS" : "TPC");
}
-void EffectOfModifiedTrackingEfficiency(const char* fileNameData, Int_t id, TH1* trackingEff, Int_t axis1, Int_t axis2 = -1, const char* name = "EffectOfModifiedTrackingEfficiency")
+void EffectOfModifiedTrackingEfficiency(const char* fileNameData, Int_t id, Int_t region, TH1* trackingEff, Int_t axis1, Int_t axis2 = -1, const char* name = "EffectOfModifiedTrackingEfficiency")
{
// trackingEff should contain the change in tracking efficiency, i.e. between before and after in the eta-pT plane
SetupRanges(corrected);
AliUEHist* ueHist = corrected->GetUEHist(id);
- Int_t region = 2;
+ Float_t ptLeadMin = -1;
+ Float_t ptLeadMax = -1;
+ if (id == 2)
+ {
+ // the uncertainty is flat in delta phi, so use this trick to get directly the uncertainty as function of leading pT
+ //ptLeadMin = 1.01;
+ //ptLeadMax = 1.99;
+ }
+
// histogram before
- TH1* before = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
+ TH1* before = ueHist->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
// copy histogram
// the CFStepTriggered step is overwritten here and cannot be used for comparison afterwards anymore
ueHist->CorrectTracks(AliUEHist::kCFStepTriggered, AliUEHist::kCFStepAll, trackingEff, axis1, axis2);
// histogram after
- TH1* after = ueHist->GetUEHist(AliUEHist::kCFStepAll, region);
+ TH1* after = ueHist->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
DrawRatio(before, after, name);
gPad->GetCanvas()->SaveAs(Form("%s.png", name));
//new TCanvas; systEffect->Draw("COLZ"); new TCanvas; effHist->Draw("COLZ");
- EffectOfModifiedTrackingEfficiency(fileNameData, id, effHist, 0, 1, histName);
+ EffectOfModifiedTrackingEfficiency(fileNameData, id, (id == 2) ? 0 : 2, effHist, 0, 1, histName);
//return;
}
for (Int_t region=0; region<maxRegion; region++)
before[region] = ueHistData->GetUEHist(AliUEHist::kCFStepAll, region, ptLeadMin, ptLeadMax);
- defaultEff = ueHistCorrections->GetTrackingEfficiency();
- defaultEffpT = ueHistCorrections->GetTrackingEfficiency(1);
+ //defaultEff = ueHistCorrections->GetTrackingEfficiency();
+ defaultEff = ueHistCorrections->GetTrackingCorrection();
+ //defaultEffpT = ueHistCorrections->GetTrackingEfficiency(1);
+ defaultEffpT = ueHistCorrections->GetTrackingCorrection(1);
defaultContainer = ueHistCorrections->GetTrackHistEfficiency();
c = new TCanvas;
ueHistCorrections->SetTrackHistEfficiency(newContainer);
// ratio
- modifiedEff = ueHistCorrections->GetTrackingEfficiency();
+ //modifiedEff = ueHistCorrections->GetTrackingEfficiency();
+ modifiedEff = ueHistCorrections->GetTrackingCorrection();
modifiedEff->Divide(modifiedEff, defaultEff);
//modifiedEff->Draw("COLZ");
c->cd();
- modifiedEffpT = ueHistCorrections->GetTrackingEfficiency(1);
+ //modifiedEffpT = ueHistCorrections->GetTrackingEfficiency(1);
+ modifiedEffpT = ueHistCorrections->GetTrackingCorrection(1);
modifiedEffpT->SetLineColor(i+1);
modifiedEffpT->Draw("SAME");
for (Int_t i=0; i<maxRegion; i++)
Printf("Largest deviation in region %d is %f", i, largestDeviation[i]);
}
+
+void MergeList()
+{
+ loadlibs();
+
+ ifstream in;
+ in.open("list");
+
+ TFileMerger m;
+
+ TString line;
+ while (in.good())
+ {
+ in >> line;
+
+ if (line.Length() == 0)
+ continue;
+
+ TString fileName;
+ fileName.Form("%s/%s/PWG4_JetTasksOutput.root", "maps", line.Data());
+ Printf("%s", fileName.Data());
+
+ m.AddFile(fileName);
+ }
+
+ m.SetFastMethod();
+ m.OutputFile("merged.root");
+ m.Merge();
+}
+
+void MergeList2(const char* listFile, Bool_t onlyPrintEvents = kFALSE)
+{
+ loadlibs();
+
+ ifstream in;
+ in.open(listFile);
+
+ AliUEHistograms* final = 0;
+ TList* finalList = 0;
+
+ TString line;
+ while (in.good())
+ {
+ in >> line;
+
+ if (line.Length() == 0)
+ continue;
+
+ TString fileName;
+ fileName.Form("%s/%s/PWG4_JetTasksOutput.root", "maps", line.Data());
+ Printf("%s", fileName.Data());
+
+ file = TFile::Open(fileName);
+ if (!file)
+ continue;
+
+ list = (TList*) file->Get(objectName);
+
+ if (!final)
+ {
+ final = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+ //final->GetEventCount()->Draw(); return;
+ Printf("Events: %d", (Int_t) final->GetEventCount()->ProjectionX()->GetBinContent(4));
+ finalList = list;
+ }
+ else
+ {
+ additional = (AliUEHistograms*) list->FindObject("AliUEHistograms");
+ Printf("Events: %d", (Int_t) additional->GetEventCount()->ProjectionX()->GetBinContent(4));
+
+ if (!onlyPrintEvents)
+ {
+ TList list2;
+ list2.Add(additional);
+ final->Merge(&list2);
+ }
+ delete additional;
+ file->Close();
+ }
+ }
+
+ if (onlyPrintEvents)
+ return;
+
+ Printf("Total events (at step 0): %d", (Int_t) final->GetEventCount()->ProjectionX()->GetBinContent(4));
+
+ file3 = TFile::Open("merged.root", "RECREATE");
+ file3->mkdir("PWG4_LeadingTrackUE");
+ file3->cd("PWG4_LeadingTrackUE");
+ finalList->Write(0, TObject::kSingleKey);
+ file3->Close();
+}
+
+void PlotAll(const char* correctedFile, const char* mcFile)
+{
+ gCache = 1;
+
+ if (gEnergy == 900)
+ {
+ Float_t range[] = { 1.5, 2 };
+ }
+ else
+ {
+ Float_t range[] = { 3, 10 };
+ }
+
+ for (Int_t id=2; id<3; id++)
+ {
+ if (id < 2)
+ gForceRange = range[id];
+ else
+ gForceRange = -1;
+
+ if (id < 2)
+ {
+ for (Int_t region=0; region<3; region++)
+ {
+ CompareStep(correctedFile, mcFile, id, 0, region);
+ gPad->GetCanvas()->SaveAs(Form("%s_%s_%d_%d.png", TString(correctedFile).Tokenize(".")->First()->GetName(), TString(mcFile).Tokenize(".")->First()->GetName(), id, region));
+ }
+ }
+ else
+ {
+ Float_t leadingPtArr[] = { 0.50, 2.0, 4.0, 6.0, 10.0, 20.0, 50.0 };
+ for (Int_t leadingPtID=0; leadingPtID<6; leadingPtID++)
+ {
+ CompareStep(correctedFile, mcFile, id, 0, 0, leadingPtArr[leadingPtID] + 0.01, leadingPtArr[leadingPtID+1] - 0.01);
+ 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]));
+ }
+ }
+ }
+}
+
\ No newline at end of file