diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h index 36985558415d4ef0a6b423dce5f527a897d16522..bfccb19a53ecd3eaa1b440689eadfac091744444 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmCuts.h @@ -57,11 +57,14 @@ public: double cutD = (stationNum == 1) ? fMvd1CutD : fMvd2CutD; double cutP = (stationNum == 1) ? fMvd1CutP : fMvd2CutP; double val = -1. * (cutP / cutD) * dmvd + cutP; - 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 + 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 e87d289669415e570e1d0d6c0b5fbf5a394ec86f..12ff21e3e9d75e7e18753ab5a1e5039606beb7e3 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx @@ -84,7 +84,7 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir, DrawMvdAndStsHist(); DrawAccRecMom(); DrawPmtXY(); - DrawMinvBg(); // TODO: do not extra method + DrawMinvBg(); // TODO: do not extra method DrawCombinatorialPairs(); SaveCanvasToImage(); @@ -301,9 +301,10 @@ void LmvmDraw::DrawCombinatorialPairs() { // Draw pair yields // TODO: scale to bin width string cbName = "CombPairsPluto/"; - + { - TCanvas* c3 = fH.fHM.CreateCanvas((cbName + "pairYields_steps").c_str(), (cbName + "pairYields_steps").c_str(), 1800, 600); + TCanvas* c3 = + fH.fHM.CreateCanvas((cbName + "pairYields_steps").c_str(), (cbName + "pairYields_steps").c_str(), 1800, 600); c3->Divide(3, 1); c3->cd(1); DrawAnaStepH1("hMinvCombPM_pluto_mixedEv", true); @@ -317,24 +318,26 @@ void LmvmDraw::DrawCombinatorialPairs() } { - TCanvas* c = fH.fHM.CreateCanvas((cbName + "pairYields").c_str(), (cbName + "pairYields").c_str(), 1800, 1800); - c->Divide(3,3); + TCanvas* c = fH.fHM.CreateCanvas((cbName + "pairYields").c_str(), (cbName + "pairYields").c_str(), 1800, 1800); + c->Divide(3, 3); int i = 1; - for (auto step : fH.fAnaSteps) { + for (auto step : fH.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); TH1* pm = fH.H1Clone("hMinvCombPM_pluto_mixedEv", step); TH1* pp = fH.H1Clone("hMinvCombPP_pluto_mixedEv", step); TH1* mm = fH.H1Clone("hMinvCombMM_pluto_mixedEv", step); - DrawH1({pm, pp, mm}, {"e^{+}e^{-} pairs (mixed ev.)", "e^{+}e^{+} pairs (mixed ev.)", "e^{-}e^{-} pairs (mixed ev.)"}, kLinear, kLog, true, 0.57, 0.79, 0.99, 0.99, "hist"); + DrawH1({pm, pp, mm}, + {"e^{+}e^{-} pairs (mixed ev.)", "e^{+}e^{+} pairs (mixed ev.)", "e^{-}e^{-} pairs (mixed ev.)"}, kLinear, + kLog, true, 0.57, 0.79, 0.99, 0.99, "hist"); fH.DrawAnaStepOnPad(step); } } // Draw ratio of e+e+/e-e- pairs - { - TCanvas* c = fH.fHM.CreateCanvas((cbName + "ratio_PPMM").c_str(), (cbName + "ratio_PPMM").c_str(), 1800, 1800); - c->Divide(3,3); + { + TCanvas* c = fH.fHM.CreateCanvas((cbName + "ratio_PPMM").c_str(), (cbName + "ratio_PPMM").c_str(), 1800, 1800); + c->Divide(3, 3); int i = 1; for (auto step : fH.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; @@ -350,8 +353,8 @@ void LmvmDraw::DrawCombinatorialPairs() // Draw geometric mean { - TCanvas* c = fH.fHM.CreateCanvas((cbName + "geomMean").c_str(), (cbName + "geomMean").c_str(), 1800, 1800); - c->Divide(3,3); + TCanvas* c = fH.fHM.CreateCanvas((cbName + "geomMean").c_str(), (cbName + "geomMean").c_str(), 1800, 1800); + c->Divide(3, 3); int i = 1; for (auto step : fH.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; @@ -369,9 +372,9 @@ void LmvmDraw::DrawCombinatorialPairs() } // draw k Factor - { - TCanvas* c = fH.fHM.CreateCanvas((cbName + "kFactor").c_str(), (cbName + "kFactor").c_str(), 1800, 1800); - c->Divide(3,3); + { + TCanvas* c = fH.fHM.CreateCanvas((cbName + "kFactor").c_str(), (cbName + "kFactor").c_str(), 1800, 1800); + c->Divide(3, 3); int i = 1; for (auto step : fH.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; @@ -379,7 +382,8 @@ void LmvmDraw::DrawCombinatorialPairs() TH1* k = fH.H1Clone("hMinvCombPM_pluto_mixedEv", step); TH1* gm = fH.H1Clone("hMinvCombPP_pluto_mixedEv", step); for (int iB = 1; iB <= k->GetNbinsX(); iB++) { - double con = std::sqrt(fH.H1("hMinvCombPP_pluto_mixedEv", step)->GetBinContent(iB) * fH.H1("hMinvCombMM_pluto_mixedEv", step)->GetBinContent(iB)); + double con = std::sqrt(fH.H1("hMinvCombPP_pluto_mixedEv", step)->GetBinContent(iB) + * fH.H1("hMinvCombMM_pluto_mixedEv", step)->GetBinContent(iB)); gm->SetBinContent(iB, con); } k->Divide(gm); @@ -449,7 +453,7 @@ void LmvmDraw::DrawCuts() Draw2DCut("hTrdLike_Pi"); Draw2DCut("hAnnRichVsMom"); Draw2DCut("hTofM2"); - + Draw1DCut("hChi2PrimVertex", "right", fCuts.fChi2PrimCut); //Draw1DCut("hPt", "left", fCuts.fPtCut); //Draw1DCut("hMom", "left"); @@ -544,7 +548,7 @@ void LmvmDraw::DrawMinvBg() 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; @@ -555,7 +559,7 @@ void LmvmDraw::DrawMinvBg() cout << "Entries pipi0: " << pipi0->GetEntries() << endl; cout << "Entries pio: " << pio->GetEntries() << endl; cout << "Entries pi0o: " << pi0o->GetEntries() << endl; - + int reb = 50; gg->Rebin(reb); @@ -563,13 +567,13 @@ void LmvmDraw::DrawMinvBg() 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."}; @@ -583,41 +587,45 @@ 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); + 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"); + DrawH1(pdgZoom, kLinear, kLog, "hist text40"); } // PID vs momentum - vector<string> yLabel = {"e^{#pm}_{PLUTO}", "e^{#pm}_{UrQMD}", "#pi^{#pm}", "p", "K^{+}", "o." }; + 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); + 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()); + 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; + 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(); + 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.; + 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); @@ -627,13 +635,13 @@ void LmvmDraw::DrawElPurity() // 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}"); - + {"#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); + 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"}) { @@ -645,8 +653,8 @@ void LmvmDraw::DrawElPurity() 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? + + 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.); @@ -683,7 +691,8 @@ 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((hist + fH.fSrcLatex[static_cast<int>(src)]).c_str(), 1, fH.H2(hist, src)->GetYaxis()->GetNbins(), "")); + 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]); } @@ -823,7 +832,10 @@ void LmvmDraw::DrawMinvBgPairSrc(ELmvmAnaStep step) vector<TH1*> hists; vector<string> latex; for (int i = 0; i < fH.fNofBgPairSrc; i++) { - 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 + 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 + "%)"); } @@ -853,20 +865,21 @@ void LmvmDraw::DrawAccRecMom() // 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> 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); + 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); + 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.)"; @@ -874,7 +887,7 @@ void LmvmDraw::DrawAccRecMom() i++; } - double y1 = 0.17; //(pdg == 211) ? 0.20 : 0.74; + 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"); @@ -908,7 +921,7 @@ void LmvmDraw::DrawAccRecMom() } void LmvmDraw::DrawSource2D(const string& cName, const string& hName, const vector<string>& yLabels, double scale, - const string& zTitle) + const string& zTitle) { fH.fHM.CreateCanvas((cName + "_abs").c_str(), (cName + "_abs").c_str(), 900, 600); TH2D* habs = fH.H2Clone(hName); @@ -975,8 +988,8 @@ void LmvmDraw::DrawBgSourceTracks() SetAnalysisStepAxis(purity); DrawSource2D("bg/SrcTracksBg_2d", "hBgSrcTracks", - {"#gamma", "#pi^{0}", "#pi^{#pm}", "p", "K", "e^{#pm}_{sec}", "oth.", "signal"}, 1000., - "Tracks per event x10^{3}"); + {"#gamma", "#pi^{0}", "#pi^{#pm}", "p", "K", "e^{#pm}_{sec}", "oth.", "signal"}, 1000., + "Tracks per event x10^{3}"); TCanvas* c = fH.fHM.CreateCanvas("nofTopoPairs", "nofTopoPairs", 1600, 800); c->Divide(2, 1); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h index 1ed6d24b2e01e2d82674422e4e6e308e9e4a8b73..362907dde02d17ec1929913d487c0626fadb83cb 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h @@ -117,7 +117,7 @@ private: void DrawMinvBg(); void DrawSource2D(const std::string& cName, const std::string& hName, const std::vector<std::string>& yLabels, - double scale, const std::string& zTitle); + double scale, const std::string& zTitle); void DrawBgSourceTracks(); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx index b6025606a479d333dac619b51e891d2d28a84002..dea2df01f3c329a94d6fffdd02e61c7093574ec9 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx @@ -168,7 +168,7 @@ T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step) void LmvmDrawAll::DrawMinvAll() { for (const ELmvmAnaStep step : fHMean.fAnaSteps) { - if (step < ELmvmAnaStep::ElId) continue; // TODO: activate also earlier histos in Task.cxx + if (step < ELmvmAnaStep::ElId) continue; // TODO: activate also earlier histos in Task.cxx string name = fHMean.GetName("lmvmAll_minv_", step); fHMean.fHM.CreateCanvas(name.c_str(), name.c_str(), 1000, 1000); DrawMinv(step); @@ -184,9 +184,9 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) TH1D* sbg = static_cast<TH1D*>(bg->Clone()); sbg->Add(cocktail); TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); - TH1D* cbSig = fHMean.H1Clone("hMinvCombSignalMc", step); + TH1D* cbSig = fHMean.H1Clone("hMinvCombSignalMc", step); TH1D* cbU = fHMean.H1Clone("hMinvCombBg_urqmd", step); - TH1D* cbSigU = fHMean.H1Clone("hMinvCombSignalMc_urqmd", step); + TH1D* cbSigU = fHMean.H1Clone("hMinvCombSignalMc_urqmd", step); double binWidth = bg->GetBinWidth(1); vector<TH1D*> sHists(fHMean.fNofSignals); @@ -252,7 +252,7 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) } void LmvmDrawAll::DrawMinvBgSourcesAll() -{ +{ TH1D* gg = fHMean.H1Clone("hMinvBgSource2_elid_gg"); TH1D* pipi = fHMean.H1Clone("hMinvBgSource2_elid_pipi"); TH1D* pi0pi0 = fHMean.H1Clone("hMinvBgSource2_elid_pi0pi0"); @@ -267,7 +267,7 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() TH1D* physBg = fHMean.H1Clone("hMinvBgSource2_elid_gg"); physBg->Add(fHMean.H1("hMinvBgSource2_elid_gpi0")); physBg->Add(fHMean.H1("hMinvBgSource2_elid_pi0pi0")); - + TH1D* cPi = fHMean.H1Clone("hMinvBgSource2_elid_pipi"); cPi->Add(fHMean.H1Clone("hMinvBgSource2_elid_gpi")); cPi->Add(fHMean.H1Clone("hMinvBgSource2_elid_pipi0")); @@ -276,7 +276,7 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() 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); @@ -286,7 +286,7 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() 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(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."); @@ -297,8 +297,8 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() 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); + 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()->SetLabelSize(0.05); @@ -315,17 +315,19 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() 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"); + 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); + 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); @@ -413,7 +415,9 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() mixed->Scale(scale); mixedU->Scale(scaleU); same->SetMinimum(1e-8); - DrawH1({same, mixed, sameU, mixedU}, {"Same " + comb, "Mixed " + comb, "Same " + comb + " (UrQMD-El)", "Mixed " + comb + " (UrQMD-El)"}, kLinear, kLog, true, 0.65, 0.7, 0.99, 0.99, "p"); + DrawH1({same, mixed, sameU, mixedU}, + {"Same " + comb, "Mixed " + comb, "Same " + comb + " (UrQMD-El)", "Mixed " + comb + " (UrQMD-El)"}, + kLinear, kLog, true, 0.65, 0.7, 0.99, 0.99, "p"); fHMean.DrawAnaStepOnPad(step); } } @@ -429,7 +433,8 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); - TH1D* same = fHMean.H1Clone("hMinvComb" + comb + "_sameEv", step); // TODO: check if this does not change origin histo + TH1D* same = + fHMean.H1Clone("hMinvComb" + comb + "_sameEv", step); // TODO: check if this does not change origin histo TH1D* mixedOverSame = fHMean.H1Clone("hMinvComb" + comb + "_mixedEv", step); mixedOverSame->Divide(mixedOverSame, same, 1., 1., "B"); DrawH1(mixedOverSame, kLinear, kLog, "hist"); @@ -480,7 +485,9 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() same->GetXaxis()->SetTitle("M_{ee}"); same->SetMinimum(1e-8); mixed->SetMinimum(1e-8); - DrawH1({same, mixed, sameU, mixedU}, {"geom. mean (same)", "geom. mean (mixed)", "geom. mean (same, UrQMD-El)", "geom. mean (mixed, UrQMD-El)"}, kLinear, kLog, true, 0.52, 0.8, 0.99, 0.99, "p"); + DrawH1({same, mixed, sameU, mixedU}, + {"geom. mean (same)", "geom. mean (mixed)", "geom. mean (same, UrQMD-El)", "geom. mean (mixed, UrQMD-El)"}, + kLinear, kLog, true, 0.52, 0.8, 0.99, 0.99, "p"); fHMean.DrawAnaStepOnPad(step); } } @@ -493,8 +500,8 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); - TH1D* k = fHMean.H1Clone("hMinvCombK", step); - TH1D* kU = fHMean.H1Clone("hMinvCombK_urqmd", step); + TH1D* k = fHMean.H1Clone("hMinvCombK", step); + TH1D* kU = fHMean.H1Clone("hMinvCombK_urqmd", step); double min = 0.6; double max = 1.4; k->SetMinimum(min); @@ -505,7 +512,7 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() fHMean.DrawAnaStepOnPad(step); } } - + // Draw CB vs MC-BG { TCanvas* c = fHMean.fHM.CreateCanvas("cb/cbVsMcBg", "cb/cbVsMcBg", 1800, 1800); @@ -529,8 +536,8 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); - TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); - TH1D* bg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step); + TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); + TH1D* bg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step); TH1D* ratio = static_cast<TH1D*>(bg->Clone()); ratio->Divide(cb); DrawH1({ratio}, {"BG_{MC}/B_{C}"}, kLinear, kLinear, true, 0.65, 0.75, 0.92, 0.9, "hist"); @@ -568,8 +575,8 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() pmSame->SetMaximum(1e-2); pmSame->SetMinimum(4e-9); DrawH1({pmSame, combBg, combSignalNpm, cocktail}, - {"N_{same}^{+-}", "B_{c}", "Signal (N_{same}^{+-} - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.53, - 0.75, 0.99, 0.92, "hist"); + {"N_{same}^{+-}", "B_{c}", "Signal (N_{same}^{+-} - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.53, 0.75, + 0.99, 0.92, "hist"); fHMean.DrawAnaStepOnPad(step); } } @@ -589,9 +596,8 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() sbg2->SetMaximum(2e-2); sbg2->SetMinimum(4e-9); TH1D* cocktail = GetCocktailMinvH1(step); - DrawH1({sbg2, combBg, sBgSignal, cocktail}, - {"Cocktail + BG", "B_{C}", "Signal (Cock+BG - B_{C})", "Cocktail"}, kLinear, kLog, true, 0.53, - 0.72, 0.99, 0.92, "p"); + DrawH1({sbg2, combBg, sBgSignal, 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); } } @@ -618,18 +624,19 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; string name = fHMean.GetName("cb/cbAndSignal_urqmd", step); - TCanvas* c = fHMean.fHM.CreateCanvas(name.c_str(), name.c_str(), 1800, 900); - TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); - TH1D* cbU = fHMean.H1Clone("hMinvCombBg_urqmd", step); - TH1D* sig = fHMean.H1Clone("hMinvCombSignalMc", step); - TH1D* sigU = fHMean.H1Clone("hMinvCombSignalMc_urqmd", step); + TCanvas* c = fHMean.fHM.CreateCanvas(name.c_str(), name.c_str(), 1800, 900); + TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); + TH1D* cbU = fHMean.H1Clone("hMinvCombBg_urqmd", step); + TH1D* sig = fHMean.H1Clone("hMinvCombSignalMc", step); + TH1D* sigU = fHMean.H1Clone("hMinvCombSignalMc_urqmd", step); cb->SetMinimum(1e-8); c->Divide(2, 1); c->cd(1); DrawH1({cb, cbU}, {"B_{C}", "B_{C} (UrQMD electrons)"}, kLinear, kLog, true, 0.6, 0.85, 0.95, 0.95, "p"); fHMean.DrawAnaStepOnPad(step); c->cd(2); - DrawH1({sig, sigU}, {"Signal (BG+Coc - B_{C})", "Signal (UrQMD electrons)"}, kLinear, kLog, true, 0.55, 0.85, 0.95, 0.95, "p"); + DrawH1({sig, sigU}, {"Signal (BG+Coc - B_{C})", "Signal (UrQMD electrons)"}, kLinear, kLog, true, 0.55, 0.85, + 0.95, 0.95, "p"); fHMean.DrawAnaStepOnPad(step); } } @@ -664,9 +671,9 @@ void LmvmDrawAll::SaveHist() void LmvmDrawAll::CalcCombBGHistos() { for (ELmvmAnaStep step : fHMean.fAnaSteps) { - for (const string& cat : {"", "_urqmd"}) { // common CB or from pure UrQMD electrons + for (const string& cat : {"", "_urqmd"}) { // common CB or from pure UrQMD electrons if (step < ELmvmAnaStep::ElId) continue; - + // Geometrical mean for (const string& ev : {"_sameEv", "_mixedEv"}) { TH1D* pp = fHMean.H1Clone("hMinvCombPP" + cat + ev, step); @@ -685,38 +692,38 @@ void LmvmDrawAll::CalcCombBGHistos() TH1D* k = fHMean.CreateHByClone<TH1D>("hMinvCombPM" + cat + "_mixedEv", "hMinvCombK" + cat, step); k->Divide(fHMean.H1("hMinvCombGeom" + cat + "_mixedEv", step)); k->Scale(0.5); - + // calculate combinatorial BG from same (<= 0.3 GeV) and mixed (> 0.3 GeV) events data - // at first, calculate norm. factor between same and mixed events in 300 - 1000 MeV/c2 interval + // at first, calculate norm. factor between same and mixed events in 300 - 1000 MeV/c2 interval TH1D* hGeomSame = fHMean.H1("hMinvCombGeom" + cat + "_sameEv", step); TH1D* hGeomMixed = fHMean.H1("hMinvCombGeom" + cat + "_mixedEv", step); int minBin = hGeomSame->FindBin(0.3); int maxBin = hGeomSame->FindBin(1.0); double normGM = hGeomSame->Integral(minBin, maxBin) / hGeomMixed->Integral(minBin, maxBin); - // combinatorial BG - TH1D* hBc = fHMean.CreateHByClone<TH1D>("hMinvCombGeom" + cat + "_sameEv", "hMinvCombBg" + cat, step); + // combinatorial BG + TH1D* hBc = fHMean.CreateHByClone<TH1D>("hMinvCombGeom" + cat + "_sameEv", "hMinvCombBg" + cat, step); for (int i = 0; i <= hBc->GetNbinsX(); i++) { - double val = (i <= minBin) ? 2 * fHMean.H1("hMinvCombK" + cat, step)->GetBinContent(i) - * fHMean.H1("hMinvCombGeom" + cat + "_sameEv", step)->GetBinContent(i) - : normGM * fHMean.H1("hMinvCombPM" + cat + "_mixedEv", step)->GetBinContent(i); + double val = (i <= minBin) ? 2 * fHMean.H1("hMinvCombK" + cat, step)->GetBinContent(i) + * fHMean.H1("hMinvCombGeom" + cat + "_sameEv", step)->GetBinContent(i) + : normGM * fHMean.H1("hMinvCombPM" + cat + "_mixedEv", step)->GetBinContent(i); hBc->SetBinContent(i, val); } - + // calculate signal from CB subtraction - // from 'N+-same': Signal = N+-same - B_{C} + // from 'N+-same': Signal = N+-same - B_{C} TH1D* hSigNpm = fHMean.CreateHByClone<TH1D>("hMinvCombPM" + cat + "_sameEv", "hMinvCombSignalNpm" + cat, step); hSigNpm->Add(fHMean.H1("hMinvCombBg" + cat, step), -1.); - // from MC 'Cocktail + BG': Signal = Coc+BG - B_{C} + // from MC 'Cocktail + BG': Signal = Coc+BG - B_{C} TH1D* hSigCoc = fHMean.CreateHByClone<TH1D>("hMinv_bg", "hMinvCombSignalMc" + cat, step); hSigCoc->Add(GetCocktailMinvH1(step)); hSigCoc->Add(fHMean.H1("hMinvCombBg" + cat, step), -1.); // calculate errors via error propagation formula int nofEvents = GetNofTotalEvents(); - double bW = fHMean.H1("hMinvCombPM" + cat + "_sameEv", step)->GetBinWidth(1); - int nofBins = fHMean.H1("hMinvCombPM" + cat + "_sameEv", step)->GetNbinsX(); + double bW = fHMean.H1("hMinvCombPM" + cat + "_sameEv", step)->GetBinWidth(1); + int nofBins = fHMean.H1("hMinvCombPM" + cat + "_sameEv", step)->GetNbinsX(); for (int i = 1; i <= nofBins; i++) { //s_ for same, m_ for mixed double s_pm = fHMean.H1("hMinvCombPM" + cat + "_sameEv", step)->GetBinContent(i) * nofEvents * bW; @@ -725,36 +732,36 @@ void LmvmDrawAll::CalcCombBGHistos() double m_pm = fHMean.H1("hMinvCombPM" + cat + "_mixedEv", step)->GetBinContent(i) * nofEvents * bW; double m_pp = fHMean.H1("hMinvCombPP" + cat + "_mixedEv", step)->GetBinContent(i) * nofEvents * bW; double m_mm = fHMean.H1("hMinvCombMM" + cat + "_mixedEv", step)->GetBinContent(i) * nofEvents * bW; - - // derivatives of B_{C} w.r. to according value + + // derivatives of B_{C} w.r. to according value 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 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); - // contributions to error propagation + // contributions to error propagation 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); - // final error propagation values + // final error propagation values 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); if (i <= minBin) fHMean.H1("hMinvCombBg" + cat, step)->SetBinError(i, errorBc); - if (i > minBin) fHMean.H1("hMinvCombBg" + cat, step)->SetBinError(i, errorBc2); + if (i > minBin) fHMean.H1("hMinvCombBg" + cat, step)->SetBinError(i, errorBc2); if (i <= minBin) fHMean.H1("hMinvCombSignalNpm" + cat, step)->SetBinError(i, errorSig); - if (i > minBin) fHMean.H1("hMinvCombSignalNpm" + cat, step)->SetBinError(i, errorSig2); + if (i > minBin) fHMean.H1("hMinvCombSignalNpm" + cat, step)->SetBinError(i, errorSig2); if (i <= minBin) fHMean.H1("hMinvCombSignalMc" + cat, step)->SetBinError(i, errorSig); - if (i > minBin) fHMean.H1("hMinvCombSignalMc" + cat, step)->SetBinError(i, errorSig2); + if (i > minBin) fHMean.H1("hMinvCombSignalMc" + cat, step)->SetBinError(i, errorSig2); } // error calculation - } // cat ("" or "_urqmd") - } // steps + } // cat ("" or "_urqmd") + } // steps } void LmvmDrawAll::CalcCutEffRange(double minMinv, double maxMinv) diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h b/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h index 7b8f55935fe94ccc566bbdeeb22978b8a4863b1a..29df08cf231f8cc0f93beadc7f8801a445a5c373 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h @@ -23,7 +23,8 @@ 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 10 * 4.52941e-4; // TODO urgent: delete factor 10 and set right parameter!!! + 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 54fa1e44046ad5e4f2dbeb733f0a1cf97247a03e..092452b9b99b5f6de8ef0c4ea5773023e5e4abf9 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx @@ -59,15 +59,18 @@ LmvmTask::~LmvmTask() {} void LmvmTask::InitHists() { string ax = "Yield"; - string axMinv = "Yield"; // "dN/dM_{ee}/N [GeV/c^{2}]^{-1}"; TODO: when in Draw.cxx the scaling to bin width is implemented this can be changed back to "dN/dM..." + string axMinv = + "Yield"; // "dN/dM_{ee}/N [GeV/c^{2}]^{-1}"; TODO: when in Draw.cxx the scaling to bin width is implemented this can be changed back to "dN/dM..." fH.CreateH2("hMomVsAnglePairSignalMc", "#sqrt{P_{e^{#pm}} P_{e^{#mp}}} [GeV/c]", "#theta_{e^{+},e^{-}} [deg]", "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("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); @@ -104,8 +107,7 @@ void LmvmTask::InitHists() 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.); + 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.); @@ -143,8 +145,9 @@ void LmvmTask::InitHists() 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}]", axMinv, 2500, 0., 2.5); - fH.CreateH1("hMinvBgSource2_elid", {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"}, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0., 2.5); // "pi" are misid. charged pions - + fH.CreateH1("hMinvBgSource2_elid", {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"}, + "M_{ee} [GeV/c^{2}]", axMinv, 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); @@ -404,10 +407,10 @@ bool LmvmTask::IsMcTrackAccepted(int mcTrackInd) void LmvmTask::FillAccRecVsMomHist() { - for (const int& pdg : {11, 211, 2212} ) { // TODO: restructure this more efficiently (wihtout these pdg-loops) + 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"; - + 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); @@ -461,14 +464,17 @@ void LmvmTask::FillAccRecVsMomHist() } } } - + // 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"; + string hName = (pdg == 211) ? "hCandMisIdAsEl_pi" + : (pdg == 111) ? "hCandMisIdAsEl_pi0" + : (pdg == 2212) ? "hCandMisIdAsEl_proton" + : "hCandMisIdAsEl_kaon"; fH.FillH1(hName, cand.fMomentum.Mag()); } } @@ -620,62 +626,74 @@ void LmvmTask::CombinatorialPairs() const auto& cand1 = fCandsTotal[iC1]; for (size_t iC2 = iC1 + 1; iC2 < nCand; iC2++) { - const auto& cand2 = fCandsTotal[iC2]; + const auto& cand2 = fCandsTotal[iC2]; 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) { + + for (auto step : fH.fAnaSteps) { // seperate PLUTO and UrQMD electrons - if (cand1.IsCutTill(step) && cand2.IsCutTill(step) && std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11) { + if (cand1.IsCutTill(step) && cand2.IsCutTill(step) && std::abs(cand1.fMcPdg) == 11 + && std::abs(cand2.fMcPdg) == 11) { // only PLUTO electrons if (cand1.IsMcSignal() && cand2.IsMcSignal()) { if (isSameEvent) { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_pluto_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_pluto_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_pluto_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_pluto_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_pluto_sameEv", step, pRec.fMinv, w); } else { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_pluto_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_pluto_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_pluto_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_pluto_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_pluto_mixedEv", step, pRec.fMinv, w); } } // only UrQMD electrons else if (!cand1.IsMcSignal() && !cand2.IsMcSignal()) { if (isSameEvent) { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmd_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmd_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmd_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_urqmd_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_urqmd_sameEv", step, pRec.fMinv, w); } else { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmd_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmd_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmd_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_urqmd_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_urqmd_mixedEv", step, pRec.fMinv, w); } } } - + // for common CB, don't use PLUTO signals - if (cand1.IsMcSignal() or cand2.IsMcSignal()) continue; + if (cand1.IsMcSignal() or cand2.IsMcSignal()) continue; if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue; if (cand1.IsCutTill(step) && cand2.IsCutTill(step)) { if (cand1.fCharge * cand2.fCharge < 0) { if (isSameEvent) fH.FillH1("hMinvCombPM_sameEv", step, pRec.fMinv, w); - else fH.FillH1("hMinvCombPM_mixedEv", step, pRec.fMinv, w); + else + fH.FillH1("hMinvCombPM_mixedEv", step, pRec.fMinv, w); } if (cand1.fCharge < 0 && cand2.fCharge < 0) { if (isSameEvent) fH.FillH1("hMinvCombMM_sameEv", step, pRec.fMinv, w); - else fH.FillH1("hMinvCombMM_mixedEv", step, pRec.fMinv, w); + else + fH.FillH1("hMinvCombMM_mixedEv", step, pRec.fMinv, w); } if (cand1.fCharge > 0 && cand2.fCharge > 0) { if (isSameEvent) fH.FillH1("hMinvCombPP_sameEv", step, pRec.fMinv, w); - else fH.FillH1("hMinvCombPP_mixedEv", step, pRec.fMinv, w); + else + fH.FillH1("hMinvCombPP_mixedEv", step, pRec.fMinv, w); } } - } // steps - } // cand2 - } // cand1 + } // steps + } // cand2 + } // cand1 } void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands) @@ -751,56 +769,71 @@ void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAna ELmvmBgPairSrc bgSrc = LmvmUtils::GetBgPairSrc(candP, candM); //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); - } + 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 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 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 neither electron nor misid. charged pion + // 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 + "pipi0", parRec.fMinv); - else fH.FillH1(hName + "oo", 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); + } + } //} } } @@ -815,12 +848,17 @@ void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg) 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; + 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 { @@ -845,12 +883,18 @@ 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); } @@ -861,8 +905,22 @@ void LmvmTask::BgPairPdg(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaS 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; + 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); } @@ -917,7 +975,7 @@ void LmvmTask::SignalAndBgReco() for (auto step : fH.fAnaSteps) { if (cand.IsCutTill(step)) { TrackSource(cand, step, pdg); - fH.FillH2("hPtYCandidate", step, cand.fRapidity, cand.fMomentum.Perp(), fW); + fH.FillH2("hPtYCandidate", step, cand.fRapidity, cand.fMomentum.Perp(), fW); } } } @@ -1163,7 +1221,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 + 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) { @@ -1173,9 +1231,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 - + 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.; @@ -1220,8 +1278,8 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s else if (mvdStationNum == 2) cand.fIsMvd2Cut = isMvdCut; - cout << "cand.fIsMvd1Cut = " <<cand.fIsMvd1Cut << endl; // TODO: delete these cout lines - cout << "cand.fIsMvd2Cut = " <<cand.fIsMvd2Cut << endl; + cout << "cand.fIsMvd1Cut = " << cand.fIsMvd1Cut << endl; // TODO: delete these cout lines + cout << "cand.fIsMvd2Cut = " << cand.fIsMvd2Cut << endl; } } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx index 347c802c19d24ae7909c1ce877707fa6aa9193de..2d7dde947dc2b37e0039ca2f94fddbaa2694f2b8 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx @@ -241,13 +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 = (momentum < 1.) ? true : CbmLitGlobalElectronId::GetInstance().IsTrdElectron(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/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx b/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx index f949a87d25bcdf325353d5bb43d2f2f5f4dde96b..f8a244e9e518daf2134999b800ccaa7a8bc20527 100644 --- a/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx +++ b/reco/littrack/cbm/elid/CbmLitGlobalElectronId.cxx @@ -30,8 +30,9 @@ #include <cmath> CbmLitGlobalElectronId::CbmLitGlobalElectronId() - : 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 + : 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,10 +101,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 (momentum < 1.) return true; + if (ann > fTrdAnnCut) return true; else return false; @@ -121,7 +122,7 @@ Bool_t CbmLitGlobalElectronId::IsTofElectron(Int_t globalTrackIndex, Double_t mo else { if (mass2 < 0.01) { return true; } } - + // stricter cut /*if (momentum >= 0.8) { if (mass2 < (- 0.020 * momentum + 0.026)) { return true; } @@ -129,10 +130,10 @@ Bool_t CbmLitGlobalElectronId::IsTofElectron(Int_t globalTrackIndex, Double_t mo 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 90788c14ee3386658019a068c777c5ce050c02a6..38c1b08a6ed04d499e543ac538396dfca99709cd 100644 --- a/reco/littrack/cbm/qa/tracking/CbmLitTrackingQa.cxx +++ b/reco/littrack/cbm/qa/tracking/CbmLitTrackingQa.cxx @@ -93,7 +93,7 @@ CbmLitTrackingQa::CbmLitTrackingQa() , fTofPoints(NULL) , fTofMatches(NULL) , 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 + , fTrdAnnCut(0.8) // at the moment, TrdAnn is likelihood value // values loose | strict: 0.2 | 0.8 { }