diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h index 7e04f6e2b9d4fc94a806ed421f2c968d608e835d..38618bd2e318f1496f558acdc93fa880f40a88f4 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h @@ -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; diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h index f61c5018b1b5b1c122271bd15efb312ea9fa3e16..5302ca3699a1a8783f8e7de4a67b444734ddcfe9 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h @@ -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; diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx index b695a27511a9fd4ef25e6090edcf4d9a36a65fad..e0116bbba216bcac376efdc4e3367cf4769fa38d 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx @@ -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"); } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx index b60d087f20a597d0936a6f9e3922021b9138ac3f..b57627d709c13e3cf2abbe1a823658264202b465 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx @@ -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); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h index f6059c6b4b33f56b6daa7b4e2389578b9fec1daf..56415b5205882100be261486a21e222b4fbe40ae 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h @@ -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(); /** diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx index 42d5cc6b27919c831ba78427289f2f7ea808c5c1..078fa62735311acb1bf99932bec54e9fe90a45fc 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx @@ -1,6 +1,6 @@ -/* 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" diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx index 0e73e67a01432a94525567973160d88ab7f92dd1..bb05b6a1625644da1319f18feb69b5872bdc3c5f 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx @@ -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 } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx index ef6f2028f750bc835a73b2e7532802d8c3a365c4..4980b92a3a692a6a40405953e89a34d49e29537e 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx @@ -1,4 +1,4 @@ -/* 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]*/ diff --git a/macro/analysis/dielectron/batch_job.py b/macro/analysis/dielectron/batch_job.py index 07563982797cc0467eb6fae1b054f84631207db8..23282b2c5bb28afd6f9d3e93601302e51ed539a0 100755 --- a/macro/analysis/dielectron/batch_job.py +++ b/macro/analysis/dielectron/batch_job.py @@ -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)) diff --git a/macro/analysis/dielectron/batch_send.py b/macro/analysis/dielectron/batch_send.py index b1817cca2dd7a78762194cf8aea5066e34251da2..76305c1719540d2418b26486b6dd89fb05167886 100755 --- a/macro/analysis/dielectron/batch_send.py +++ b/macro/analysis/dielectron/batch_send.py @@ -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" diff --git a/macro/analysis/dielectron/draw_all.py b/macro/analysis/dielectron/draw_all.py index bfa69559f744dfb9f29373cc2d20970e76b836a6..1c57c45021138b055e1d88d7c86b619785b9dce2 100755 --- a/macro/analysis/dielectron/draw_all.py +++ b/macro/analysis/dielectron/draw_all.py @@ -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__': diff --git a/macro/analysis/dielectron/draw_analysis.C b/macro/analysis/dielectron/draw_analysis.C index 29a4e7efdba6f10020b393dcebc023a187f18eaf..f5369bc0beb3e87b020102582cc3fc07ebcc36fb 100644 --- a/macro/analysis/dielectron/draw_analysis.C +++ b/macro/analysis/dielectron/draw_analysis.C @@ -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); diff --git a/macro/analysis/dielectron/draw_analysis_all.C b/macro/analysis/dielectron/draw_analysis_all.C index e7bc5c2efc3e124af82157bfd2fa8dc110014450..f8f487bc56ca058547fb28e7be8db92a24c41cca 100644 --- a/macro/analysis/dielectron/draw_analysis_all.C +++ b/macro/analysis/dielectron/draw_analysis_all.C @@ -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", + const string& fileQgp = "/lustre/cbm/users/criesen/data/lmvm/qgp/analysis.all.root", + const string& fileOmega = "/lustre/cbm/users/criesen/data/lmvm/omegaepem/analysis.all.root", + const string& filePhi = "/lustre/cbm/users/criesen/data/lmvm/phi/analysis.all.root", + const string& fileOmegaD = "/lustre/cbm/users/criesen/data/lmvm/omegadalitz/analysis.all.root", + const string& resultDir = "/lustre/cbm/users/criesen/data/lmvm/results/", Bool_t useMvd = false) { LmvmDrawAll* draw = new LmvmDrawAll(); diff --git a/macro/analysis/dielectron/draw_litqa.C b/macro/analysis/dielectron/draw_litqa.C index d11a2438964ef113270c167982d4223735965fc5..67d8e25e9a9ccc78652500fe5903e3fa09536333 100644 --- a/macro/analysis/dielectron/draw_litqa.C +++ b/macro/analysis/dielectron/draw_litqa.C @@ -2,17 +2,17 @@ SPDX-License-Identifier: GPL-3.0-only Authors: Elena Lebedeva [committer], Semen Lebedev */ -void draw_litqa(const string& histRootFile = "/Users/slebedev/Development/cbm/data/lmvm/nov21/litqa.phi100k.all.root", - const string& resultDir = "/Users/slebedev/Development/cbm/data/lmvm/nov21/litqa/") +void draw_litqa(const string& histRootFile = "/lustre/cbm/users/criesen/data/lmvm/inmed/litqa.all.root", + const string& resultDir = "/lustre/cbm/users/criesen/data/lmvm/results/litqa/") { - string outputDirTracking = resultDir + "/tracking/"; + string outputDirTracking = resultDir + "tracking/"; gSystem->mkdir(outputDirTracking.c_str(), true); 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(); diff --git a/macro/analysis/dielectron/hadd_many.py b/macro/analysis/dielectron/hadd_many.py index a4c6a0a23eba7ad33a35f3067980b43133756ab4..bc6f8aaca5ecf9b462f6aff31e0aac44c13bd6d9 100755 --- a/macro/analysis/dielectron/hadd_many.py +++ b/macro/analysis/dielectron/hadd_many.py @@ -3,7 +3,7 @@ import shutil def main(): plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"] - dataDir = "/lustre/nyx/cbm/users/criesen/data/lmvm/" + dataDir = "/lustre/nyx/cbm/users/criesen/data/lmvm" for plutoParticle in plutoParticles: dataDirPluto = dataDir + "/" + plutoParticle diff --git a/macro/analysis/dielectron/run_analysis.C b/macro/analysis/dielectron/run_analysis.C index 37e6a65860eb0af2576729cbecb35db9759174ba..45a055edc98f47194d43e47519cfd07e41520a3a 100644 --- a/macro/analysis/dielectron/run_analysis.C +++ b/macro/analysis/dielectron/run_analysis.C @@ -26,9 +26,9 @@ void run_analysis(const std::string& traFile, const std::string& parFile, const gDebug = 0; FairRunAna* run = new FairRunAna(); - FairFileSource* inputSource = new FairFileSource(digiFile.c_str()); + FairFileSource* inputSource = new FairFileSource(recoFile.c_str()); inputSource->AddFriend(traFile.c_str()); - inputSource->AddFriend(recoFile.c_str()); + inputSource->AddFriend(digiFile.c_str()); run->SetSource(inputSource); run->SetOutputFile(analysisFile.c_str()); run->SetGenerateRunInfo(kFALSE); diff --git a/macro/analysis/dielectron/run_litqa.C b/macro/analysis/dielectron/run_litqa.C index c3f1a6a8cb21d5e7cdb192cf94c9647914529964..1fec1e796b196289ab627cdee56c921dc7be5daa 100644 --- a/macro/analysis/dielectron/run_litqa.C +++ b/macro/analysis/dielectron/run_litqa.C @@ -42,9 +42,9 @@ void run_litqa(const std::string& traFile, const std::string& parFile, const std timer.Start(); FairRunAna* run = new FairRunAna(); - FairFileSource* inputSource = new FairFileSource(digiFile.c_str()); + FairFileSource* inputSource = new FairFileSource(recoFile.c_str()); inputSource->AddFriend(traFile.c_str()); - inputSource->AddFriend(recoFile.c_str()); + inputSource->AddFriend(digiFile.c_str()); run->SetSource(inputSource); run->SetOutputFile(qaFile.c_str()); run->SetGenerateRunInfo(kFALSE);