diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h
index ff1383d03c74aa9165c6d1364b4ab5025a7b7921..36985558415d4ef0a6b423dce5f527a897d16522 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 4db42c028bdc85bc67c1a5f3718d7cb161d4a189..520bc4ba1092587f741354361748cb4218bc4d68 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 1ac829f5b59f00b8568c0e9aa43ac40de66f9867..1ed6d24b2e01e2d82674422e4e6e308e9e4a8b73 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 747228bcea49c2b2b40420116349ef6ef47941f1..ef171b6af7a10aad0eaa049a73c835d50084ec11 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 4a2871b6916b646f2b1fd2c6bb0d8fadcdebe432..9d138f1c309225feeb93fee48d2919fc9f50c02f 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 414830101a89c32b8c994a0dacc161f2558b358a..7b8f55935fe94ccc566bbdeeb22978b8a4863b1a 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 3e09d218f3639d2190cdcb3dc90ad04acf7660f0..f67f6a618375b5e42541f11ed27080de4a9ba18e 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 688ccc3a1120fdc6222f5bee739e7a78acd85436..bb60cd1c905b761bf55295a89d8efbeb067add4f 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 4980b92a3a692a6a40405953e89a34d49e29537e..347c802c19d24ae7909c1ce877707fa6aa9193de 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 42d42f206210322cb26d85e30a92fb04db70fe18..7cf03907153e626ecf133fcf06f4b42d0e69c8fe 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 76305c1719540d2418b26486b6dd89fb05167886..df4835a343f8c5868837e00343fa3dbb26f0b5c0 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 ef19322f7c8180fd63fb91a309a52fb924bac1d2..6ffaab8e60438c9656f5d981dbc9ff60de97db2b 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 67d8e25e9a9ccc78652500fe5903e3fa09536333..bd092fe91d5580243b3bb48cc861ff277da5d0fd 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 bc6f8aaca5ecf9b462f6aff31e0aac44c13bd6d9..c16b6f56c7971410b6fee1684be70f456635492d 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 7eea67966c3e9badae5af6ad9917de1480188ddc..32c3ca54828a066f3724434437a435b81122b496 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 447015d501fe1fbf9a43b46c1d843c5db5553011..013ece75f1d58ef6266e5b2a5c3b95376041dc5c 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 829bc4fcd047396b5e918dba878056d4824e4ebc..f949a87d25bcdf325353d5bb43d2f2f5f4dde96b 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 88f7c53653abe953116ac2aa6fc1631e88161f9c..90788c14ee3386658019a068c777c5ce050c02a6 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
{
}