From 231aa7abbd6ff5721f26ead29f62331a5ac1064e Mon Sep 17 00:00:00 2001
From: Cornelius Feier-Riesen <cornelius.riesen@physik.uni-giessen.de>
Date: Thu, 3 Mar 2022 16:59:18 +0100
Subject: [PATCH] LMVM: add new PLUTO path (Tetyanas folder) add purity
histograms add histograms to investigate diff. between PLUTO/UrQMD and
electron/positron add CB histograms in LmvmDraw.cxx add histogram to
investigate BG and its sources add histosgrams for acc and rec. pions and
electrons for var. detector combinations some more small changes
---
analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h | 8 +-
analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx | 424 +++++++++++++-----
analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h | 16 +-
.../PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx | 312 ++++++++-----
analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h | 3 +-
.../PWGDIL/dielectron/lmvm/LmvmSimParam.h | 2 +-
analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx | 298 ++++++++----
analysis/PWGDIL/dielectron/lmvm/LmvmTask.h | 4 +-
analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx | 9 +-
macro/analysis/dielectron/batch_job.py | 18 +
macro/analysis/dielectron/batch_send.py | 6 +-
macro/analysis/dielectron/draw_all.py | 35 +-
macro/analysis/dielectron/draw_litqa.C | 4 +-
macro/analysis/dielectron/hadd_many.py | 6 +-
macro/analysis/dielectron/run_local.py | 14 +-
macro/analysis/dielectron/run_transport.C | 2 +-
.../cbm/elid/CbmLitGlobalElectronId.cxx | 25 +-
.../cbm/qa/tracking/CbmLitTrackingQa.cxx | 5 +-
18 files changed, 825 insertions(+), 366 deletions(-)
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h
index ff1383d03c..3698555841 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h
@@ -53,11 +53,15 @@ public:
return false;
}
// it is assumed that stationNum can be 1 or 2
+ LOG(info) << "I am in LmvmCuts::IsMvdCutOk"; // TODO: delete this line
double cutD = (stationNum == 1) ? fMvd1CutD : fMvd2CutD;
double cutP = (stationNum == 1) ? fMvd1CutP : fMvd2CutP;
double val = -1. * (cutP / cutD) * dmvd + cutP;
- if (!(dmvd < cutD && val > mom)) return true;
- return false;
+ if (!(dmvd < cutD && val > mom)) {LOG(info) << "MVD cut passed."; return true;} // TODO: delete cout
+ else { // TODO: delete cout and set back 'return' without else bracket
+ LOG(info) << "MVD cut not passed.";
+ return false;
+ }
}
std::string ToString()
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx
index 4db42c028b..520bc4ba10 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx
@@ -55,31 +55,36 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir,
fH.fHM.ScaleByPattern(".*", 1. / fNofEvents);
RebinMinvHist();
- DrawAnaStepMany("lmvm_pair_pty", [this](ELmvmAnaStep step) { DrawPtY(step); });
- DrawAnaStepMany("lmvm_pair_rapidity", [this](ELmvmAnaStep step) { DrawRapidity(step); });
- //DrawAnaStepMany("lmvm_pair_pty_efficiency", [this](ELmvmAnaStep step) { DrawPtYEfficiency(step); }); // TODO: causes segmentation violation error
- DrawAnaStepMany("lmvm_minv_sbg", [this](ELmvmAnaStep step) { DrawMinvSBg(step); });
- //DrawAnaStepMany("lmvm_minv_bgPairSrc", [this](ELmvmAnaStep step) { DrawMinvBgPairSrc(step); }); // TODO: causes segmentation violation error
- //DrawAnaStepMany("lmvm_minv_matching", [this](ELmvmAnaStep step) { DrawMinvMatching(step); }); // TODO: causes segmentation violation error
- DrawAnaStepMany("lmvm_minv_pt", [this](ELmvmAnaStep step) { DrawMinvPt(step); });
- DrawAnaStepMany("lmvm_anglePair", [this](ELmvmAnaStep step) { DrawSrcAnaStepH1("hAnglePair", step); });
+ //DrawAnaStepMany("pair_pty", [this](ELmvmAnaStep step) { DrawPtY(step); });
+ DrawAnaStepMany("pty_pair_signal", [this](ELmvmAnaStep step) { DrawPtY("hPtYPairSignal", step); });
+ DrawAnaStepMany("pty_cand", [this](ELmvmAnaStep step) { DrawPtY("hPtYCandidate", step); });
+ DrawAnaStepMany("pair_rapidity", [this](ELmvmAnaStep step) { DrawRapidity(step); });
+ //DrawAnaStepMany("pair_pty_efficiency", [this](ELmvmAnaStep step) { DrawPtYEfficiency(step); }); // TODO: causes segmentation violation error
+ DrawAnaStepMany("minv_sbg", [this](ELmvmAnaStep step) { DrawMinvSBg(step); });
+ //DrawAnaStepMany("minv_bgPairSrc", [this](ELmvmAnaStep step) { DrawMinvBgPairSrc(step); }); // TODO: causes segmentation violation error
+ //DrawAnaStepMany("minv_matching", [this](ELmvmAnaStep step) { DrawMinvMatching(step); }); // TODO: causes segmentation violation error
+ DrawAnaStepMany("minv_pt", [this](ELmvmAnaStep step) { DrawMinvPt(step); });
+ DrawAnaStepMany("anglePair", [this](ELmvmAnaStep step) { DrawSrcAnaStepH1("hAnglePair", step); });
// draw momentum histograms
for (const string& hName : {"hMom", "hMomPx", "hMomPy", "hMomPz", "hPt", "hRapidity"}) {
- DrawAnaStepMany("lmvm_" + hName, [this, hName](ELmvmAnaStep step) { DrawSrcAnaStepH1(hName, step); });
- DrawAnaStepMany("lmvm_" + hName + "EpEm", [this, hName](ELmvmAnaStep step) { DrawSrcAnaStepEpEmH1(hName, step); });
+ DrawAnaStepMany(hName, [this, hName](ELmvmAnaStep step) { DrawSrcAnaStepH1(hName, step); });
+ DrawAnaStepMany(hName + "EpEm", [this, hName](ELmvmAnaStep step) { DrawSrcAnaStepEpEmH1(hName, step); });
}
- DrawMomAccEpEm();
+ DrawElPurity();
+ DrawCombinatorialPairs();
DrawMisc();
DrawGammaVertex();
DrawCuts();
DrawMinvAll();
DrawBgSourceTracks();
+ DrawBgSourcePairsAll();
DrawMismatchesAndGhosts();
DrawMvdCutQa();
DrawMvdAndStsHist();
- DrawPiMom();
+ DrawAccRecMom();
DrawPmtXY();
+ DrawMinvBg(); // TODO: do not extra method
SaveCanvasToImage();
/// Restore old global file and folder pointer to avoid messing with FairRoot
@@ -95,13 +100,16 @@ bool LmvmDraw::SkipMvd(ELmvmAnaStep step)
void LmvmDraw::RebinMinvHist()
{
int nGroup = 20;
- int nGroupCB = 10; // rebin for CB histos
+ int nGroupCB = 100; // rebin for CB histos
int nGroupMatch = 50;
int nGroupBgSrc = 50;
fH.Rebin("hMinv", fH.fSrcNames, fH.fAnaStepNames, nGroup);
fH.Rebin("hMinvCombPM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, nGroupCB);
fH.Rebin("hMinvCombPP", {"sameEv", "mixedEv"}, fH.fAnaStepNames, nGroupCB);
fH.Rebin("hMinvCombMM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, nGroupCB);
+ fH.Rebin("hMinvCombPM", {"pluto", "urqmd"}, fH.fAnaStepNames, nGroupCB);
+ fH.Rebin("hMinvCombPP", {"pluto", "urqmd"}, fH.fAnaStepNames, nGroupCB);
+ fH.Rebin("hMinvCombMM", {"pluto", "urqmd"}, fH.fAnaStepNames, nGroupCB);
fH.Rebin("hMinvBgMatch", {"trueMatch", "trueMatchEl", "trueMatchNotEl", "mismatch"}, fH.fAnaStepNames, nGroupMatch);
fH.Rebin("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, nGroupBgSrc);
}
@@ -197,7 +205,6 @@ void LmvmDraw::DrawCutEffH1(const string& hist, const string& option)
DrawH1(effHist, fH.fSrcLatex, kLinear, kLinear, true, 0.8, 0.8, 0.99, 0.99, "hist");
}
-
void LmvmDraw::DrawAnaStepMany(const string& cName, function<void(ELmvmAnaStep)> drawFunc)
{
int hi = 1;
@@ -218,10 +225,10 @@ void LmvmDraw::DrawAnaStepMany(const string& cName, function<void(ELmvmAnaStep)>
}
}
-void LmvmDraw::DrawPtY(ELmvmAnaStep step)
+void LmvmDraw::DrawPtY(const string& hist, ELmvmAnaStep step)
{
- TH2D* h = fH.H2("hPtYPairSignal", step);
- TH2D* hmc = fH.H2("hPtYPairSignal", ELmvmAnaStep::Mc);
+ TH2D* h = fH.H2(hist.c_str(), step);
+ TH2D* hmc = fH.H2(hist.c_str(), ELmvmAnaStep::Mc);
DrawH2(h, kLinear, kLinear, kLinear, "COLZ");
bool drawAnaStep = true;
if (drawAnaStep) fH.DrawEfficiency(h, hmc, 0.2, 1.8);
@@ -273,39 +280,62 @@ void LmvmDraw::DrawSrcAnaStepEpEmH1(const string& hName, ELmvmAnaStep step)
fH.DrawAnaStepOnPad(step);
}
-void LmvmDraw::DrawMomAccEpEm()
-{
- for (const string& det : {"sts", "rich", "trd", "tof"}) {
- vector<TH1*> hVec;
- vector<string> latex;
- for (const ELmvmSrc src : {ELmvmSrc::Signal, ELmvmSrc::Pi0, ELmvmSrc::Gamma}) {
- for (const string& pm : {"+", "-"}) {
- hVec.push_back(fH.H1("hMomAcc" + pm + "_" + det, src));
- latex.push_back(fH.fSrcLatex[static_cast<int>(src)] + " (e" + pm + ")");
- }
- }
- fH.fHM.CreateCanvas("lmvm_momDetAcc/lmvm_momDetAcc_" + det, "lmvm_momDetAcc/lmvm_momDetAcc_" + det, 800, 800);
- DrawH1(hVec, latex, kLinear, kLog, true, 0.90, 0.7, 0.99, 0.99, "hist");
- DrawTextOnPad(det, 0.4, 0.9, 0.6, 0.999);
- }
-}
-
void LmvmDraw::DrawMisc()
{
- fH.fHM.CreateCanvas("lmvm_mom_pairSignal", "lmvm_mom_pairSignal", 800, 800);
+ fH.fHM.CreateCanvas("mom_pairSignal", "mom_pairSignal", 800, 800);
DrawAnaStepH1("hMomPairSignal", true);
- fH.fHM.CreateCanvas("lmvm_mother_pdg", "lmvm_mother_pdg", 800, 800);
+ fH.fHM.CreateCanvas("mother_pdg", "mother_pdg", 800, 800);
DrawH1({fH.H1("hMotherPdg_mc"), fH.H1("hMotherPdg_acc")}, {"MC", "acc"}, kLinear, kLog, true, 0.7, 0.7, 0.99, 0.99,
"hist");
- fH.fHM.CreateCanvas("lmvm_momVsAngle_pairSignal", "lmvm_momVsAngle_pairSignal", 800, 800);
+ fH.fHM.CreateCanvas("momVsAngle_pairSignal", "momVsAngle_pairSignal", 800, 800);
DrawH2(fH.H2("hMomVsAnglePairSignalMc"));
}
+void LmvmDraw::DrawCombinatorialPairs()
+{
+ // Draw same sign pairs, geometric mean and k factor seperate for PLUTO and UrQMD
+ for (const string& type :{"pluto", "urqmd"}) {
+ TH1* pm = fH.H1Clone("hMinvCombPM_" + type + "_ttcut");
+ TH1* pp = fH.H1Clone("hMinvCombPP_" + type + "_ttcut");
+ TH1* mm = fH.H1Clone("hMinvCombMM_" + type + "_ttcut");
+
+ TH1* rat = fH.H1Clone("hMinvCombPP_" + type + "_ttcut");
+ rat->Divide(mm);
+
+ fH.fHM.CreateCanvas("CombPairs/pairs_" + type, "CombPairs/pairs_" + type, 800, 800);
+ DrawH1({pm, pp, mm}, {"e^{+}e^{-} pairs", "e^{+}e^{+} pairs", "e^{-}e^{-} pairs"}, kLinear, kLog, true, 0.7, 0.7, 0.99, 0.99, "hist");
+ DrawTextOnPad(type, 0.4, 0.87, 0.6, 0.97);
+
+ fH.fHM.CreateCanvas("CombPairs/ratio_ppmm_" + type, "CombPairs/ratio_ppmm_" + type, 800, 800);
+ DrawH1(rat, kLinear, kLinear, "hist");
+ DrawTextOnPad(type, 0.4, 0.87, 0.6, 0.97);
+
+ TH1* gm = fH.H1Clone("hMinvCombPP_" + type + "_ttcut");
+
+ for (int i = 1; i <= gm->GetNbinsX(); i++) {
+ double con = std::sqrt(pp->GetBinContent(i) * mm->GetBinContent(i));
+ gm->SetBinContent(i, con);
+ }
+
+ TH1D* hK = (TH1D*) pm->Clone("hK");
+ hK->Divide(gm);
+ hK->Scale(0.5);
+
+ fH.fHM.CreateCanvas("CombPairs/kFactor_" + type, "CombPairs/kFactor_" + type, 800, 800);
+ DrawH1(hK, kLinear, kLinear, "p");
+ DrawTextOnPad(type, 0.4, 0.87, 0.6, 0.97);
+
+ fH.fHM.CreateCanvas("CombPairs/geomMean_" + type, "CombPairs/geomMean_" + type, 800, 800);
+ DrawH1(gm, kLinear, kLog, "hist");
+ DrawTextOnPad(type, 0.4, 0.87, 0.6, 0.97);
+ }
+}
+
void LmvmDraw::DrawPmtXY()
{
- TCanvas* c = fH.fHM.CreateCanvas("lmvm_pmtXY", "lmvm_pmtXY", 1800, 600);
+ TCanvas* c = fH.fHM.CreateCanvas("pmtXY", "pmtXY", 1800, 600);
c->Divide(3, 1);
vector<ELmvmSrc> src {ELmvmSrc::Signal, ELmvmSrc::Pi0, ELmvmSrc::Gamma};
for (size_t i = 0; i < src.size(); i++) {
@@ -333,7 +363,7 @@ void LmvmDraw::Draw1DCut(const string& hist, const string& sigOption, double cut
{
int w = 2400;
int h = 800;
- TCanvas* c = fH.fHM.CreateCanvas(("lmvm_cuts/lmvm_" + hist).c_str(), ("lmvm_cuts/lmvm_" + hist).c_str(), w, h);
+ TCanvas* c = fH.fHM.CreateCanvas(("cuts/" + hist).c_str(), ("cuts/" + hist).c_str(), w, h);
c->Divide(3, 1);
c->cd(1);
@@ -356,19 +386,17 @@ void LmvmDraw::Draw1DCut(const string& hist, const string& sigOption, double cut
void LmvmDraw::DrawCuts()
{
Draw1DCut("hAnnRich", "left", -0.4); // CbmLitGlobalElectronId::GetInstance().GetRichAnnCut()
- Draw1DCut("hAnnTrd", "left", 0.1); // CbmLitGlobalElectronId::GetInstance().GetTrdAnnCut()
+ //Draw1DCut("hAnnTrd", "left", 0.1); // CbmLitGlobalElectronId::GetInstance().GetTrdAnnCut() // TODO: uncomment when Trd Ann works again
+ Draw2DCut("hTrdLike_El");
+ Draw2DCut("hTrdLike_Pi");
Draw2DCut("hAnnRichVsMom");
Draw2DCut("hTofM2");
- DrawTofM2Cut();
-
+
Draw1DCut("hChi2PrimVertex", "right", fCuts.fChi2PrimCut);
//Draw1DCut("hPt", "left", fCuts.fPtCut);
//Draw1DCut("hMom", "left");
Draw1DCut("hChi2Sts", "right");
- Draw2DCut("hTrdLike_El");
- Draw2DCut("hTrdLike_Pi");
-
for (const string& type : {"all", "pion", "truePair"}) {
Draw2DCut("hStCut_" + type, fCuts.fStCutPP, fCuts.fStCutAngle);
Draw2DCut("hTtCut_" + type, fCuts.fTtCutPP, fCuts.fTtCutAngle);
@@ -382,16 +410,16 @@ void LmvmDraw::DrawCuts()
}
-void LmvmDraw::DrawSrcBgPairs(ELmvmAnaStep step, bool inPercent, bool drawAnaStep)
+void LmvmDraw::DrawBgSourcePairs(ELmvmAnaStep step, bool inPercent, bool drawAnaStep)
{
TH2D* h = fH.H2Clone("hSrcBgPairsEpEm", step);
gStyle->SetPaintTextFormat("4.1f");
- string labels[3] = {"#gamma", "#pi^{0}", "oth"};
- for (int i = 1; i <= 3; i++) {
+ string labels[4] = {"#gamma", "#pi^{0}", "#pi^{#pm}", "oth"};
+ for (int i = 1; i <= 4; i++) {
h->GetYaxis()->SetBinLabel(i, labels[i - 1].c_str());
h->GetXaxis()->SetBinLabel(i, labels[i - 1].c_str());
}
- h->SetMarkerSize(3);
+ h->SetMarkerSize(2.5);
if (inPercent) {
h->Scale(100. / h->Integral());
h->GetZaxis()->SetTitle("[%]");
@@ -406,34 +434,34 @@ void LmvmDraw::DrawSrcBgPairs(ELmvmAnaStep step, bool inPercent, bool drawAnaSte
if (drawAnaStep) fH.DrawAnaStepOnPad(step);
}
-void LmvmDraw::DrawSrcBgPairsAll()
+void LmvmDraw::DrawBgSourcePairsAll()
{
int hi = 1;
- TCanvas* c1 = fH.fHM.CreateCanvas("lmvm_srcBgPairs_abs", "lmvm_srcBgPairs_abs", 1500, 1500);
+ TCanvas* c1 = fH.fHM.CreateCanvas("bg/srcPairs_abs", "bg/srcPairs_abs", 1500, 1500);
c1->Divide(3, 3);
for (ELmvmAnaStep step : fH.fAnaSteps) {
if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc || SkipMvd(step)) continue;
c1->cd(hi++);
- DrawSrcBgPairs(step, false);
+ DrawBgSourcePairs(step, false);
}
hi = 1;
- TCanvas* c2 = fH.fHM.CreateCanvas("lmvm_srcBgPairs_perc", "lmvm_srcBgPairs_perc", 1500, 1500);
+ TCanvas* c2 = fH.fHM.CreateCanvas("bg/srcPairs_perc", "bg/srcPairs_perc", 1500, 1500);
c2->Divide(3, 3);
for (ELmvmAnaStep step : fH.fAnaSteps) {
if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc || SkipMvd(step)) continue;
- c1->cd(hi++);
- DrawSrcBgPairs(step, true);
+ c2->cd(hi++);
+ DrawBgSourcePairs(step, true);
}
- string ptCutName = fH.fAnaStepNames[static_cast<int>(ELmvmAnaStep::PtCut)];
- fH.fHM.CreateCanvas("lmvm_srcBgPairs_abs_" + ptCutName, "lmvm_srcBgPairs_abs_" + ptCutName, 900, 900);
- DrawSrcBgPairs(ELmvmAnaStep::PtCut, false, false);
+ string elidCutName = fH.fAnaStepNames[static_cast<int>(ELmvmAnaStep::ElId)];
+ fH.fHM.CreateCanvas("bg/srcPairs_abs_" + elidCutName, "bg/srcPairs_abs_" + elidCutName, 1100, 800);
+ DrawBgSourcePairs(ELmvmAnaStep::ElId, false, false);
- fH.fHM.CreateCanvas("lmvm_srcBgPairs_perc_" + ptCutName, "lmvm_srcBgPairs_perc_" + ptCutName, 900, 900);
- DrawSrcBgPairs(ELmvmAnaStep::PtCut, true, false);
+ fH.fHM.CreateCanvas("bg/srcPairs_perc_" + elidCutName, "bg/srcPairs_perc_" + elidCutName, 1100, 800);
+ DrawBgSourcePairs(ELmvmAnaStep::ElId, true, false);
- DrawBgSource2D("lmvm_srcBgPairs_2d", "hSrcBgPairs", fH.fBgPairSrcLatex, 1000., "Pairs per event x10^{3}");
+ DrawSource2D("bg/srcPairs_2d", "hSrcBgPairs", fH.fBgPairSrcLatex, 1000., "Pairs per event x10^{3}");
}
void LmvmDraw::Draw2DCutTriangle(double xCr, double yCr)
@@ -446,32 +474,148 @@ void LmvmDraw::Draw2DCutTriangle(double xCr, double yCr)
}
}
-void LmvmDraw::DrawTofM2Cut()
+void LmvmDraw::DrawMinvBg()
{
- vector<TLine*> lines {new TLine(0., 0.01, 1.3, 0.01), new TLine(1.3, 0.01, 2.5, 0.022)}; // set by hand
- TCanvas* c = fH.fHM.CreateCanvas("lmvm_cuts/lmvm_hTofM2_zoom", "lmvm_cuts/lmvm_hTofM2_zoom", 1600, 800);
- c->Divide(2, 1);
- int hi = 1;
- for (ELmvmSrc src : {ELmvmSrc::Signal, ELmvmSrc::Bg}) {
- TH2D* hist = fH.H2Clone("hTofM2", src);
+ TH1D* gg = fH.H1Clone("hMinvBgSource2_elid_gg");
+ TH1D* pipi = fH.H1Clone("hMinvBgSource2_elid_pipi");
+ TH1D* pi0pi0 = fH.H1Clone("hMinvBgSource2_elid_pi0pi0");
+ TH1D* oo = fH.H1Clone("hMinvBgSource2_elid_oo");
+ TH1D* gpi = fH.H1Clone("hMinvBgSource2_elid_gpi");
+ TH1D* gpi0 = fH.H1Clone("hMinvBgSource2_elid_gpi0");
+ TH1D* go = fH.H1Clone("hMinvBgSource2_elid_go");
+ TH1D* pipi0 = fH.H1Clone("hMinvBgSource2_elid_pipi0");
+ TH1D* pio = fH.H1Clone("hMinvBgSource2_elid_pio");
+ TH1D* pi0o = fH.H1Clone("hMinvBgSource2_elid_pi0o");
+
+ cout << "Entries gg: " << gg->GetEntries() << endl;
+ cout << "Entries pipi: " << pipi->GetEntries() << endl;
+ cout << "Entries pi0pi0: " << pi0pi0->GetEntries() << endl;
+ cout << "Entries oo: " << oo->GetEntries() << endl;
+ cout << "Entries gpi: " << gpi->GetEntries() << endl;
+ cout << "Entries gpi0: " << gpi0->GetEntries() << endl;
+ cout << "Entries go: " << go->GetEntries() << endl;
+ cout << "Entries pipi0: " << pipi0->GetEntries() << endl;
+ cout << "Entries pio: " << pio->GetEntries() << endl;
+ cout << "Entries pi0o: " << pi0o->GetEntries() << endl;
+
+ int reb = 50;
+
+ gg->Rebin(reb);
+ pi0pi0->Rebin(reb);
+ gpi0->Rebin(reb);
+ go->Rebin(reb);
+ pi0o->Rebin(reb);
+
+ gg->Scale(1. / reb);
+ pi0pi0->Scale(1. / reb);
+ gpi0->Scale(1. / reb);
+ go->Scale(1. / reb);
+ pi0o->Scale(1. / reb);
+
+ string cName = "minvBgSrc/minvBgSrc";
+ //vector<string> names = {"#gamma-#gamma", "#pi^{#pm}-#pi^{#pm}", "#pi^{0}-#pi^{0}", "o.-o.", "#gamma-#pi^{#pm}", "#gamma-#pi^{0}", "#gamma-o.", "#pi^{#pm}-#pi^{0}", "#pi^{#pm}-o.", "#pi^{0}-o.", "misid. #pi^{#pm}"};
+ vector<string> names = {"#gamma-#gamma", "#pi^{0}-#pi^{0}", "#gamma-#pi^{0}", "#gamma-o.", "#pi^{0}-o."};
+ fH.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1000, 1000);
+ //DrawH1({gg, pipi, pi0pi0, oo, gpi, gpi0, go, pipi0, pio, pi0o}, names, kLinear, kLog, true, 0.85, 0.7, 0.99, 0.99, "hist");
+ DrawH1({gg, pi0pi0, gpi0, go, pi0o}, names, kLinear, kLog, true, 0.85, 0.7, 0.99, 0.99, "hist");
+}
+
+void LmvmDraw::DrawElPurity()
+{
+ // All occuring PIDs
+ for (ELmvmAnaStep step : fH.fAnaSteps) {
+ if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue;
+ TCanvas* c = fH.fHM.CreateCanvas("purity/pid_"+ fH.fAnaStepNames[static_cast<int>(step)], "purity/pid_" + fH.fAnaStepNames[static_cast<int>(step)], 1600, 800);
+ c->Divide(2,1);
+ c->cd(1);
+ DrawH1(fH.H1("hCandPdg_" + fH.fAnaStepNames[static_cast<int>(step)]), kLinear, kLog, "hist text40");
+ c->cd(2);
+ TH1D* pdgZoom = fH.H1Clone("hCandPdg_" + fH.fAnaStepNames[static_cast<int>(step)]);
+ pdgZoom->GetXaxis()->SetRangeUser(-20., 20.);
+ DrawH1(pdgZoom, kLinear, kLog, "hist text40");
+ }
+
+ // PID vs momentum
+ vector<string> yLabel = {"e^{#pm}_{PLUTO}", "e^{#pm}_{UrQMD}", "#pi^{#pm}", "p", "K^{+}", "o." };
+ for (ELmvmAnaStep step : fH.fAnaSteps) {
+ if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue;
+ fH.fHM.CreateCanvas("purity/PidVsMom_" + fH.fAnaStepNames[static_cast<int>(step)], "purity/PidVsMom_" + fH.fAnaStepNames[static_cast<int>(step)], 800, 600);
+ TH2D* hPidMom = fH.H2("hCandPdgVsMom_" + fH.fAnaStepNames[static_cast<int>(step)]);
+ hPidMom->SetMinimum(5e-7);
+ DrawH2(hPidMom, kLinear, kLinear, kLog, "COLZ");
+ for (size_t y = 1; y <= yLabel.size(); y++) {
+ hPidMom->GetYaxis()->SetBinLabel(y, yLabel[y-1].c_str());
+ }
+ double nEl = hPidMom->Integral(1, hPidMom->GetXaxis()->GetNbins(), 2, 2); // do not count PLUTO particles
+ double purity = (nEl / hPidMom->Integral(1, hPidMom->GetXaxis()->GetNbins(), 2, hPidMom->GetYaxis()->GetNbins())) * 100;
+ DrawTextOnPad("Purity: " + Cbm::NumberToString(purity, 1) + " %", 0.1, 0.9, 0.45, 0.99);
+ }
+
+ // Purity vs momentum
+ int nBins = fH.H2("hCandPdgVsMom_elid")->GetXaxis()->GetNbins();
+ double xMin = fH.H2("hCandPdgVsMom_elid")->GetXaxis()->GetXmin();
+ double xMax = fH.H2("hCandPdgVsMom_elid")->GetXaxis()->GetXmax();
+ TH1D* purity = new TH1D("purity_Mom", "purity_Mom; P [GeV/c]; Purity [%]", nBins, xMin, xMax);
+ for (int i = 1; i <= purity->GetNbinsX(); i++) {
+ double nEl = fH.H2("hCandPdgVsMom", ELmvmAnaStep::ElId)->GetBinContent(i, 2);
+ double nAll = fH.H2("hCandPdgVsMom", ELmvmAnaStep::ElId)->Integral(i, i, 2, fH.H2("hCandPdgVsMom_elid")->GetYaxis()->GetNbins());
+ double val = (nAll != 0) ? 100 * nEl/nAll : 0.;
+ purity->SetBinContent(i, val);
+ }
+ purity->Rebin(5);
+ purity->Scale(1. / 5.);
+ fH.fHM.CreateCanvas("purity/purity_mom_elid", "purity/purity_mom_elid", 800, 800);
+ DrawH1(purity, kLinear, kLinear, "pe");
+
+ // Source of electron (PDG = +-11) candidates
+ DrawSource2D("purity/SrcTracksEl_2d", "hCandElSrc",
+ {"#gamma", "#pi^{0}", "#pi^{#pm}", "p", "K", "e^{#pm}_{sec}", "oth.", "signal"}, 1000.,
+ "Tracks per event x10^{3}");
+
+ // Occurency of Electrons and not-Electrons for various cut categories
+ for (const string& hName : {"hAnnRichVsMomPur", "hTrdElLikePur", "hTofM2Pur"}) {
+ string cName = "purity/cuts_" + hName;
+ TCanvas* c = fH.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 800);
+ c->Divide(3, 1);
+ int hi = 1;
+ for (const string& id : {"El", "Bg"}) {
+ TH2D* hist = fH.H2Clone(hName + "_" + id);
+ if (hName == "hTofM2Pur") {
+ hist->GetXaxis()->SetRangeUser(0., 3.);
+ hist->GetYaxis()->SetRangeUser(-0.1, 0.05);
+ }
+ c->cd(hi);
+ DrawH2(hist, kLinear, kLinear, kLog, "COLZ");
+ DrawTextOnPad(id.c_str(), 0.6, 0.89, 0.7, 0.99);
+
+ if (hName == "hTofM2Pur") { // TODO: these drawn lines needed?
+ vector<TLine*> lines {new TLine(0., 0.01, 1.3, 0.01), new TLine(1.3, 0.01, 2.5, 0.022)}; // set by hand
+ for (size_t i = 0; i < lines.size(); i++) {
+ lines[i]->SetLineWidth(2.);
+ lines[i]->Draw();
+ }
+ }
+ hi++;
+ }
c->cd(hi);
- hist->GetXaxis()->SetRangeUser(0., 2.5);
- hist->GetYaxis()->SetRangeUser(-0.08, 0.05);
- DrawH2(hist);
- DrawTextOnPad(fH.fSrcLatex[static_cast<int>(src)], 0.6, 0.89, 0.7, 0.99);
- for (size_t i = 0; i < lines.size(); i++) {
- lines[i]->SetLineWidth(2.);
- lines[i]->Draw();
+ TH2D* ratio = fH.H2Clone(hName + "_El");
+ ratio->Divide(fH.H2(hName + "_Bg"));
+ if (hName == "hTofM2Pur") {
+ ratio->GetXaxis()->SetRangeUser(0., 3.);
+ ratio->GetYaxis()->SetRangeUser(-0.1, 0.05);
}
- hi++;
+ DrawH2(ratio, kLinear, kLinear, kLog, "COLZ");
+ DrawTextOnPad("Ratio El/Bg", 0.4, 0.85, 0.8, 0.99);
}
}
void LmvmDraw::Draw2DCut(const string& hist, double cutCrossX, double cutCrossY)
{
- TCanvas* c = fH.fHM.CreateCanvas(("lmvm_cuts/lmvm_" + hist).c_str(), ("lmvm_cuts/lmvm_" + hist).c_str(), 1000, 1500);
+ TCanvas* c = fH.fHM.CreateCanvas(("cuts/" + hist).c_str(), ("cuts/" + hist).c_str(), 1000, 1500);
c->Divide(2, 3);
vector<TH1*> projX, projY;
+ projX.clear(); // TODO: clearing needed?
+ projY.clear();
vector<string> latex;
for (ELmvmSrc src : {ELmvmSrc::Signal, ELmvmSrc::Bg, ELmvmSrc::Pi0, ELmvmSrc::Gamma}) {
int srcInt = static_cast<int>(src);
@@ -481,7 +625,7 @@ void LmvmDraw::Draw2DCut(const string& hist, double cutCrossX, double cutCrossY)
DrawTextOnPad((Cbm::NumberToString(nofPerEvent, 2) + "/ev."), 0.1, 0.9, 0.5, 0.99);
DrawTextOnPad(fH.fSrcLatex[srcInt], 0.6, 0.89, 0.7, 0.99);
Draw2DCutTriangle(cutCrossX, cutCrossY);
- projX.push_back(fH.H2(hist, src)->ProjectionX());
+ projX.push_back(fH.H2(hist, src)->ProjectionX((hist + fH.fSrcLatex[static_cast<int>(src)]).c_str(), 1, fH.H2(hist, src)->GetYaxis()->GetNbins(), ""));
projY.push_back(fH.H2(hist, src)->ProjectionY());
latex.push_back(fH.fSrcLatex[srcInt]);
}
@@ -498,7 +642,7 @@ void LmvmDraw::DrawGammaVertex()
vector<TH2D*> xyz {fH.H2("hVertexGammaXZ", ELmvmAnaStep::Mc), fH.H2("hVertexGammaYZ", ELmvmAnaStep::Mc),
fH.H2("hVertexGammaXY", ELmvmAnaStep::Mc)};
- TCanvas* c = fH.fHM.CreateCanvas("lmvm_vertexGamma_mc", "lmvm_vertexGamma_mc", 1800, 600);
+ TCanvas* c = fH.fHM.CreateCanvas("vertexGamma_mc", "vertexGamma_mc", 1800, 600);
c->Divide(3, 1);
for (size_t i = 0; i < xyz.size(); i++) {
c->cd(i + 1);
@@ -507,7 +651,7 @@ void LmvmDraw::DrawGammaVertex()
gPad->SetLogz(true);
}
- TCanvas* cZ = fH.fHM.CreateCanvas("lmvm_vertexGamma_z", "lmvm_vertexGamma_z", 1500, 750);
+ TCanvas* cZ = fH.fHM.CreateCanvas("vertexGamma_z", "vertexGamma_z", 1500, 750);
cZ->Divide(2, 1);
int counter = 1;
for (ELmvmAnaStep step : {ELmvmAnaStep::Mc, ELmvmAnaStep::PtCut}) {
@@ -520,7 +664,7 @@ void LmvmDraw::DrawGammaVertex()
fH.DrawAnaStepOnPad(step);
}
- TCanvas* cZoom = fH.fHM.CreateCanvas("lmvm_vertexGamma_mc_zoom", "lmvm_vertexGamma_mc_zoom", 1800, 600);
+ TCanvas* cZoom = fH.fHM.CreateCanvas("vertexGamma_mc_zoom", "vertexGamma_mc_zoom", 1800, 600);
cZoom->Divide(3, 1);
for (size_t i = 0; i < xyz.size(); i++) {
TH2D* hZoom = (TH2D*) xyz[i]->Clone();
@@ -532,7 +676,7 @@ void LmvmDraw::DrawGammaVertex()
gPad->SetLogz(true);
}
- fH.fHM.CreateCanvas("lmvm_vertexGamma_rz_mc", "lmvm_vertexGamma_rz_mc", 900, 900);
+ fH.fHM.CreateCanvas("vertexGamma_rz_mc", "vertexGamma_rz_mc", 900, 900);
DrawH2(fH.H2("hVertexGammaRZ", ELmvmAnaStep::Mc));
fH.H2("hVertexGammaRZ", ELmvmAnaStep::Mc)->SetMinimum(1e-3);
gPad->SetLogz(true);
@@ -570,14 +714,14 @@ void LmvmDraw::DrawAnaStepH1(const string& name, bool logy)
void LmvmDraw::DrawMinvAll()
{
- TCanvas* c1 = fH.fHM.CreateCanvas("lmvm_minv_sbg_anaStep", "lmvm_minv_sbg_anaStep", 1200, 600);
+ TCanvas* c1 = fH.fHM.CreateCanvas("minv_sbg_anaStep", "minv_sbg_anaStep", 1200, 600);
c1->Divide(2, 1);
c1->cd(1);
DrawAnaStepH1(fH.GetName("hMinv", ELmvmSrc::Signal), true);
c1->cd(2);
DrawAnaStepH1(fH.GetName("hMinv", ELmvmSrc::Bg), true);
- TCanvas* c2 = fH.fHM.CreateCanvas("lmvm_minv_pi0_eta_gamma_anaStep", "lmvm_minv_pi0_eta_gamma_anaStep", 1800, 600);
+ TCanvas* c2 = fH.fHM.CreateCanvas("minv_pi0_eta_gamma_anaStep", "minv_pi0_eta_gamma_anaStep", 1800, 600);
c2->Divide(3, 1);
c2->cd(1);
DrawAnaStepH1(fH.GetName("hMinv", ELmvmSrc::Pi0), true);
@@ -587,15 +731,18 @@ void LmvmDraw::DrawMinvAll()
DrawAnaStepH1(fH.GetName("hMinv", ELmvmSrc::Gamma), true);
for (const string& name : {"sameEv", "mixedEv"}) {
- string cName = "lmvm_minv_combPairs_" + name;
+ string cName = "minv_combPairs_" + name;
TCanvas* c3 = fH.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 600);
c3->Divide(3, 1);
c3->cd(1);
DrawAnaStepH1("hMinvCombPM_" + name, true);
+ DrawTextOnPad("e^{+}e^{-} pairs " + name, 0.35, 0.9, 0.65, 0.99);
c3->cd(2);
DrawAnaStepH1("hMinvCombPP_" + name, true);
+ DrawTextOnPad("e^{+}e^{+} pairs " + name, 0.35, 0.9, 0.65, 0.99);
c3->cd(3);
DrawAnaStepH1("hMinvCombMM_" + name, true);
+ DrawTextOnPad("e^{-}e^{-} pairs " + name, 0.35, 0.9, 0.65, 0.99);
}
}
@@ -633,10 +780,7 @@ void LmvmDraw::DrawMinvBgPairSrc(ELmvmAnaStep step)
vector<TH1*> hists;
vector<string> latex;
for (int i = 0; i < fH.fNofBgPairSrc; i++) {
- stringstream ss;
- ss << "hMinvBgSource_" << fH.fBgPairSrcNames[i] << "_" << fH.fAnaStepNames[static_cast<int>(step)];
- hists.push_back(fH.H1(ss.str()));
- hists[i]->SetMinimum(1e-8);
+ hists.push_back(fH.H1("hMinvBgSource_" + fH.fBgPairSrcNames[i] + "_" + fH.fAnaStepNames[static_cast<int>(step)])); // segmentation violation error is caused by this specific histogram; works with others
string perc = Cbm::NumberToString(100. * hists[i]->GetEntries() / nofBg, 1);
latex.push_back(fH.fBgPairSrcLatex[i] + "(" + perc + "%)");
}
@@ -661,39 +805,73 @@ void LmvmDraw::DrawMinvMatching(ELmvmAnaStep step)
fH.DrawAnaStepOnPad(step);
}
-void LmvmDraw::DrawPiMom()
+void LmvmDraw::DrawAccRecMom()
{
- vector<string> subNames {"mc", "acc", "rec", "recOnlySts", "recStsRichTrd", "recStsRichTrdTof"};
- vector<string> latex {"MC", "Acc", "Rec", "Rec only STS", "Rec STS-RICH-TRD", "Rec STS-RICH-TRD-TOF"};
- vector<string> latexAll(latex.size()), latexPrim(latex.size());
- vector<TH1*> histsAll, histsPrim;
+ // Acceptance and reconstruction yields for various detector combinations
+ for (const int& pdg : {11, 211, 2212}) {
+ vector<string> subNames {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"};
+ vector<string> latex {"MC", "Acc", "Rec in STS", "Rec in STS-RICH", "Rec in STS-RICH-TRD", "Rec in STS-RICH-TRD-TOF"};
+ vector<string> latexAll(latex.size()), latexPrim(latex.size());
+ string hName = (pdg == 11) ? "hElMom" : (pdg == 211) ? "hPiMom" : "hProtonMom";
+ string cName = "AccRecMom/" + hName;
+ vector<TH1*> histsAll, histsPrim;
+ int i = 0;
+ for (const string& subName : subNames) {
+ TH1D* hAll = fH.H1(hName + "_all_" + subName);
+ hAll->SetMinimum(3e-6);
+ hAll->SetMaximum(50);
+ latexAll[i] = latex[i] + " (" + Cbm::NumberToString(hAll->GetEntries() / fNofEvents, 2) + "/ev.)";
+ histsAll.push_back(hAll);
+
+ TH1D* hPrim = fH.H1(hName + "_prim_" + subName);
+ hPrim->SetMinimum(3e-6);
+ hPrim->SetMaximum(50);
+ latexPrim[i] = latex[i] + " (" + Cbm::NumberToString(hPrim->GetEntries() / fNofEvents, 2) + "/ev.)";
+ histsPrim.push_back(hPrim);
+ i++;
+ }
- int i = 0;
- for (const string& subName : subNames) {
- TH1D* hAll = fH.H1("hPiMom_all_" + subName);
- latexAll[i] = latex[i] + " (" + Cbm::NumberToString(hAll->GetEntries() / fNofEvents, 2) + "/ev.)";
- histsAll.push_back(hAll);
-
- TH1D* hPrim = fH.H1("hPiMom_prim_" + subName);
- latexPrim[i] = latex[i] + " (" + Cbm::NumberToString(hPrim->GetEntries() / fNofEvents, 2) + "/ev.)";
- histsPrim.push_back(hPrim);
- i++;
+ double y1 = 0.17; //(pdg == 211) ? 0.20 : 0.74;
+ double y2 = 0.42; //(pdg == 211) ? 0.45 : 0.99;
+ fH.fHM.CreateCanvas(cName, cName, 900, 900);
+ DrawH1(histsAll, latexAll, kLinear, kLog, true, 0.4, y1, 0.95, y2, "hist");
+
+ fH.fHM.CreateCanvas(cName + "_prime", cName + "_prime", 900, 900);
+ DrawH1(histsPrim, latexPrim, kLinear, kLog, true, 0.4, y1, 0.95, y2, "hist");
}
- fH.fHM.CreateCanvas("lmvm_piMom", "lmvm_piMom", 900, 900);
- DrawH1(histsAll, latexAll, kLinear, kLog, true, 0.45, 0.75, 0.99, 0.99, "hist");
+ // Pions, protons and kaons that are misidentified as electrons
+ TH1* misPi = fH.H1Clone("hCandMisIdAsEl_pi");
+ TH1* misProt = fH.H1Clone("hCandMisIdAsEl_proton");
+ TH1* misKa = fH.H1Clone("hCandMisIdAsEl_kaon");
+
+ fH.fHM.CreateCanvas("AccRecMom/MisidCands", "AccRecMom/MisidCands", 800, 800);
+ DrawH1({misPi, misProt, misKa}, {"#pi^{#pm}", "proton", "kaon"}, kLinear, kLog, true, 0.89, 0.75, 0.99, 0.92, "p");
- fH.fHM.CreateCanvas("lmvm_piMomPrim", "lmvm_piMomPrim", 900, 900);
- DrawH1(histsPrim, latexPrim, kLinear, kLog, true, 0.45, 0.75, 0.99, 0.99, "hist");
+ // Acceptance in single detectors
+ for (const string& det : {"sts", "rich", "trd", "tof"}) {
+ vector<TH1*> hVec;
+ vector<string> latex;
+ for (const ELmvmSrc src : {ELmvmSrc::Signal, ELmvmSrc::Pi0, ELmvmSrc::Gamma}) {
+ for (const string& pm : {"+", "-"}) {
+ hVec.push_back(fH.H1("hMomAcc" + pm + "_" + det, src));
+ latex.push_back(fH.fSrcLatex[static_cast<int>(src)] + " (e" + pm + ")");
+ }
+ }
+ fH.fHM.CreateCanvas("AccRecMom/momDetAcc_" + det, "AccRecMom/momDetAcc_" + det, 800, 800);
+ DrawH1(hVec, latex, kLinear, kLog, true, 0.90, 0.7, 0.99, 0.99, "hist");
+ DrawTextOnPad(det, 0.4, 0.9, 0.6, 0.999);
+ }
}
-void LmvmDraw::DrawBgSource2D(const string& cName, const string& hName, const vector<string>& yLabels, double scale,
+void LmvmDraw::DrawSource2D(const string& cName, const string& hName, const vector<string>& yLabels, double scale,
const string& zTitle)
{
fH.fHM.CreateCanvas((cName + "_abs").c_str(), (cName + "_abs").c_str(), 900, 600);
TH2D* habs = fH.H2Clone(hName);
habs->SetStats(false);
habs->Scale(scale);
+ habs->SetMinimum(1e-2);
habs->GetZaxis()->SetTitle(zTitle.c_str());
habs->SetMarkerSize(1.4);
DrawH2(habs, kLinear, kLinear, kLog, "text COLZ");
@@ -716,7 +894,7 @@ void LmvmDraw::DrawBgSource2D(const string& cName, const string& hName, const ve
hperc->GetZaxis()->SetTitle("[%]");
hperc->GetYaxis()->SetLabelSize(0.06);
hperc->SetMarkerColor(kBlack);
- hperc->SetMarkerSize(1.8);
+ hperc->SetMarkerSize(1.4);
DrawH2(hperc, kLinear, kLinear, kLinear, "text COLZ");
for (size_t y = 1; y <= yLabels.size(); y++) {
@@ -732,18 +910,18 @@ void LmvmDraw::DrawBgSourceTracks()
{
gStyle->SetPaintTextFormat("4.1f");
- fH.fHM.CreateCanvas("lmvm_nofBgTracks", "lmvm_nofBgTracks", 900, 900);
+ fH.fHM.CreateCanvas("bg/nofBgTracks", "bg/nofBgTracks", 900, 900);
TH1D* hbg = fH.H1Clone("hNofBgTracks");
hbg->Scale(10);
hbg->GetYaxis()->SetTitle("Tracks/event x10");
DrawH1(hbg, kLinear, kLog, "hist text0");
hbg->SetMarkerSize(2.);
- fH.fHM.CreateCanvas("lmvm_nofSignalTracks", "lmvm_nofSignalTracks", 900, 900);
+ fH.fHM.CreateCanvas("signal/nofSignalTracks", "signal/nofSignalTracks", 900, 900);
TH1D* hel = fH.H1("hNofSignalTracks");
DrawH1(hel, kLinear, kLog, "hist");
- fH.fHM.CreateCanvas("lmvm_purity", "lmvm_purity", 900, 900);
+ fH.fHM.CreateCanvas("purity", "purity", 900, 900);
TH1D* purity = new TH1D("purity", "Purity;Analysis steps;Purity", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
purity->Divide(fH.H1("hNofBgTracks"), fH.H1("hNofSignalTracks"));
DrawH1(purity, kLinear, kLog, "hist text30");
@@ -753,11 +931,11 @@ void LmvmDraw::DrawBgSourceTracks()
SetAnalysisStepAxis(hel);
SetAnalysisStepAxis(purity);
- DrawBgSource2D("lmvm_bgSrcTracks_2d", "hBgSrcTracks",
- {"#gamma", "#pi^{0}", "#pi^{#pm}", "p", "K", "e^{#pm}_{sec}", "oth."}, 100.,
- "Tracks per event x10^{2}");
+ DrawSource2D("bg/SrcTracksBg_2d", "hBgSrcTracks",
+ {"#gamma", "#pi^{0}", "#pi^{#pm}", "p", "K", "e^{#pm}_{sec}", "oth.", "signal"}, 1000.,
+ "Tracks per event x10^{3}");
- TCanvas* c = fH.fHM.CreateCanvas("lmvm_nofTopoPairs", "lmvm_nofTopoPairs", 1600, 800);
+ TCanvas* c = fH.fHM.CreateCanvas("nofTopoPairs", "nofTopoPairs", 1600, 800);
c->Divide(2, 1);
int i = 1;
for (const string& p : {"gamma", "pi0"}) {
@@ -772,7 +950,7 @@ void LmvmDraw::DrawBgSourceTracks()
void LmvmDraw::DrawMismatchesAndGhosts()
{
gStyle->SetPaintTextFormat("4.1f");
- TCanvas* c1 = fH.fHM.CreateCanvas("lmvm_nofMismatches", "lmvm_nofMismatches", 1500, 1500);
+ TCanvas* c1 = fH.fHM.CreateCanvas("nofMismatches", "nofMismatches", 1500, 1500);
c1->Divide(2, 2);
vector<string> dets {"all", "rich", "trd", "tof"};
for (size_t i = 0; i < dets.size(); i++) {
@@ -785,7 +963,7 @@ void LmvmDraw::DrawMismatchesAndGhosts()
SetAnalysisStepAxis(h);
}
- fH.fHM.CreateCanvas("lmvm_nofGhosts", "lmvm_nofGhosts", 900, 900);
+ fH.fHM.CreateCanvas("nofGhosts", "nofGhosts", 900, 900);
DrawH1(fH.H1("hNofGhosts"), kLinear, kLog, "hist");
SetAnalysisStepAxis(fH.H1("hNofGhosts"));
}
@@ -821,7 +999,7 @@ void LmvmDraw::SetAnalysisStepAxis(TH1* h)
void LmvmDraw::DrawMvdCutQa()
{
if (!fUseMvd) return;
- TCanvas* c = fH.fHM.CreateCanvas("lmvm_cuts/lmvm_mvdCutQa", "lmvm_cuts/lmvm_mvd1cut_qa", 1600, 800);
+ TCanvas* c = fH.fHM.CreateCanvas("cuts/mvdCutQa", "cuts/mvd1cut_qa", 1600, 800);
c->Divide(2, 1);
int i = 1;
for (const string& num : {"1", "2"}) {
@@ -839,7 +1017,7 @@ void LmvmDraw::DrawMvdCutQa()
void LmvmDraw::DrawMvdAndStsHist()
{
if (!fUseMvd) return;
- TCanvas* c1 = fH.fHM.CreateCanvas("lmvm_nofHitsMvdSts", "lmvm_nofHitsMvdSts", 1600, 800);
+ TCanvas* c1 = fH.fHM.CreateCanvas("nofHitsMvdSts", "nofHitsMvdSts", 1600, 800);
c1->Divide(2, 1);
c1->cd(1);
DrawSrcH1("hNofMvdHits");
@@ -847,11 +1025,11 @@ void LmvmDraw::DrawMvdAndStsHist()
DrawSrcH1("hNofStsHits");
Draw2DCut("hMvdXY_1");
- fH.fHM.CreateCanvas("lmvm_mvd1", "lmvm_mvd1", 900, 900);
+ fH.fHM.CreateCanvas("mvd1", "mvd1", 900, 900);
DrawSrcH1("hMvdR_1");
Draw2DCut("hMvdXY_2");
- fH.fHM.CreateCanvas("lmvm_mvd2", "lmvm_mvd2", 900, 900);
+ fH.fHM.CreateCanvas("mvd2", "mvd2", 900, 900);
DrawSrcH1("hMvdR_2");
}
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h
index 1ac829f5b5..1ed6d24b2e 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h
@@ -81,7 +81,7 @@ private:
void DrawAnaStepMany(const std::string& cName, std::function<void(ELmvmAnaStep)> drawFunc);
- void DrawPtY(ELmvmAnaStep step);
+ void DrawPtY(const std::string& hist, ELmvmAnaStep step);
void DrawRapidity(ELmvmAnaStep step);
void DrawPtYEfficiency(ELmvmAnaStep step);
void DrawMinvSBg(ELmvmAnaStep step);
@@ -97,14 +97,14 @@ private:
void Draw1DCut(const std::string& hist, const std::string& sigOption, double cut = -999999.);
void Draw2DCut(const std::string& hist, double cutCrossX = -999999., double cutCrossY = -999999.);
void DrawTofM2Cut();
+ void DrawElPurity();
void DrawCuts();
+ void DrawCombinatorialPairs();
- void DrawSrcBgPairs(ELmvmAnaStep step, bool inPercent, bool drawAnaStep = true);
- void DrawSrcBgPairsAll();
+ void DrawBgSourcePairs(ELmvmAnaStep step, bool inPercent, bool drawAnaStep = true);
+ void DrawBgSourcePairsAll();
- void DrawMomAccEpEm();
-
- void DrawPiMom();
+ void DrawAccRecMom();
void Draw2DCutTriangle(double xCr, double yCr);
@@ -114,7 +114,9 @@ private:
void DrawMinvAll();
- void DrawBgSource2D(const std::string& cName, const std::string& hName, const std::vector<std::string>& yLabels,
+ void DrawMinvBg();
+
+ void DrawSource2D(const std::string& cName, const std::string& hName, const std::vector<std::string>& yLabels,
double scale, const std::string& zTitle);
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx
index 747228bcea..ef171b6af7 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx
@@ -70,6 +70,7 @@ void LmvmDrawAll::DrawHistFromFile(const string& fileInmed, const string& fileQg
DrawMinvAll();
DrawMinvCombSignalAndBg();
DrawMinvPtAll();
+ DrawMinvBgSourcesAll();
DrawSBgVsMinv();
SaveHist();
SaveCanvasToImage();
@@ -120,6 +121,16 @@ void LmvmDrawAll::CreateMeanHistAll()
}
}
}
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_gg", nofEvents, fRebMinv); // TODO: automatize this
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_pipi", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_pi0pi0", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_oo", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_gpi", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_gpi0", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_go", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_pipi0", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_pio", nofEvents, fRebMinv);
+ CreateMeanHist<TH1D>("hMinvBgSource2_elid_pi0o", nofEvents, fRebMinv);
}
TH1D* LmvmDrawAll::GetCocktailMinvH1(ELmvmAnaStep step) { return GetCocktailMinv<TH1D>("hMinv", step); }
@@ -168,6 +179,8 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
TH1D* cocktail = GetCocktailMinvH1(step);
TH1D* sbg = static_cast<TH1D*>(bg->Clone());
sbg->Add(cocktail);
+ TH1D* cb = fHMean.H1Clone("hMinvCombBgAsm", step);
+ TH1D* cbSig = fHMean.H1Clone("hMinvSBgSignalAsm", step);
double binWidth = bg->GetBinWidth(1);
vector<TH1D*> sHists(fHMean.fNofSignals);
@@ -184,10 +197,11 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
drawData.emplace_back(sbg, kBlack, kBlack, 1, -1, "Cocktail+BG");
drawData.emplace_back(bg, kGray, kBlack, 1, -1, "Background");
}
- drawData.emplace_back(sHists[static_cast<int>(ELmvmSignal::OmegaD)], kCyan + 2, kCyan + 4, 2, -1,
- "#omega #rightarrow #pi^{0}e^{+}e^{-}");
+
drawData.emplace_back(pi0, kGreen - 3, kGreen + 3, 2, -1, "#pi^{0} #rightarrow #gammae^{+}e^{-}");
drawData.emplace_back(eta, kRed - 4, kRed + 2, 2, -1, "#eta #rightarrow #gammae^{+}e^{-}");
+ drawData.emplace_back(sHists[static_cast<int>(ELmvmSignal::OmegaD)], kCyan + 2, kCyan + 4, 2, -1,
+ "#omega #rightarrow #pi^{0}e^{+}e^{-}");
drawData.emplace_back(sHists[static_cast<int>(ELmvmSignal::Omega)], kOrange + 7, kOrange + 4, 2, -1,
"#omega #rightarrow e^{+}e^{-}");
drawData.emplace_back(sHists[static_cast<int>(ELmvmSignal::Phi)], kAzure + 2, kAzure + 3, 2, -1,
@@ -196,12 +210,14 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
drawData.emplace_back(sHists[static_cast<int>(ELmvmSignal::Inmed)], kMagenta - 3, kMagenta - 2, 4, 3018,
"in-medium #rho");
drawData.emplace_back(cocktail, -1, kRed + 2, 4, -1, "Cocktail");
+ drawData.emplace_back(cb, -1, kCyan + 2, 4, -1, "comb. BG (CB)");
+ drawData.emplace_back(cbSig, -1, kCyan, 3, -1, "signal (CB subtr.)");
double min = std::numeric_limits<Double_t>::max();
double max = std::numeric_limits<Double_t>::min();
TH1D* h0 = nullptr;
- TLegend* leg = new TLegend(0.7, 0.6, 0.99, 0.99);
+ TLegend* leg = new TLegend(0.7, 0.57, 0.99, 0.97);
for (size_t i = 0; i < drawData.size(); i++) {
const auto& d = drawData[i];
d.fH->GetYaxis()->SetTitle("dN/dM_{ee} [GeV/c^{2}]^{-1}");
@@ -212,8 +228,11 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
leg->AddEntry(d.fH, d.fLegend.c_str(), "f");
DrawH1(d.fH, kLinear, kLinear, (h0 == nullptr) ? "hist" : "hist,same", d.fLineColor, d.fLineWidth, 0);
if (h0 == nullptr) h0 = d.fH;
- min = std::min(d.fH->GetMinimum(), min);
+ //min = std::min(d.fH->GetMinimum(), min);
+ min = 1e-8;
max = std::max(d.fH->GetMaximum(), max);
+ //h0->SetMinimum(1e-7);
+ //h0->SetMaximum(50);
}
if (min == 0.) min = std::min(1e-4, max * 1e-7);
if (h0 != nullptr) h0->SetMinimum(min);
@@ -225,6 +244,96 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
//DrawH1(hists, legendStr, kLinear, kLog, true, 0.7, 0.6, 0.99, 0.99, "HIST L");
}
+void LmvmDrawAll::DrawMinvBgSourcesAll()
+{
+ TH1D* gg = fHMean.H1Clone("hMinvBgSource2_elid_gg");
+ TH1D* pipi = fHMean.H1Clone("hMinvBgSource2_elid_pipi");
+ TH1D* pi0pi0 = fHMean.H1Clone("hMinvBgSource2_elid_pi0pi0");
+ TH1D* oo = fHMean.H1Clone("hMinvBgSource2_elid_oo");
+ TH1D* gpi = fHMean.H1Clone("hMinvBgSource2_elid_gpi");
+ TH1D* gpi0 = fHMean.H1Clone("hMinvBgSource2_elid_gpi0");
+ TH1D* go = fHMean.H1Clone("hMinvBgSource2_elid_go");
+ TH1D* pipi0 = fHMean.H1Clone("hMinvBgSource2_elid_pipi0");
+ TH1D* pio = fHMean.H1Clone("hMinvBgSource2_elid_pio");
+ TH1D* pi0o = fHMean.H1Clone("hMinvBgSource2_elid_pi0o");
+
+ TH1D* physBg = fHMean.H1Clone("hMinvBgSource2_elid_gg");
+ physBg->Add(fHMean.H1("hMinvBgSource2_elid_gpi0"));
+ physBg->Add(fHMean.H1("hMinvBgSource2_elid_pi0pi0"));
+ physBg->GetYaxis()->SetTitle("dN/dM_{ee} [GeV/c^{2}]^{-1}");
+
+ TH1D* cPi = fHMean.H1Clone("hMinvBgSource2_elid_pipi");
+ cPi->Add(fHMean.H1Clone("hMinvBgSource2_elid_gpi"));
+ cPi->Add(fHMean.H1Clone("hMinvBgSource2_elid_pipi0"));
+ cPi->Add(fHMean.H1Clone("hMinvBgSource2_elid_pio"));
+
+ TH1D* rest = fHMean.H1Clone("hMinvBgSource2_elid_oo");
+ rest->Add(fHMean.H1Clone("hMinvBgSource2_elid_go"));
+ rest->Add(fHMean.H1Clone("hMinvBgSource2_elid_pi0o"));
+
+
+ //physBg->SetMinimum(1e-6);
+ //physBg->SetMaximum(1e-2);
+
+ vector<LmvmDrawMinvData> drawData;
+
+ drawData.emplace_back(gpi0, kBlack, kBlack, 1, -1, "#gamma - #pi^{0}");
+ drawData.emplace_back(pi0pi0, kGray + 1, kGray + 1, 1, -1, "#pi^{0} - #pi^{0}");
+ drawData.emplace_back(gg, kCyan + 2, kCyan + 4, 2, -1, "#gamma - #gamma");
+ drawData.emplace_back(pipi0, kGreen -3 , kGreen + 3, 2, -1, "#pi^{#pm}_{misid} - #pi^{0}");
+ drawData.emplace_back(pi0o, kRed - 4, kRed + 2, 1, -1, "#pi^{0} - o.");
+ drawData.emplace_back(gpi, kOrange + 7, kOrange + 4, 2, -1, "#gamma - #pi^{#pm}_{misid}");
+ drawData.emplace_back(go, kAzure + 2, kAzure + 3, 2, -1, "#gamma - o.");
+ drawData.emplace_back(pipi, kOrange - 2, kOrange - 3, 4, 3112, "#pi^{#pm}_{misid} - #pi^{#pm}_{misid}");
+ drawData.emplace_back(pio, kMagenta - 3, kMagenta - 2, 4, 3018, "#pi^{#pm}_{misid} - o.");
+ drawData.emplace_back(oo, -1, kRed + 2, 4, -1, "o. - o.");
+
+ string cName = "minvBgSrc/minvBgSrc";
+ fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1000, 1000);
+
+ TH1D* h0 = nullptr;
+ TLegend* leg = new TLegend(0.65, 0.55, 0.95, 0.95);
+ for (size_t i = 0; i < drawData.size(); i++) {
+ const auto& d = drawData[i];
+ d.fH->GetYaxis()->SetTitle("dN/dM_{ee} [GeV/c^{2}]^{-1}");
+ d.fH->GetYaxis()->SetLabelSize(0.05);
+
+ if (d.fFillColor != -1) d.fH->SetFillColor(d.fFillColor);
+ if (d.fFillStyle != -1) d.fH->SetFillStyle(d.fFillStyle);
+ leg->AddEntry(d.fH, d.fLegend.c_str(), "f");
+ DrawH1(d.fH, kLinear, kLinear, (h0 == nullptr) ? "hist" : "hist,same", d.fLineColor, d.fLineWidth, 0);
+ if (h0 == nullptr) h0 = d.fH;
+ }
+
+ leg->SetFillColor(kWhite);
+ leg->Draw();
+ gPad->SetLogy(true);
+
+ fHMean.fHM.CreateCanvas((cName + "_misidPi").c_str(), (cName + "_misidPi").c_str(), 1000, 1000);
+ DrawH1({physBg, cPi, rest}, {"Phys. BG", "BG w. misid. #pi^{#pm}", "Rest"}, kLinear, kLog, true, 0.7, 0.8, 0.95, 0.91, "p");
+
+
+ //check if all bg pair combinations are considered
+ TH1D* bgRat = static_cast<TH1D*>(physBg->Clone());
+ bgRat->Add(cPi);
+ bgRat->Add(rest);
+
+ int nBins = bgRat->GetNbinsX();
+ for (int i = 1; i <= nBins; i++) {
+ if (fHMean.H1("hMinv", ELmvmSrc::Bg, ELmvmAnaStep::ElId)->GetBinContent(i) == 0 && bgRat->GetBinContent(i) == 0) bgRat->SetBinContent(i, 1);
+ else {
+ double r = bgRat->GetBinContent(i) / fHMean.H1("hMinv", ELmvmSrc::Bg, ELmvmAnaStep::ElId)->GetBinContent(i);
+ bgRat->SetBinContent(i, r);
+ }
+ }
+
+ bgRat->SetMinimum(0.9);
+ bgRat->SetMaximum(1.1);
+
+ fHMean.fHM.CreateCanvas((cName + "_compWithBg").c_str(), (cName + "_compWithBg").c_str(), 800, 800);
+ DrawH1(bgRat, kLinear, kLinear, "hist");
+}
+
void LmvmDrawAll::DrawMinvPtAll()
{
for (const ELmvmAnaStep step : fHMean.fAnaSteps) {
@@ -243,7 +352,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
// draw MM/PP ratio for same events
{
- TCanvas* c = fHMean.fHM.CreateCanvas("lmvmAll_CBsameEv_RatioMMPP", "lmvmAll_CBsameEv_RatioMMPP", 1800, 1800);
+ TCanvas* c = fHMean.fHM.CreateCanvas("cb/minv_RatioMMPP_same", "cb/minv_RatioMMPP_same", 1800, 1800);
c->Divide(3, 3);
int i = 1;
for (auto step : fHMean.fAnaSteps) {
@@ -261,7 +370,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
// Draw PM, MM, PP in one plot
{
for (const string& ev : {"sameEv", "mixedEv"}) {
- string cName = "lmvmAll_CB" + ev;
+ string cName = "cb/minv_pairYield_" + ev;
TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800);
c->Divide(3, 3);
int i = 1;
@@ -281,7 +390,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
// Draw same and mixed for PP, MM, PM cases; normalize to 400 - 700 MeV/c2 interval
{
for (const string& comb : {"PM", "PP", "MM"}) {
- string cName = "lmvmAll_CB_sameMixed" + comb;
+ string cName = "cb/minv_pairYield_" + comb + "_sameMixed";
TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800);
c->Divide(3, 3);
int i = 1;
@@ -294,12 +403,10 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
int maxBin = same->FindBin(0.7);
double scale = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin);
mixed->Scale(scale);
- DrawH1({same, mixed}, {"Same " + comb, "Mixed " + comb}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99, "HIST");
+ same->SetMinimum(2e-7);
+ mixed->SetMinimum(2e-7);
+ DrawH1({same, mixed}, {"Same " + comb, "Mixed " + comb}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99, "p");
fHMean.DrawAnaStepOnPad(step);
- int nofSamePairs = same->GetEntries();
- int nofMixedPairs = mixed->GetEntries();
- DrawTextOnPad("same: " + Cbm::NumberToString(nofSamePairs, 5), 0.4, 0.7, 0.8, 0.9);
- DrawTextOnPad("mixed: " + Cbm::NumberToString(nofMixedPairs, 5), 0.4, 0.5, 0.8, 0.7);
}
}
}
@@ -307,7 +414,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
// Draw ratio mixed/same for PP, MM, PM cases
{
for (const string& comb : {"PM", "PP", "MM"}) {
- string cName = "lmvmAll_CB_mixedOverSame" + comb;
+ string cName = "cb/minv_mixedOverSame" + comb;
TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800);
c->Divide(3, 3);
int i = 1;
@@ -323,6 +430,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
}
}
+
// compare PM with sbg (Cocktail + BG)
{
for (auto step : fHMean.fAnaSteps) {
@@ -332,7 +440,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
sbg->Add(GetCocktailMinvH1(step));
TH1D* ratio = (TH1D*) pmSame->Clone();
- string cName = "lmvmAll_CB_SameVsSbg_" + fHMean.fAnaStepNames[static_cast<int>(step)];
+ string cName = "cb/minv_N+-VsCocBg_" + fHMean.fAnaStepNames[static_cast<int>(step)];
TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 900);
c->Divide(2, 1);
c->cd(1);
@@ -348,7 +456,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
//Draw Geometric Mean
{
- TCanvas* c = fHMean.fHM.CreateCanvas("lmvmAll_minvCombGeom", "lmvmAll_minvCombGeom", 1800, 1800);
+ TCanvas* c = fHMean.fHM.CreateCanvas("cb/minv_geomMean", "cb/minv_geomMean", 1800, 1800);
c->Divide(3, 3);
int i = 1;
for (auto step : fHMean.fAnaSteps) {
@@ -361,80 +469,47 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
double scale = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin);
mixed->Scale(scale);
same->GetXaxis()->SetTitle("M_{ee}");
- DrawH1({same, mixed}, {"same", "mixed"}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99, "HIST");
+ same->SetMinimum(1e-8);
+ mixed->SetMinimum(1e-8);
+ DrawH1({same, mixed}, {"same", "mixed"}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99, "p");
fHMean.DrawAnaStepOnPad(step);
}
}
- // Draw CombBg vs Cocktail+Bg vs combBg(Asm)
+ // Draw k factor (deleted old draw procedure)
+ //{
+ //}
+
+
+ // Draw CombBg vs BG from MC input
{
for (auto step : fHMean.fAnaSteps) {
if (step < ELmvmAnaStep::ElId) continue;
- TH1D* cbg = fHMean.H1Clone("hMinvCombBg", step);
- TH1D* cbgAsm = fHMean.H1Clone("hMinvCombBgAsm", step);
- TH1D* sbg = GetCocktailMinvH1(step);
- sbg->Add(fHMean.H1("hMinv", ELmvmSrc::Bg, step));
-
- TH1D* ratio1 = static_cast<TH1D*>(sbg->Clone());
- TH1D* ratio2 = static_cast<TH1D*>(sbg->Clone());
-
- string cName = "lmvmAll_minvCombBgVsSBgVsCombBgAsm_" + fHMean.fAnaStepNames[static_cast<int>(step)];
+ TH1D* cb = fHMean.H1Clone("hMinvCombBgAsm", step);
+ TH1D* bg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step);
+
+ TH1D* ratio = static_cast<TH1D*>(bg->Clone());
+ ratio->Divide(cb);
+
+ string cName = "cb/minv_bgVsCb_" + fHMean.fAnaStepNames[static_cast<int>(step)];
TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 900);
c->Divide(2, 1);
c->cd(1);
- DrawH1({cbg, sbg, cbgAsm}, {"CombBG", "Cocktail + BG", "CombBG(Asm)"}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99,
- "HIST");
+ DrawH1({bg, cb}, {"BG from MC input", "comb. BG (B_{C})"}, kLinear, kLog, true, 0.7, 0.75, 0.99, 0.90,
+ "hist");
c->cd(2);
- ratio1->Divide(cbg);
- ratio2->Divide(cbgAsm);
- DrawH1({ratio1, ratio2}, {"Cocktail+BG/CombBG", "Cocktail+BG/CombBG(Asm)"}, kLinear, kLog, true, 0.65, 0.85, 0.99,
- 0.99, "HIST");
+ DrawH1({ratio}, {"BG_{MC}/B_{C}"}, kLinear, kLinear, true, 0.65, 0.75, 0.92,
+ 0.9, "hist");
fHMean.DrawAnaStepOnPad(step);
}
}
- //Draw k factor
- {
- fHMean.fHM.CreateCanvas("lmvmAll_minvCombK", "lmvmAll_minvCombK", 1000, 1000);
- vector<string> legend;
- vector<TH1*> hists;
- for (auto step : fHMean.fAnaSteps) {
- if (step < ELmvmAnaStep::Reco) continue;
- hists.push_back(fHMean.H1Clone("hMinvCombK", step));
- legend.push_back(fHMean.fAnaStepLatex[static_cast<int>(step)]);
- }
- hists[0]->GetXaxis()->SetTitle("M_{ee}");
- hists[0]->GetYaxis()->SetTitle("k factor");
- DrawH1(hists, legend, kLinear, kLinear, true, 0.8, 0.6, 0.99, 0.99, "HIST");
- }
-
- //Draw Combinatorial Signal from N+-same (without mixed geom. mean) // TODO: without mixed: needed or delete?
- {
- TCanvas* c = fHMean.fHM.CreateCanvas("lmvmAll_minvCombPMSignal", "lmvmAll_minvCombPMSignal", 1800, 1800);
- c->Divide(3, 3);
- int i = 1;
- for (auto step : fHMean.fAnaSteps) {
- if (step < ELmvmAnaStep::ElId) continue;
- c->cd(i++);
- TH1D* pmSame = fHMean.H1Clone("hMinvCombPM_sameEv", step);
- TH1D* combPMSignalAsm = fHMean.H1Clone("hMinvCombPMSignal", step);
- TH1D* combBgAsm = fHMean.H1Clone("hMinvCombBg", step);
- TH1D* cocktail = GetCocktailMinvH1(step);
- pmSame->SetMaximum(3e-1);
- pmSame->SetMinimum(4e-8);
- DrawH1({pmSame, combBgAsm, combPMSignalAsm, cocktail},
- {"N_{same}^{+-}", "B_{c}", "Signal (N_{same}^{+-} - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.8, 0.8,
- 0.99, 0.99, "hist");
- fHMean.DrawAnaStepOnPad(step);
- }
- }
-
//Draw Combinatorial Signal from N+-same
{
- TCanvas* c = fHMean.fHM.CreateCanvas("lmvmAll_minvCombPMSignalAsm", "lmvmAll_minvCombPMSignalAsm", 1800, 1800);
+ TCanvas* c = fHMean.fHM.CreateCanvas("cb/minv_cbVsInput_PM", "cb/minv_cbVsInput_PM", 1800, 1800);
c->Divide(3, 3);
int i = 1;
for (auto step : fHMean.fAnaSteps) {
@@ -453,52 +528,48 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
}
}
+
//Draw Combinatorial Signal from Cocktail+BG
{
- TCanvas* c = fHMean.fHM.CreateCanvas("lmvmAll_minvSbgSignalAsm", "lmvmAll_minvSbgSignalAsm", 1800, 1800);
+ TCanvas* c = fHMean.fHM.CreateCanvas("cb/minv_cbVsInput_cocBg", "cb/minv_cbVsInput_cocBg", 1800, 1800);
c->Divide(3, 3);
int i = 1;
for (auto step : fHMean.fAnaSteps) {
if (step < ELmvmAnaStep::ElId) continue;
c->cd(i++);
- TH1D* sbg = fHMean.H1Clone("hMinv_bg", step);
- sbg->Add(GetCocktailMinvH1(step));
+ TH1D* sbg2 = fHMean.H1Clone("hMinv_bg", step);
+ sbg2->Add(GetCocktailMinvH1(step));
TH1D* combBgAsm = fHMean.H1Clone("hMinvCombBgAsm", step);
TH1D* sBgSignalAsm = fHMean.H1Clone("hMinvSBgSignalAsm", step);
- sbg->SetMaximum(3e-1);
- sbg->SetMinimum(4e-8);
+ sbg2->GetYaxis()->SetTitle("dN/dM_{ee} [GeV/c^{2}]^{-1}");
+ sbg2->SetMaximum(1e-2);
+ sbg2->SetMinimum(4e-9);
TH1D* cocktail = GetCocktailMinvH1(step);
- DrawH1({sbg, combBgAsm, sBgSignalAsm, cocktail},
- {"Cocktail + BG", "B_{c} (asm)", "Signal (Cocktail+BG - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.8,
- 0.8, 0.99, 0.99, "hist");
+ DrawH1({sbg2, combBgAsm, sBgSignalAsm, cocktail},
+ {"Cocktail + BG", "B_{C}", "Signal (Cock+BG - B_{C})", "Cocktail"}, kLinear, kLog, true, 0.53,
+ 0.72, 0.99, 0.92, "p");
fHMean.DrawAnaStepOnPad(step);
}
}
-
- //Draw Combinatorial Signal from Cocktail+BG (without mixed geom. mean) // TODO: without mixed: needed or delete?
+ // Draw ratio of CB-signal and cocktail
{
- TCanvas* c = fHMean.fHM.CreateCanvas("lmvmAll_minvSbgSignal", "lmvmAll_minvSbgSignal", 1800, 1800);
- c->Divide(3, 3);
- int i = 1;
- for (auto step : fHMean.fAnaSteps) {
- if (step < ELmvmAnaStep::ElId) continue;
- c->cd(i++);
- TH1D* sbg = fHMean.H1Clone("hMinv_bg", step);
- sbg->Add(GetCocktailMinvH1(step));
- TH1D* combBgAsm = fHMean.H1Clone("hMinvCombBg", step);
- TH1D* sBgSignalAsm = fHMean.H1Clone("hMinvSBgSignal", step);
- sbg->SetMaximum(3e-1);
- sbg->SetMinimum(4e-8);
- TH1D* cocktail = GetCocktailMinvH1(step);
- DrawH1({sbg, combBgAsm, sBgSignalAsm, cocktail},
- {"Cocktail + BG", "B_{c}", "Signal (Cocktail+BG - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.8, 0.8,
- 0.99, 0.99, "hist");
- fHMean.DrawAnaStepOnPad(step);
- }
+ ELmvmAnaStep step = ELmvmAnaStep::TtCut;
+
+ fHMean.fHM.CreateCanvas("cb/minv_ratio_cock-CBSignal", "cb/minv_ratio_cock-CBSignal", 1800, 1800);
+
+ TH1D* sig = fHMean.H1("hMinvSBgSignalAsm", step);
+ TH1D* rat = GetCocktailMinvH1(step);
+ rat->Divide(sig);
+
+ rat->SetMaximum(3.);
+
+ DrawH1({rat}, {"Cocktail/S_{CB}"}, kLinear, kLinear, true, 0.5, 0.8, 0.85, 0.86, "hist");
+ fHMean.DrawAnaStepOnPad(step);
}
}
+
void LmvmDrawAll::DrawSBgVsMinv()
{
vector<TH1*> hists;
@@ -600,18 +671,22 @@ void LmvmDrawAll::CalcCombBGHistos()
hBCSignalAsm->Add(GetCocktailMinvH1(step));
hBCSignalAsm->Add(fHMean.H1("hMinvCombBgAsm", step), -1.);
- // TODO: @Cornelius check if calculations are correct
+ // calculate errors by error propagation formula
+ // TODO: @Cornelius check if calculations are correct
+
+ int nofEvents = GetNofTotalEvents();
+ double bW = fHMean.H1("hMinvCombPM_sameEv", step)->GetBinWidth(1);
int nofBins = fHMean.H1("hMinvCombPM_sameEv", step)->GetNbinsX();
for (int i = 1; i <= nofBins; i++) {
//s_ for same, m_ for mixed
- double s_pm = fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(i);
- double s_pp = fHMean.H1("hMinvCombPP_sameEv", step)->GetBinContent(i);
- double s_mm = fHMean.H1("hMinvCombMM_sameEv", step)->GetBinContent(i);
- double m_pm = fHMean.H1("hMinvCombPM_mixedEv", step)->GetBinContent(i);
- double m_pp = fHMean.H1("hMinvCombPP_mixedEv", step)->GetBinContent(i);
- double m_mm = fHMean.H1("hMinvCombMM_mixedEv", step)->GetBinContent(i);
-
- double s_dpm = 1.; // derivatives of B_c and signal resp. to single contributions
+ double s_pm = fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(i) * nofEvents * bW;
+ double s_pp = fHMean.H1("hMinvCombPP_sameEv", step)->GetBinContent(i) * nofEvents * bW;
+ double s_mm = fHMean.H1("hMinvCombMM_sameEv", step)->GetBinContent(i) * nofEvents * bW;
+ double m_pm = fHMean.H1("hMinvCombPM_mixedEv", step)->GetBinContent(i) * nofEvents * bW;
+ double m_pp = fHMean.H1("hMinvCombPP_mixedEv", step)->GetBinContent(i) * nofEvents * bW;
+ double m_mm = fHMean.H1("hMinvCombMM_mixedEv", step)->GetBinContent(i) * nofEvents * bW;
+
+ /*double s_dpm = 1.; // derivatives of B_c and signal resp. to single contributions
double d1 = 1. / std::sqrt(s_pp * s_mm * m_pp * m_mm);
double s_dpp = 0.5 * m_pm * s_mm * d1;
double s_dmm = 0.5 * m_pm * s_pp * d1;
@@ -630,10 +705,27 @@ void LmvmDrawAll::CalcCombBGHistos()
double m_fmm = std::pow(std::sqrt(m_mm) * m_dmm, 2); //TODO: @Cornelius check, changed from s_mm2 to m_mm2
double m_fpm2 = std::pow(std::sqrt(m_pm) * m_dpm2, 2); // for > 300 MeV
- double errorBc = std::sqrt(s_fpp + s_fmm + m_fpm + m_fpp + m_fmm); // final error propagation value
- double errorBc2 = std::sqrt(m_fpm2);
- double errorSig = std::sqrt(s_fpm + s_fpp + s_fmm + m_fpm + m_fpp + m_fmm);
- double errorSig2 = std::sqrt(s_fpm + m_fpm2);
+ double errorBc = std::sqrt(s_fpp + s_fmm + m_fpm + m_fpp + m_fmm) / (nofEvents * bW); // final error propagation value
+ double errorBc2 = std::sqrt(m_fpm2) / (nofEvents * bW);
+ double errorSig = std::sqrt(s_fpm + s_fpp + s_fmm + m_fpm + m_fpp + m_fmm) / (nofEvents * bW);
+ double errorSig2 = std::sqrt(s_fpm + m_fpm2) / (nofEvents * bW);*/
+
+ double d_m_pm = std::sqrt(s_pp * s_mm / (m_pp * m_mm));
+ double d_m_pp = -0.5 * m_pm * m_mm * std::sqrt(s_pp * s_mm) / (std::pow(std::sqrt(m_pp * m_mm), 3));
+ double d_m_mm = -0.5 * m_pm * m_pp * std::sqrt(s_pp * s_mm) / (std::pow(std::sqrt(m_pp * m_mm), 3));
+ double d_s_pp = 0.5 * m_pm * s_mm / std::sqrt(m_pp * m_mm * s_pp * s_mm);
+ double d_s_mm = 0.5 * m_pm * s_pp / std::sqrt(m_pp * m_mm * s_pp * s_mm);
+
+ double f_m_pm = std::pow(d_m_pm * std::sqrt(m_pm), 2);
+ double f_m_pp = std::pow(d_m_pp * std::sqrt(m_pp), 2);
+ double f_m_mm = std::pow(d_m_mm * std::sqrt(m_mm), 2);
+ double f_s_pp = std::pow(d_s_pp * std::sqrt(s_pp), 2);
+ double f_s_mm = std::pow(d_s_mm * std::sqrt(s_mm), 2);
+
+ double errorBc = std::sqrt(f_m_pm + f_m_pp + f_m_mm + f_s_pp + f_s_mm) / (nofEvents * bW);
+ double errorBc2 = normGM * std::sqrt(m_pm) / (nofEvents * bW);// for range > 0.3 GeV
+ double errorSig = std::sqrt(s_pm + std::pow(errorBc, 2)) / (nofEvents * bW);
+ double errorSig2 = std::sqrt(s_pm + std::pow(errorBc2, 2)) / (nofEvents * bW);
fHMean.H1("hMinvCombBg", step)->SetBinError(i, errorBc);
if (i <= minBin) fHMean.H1("hMinvCombBgAsm", step)->SetBinError(i, errorBc);
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h
index 4a2871b691..9d138f1c30 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h
@@ -37,7 +37,7 @@ private:
std::vector<LmvmHist*> fH;
LmvmHist fHMean;
- int fRebMinv = 40; // Rebin for minv histograms
+ int fRebMinv = 50; // Rebin for minv histograms
std::string fOutputDir; // output directory for figures
@@ -58,6 +58,7 @@ private:
void DrawMinvAll();
void DrawMinv(ELmvmAnaStep step);
void DrawMinvPtAll();
+ void DrawMinvBgSourcesAll();
/**
* \brief Draw invariant mass spectra for all signal types for specified analysis step with BG reduced by combinatorial BG.
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h b/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h
index 414830101a..7b8f55935f 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h
@@ -23,7 +23,7 @@ public:
if (particle == "omegadalitz" || particle == "wdalitz") return 2.28721 * 7.7e-4;
if (particle == "phi") return 0.311619 * 2.97e-4;
if (particle == "inmed" || particle == "rho0") return 0.0304706;
- if (particle == "qgp" || particle == "qgp_epem") return 4.52941e-4;
+ if (particle == "qgp" || particle == "qgp_epem") return 10 * 4.52941e-4; // TODO urgent: delete factor 10 and set right parameter!!!
}
else if (energy == "25gev") {
if (particle == "rho0") return 23 * 4.7e-5;
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx
index 3e09d218f3..f67f6a6183 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx
@@ -64,6 +64,10 @@ void LmvmTask::InitHists()
"Counter", 100, 0., 5., 1000, 0., 50.);
fH.CreateH1("hMotherPdg", {"mc", "acc"}, "Pdg code", "Particles/event", 7000, -3500., 3500.);
+ fH.CreateH1("hCandPdg", fH.fAnaStepNames, "Pdg code", "Particles/event", 7001, -3500., 3500.);
+ fH.CreateH2("hCandPdgVsMom", fH.fAnaStepNames, "P [GeV/c]", "Particle ID", "Yield/(Event * Bin)", 120, 0., 6., 6, 0., 6.);
+ fH.CreateH2("hCandElSrc", "Analysis step", "Mother of Electron Candidate", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.);
+ fH.CreateH2("hBgPairPdg", fH.fAnaStepNames, "PDG of Candidate 1", "PDG of Candidate 2", ax, 8, 0., 8., 8, 0., 8.);
fH.CreateH2("hPmtXY", fH.fSrcNames, "X [cm]", "Y [cm]", "Counter", 110, -110, 110, 200, -200, 200);
@@ -75,7 +79,7 @@ void LmvmTask::InitHists()
fH.CreateH1("hNofBgTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
fH.CreateH1("hNofSignalTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
- fH.CreateH2("hBgSrcTracks", "Analysis step", "Particle", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 7, 0., 7.);
+ fH.CreateH2("hBgSrcTracks", "Analysis step", "Candidate Source", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.);
fH.CreateH1("hNofTopoPairs", {"gamma", "pi0"}, "Pair type", "Pairs/event", 8, 0., 8);
fH.CreateH1("hNofMismatches", {"all", "rich", "trd", "tof"}, "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0.,
@@ -85,8 +89,8 @@ void LmvmTask::InitHists()
fH.CreateH2("hSrcBgPairs", "Analysis step", "Pair", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, fH.fNofBgPairSrc, 0.,
fH.fNofBgPairSrc);
- fH.CreateH2("hSrcBgPairsEpEm", fH.fAnaStepNames, "mother particle e+", "mother particle e-", ax, 3, 0., 3., 3, 0.,
- 3.);
+ fH.CreateH2("hSrcBgPairsEpEm", fH.fAnaStepNames, "mother particle e+", "mother particle e-", ax, 4, 0., 4., 4, 0.,
+ 4.);
fH.CreateH1("hEventNumber", "", "", 1, 0, 1.);
fH.CreateH1("hEventNumberMixed", "", "", 1, 0, 1.);
@@ -94,13 +98,17 @@ void LmvmTask::InitHists()
fH.CreateH1("hAnnRich", fH.fSrcNames, "RICH ANN output", ax, 100, -1.1, 1.1);
fH.CreateH2("hAnnRichVsMom", fH.fSrcNames, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 6., 100, -1.1, 1.1);
fH.CreateH1("hAnnTrd", fH.fSrcNames, "Likelihood output", ax, 100, -.1, 1.1); // TODO: change back to "TRD ANN"
- fH.CreateH2("hTofM2", fH.fSrcNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 4., 500, -0.1, 1.0);
+ fH.CreateH2("hTofM2", fH.fSrcNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 150, 0., 6., 500, -0.1, 1.0);
fH.CreateH1("hChi2Sts", fH.fSrcNames, "#chi^{2}", ax, 200, 0., 20.);
fH.CreateH1("hChi2PrimVertex", fH.fSrcNames, "#chi^{2}_{prim}", ax, 200, 0., 20.);
fH.CreateH1("hNofMvdHits", fH.fSrcNames, "Number of hits in MVD", ax, 5, -0.5, 4.5);
fH.CreateH1("hNofStsHits", fH.fSrcNames, "Number of hits in STS", ax, 9, -0.5, 8.5);
- fH.CreateH2("hTrdLike", {"El", "Pi"}, fH.fSrcNames, "P [GeV/c]", "Likelihood output", ax, 100, 0., 6., 100, -0.1,
- 1.1);
+ fH.CreateH2("hTrdLike", {"El", "Pi"}, fH.fSrcNames, "P [GeV/c]", "Likelihood output", ax, 100, 0., 6., 100, 0.,
+ 1.);
+
+ fH.CreateH2("hAnnRichVsMomPur", {"El", "Bg"}, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 6., 100, -1.1, 1.1);
+ fH.CreateH2("hTrdElLikePur", {"El", "Bg"}, "P [GeV/c]", "Likelihood output", ax, 100, 0., 6., 100, 0., 1.);
+ fH.CreateH2("hTofM2Pur", {"El", "Bg"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 150, 0., 6., 500, -0.1, 1.0);
fH.CreateH2("hTtCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
"#theta_{e^{+},e^{-}} [deg]", ax, 100, 0., 5., 100, 0., 5.);
@@ -117,18 +125,25 @@ void LmvmTask::InitHists()
fH.CreateH1("hMvdMcDist", {"1", "2"}, fH.fSrcNames, "Track-Hit distance [cm]", ax, 100, 0., 10.);
fH.CreateH1("hMinv", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0., 2.5);
- fH.CreateH1("hMinvCombPM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{e+e-} [GeV/c^{2}]", ax, 2500, 0., 2.5);
- fH.CreateH1("hMinvCombPP", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{e+e+} [GeV/c^{2}]", ax, 2500, 0., 2.5);
- fH.CreateH1("hMinvCombMM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{e-e-} [GeV/c^{2}]", ax, 2500, 0., 2.5);
+ fH.CreateH1("hMinvCombPM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0., 2.5);
+ fH.CreateH1("hMinvCombPP", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0., 2.5);
+ fH.CreateH1("hMinvCombMM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0., 2.5);
+
+ fH.CreateH1("hMinvCombPM", {"pluto", "urqmd"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0, 2.5); // mind: changed (18.05.22); was '"hMinvCombPM_elid", {"pluto", "urqmd"}' before
+ fH.CreateH1("hMinvCombPP", {"pluto", "urqmd"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0, 2.5);
+ fH.CreateH1("hMinvCombMM", {"pluto", "urqmd"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0, 2.5);
+
fH.CreateH1("hMinvBgMatch", {"trueMatch", "trueMatchEl", "trueMatchNotEl", "mismatch"}, fH.fAnaStepNames,
"M_{ee} [GeV/c^{2}]", ax, 2000, 0., 2.5);
fH.CreateH1("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 2500, 0., 2.5);
-
+ fH.CreateH1("hMinvBgSource2_elid", {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"}, "M_{ee} [GeV/c^{2}]", ax, 2500, 0., 2.5); // "pi" are misid. charged pions
+
fH.CreateH2("hMinvPt", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", ax, 100, 0., 2., 25, 0.,
2.5);
fH.CreateH1("hMomPairSignal", fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 15.);
fH.CreateH2("hPtYPairSignal", fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
+ fH.CreateH2("hPtYCandidate", fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
fH.CreateH1("hAnglePair", fH.fSrcNames, fH.fAnaStepNames, "#Theta_{1,2} [deg]", ax, 160, 0., 80.);
for (const string& suff : {"", "+", "-"}) {
@@ -143,8 +158,13 @@ void LmvmTask::InitHists()
fH.CreateH1("hMomAcc+", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.);
fH.CreateH1("hMomAcc-", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.);
- fH.CreateH1("hPiMom", {"all", "prim"}, {"mc", "acc", "rec", "recOnlySts", "recStsRichTrd", "recStsRichTrdTof"},
+ fH.CreateH1("hPiMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
"P [GeV/c]", ax, 30, 0., 3.);
+ fH.CreateH1("hElMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
+ "P [GeV/c]", ax, 30, 0., 3.);
+ fH.CreateH1("hProtonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
+ "P [GeV/c]", ax, 30, 0., 3.);
+ fH.CreateH1("hCandMisIdAsEl", {"pi", "pi0", "proton", "kaon"}, "P [GeV/c]", ax, 30, 0., 3.);
}
InitStatus LmvmTask::Init()
@@ -215,7 +235,7 @@ void LmvmTask::Exec(Option_t*)
CalculateNofTopologyPairs("hNofTopoPairs_pi0", ELmvmSrc::Pi0);
DifferenceSignalAndBg();
SignalAndBgReco();
- FillPionsHist();
+ FillAccRecVsMomHist();
fCandsTotal.insert(fCandsTotal.end(), fCands.begin(), fCands.end());
LOG(info) << "fCandsTotal.size = " << fCandsTotal.size();
@@ -375,61 +395,75 @@ bool LmvmTask::IsMcTrackAccepted(int mcTrackInd)
&& tr->GetNPoints(ECbmModuleId::kTrd) >= 2 && tr->GetNPoints(ECbmModuleId::kTof) > 1);
}
-void LmvmTask::FillPionsHist()
+void LmvmTask::FillAccRecVsMomHist()
{
- int nofMcTracks = fMCTracks->GetEntriesFast();
- for (int i = 0; i < nofMcTracks; i++) {
- CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i));
- bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4);
- TVector3 vertex;
- mct->GetStartVertex(vertex);
-
- if (std::abs(mct->GetPdgCode()) == 211) {
- fH.FillH1("hPiMom_all_mc", mct->GetP());
- if (isAccSts) fH.FillH1("hPiMom_all_acc", mct->GetP());
-
- if (vertex.Mag() < 0.1) {
- fH.FillH1("hPiMom_prim_mc", mct->GetP());
- if (isAccSts) fH.FillH1("hPiMom_prim_acc", mct->GetP());
+ for (const int& pdg : {11, 211, 2212} ) { // TODO: restructure this more efficiently (wihtout these pdg-loops)
+ int nofMcTracks = fMCTracks->GetEntriesFast();
+ string hName = (pdg == 11) ? "hElMom" : (pdg == 211) ? "hPiMom" : "hProtonMom";
+
+ for (int i = 0; i < nofMcTracks; i++) {
+ CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i));
+ bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4);
+ TVector3 vertex;
+ mct->GetStartVertex(vertex);
+
+ if (std::abs(mct->GetPdgCode()) == pdg) {
+ fH.FillH1(hName + "_all_mc", mct->GetP());
+ if (isAccSts) fH.FillH1(hName + "_all_acc", mct->GetP());
+
+ if (vertex.Mag() < 0.1) {
+ fH.FillH1(hName + "_prim_mc", mct->GetP());
+ if (isAccSts) fH.FillH1(hName + "_prim_acc", mct->GetP());
+ }
}
}
- }
- int ngTracks = fGlobalTracks->GetEntriesFast();
- for (int i = 0; i < ngTracks; i++) {
- CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(i));
- if (gTrack == nullptr) continue;
- int stsInd = gTrack->GetStsTrackIndex();
- bool isRich = (gTrack->GetRichRingIndex() >= 0);
- bool isTrd = (gTrack->GetTrdTrackIndex() >= 0);
- bool isTof = (gTrack->GetTofHitIndex() >= 0);
+ int ngTracks = fGlobalTracks->GetEntriesFast();
+ for (int i = 0; i < ngTracks; i++) {
+ CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(i));
+ if (gTrack == nullptr) continue;
+ int stsInd = gTrack->GetStsTrackIndex();
+ bool isRich = (gTrack->GetRichRingIndex() >= 0);
+ bool isTrd = (gTrack->GetTrdTrackIndex() >= 0);
+ bool isTof = (gTrack->GetTofHitIndex() >= 0);
- if (stsInd < 0) continue;
- CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
- if (stsTrack == nullptr) continue;
- CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
- if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
- int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
- if (stsMcTrackId < 0) continue;
- CbmMCTrack* mct = (CbmMCTrack*) fMCTracks->At(stsMcTrackId);
- if (mct == nullptr) continue;
- TVector3 vertex;
- mct->GetStartVertex(vertex);
-
- if (std::abs(mct->GetPdgCode()) == 211) {
- fH.FillH1("hPiMom_all_rec", mct->GetP());
- if (!isRich && !isTrd && !isTof) { fH.FillH1("hPiMom_all_recOnlySts", mct->GetP()); }
- if (isRich && isTrd) { fH.FillH1("hPiMom_all_recStsRichTrd", mct->GetP()); }
- if (isRich && isTrd && isTof) { fH.FillH1("hPiMom_all_recStsRichTrdTof", mct->GetP()); }
-
- if (vertex.Mag() < 0.1) {
- fH.FillH1("hPiMom_prim_rec", mct->GetP());
- if (!isRich && !isTrd && !isTof) { fH.FillH1("hPiMom_prim_recOnlySts", mct->GetP()); }
- if (isRich && isTrd) { fH.FillH1("hPiMom_prim_recStsRichTrd", mct->GetP()); }
- if (isRich && isTrd && isTof) { fH.FillH1("hPiMom_prim_recStsRichTrdTof", mct->GetP()); }
+ if (stsInd < 0) continue;
+ CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
+ if (stsTrack == nullptr) continue;
+ CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
+ if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
+ int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
+ if (stsMcTrackId < 0) continue;
+ CbmMCTrack* mct = (CbmMCTrack*) fMCTracks->At(stsMcTrackId);
+ if (mct == nullptr) continue;
+ TVector3 vertex;
+ mct->GetStartVertex(vertex);
+
+ if (std::abs(mct->GetPdgCode()) == pdg) {
+ fH.FillH1(hName + "_all_recSts", mct->GetP());
+ if (isRich) { fH.FillH1(hName + "_all_recStsRich", mct->GetP()); }
+ if (isRich && isTrd) { fH.FillH1(hName + "_all_recStsRichTrd", mct->GetP()); }
+ if (isRich && isTrd && isTof) { fH.FillH1(hName + "_all_recStsRichTrdTof", mct->GetP()); }
+
+ if (vertex.Mag() < 0.1) {
+ fH.FillH1(hName + "_prim_recSts", mct->GetP());
+ if (isRich) { fH.FillH1(hName + "_prim_recStsRich", mct->GetP()); }
+ if (isRich && isTrd) { fH.FillH1(hName + "_prim_recStsRichTrd", mct->GetP()); }
+ if (isRich && isTrd && isTof) { fH.FillH1(hName + "_prim_recStsRichTrdTof", mct->GetP()); }
+ }
}
}
}
+
+ // Fill histo with misidentified pions, protons, kaons // TODO: move this to SignalAndBgReco()?
+ for (const auto& cand : fCands) {
+ int pdg = cand.fMcPdg;
+ if (!(pdg == 211 || pdg == 111 || pdg == 2212 || pdg == 321)) continue;
+ if (!cand.fIsElectron) continue;
+ cout << "LmvmTask::FillAccRecVsMomHist(): pdg = " << pdg << endl;
+ string hName = (pdg == 211) ? "hCandMisIdAsEl_pi" : (pdg == 111) ? "hCandMisIdAsEl_pi0" : (pdg == 2212) ? "hCandMisIdAsEl_proton" : "hCandMisIdAsEl_kaon";
+ fH.FillH1(hName, cand.fMomentum.Mag());
+ }
}
@@ -580,12 +614,29 @@ void LmvmTask::CombinatorialPairs()
for (size_t iC2 = iC1 + 1; iC2 < nCand; iC2++) {
const auto& cand2 = fCandsTotal[iC2];
- if (cand1.IsMcSignal() xor cand2.IsMcSignal()) continue;
+ if (cand1.IsMcSignal() or cand2.IsMcSignal()) continue; // TODO: 'OR' or 'XOR'?
LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2);
double w = (cand1.IsMcSignal() && cand2.IsMcSignal()) ? fW : 1.;
bool isSameEvent = (cand1.fEventNumber == cand2.fEventNumber);
for (auto step : fH.fAnaSteps) {
if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue;
+
+ // seperate PLUTO and UrQMD candidates
+ if (cand1.IsCutTill(step) && cand2.IsCutTill(step) && std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11) {
+ // only PLUTO
+ if (cand1.IsMcSignal()) {
+ if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_pluto", step, pRec.fMinv, w);
+ else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_pluto", step, pRec.fMinv, w);
+ else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_pluto", step, pRec.fMinv, w);
+ }
+ // only UrQMD
+ else if (!cand1.IsMcSignal()) {
+ if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmd", step, pRec.fMinv, w);
+ else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmd", step, pRec.fMinv, w);
+ else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmd", step, pRec.fMinv, w);
+ }
+ }
+
if (cand1.IsCutTill(step) && cand2.IsCutTill(step)) {
if (cand1.fCharge * cand2.fCharge < 0) {
if (isSameEvent) fH.FillH1("hMinvCombPM_sameEv", step, pRec.fMinv, w);
@@ -674,17 +725,64 @@ void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAna
fH.FillH1("hAnglePair", src, step, parRec.fAngle, fW);
if (src == ELmvmSrc::Bg) {
- // gamma=0.5, pi0=1.5, other=2.5
- double indM = candM.IsMcGamma() ? 0.5 : (candM.IsMcPi0() ? 1.5 : 2.5);
- double indP = candP.IsMcGamma() ? 0.5 : (candP.IsMcPi0() ? 1.5 : 2.5);
+ // gamma=0.5, pi0=1.5, pions=2.5, other=3.5
+ double indM = candM.IsMcGamma() ? 0.5 : (candM.IsMcPi0() ? 1.5 : (std::abs(candM.fMcPdg) == 211) ? 2.5 : 3.5);
+ double indP = candP.IsMcGamma() ? 0.5 : (candP.IsMcPi0() ? 1.5 : (std::abs(candP.fMcPdg) == 211) ? 2.5 : 3.5);
fH.FillH2("hSrcBgPairsEpEm", step, indP, indM);
ELmvmBgPairSrc bgSrc = LmvmUtils::GetBgPairSrc(candP, candM);
- if (bgSrc != ELmvmBgPairSrc::Undefined) {
+ //if (bgSrc != ELmvmBgPairSrc::Undefined) {
string name = fH.GetName("hMinvBgSource_" + fH.fBgPairSrcNames[static_cast<int>(bgSrc)], step);
fH.FillH1(name, parRec.fMinv);
fH.FillH2("hSrcBgPairs", static_cast<int>(step) + 0.5, static_cast<double>(bgSrc) + 0.5);
- }
+
+ if (step == ELmvmAnaStep::ElId) {
+ string hName = "hMinvBgSource2_elid_";
+
+ if (std::abs(candP.fMcPdg) == 11) {
+
+ // cand1 is El and from Gamma
+ if (candP.IsMcGamma()) {
+ if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gg", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0()) fH.FillH1(hName + "gpi0", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 211) fH.FillH1(hName + "gpi", parRec.fMinv);
+ else fH.FillH1(hName + "go", parRec.fMinv);
+ }
+
+ // cand1 is El and from Pi0
+ else if (candP.IsMcPi0()) {
+ if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gpi0", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0()) fH.FillH1(hName + "pi0pi0", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 211) fH.FillH1(hName + "pipi0", parRec.fMinv);
+ else fH.FillH1(hName + "pi0o", parRec.fMinv);
+ }
+
+ // cand1 is El but not from Gamma or Pi0
+ else {
+ if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "go", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0()) fH.FillH1(hName + "pi0o", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 211) fH.FillH1(hName + "pio", parRec.fMinv);
+ else fH.FillH1(hName + "oo", parRec.fMinv);
+ }
+ }
+
+ // cand1 is misid. charged pion
+ else if (std::abs(candP.fMcPdg) == 211) {
+ if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gpi", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0()) fH.FillH1(hName + "pipi0", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 211) fH.FillH1(hName + "pipi", parRec.fMinv);
+ else fH.FillH1(hName + "pio", parRec.fMinv);
+ }
+
+ // cand1 is neither electron nor misid. charged pion
+ else {
+ if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "go", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0()) fH.FillH1(hName + "pi0o", parRec.fMinv);
+ else if (std::abs(candM.fMcPdg) == 211) fH.FillH1(hName + "pipi0", parRec.fMinv);
+ else fH.FillH1(hName + "oo", parRec.fMinv);
+ }
+ }
+ //}
}
}
@@ -696,8 +794,16 @@ void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg)
double stepBin = static_cast<double>(step) + 0.5;
FillMomHists(nullptr, &cand, cand.fMcSrc, step);
+ fH.FillH1("hCandPdg", step, cand.fMcPdg);
+
+ int absPdg = std::abs(pdg);
+ double pdgBin = (absPdg == 11 && cand.IsMcSignal()) ? 0.5 : (absPdg == 11 && !cand.IsMcSignal()) ? 1.5 : (absPdg == 211) ? 2.5 : (absPdg == 2212) ? 3.5 : (absPdg == 321) ? 4.5 : 5.5;
+ fH.FillH2("hCandPdgVsMom", step, cand.fMomentum.Mag(), pdgBin);
- if (cand.IsMcSignal()) { fH.FillH1("hNofSignalTracks", stepBin, fW); }
+ if (cand.IsMcSignal()) {
+ fH.FillH1("hNofSignalTracks", stepBin, fW);
+ fH.FillH2("hCandElSrc", stepBin, 7.5, fW);
+ }
else {
if (LmvmUtils::IsMismatch(cand)) fH.FillH1("hNofMismatches_all", stepBin);
if (cand.fStsMcTrackId != cand.fRichMcTrackId) fH.FillH1("hNofMismatches_rich", stepBin);
@@ -720,22 +826,28 @@ void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg)
double srcBin = 0.0;
if (cand.IsMcGamma()) srcBin = 0.5;
- else if (cand.IsMcPi0())
- srcBin = 1.5;
- else if (std::abs(pdg) == 211)
- srcBin = 2.5;
- else if (pdg == 2212)
- srcBin = 3.5;
- else if (std::abs(pdg) == 321)
- srcBin = 4.5;
- else if ((std::abs(pdg) == 11) && !cand.IsMcGamma() && !cand.IsMcPi0() && !cand.IsMcSignal())
- srcBin = 5.5;
- else
- srcBin = 6.5;
+ else if (cand.IsMcPi0()) srcBin = 1.5;
+ else if (std::abs(pdg) == 211) srcBin = 2.5;
+ else if (pdg == 2212) srcBin = 3.5;
+ else if (std::abs(pdg) == 321) srcBin = 4.5;
+ else if ((std::abs(pdg) == 11) && !cand.IsMcGamma() && !cand.IsMcPi0() && !cand.IsMcSignal())srcBin = 5.5;
+ else srcBin = 6.5;
fH.FillH2("hBgSrcTracks", stepBin, srcBin);
+ if (std::abs(cand.fMcPdg) == 11) fH.FillH2("hCandElSrc", stepBin, srcBin);
}
}
+void LmvmTask::BgPairPdg(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step)
+{
+ int pdgX = candP.fMcPdg;
+ int pdgY = candM.fMcPdg;
+
+ double pdgBinX = (std::abs(pdgX) == 11 && candP.IsMcSignal()) ? 0.5 : (std::abs(pdgX) == 11 && !candP.IsMcSignal()) ? 1.5 : (std::abs(pdgX) == 211) ? 2.5 : (pdgX == 2212) ? 3.5 : (pdgX == 321) ? 4.5 : (pdgX == 3112 or pdgX == 3222) ? 5.5 : (std::abs(pdgX) == 13) ? 6.5 : 7.5;
+ double pdgBinY = (std::abs(pdgY) == 11 && candM.IsMcSignal()) ? 0.5 : (std::abs(pdgY) == 11 && !candM.IsMcSignal()) ? 1.5 : (std::abs(pdgY) == 211) ? 2.5 : (pdgY == 2212) ? 3.5 : (pdgY == 321) ? 4.5 : (pdgY == 3112 or pdgY == 3222) ? 5.5 : (std::abs(pdgY) == 13) ? 6.5 : 7.5;
+
+ fH.FillH2("hBgPairPdg", step, pdgBinX, pdgBinY);
+}
+
void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const LmvmKinePar& parMc,
const LmvmKinePar& parRec, ELmvmAnaStep step)
{
@@ -754,6 +866,7 @@ void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const
fH.FillH1("hMomPairSignal", step, parMc.fMomentumMag, fW);
}
if (src == ELmvmSrc::Bg) {
+ BgPairPdg(candP, candM, step);
if (isMismatch) { fH.FillH1("hMinvBgMatch_mismatch", step, parRec.fMinv); }
else {
fH.FillH1("hMinvBgMatch_trueMatch", step, parRec.fMinv);
@@ -772,6 +885,7 @@ void LmvmTask::SignalAndBgReco()
CheckTopologyCut(ELmvmTopologyCut::TT, "hTtCut");
CheckTopologyCut(ELmvmTopologyCut::RT, "hRtCut");
if (fUseMvd) {
+ cout << "I am in LmvmTask::SignalAndBgReco() -> fUseMvd" << endl; // TODO: delete this line
CheckClosestMvdHit(1, "hMvdCut_1", "hMvdCutQa_1");
CheckClosestMvdHit(2, "hMvdCut_2", "hMvdCutQa_2");
}
@@ -783,7 +897,10 @@ void LmvmTask::SignalAndBgReco()
if (mcTrack != nullptr) pdg = mcTrack->GetPdgCode();
}
for (auto step : fH.fAnaSteps) {
- if (cand.IsCutTill(step)) TrackSource(cand, step, pdg);
+ if (cand.IsCutTill(step)) {
+ TrackSource(cand, step, pdg);
+ fH.FillH2("hPtYCandidate", step, cand.fRapidity, cand.fMomentum.Perp(), fW);
+ }
}
}
@@ -960,10 +1077,22 @@ void LmvmTask::DifferenceSignalAndBg()
if (!cand.fIsChi2Prim) continue;
fH.FillH1("hAnnRich", cand.fMcSrc, cand.fRichAnn, fW);
fH.FillH2("hAnnRichVsMom", cand.fMcSrc, cand.fMomentum.Mag(), cand.fRichAnn, fW);
- fH.FillH1("hAnnTrd", cand.fMcSrc, cand.fTrdAnn, fW);
- fH.FillH2("hTofM2", cand.fMcSrc, cand.fMomentum.Mag(), cand.fMass2, fW);
+ //fH.FillH1("hAnnTrd", cand.fMcSrc, cand.fTrdAnn, fW); // TODO: uncomment when TRD ANN is working (CbmLitGlobalElectronId::GetTrdAnn() gives back El-Likelihood)
fH.FillH2("hTrdLike_El", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikeEl, fW);
fH.FillH2("hTrdLike_Pi", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikePi, fW);
+ fH.FillH2("hTofM2", cand.fMcSrc, cand.fMomentum.Mag(), cand.fMass2, fW);
+
+ // electron purity
+ if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) == 11) {
+ fH.FillH2("hAnnRichVsMomPur_El", cand.fMomentum.Mag(), cand.fRichAnn, fW);
+ fH.FillH2("hTrdElLikePur_El", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW);
+ fH.FillH2("hTofM2Pur_El", cand.fMomentum.Mag(), cand.fMass2, fW);
+ }
+ else if (!cand.IsMcSignal() && std::abs(cand.fMcPdg != 11)) {
+ fH.FillH2("hAnnRichVsMomPur_Bg", cand.fMomentum.Mag(), cand.fRichAnn, fW);
+ fH.FillH2("hTrdElLikePur_Bg", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW);
+ fH.FillH2("hTofM2Pur_Bg", cand.fMomentum.Mag(), cand.fMass2, fW);
+ }
if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue;
//fH.FillSourceH1("hPt", cand.fMcSrc, cand.fMomentum.Perp(), fW);
@@ -1016,6 +1145,7 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s
if (fCands[iC].IsCutTill(ELmvmAnaStep::ElId)) {
CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(fCands[iC].fStsInd));
if (stsTrack == nullptr) continue;
+ cout << "NofMvdHits = " << stsTrack->GetNofMvdHits() << endl; // TODO: delete
for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
CbmMvdHit* candHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
if (candHit != nullptr && candHit->GetStationNr() == mvdStationNum) {
@@ -1025,6 +1155,9 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s
}
}
+ cout << "mvdV.size() = " << mvdV.size() << endl; // TODO: delete this line
+ cout << "candV.size() = " << candV.size() << endl; // TODO: delete this line
+
for (size_t iC = 0; iC < candV.size(); iC++) {
LmvmCand& cand = fCands[candV[iC].fInd];
double minD = 9999999.;
@@ -1068,6 +1201,9 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s
if (mvdStationNum == 1) cand.fIsMvd1Cut = isMvdCut;
else if (mvdStationNum == 2)
cand.fIsMvd2Cut = isMvdCut;
+
+ cout << "cand.fIsMvd1Cut = " <<cand.fIsMvd1Cut << endl; // TODO: delete these cout lines
+ cout << "cand.fIsMvd2Cut = " <<cand.fIsMvd2Cut << endl;
}
}
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h
index 688ccc3a11..bb60cd1c90 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h
@@ -87,6 +87,8 @@ public:
void TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg);
+ void BgPairPdg(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step);
+
void DoMcTrack();
void DoMcPair();
@@ -137,7 +139,7 @@ public:
virtual void Finish();
- void FillPionsHist();
+ void FillAccRecVsMomHist();
ClassDef(LmvmTask, 1);
diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx
index 4980b92a3a..347c802c19 100644
--- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx
@@ -241,8 +241,13 @@ void LmvmUtils::IsElectron(int globalTrackIndex, double momentum, double momentu
{
bool richEl = CbmLitGlobalElectronId::GetInstance().IsRichElectron(globalTrackIndex, momentum);
cand->fRichAnn = CbmLitGlobalElectronId::GetInstance().GetRichAnn(globalTrackIndex, momentum);
- bool trdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum);
- cand->fTrdAnn = CbmLitGlobalElectronId::GetInstance().GetTrdAnn(globalTrackIndex, momentum);
+
+ bool trdEl = (momentum < 1.) ? true : CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum);
+ //bool trdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum);
+
+ //cand->fTrdAnn = CbmLitGlobalElectronId::GetInstance().GetTrdAnn(globalTrackIndex, momentum); //TODO: uncomment this and delete next line when TrdAnn is working again
+ cand->fTrdLikeEl = CbmLitGlobalElectronId::GetInstance().GetTrdAnn(globalTrackIndex, momentum);
+
bool tofEl = CbmLitGlobalElectronId::GetInstance().IsTofElectron(globalTrackIndex, momentum);
cand->fMass2 = CbmLitGlobalElectronId::GetInstance().GetTofM2(globalTrackIndex, momentum);
bool isMomCut = (momentumCut > 0.) ? (momentum < momentumCut) : true;
diff --git a/macro/analysis/dielectron/batch_job.py b/macro/analysis/dielectron/batch_job.py
index 42d42f2062..7cf0390715 100755
--- a/macro/analysis/dielectron/batch_job.py
+++ b/macro/analysis/dielectron/batch_job.py
@@ -28,6 +28,7 @@ def main():
os.chdir(workDir)
plutoFile = getPlutoPath(colSystem, colEnergy, plutoParticle, taskId)
+ #plutoFile = getPlutoPathTetyana(colSystem, colEnergy, plutoParticle, taskId)
urqmdFile = "/lustre/nyx/cbm/prod/gen/urqmd/auau/" + colEnergy + "/centr/urqmd.auau."+ colEnergy + ".centr." + str(taskId).zfill(5) + ".root"
@@ -69,5 +70,22 @@ def getPlutoPath(colSystem, colEnergy, plutoParticle, taskId):
elif plutoParticle == "urqmd":
return ""
+def getPlutoPathTetyana(colSystem, colEnergy, plutoParticle, taskId):
+ if plutoParticle == "rho0":
+ return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/rho0_pluto/sum_rho0_1000_001.zip#rho0_" + str(taskId).zfill(5) + ".root"
+ elif plutoParticle == "omegaepem":
+ return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/w_pluto/sum_w_1000_001.zip#w_" + str(taskId).zfill(5) + ".root"
+ elif plutoParticle == "omegadalitz":
+ return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/wdalitz_pluto/sum_wdalitz_1000_001.zip#wdalitz_" + str(taskId).zfill(5) + ".root"
+ elif plutoParticle == "phi":
+ return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/phi_pluto/sum_phi_1000_001.zip#phi_" + str(taskId).zfill(5) + ".root"
+ elif plutoParticle == "inmed":
+ return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/inmed_had_epem_" + colEnergy + "_pluto/sum_inmed_had_6800_001.zip#inmed_had_epem_8gev_" + str(taskId).zfill(5) + ".root"
+ elif plutoParticle == "qgp":
+ return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/qgp_epem_8gev_pluto/sum_qgp_6800_001.zip#qgp_epem_8gev_" + str(taskId).zfill(5) + ".root"
+ elif plutoParticle == "urqmd":
+ return ""
+
+
if __name__ == '__main__':
main()
diff --git a/macro/analysis/dielectron/batch_send.py b/macro/analysis/dielectron/batch_send.py
index 76305c1719..df4835a343 100755
--- a/macro/analysis/dielectron/batch_send.py
+++ b/macro/analysis/dielectron/batch_send.py
@@ -4,10 +4,10 @@ import os
import shutil
def main():
- nofJobs = 1000
+ nofJobs = 1000
timeLimit = "08:00:00"
geoSetup = "sis100_electron"
- plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"]
+ plutoParticles = ["inmed", "omegadalitz", "omegaepem", "phi", "qgp"]
dataDir = "/lustre/nyx/cbm/users/criesen/data/lmvm/"
@@ -21,7 +21,7 @@ def main():
os.makedirs(dataDir,exist_ok = True)
for plutoParticle in plutoParticles:
- jobName = "LMVM_" + plutoParticle
+ jobName = plutoParticle
dataDirPluto = dataDir + plutoParticle
logFile = dataDirPluto + "/log/log_slurm-%A_%a.out"
errorFile = dataDirPluto + "/error/error_slurm-%A_%a.out"
diff --git a/macro/analysis/dielectron/draw_all.py b/macro/analysis/dielectron/draw_all.py
index ef19322f7c..6ffaab8e60 100755
--- a/macro/analysis/dielectron/draw_all.py
+++ b/macro/analysis/dielectron/draw_all.py
@@ -4,34 +4,41 @@ import os, shutil
def main():
plutoParticles = ["inmed", "omegadalitz", "omegaepem", "phi", "qgp"]
- dataDir = "/lustre/nyx/cbm/users/criesen/data/lmvm/"
- dataDirOut = dataDir + "/results/"
+
+ # paths on lustre
macroDir = "/lustre/cbm/users/criesen/cbmroot/macro/analysis/dielectron/"
- useMvd = True
+ dataDir = "/lustre/cbm/users/criesen/data/lmvm/"
+
+ # local paths
+ #macroDir = "/home/aghoehne/soft/cbm/cbmroot/macro/analysis/dielectron/"
+ #dataDir = "/home/aghoehne/soft/cbm/data/output/"
- if os.path.exists(dataDirOut):
- shutil.rmtree(dataDirOut)
- os.mkdir(dataDirOut)
+ resultDir = dataDir + "results/"
+
+ useMvd = False
+ if os.path.exists(resultDir):
+ shutil.rmtree(resultDir)
+ os.mkdir(resultDir)
for plutoParticle in plutoParticles:
- resultDir = dataDirOut + plutoParticle
+ resultDirPP = resultDir + plutoParticle
- resultDirAna = resultDir + "/lmvm/"
+ resultDirAna = resultDirPP + "/lmvm/"
inRootAna = dataDir + plutoParticle + "/analysis.all.root"
os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, inRootAna, resultDirAna, useMvd))
- resultDirLitqa = resultDir + "/litqa/"
+ resultDirLitqa = resultDirPP + "/litqa/"
inRootLitqa = dataDir + plutoParticle + "/litqa.all.root"
- os.system(('root -l -b -q {}/draw_litqa.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inRootLitqa, resultDirLitqa))
+ #os.system(('root -l -b -q {}/draw_litqa.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inRootLitqa, resultDirLitqa))
allInmed = dataDir + "/inmed/analysis.all.root"
- allQgp = dataDir + "/qgp/analysis.all.root"
- allOmega = dataDir + "/omegaepem/analysis.all.root"
- allPhi = dataDir + "/phi/analysis.all.root"
+ allQgp = dataDir + "/qgp/analysis.all.root"
+ allOmega = dataDir + "/omegaepem/analysis.all.root"
+ allPhi = dataDir + "/phi/analysis.all.root"
allOmegaD = dataDir + "/omegadalitz/analysis.all.root"
- resultDirAll = dataDirOut + "/all/"
+ resultDirAll = resultDir + "/all/"
os.system(('root -l -b -q {}/draw_analysis_all.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, allInmed, allQgp, allOmega, allPhi, allOmegaD, resultDirAll, useMvd))
diff --git a/macro/analysis/dielectron/draw_litqa.C b/macro/analysis/dielectron/draw_litqa.C
index 67d8e25e9a..bd092fe91d 100644
--- a/macro/analysis/dielectron/draw_litqa.C
+++ b/macro/analysis/dielectron/draw_litqa.C
@@ -12,11 +12,11 @@ void draw_litqa(const string& histRootFile = "/lustre/cbm/users/criesen/data/lmv
CbmSimulationReport* trackingQaReport = new CbmLitTrackingQaReport();
trackingQaReport->Create(histRootFile, outputDirTracking);
- string outputDirClustering = resultDir + "clustering/";
+ /*string outputDirClustering = resultDir + "clustering/";
gSystem->mkdir(outputDirClustering.c_str(), true);
CbmSimulationReport* clusteringQaReport = new CbmLitClusteringQaReport();
- clusteringQaReport->Create(histRootFile, outputDirClustering);
+ clusteringQaReport->Create(histRootFile, outputDirClustering);*/
// CbmSimulationReport* fitQaReport = new CbmLitFitQaReport();
// fitQaReport->Create(fileName, outputDir);
diff --git a/macro/analysis/dielectron/hadd_many.py b/macro/analysis/dielectron/hadd_many.py
index bc6f8aaca5..c16b6f56c7 100755
--- a/macro/analysis/dielectron/hadd_many.py
+++ b/macro/analysis/dielectron/hadd_many.py
@@ -10,9 +10,9 @@ def main():
inFilesLitQa = dataDirPluto + "/litqa.*.root"
outFileLitQa = dataDirPluto + "/litqa.all.root"
- if os.path.exists(outFileLitQa):
- os.remove(outFileLitQa)
- os.system(('hadd -j -T -f {} {}').format(outFileLitQa, inFilesLitQa))
+ #if os.path.exists(outFileLitQa):
+ # os.remove(outFileLitQa)
+ #os.system(('hadd -j -T -f {} {}').format(outFileLitQa, inFilesLitQa))
inFilesAna = dataDirPluto + "/analysis.*.root"
outFileAna = dataDirPluto + "/analysis.all.root"
diff --git a/macro/analysis/dielectron/run_local.py b/macro/analysis/dielectron/run_local.py
index 7eea67966c..32c3ca5482 100755
--- a/macro/analysis/dielectron/run_local.py
+++ b/macro/analysis/dielectron/run_local.py
@@ -18,13 +18,13 @@ def main():
colSystem = "auau"
plutoParticle = "phi"
- traFile = dataDir + "/tra." + runId + ".root"
- parFile = dataDir + "/param." + runId + ".root"
- digiFile = dataDir + "/digi." + runId + ".root"
- recoFile = dataDir + "/reco." + runId + ".root"
- litQaFile = dataDir + "/litqa." + runId + ".root"
- analysisFile = dataDir + "/analysis." + runId + ".root"
- geoSimFile = dataDir + "/geosim." + runId + ".root"
+ traFile = dataDir + "tra." + runId + ".root"
+ parFile = dataDir + "param." + runId + ".root"
+ digiFile = dataDir + "digi." + runId + ".root"
+ recoFile = dataDir + "reco." + runId + ".root"
+ litQaFile = dataDir + "litqa." + runId + ".root"
+ analysisFile = dataDir + "analysis." + runId + ".root"
+ geoSimFile = dataDir + "geosim." + runId + ".root"
traCmd = ('root -l -b -q run_transport.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(urqmdFile, plutoFile, traFile, parFile, geoSimFile, geoSetup, nEvents)
digiCmd = ('root -l -b -q run_digi.C\(\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(traFile, parFile, digiFile, nEvents)
diff --git a/macro/analysis/dielectron/run_transport.C b/macro/analysis/dielectron/run_transport.C
index 447015d501..013ece75f1 100644
--- a/macro/analysis/dielectron/run_transport.C
+++ b/macro/analysis/dielectron/run_transport.C
@@ -27,7 +27,7 @@ void run_transport(const std::string& urqmdFile, // if "", no urqmd
run.SetParFileName(parFile.c_str());
run.SetGeoFileName(geoFile.c_str());
run.LoadSetup(geoSetup.c_str());
- run.SetTarget("Gold", 0.0025, 2.5); // for lmvm thickness = 0.0025; // 25 mum
+ run.SetTarget("Gold", 0.0025, 2.5, 0., 0., -44., 0.); // for lmvm thickness = 0.0025; // 25 mum
run.SetBeamPosition(0., 0., 0.1, 0.1);
//run.SetEngine(kGeant4);
//run.SetRandomEventPlane();
diff --git a/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx b/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx
index 829bc4fcd0..f949a87d25 100644
--- a/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx
+++ b/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx
@@ -30,8 +30,8 @@
#include <cmath>
CbmLitGlobalElectronId::CbmLitGlobalElectronId()
- : fTrdAnnCut(0.1)
- , fRichAnnCut(-0.4)
+ : fTrdAnnCut(0.8) // TODO: at the moment, TrdAnn is Likelilhood value; change this if TrdAnn works again // values loose | strict: 0.2 | 0.8
+ , fRichAnnCut(-0.0) // values loose | strict: -0.4 | 0.0
, fRichUseAnn(true)
, fRichMeanA(-1.)
, fRichMeanB(-1.)
@@ -100,6 +100,10 @@ Bool_t CbmLitGlobalElectronId::IsTrdElectron(Int_t globalTrackIndex, Double_t mo
if (NULL == trdTrack) return false;
Double_t ann = trdTrack->GetPidLikeEL();
+
+ // TODO @Cornelius: remove following P<1GeV condition
+ if (momentum < 1.) return true;
+
if (ann > fTrdAnnCut) return true;
else
return false;
@@ -110,14 +114,25 @@ Bool_t CbmLitGlobalElectronId::IsTofElectron(Int_t globalTrackIndex, Double_t mo
Double_t mass2 = GetTofM2(globalTrackIndex, momentum, eventTime);
if (mass2 == -1.) return false;
+ // loose cut
if (momentum >= 1.3) {
if (mass2 < (0.010 * momentum - 0.003)) { return true; }
}
else {
- if (mass2 < 0.01) {
- return true; //fTofM2
- }
+ if (mass2 < 0.01) { return true; }
+ }
+
+ // stricter cut
+ /*if (momentum >= 0.8) {
+ if (mass2 < (- 0.020 * momentum + 0.026)) { return true; }
}
+ else {
+ if (mass2 < 0.01) { return true; }
+ }*/
+
+ // simple linear cut
+ //if (mass2 < 0.01) return true;
+
return false;
}
diff --git a/reco/littrack/cbm/qa/tracking/CbmLitTrackingQa.cxx b/reco/littrack/cbm/qa/tracking/CbmLitTrackingQa.cxx
index 88f7c53653..90788c14ee 100644
--- a/reco/littrack/cbm/qa/tracking/CbmLitTrackingQa.cxx
+++ b/reco/littrack/cbm/qa/tracking/CbmLitTrackingQa.cxx
@@ -92,9 +92,8 @@ CbmLitTrackingQa::CbmLitTrackingQa()
, fTrdMatches(NULL)
, fTofPoints(NULL)
, fTofMatches(NULL)
- , fMcToRecoMap()
- , fTrdAnnCut(0.85)
- , fRichAnnCut(0.25)
+ , fRichAnnCut(0.0) // values loose | strict: -0.4 | 0.0
+ , fTrdAnnCut(0.8) // at the moment, TrdAnn is likelihood value // values loose | strict: 0.2 | 0.8
{
}
--
GitLab