Skip to content

GitLab

Commit 5dd72067 authored by Cornelius Feier-Riesen's avatar Cornelius Feier-Riesen
Browse files

Implement generic rebinning/scaling and adjust topo cut values

parent a8601a15
Hide whitespace changes
Inline Side-by-side
Showing with 129 additions and 118 deletions
analysis/PWGDIL/dielectron/lmvm/LmvmCand.h
View file @ 5dd72067
......@@ -73,7 +73,7 @@ public:
// Cuts. If true then cut is passed
bool fIsChi2Prim = false;
bool fIsElectron = false;
bool fIsGammaCut = false;
bool fIsGammaCut = true; // Will be set to 'false' as soon as a partner with minv < 25 MeV is found
bool fIsMvd1Cut = false;
bool fIsMvd2Cut = false;
bool fIsTtCut = false;
......
analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h
View file @ 5dd72067
......@@ -86,12 +86,12 @@ public:
double fAngleCut = 1.;
double fChi2PrimCut = 3.;
double fGammaCut = 0.025;
double fStCutAngle = 2.4;
double fStCutPP = 1.;
double fTtCutAngle = 2.0;
double fTtCutPP = 3.0;
double fRtCutAngle = 1.2;
double fRtCutPP = 1.6;
double fStCutAngle = 2.3;
double fStCutPP = 2.9;
double fTtCutAngle = 2.2;
double fTtCutPP = 3.2;
double fRtCutAngle = 2.4;
double fRtCutPP = 3.0;
double fMvd1CutP = 1.2;
double fMvd1CutD = 0.4;
double fMvd2CutP = 1.5;
......
analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx
View file @ 5dd72067
......@@ -58,10 +58,10 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir,
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); });
//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); });
DrawAnaStepMany("lmvm_minv_matching", [this](ELmvmAnaStep step) { DrawMinvMatching(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); });
......@@ -71,7 +71,6 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir,
DrawAnaStepMany("lmvm_" + hName + "EpEm", [this, hName](ELmvmAnaStep step) { DrawSrcAnaStepEpEmH1(hName, step); });
}
DrawMomAccEpEm();
DrawMisc();
DrawGammaVertex();
DrawCuts();
......@@ -804,5 +803,5 @@ void LmvmDraw::DrawMvdAndStsHist()
void LmvmDraw::SaveCanvasToImage()
{
fH.fHM.SaveCanvasToImage(fOutputDir, "png;eps"); // fHM->SaveCanvasToImage(fOutputDir, "png;eps");
fH.fHM.SaveCanvasToImage(fOutputDir, "png"); // fHM->SaveCanvasToImage(fOutputDir, "png;eps");
}
analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx
View file @ 5dd72067
......@@ -57,7 +57,6 @@ void LmvmDrawAll::DrawHistFromFile(const string& fileInmed, const string& fileQg
TFile* file = new TFile(fileNames[i].c_str());
fH[i]->fHM.ReadFromFile(file);
int nofEvents = (int) fH[i]->H1("hEventNumber")->GetEntries();
//fHM[i]->ScaleByPattern(".*", 1. / nofEvents); // TODO: keep this commented?
LOG(info) << "Signal:" << fHMean.fSignalNames[i] << " nofEvents:" << nofEvents << endl;
}
......@@ -90,7 +89,7 @@ int LmvmDrawAll::GetNofTotalEvents()
}
template<class T>
void LmvmDrawAll::CreateMeanHist(const string& name, int nofEvents)
void LmvmDrawAll::CreateMeanHist(const string& name, int nofEvents, int nofRebins)
{
for (ELmvmSignal sig : fHMean.fSignals) {
if (static_cast<int>(sig) == 0) fHMean.fHM.Add(name, static_cast<T*>(H(sig)->GetObject(name)->Clone()));
......@@ -98,6 +97,11 @@ void LmvmDrawAll::CreateMeanHist(const string& name, int nofEvents)
static_cast<T*>(fHMean.GetObject(name))->Add(static_cast<T*>(H(sig)->GetObject(name)->Clone()));
}
static_cast<T*>(fHMean.GetObject(name))->Scale(1. / (double) nofEvents);
if (nofRebins > 0) {
static_cast<T*>(fHMean.GetObject(name))->Rebin(nofRebins);
double binWidth = static_cast<T*>(fHMean.GetObject(name))->GetBinWidth(1);
static_cast<T*>(fHMean.GetObject(name))->Scale(1. / binWidth);
}
}
void LmvmDrawAll::CreateMeanHistAll()
......@@ -106,14 +110,13 @@ void LmvmDrawAll::CreateMeanHistAll()
for (auto step : fHMean.fAnaSteps) {
for (auto src : {ELmvmSrc::Bg, ELmvmSrc::Eta, ELmvmSrc::Pi0}) {
CreateMeanHist<TH1D>(fHMean.GetName("hMinv", src, step), nofEvents);
CreateMeanHist<TH1D>(fHMean.GetName("hMinv", src, step), nofEvents, 20);
CreateMeanHist<TH2D>(fHMean.GetName("hMinvPt", src, step), nofEvents);
}
for (const string& comb : {"PM", "PP", "MM"}) {
for (const string& ev : {"sameEv", "mixedEv"}) {
// TODO: @Cornelius implement Proper scaling, nofEvents or nofMixedEvents
CreateMeanHist<TH2D>(fHMean.GetName("hMinvComb" + comb + "_" + ev, step), nofEvents);
CreateMeanHist<TH1D>(fHMean.GetName("hMinvComb" + comb + "_" + ev, step), nofEvents, 20);
}
}
}
......@@ -126,21 +129,26 @@ T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step)
{
T* sEta = dynamic_cast<T*>(fHMean.GetObject(fHMean.GetName(name, ELmvmSrc::Eta, step)));
T* sPi0 = dynamic_cast<T*>(fHMean.GetObject(fHMean.GetName(name, ELmvmSrc::Pi0, step)));
double binWidth = sEta->GetBinWidth(1);
T* coctail = nullptr;
T* cocktail = nullptr;
for (ELmvmSignal signal : fHMean.fSignals) {
string nameFull = fHMean.GetName(name, ELmvmSrc::Signal, step);
T* sHist = dynamic_cast<T*>(H(signal)->GetObject(nameFull)->Clone());
int nofEvents = (int) H(signal)->H1("hEventNumber")->GetEntries();
sHist->Scale(1. / nofEvents);
if (coctail == nullptr) coctail = sHist;
if (name != "hMinvPt") { // TODO: find another solution for this; rebin here has to be exactly the same as before
sHist->Rebin(20);
sHist->Scale(1. / binWidth);
}
if (cocktail == nullptr) cocktail = sHist;
else
coctail->Add(sHist);
cocktail->Add(sHist);
}
coctail->Add(sEta);
coctail->Add(sPi0);
cocktail->Add(sEta);
cocktail->Add(sPi0);
return coctail;
return cocktail;
}
void LmvmDrawAll::DrawMinvAll()
......@@ -160,11 +168,14 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
TH1D* cocktail = GetCocktailMinvH1(step);
TH1D* sbg = static_cast<TH1D*>(bg->Clone());
sbg->Add(cocktail);
double binWidth = bg->GetBinWidth(1);
vector<TH1D*> sHists(fHMean.fNofSignals);
for (ELmvmSignal signal : fHMean.fSignals) {
TH1D* sHist = H(signal)->H1Clone("hMinv", ELmvmSrc::Signal, step);
sHist->Scale(1. / H(signal)->H1("hEventNumber")->GetEntries());
sHist->Rebin(20); // TODO urgent: find solution to rebin w. same value as before!
sHist->Scale(1. / binWidth);
sHists[static_cast<int>(signal)] = sHist;
}
......@@ -187,16 +198,12 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step)
drawData.emplace_back(cocktail, -1, kRed + 2, 4, -1, "Cocktail");
double binWidth = drawData[0].fH->GetBinWidth(1);
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);
for (size_t i = 0; i < drawData.size(); i++) {
const auto& d = drawData[i];
// TODO: Think about rebin
d.fH->Rebin(20);
d.fH->Scale(1. / binWidth);
d.fH->GetYaxis()->SetTitle("dN/dM_{ee} [GeV/c^{2}]^{-1}");
d.fH->GetYaxis()->SetLabelSize(0.05);
......@@ -244,8 +251,6 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
c->cd(i++);
TH1D* pp = fHMean.H1Clone("hMinvCombPP_sameEv", step);
TH1D* mm = fHMean.H1Clone("hMinvCombMM_sameEv", step);
pp->Rebin(20);
mm->Rebin(20);
mm->Divide(pp);
mm->GetYaxis()->SetTitle("Ratio e^{-}e^{-}/e^{+}e^{+}");
DrawH1(mm, kLinear, kLinear, "hist");
......@@ -266,9 +271,6 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
TH1D* pp = fHMean.H1Clone("hMinvCombPP_" + ev, step);
TH1D* mm = fHMean.H1Clone("hMinvCombMM_" + ev, step);
TH1D* pm = fHMean.H1Clone("hMinvCombPM_" + ev, step);
pp->Rebin(20);
mm->Rebin(20);
pm->Rebin(20);
pm->GetYaxis()->SetTitle("M_{ee}");
DrawH1({pm, pp, mm}, {"e+e-", "e+e+", "e-e-"}, kLinear, kLog, true, 0.85, 0.7, 0.99, 0.99, "HIST");
fHMean.DrawAnaStepOnPad(step);
......@@ -292,10 +294,31 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
int maxBin = same->FindBin(0.7);
double scale = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin);
mixed->Scale(scale);
same->Rebin(20);
mixed->Rebin(20);
DrawH1({same, mixed}, {"Same " + comb, "Mixed " + comb}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99, "HIST");
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);
}
}
}
// Draw ratio mixed/same for PP, MM, PM cases
{
for (const string& comb : {"PM", "PP", "MM"}) {
string cName = "lmvmAll_CB_mixedOverSame" + comb;
TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800);
c->Divide(3, 3);
int i = 1;
for (auto step : fHMean.fAnaSteps) {
if (step < ELmvmAnaStep::ElId) continue;
c->cd(i++);
TH1D* same = fHMean.H1Clone("hMinvComb" + comb + "_sameEv", step);
TH1D* mixedOverSame = fHMean.CreateHByClone<TH1D>("hMinvComb" + comb + "_mixedEv", "hMinvMixedOverSame", step);
mixedOverSame->Divide(mixedOverSame, same, 1., 1., "B");
DrawH1(mixedOverSame, kLinear, kLog, "hist");
fHMean.DrawAnaStepOnPad(step);
}
}
}
......@@ -307,8 +330,6 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
TH1D* pmSame = fHMean.H1Clone("hMinvCombPM_sameEv", step);
TH1D* sbg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step);
sbg->Add(GetCocktailMinvH1(step));
pmSame->Rebin(20);
sbg->Rebin(20);
TH1D* ratio = (TH1D*) pmSame->Clone();
string cName = "lmvmAll_CB_SameVsSbg_" + fHMean.fAnaStepNames[static_cast<int>(step)];
......@@ -339,9 +360,7 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
int maxBin = same->FindBin(0.7);
double scale = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin);
mixed->Scale(scale);
same->Rebin(20);
mixed->Rebin(20);
same->GetYaxis()->SetTitle("M_{ee}");
same->GetXaxis()->SetTitle("M_{ee}");
DrawH1({same, mixed}, {"same", "mixed"}, kLinear, kLog, true, 0.8, 0.8, 0.99, 0.99, "HIST");
fHMean.DrawAnaStepOnPad(step);
}
......@@ -356,9 +375,6 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
TH1D* cbgAsm = fHMean.H1Clone("hMinvCombBgAsm", step);
TH1D* sbg = GetCocktailMinvH1(step);
sbg->Add(fHMean.H1("hMinv", ELmvmSrc::Bg, step));
cbg->Rebin(20);
cbgAsm->Rebin(20);
sbg->Rebin(20);
TH1D* ratio1 = static_cast<TH1D*>(sbg->Clone());
TH1D* ratio2 = static_cast<TH1D*>(sbg->Clone());
......@@ -388,7 +404,6 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
for (auto step : fHMean.fAnaSteps) {
if (step < ELmvmAnaStep::Reco) continue;
hists.push_back(fHMean.H1Clone("hMinvCombK", step));
hists.back()->Rebin(20);
legend.push_back(fHMean.fAnaStepLatex[static_cast<int>(step)]);
}
hists[0]->GetXaxis()->SetTitle("M_{ee}");
......@@ -407,12 +422,8 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
TH1D* pmSame = fHMean.H1Clone("hMinvCombPM_sameEv", step);
TH1D* combPMSignalAsm = fHMean.H1Clone("hMinvCombPMSignalAsm", step);
TH1D* combBgAsm = fHMean.H1Clone("hMinvCombBgAsm", step);
TH1D* coctail = GetCocktailMinvH1(step);
pmSame->Rebin(20);
combBgAsm->Rebin(20);
combPMSignalAsm->Rebin(20);
coctail->Rebin(20);
DrawH1({pmSame, combBgAsm, combPMSignalAsm, coctail},
TH1D* cocktail = GetCocktailMinvH1(step);
DrawH1({pmSame, combBgAsm, combPMSignalAsm, cocktail},
{"N_{same}^{+-}", "B_{c} (asm)", "Signal (N_{same}^{+-} - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.8,
0.8, 0.99, 0.99, "hist");
fHMean.DrawAnaStepOnPad(step);
......@@ -431,12 +442,8 @@ void LmvmDrawAll::DrawMinvCombSignalAndBg()
sbg->Add(GetCocktailMinvH1(step));
TH1D* combBgAsm = fHMean.H1Clone("hMinvCombBgAsm", step);
TH1D* sBgSignalAsm = fHMean.H1Clone("hMinvSBgSignalAsm", step);
TH1D* coctail = GetCocktailMinvH1(step);
sbg->Rebin(20);
combBgAsm->Rebin(20);
sBgSignalAsm->Rebin(20);
coctail->Rebin(20);
DrawH1({sbg, combBgAsm, sBgSignalAsm, coctail},
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");
fHMean.DrawAnaStepOnPad(step);
......@@ -450,17 +457,14 @@ void LmvmDrawAll::DrawSBgVsMinv()
for (ELmvmAnaStep step : {ELmvmAnaStep::TtCut, ELmvmAnaStep::PtCut}) {
TH1D* bg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step);
//TH1D* combBg = fHMean.H1Clone("hMinvCombBg", step);
TH1D* coctail = GetCocktailMinvH1(step);
TH1D* cocktailOverBg = fHMean.CreateHByClone<TH1D>("hMinv_bg", "hMinvCoctailOverBg", step);
bg->Rebin(20);
coctail->Rebin(20);
cocktailOverBg->Rebin(20);
TH1D* cocktail = GetCocktailMinvH1(step);
TH1D* cocktailOverBg = fHMean.CreateHByClone<TH1D>("hMinv_bg", "hMinvCocktailOverBg", step);
cocktailOverBg->GetYaxis()->SetTitle("Cocktail/Background");
cocktailOverBg->Divide(coctail, bg, 1., 1., "B");
cocktailOverBg->Divide(cocktail, bg, 1., 1., "B");
hists.push_back(cocktailOverBg);
}
fHMean.fHM.CreateCanvas("lmvmAll_minvCoctailOverBg", "lmvmAll_minvCoctailOverBg", 1000, 1000);
fHMean.fHM.CreateCanvas("lmvmAll_minvCocktailOverBg", "lmvmAll_minvCocktailOverBg", 1000, 1000);
DrawH1(hists, {"Without Pt cut", "With Pt cut"}, kLinear, kLog, true, 0.6, 0.85, 0.99, 0.99, "hist");
}
......@@ -485,13 +489,13 @@ void LmvmDrawAll::CalcCombBGHistos()
TH1D* pp = fHMean.H1Clone("hMinvCombPP_" + ev, step);
TH1D* mm = fHMean.H1Clone("hMinvCombMM_" + ev, step);
TH1D* geom = fHMean.CreateHByClone<TH1D>("hMinvCombMM_" + ev, "hMinvCombGeom_" + ev, step);
for (int i = 1; i <= geom->GetNbinsX(); i++) {
double cpp = pp->GetBinContent(i);
double cmm = mm->GetBinContent(i);
double content = std::sqrt(cpp * cmm);
if (cpp == 0 && cmm != 0) content = cmm;
else if (cpp != 0 && cmm == 0)
content = cpp;
double content = std::sqrt(
cpp
* cmm); // TODO: what, if one of these values is 0, even after rebinning? Should then last result be applied with a warning?
geom->SetBinContent(i, content);
}
}
......@@ -617,7 +621,7 @@ void LmvmDrawAll::CalcCutEffRange(double minMinv, double maxMinv)
stringstream ss;
ss << "lmvmAll_cutEff_" << minMinv << "to" << maxMinv;
fHMean.fHM.CreateCanvas(ss.str(), ss.str(), 1000, 1000);
DrawH1({hBg, hS}, {"BG", "Coctail"}, kLinear, kLinear, true, 0.75, 0.85, 1.0, 1.0, "hist");
DrawH1({hBg, hS}, {"BG", "Cocktail"}, kLinear, kLinear, true, 0.75, 0.85, 1.0, 1.0, "hist");
hBg->SetLineWidth(4);
hS->SetLineWidth(4);
hBg->SetMinimum(1);
......
analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h
View file @ 5dd72067
......@@ -64,7 +64,7 @@ private:
void DrawMinvCombSignalAndBg();
template<class T>
void CreateMeanHist(const std::string& name, int nofEvents);
void CreateMeanHist(const std::string& name, int nofEvents, int nofRebins = -1); // if nRebin = -1, no rebin
void CreateMeanHistAll();
/**
......
analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx
View file @ 5dd72067
/* Copyright (C) 2012-2021 UGiessen, JINR-LIT
/* Copyright (C) 2021 Justus-Liebig-Universitaet Giessen, Giessen
SPDX-License-Identifier: GPL-3.0-only
Authors: Elena Lebedeva, Semen Lebedev [committer] */
Authors: Semen Lebedev [committer] */
#include "LmvmHist.h"
......
analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx
View file @ 5dd72067
......@@ -190,7 +190,8 @@ void LmvmTask::Exec(Option_t*)
// if (impactPar <= 7.7) isCentralCollision = true;
// }
LOG(info) << "LmvmTask event number " << fEventNumber;
cout << "========================================================" << endl;
LOG(info) << "LmvmTask event number " << fEventNumber;
LOG(info) << "fPionMisidLevel = " << fPionMisidLevel;
LOG(info) << fCuts.ToString();
LOG(info) << "fW = " << fW;
......@@ -215,7 +216,6 @@ void LmvmTask::Exec(Option_t*)
fCandsTotal.insert(fCandsTotal.end(), fCands.begin(), fCands.end());
LOG(info) << "fCandsTotal.size = " << fCandsTotal.size();
//}
} // Exec
......@@ -259,6 +259,8 @@ void LmvmTask::FillMomHists(const CbmMCTrack* mct, const LmvmCand* cand, ELmvmSr
void LmvmTask::DoMcTrack()
{
int nMcTracks = fMCTracks->GetEntriesFast();
LOG(info) << "nMcTracks = " << nMcTracks;
for (int i = 0; i < nMcTracks; i++) {
CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i));
if (mct == nullptr) continue;
......@@ -276,16 +278,6 @@ void LmvmTask::DoMcTrack()
if (mct->GetNPoints(ECbmModuleId::kTrd) >= 2) fH.FillH1("hMomAcc" + chargeStr + "_trd", src, mom, fW);
if (mct->GetNPoints(ECbmModuleId::kTof) >= 1) fH.FillH1("hMomAcc" + chargeStr + "_tof", src, mom, fW);
if (std::abs(mct->GetPdgCode()) == 11) {
int mcMotherPdg = 0;
if (mct->GetMotherId() != -1) {
CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMCTracks->At(mct->GetMotherId()));
if (mother != nullptr) mcMotherPdg = mother->GetPdgCode();
}
fH.FillH1("hMotherPdg_mc", mcMotherPdg);
if (isAcc) fH.FillH1("hMotherPdg_acc", mcMotherPdg);
}
if (LmvmUtils::IsMcGammaEl(mct, fMCTracks)) {
TVector3 v;
mct->GetStartVertex(v);
......@@ -297,6 +289,20 @@ void LmvmTask::DoMcTrack()
fH.FillH2("hVertexGammaRZ", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
}
}
// Fill PDG histos
int mcPdg = mct->GetPdgCode();
if (std::abs(mcPdg) == 11 || mcPdg == 99009911) {
int mcMotherPdg = 0;
if (mct->GetMotherId() != -1) {
CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMCTracks->At(mct->GetMotherId()));
if (mother != nullptr) mcMotherPdg = mother->GetPdgCode();
}
if (std::abs(mcPdg) == 11) {
fH.FillH1("hMotherPdg_mc", mcMotherPdg);
if (isAcc) fH.FillH1("hMotherPdg_acc", mcMotherPdg);
}
}
}
}
......@@ -307,7 +313,8 @@ void LmvmTask::DoMcPair()
CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc1));
ELmvmSrc src = LmvmUtils::GetMcSrc(mct1, fMCTracks);
// To speed up: select only signal, eta and pi0 electrons
if (!(src == ELmvmSrc::Signal || src == ELmvmSrc::Pi0 || src == ELmvmSrc::Eta)) continue;
if (!(src == ELmvmSrc::Signal || src == ELmvmSrc::Pi0 || src == ELmvmSrc::Eta)) continue; // TODO: un-/comment?
bool isAcc1 = IsMcTrackAccepted(iMc1);
for (int iMc2 = iMc1 + 1; iMc2 < nMcTracks; iMc2++) {
......@@ -572,9 +579,9 @@ void LmvmTask::CombinatorialPairs()
for (size_t iC2 = iC1 + 1; iC2 < nCand; iC2++) {
const auto& cand2 = fCandsTotal[iC2];
if (cand2.IsMcSignal()) continue;
//double weight = (cand1.IsMcSignal() || cand2.IsMcSignal()) ? fW : 1.;
LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2);
//double weight = (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;
......@@ -682,7 +689,7 @@ void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAna
void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg)
{
// mo need to fill histograms for MC and Acc steps
// no need to fill histograms for MC and Acc steps
if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) return;
double stepBin = static_cast<double>(step) + 0.5;
......@@ -731,7 +738,7 @@ void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg)
void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const LmvmKinePar& parMc,
const LmvmKinePar& parRec, ELmvmAnaStep step)
{
// mo need to fill histograms for MC and Acc steps
// no need to fill histograms for MC and Acc steps
if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) return;
bool isMismatch = (LmvmUtils::IsMismatch(candP) || LmvmUtils::IsMismatch(candM));
ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM);
......@@ -785,6 +792,7 @@ void LmvmTask::SignalAndBgReco()
(candP.fStsMcTrackId >= 0) ? static_cast<CbmMCTrack*>(fMCTracks->At(candP.fStsMcTrackId)) : nullptr;
for (const auto& candM : fCands) {
if (candM.fCharge > 0) continue;
CbmMCTrack* mctrackM =
(candM.fStsMcTrackId >= 0) ? static_cast<CbmMCTrack*>(fMCTracks->At(candM.fStsMcTrackId)) : nullptr;
......@@ -806,9 +814,9 @@ void LmvmTask::CheckGammaConvAndPi0()
if (candM.fCharge > 0) continue;
if (candP.IsCutTill(ELmvmAnaStep::ElId) && candM.IsCutTill(ELmvmAnaStep::ElId)) {
LmvmKinePar pRec = LmvmKinePar::Create(&candP, &candM);
if (fCuts.IsGammaCutOk(pRec.fMinv)) {
candM.fIsGammaCut = true;
candP.fIsGammaCut = true;
if (!fCuts.IsGammaCutOk(pRec.fMinv)) {
candM.fIsGammaCut = false;
candP.fIsGammaCut = false;
}
}
}
......@@ -1041,7 +1049,7 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s
stsMotherId = (mct2 != nullptr) ? mct2->GetMotherId() : -2;
}
bin = (mvdMotherId != -1 && mvdMotherId == stsMotherId) ? 0.5 : bin = 1.5; // correct or wrong assignment
bin = (mvdMotherId != -1 && mvdMotherId == stsMotherId) ? 0.5 : 1.5; // correct or wrong assignment
if (cand.IsMcSignal()) {
bin = (mvdMotherId == stsMotherId && mcMvdHitPdg == 11) ? 0.5 : 1.5; // correct or wrong assignment
}
......
analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx
View file @ 5dd72067
/* Copyright (C) 2021 Justus-Liebig-Universitaet Giessen, Giessen
/* Copyright (C) 2010-2021 Justus-Liebig-Universitaet Giessen, Giessen
SPDX-License-Identifier: GPL-3.0-only
Authors: Semen Lebedev [committer]*/
......
macro/analysis/dielectron/batch_job.py
View file @ 5dd72067
......@@ -10,7 +10,7 @@ def main():
plutoParticle = sys.argv[3]
cbmrootConfigPath = "/lustre/nyx/cbm/users/criesen/build/config.sh"
macroDir = "/lustre/nyx/cbm/users/criesen/cbmroot/macro/analysis/dielectron/"
nofEvents = 1000
nofEvents = 500
taskId = os.environ.get('SLURM_ARRAY_TASK_ID')
jobId = os.environ.get('SLURM_ARRAY_JOB_ID')
......@@ -45,9 +45,9 @@ def main():
qaCmd = ('root -l -b -q {}/run_litqa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, traFile, parFile, digiFile, recoFile, litQaFile, geoSetup, nofEvents)
anaCmd = ('root -l -b -q {}/run_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, traFile, parFile, digiFile, recoFile, analysisFile, plutoParticle, colSystem, colEnergy, geoSetup, nofEvents)
os.system((". /{} -a; {}").format(cbmrootConfigPath, traCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, digiCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, recoCmd))
#os.system((". /{} -a; {}").format(cbmrootConfigPath, traCmd))
#os.system((". /{} -a; {}").format(cbmrootConfigPath, digiCmd))
#os.system((". /{} -a; {}").format(cbmrootConfigPath, recoCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, qaCmd))
os.system((". /{} -a; {}").format(cbmrootConfigPath, anaCmd))
......
macro/analysis/dielectron/batch_send.py
View file @ 5dd72067
......@@ -4,7 +4,7 @@ import os
import shutil
def main():
nofJobs = 100
nofJobs = 1000
timeLimit = "08:00:00"
geoSetup = "sis100_electron"
plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"]
......@@ -21,7 +21,7 @@ def main():
os.makedirs(dataDir,exist_ok = True)
for plutoParticle in plutoParticles:
jobName = "LMVM" + plutoParticle
jobName = "LMVM_" + plutoParticle
dataDirPluto = dataDir + plutoParticle
logFile = dataDirPluto + "/log/log_slurm-%A_%a.out"
errorFile = dataDirPluto + "/error/error_slurm-%A_%a.out"
......
macro/analysis/dielectron/draw_all.py
View file @ 5dd72067
......@@ -17,6 +17,7 @@ def main():
for plutoParticle in plutoParticles:
resultDir = dataDirOut + plutoParticle
resultDirAna = resultDir + "/lmvm/"
inRootAna = dataDir + plutoParticle + "/analysis.all.root"
os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, inRootAna, resultDirAna, useMvd, drawSig))
......@@ -32,6 +33,7 @@ def main():
allPhi = dataDir + "/phi/analysis.all.root"
allOmegaD = dataDir + "/omegadalitz/analysis.all.root"
resultDirAll = dataDirOut + "/lmvm_all/"
os.system(('root -l -b -q {}/draw_analysis_all.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, allInmed, allQgp, allOmega, allPhi, allOmegaD, resultDirAll, useMvd))
if __name__ == '__main__':
......
macro/analysis/dielectron/draw_analysis.C
View file @ 5dd72067
......@@ -4,10 +4,9 @@
//#include <experimental/filesystem>
void draw_analysis(
const string& histRootFile = "/Users/slebedev/Development/cbm/data/lmvm/nov21/analysis.inmed.all.root",
const string& resultDir = "/Users/slebedev/Development/cbm/data/lmvm/nov21/results/", Bool_t useMvd = true,
Bool_t drawSignificance = true)
void draw_analysis(const string& histRootFile = "/lustre/cbm/users/criesen/data/lmvm/inmed/analysis.all.root",
const string& resultDir = "/lustre/cbm/users/criesen/data/lmvm/results/", Bool_t useMvd = true,
Bool_t drawSignificance = true)
{
gSystem->mkdir(resultDir.c_str(), true);
......
macro/analysis/dielectron/draw_analysis_all.C
View file @ 5dd72067
......@@ -2,13 +2,12 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Elena Lebedeva [committer], Semen Lebedev */
void draw_analysis_all(
const string& fileInmed = "/Users/slebedev/Development/cbm/data/lmvm/nov21/analysis.inmed.all.root",
const string& fileQgp = "/Users/slebedev/Development/cbm/data/lmvm/nov21/analysis.qgp.all.root",
const string& fileOmega = "/Users/slebedev/Development/cbm/data/lmvm/nov21/analysis.omega.all.root",
const string& filePhi = "/Users/slebedev/Development/cbm/data/lmvm/nov21/analysis.phi.all.root",
const string& fileOmegaD = "/Users/slebedev/Development/cbm/data/lmvm/nov21/analysis.omegaD.all.root",
const string& resultDir = "/Users/slebedev/Development/cbm/data/lmvm/nov21/results_all/", Bool_t useMvd = false)
void draw_analysis_all(const string& fileInmed = "/lustre/cbm/users/criesen/data/lmvm/inmed/analysis.all.root",