From 8eb40c5a8c7a798a3b805464409688045e2242b1 Mon Sep 17 00:00:00 2001 From: Cornelius Feier-Riesen <cornelius.riesen@physik.uni-giessen.de> Date: Fri, 21 Oct 2022 15:09:55 +0200 Subject: [PATCH] Update LMVM analysis --- .../PWGDIL/dielectron/lmvm/CMakeLists.txt | 1 + analysis/PWGDIL/dielectron/lmvm/LmvmCand.h | 20 +- analysis/PWGDIL/dielectron/lmvm/LmvmDef.h | 13 + analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx | 133 +- analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h | 20 +- .../PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx | 1684 ++++++++++++++--- analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h | 46 +- analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx | 269 ++- analysis/PWGDIL/dielectron/lmvm/LmvmHist.h | 25 +- .../PWGDIL/dielectron/lmvm/LmvmSimParam.h | 3 +- analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx | 1166 +++++++++--- analysis/PWGDIL/dielectron/lmvm/LmvmTask.h | 35 +- analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx | 82 +- analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h | 71 +- macro/analysis/dielectron/batch_job.py | 6 +- macro/analysis/dielectron/batch_send.py | 2 +- macro/analysis/dielectron/draw_all.py | 4 +- 17 files changed, 2847 insertions(+), 733 deletions(-) mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmHist.h mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmTask.h mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx mode change 100644 => 100755 analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h diff --git a/analysis/PWGDIL/dielectron/lmvm/CMakeLists.txt b/analysis/PWGDIL/dielectron/lmvm/CMakeLists.txt index 012ae48099..89fd4f5206 100644 --- a/analysis/PWGDIL/dielectron/lmvm/CMakeLists.txt +++ b/analysis/PWGDIL/dielectron/lmvm/CMakeLists.txt @@ -30,6 +30,7 @@ set(LIBRARY_NAME CbmLmvm) set(LINKDEF LmvmLinkDef.h) set(PUBLIC_DEPENDENCIES CbmBase + CbmRichBase CbmData KF NicaCbmFormat diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h index 65a54829c5..575bef1faa 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h @@ -49,11 +49,23 @@ public: TVector3 fPosition; TVector3 fMomentum; double fMass = 0.; + double fMassSig = -1.; // mass of signal mother (remains '-1' if candidate is not signal electron) double fEnergy = 0.; double fRapidity = 0.; int fCharge = 0; double fChi2Prim = 0.; - double fChi2sts = 0.; + double fChi2Sts = 0.; + double fChi2Rich = 0.; + double fChi2Trd = 0.; + double fChi2Tof = 0.; + double fLength = 0.; // length of according global track + double fTime = 0.; // + double fTofDist = -1.; // Hit-Track-Distance in ToF + + // To investigate misidentifications in single sub detectors + bool fIsRichElectron = false; + bool fIsTrdElectron = false; + bool fIsTofElectron = false; int fMcMotherId = -1; int fEventNumber = 0; @@ -64,8 +76,10 @@ public: int fStsInd = -1; int fRichInd = -1; int fTrdInd = -1; - int fTofInd = -1; + int fTofInd = -1;// TODO: change to "fTofHitInd" + int fTofTrackInd = -1; int fMcPdg = -1; + int fGTrackInd = -1; double fRichAnn = 0.; double fTrdAnn = 0.; double fMass2 = 0.; @@ -75,7 +89,7 @@ public: // Cuts. If true then cut is passed bool fIsChi2Prim = false; bool fIsElectron = false; - bool fIsGammaCut = true; // Will be set to 'false' as soon as a partner with minv < 25 MeV is found + 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/LmvmDef.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h index 1fe0053cc8..809283fd9a 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h @@ -123,5 +123,18 @@ public: double fFitSigma = 0; }; +class LmvmLegend { +public: + LmvmLegend() {} + LmvmLegend(TH1D* h, const char* name, Option_t* opt) + : fH(h) + , fName(name) + , fOpt(opt) + { + } + TH1D* fH = nullptr; + const char* fName = ""; + Option_t* fOpt = "l"; +}; #endif diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx index 12ff21e3e9..a4b7d5173c 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx @@ -56,9 +56,12 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir, fH.fHM.ScaleByPattern(".*", 1. / fNofEvents); RebinMinvHist(); - //DrawAnaStepMany("pair_pty", [this](ELmvmAnaStep step) { DrawPtY(step); }); DrawAnaStepMany("pty_pair_signal", [this](ELmvmAnaStep step) { DrawPtY("hPtYPairSignal", step); }); - DrawAnaStepMany("pty_cand", [this](ELmvmAnaStep step) { DrawPtY("hPtYCandidate", step); }); + for (const auto& cand : fH.fCandNames) { + const string& cName = "hPtY_cands/" + cand; + const string& hName = "hPtY_cands_" + cand; + DrawAnaStepMany(cName, [this, hName](ELmvmAnaStep step) { DrawPtY(hName, step); }); + } DrawAnaStepMany("pair_rapidity", [this](ELmvmAnaStep step) { DrawRapidity(step); }); //DrawAnaStepMany("pair_pty_efficiency", [this](ELmvmAnaStep step) { DrawPtYEfficiency(step); }); // TODO: causes segmentation violation error DrawAnaStepMany("minv_sbg", [this](ELmvmAnaStep step) { DrawMinvSBg(step); }); @@ -82,9 +85,10 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir, DrawMismatchesAndGhosts(); DrawMvdCutQa(); DrawMvdAndStsHist(); - DrawAccRecMom(); + DrawAccRecVsMom(); // TODO: finish this method! DrawPmtXY(); DrawMinvBg(); // TODO: do not extra method + DrawBetaMomSpectra(); DrawCombinatorialPairs(); SaveCanvasToImage(); @@ -118,70 +122,16 @@ void LmvmDraw::RebinMinvHist() fH.Rebin("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, nGroupBgSrc); } -TH1D* LmvmDraw::CreateSignificanceH1(TH1D* s, TH1D* bg, const string& name, const string& option) +void LmvmDraw::DrawBetaMomSpectra() { - int nBins = s->GetNbinsX(); - TH1D* hsig = new TH1D(name.c_str(), name.c_str(), nBins, s->GetXaxis()->GetXmin(), s->GetXaxis()->GetXmax()); - hsig->GetXaxis()->SetTitle(s->GetXaxis()->GetTitle()); - - // "right" - reject right part of the histogram. value > cut -> reject - if (option == "right") { - for (int i = 1; i <= nBins; i++) { - double sumSignal = s->Integral(1, i, "width"); - double sumBg = bg->Integral(1, i, "width"); - double sign = (sumSignal + sumBg != 0.) ? sumSignal / std::sqrt(sumSignal + sumBg) : 0.; - hsig->SetBinContent(i, sign); - hsig->GetYaxis()->SetTitle("Significance"); - } - } - // "left" - reject left part of the histogram. value < cut -> reject - else if (option == "left") { - for (int i = nBins; i >= 1; i--) { - double sumSignal = s->Integral(i, nBins, "width"); - double sumBg = bg->Integral(i, nBins, "width"); - double sign = (sumSignal + sumBg != 0.) ? sumSignal / std::sqrt(sumSignal + sumBg) : 0.; - hsig->SetBinContent(i, sign); - hsig->GetYaxis()->SetTitle("Significance"); - } - } - return hsig; -} - -TH2D* LmvmDraw::CreateSignificanceH2(TH2D* signal, TH2D* bg, const string& name, const string& title) -{ - double xmin = 1.0; - double xmax = 5.0; - double ymin = 1.0; - double ymax = 5.0; - double delta = 0.1; - int nStepsX = (int) ((xmax - xmin) / delta); - int nStepsY = (int) ((ymax - ymin) / delta); - - TH2D* hsig = new TH2D(name.c_str(), title.c_str(), nStepsX, xmin, xmax, nStepsY, ymin, ymax); - - int binX = 1; - for (double xcut = xmin; xcut <= xmax; xcut += delta, binX++) { - int binY = 1; - for (double ycut = ymin; ycut <= ymax; ycut += delta, binY++) { - double sumSignal = 0; - double sumBg = 0; - for (int ix = 1; ix <= signal->GetNbinsX(); ix++) { - for (int iy = 1; iy <= signal->GetNbinsY(); iy++) { - double xc = signal->GetXaxis()->GetBinCenter(ix); - double yc = signal->GetYaxis()->GetBinCenter(iy); - double val = -1 * (ycut / xcut) * xc + ycut; - - if (!(xc < xcut && yc < val)) { - sumSignal += signal->GetBinContent(ix, iy); - sumBg += bg->GetBinContent(ix, iy); - } - } - } - double sign = (sumSignal + sumBg != 0.) ? sumSignal / std::sqrt(sumSignal + sumBg) : 0.; - hsig->SetBinContent(binX, binY, sign); - } - } - return hsig; + TH2D* hPos = fH.H2Clone("hBetaMom_cands_plutoEl+"); + TH2D* hEl = fH.H2Clone("hBetaMom_cands_plutoEl-"); + TCanvas* c = fH.fHM.CreateCanvas("betaMom/", "betaMom/", 1600, 800); + c->Divide(2,1); + c->cd(1); + DrawH2(hEl, kLinear, kLinear, kLog, "colz"); + c->cd(2); + DrawH2(hPos, kLinear, kLinear, kLog, "colz"); } void LmvmDraw::DrawCutEffH1(const string& hist, const string& option) @@ -440,8 +390,7 @@ void LmvmDraw::Draw1DCut(const string& hist, const string& sigOption, double cut DrawCutEffH1(hist, sigOption); c->cd(3); - TH1D* sign = - CreateSignificanceH1(fH.H1(hist, ELmvmSrc::Signal), fH.H1(hist, ELmvmSrc::Bg), hist + "_significance", sigOption); + TH1D* sign = fH.CreateSignificanceH1(fH.H1(hist, ELmvmSrc::Signal), fH.H1(hist, ELmvmSrc::Bg), hist + "_significance", sigOption); DrawH1(sign, kLinear, kLinear, "hist"); } @@ -639,37 +588,21 @@ void LmvmDraw::DrawElPurity() "Tracks per event x10^{3}"); // Occurency of Electrons and not-Electrons for various cut categories - for (const string& hName : {"hAnnRichVsMomPur", "hTrdElLikePur", "hTofM2Pur"}) { + for (const string& hName : {"hAnnRichVsMomPur", "hTrdElLikePur"}) { string cName = "purity/cuts_" + hName; TCanvas* c = fH.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 800); c->Divide(3, 1); int hi = 1; for (const string& id : {"El", "Bg"}) { TH2D* hist = fH.H2Clone(hName + "_" + id); - if (hName == "hTofM2Pur") { - hist->GetXaxis()->SetRangeUser(0., 3.); - hist->GetYaxis()->SetRangeUser(-0.1, 0.05); - } c->cd(hi); DrawH2(hist, kLinear, kLinear, kLog, "COLZ"); DrawTextOnPad(id.c_str(), 0.6, 0.89, 0.7, 0.99); - - if (hName == "hTofM2Pur") { // TODO: these drawn lines needed? - vector<TLine*> lines {new TLine(0., 0.01, 1.3, 0.01), new TLine(1.3, 0.01, 2.5, 0.022)}; // set by hand - for (size_t i = 0; i < lines.size(); i++) { - lines[i]->SetLineWidth(2.); - lines[i]->Draw(); - } - } hi++; } c->cd(hi); TH2D* ratio = fH.H2Clone(hName + "_El"); ratio->Divide(fH.H2(hName + "_Bg")); - if (hName == "hTofM2Pur") { - ratio->GetXaxis()->SetRangeUser(0., 3.); - ratio->GetYaxis()->SetRangeUser(-0.1, 0.05); - } DrawH2(ratio, kLinear, kLinear, kLog, "COLZ"); DrawTextOnPad("Ratio El/Bg", 0.4, 0.85, 0.8, 0.99); } @@ -860,26 +793,26 @@ void LmvmDraw::DrawMinvMatching(ELmvmAnaStep step) fH.DrawAnaStepOnPad(step); } -void LmvmDraw::DrawAccRecMom() +void LmvmDraw::DrawAccRecVsMom() { - // Acceptance and reconstruction yields for various detector combinations - for (const int& pdg : {11, 211, 2212}) { + // Acceptance and reconstruction yields cs. momentum for various detector combinations + for (const int& pdg : {11, 211, 2212, 321}) { 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; + string ptcl = (pdg == 11) ? "hEl" : (pdg == 211) ? "hPi" : (pdg == 2212) ? "hProton" : "hKaon"; + vector<TH1*> histsAll, histsPrim; int i = 0; + for (const string& subName : subNames) { - TH1D* hAll = fH.H1(hName + "_all_" + subName); + TH1D* hAll = fH.H1(ptcl + "Mom_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(ptcl + "Mom_prim_" + subName); hPrim->SetMinimum(3e-6); hPrim->SetMaximum(50); latexPrim[i] = latex[i] + " (" + Cbm::NumberToString(hPrim->GetEntries() / fNofEvents, 2) + "/ev.)"; @@ -887,8 +820,10 @@ void LmvmDraw::DrawAccRecMom() i++; } - double y1 = 0.17; //(pdg == 211) ? 0.20 : 0.74; + //if (pdg == 321) continue; TODO: with kaons? + double y1 = 0.17; //(pdg == 211) ? 0.20 : 0.74; double y2 = 0.42; //(pdg == 211) ? 0.45 : 0.99; + string cName = "AccRecMom/" + ptcl + "Mom"; fH.fHM.CreateCanvas(cName, cName, 900, 900); DrawH1(histsAll, latexAll, kLinear, kLog, true, 0.4, y1, 0.95, y2, "hist"); @@ -896,14 +831,6 @@ void LmvmDraw::DrawAccRecMom() DrawH1(histsPrim, latexPrim, kLinear, kLog, true, 0.4, y1, 0.95, y2, "hist"); } - // Pions, protons and kaons that are misidentified as electrons - TH1* misPi = fH.H1Clone("hCandMisIdAsEl_pi"); - TH1* misProt = fH.H1Clone("hCandMisIdAsEl_proton"); - TH1* misKa = fH.H1Clone("hCandMisIdAsEl_kaon"); - - fH.fHM.CreateCanvas("AccRecMom/MisidCands", "AccRecMom/MisidCands", 800, 800); - DrawH1({misPi, misProt, misKa}, {"#pi^{#pm}", "proton", "kaon"}, kLinear, kLog, true, 0.89, 0.75, 0.99, 0.92, "p"); - // Acceptance in single detectors for (const string& det : {"sts", "rich", "trd", "tof"}) { vector<TH1*> hVec; diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h index 362907dde0..ec5f134242 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.h @@ -45,7 +45,7 @@ public: private: Int_t fNofEvents = 0; // number of events of current job - bool fUseMvd = false; // do you want to draw histograms related to the MVD detector? + bool fUseMvd = false; // do you want to draw histograms related to the MVD detector? LmvmCuts fCuts; // electron identification and analysis cuts @@ -63,20 +63,6 @@ private: */ void SaveCanvasToImage(); - /** - * Produce 1D significance histogram Significance=S/sqrt(S+BG). - * \param[in] s Histogram with signal. - * \param[in] bg Histogram eith background. - * \param[in] name Name of new significance histogram. - * \param[in] option Could be "right" or "left". - */ - TH1D* CreateSignificanceH1(TH1D* s, TH1D* bg, const std::string& name, const std::string& option); - - /** - * Produce 2D significance histogram Significance=S/sqrt(S+BG). - */ - TH2D* CreateSignificanceH2(TH2D* signal, TH2D* bg, const std::string& name, const std::string& title); - void DrawCutEffH1(const std::string& hist, const std::string& option); void DrawAnaStepMany(const std::string& cName, std::function<void(ELmvmAnaStep)> drawFunc); @@ -104,7 +90,9 @@ private: void DrawBgSourcePairs(ELmvmAnaStep step, bool inPercent, bool drawAnaStep = true); void DrawBgSourcePairsAll(); - void DrawAccRecMom(); + void DrawAccRecVsMom(); + + void DrawBetaMomSpectra(); void Draw2DCutTriangle(double xCr, double yCr); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx old mode 100644 new mode 100755 index dea2df01f3..9828b8866b --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx @@ -19,6 +19,8 @@ #include "TH1D.h" #include "TH2D.h" #include "TKey.h" +#include "TLine.h" +#include "TGraph.h" #include "TLatex.h" #include "TMath.h" #include "TStyle.h" @@ -31,6 +33,7 @@ #include <string> #include "LmvmDef.h" +#include "LmvmUtils.h" using namespace std; using namespace Cbm; @@ -64,14 +67,27 @@ void LmvmDrawAll::DrawHistFromFile(const string& fileInmed, const string& fileQg CalcCutEffRange(0.0, 0.2); CalcCutEffRange(0.2, 0.6); CalcCutEffRange(0.6, 1.2); + DrawCutEffSignal(); + DrawPionSuppression(); CalcCombBGHistos(); SBgRangeAll(); DrawSBgResults(); DrawMinvAll(); DrawMinvCombBgAndSignal(); + DrawMinvOfficialStyle(); DrawMinvPtAll(); DrawMinvBgSourcesAll(); DrawSBgVsMinv(); + InvestigateMisid(); // TODO: move some procedures in here into seperate methods? + DrawBetaMomSpectra(); + DrawMomRecoPrecision(); + DrawMomPluto(); + DrawPurity(); + DrawSignificancesAll(); + DrawTofM2(); + DrawGTrackVertex(); + DrawTofHitXY(); + DrawMinvScaleValues(); SaveHist(); SaveCanvasToImage(); @@ -109,33 +125,131 @@ void LmvmDrawAll::CreateMeanHistAll() { int nofEvents = GetNofTotalEvents(); + // Global Track Loop + for (const string& ptcl : fHMean.fGTrackNames) { + CreateMeanHist<TH1D>("hMom_gTracks_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hPtY_gTracks_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hTofM2_gTracks_" + ptcl, nofEvents); + /*CreateMeanHist<TH2D>("hIndexStsRich_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hIndexStsTrd_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hIndexStsTof_" + ptcl, nofEvents);*/ + CreateMeanHist<TH2D>("hRichRingTrackDist_gTracks_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hMatchId_gTracks_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofM2Calc_gTracks_" + ptcl, nofEvents); + + for (const string& cat : {"trueid", "misid", "truematch", "mismatch"}) { + CreateMeanHist<TH2D>("hTofHitXY_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofPointXY_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hTofHitPointDist_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofTimeVsMom_gTracks_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofHitTrackDist_gTracks_" + cat + "_" + ptcl, nofEvents); + + } + + for (const string& cat : {"all", "complete"}) { + CreateMeanHist<TH1D>("hNofMismatches_gTracks_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hNofMismatchedTrackSegments_" + cat + "_" + ptcl, nofEvents); + for (const string& match : {"truematch", "mismatch"}) { + for (const string& det : {"rich", "trd", "tof"}) { + CreateMeanHist<TH1D>("hChi2_" + match + "_" + cat + "_" + det + "_" + ptcl, nofEvents); + } + } + } + } // fGTrackNames + + // Candidate Loop + for (const string& ptcl : fHMean.fCandNames) { + CreateMeanHist<TH2D>("hTofPileHitXY_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofPilePointXY_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hTofPileHitPointDist_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hTofPilePty_cands_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofTimeVsMom_cands_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofHitTrackDist_cands_" + ptcl, nofEvents); + + for (auto step : fHMean.fAnaSteps) { + CreateMeanHist<TH1D>(fHMean.GetName("hMom_cands_" + ptcl, step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hMomRatio_cands_" + ptcl, step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hMomRatioVsMom_cands_" + ptcl, step), 1); + CreateMeanHist<TH2D>(fHMean.GetName("hPtY_cands_" + ptcl, step), nofEvents); + CreateMeanHist<TH2D>(fHMean.GetName("hTofM2_cands_" + ptcl, step), nofEvents); + CreateMeanHist<TH2D>(fHMean.GetName("hRichRingTrackDist_cands_" + ptcl, step), nofEvents); + } // steps + + for (const string& cat : {"gTracks", "cands"}) { + CreateMeanHist<TH1D>("hBetaMom_" + cat + "_" + ptcl, nofEvents); + } + + + for (const string det : {"sts", "rich", "trd", "tof"}) { + CreateMeanHist<TH2D>("hChi2VsMom_" + det + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofTimeVsChi2_" + det + "_" + ptcl, nofEvents); + } + + for (const string det : {"StsRich", "StsTrd", "RichTrd"}) { + CreateMeanHist<TH2D>("hChi2Comb_" + det + "_" + ptcl, nofEvents); + } + + + } // fCandNames + + // Step Loop for (auto step : fHMean.fAnaSteps) { for (auto src : {ELmvmSrc::Bg, ELmvmSrc::Eta, ELmvmSrc::Pi0}) { CreateMeanHist<TH1D>(fHMean.GetName("hMinv", src, step), nofEvents, fRebMinv); + CreateMeanHist<TH1D>(fHMean.GetName("hMinv_urqmdAll", src, step), nofEvents, fRebMinv); + CreateMeanHist<TH1D>(fHMean.GetName("hMinv_urqmdEl", src, step), nofEvents, fRebMinv); CreateMeanHist<TH2D>(fHMean.GetName("hMinvPt", src, step), nofEvents); } for (const string& comb : {"PM", "PP", "MM"}) { - for (const string& cat : {"", "_urqmd"}) { + for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { for (const string& ev : {"_sameEv", "_mixedEv"}) { CreateMeanHist<TH1D>(fHMean.GetName("hMinvComb" + comb + cat + ev, step), nofEvents, fRebMinv); } } } + + CreateMeanHist<TH1D>(fHMean.GetName("hPionSupp_idEl", step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hPionSupp_pi", step), nofEvents); + CreateMeanHist<TH2D>(fHMean.GetName("hCandPdgVsMom", step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hTofPilePdgs_cands", step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hVertexXZ_misidTof", step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hVertexYZ_misidTof", step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hVertexXY_misidTof", step), nofEvents); + CreateMeanHist<TH1D>(fHMean.GetName("hVertexRZ_misidTof", step), nofEvents); + } // steps + + string hBgSrc = "hMinvBgSource2_elid_"; + for (const string& pair : {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"}) { + string hPairFull = hBgSrc + pair; + CreateMeanHist<TH1D>(hPairFull, nofEvents, fRebMinv); + } + + for (const string& cat : {"true", "misid"}) { + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { + CreateMeanHist<TH1D>("hMom_" + fHMean.fCandNames[iP] + "_" + cat, nofEvents); + CreateMeanHist<TH1D>("hPtY_" + fHMean.fCandNames[iP] + "_" + cat, nofEvents); + CreateMeanHist<TH1D>("hTofM2_" + fHMean.fCandNames[iP] + "_" + cat, nofEvents); + } + } + + for (const string& det : {"rich", "trd", "tof"}) { + for (const string& cat : {"all", "complete"}) { + CreateMeanHist<TH2D>("hPdgVsMom_gTracks_" + det + "_" + cat, nofEvents); + } + } + + for (const string& cat : {"El", "Bg"}) { + CreateMeanHist<TH2D>("hAnnRichVsMomPur_" + cat, nofEvents); + CreateMeanHist<TH2D>("hTrdElLikePur_" + cat, nofEvents); } - CreateMeanHist<TH1D>("hMinvBgSource2_elid_gg", nofEvents, fRebMinv); // TODO: automatize this - CreateMeanHist<TH1D>("hMinvBgSource2_elid_pipi", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_pi0pi0", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_oo", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_gpi", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_gpi0", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_go", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_pipi0", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_pio", nofEvents, fRebMinv); - CreateMeanHist<TH1D>("hMinvBgSource2_elid_pi0o", nofEvents, fRebMinv); + + CreateMeanHist<TH1D>("hBetaMom_allGTracks", nofEvents); + CreateMeanHist<TH1D>("hTofPilePdgs_gTracks", nofEvents); + CreateMeanHist<TH2D>("hVertexGTrackRZ", nofEvents); } -TH1D* LmvmDrawAll::GetCocktailMinvH1(ELmvmAnaStep step) { return GetCocktailMinv<TH1D>("hMinv", step); } +TH1D* LmvmDrawAll::GetCocktailMinvH1(const string& name, ELmvmAnaStep step) { return GetCocktailMinv<TH1D>(name, step); } template<class T> T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step) @@ -156,8 +270,7 @@ T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step) sHist->Scale(1. / binWidth); } if (cocktail == nullptr) cocktail = sHist; - else - cocktail->Add(sHist); + else cocktail->Add(sHist); } cocktail->Add(sEta); cocktail->Add(sPi0); @@ -165,13 +278,896 @@ T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step) return cocktail; } +void LmvmDrawAll::DrawMinvScaleValues() +{ + TH1D* inmed = new TH1D("inmed", "inmed; M_{ee} [GeV/c^{2}]; Scale Value", 3400., 0., 3.4); + TH1D* qgp = new TH1D("qgp", "qgp; M_{ee} [GeV/c^{2}]; Scale Value", 3400., 0., 3.4); + for (int iB = 1; iB <= 3400; iB++) { + double binCenter = inmed->GetBinCenter(iB); + double sInmed = LmvmUtils::GetMassScaleInmed(binCenter); + double sQgp = LmvmUtils::GetMassScaleQgp(binCenter); + inmed->SetBinContent(iB, sInmed); + qgp->SetBinContent(iB, sQgp); + } + fHMean.fHM.CreateCanvas("MinvScaleValues", "MinvScaleValues", 2400, 800); + DrawH1({inmed, qgp}, {"inmed", "QGP"}, kLinear, kLog, true, 0.7, 0.7, 0.9, 0.9, "hist"); + + TH1D* inmed2 = (TH1D*) inmed->Clone(); + TH1D* qgp2 = (TH1D*) qgp->Clone(); + inmed2->GetXaxis()->SetRangeUser(1., 1.1); + qgp2->GetXaxis()->SetRangeUser(1., 1.1); + inmed2->GetYaxis()->SetRangeUser(3e-3, 5e-1); + qgp2->GetYaxis()->SetRangeUser(3e-3, 5e-1); + + fHMean.fHM.CreateCanvas("MinvScaleValues2", "MinvScaleValues2", 2400, 800); + DrawH1({inmed2, qgp2}, {"inmed", "QGP"}, kLinear, kLog, true, 0.7, 0.7, 0.9, 0.9, "p"); + +} + +void LmvmDrawAll::DrawSignificance(TH2D* hEl, TH2D* hBg, const string& name, double minZ, double maxZ, const string& option) +{ + string hElProjName = name + "_yProjEl"; + string hBgProjName = name + "_yProjBg"; + TH1D* el = hEl->ProjectionY(hElProjName.c_str(), 1, hEl->GetXaxis()->GetNbins(), ""); + TH1D* bg = hBg->ProjectionY(hBgProjName.c_str(), 1, hBg->GetXaxis()->GetNbins(), ""); + TH2D* rat = (TH2D*) hEl->Clone(); + rat->Divide(hBg); + + const string hName = name + "_sign"; + const string cName = "Significance/" + name; + + TH1D* sign = fHMean.CreateSignificanceH1(el, bg, hName.c_str(), option); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 1600); + c->Divide(3, 2); + c->cd(1); + DrawH2(hEl, kLinear, kLinear, kLog, "colz"); + hEl->GetZaxis()->SetRangeUser(minZ, maxZ); + DrawTextOnPad("Electrons", 0.25, 0.9, 0.75, 0.99); + c->cd(2); + DrawH2(hBg, kLinear, kLinear, kLog, "colz"); + hBg->GetZaxis()->SetRangeUser(minZ, maxZ); + DrawTextOnPad("Other", 0.25, 0.9, 0.75, 0.99); + c->cd(3); + DrawH2(rat, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("Ratio Electrons/Other", 0.2, 0.9, 0.8, 0.99); + c->cd(4); + DrawH1({el, bg}, {"Electrons", "Other"}, kLinear, kLog, true, 0.75, 0.75, 0.91, 0.91); + DrawTextOnPad("Yield / Event", 0.25, 0.9, 0.75, 0.99); + c->cd(5); + DrawH1(sign, kLinear, kLinear, "hist"); + double maxX = -999999.; + double maxY = -999999.; + for (int iB = 1; iB <= sign->GetNbinsX(); iB++) { + if (option == "right" && sign->GetBinContent(iB) >= maxY) { + maxX = sign->GetBinCenter(iB); + maxY = sign->GetBinContent(iB); + } + if (option == "left" && sign->GetBinContent(iB) > maxY) { + maxX = sign->GetBinCenter(iB); + maxY = sign->GetBinContent(iB); + } + } + DrawTextOnPad("max. at: " + Cbm::NumberToString(maxX, 1), 0.55, 0.2, 0.85, 0.5); + DrawTextOnPad("Significance", 0.25, 0.9, 0.75, 0.99); +} + +void LmvmDrawAll::DrawSignificancesAll() // TODO: implement automatization to seperate electrons from not-electrons in gTracks and cands +{ + // RICH ANN and Electron Likelihood + { + for (const string& hist : {"hAnnRichVsMomPur", "hTrdElLikePur"}) { + TH2D* el = fHMean.H2Clone(hist + "_El"); + TH2D* bg = fHMean.H2Clone(hist + "_Bg"); + DrawSignificance(el, bg, hist, 1e-9, 1e-3, "left"); + } + } + + // X2 in TRD (filled after ElID step) + { + TH2D* el = fHMean.H2Clone("hChi2VsMom_trd_plutoEl+"); + el->Add(fHMean.H2("hChi2VsMom_trd_plutoEl-")); + + TH2D* bg = fHMean.H2Clone("hChi2VsMom_trd_pion+"); + for (size_t iP = 2; iP < fHMean.fCandNames.size(); iP++) { // consider UrQMD-electrons as "BG" + string ptcl = fHMean.fCandNames[iP]; + string hName1 = "hChi2VsMom_trd_" + ptcl; + bg->Add(fHMean.H2(hName1.c_str())); + } + + DrawSignificance(el, bg, "hChi2VsMom_trd", 1e-9, 1e-3, "right"); + } + + // Time in ToF (filled after ElID step) + { + TH2D* el = fHMean.H2Clone("hTofTimeVsMom_cands_plutoEl+"); + TH2D* bg = fHMean.H2Clone("hTofTimeVsMom_cands_pion+"); + + for (size_t iP = 1; iP < fHMean.fCandNames.size(); iP++) { + string hName = "hTofTimeVsMom_cands_" + fHMean.fCandNames[iP]; + if (iP <= 3) el->Add(fHMean.H2(hName.c_str())); + else if (iP == 4) continue; + else bg->Add(fHMean.H2(hName.c_str())); + } + + DrawSignificance(el, bg, "hTofTimeVsMom", 1e-9, 1e-2, "right"); + } + + // Hit-Track-Distance in ToF (filled after ElID step) + { + TH2D* el = fHMean.H2Clone("hTofHitTrackDist_cands_plutoEl+"); + TH2D* bg = fHMean.H2Clone("hTofHitTrackDist_cands_pion+"); + + for (size_t iP = 1; iP < fHMean.fCandNames.size(); iP++) { + string hName = "hTofHitTrackDist_cands_" + fHMean.fCandNames[iP]; + if (iP <= 1) el->Add(fHMean.H2(hName.c_str())); // consider UrQMD-electrons as "BG" + else if (iP == 4) continue; + else bg->Add(fHMean.H2(hName.c_str())); + } + + DrawSignificance(el, bg, "hTofHitTrackDist", 1e-9, 1e-2, "right"); + } +} + +void LmvmDrawAll::DrawGTrackVertex() +{ + TCanvas* c = fHMean.fHM.CreateCanvas("Vertex/globalTrackRZ", "Vertex/globalTrackRZ", 1600, 800); + c->Divide(2, 1); + c->cd(1); + TH2D* h = fHMean.H2Clone("hVertexGTrackRZ"); + h->GetYaxis()->SetRangeUser(0., 50); + h->GetZaxis()->SetRangeUser(1e-7, 10); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + c->cd(2); + TH2D* h2 = fHMean.H2Clone("hVertexGTrackRZ"); + h2->GetXaxis()->SetRangeUser(-49., -39.); + h2->GetYaxis()->SetRangeUser(0., 10.); + h2->GetZaxis()->SetRangeUser(1e-7, 10); + DrawH2(h2, kLinear, kLinear, kLog, "colz"); +} + +void LmvmDrawAll::DrawTofHitXY() +{ + //double min = 1e-7; + //double max = 2e-3; + double minD = 1e-7; + double maxD = 2e-2; + + /*TCanvas* c = fHMean.fHM.CreateCanvas("ToF/point-hit-dist_gTracks", "ToF/point-hit-dist_gTracks", 2400, 2400); // TODO: do these as loop for pile and not-pile histograms + c->Divide(4, 4); + int i = 1; + for (const string& ptcl : fHMean.fGTrackNames) { + c->cd(i++); + string hName = "hTofHitPointDist_" + ptcl; + TH1D* h = fHMean.H1Clone(hName.c_str()); + h->Fit("gaus", "Q", "", 0., 2.); + h->GetYaxis()->SetRangeUser(minD, maxD); + DrawH1(h, kLinear, kLog, "hist"); + DrawTextOnPad(fHMean.fGTrackLatex[i-2], 0.25, 0.9, 0.75, 0.999); + TF1* func = h->GetFunction("gaus"); + double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; + DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 2), 0.2, 0.65, 0.6, 0.89); + } + + for (const string& ptcl : fHMean.fGTrackNames) { + TCanvas* c2 = fHMean.fHM.CreateCanvas("ToF/HitXY/" + ptcl, "ToF/HitXY/" + ptcl, 2400, 800); + c2->Divide(3,1); + c2->cd(1); + TH2D* point = fHMean.H2Clone("hTofPointXY_" + ptcl); + point->GetZaxis()->SetRangeUser(min, max); + DrawH2(point, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("ToF Points", 0.3, 0.9, 0.7, 0.99); + c2->cd(2); + TH2D* hit = fHMean.H2Clone("hTofHitXY_" + ptcl); + hit->GetZaxis()->SetRangeUser(min, max); + DrawH2(hit, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("ToF Hits", 0.3, 0.9, 0.7, 0.99); + c2->cd(3); + TH1D* dist = fHMean.H1Clone("hTofHitPointDist_" + ptcl); + dist->GetYaxis()->SetRangeUser(minD, maxD); + dist->Fit("gaus", "Q", "", 0., 2.); + DrawH1(dist, kLinear, kLog, "hist"); + DrawTextOnPad("Distance Hit-Point", 0.3, 0.9, 0.7, 0.99); + TF1* func = dist->GetFunction("gaus"); + double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; + DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 1), 0.2, 0.65, 0.5, 0.89); + }*/ + + // draw only for particles in ToF pile that are identified as electrons in ToF + double minH = 2e-6; + double maxH = 1e-4; + + TCanvas* c3 = fHMean.fHM.CreateCanvas("ToF/HitXY/MisidsInTofPile/point-hit-dist_all", "ToF/HitXY/MisidsInTofPile/point-hit-dist_all", 2400, 1600); + c3->Divide(3, 2); + int iC = 1; + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + c3->cd(iC++); + string ptcl = fHMean.fCandNames[iP]; + string hName = "hTofPileHitPointDist_" + ptcl; + TH1D* h = fHMean.H1Clone(hName.c_str()); + h->GetYaxis()->SetRangeUser(minD, maxD); + h->Fit("gaus", "Q", "", 0., 2.); + DrawH1(h, kLinear, kLog, "hist"); + DrawTextOnPad(fHMean.fCandLatex[iP], 0.4, 0.9, 0.54, 0.99); + TF1* func = h->GetFunction("gaus"); + double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; + DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 2), 0.2, 0.65, 0.6, 0.89); + } + + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + string ptcl = fHMean.fCandNames[iP]; + TCanvas* c4 = fHMean.fHM.CreateCanvas("ToF/HitXY/MisidsInTofPile/" + ptcl, "ToF/HitXY/MisidsInTofPile/" + ptcl, 2400, 800); + c4->Divide(3,1); + c4->cd(1); + TH2D* hit = fHMean.H2Clone("hTofPileHitXY_" + ptcl); + hit->GetZaxis()->SetRangeUser(minH, maxH); + DrawH2(hit, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("ToF Hits", 0.3, 0.9, 0.7, 0.99); + c4->cd(2); + TH2D* point = fHMean.H2Clone("hTofPilePointXY_" + ptcl); + point->GetZaxis()->SetRangeUser(minH, maxH); + DrawH2(point, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("ToF Points", 0.3, 0.9, 0.7, 0.99); + c4->cd(3); + TH1D* dist = fHMean.H1Clone("hTofPileHitPointDist_" + ptcl); + dist->Fit("gaus", "Q", "", 0., 2.); + dist->GetYaxis()->SetRangeUser(minD, maxD); + DrawH1(dist, kLinear, kLog, "hist"); + DrawTextOnPad("Distance Hit-Point", 0.3, 0.9, 0.7, 0.99); + TF1* func = dist->GetFunction("gaus"); + double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; + DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 1), 0.2, 0.65, 0.5, 0.89); + } +} + +void LmvmDrawAll::DrawBetaMomSpectra() +{ + double min = 1e-5; + double max = 10; + for (const string& cat : {"gTracks", "cands"}) { + for (size_t i = 0; i < fHMean.fCandNames.size(); i++ ) { + string cName = "BetaMom/" + cat + "/" + fHMean.fCandNames[i]; + string hName = "hBetaMom_" + cat + "_" + fHMean.fCandNames[i]; + fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 800, 800); + TH2D* h = fHMean.H2Clone(hName.c_str()); + h->GetZaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad((fHMean.fCandLatex[i]).c_str(), 0.25, 0.9, 0.75, 0.999); + } + } + fHMean.fHM.CreateCanvas("BetaMom/gTracks/all", "BetaMom/gTracks/all", 800, 800); + TH2D* h = fHMean.H2Clone("hBetaMom_allGTracks"); + h->GetZaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); +} + +void LmvmDrawAll::DrawMomPluto() +{ + for (const string& pi : {"Px", "Py", "Pz"}) { + string cName = "MomDistPluto/" + pi; + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 1600); + c->Divide(3,2); + int i = 1; + for (ELmvmSignal signal : fHMean.fSignals) { + c->cd(i++); + string hName = "hMom" + pi; + TH1D* ep = H(signal)->H1Clone((hName + "+_signal_mc").c_str()); + TH1D* em = H(signal)->H1Clone((hName + "-_signal_mc").c_str()); + double bW = ep->GetBinWidth(1); + ep->Scale(1. / (H(signal)->H1("hEventNumber")->GetEntries() * bW)); + em->Scale(1. / (H(signal)->H1("hEventNumber")->GetEntries() * bW)); + ep->GetYaxis()->SetTitle("1/N dN/dP [GeV/c]^{-1}"); + DrawH1({ep, em}, {"e^{+}", "e^{-}"}, kLinear, kLog, true, 0.75, 0.75, 0.91, 0.91); + DrawTextOnPad(fHMean.fSignalNames[static_cast<int>(signal)].c_str(), 0.23, 0.8, 0.48, 0.89); + } + } +} + +void LmvmDrawAll::DrawMomentum() +{ + // draw momentum of all global tracks + fHMean.DrawAllGTracks(1, "hMom/gTracks", "hMom_gTracks", {""}, {""}, 1e-7, 10.); + + // draw momentum for misidentified and not-misidentified particles (from global tracks) + { + double min = 5e-8; + double max = 10; + TCanvas* c = fHMean.fHM.CreateCanvas("hMom/misid_vs_True", "hMom/misid_vs_True", 2400, 1600); + c->Divide(3, 2); + int iC = 1; + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + c->cd(iC++); + string ptcl = fHMean.fCandNames[iP]; + TH1D* hMis = fHMean.H1Clone("hMom_" + ptcl + "_misid"); + TH1D* hTrue = fHMean.H1Clone("hMom_" + ptcl + "_true"); + hMis->GetYaxis()->SetRangeUser(min, max); + hTrue->GetYaxis()->SetRangeUser(min, max); + DrawH1({hTrue, hMis}, {"not misidentified", "misidentified as electron"}, kLinear, kLog, true, 0.55, 0.8, 0.95, 0.91 ,"hist"); + DrawTextOnPad(fHMean.fCandLatex[iP], 0.4, 0.9, 0.54, 0.99); + } + } +} + +void LmvmDrawAll::DrawCutEffSignal() +{ + TH1D* mc = fHMean.H1Clone("hMom_cands_" + fHMean.fCandNames[0], ELmvmAnaStep::Mc); + mc->Add(fHMean.H1("hMom_cands_" + fHMean.fCandNames[1], ELmvmAnaStep::Mc)); + + TH1D* elid = fHMean.H1Clone("hMom_cands_" + fHMean.fCandNames[0], ELmvmAnaStep::ElId); + elid->Add(fHMean.H1("hMom_cands_" + fHMean.fCandNames[1], ELmvmAnaStep::ElId)); + elid->Divide(mc); + + TH1D* gamma = fHMean.H1Clone("hMom_cands_" + fHMean.fCandNames[0], ELmvmAnaStep::GammaCut); + gamma->Add(fHMean.H1("hMom_cands_" + fHMean.fCandNames[1], ELmvmAnaStep::GammaCut)); + gamma->Divide(mc); + + TH1D* tt = fHMean.H1Clone("hMom_cands_" + fHMean.fCandNames[0], ELmvmAnaStep::TtCut); + tt->Add(fHMean.H1("hMom_cands_" + fHMean.fCandNames[1], ELmvmAnaStep::TtCut)); + tt->Divide(mc); + + TH1D* pt = fHMean.H1Clone("hMom_cands_" + fHMean.fCandNames[0], ELmvmAnaStep::PtCut); + pt->Add(fHMean.H1("hMom_cands_" + fHMean.fCandNames[1], ELmvmAnaStep::PtCut)); + pt->Divide(mc); + + fHMean.fHM.CreateCanvas("CutEff_signal", "CutEff_signal", 900, 900); + elid->GetYaxis()->SetTitle("Efficiency"); + DrawH1({elid, gamma, tt, pt}, {"El-ID", "Gamma cut", "TT cut", "P_{t} cut"}, kLinear, kLinear, true, 0.6, 0.75, 0.95, 0.91, "hist"); + DrawTextOnPad("Signal Efficiency", 0.3, 0.9, 0.7, 0.99); +} + +void LmvmDrawAll::DrawPionSuppression() +{ + for (auto step : {ELmvmAnaStep::Mc, ELmvmAnaStep::Acc}) { + TH1D* elidMc = fHMean.H1Clone("hPionSupp_pi", step); + elidMc->Divide(fHMean.H1("hPionSupp_idEl", ELmvmAnaStep::ElId)); + + TH1D* gammaMc = fHMean.H1Clone("hPionSupp_pi", step); + gammaMc->Divide(fHMean.H1("hPionSupp_idEl", ELmvmAnaStep::GammaCut)); + + TH1D* ttMc = fHMean.H1Clone("hPionSupp_pi", step); + ttMc->Divide(fHMean.H1("hPionSupp_idEl", ELmvmAnaStep::TtCut)); + + TH1D* ptMc = fHMean.H1Clone("hPionSupp_pi", step); + ptMc->Divide(fHMean.H1("hPionSupp_idEl", ELmvmAnaStep::PtCut)); + + string text = fHMean.fAnaStepNames[static_cast<int>(step)]; + fHMean.fHM.CreateCanvas("PionSuppression_" + text, "PionSuppression_" + text, 900, 900); + elidMc->GetYaxis()->SetTitle("Pion Suppression"); + DrawH1({elidMc, gammaMc, ttMc, ptMc}, {"El-ID", "Gamma cut", "TT cut", "P_{t} cut"}, kLinear, kLog, true, 0.7, 0.8, 0.95, 0.91, "hist"); + DrawTextOnPad("Pion Suppression (norm. to " + text + ")", 0.25, 0.88, 0.75, 0.99); + } +} + +void LmvmDrawAll::DrawTofM2() +{ + // for all global tracks + { + TH2D* hEl = fHMean.H2Clone("hTofM2_gTracks_" + fHMean.fGTrackNames[0]); + for (size_t iP = 1; iP < 4; iP++) { + string hName = "hTofM2_gTracks_" + fHMean.fGTrackNames[iP]; + hEl->Add(fHMean.H2(hName.c_str())); + } + TH2D* hBg = fHMean.H2Clone("hTofM2_gTracks_"+ fHMean.fGTrackNames[4]); + for (size_t iP = 5; iP < fHMean.fGTrackNames.size(); iP++) { + string hName = "hTofM2_gTracks_" + fHMean.fGTrackNames[iP]; + hBg->Add(fHMean.H2(hName.c_str())); + } + + double minX = 0.; + double maxX = 2.5; + double minY = -0.05; + double maxY = 0.05; + double minZ = 1e-8; + double maxZ = 1; + vector<TLine*> lines {new TLine(0., 0.01, 1.3, 0.01), new TLine(1.3, 0.01, 2.5, 0.022)}; // set by hand + TCanvas* c = fHMean.fHM.CreateCanvas("ToF/Purity/gTracks", "ToF/Purity/gTracks", 1600, 800); // TODO: check: is this all gTracks? + c->Divide(2, 1); + c->cd(1); + hEl->GetXaxis()->SetRangeUser(minX, maxX); + hEl->GetYaxis()->SetRangeUser(minY, maxY); + hEl->GetZaxis()->SetRangeUser(minZ, maxZ); + DrawH2(hEl, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("Electrons", 0.35, 0.89, 0.65, 0.99); + for (size_t i = 0; i < lines.size(); i++) { + lines[i]->SetLineWidth(2.); + lines[i]->Draw(); + } + c->cd(2); + hBg->GetXaxis()->SetRangeUser(minX, maxX); + hBg->GetYaxis()->SetRangeUser(minY, maxY); + hBg->GetZaxis()->SetRangeUser(minZ, maxZ); + DrawH2(hBg, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("Other", 0.35, 0.89, 0.65, 0.99); + for (size_t i = 0; i < lines.size(); i++) { + lines[i]->SetLineWidth(2.); + lines[i]->Draw(); + } + } + + // for candidates + { + double minX = 0.; + double maxX = 2.5; + double minY = -0.05; + double maxY = 0.05; + double minZ = 1e-8; + double maxZ = 10; + 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 (auto step : fHMean.fAnaSteps) { + TH2D* hEl = fHMean.H2Clone(("hTofM2_cands_" + fHMean.fCandNames[0]), step); + for (size_t iP = 1; iP < 4; iP++) { + string hName = "hTofM2_cands_" + fHMean.fCandNames[iP] + "_"+ fHMean.fAnaStepNames[static_cast<int>(step)]; + hEl->Add(fHMean.H2(hName.c_str())); + } + TH2D* hBg = fHMean.H2Clone(("hTofM2_cands_" + fHMean.fCandNames[4]), step); + for (size_t iP = 5; iP < fHMean.fCandNames.size(); iP++) { + string hName = "hTofM2_cands_" + fHMean.fCandNames[iP] + "_"+ fHMean.fAnaStepNames[static_cast<int>(step)]; + hBg->Add(fHMean.H2(hName.c_str())); + } + string cName = "ToF/Purity/" + fHMean.fAnaStepNames[static_cast<int>(step)]; + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1600, 800); + c->Divide(2, 1); + c->cd(1); + hEl->GetXaxis()->SetRangeUser(minX, maxX); + hEl->GetYaxis()->SetRangeUser(minY, maxY); + hEl->GetZaxis()->SetRangeUser(minZ, maxZ); + DrawH2(hEl, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("Electrons (" + fHMean.fAnaStepNames[static_cast<int>(step)] + ")", 0.3, 0.89, 0.6, 0.99); + for (size_t i = 0; i < lines.size(); i++) { + lines[i]->SetLineWidth(2.); + lines[i]->Draw(); + } + c->cd(2); + hBg->GetXaxis()->SetRangeUser(minX, maxX); + hBg->GetYaxis()->SetRangeUser(minY, maxY); + hBg->GetZaxis()->SetRangeUser(minZ, maxZ); + DrawH2(hBg, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("Other (" + fHMean.fAnaStepNames[static_cast<int>(step)] + ")", 0.3, 0.89, 0.55, 0.99); + for (size_t i = 0; i < lines.size(); i++) { + lines[i]->SetLineWidth(2.); + lines[i]->Draw(); + } + } + } +} + +void LmvmDrawAll::DrawPtYAndTofM2Elid() +{ + double x[200], y[200]; + int n = 200; + double m = 511e-6; + double m2 = m*m; + double p = 6; + double p2 = p*p; + for (int i=0; i < n; i++) { + x[i] = 0.02 * i; + double counter = std::sqrt(2 * m2 * TMath::Exp(2 * x[i]) - m2 * TMath::Exp(4 * x[i]) + 4 * p2 * TMath::Exp(2 * x[i]) - m2); + double denom = std::sqrt(1 + 2 * TMath::Exp(2 * x[i]) + TMath::Exp(4 * x[i])); + double r = counter / denom; + y[i] = r; + } + auto pLine = new TGraph(n, x, y); // draw momentum line into Pt-Y histogram + + for (const string& cat : {"hPtY", "hTofM2"}) { + // draw all global tracks seperate for diff. particles + { + double min = 1e-8; + double max = 10; + string cName = cat + "/globalTracks/all"; + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 2400); + c->Divide(4, 4); + int i = 1; + for (auto ptcl : fHMean.fGTrackNames) { + c->cd(i++); + string hName = cat + "_gTracks_" + ptcl; + TH2D* h = fHMean.H2Clone(hName.c_str()); + h->GetZaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad(fHMean.fGTrackLatex[i-2], 0.25, 0.9, 0.75, 0.999); + if (cat == "hPtY") pLine->Draw("C"); + } + } + + // draw zoomed in TofM2 + { + if (cat == "hTofM2") { + double min = 1e-8; + double max = 10; + string cName = cat + "/globalTracks/all_zoom"; + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 2400); + c->Divide(4, 4); + int i = 1; + for (auto ptcl : fHMean.fGTrackNames) { + c->cd(i++); + string hName = cat + "_gTracks_" + ptcl; + TH2D* h = fHMean.H2Clone(hName.c_str()); + h->GetXaxis()->SetRangeUser(0., 2.5); + h->GetYaxis()->SetRangeUser(-0.05, 0.05); + h->GetZaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad(fHMean.fGTrackLatex[i-2], 0.25, 0.9, 0.75, 0.999); + } + } + } + + // draw misidentified and not-misidentified particles and ratio of both (all from global tracks) + { + double min = (cat == "hPtY") ? 1e-7 : 1e-8; + double max = (cat == "hPtY") ? 10 : 10; + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + string ptcl = fHMean.fCandNames[iP]; + string cName = cat + "/candidates/" + ptcl + "_misid"; + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 800); + c->Divide(3, 1); + c->cd(1); + TH2D* hTrue = fHMean.H2Clone(cat + "_" + ptcl + "_true"); + if (cat == "hTofM2") { + hTrue->GetXaxis()->SetRangeUser(0., 2.5); + hTrue->GetYaxis()->SetRangeUser(-0.05, 0.05); + } + hTrue->GetZaxis()->SetRangeUser(min, max); + DrawH2(hTrue, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad(fHMean.fCandLatex[iP] + " (not misidentified)" , 0.2, 0.9, 0.75, 0.99); + if (cat == "hPtY") pLine->Draw("C"); + + c->cd(2); + TH2D* hMis = fHMean.H2Clone(cat + "_" + ptcl + "_misid"); + if (cat == "hTofM2") { + hMis->GetXaxis()->SetRangeUser(0., 2.5); + hMis->GetYaxis()->SetRangeUser(-0.05, 0.05); + } + hMis->GetZaxis()->SetRangeUser(min, max); + DrawH2(hMis, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad(fHMean.fCandLatex[iP] + " (misidentified as electron)" , 0.15, 0.9, 0.75, 0.99); + if (cat == "hPtY") pLine->Draw("C"); + + c->cd(3); + TH2D* hRat = fHMean.H2Clone(cat + "_" + ptcl + "_misid"); + hRat->Divide(hTrue); + if (cat == "hTofM2") { + hRat->GetXaxis()->SetRangeUser(0., 2.5); + hRat->GetYaxis()->SetRangeUser(-0.05, 0.05); + } + hRat->GetZaxis()->SetRangeUser(1e-5, 1); // TODO: check best min value + hRat->GetZaxis()->SetTitle("Ratio #misid./#not misid."); + DrawH2(hRat, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad(fHMean.fCandLatex[iP] + ": Ratio #misid./#not misid.", 0.1, 0.9, 0.6, 0.99); + if (cat == "hPtY") pLine->Draw("C"); + } + } + + // draw misidentified candidates for each step + { + double min = (cat == "hPtY") ? 2e-8 : 2e-8; + double max = (cat == "hPtY") ? 5e-3 : 5e-3; + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { + string ptcl = fHMean.fCandNames[iP]; + string cName = cat + "/candidates/" + ptcl + "_misid_steps"; + TCanvas* cPty = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); + cPty->Divide(3, 3); + int i = 1; + for (auto step : fHMean.fAnaSteps) { + if (step < ELmvmAnaStep::ElId) continue; + cPty->cd(i++); + string hName = cat + "_cands_" + ptcl; + TH2D* h = fHMean.H2Clone(hName.c_str(), step); + if (cat == "hTofM2") { + h->GetXaxis()->SetRangeUser(0., 2.5); + h->GetYaxis()->SetRangeUser(-0.05, 0.05); + } + h->GetZaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + fHMean.DrawAnaStepOnPad(step); + if (cat == "hPtY") pLine->Draw("C"); + } + } + } + } // cat: "hPtY", "hTofM2" +} + +void LmvmDrawAll::DrawTofPilePids() +{ + // draw particles that occur in "ToF Pile" + vector<string> xLabel = {"e^{-}_{PLUTO}", "e^{+}_{PLUTO}", "e^{-}_{UrQMD_prim}", "e^{+}_{UrQMD_prim}", "e^{-}_{UrQMD_sec}", "e^{+}_{UrQMD_sec}", "#pi^{+}", "#pi^{-}", "p", "K^{+}", "K^{-}", "o."}; + double max = 3; + double min = 1e-5; + TCanvas* cTofPile = fHMean.fHM.CreateCanvas("ToF/TofPilePids", "ToF/TofPilePids", 1800, 1800); + cTofPile->Divide(4, 4); + + cTofPile->cd(1); + TH1D* h0 = fHMean.H1Clone("hTofPilePdgs_gTracks"); + for (size_t iL = 0; iL < xLabel.size(); iL++) { + h0->GetXaxis()->SetBinLabel(iL + 1, xLabel[iL].c_str()); + } + h0->GetYaxis()->SetRangeUser(min, max); + DrawH1(h0, kLinear, kLog, "hist"); + DrawTextOnPad("Global Tracks", 0.35, 0.9, 0.65, 0.999); + + int i = 2; + for (auto step : fHMean.fAnaSteps) { + cTofPile->cd(i++); + TH1D* h = fHMean.H1Clone("hTofPilePdgs_cands", step); + for (size_t iL = 0; iL < xLabel.size(); iL++) { + h->GetXaxis()->SetBinLabel(iL + 1, xLabel[iL].c_str()); + } + h->GetYaxis()->SetRangeUser(min, max); + DrawH1(h, kLinear, kLog, "hist"); + fHMean.DrawAnaStepOnPad(step); + } +} + +void LmvmDrawAll::DrawPurity() +{ + // Misidentified Particles vs. Momentum + { + vector<string> yLabel = {"e^{#pm}_{PLUTO}", "e^{#pm}_{UrQMD}", "#pi^{#pm}", "p", "K^{+}", "o."}; + for (ELmvmAnaStep step : fHMean.fAnaSteps) { + if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue; + string cName = "Purity/purity_" + fHMean.fAnaStepNames[static_cast<int>(step)]; + string hName = "hCandPdgVsMom_" + fHMean.fAnaStepNames[static_cast<int>(step)]; + fHMean.fHM.CreateCanvas(cName, cName, 800, 600); + TH2D* h = fHMean.H2Clone(hName); + h->GetZaxis()->SetRangeUser(5e-7, 1e-2); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + for (size_t y = 1; y <= yLabel.size(); y++) { + h->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + } + double nEl = h->Integral(1, h->GetXaxis()->GetNbins(), 2, 2); // do not count PLUTO particles + double purity = (nEl / h->Integral(1, h->GetXaxis()->GetNbins(), 2, h->GetYaxis()->GetNbins())) * 100; + DrawTextOnPad("Purity: " + Cbm::NumberToString(purity, 1) + " %", 0.1, 0.9, 0.45, 0.99); + } + } + + // Purity vs. Momentum after ElID step + { + int nBins = fHMean.H2("hCandPdgVsMom_elid")->GetXaxis()->GetNbins(); + double xMin = fHMean.H2("hCandPdgVsMom_elid")->GetXaxis()->GetXmin(); + double xMax = fHMean.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 = fHMean.H2("hCandPdgVsMom", ELmvmAnaStep::ElId)->GetBinContent(i, 2); + double nAll = fHMean.H2("hCandPdgVsMom", ELmvmAnaStep::ElId)->Integral(i, i, 2, fHMean.H2("hCandPdgVsMom_elid")->GetYaxis()->GetNbins()); + double val = (nAll != 0) ? 100 * nEl / nAll : 0.; + purity->SetBinContent(i, val); + } + purity->Rebin(5); + purity->Scale(1. / 5.); + fHMean.fHM.CreateCanvas("Purity/purityVsMom_elid", "Purity/purityVsMom_elid", 800, 800); + DrawH1(purity, kLinear, kLinear, "p"); + } + + // Purity seperate for ID Detectors + { + vector<string> yLabel = {"#pi^{+}", "#pi^{-}", "p", "K^{+}", "K^{-}", "o."}; + double min = 1e-7; + double max = 10; + + TH2D* rich1 = fHMean.H2Clone("hPdgVsMom_gTracks_rich_all"); + TH2D* trd1 = fHMean.H2Clone("hPdgVsMom_gTracks_trd_all"); + TH2D* tof1 = fHMean.H2Clone("hPdgVsMom_gTracks_tof_all"); + TH2D* rich2 = fHMean.H2Clone("hPdgVsMom_gTracks_rich_complete"); + TH2D* trd2 = fHMean.H2Clone("hPdgVsMom_gTracks_trd_complete"); + TH2D* tof2 = fHMean.H2Clone("hPdgVsMom_gTracks_tof_complete"); + + rich1->GetZaxis()->SetRangeUser(min, max); + trd1 ->GetZaxis()->SetRangeUser(min, max); + tof1 ->GetZaxis()->SetRangeUser(min, max); + rich2->GetZaxis()->SetRangeUser(min, max); + trd2 ->GetZaxis()->SetRangeUser(min, max); + tof2 ->GetZaxis()->SetRangeUser(min, max); + + for (size_t y = 1; y <= yLabel.size(); y++) { + rich1->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + trd1 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + tof1 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + rich2->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + trd2 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + tof2 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + } + + TCanvas* c1 = fHMean.fHM.CreateCanvas("Misid/misid_gTracks_all", "Misid/misid_gTracks_all", 2400, 800); + c1->Divide(3,1); + c1->cd(1); + DrawH2(rich1, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("RICH", 0.25, 0.9, 0.75, 0.999); + DrawPurityHistText(rich1); + c1->cd(2); + DrawH2(trd1, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("TRD", 0.25, 0.9, 0.75, 0.999); + DrawPurityHistText(trd1); + c1->cd(3); + DrawH2(tof1, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("ToF", 0.25, 0.9, 0.75, 0.999); + DrawPurityHistText(tof1); + + TCanvas* c2 = fHMean.fHM.CreateCanvas("Misid/misid_gTracks_complete", "Misid/misid_gTracks_complete", 2400, 800); + c2->Divide(3,1); + c2->cd(1); + DrawH2(rich2, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("RICH", 0.25, 0.9, 0.75, 0.999); + DrawPurityHistText(rich2); + + c2->cd(2); + DrawH2(trd2, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("TRD", 0.25, 0.9, 0.75, 0.999); + DrawPurityHistText(trd2); + + c2->cd(3); + DrawH2(tof2, kLinear, kLinear, kLog, "colz"); + DrawTextOnPad("ToF", 0.25, 0.9, 0.75, 0.999); + DrawPurityHistText(tof2); + } +} + +void LmvmDrawAll::DrawPurityHistText(TH2D* h) +{ + int nX = h->GetXaxis()->GetNbins(); + + double xMin = 0.4; + double xMax = 0.75; + double piP = h->Integral(1, nX, 1, 1); + double piM = h->Integral(1, nX, 2, 2); + double p = h->Integral(1, nX, 3, 3); + double KP = h->Integral(1, nX, 4, 4); + double KM = h->Integral(1, nX, 5, 5); + double o = h->Integral(1, nX, 6, 6); + + DrawTextOnPad("#pi^{+}: " + Cbm::NumberToString(piP, 1) + " /Ev", xMin, 0.12, xMax, 0.3); + DrawTextOnPad("#pi^{-}: " + Cbm::NumberToString(piM, 1) + " /Ev", xMin, 0.28, xMax, 0.37); + DrawTextOnPad("p: " + Cbm::NumberToString(p, 1) + " /Ev", xMin, 0.43, xMax, 0.48); + DrawTextOnPad("K^{+}: " + Cbm::NumberToString(KP, 1) + " /Ev", xMin, 0.53, xMax, 0.65); + DrawTextOnPad("K^{-}: " + Cbm::NumberToString(KM, 1) + " /Ev", xMin, 0.66, xMax, 0.76); + DrawTextOnPad("o.: " + Cbm::NumberToString(o, 1) + " /Ev", xMin, 0.77, xMax, 0.9); +} + +void LmvmDrawAll::InvestigateMisid() +{ + DrawMomentum(); + DrawPtYAndTofM2Elid(); + DrawTofPilePids(); + DrawChi2(); + + // Connection between Misidentification and Mismatches + { + vector<string> xLabel = {"0", "1", "2", "3"}; + vector<string> yLabel = {"true-ID", "mis-ID"}; + fHMean.DrawAllGTracks(2, "hMatchId_gTracks", "hMatchId_gTracks", {xLabel}, {yLabel}, 1e-7, 1e3); + } + + // RICH: Ring-Track Distance + fHMean.DrawAllGTracks(2, "hRichRingTrackDist/gTracks_all", "hRichRingTrackDist_gTracks", {""}, {""}, 1e-7, 1e-2); // TODO: move this anywhere else? + fHMean.DrawAllCandsAndSteps(2, "hRichRingTrackDist/", "hRichRingTrackDist_cands", {""}, {""}, 1e-7, 1e-2); + + // ToF + { + double minXY = 1e-6; + double maxXY = 1e-3; + double minD = 3e-6; + double maxD = 3e-1; + for (const string& cat : {"trueid", "misid", "truematch", "mismatch"}) { + fHMean.DrawAllGTracks(2, "ToF/HitTrackDist/HitPointDist_" + cat, "hTofHitTrackDist_gTracks_" + cat, {""}, {""}, 1e-9, 1.); + fHMean.DrawAllGTracks(2, "ToF/HitXY/PointXY_" + cat, "hTofPointXY_" + cat, {""}, {""}, minXY, maxXY); + fHMean.DrawAllGTracks(2, "ToF/HitXY/HitXY_" + cat, "hTofHitXY_" + cat, {""}, {""}, minXY, maxXY); + fHMean.DrawAllGTracks(1, "ToF/HitXY/HitPointDist_" + cat, "hTofHitPointDist_" + cat, {""}, {""}, minD, maxD); + } + fHMean.DrawAllCands(2, "ToF/Time/HitTrackDist_cands", "hTofHitTrackDist_cands", {""}, {""}, 1e-9, 1.); + fHMean.DrawAllCands(2, "ToF/Time/TimeVsMom_", "hTofTimeVsMom_cands", {""}, {""}, 1e-9, 1e-2); + fHMean.DrawAllGTracks(2, "ToF/Time/IdAndMatch", "hTofM2Calc_gTracks", {""}, {"true-ID", "mis-ID", "truematch", "mismatch"}, 1e-9, 10); + + TCanvas* c1 = fHMean.fHM.CreateCanvas("ToF/HitXY/HitPointDist_all", "ToF/HitXY/HitPointDist_all", 2400, 2400); + c1->Divide(4, 4); + int iC = 1; + for (size_t iP = 0; iP < fHMean.fGTrackNames.size(); iP++) { + c1->cd(iC++); + string hName = (iP <= 3) ? "hTofHitPointDist_trueid_" + fHMean.fGTrackNames[iP] : "hTofHitPointDist_misid_" + fHMean.fGTrackNames[iP]; + string text = (iP <= 3) ? "true-ID" : "mis-ID"; + TH1D* h = fHMean.H1Clone(hName); + h->GetYaxis()->SetRangeUser(minD, maxD); + h->Fit("gaus", "Q", "", 0., 2.); + DrawH1(h, kLinear, kLog, "hist"); + DrawTextOnPad(text, 0.5, 0.7, 0.8, 0.9); + DrawTextOnPad(fHMean.fGTrackLatex[iP], 0.25, 0.9, 0.75, 0.99); + TF1* func = h->GetFunction("gaus"); + double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; + DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 2), 0.2, 0.65, 0.6, 0.89); + } + } + + // draw chi2 of true-matched vs. mismatched + for (const string& match : {"truematch", "mismatch"}) { + for (const string& cat : {"all", "complete"}) { + for (const string& det : {"rich", "trd", "tof"}) { + string cName = "hChi2/" + match + "_" + det + "_" + cat; + string hName = "hChi2_" + match + "_" + cat + "_" + det; + fHMean.DrawAllGTracks(1, cName, hName, {""}, {""}, 2e-7, 10); + } + } + } + + // draw vertex of ToF-misidentifications + { + double min = 1e-7; + for (auto step : fHMean.fAnaSteps) { + string cName = "ToF/Vertex/tofMisid_" + fHMean.fAnaStepNames[static_cast<int>(step)]; + vector<TH2D*> xyz {fHMean.H2("hVertexXZ_misidTof", step), fHMean.H2("hVertexYZ_misidTof", step), fHMean.H2("hVertexXY_misidTof", step)}; + + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 600); + c->Divide(3, 1); + for (size_t i = 0; i < xyz.size(); i++) { + c->cd(i + 1); + DrawH2(xyz[i]); + xyz[i]->SetMinimum(min); + gPad->SetLogz(true); + } + + TCanvas* cZoom = fHMean.fHM.CreateCanvas((cName + "_zoom").c_str(), (cName + "_zoom").c_str(), 1800, 600); + cZoom->Divide(3, 1); + for (size_t i = 0; i < xyz.size(); i++) { + TH2D* hZoom = (TH2D*) xyz[i]->Clone(); + cZoom->cd(i + 1); + if (i == 2) { + hZoom->GetXaxis()->SetRangeUser(-10., 10.); + hZoom->GetYaxis()->SetRangeUser(-10., 10.); + } + else { + hZoom->GetXaxis()->SetRangeUser(fZ - 1., fZ + 10.); + hZoom->GetYaxis()->SetRangeUser(-10., 10.); + } + DrawH2(hZoom); + hZoom->SetMinimum(min); + gPad->SetLogz(true); + } + + fHMean.fHM.CreateCanvas((cName + "_RZ").c_str(), (cName + "_RZ").c_str(), 900, 900); + TH2D* hRZ = fHMean.H2Clone("hVertexRZ_misidTof", step); + hRZ->SetMinimum(min); + DrawH2(hRZ, kLinear, kLinear, kLog); + } + } + + // draw RICH/TRD/ToF index vs. STS index + /*{ + for (const string& det : {"Rich", "Trd", "Tof"}) { + string cName = "hIndex/" + det; + string hName = "hIndexSts" + det; + fHMean.DrawAllGTracks(2, cName, hName, {""}, {""}, 1e-7, 2e-3); + } + }*/ + + double minMism = 1e-4; + double maxMism = 100.; + + vector<string> xLabelDet = {"STS_{all}", "RICH_{all}", "RICH_{mis}", "TRD_{all}", "TRD_{mis}", "ToF_{all}", "ToF_{mis}"}; + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatches_all", "hNofMismatches_gTracks_all", xLabelDet, {""}, minMism, maxMism); + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatches_complete", "hNofMismatches_gTracks_complete", xLabelDet, {""}, minMism, maxMism); + + vector<string> xLabelSeg = {"0", "1", "2", "3"}; + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatchedTrackSegments_all", "hNofMismatchedTrackSegments_all", xLabelSeg, {""}, minMism, maxMism); + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatchedTrackSegments_complete", "hNofMismatchedTrackSegments_complete", xLabelSeg, {""}, minMism, maxMism); +} + +void LmvmDrawAll::DrawChi2() +{ + double min = 1e-9; + double max = 5e-3; + for (const string det : {"sts", "rich", "trd", "tof"}) { + fHMean.DrawAllCands(2, "hChi2/Detectors/hChi2VsMom_" + det, "hChi2VsMom_" + det, {""}, {""}, min, max); + fHMean.DrawAllCands(2, "hChi2/Detectors/hTofTimeVsChi2_" + det, "hTofTimeVsChi2_" + det, {""}, {""}, min, max); + } + + for (const string det : {"StsRich", "StsTrd", "RichTrd"}) { + fHMean.DrawAllCands(2, "hChi2/Detectors/hChi2Comb" + det, "hChi2Comb_" + det, {""}, {""}, min, max); + } +} + void LmvmDrawAll::DrawMinvAll() { for (const ELmvmAnaStep step : fHMean.fAnaSteps) { - 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); + if (step < ELmvmAnaStep::ElId) continue; + string cName = "Minv/" + fHMean.fAnaStepNames[static_cast<int>(step)]; + fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1000, 1000); DrawMinv(step); + fHMean.DrawAnaStepOnPad(step); } } @@ -180,13 +1176,10 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) TH1D* bg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step); TH1D* pi0 = fHMean.H1Clone("hMinv", ELmvmSrc::Pi0, step); TH1D* eta = fHMean.H1Clone("hMinv", ELmvmSrc::Eta, step); - TH1D* cocktail = GetCocktailMinvH1(step); + TH1D* cocktail = GetCocktailMinvH1("hMinv", step); TH1D* sbg = static_cast<TH1D*>(bg->Clone()); sbg->Add(cocktail); TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); - TH1D* cbSig = fHMean.H1Clone("hMinvCombSignalMc", step); - TH1D* cbU = fHMean.H1Clone("hMinvCombBg_urqmd", step); - TH1D* cbSigU = fHMean.H1Clone("hMinvCombSignalMc_urqmd", step); double binWidth = bg->GetBinWidth(1); vector<TH1D*> sHists(fHMean.fNofSignals); @@ -201,7 +1194,7 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) vector<LmvmDrawMinvData> drawData; if (step != ELmvmAnaStep::Mc) { drawData.emplace_back(sbg, kBlack, kBlack, 1, -1, "Cocktail+BG"); - drawData.emplace_back(bg, kGray, kBlack, 1, -1, "Background"); + drawData.emplace_back(bg, kGray, kBlack, 1, 3001, "Background"); } drawData.emplace_back(pi0, kGreen - 3, kGreen + 3, 2, -1, "#pi^{0} #rightarrow #gammae^{+}e^{-}"); @@ -217,11 +1210,7 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) "in-medium #rho"); drawData.emplace_back(cocktail, -1, kRed + 2, 4, -1, "Cocktail"); drawData.emplace_back(cb, -1, kCyan + 2, 4, -1, "CB"); - drawData.emplace_back(cbSig, -1, kCyan, 3, -1, "CB-Signal"); - drawData.emplace_back(cbU, -1, kMagenta + 3, 4, -1, "CB (UrQMD-El)"); - drawData.emplace_back(cbSigU, -1, kMagenta, 3, -1, "CB-Signal (UrQMD-El)"); - - + double min = std::numeric_limits<Double_t>::max(); double max = std::numeric_limits<Double_t>::min(); TH1D* h0 = nullptr; @@ -235,11 +1224,8 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) leg->AddEntry(d.fH, d.fLegend.c_str(), "f"); DrawH1(d.fH, kLinear, kLinear, (h0 == nullptr) ? "hist" : "hist,same", d.fLineColor, d.fLineWidth, 0); if (h0 == nullptr) h0 = d.fH; - //min = std::min(d.fH->GetMinimum(), min); - min = 1e-8; + min = 1e-9; max = std::max(d.fH->GetMaximum(), max); - //h0->SetMinimum(1e-7); - //h0->SetMaximum(50); } if (min == 0.) min = std::min(1e-4, max * 1e-7); if (h0 != nullptr) h0->SetMinimum(min); @@ -248,7 +1234,233 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) leg->SetFillColor(kWhite); leg->Draw(); gPad->SetLogy(true); - //DrawH1(hists, legendStr, kLinear, kLog, true, 0.7, 0.6, 0.99, 0.99, "HIST L"); +} + +void LmvmDrawAll::DrawMomRecoPrecision() +{ + // draw ratio P_rec/P_MC + { + for (const string& cat : {"", "_zoom", "_zoom2"}) { + string cName = "hMom/precisionReco/mom_ratioRecOverMc" + cat; + 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; + vector<TH1*> hists; + vector<string> legend; + c->cd(i++); + for (size_t iP = 0; iP < fHMean.fCandNames.size(); iP++ ) { + TH1D* h = fHMean.H1Clone("hMomRatio_cands_" + fHMean.fCandNames[iP], step); + if (cat == "_zoom") h->GetXaxis()->SetRangeUser(0.9, 1.1); + if (cat == "_zoom2") h->GetXaxis()->SetRangeUser(7., 7.5); + hists.push_back(h); + legend.push_back((fHMean.fCandLatex[iP]).c_str()); + } + DrawH1(hists, legend, kLinear, kLog, true, 0.8, 0.65, 0.99, 0.99, "hist"); + fHMean.DrawAnaStepOnPad(step); + } + } + } + + // draw ratio P_rec/P_MC vs P + { + ELmvmAnaStep step = ELmvmAnaStep::ElId; + double min = 0.5; + double max = 1.5; + TCanvas* c = fHMean.fHM.CreateCanvas("hMom/precisionReco/mom_ratioRecOverMcVsMom", "hMom/precisionReco/mom_ratioRecOverMcVsMom", 2700, 1800); + c->Divide(4,3); + int i = 1; + for (auto ptcl : fHMean.fCandNames) { + int iL = i-1; + c->cd(i++); + string hName = "hMomRatioVsMom_cands_" + ptcl; + TH2D* h = fHMean.H2Clone(hName.c_str(), step); + h->GetYaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + string text = fHMean.fCandLatex[iL] + ", "+ fHMean.fAnaStepLatex[static_cast<int>(step)]; + DrawTextOnPad(text.c_str(), 0.25, 0.9, 0.75, 0.999); + } + } +} + +void LmvmDrawAll::DrawMinvOfficialStyle() +{ + // variable binning + TH1D* h = fHMean.H1("hMinv", ELmvmSrc::Bg, ELmvmAnaStep::Mc); // template for values + double ChangeValue = 1.2; + int f = 2; // factor larger/smaller binsize + int binChange = h->GetXaxis()->FindBin(ChangeValue); + double bW = h->GetXaxis()->GetBinWidth(1); + const int nBins = (int) binChange + (h->GetNbinsX() - binChange)/f; + vector<double> binEdges(nBins+1, 0.); + + for (int iB = 1; iB <= nBins; iB++) { + binEdges[iB] = (iB < binChange) ? binEdges[iB-1] + bW : binEdges[iB-1] + f * bW; + } + binEdges[nBins+1] = h->GetBinCenter(h->GetNbinsX()) + bW/f; + + int nofInmed = H(ELmvmSignal::Inmed) ->H1("hEventNumber")->GetEntries(); + int nofQgp = H(ELmvmSignal::Qgp) ->H1("hEventNumber")->GetEntries(); + int nofOmega = H(ELmvmSignal::Omega) ->H1("hEventNumber")->GetEntries(); + int nofOmegaD = H(ELmvmSignal::OmegaD)->H1("hEventNumber")->GetEntries(); + int nofPhi = H(ELmvmSignal::Phi) ->H1("hEventNumber")->GetEntries(); + + for (auto step : fHMean.fAnaSteps) { + if (step < ELmvmAnaStep::ElId) continue; + TH1D* h_npm = new TH1D("h_npm", "h_npm", nBins, binEdges.data()); + TH1D* h_bc = new TH1D("h_bc", "h_bc", nBins, binEdges.data()); + TH1D* h_coc = new TH1D("h_coc", "h_coc", nBins, binEdges.data()); + TH1D* h_eta = new TH1D("h_eta", "h_eta", nBins, binEdges.data()); + TH1D* h_pi0 = new TH1D("h_pi0", "h_pi0", nBins, binEdges.data()); + TH1D* h_inmed = new TH1D("h_inmed", "h_inmed", nBins, binEdges.data()); + TH1D* h_qgp = new TH1D("h_qgp", "h_qgp", nBins, binEdges.data()); + TH1D* h_omega = new TH1D("h_omega", "h_omega", nBins, binEdges.data()); + TH1D* h_omegaD = new TH1D("h_omegaD", "h_omegaD", nBins, binEdges.data()); + TH1D* h_phi = new TH1D("h_phi", "h_phi", nBins, binEdges.data()); + + TH1D* hInmed0 = H(ELmvmSignal::Inmed)->H1Clone("hMinv", ELmvmSrc::Signal, step); + TH1D* hQgp0 = H(ELmvmSignal::Qgp)->H1Clone("hMinv", ELmvmSrc::Signal, step); + TH1D* hOmega0 = H(ELmvmSignal::Omega)->H1Clone("hMinv", ELmvmSrc::Signal, step); + TH1D* hOmegaD0 = H(ELmvmSignal::OmegaD)->H1Clone("hMinv", ELmvmSrc::Signal, step); + TH1D* hPhi0 = H(ELmvmSignal::Phi)->H1Clone("hMinv", ELmvmSrc::Signal, step); + + hInmed0->Rebin(fRebMinv); + hQgp0->Rebin(fRebMinv); + hOmega0->Rebin(fRebMinv); + hOmegaD0->Rebin(fRebMinv); + hPhi0->Rebin(fRebMinv); + + hInmed0->Scale(1. / (nofInmed * bW)); + hQgp0->Scale(1. / (nofQgp * bW)); + hOmega0->Scale(1. / (nofOmega * bW)); + hOmegaD0->Scale(1. / (nofOmegaD * bW)); + hPhi0->Scale(1. / (nofPhi * bW)); + + for (int iB = 1; iB <= h->GetNbinsX(); iB++) { + if (iB < binChange) { + h_npm ->SetBinContent(iB, fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(iB)); + h_bc ->SetBinContent(iB, fHMean.H1("hMinvCombBg", step)->GetBinContent(iB)); + h_coc ->SetBinContent(iB, GetCocktailMinvH1("hMinv", step)->GetBinContent(iB)); + h_eta ->SetBinContent(iB, fHMean.H1("hMinv", ELmvmSrc::Eta, step)->GetBinContent(iB)); + h_pi0 ->SetBinContent(iB, fHMean.H1("hMinv", ELmvmSrc::Pi0, step)->GetBinContent(iB)); + h_inmed ->SetBinContent(iB, hInmed0 ->GetBinContent(iB)); + h_qgp ->SetBinContent(iB, hQgp0 ->GetBinContent(iB)); + h_omega ->SetBinContent(iB, hOmega0 ->GetBinContent(iB)); + h_omegaD->SetBinContent(iB, hOmegaD0->GetBinContent(iB)); + h_phi ->SetBinContent(iB, hPhi0 ->GetBinContent(iB)); + } + else { + int iB2 = (int) binChange + (iB-binChange)/f; + h_npm ->SetBinContent(iB2, h_npm->GetBinContent(iB2) + fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(iB)/f); + h_bc ->SetBinContent(iB2, h_bc->GetBinContent(iB2) + fHMean.H1("hMinvCombBg", step)->GetBinContent(iB)/f); + h_coc ->SetBinContent(iB2, h_coc->GetBinContent(iB2) + GetCocktailMinvH1("hMinv", step)->GetBinContent(iB)/f); + h_eta ->SetBinContent(iB2, h_eta->GetBinContent(iB2) + fHMean.H1("hMinv", ELmvmSrc::Eta, step)->GetBinContent(iB)/f); + h_pi0 ->SetBinContent(iB2, h_pi0->GetBinContent(iB2) + fHMean.H1("hMinv", ELmvmSrc::Pi0, step)->GetBinContent(iB)/f); + h_inmed ->SetBinContent(iB2, h_inmed->GetBinContent(iB2) + hInmed0 ->GetBinContent(iB)/f); + h_qgp ->SetBinContent(iB2, h_qgp->GetBinContent(iB2) + hQgp0 ->GetBinContent(iB)/f); + h_omega ->SetBinContent(iB2, h_omega->GetBinContent(iB2) + hOmega0 ->GetBinContent(iB)/f); + h_omegaD->SetBinContent(iB2, h_omegaD->GetBinContent(iB2) + hOmegaD0->GetBinContent(iB)/f); + h_phi ->SetBinContent(iB2, h_phi->GetBinContent(iB2) + hPhi0 ->GetBinContent(iB)/f); + } + } + + // set errors + int nEv = GetNofTotalEvents(); + for (int iB = 1; iB <= h_npm->GetNbinsX(); iB++) { + double bW2 = h_npm->GetBinLowEdge(iB+1) - h_npm->GetBinLowEdge(iB); + //cout << "iB = " << iB << " | bW2 = " << bW2 << endl; + h_npm->SetBinError(iB, std::sqrt(h_npm->GetBinContent(iB) * bW2 * nEv) / (bW2 * nEv)); + h_bc ->SetBinError(iB, std::sqrt(h_bc ->GetBinContent(iB) * bW2 * nEv) / (bW2 * nEv)); + } + + TH1D* h_sum = (TH1D*) h_eta->Clone(); + h_sum->Add(h_pi0); + h_sum->Add(h_inmed); + h_sum->Add(h_qgp); + h_sum->Add(h_omega); + h_sum->Add(h_omegaD); + h_sum->Add(h_phi); + + TH1D* h_sum2 = (TH1D*) h_sum->Clone(); + + fHMean.SetOptH1(h_npm, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 21, 1.3, kPink, "marker"); + fHMean.SetOptH1(h_bc, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 22, 1.3, kBlue-2, "marker"); + fHMean.SetOptH1(h_sum, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 20, 1.1, kBlack, "marker"); + + fHMean.SetOptH1(h_sum2, "", "",510, 1, 1, kBlack, "line"); + fHMean.SetOptH1(h_eta, "", "",510, 1, 1, kBlue+2, "line"); + fHMean.SetOptH1(h_pi0, "", "",510, 1, 1, kAzure-4, "line"); + fHMean.SetOptH1(h_inmed, "", "",510, 1, 1, kPink, "line"); + fHMean.SetOptH1(h_qgp, "", "",510, 1, 1, kOrange, "line"); + fHMean.SetOptH1(h_omega, "", "",510, 1, 1, kGreen+1, "line"); + fHMean.SetOptH1(h_omegaD, "", "",510, 1, 1, kGreen+3, "line"); + fHMean.SetOptH1(h_phi, "", "",510, 1, 1, kMagenta+2, "line"); + + //h_npm->GetXaxis()->SetRangeUser(0.,2.); + h_npm->GetYaxis()->SetRangeUser(1e-9,10); + + string cName = "MinvOfficialStyle/minv_" + fHMean.fAnaStepNames[static_cast<int>(step)]; + TCanvas* can = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 900, 918); + fHMean.SetOptCanvas(can); + can->SetLogy(); + can->cd(); + + h_npm->Draw("peist"); + h_bc ->Draw("peistsame"); + h_sum->Draw("phistsame"); + + h_sum2 ->Draw("histsame"); + h_eta ->Draw("histsame"); + h_pi0 ->Draw("histsame"); + h_inmed ->Draw("histsame"); + h_qgp ->Draw("histsame"); + h_omega ->Draw("histsame"); + h_omegaD->Draw("histsame"); + h_phi ->Draw("histsame"); + + vector<LmvmLegend> legend3; + legend3.emplace_back(h_npm, "N^{#pm}_{same}", "p"); + legend3.emplace_back(h_bc, "CB_{calc}", "p"); + legend3.emplace_back(h_sum, "MC Cocktail", "p"); + + vector<LmvmLegend> legend8; + legend8.emplace_back(h_sum2, "MC Cocktail", "l"); + legend8.emplace_back(h_eta, "#eta#rightarrow#gammae^{+}e^{-}", "l"); + legend8.emplace_back(h_pi0, "#pi^{0}#rightarrow#gammae^{+}e^{-}", "l"); + legend8.emplace_back(h_omegaD, "#omega#rightarrow#pi^{0}e^{+}e^{-}", "l"); + legend8.emplace_back(h_omega, "#omega#rightarrowe^{+}e^{-}", "l"); + legend8.emplace_back(h_phi, "#phi#rightarrowe^{+}e^{-}", "l"); + legend8.emplace_back(h_inmed, "Rapp in-medium SF", "l"); + legend8.emplace_back(h_qgp, "Rapp QGP", "l"); + + fHMean.SetLegend(legend3, 0.035, 0.37, 0.64, 0.57, 0.77); + fHMean.SetLegend(legend8, 0.025, 0.66, 0.57, 0.82, 0.87); + + TLatex* tex = new TLatex(0.15, 2.5,"CBM Simulations"); + tex->SetTextColor(1); //kAzure+6); + tex->SetTextSize(0.035); + tex->SetTextFont(42); + tex->Draw(); + TLatex* tex2 = new TLatex(0.15, 0.7,"0 fm Au+Au, #sqrt{s_{NN}}=4.1 GeV"); + tex2->SetTextColor(1); //kAzure+6); + tex2->SetTextSize(0.035); + tex2->SetTextFont(42); + tex2->Draw(); + + /*delete h_npm; + delete h_bc; + delete h_coc; + delete h_eta; + delete h_pi0; + delete h_inmed; + delete h_qgp; + delete h_omega; + delete h_omegaD; + delete h_phi; + delete h_sum; + delete h_sum2;*/ + } // steps } void LmvmDrawAll::DrawMinvBgSourcesAll() @@ -326,8 +1538,7 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() 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->SetBinContent(i, 2); else { double r = bgRat->GetBinContent(i) / fHMean.H1("hMinv", ELmvmSrc::Bg, ELmvmAnaStep::ElId)->GetBinContent(i); bgRat->SetBinContent(i, r); @@ -335,7 +1546,7 @@ void LmvmDrawAll::DrawMinvBgSourcesAll() } bgRat->SetMinimum(0.9); - bgRat->SetMaximum(1.1); + bgRat->SetMaximum(2.1); fHMean.fHM.CreateCanvas((cName + "_compWithBg").c_str(), (cName + "_compWithBg").c_str(), 800, 800); DrawH1(bgRat, kLinear, kLinear, "hist"); @@ -346,298 +1557,229 @@ void LmvmDrawAll::DrawMinvPtAll() for (const ELmvmAnaStep step : fHMean.fAnaSteps) { string name = fHMean.GetName("minvPt/lmvmAll_minvPt", step); fHMean.fHM.CreateCanvas(name.c_str(), name.c_str(), 1000, 1000); - DrawH2(GetCocktailMinv<TH2D>("hMinvPt", step)); + DrawH2(GetCocktailMinv<TH2D>("hMinvPt", step), kLinear, kLinear, kLog, "colz"); } } void LmvmDrawAll::DrawMinvCombBgAndSignal() { - // double yMin = 1e-7; - // double yMax = 5e-2; - // string yTitle = "dN/dM_{ee}/N [GeV/c^{2}]^{-1}"; - // string xTitle = "M_{ee} [GeV/c^2]"; + string folder = "CombinatorialBackground/"; + string folderUAll = folder + "UrQMD/All/"; + string folderUEl = folder + "UrQMD/Electrons/"; - // draw MM/PP ratio for same events + // Draw PM, MM, PP in one plot { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/RatioMMPP_same", "cb/RatioMMPP_same", 1800, 1800); - c->Divide(3, 3); - int i = 1; - for (auto step : fHMean.fAnaSteps) { - if (step < ELmvmAnaStep::ElId) continue; - c->cd(i++); - TH1D* pp = fHMean.H1Clone("hMinvCombPP_sameEv", step); - TH1D* mm = fHMean.H1Clone("hMinvCombMM_sameEv", step); - mm->Divide(pp); - mm->GetYaxis()->SetTitle("Ratio e^{-}e^{-}/e^{+}e^{+}"); - DrawH1(mm, kLinear, kLinear, "hist"); - fHMean.DrawAnaStepOnPad(step); + for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { + for (const string& ev : {"sameEv", "mixedEv"}) { + string cName = (cat == "_urqmdAll") ? folderUAll + "PairYields_" + ev :(cat == "_urqmdEl") ? folderUEl + "PairYields_" + ev : folder + "PairYields_" + ev; + 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* pp = fHMean.H1Clone("hMinvCombPP" + cat + "_"+ ev, step); + TH1D* mm = fHMean.H1Clone("hMinvCombMM" + cat + "_" + ev, step); + TH1D* pm = fHMean.H1Clone("hMinvCombPM" + cat + "_" + ev, step); + pm->GetYaxis()->SetTitle("Yield"); + if (ev == "sameEv") pm->GetYaxis()->SetRangeUser(1e-6, 2e-2); + 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); + } + } } } - // Draw PM, MM, PP in one plot + // draw MM/PP ratio for same events { - for (const string& ev : {"sameEv", "mixedEv"}) { - string cName = "cb/pairYield_" + ev; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); + for (const string& cat : {"", "_urqmdEl"}) { + string cName = (cat == "_urqmdEl") ? folderUEl : folder; + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "RatioMMPP_same", cName + "RatioMMPP_same", 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); - TH1D* pp = fHMean.H1Clone("hMinvCombPP_" + ev, step); - TH1D* mm = fHMean.H1Clone("hMinvCombMM_" + ev, step); - TH1D* pm = fHMean.H1Clone("hMinvCombPM_" + ev, step); - 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"); + TH1D* pp = fHMean.H1Clone("hMinvCombPP" + cat + "_sameEv", step); + TH1D* mm = fHMean.H1Clone("hMinvCombMM" + cat + "_sameEv", step); + mm->Divide(pp); + mm->GetYaxis()->SetTitle("Ratio e^{-}e^{-}/e^{+}e^{+}"); + DrawH1(mm, kLinear, kLinear, "hist"); fHMean.DrawAnaStepOnPad(step); } } } - // Draw same and mixed for PP, MM, PM cases; normalize to 300 - 1000 MeV/c2 interval + //Draw Geometric Mean for same and mixed (mixed normalized) { - for (const string& comb : {"PM", "PP", "MM"}) { - string cName = "cb/pairYield_" + comb + "_sameMixed"; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); + for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { + string cName = (cat == "_urqmdAll") ? folderUAll : (cat == "_urqmdEl") ? folderUEl : folder; + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "geomMean", cName + "geomMean", 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* mixed = fHMean.H1Clone("hMinvComb" + comb + "_mixedEv", step); - TH1D* sameU = fHMean.H1Clone("hMinvComb" + comb + "_urqmd_sameEv", step); - TH1D* mixedU = fHMean.H1Clone("hMinvComb" + comb + "_urqmd_mixedEv", step); - int minBin = same->FindBin(0.3); - int maxBin = same->FindBin(1.0); + TH1D* same = fHMean.H1Clone("hMinvCombGeom" + cat + "_sameEv", step); + TH1D* mixed = fHMean.H1Clone("hMinvCombGeom" + cat + "_mixedEv", step); + int minBin = same->FindBin(0.4); + int maxBin = same->FindBin(0.7); double scale = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin); - double scaleU = sameU->Integral(minBin, maxBin) / mixedU->Integral(minBin, maxBin); mixed->Scale(scale); - mixedU->Scale(scaleU); + same->GetXaxis()->SetTitle("M_{ee} [GeV/c^{2}]"); 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"); + mixed->SetMinimum(1e-8); + DrawH1({same, mixed}, {"geom. mean (same)", "geom. mean (mixed)"}, kLinear, kLog, true, 0.52, 0.8, 0.99, 0.99, "p"); fHMean.DrawAnaStepOnPad(step); } } } - // Draw ratio mixed/same for PP, MM, PM cases - { - for (const string& comb : {"PM", "PP", "MM"}) { - string cName = "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); // 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"); - fHMean.DrawAnaStepOnPad(step); - } - } + // Draw k factor + { // draw all three k factors in one histo // TODO: keep this or next method or both? + ELmvmAnaStep step = ELmvmAnaStep::ElId; + fHMean.fHM.CreateCanvas(folder + "k_all", folder + "k_all", 1800, 1800); + TH1D* kAll = fHMean.H1Clone("hMinvCombK", step); + TH1D* kUAll = fHMean.H1Clone("hMinvCombK_urqmdAll", step); + TH1D* kUEl = fHMean.H1Clone("hMinvCombK_urqmdEl", step); + DrawH1({kAll, kUAll, kUEl}, {"k (all)", "k (UrQMD all)", "k (UrQMD El"}, kLinear, kLinear, true, 0.65, 0.8, 0.99, 0.95, "p"); + kAll->GetYaxis()->SetTitle("k Factor"); + kAll->GetYaxis()->SetRangeUser(0.8, 1.2); + fHMean.DrawAnaStepOnPad(step); } - // compare PM with sbg (Cocktail + BG) + // Compare N+-, BG+Coc, BG and CB after ElID step { - for (auto step : fHMean.fAnaSteps) { - if (step < ELmvmAnaStep::ElId) continue; - TH1D* pmSame = fHMean.H1Clone("hMinvCombPM_sameEv", step); - TH1D* sbg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step); - sbg->Add(GetCocktailMinvH1(step)); - TH1D* ratio = (TH1D*) pmSame->Clone(); - string cName = "cb/N+-VsCocBg_" + fHMean.fAnaStepNames[static_cast<int>(step)]; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 900); + ELmvmAnaStep step = ELmvmAnaStep::ElId; + + for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { + string cName = (cat == "_urqmdAll") ? folderUAll : (cat == "_urqmdEl") ? folderUEl : folder; + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "N+-VsBgCoc", cName + "N+-VsBgCoc", 1800, 900); c->Divide(2, 1); + TH1D* npm = fHMean.H1Clone("hMinvCombPM" + cat + "_sameEv", step); + TH1D* cb = fHMean.H1Clone("hMinvCombBg" + cat, step); + TH1D* bg = fHMean.H1Clone("hMinv" + cat, ELmvmSrc::Bg, step); + TH1D* cocktail = GetCocktailMinvH1("hMinv" + cat, step); + TH1D* sbg = static_cast<TH1D*>(bg->Clone()); + sbg->Add(cocktail); + TH1D* ratS = (TH1D*) npm->Clone(); + ratS->Divide(sbg); + TH1D* ratB = (TH1D*) cb->Clone(); + ratB->Divide(bg); c->cd(1); - DrawH1({pmSame, sbg}, {"N_{same}^{+-}", "Cocktail + BG"}, kLinear, kLog, true, 0.7, 0.8, 0.99, 0.99, "HIST"); - pmSame->SetMinimum(1e-5 * sbg->GetMaximum()); - c->cd(2); - ratio->Divide(sbg); - DrawH1(ratio, kLinear, kLog, "hist"); - ratio->SetMinimum(0.01); - } - } - - //Draw Geometric Mean - { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/geomMean", "cb/geomMean", 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("hMinvCombGeom_sameEv", step); - TH1D* mixed = fHMean.H1Clone("hMinvCombGeom_mixedEv", step); - TH1D* sameU = fHMean.H1Clone("hMinvCombGeom_urqmd_sameEv", step); - TH1D* mixedU = fHMean.H1Clone("hMinvCombGeom_urqmd_mixedEv", step); - int minBin = same->FindBin(0.4); - int maxBin = same->FindBin(0.7); - double scale = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin); - double scaleU = same->Integral(minBin, maxBin) / mixed->Integral(minBin, maxBin); - mixed->Scale(scale); - mixedU->Scale(scaleU); - 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"); - fHMean.DrawAnaStepOnPad(step); - } - } - - // Draw k factor - { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/k", "cb/k", 1800, 1800); - c->Divide(3, 3); - int i = 1; - 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); - double min = 0.6; - double max = 1.4; - k->SetMinimum(min); - kU->SetMinimum(min); - k->SetMaximum(max); - kU->SetMaximum(max); - DrawH1({k, kU}, {"k Factor", "k Factor (UrQMD electrons)"}, kLinear, kLinear, true, 0.65, 0.82, 0.99, 0.95, "p"); - fHMean.DrawAnaStepOnPad(step); - } - } - - // Draw CB vs MC-BG - { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/cbVsMcBg", "cb/cbVsMcBg", 1800, 1800); - c->Divide(3, 3); - int i = 1; - 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); - DrawH1({bg, cb}, {"BG from MC input", "comb. BG (B_{C})"}, kLinear, kLog, true, 0.7, 0.75, 0.99, 0.90, "hist"); + npm->GetYaxis()->SetRangeUser(1e-6, 2e-2); + DrawH1({npm, sbg, cb, bg}, {"N^{+-}_{same}", "BG + Cocktail", "CB_{calc}", "Background"}, kLinear, kLog, true, 0.7, 0.7, 0.99, 0.95, "p"); fHMean.DrawAnaStepOnPad(step); + c->cd(2); + ratS->GetYaxis()->SetTitle("Ratio"); + ratS->GetYaxis()->SetRangeUser(0., 2.); + ratB->GetYaxis()->SetRangeUser(0., 2.); + DrawH1({ratS, ratB}, {"N^{+-}_{same} / BG+Coc", "CB_{calc} / BG"}, kLinear, kLinear, true, 0.7, 0.8, 0.99, 0.95, "p"); } } - // Draw ratio of CB and MC-BG + // Draw Ratios N+-/CocBg and CB/BG for all three categories { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/cbVsMcBg_ratio", "cb/cbVsMcBg_ratio", 1800, 1800); - c->Divide(3, 3); - int i = 1; - 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* 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"); - fHMean.DrawAnaStepOnPad(step); - } + ELmvmAnaStep step = ELmvmAnaStep::ElId; + + TH1D* rsAll = fHMean.H1Clone("hMinvCombPM_sameEv", step); + TH1D* rsUAll = fHMean.H1Clone("hMinvCombPM_urqmdAll_sameEv", step); + TH1D* rsUEl = fHMean.H1Clone("hMinvCombPM_urqmdEl_sameEv", step); + + TH1D* cbgAll = GetCocktailMinvH1("hMinv", step); + cbgAll->Add(fHMean.H1("hMinv", ELmvmSrc::Bg, step)); + TH1D* cbgUAll = GetCocktailMinvH1("hMinv_urqmdAll", step); + cbgUAll->Add(fHMean.H1("hMinv_urqmdAll", ELmvmSrc::Bg, step)); + TH1D* cbgUEl = GetCocktailMinvH1("hMinv_urqmdEl", step); + cbgUEl->Add(fHMean.H1("hMinv_urqmdEl", ELmvmSrc::Bg, step)); + + rsAll->Divide(cbgAll); + rsUAll->Divide(cbgUAll); + rsUEl->Divide(cbgUEl); + + TH1D* rbAll = fHMean.H1Clone("hMinvCombBg", step); + TH1D* rbUAll = fHMean.H1Clone("hMinvCombBg_urqmdAll", step); + TH1D* rbUEl = fHMean.H1Clone("hMinvCombBg_urqmdEl", step); + + rbAll->Divide(fHMean.H1("hMinv", ELmvmSrc::Bg, step)); + rbUAll->Divide(fHMean.H1("hMinv_urqmdAll", ELmvmSrc::Bg, step)); + rbUEl->Divide(fHMean.H1("hMinv_urqmdEl", ELmvmSrc::Bg, step)); + + rsAll->GetYaxis()->SetTitle("N^{+-}_{same} / BG+Coc"); + rsAll->GetYaxis()->SetRangeUser(0., 2.); + rbAll->GetYaxis()->SetTitle("CB_{calc} / BG"); + rbAll->GetYaxis()->SetRangeUser(0., 2.); + + TCanvas* c = fHMean.fHM.CreateCanvas(folder + "CbVsMc_ratio", folder + "CbVsMc_ratio", 1800, 900); + c->Divide(2, 1); + c->cd(1); + DrawH1({rsAll, rsUAll, rsUEl}, {"UrQMD + PLUTO", "UrQMD", "UrQMD electrons"}, kLinear, kLinear, true, 0.55, 0.78, 0.99, 0.95, "hist"); + fHMean.DrawAnaStepOnPad(step); + c->cd(2); + DrawH1({rbAll, rbUAll, rbUEl}, {"UrQMD + PLUTO", "UrQMD", "UrQMD electrons"}, kLinear, kLinear, true, 0.55, 0.78, 0.99, 0.95, "hist"); + fHMean.DrawAnaStepOnPad(step); } // Draw common CB vs. CB from pure UrQMD electrons { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/cbVsUrqmdCb", "cb/cbVsUrqmdCb", 1800, 1800); + TCanvas* c = fHMean.fHM.CreateCanvas(folder + "CbVsUrqmdCb", folder + "CbVsUrqmdCb", 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); TH1D* cb = fHMean.H1Clone("hMinvCombBg", step); - TH1D* cbU = fHMean.H1Clone("hMinvCombBg_urqmd", step); + TH1D* cbU = fHMean.H1Clone("hMinvCombBg_urqmdEl", step); DrawH1({cb, cbU}, {"B_{C}", "B_{C} from UrQMD electrons"}, kLinear, kLog, true, 0.55, 0.8, 0.99, 0.95); fHMean.DrawAnaStepOnPad(step); } } - //Draw Combinatorial Signal from N+-same + //Draw Combinatorial Signal with N+-same, CB, Cocktail { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/cbVsInput_Npm", "cb/cbVsInput_Npm", 1800, 1800); - c->Divide(3, 3); - int i = 1; - for (auto step : fHMean.fAnaSteps) { - if (step < ELmvmAnaStep::ElId) continue; - c->cd(i++); - TH1D* pmSame = fHMean.H1Clone("hMinvCombPM_sameEv", step); - TH1D* combBg = fHMean.H1Clone("hMinvCombBg", step); - TH1D* combSignalNpm = fHMean.H1Clone("hMinvCombSignalNpm", step); - TH1D* cocktail = GetCocktailMinvH1(step); - pmSame->SetMaximum(1e-2); - pmSame->SetMinimum(4e-9); - DrawH1({pmSame, combBg, combSignalNpm, cocktail}, + for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { + string cName = (cat == "_urqmdEl") ? folderUEl : (cat == "_urqmdAll") ? folderUAll : folder; + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "CbVsInput_Npm", cName + "CbVsInput_Npm", 1800, 1800); + c->Divide(3, 3); + int i = 1; + for (auto step : fHMean.fAnaSteps) { + if (step < ELmvmAnaStep::ElId) continue; + c->cd(i++); + TH1D* pmSame = fHMean.H1Clone("hMinvCombPM" + cat + "_sameEv", step); + TH1D* combBg = fHMean.H1Clone("hMinvCombBg" + cat, step); + TH1D* combSignalNpm = fHMean.H1Clone("hMinvCombSignalNpm" + cat, step); + TH1D* cocktail = GetCocktailMinvH1("hMinv" + cat, step); + 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"); - fHMean.DrawAnaStepOnPad(step); + 0.99, 0.92, "p"); + fHMean.DrawAnaStepOnPad(step); + } } } //Draw Combinatorial Signal from Cocktail+BG { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/cbVsInput_cocBg", "cb/cbVsInput_cocBg", 1800, 1800); - c->Divide(3, 3); - int i = 1; - for (auto step : fHMean.fAnaSteps) { - if (step < ELmvmAnaStep::ElId) continue; - c->cd(i++); - TH1D* sbg2 = fHMean.H1Clone("hMinv_bg", step); - sbg2->Add(GetCocktailMinvH1(step)); - TH1D* combBg = fHMean.H1Clone("hMinvCombBg", step); - TH1D* sBgSignal = fHMean.H1Clone("hMinvCombSignalMc", step); - 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"); - fHMean.DrawAnaStepOnPad(step); - } - } - - // Draw ratio of CB-signal and cocktail - { - TCanvas* c = fHMean.fHM.CreateCanvas("cb/SigVsCoc", "cb/SigVsCoc", 1800, 1800); - c->Divide(3, 3); - int i = 1; - for (auto step : fHMean.fAnaSteps) { - if (step < ELmvmAnaStep::ElId) continue; - c->cd(i++); - TH1D* sig = fHMean.H1("hMinvCombSignalMc", step); - TH1D* rat = GetCocktailMinvH1(step); - rat->Divide(sig); - rat->SetMaximum(3.); - DrawH1({rat}, {"Cocktail/S_{CB}"}, kLinear, kLinear, true, 0.5, 0.8, 0.85, 0.86, "hist"); - fHMean.DrawAnaStepOnPad(step); - } - } - - // Compare common CB and signal with versions with pure UrQMD electrons - { - 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); - 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"); - fHMean.DrawAnaStepOnPad(step); + for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { + string cName = (cat == "_urqmdEl") ? folderUEl : (cat == "_urqmdAll") ? folderUAll : folder; + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "CbVsInput_cocBg", cName + "CbVsInput_cocBg", 1800, 1800); + c->Divide(3, 3); + int i = 1; + for (auto step : fHMean.fAnaSteps) { + if (step < ELmvmAnaStep::ElId) continue; + c->cd(i++); + TH1D* sbg2 = fHMean.H1Clone("hMinv" + cat, ELmvmSrc::Bg, step); + sbg2->Add(GetCocktailMinvH1("hMinv" + cat, step)); + TH1D* combBg = fHMean.H1Clone("hMinvCombBg" + cat, step); + TH1D* sBgSignal = fHMean.H1Clone("hMinvCombSignalMc" + cat, step); + sbg2->SetMaximum(2e-2); + sbg2->SetMinimum(4e-9); + TH1D* cocktail = GetCocktailMinvH1("hMinv" + cat, 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"); + fHMean.DrawAnaStepOnPad(step); + } } } } @@ -647,7 +1789,7 @@ void LmvmDrawAll::DrawSBgVsMinv() vector<TH1*> hists; for (ELmvmAnaStep step : {ELmvmAnaStep::TtCut, ELmvmAnaStep::PtCut}) { TH1D* bg = fHMean.H1Clone("hMinv", ELmvmSrc::Bg, step); - TH1D* cocktail = GetCocktailMinvH1(step); + TH1D* cocktail = GetCocktailMinvH1("hMinv", step); TH1D* cocktailOverBg = fHMean.CreateHByClone<TH1D>("hMinv_bg", "hMinvCocktailOverBg", step); cocktailOverBg->GetYaxis()->SetTitle("Cocktail/Background"); cocktailOverBg->Divide(cocktail, bg, 1., 1., "B"); @@ -671,10 +1813,10 @@ 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 : {"", "_urqmdEl", "_urqmdAll"}) { // common CB or from pure UrQMD (all or electrons) if (step < ELmvmAnaStep::ElId) continue; - // Geometrical mean + // Geometric mean for (const string& ev : {"_sameEv", "_mixedEv"}) { TH1D* pp = fHMean.H1Clone("hMinvCombPP" + cat + ev, step); TH1D* mm = fHMean.H1Clone("hMinvCombMM" + cat + ev, step); @@ -683,7 +1825,7 @@ void LmvmDrawAll::CalcCombBGHistos() for (int i = 1; i <= geom->GetNbinsX(); i++) { double cpp = pp->GetBinContent(i); double cmm = mm->GetBinContent(i); - double content = std::sqrt(cpp * cmm); // TODO: mind case if one of these values is 0 + double content = std::sqrt(cpp * cmm); geom->SetBinContent(i, content); } } @@ -716,8 +1858,8 @@ void LmvmDrawAll::CalcCombBGHistos() hSigNpm->Add(fHMean.H1("hMinvCombBg" + cat, step), -1.); // from MC 'Cocktail + BG': Signal = Coc+BG - B_{C} - TH1D* hSigCoc = fHMean.CreateHByClone<TH1D>("hMinv_bg", "hMinvCombSignalMc" + cat, step); - hSigCoc->Add(GetCocktailMinvH1(step)); + TH1D* hSigCoc = fHMean.CreateHByClone<TH1D>("hMinv" + cat + "_bg", "hMinvCombSignalMc" + cat, step); + hSigCoc->Add(GetCocktailMinvH1("hMinv" + cat, step)); hSigCoc->Add(fHMean.H1("hMinvCombBg" + cat, step), -1.); // calculate errors via error propagation formula @@ -760,7 +1902,7 @@ void LmvmDrawAll::CalcCombBGHistos() if (i <= minBin) fHMean.H1("hMinvCombSignalMc" + cat, step)->SetBinError(i, errorSig); if (i > minBin) fHMean.H1("hMinvCombSignalMc" + cat, step)->SetBinError(i, errorSig2); } // error calculation - } // cat ("" or "_urqmd") + } // cat ("" or "_urqmdEl") } // steps } @@ -773,7 +1915,7 @@ void LmvmDrawAll::CalcCutEffRange(double minMinv, double maxMinv) TH1D* hS = fHMean.H1(ss1.str() + "_s"); TH1D* hBg = fHMean.H1(ss1.str() + "_bg"); int x = 1; - TH1D* cocktail = GetCocktailMinvH1(ELmvmAnaStep::ElId); + TH1D* cocktail = GetCocktailMinvH1("hMinv", ELmvmAnaStep::ElId); int binMin = cocktail->FindBin(minMinv); int binMax = cocktail->FindBin(maxMinv); double sIntElId = cocktail->Integral(binMin, binMax); @@ -782,7 +1924,7 @@ void LmvmDrawAll::CalcCutEffRange(double minMinv, double maxMinv) if (step < ELmvmAnaStep::ElId) continue; if (!fUseMvd && (step == ELmvmAnaStep::Mvd1Cut || step == ELmvmAnaStep::Mvd2Cut)) continue; - double effS = 100. * GetCocktailMinvH1(step)->Integral(binMin, binMax) / sIntElId; + double effS = 100. * GetCocktailMinvH1("hMinv", step)->Integral(binMin, binMax) / sIntElId; double effB = 100. * fHMean.H1("hMinv", ELmvmSrc::Bg, step)->Integral(binMin, binMax) / bgIntElId; hBg->GetXaxis()->SetBinLabel(x, fHMean.fAnaStepLatex[static_cast<int>(step)].c_str()); @@ -818,14 +1960,14 @@ TH1D* LmvmDrawAll::SBgRange(double minMinv, double maxMinv) TH1D* hSBg = fHMean.H1(ss.str()); hSBg->GetXaxis()->SetLabelSize(0.06); int x = 1; - TH1D* cocktail = GetCocktailMinvH1(ELmvmAnaStep::ElId); + TH1D* cocktail = GetCocktailMinvH1("hMinv", ELmvmAnaStep::ElId); int binMin = cocktail->FindBin(minMinv); int binMax = cocktail->FindBin(maxMinv); for (ELmvmAnaStep step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; if (!fUseMvd && (step == ELmvmAnaStep::Mvd1Cut || step == ELmvmAnaStep::Mvd2Cut)) continue; - double intS = 100. * GetCocktailMinvH1(step)->Integral(binMin, binMax); + double intS = 100. * GetCocktailMinvH1("hMinv", step)->Integral(binMin, binMax); double intBg = 100. * fHMean.H1("hMinv", ELmvmSrc::Bg, step)->Integral(binMin, binMax); double sbg = intS / intBg; diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h old mode 100644 new mode 100755 index d13c615132..ede49439e5 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h @@ -37,7 +37,7 @@ private: std::vector<LmvmHist*> fH; LmvmHist fHMean; - int fRebMinv = 50; // Rebin for minv histograms + int fRebMinv = 5; // Rebin for minv histograms (binSize in Task.cxx is 10 MeV) std::string fOutputDir; // output directory for figures @@ -49,7 +49,7 @@ private: template<class T> T* GetCocktailMinv(const std::string& name, ELmvmAnaStep step); - TH1D* GetCocktailMinvH1(ELmvmAnaStep step); + TH1D* GetCocktailMinvH1(const std::string& name, ELmvmAnaStep step); /** * \brief Draw S/Bg vs minv. @@ -58,14 +58,19 @@ private: void DrawMinvAll(); void DrawMinv(ELmvmAnaStep step); void DrawMinvPtAll(); + void DrawMomRecoPrecision(); void DrawMinvBgSourcesAll(); /** * \brief Draw invariant mass spectra for all signal types for specified analysis step with BG reduced by combinatorial BG. - * \param[in] step Analysis step. */ void DrawMinvCombBgAndSignal(); + /** + * \brief Draw invariant mass spectra in official style. + */ + void DrawMinvOfficialStyle(); + template<class T> void CreateMeanHist(const std::string& name, int nofEvents, int nofRebins = -1); // if nRebin = -1, no rebin void CreateMeanHistAll(); @@ -99,6 +104,40 @@ private: */ void SBgRangeAll(); + /** + * \brief Draw properties of misidentified particles. + */ + void InvestigateMisid(); + + void DrawBetaMomSpectra(); + + void DrawMomPluto(); + + void DrawTofM2(); + + void DrawGTrackVertex(); + + void DrawSignificancesAll(); + void DrawSignificance(TH2D* hEl, TH2D* hBg, const std::string& name, double minZ, double maxZ, const std::string& option); + + void DrawCutEffSignal(); + + void DrawPionSuppression(); + + // investigate misidentifications + void DrawMomentum(); + void DrawPtYAndTofM2Elid(); + void DrawTofPilePids(); + void DrawTofHitXY(); + void DrawPurity(); + void DrawPurityHistText(TH2D* h); + void DrawChi2(); + + void DrawMinvScaleValues(); + + /** + * \brief Draw properties of misidentified particles in comparison with not-misidentified. + */ void DrawSBgResults(); @@ -109,6 +148,7 @@ private: */ void SaveCanvasToImage(); + double fZ = -44.; // z-position of target LmvmDrawAll(const LmvmDrawAll&); LmvmDrawAll operator=(const LmvmDrawAll&); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx old mode 100644 new mode 100755 index 42d5cc6b27..957d44f927 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx @@ -11,6 +11,10 @@ #include "Logger.h" #include "TText.h" +#include "TStyle.h" +#include "TLegend.h" + +#include "LmvmUtils.h" using std::string; @@ -34,7 +38,7 @@ const vector<int> LmvmHist::fAnaStepColors = {kGreen + 3, kOrange + 3, kBlack, kRed, kPink - 6, kGreen, kOrange - 3, kViolet + 10, kGreen - 3, kMagenta, kYellow + 1}; -const vector<std::string> LmvmHist::fSignalNames = {"immed", "qgp", "omega", "phi", "omegaD"}; +const vector<std::string> LmvmHist::fSignalNames = {"inmed", "qgp", "omega", "phi", "omegaD"}; const vector<ELmvmSignal> LmvmHist::fSignals = {ELmvmSignal::Inmed, ELmvmSignal::Qgp, ELmvmSignal::Omega, ELmvmSignal::Phi, ELmvmSignal::OmegaD}; @@ -42,6 +46,17 @@ const vector<std::string> LmvmHist::fBgPairSrcNames = {"GG", "PP", "OO", "GP", " const vector<std::string> LmvmHist::fBgPairSrcLatex = {"#gamma-#gamma", "#pi^{0}-#pi^{0}", "o.-o.", "#gamma-#pi^{0}", "#gamma-o.", "#pi^{0}-o."}; +// the following candidates and global track vectors are mainly set up to investigate misidentifications +const vector<std::string> LmvmHist::fGTrackNames = { + "el+", "el+_prim", "el-", "el-_prim", "pion+", "pion+_prim", "pion-", "pion-_prim", + "proton", "proton_prim", "kaon+", "kaon+_prim", "kaon-", "kaon-_prim", "other"}; +const vector<std::string> LmvmHist::fGTrackLatex = { + "e^{+}", "e^{+}_{prim}", "e^{-}", "e^{-}_{prim}", "#pi^{+}", "#pi^{+}_{prim}", "#pi^{-}", "#pi^{-}_{prim}", + "p", "p_{prim}", "K^{+}", "K^{+}_{prim}", "K^{-}", "K^{-}_{prim}", "other"}; + +const vector<std::string> LmvmHist::fCandNames = {"plutoEl+", "plutoEl-", "urqmdEl+", "urqmdEl-", "pion+", "pion-", "proton", "kaon+", "kaon-", "other"}; +const vector<std::string> LmvmHist::fCandLatex = {"e^{+}_{PLUTO}", "e^{-}_{PLUTO}", "e^{+}_{UrQMD}", "e^{-}_{UrQMD}", "#pi^{+}", "#pi^{-}", "p", "K^{+}", "K^{-}", "o."}; + LmvmHist::LmvmHist() {} @@ -130,15 +145,17 @@ void LmvmHist::FillH2(const string& name, double x, double y, double w) { H2(nam void LmvmHist::FillH1(const string& name, ELmvmSrc src, double x, double wSignal) { if (src == ELmvmSrc::Undefined) return; - double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; - H1(name, src)->Fill(x, myWeight); + //double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; // TODO: delete commented lines + //H1(name, src)->Fill(x, myWeight); + H1(name, src)->Fill(x, wSignal); } void LmvmHist::FillH2(const string& name, ELmvmSrc src, double x, double y, double wSignal) { if (src == ELmvmSrc::Undefined) return; - double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; - H2(name, src)->Fill(x, y, myWeight); + //double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; // TODO: delete commented lines + //H2(name, src)->Fill(x, y, myWeight); + H2(name, src)->Fill(x, y, wSignal); } void LmvmHist::FillH1(const string& name, ELmvmAnaStep step, double x, double w) { H1(name, step)->Fill(x, w); } @@ -151,15 +168,17 @@ void LmvmHist::FillH2(const string& name, ELmvmAnaStep step, double x, double y, void LmvmHist::FillH1(const string& name, ELmvmSrc src, ELmvmAnaStep step, double x, double wSignal) { if (src == ELmvmSrc::Undefined || step == ELmvmAnaStep::Undefined) return; - double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; - FillH1(GetName(name, src, step), x, myWeight); + //double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; // TODO: delete commented lines + //FillH1(GetName(name, src, step), x, myWeight); + FillH1(GetName(name, src, step), x, wSignal); } void LmvmHist::FillH2(const string& name, ELmvmSrc src, ELmvmAnaStep step, double x, double y, double wSignal) { if (src == ELmvmSrc::Undefined || step == ELmvmAnaStep::Undefined) return; - double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; - FillH2(GetName(name, src, step), x, y, myWeight); + //double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; // TODO: delete commented lines + //FillH2(GetName(name, src, step), x, y, myWeight); + FillH2(GetName(name, src, step), x, y, wSignal); } @@ -186,6 +205,91 @@ string LmvmHist::GetName(const string& name, ELmvmSrc src, ELmvmAnaStep step) return GetName(GetName(name, src), step); } +void LmvmHist::SetOptH1(TH1D *hist, TString xAxisTitle=" ", TString yAxisTitle=" ", Int_t Ndevision=510, Int_t style=1, Float_t size=2, Int_t color=1, string opt) +{ + hist->GetXaxis()->SetTitle(xAxisTitle); + hist->GetYaxis()->SetTitle(yAxisTitle); + hist->GetXaxis()->SetTitleSize(0.05); + hist->GetYaxis()->SetTitleSize(0.05); + hist->GetXaxis()->SetTitleFont(42); + hist->GetYaxis()->SetTitleFont(42); + hist->GetXaxis()->SetNdivisions(Ndevision); + hist->GetYaxis()->SetTitleOffset(1.4); + hist->GetXaxis()->SetTitleOffset(1.); + //hist->GetXaxis()->CenterTitle(); + //hist->GetYaxis()->CenterTitle(); + hist->GetXaxis()->SetLabelFont(42); + hist->GetYaxis()->SetLabelFont(42); + hist->GetXaxis()->SetLabelSize(0.04); + hist->GetYaxis()->SetLabelSize(0.04); + if (opt == "marker") { + hist->SetMarkerStyle(style); + hist->SetMarkerSize(size); + hist->SetMarkerColor(color); + hist->SetLineColor(color); + hist->SetLineWidth(2); + hist->SetTitle(""); + //hist->SetLineWidth(2); + //hist->SetMinimum(0.1); + //hist->SetMaximum(1e5); + hist->GetXaxis()->SetLabelColor(1); + hist->GetYaxis()->SetLabelColor(1); + hist->GetXaxis()->SetTitleColor(1); + hist->GetYaxis()->SetTitleColor(1); + hist->SetFillColor(1); + } + else if (opt == "line") { + hist->SetLineStyle(style); + hist->SetLineColor(color); + hist->SetTitle(""); // TODO: with ""? + hist->SetLineWidth(size); + } + else LOG(error) << "Option '" << opt << "' undefined. Choose 'marker' or 'line'." << std::endl; +} + +void LmvmHist::SetOptCanvas(TCanvas *canvas) +{ + gStyle->SetOptStat(0); + gStyle->SetEndErrorSize(5); + //gStyle->SetErrorX(0); // X errors of the data points set to be 0 + gStyle->SetLineStyleString(22,"80 18 12 18 12 12"); // special style for the line + gStyle->SetEndErrorSize(5); // define end width of error bars + gStyle->SetCanvasColor(10); + gStyle->SetPadColor(10); + canvas->SetLeftMargin(0.15); + canvas->SetRightMargin(0.04); + canvas->SetTopMargin(0.05); + canvas->SetBottomMargin(0.12); + canvas->ToggleEventStatus(); + canvas->Range(-200,-10,1000,-2); + canvas->SetFillColor(0); + canvas->SetBorderMode(0); + canvas->SetBorderSize(0); + canvas->SetTickx(); + canvas->SetTicky(); + canvas->SetLogy(); + canvas->SetFrameLineWidth(2); + canvas->SetFrameBorderMode(0); + canvas->SetFrameBorderSize(0); +} + +void LmvmHist::SetLegend(vector<LmvmLegend> legendV, double textsize, Double_t x1, Double_t y1, Double_t x2, Double_t y2) +{ + auto leg = new TLegend(x1, y1, x2, y2); + leg->SetLineColor(10); + leg->SetLineStyle(1); + leg->SetLineWidth(1); + leg->SetFillColor(10); + leg->SetTextSize(textsize); + leg->SetTextFont(42); + leg->SetMargin(0.4); + leg->SetBorderSize(0); + for (size_t i = 0; i < legendV.size(); i++) { + const auto& l = legendV[i]; + leg->AddEntry(l.fH, l.fName, l.fOpt); + } + leg->Draw("same"); +} void LmvmHist::Rebin(const string& name, int nGroup) { fHM.Rebin(name, nGroup); } @@ -205,6 +309,153 @@ void LmvmHist::Rebin(const string& name, const vector<string>& subNames1, const } } +TH1D* LmvmHist::CreateSignificanceH1(TH1D* s, TH1D* bg, const string& name, const string& option) +{ + int nBins = s->GetNbinsX(); + TH1D* hsig = new TH1D(name.c_str(), name.c_str(), nBins, s->GetXaxis()->GetXmin(), s->GetXaxis()->GetXmax()); + hsig->GetXaxis()->SetTitle(s->GetXaxis()->GetTitle()); + + // "right" - reject right part of the histogram. value > cut -> reject + if (option == "right") { + for (int i = 1; i <= nBins; i++) { + double sumSignal = s->Integral(1, i, "width"); + double sumBg = bg->Integral(1, i, "width"); + double sign = (sumSignal + sumBg != 0.) ? sumSignal / std::sqrt(sumSignal + sumBg) : 0.; + hsig->SetBinContent(i, sign); + hsig->GetYaxis()->SetTitle("Significance"); + } + } + // "left" - reject left part of the histogram. value < cut -> reject + else if (option == "left") { + for (int i = nBins; i >= 1; i--) { + double sumSignal = s->Integral(i, nBins, "width"); + double sumBg = bg->Integral(i, nBins, "width"); + double sign = (sumSignal + sumBg != 0.) ? sumSignal / std::sqrt(sumSignal + sumBg) : 0.; + hsig->SetBinContent(i, sign); + hsig->GetYaxis()->SetTitle("Significance"); + } + } + return hsig; +} + +TH2D* LmvmHist::CreateSignificanceH2(TH2D* signal, TH2D* bg, const string& name, const string& title) +{ + double xmin = 1.0; + double xmax = 5.0; + double ymin = 1.0; + double ymax = 5.0; + double delta = 0.1; + int nStepsX = (int) ((xmax - xmin) / delta); + int nStepsY = (int) ((ymax - ymin) / delta); + + TH2D* hsig = new TH2D(name.c_str(), title.c_str(), nStepsX, xmin, xmax, nStepsY, ymin, ymax); + + int binX = 1; + for (double xcut = xmin; xcut <= xmax; xcut += delta, binX++) { + int binY = 1; + for (double ycut = ymin; ycut <= ymax; ycut += delta, binY++) { + double sumSignal = 0; + double sumBg = 0; + for (int ix = 1; ix <= signal->GetNbinsX(); ix++) { + for (int iy = 1; iy <= signal->GetNbinsY(); iy++) { + double xc = signal->GetXaxis()->GetBinCenter(ix); + double yc = signal->GetYaxis()->GetBinCenter(iy); + double val = -1 * (ycut / xcut) * xc + ycut; + + if (!(xc < xcut && yc < val)) { + sumSignal += signal->GetBinContent(ix, iy); + sumBg += bg->GetBinContent(ix, iy); + } + } + } + double sign = (sumSignal + sumBg != 0.) ? sumSignal / std::sqrt(sumSignal + sumBg) : 0.; + hsig->SetBinContent(binX, binY, sign); + } + } + return hsig; +} + +void LmvmHist::DrawAllGTracks(int dim, const string& cName, const string& hName, vector<string> xLabel, vector<string> yLabel, double min, double max) +{ + TCanvas* c = fHM.CreateCanvas(cName, cName, 2400, 2400); + c->Divide(4, 4); + int i = 1; + for (const string& ptcl : fGTrackNames) { + c->cd(i++); + string hFullname = hName + "_" + ptcl; + DrawAll(dim, hFullname, fGTrackLatex[i-2], xLabel, yLabel, min, max); + } +} + +void LmvmHist::DrawAllCands(int dim, const string& cName, const string& hName, vector<string> xLabel, vector<string> yLabel, double min, double max) +{ + TCanvas* c = fHM.CreateCanvas(cName, cName, 2400, 1800); + c->Divide(4, 3); + int i = 1; + for (const string& ptcl : fCandNames) { + c->cd(i++); + string hFullname = hName + "_" + ptcl; + DrawAll(dim, hFullname, fCandLatex[i-2], xLabel, yLabel, min, max); + } +} + +void LmvmHist::DrawAllCandsAndSteps(int dim, const string& cName, const string& hName, vector<string> xLabel, vector<string> yLabel, double min, double max) +{ + for (auto step : fAnaSteps) { + TCanvas* c = fHM.CreateCanvas(GetName(cName + "cands_all", step), GetName(cName + "cands_all", step), 2400, 1800); + c->Divide(4, 3); + int i = 1; + for (const string& ptcl : fCandNames) { + c->cd(i++); + string hFullname = GetName(hName + "_" + ptcl, step); + DrawAll(dim, hFullname.c_str(), fCandLatex[i-2], xLabel, yLabel, min, max); + } + } + + for (const string& ptcl : fCandNames) { + TCanvas* c = fHM.CreateCanvas(cName + ptcl + "_steps", cName + ptcl + "_steps", 2400, 1800); + c->Divide(4, 3); + int i = 1; + for (auto step : fAnaSteps) { + c->cd(i++); + string hFullname = GetName(hName + "_" + ptcl, step); + DrawAll(dim, hFullname.c_str(), fAnaStepNames[static_cast<int>(step)], xLabel, yLabel, min, max); + } + } +} + +void LmvmHist::DrawAll(int dim, const string& hFullname, const string& text, vector<string> xLabel, vector<string> yLabel, double min, double max) +{ + if (dim == 1) { + TH1D* h = H1Clone(hFullname.c_str()); + h->GetYaxis()->SetRangeUser(min, max); + DrawH1(h, kLinear, kLog, "hist"); + if (xLabel.size() > 1) { + for (size_t iL = 0; iL < xLabel.size(); iL++) { + h->GetXaxis()->SetBinLabel(iL + 1, xLabel[iL].c_str()); + } + } + } + else if (dim == 2) { + TH2D* h = H2Clone(hFullname.c_str()); + h->GetZaxis()->SetRangeUser(min, max); + DrawH2(h, kLinear, kLinear, kLog, "colz"); + if (xLabel.size() > 1) { + for (size_t iL = 0; iL < xLabel.size(); iL++) { + h->GetXaxis()->SetBinLabel(iL + 1, xLabel[iL].c_str()); + } + } + if (yLabel.size() > 1) { + for (size_t iL = 0; iL < yLabel.size(); iL++) { + h->GetYaxis()->SetBinLabel(iL + 1, yLabel[iL].c_str()); + } + } + } + else LOG(error) << "LmvmHist::DrawAll: Choose dimension 1 or 2"; + + DrawTextOnPad(text, 0.25, 0.9, 0.75, 0.999); +} + void LmvmHist::WriteToFile() { fHM.WriteToFile(); } void LmvmHist::DrawEfficiency(TH1* h1, TH1* h2, double xPos, double yPos) diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h old mode 100644 new mode 100755 index fd431835b7..f365103361 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h @@ -42,6 +42,15 @@ public: const static std::vector<std::string> fBgPairSrcNames; const static std::vector<std::string> fBgPairSrcLatex; + const static int fNofGTrackNames = 15; + const static std::vector<std::string> fGTrackNames; + const static std::vector<std::string> fGTrackLatex; + + const static int fNofCandNames = 10; + const static std::vector<std::string> fCandNames; + const static std::vector<std::string> fCandLatex; + + std::vector<std::string> CombineNames(const std::string& name, const std::vector<std::string>& subNames); std::vector<std::string> CombineNames(const std::string& name, const std::vector<std::string>& subNames1, @@ -125,11 +134,21 @@ public: std::string GetName(const std::string& name, ELmvmSrc src); std::string GetName(const std::string& name, ELmvmSrc src, ELmvmAnaStep step); + void SetOptH1(TH1D *hist, TString xAxisTitle, TString yAxisTitle, Int_t Ndevision, Int_t style, Float_t size, Int_t color, std::string opt = ""); // copied from Tetyanas macro + void SetOptCanvas(TCanvas *canvas); + void SetLegend(std::vector<LmvmLegend>, double textsize, double x1, double y1, double x2, double y2); - void Rebin(const std::string& name, int nGroup); + void Rebin(const std::string& name, int nGroup); // TODO: used? void Rebin(const std::string& name, const std::vector<std::string>& subNames, int nGroup); - void Rebin(const std::string& name, const std::vector<std::string>& subNames1, - const std::vector<std::string>& subNames2, int nGroup); + void Rebin(const std::string& name, const std::vector<std::string>& subNames1, const std::vector<std::string>& subNames2, int nGroup); + + TH1D* CreateSignificanceH1(TH1D* s, TH1D* bg, const std::string& name, const std::string& option); + TH2D* CreateSignificanceH2(TH2D* signal, TH2D* bg, const std::string& name, const std::string& title); + + void DrawAll(int dim, const std::string& hFullname, const std::string& padText, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); + void DrawAllGTracks(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); + void DrawAllCands(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); + void DrawAllCandsAndSteps(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); void WriteToFile(); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h b/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h index 29df08cf23..414830101a 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmSimParam.h @@ -23,8 +23,7 @@ public: if (particle == "omegadalitz" || particle == "wdalitz") return 2.28721 * 7.7e-4; if (particle == "phi") return 0.311619 * 2.97e-4; if (particle == "inmed" || particle == "rho0") return 0.0304706; - if (particle == "qgp" || particle == "qgp_epem") - return 10 * 4.52941e-4; // TODO urgent: delete factor 10 and set right parameter!!! + if (particle == "qgp" || particle == "qgp_epem") return 4.52941e-4; } 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 old mode 100644 new mode 100755 index c8d9cdfc55..149e86bac5 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx @@ -13,14 +13,17 @@ #include "CbmRichHit.h" #include "CbmRichPoint.h" #include "CbmRichRing.h" +#include "CbmRichUtil.h" + #include "CbmStsHit.h" #include "CbmStsKFTrackFitter.h" #include "CbmStsTrack.h" +#include "CbmTrdHit.h" +#include "CbmTrdTrack.h" #include "CbmTofHit.h" #include "CbmTofPoint.h" +#include "CbmTofTrack.h" #include "CbmTrackMatchNew.h" -#include "CbmTrdHit.h" -#include "CbmTrdTrack.h" #include "CbmVertex.h" #include "cbm/elid/CbmLitGlobalElectronId.h" #include "cbm/qa/mc/CbmLitMCTrackCreator.h" @@ -34,6 +37,7 @@ #include "TClonesArray.h" #include "TDatabasePDG.h" +#include "TMCProcess.h" #include "TFile.h" #include "TRandom3.h" #include "TVector3.h" @@ -43,7 +47,7 @@ #include "LmvmSimParam.h" #include "LmvmUtils.h" - +#include "LmvmHist.h" using namespace std; @@ -59,18 +63,15 @@ 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 = "dN/dM_{ee}/N [GeV/c^{2}]^{-1}"; 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.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", "Yield/(Event * Bin)", 200, 0., 10., 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); @@ -78,9 +79,15 @@ void LmvmTask::InitHists() fH.CreateH2("hVertexGammaXZ", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 200, -10., 190., 400, -130., 130.); fH.CreateH2("hVertexGammaYZ", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 200, -10., 190., 400, -130., 130.); fH.CreateH2("hVertexGammaXY", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 400, -130., 130., 400, -130., 130.); - fH.CreateH2("hVertexGammaRZ", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 300, -10., 190., 300, 0., - 150.); - + fH.CreateH2("hVertexGammaRZ", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 50, 0., + 50.); + + // vertices of not-electron candidates that are misidentified as electrons in ToF + fH.CreateH2("hVertexXZ_misidTof", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.); + fH.CreateH2("hVertexYZ_misidTof", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.); + fH.CreateH2("hVertexXY_misidTof", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 100, -50., 50., 100, -50., 50.); + fH.CreateH2("hVertexRZ_misidTof", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 100, 0., + 50.); fH.CreateH1("hNofBgTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); fH.CreateH1("hNofSignalTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); fH.CreateH2("hBgSrcTracks", "Analysis step", "Candidate Source", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.); @@ -88,6 +95,14 @@ void LmvmTask::InitHists() fH.CreateH1("hNofTopoPairs", {"gamma", "pi0"}, "Pair type", "Pairs/event", 8, 0., 8); fH.CreateH1("hNofMismatches", {"all", "rich", "trd", "tof"}, "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); + fH.CreateH1("hNofMismatches_gTracks", {"all", "complete"}, fH.fGTrackNames, "Detector", ax, 7., 0., 7.); + fH.CreateH1("hNofMismatchedTrackSegments", {"all", "complete"}, fH.fGTrackNames, "nof mism. track segments", ax, 4., 0., 4.); //"all": also partial tracks; "complete": only tracks with entries in all detectors + fH.CreateH2("hMatchId_gTracks", fH.fGTrackNames, "Nof mismatched Track Segments", "Identification", ax, 4., 0., 4., 2., 0., 2.); + + fH.CreateH1("hChi2_mismatch_all", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); + fH.CreateH1("hChi2_truematch_all", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); + fH.CreateH1("hChi2_mismatch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); + fH.CreateH1("hChi2_truematch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); fH.CreateH1("hNofGhosts", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); @@ -100,18 +115,18 @@ void LmvmTask::InitHists() fH.CreateH1("hEventNumberMixed", "", "", 1, 0, 1.); fH.CreateH1("hAnnRich", fH.fSrcNames, "RICH ANN output", ax, 100, -1.1, 1.1); - fH.CreateH2("hAnnRichVsMom", fH.fSrcNames, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 6., 100, -1.1, 1.1); + fH.CreateH2("hAnnRichVsMom", fH.fSrcNames, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 10., 100, -1.1, 1.1); fH.CreateH1("hAnnTrd", fH.fSrcNames, "Likelihood output", ax, 100, -.1, 1.1); // TODO: change back to "TRD ANN" fH.CreateH2("hTofM2", fH.fSrcNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 150, 0., 6., 500, -0.1, 1.0); + fH.CreateH1("hChi2Sts", fH.fSrcNames, "#chi^{2}", ax, 200, 0., 20.); fH.CreateH1("hChi2PrimVertex", fH.fSrcNames, "#chi^{2}_{prim}", ax, 200, 0., 20.); fH.CreateH1("hNofMvdHits", fH.fSrcNames, "Number of hits in MVD", ax, 5, -0.5, 4.5); fH.CreateH1("hNofStsHits", fH.fSrcNames, "Number of hits in STS", ax, 9, -0.5, 8.5); fH.CreateH2("hTrdLike", {"El", "Pi"}, fH.fSrcNames, "P [GeV/c]", "Likelihood output", ax, 100, 0., 6., 100, 0., 1.); - fH.CreateH2("hAnnRichVsMomPur", {"El", "Bg"}, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 6., 100, -1.1, 1.1); - fH.CreateH2("hTrdElLikePur", {"El", "Bg"}, "P [GeV/c]", "Likelihood output", ax, 100, 0., 6., 100, 0., 1.); - fH.CreateH2("hTofM2Pur", {"El", "Bg"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 150, 0., 6., 500, -0.1, 1.0); + fH.CreateH2("hAnnRichVsMomPur", {"El", "Bg"}, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 10., 100, -1.1, 1.1); + fH.CreateH2("hTrdElLikePur", {"El", "Bg"}, "P [GeV/c]", "Likelihood output", ax, 100, 0., 10., 100, 0., 1.); fH.CreateH2("hTtCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]", "#theta_{e^{+},e^{-}} [deg]", ax, 100, 0., 5., 100, 0., 5.); @@ -127,33 +142,29 @@ void LmvmTask::InitHists() 2.); // [0.5]-correct, [1.5]-wrong fH.CreateH1("hMvdMcDist", {"1", "2"}, fH.fSrcNames, "Track-Hit distance [cm]", ax, 100, 0., 10.); - fH.CreateH1("hMinv", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0., 2.5); - - // Combinatorial BG histograms - fH.CreateH1("hMinvCombPM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0., 2.5); - fH.CreateH1("hMinvCombPP", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0., 2.5); - fH.CreateH1("hMinvCombMM", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0., 2.5); - - // Combinatorial BG histograms to distinguish between PLUTO and UrQMD particles - fH.CreateH1("hMinvCombPM_pluto", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0, 2.5); - fH.CreateH1("hMinvCombPP_pluto", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0, 2.5); - fH.CreateH1("hMinvCombMM_pluto", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0, 2.5); - fH.CreateH1("hMinvCombPM_urqmd", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0, 2.5); - fH.CreateH1("hMinvCombPP_urqmd", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0, 2.5); - fH.CreateH1("hMinvCombMM_urqmd", {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0, 2.5); + fH.CreateH1("hMinv", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); + fH.CreateH1("hMinv_urqmdAll", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); // for UrQMD particles only + fH.CreateH1("hMinv_urqmdEl", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); // for UrQMD electrons only + // Pair Yield histograms for combinatorial BG calculation + for (const string& comb : {"PM", "PP", "MM"}) { + for (const string& cat : {"", "_urqmdAll", "_urqmdEl"}) { + string hName = "hMinvComb" + comb + cat; + fH.CreateH1(hName, {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); + } + } + fH.CreateH1("hMinvBgMatch", {"trueMatch", "trueMatchEl", "trueMatchNotEl", "mismatch"}, fH.fAnaStepNames, - "M_{ee} [GeV/c^{2}]", ax, 2000, 0., 2.5); - fH.CreateH1("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 2500, 0., 2.5); + "M_{ee} [GeV/c^{2}]", ax, 250, 0., 2.5); + fH.CreateH1("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 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 + "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); // "pi" are misid. charged pions fH.CreateH2("hMinvPt", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", ax, 100, 0., 2., 25, 0., 2.5); fH.CreateH1("hMomPairSignal", fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 15.); fH.CreateH2("hPtYPairSignal", fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); - fH.CreateH2("hPtYCandidate", fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); fH.CreateH1("hAnglePair", fH.fSrcNames, fH.fAnaStepNames, "#Theta_{1,2} [deg]", ax, 160, 0., 80.); for (const string& suff : {"", "+", "-"}) { @@ -168,13 +179,71 @@ void LmvmTask::InitHists() fH.CreateH1("hMomAcc+", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.); fH.CreateH1("hMomAcc-", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.); - fH.CreateH1("hPiMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}, - "P [GeV/c]", ax, 30, 0., 3.); fH.CreateH1("hElMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}, - "P [GeV/c]", ax, 30, 0., 3.); + "P [GeV/c]", ax, 100, 0., 10.); + fH.CreateH1("hPiMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}, + "P [GeV/c]", ax, 100, 0., 10.); fH.CreateH1("hProtonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}, - "P [GeV/c]", ax, 30, 0., 3.); - fH.CreateH1("hCandMisIdAsEl", {"pi", "pi0", "proton", "kaon"}, "P [GeV/c]", ax, 30, 0., 3.); + "P [GeV/c]", ax, 100, 0., 10.); + fH.CreateH1("hKaonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}, + "P [GeV/c]", ax, 100, 0., 10.); + + fH.CreateH1("hPionSupp", {"pi", "idEl"},fH.fAnaStepNames, "P [GeV/c]", ax, 20, 0., 10.); + + fH.CreateH1("hMom_gTracks", fH.fGTrackNames, "P [GeV/c]", ax, 100, 0., 10.); + fH.CreateH1("hMom_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 10.); + fH.CreateH2("hPtY_gTracks", fH.fGTrackNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); + fH.CreateH2("hPtY_cands", fH.fCandNames, fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); + fH.CreateH2("hTofM2_gTracks", fH.fGTrackNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.); + fH.CreateH2("hTofM2_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.); + + fH.CreateH1("hTofPilePdgs_cands", fH.fAnaStepNames, "Particle", ax, 12, 0., 12.); + fH.CreateH1("hTofPilePdgs_gTracks", "Particle", ax, 12, 0., 12.); + fH.CreateH2("hTofPileHitXY", fH.fCandNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); + fH.CreateH2("hTofPilePointXY", fH.fCandNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); + fH.CreateH1("hTofPileHitPointDist", fH.fCandNames, "#sqrt{dX^{2} + dY^{2}} [cm]", ax, 1000., 0., 20.); + fH.CreateH2("hTofPilePty_cands", fH.fCandNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); + + // ToF Hits + fH.CreateH2("hTofPointXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); // to see if maybe misid/mismatches stem from a certain region in ToF + fH.CreateH2("hTofHitXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); + fH.CreateH1("hTofHitPointDist", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "#sqrt{dX^{2} + dY^{2}} [cm]", ax, 400., 0., 20.); + fH.CreateH2("hTofHitTrackDist_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]", "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.); + fH.CreateH2("hTofHitTrackDist_cands", fH.fCandNames, "P [GeV/c]", "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.); + + fH.CreateH1("hMomRatio_cands", fH.fCandNames, fH.fAnaStepNames, "Ratio P_{rec}/P_{MC}", ax, 100, 0., 2.); + fH.CreateH2("hMomRatioVsMom_cands", fH.fCandNames, fH.fAnaStepNames, "P_{MC} [GeV/c]", "Ratio P_{rec}/P_{MC}", ax, 100, 0., 10., 200., 0., 2.); + + // compare misid. candidates (pions, proton, ...) with not misid. ones + for (size_t iP = 4; iP < fH.fCandNames.size(); iP++) { + fH.CreateH1("hMom_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", ax, 100, 0., 10.); + fH.CreateH2("hPtY_" + fH.fCandNames[iP], {"true", "misid"}, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 40, 0., 2.); + fH.CreateH2("hTofM2_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.); + } + + // beta-momentum spectra + fH.CreateH2("hBetaMom", {"gTracks", "cands"}, fH.fCandNames, "P * Q ", "beta", ax, 200, -10., 10., 150., 0., 1.5); + fH.CreateH2("hBetaMom_allGTracks", "P * Q ", "beta", ax, 200, -10., 10., 150., 0., 1.5); + + // z vertex of global tracks + fH.CreateH2("hVertexGTrackRZ", "Z [cm]", "R [cm]", ax, 1500, -50., 100., 1000, -50., 50.); + + /*fH.CreateH2("hIndexStsRich", fH.fGTrackNames, "STS Index", "RICH Index", ax, 400., 0., 400., 100., 0., 100.); // TODO: can be deleted + fH.CreateH2("hIndexStsTrd", fH.fGTrackNames, "STS Index", "TRD Index", ax, 400., 0., 400., 400., 0., 400.); + fH.CreateH2("hIndexStsTof", fH.fGTrackNames, "STS Index", "ToF Index", ax, 400., 0., 400., 800., 0., 800.);*/ + + fH.CreateH2("hPdgVsMom_gTracks", {"rich", "trd", "tof"}, {"all", "complete"}, "P [GeV/c]", "Misid. Particle", "Yield/(Event * Bin)", 200, 0., 10., 6, 0., 6.); + + fH.CreateH2("hTofM2Calc_gTracks", fH.fGTrackNames, "Time", "", "Yield", 100., 5., 15., 4., 0., 4.); + fH.CreateH2("hTofTimeVsMom_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]", "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., 15.); + fH.CreateH2("hTofTimeVsMom_cands", fH.fCandNames, "P [GeV/c]", "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., 15.); + fH.CreateH2("hRichRingTrackDist_gTracks", fH.fGTrackNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., 10., 200., 0., 20.); + fH.CreateH2("hRichRingTrackDist_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., 10., 200., 0., 20.); + + // check chi2 of candidates after El-ID step + fH.CreateH2("hChi2VsMom", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "P [GeV/c]", "#chi^{2}", ax, 100., 0., 10., 200, 0., 20.); + fH.CreateH2("hTofTimeVsChi2", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "#chi^{2}", "Time (ct) [m]", ax, 200., 0., 20., 200, 0., 20.); + fH.CreateH2("hChi2Comb", {"StsRich", "StsTrd", "StsTof", "RichTrd", "RichTof", "TrdTof"}, fH.fCandNames, "#chi^{2}_{1}", "#chi^{2}_{2}", ax, 200., 0., 20., 200, 0., 20.); } InitStatus LmvmTask::Init() @@ -202,6 +271,11 @@ InitStatus LmvmTask::Init() fTofHitsMatches = InitOrFatal<TClonesArray>("TofHitMatch"); fPrimVertex = InitOrFatal<CbmVertex>("PrimaryVertex."); + //fTofTracks = InitOrFatal<TClonesArray>("TofTracks"); // TODO: check this and next lines + FairRootManager* fairRootMan = FairRootManager::Instance(); + fTofTracks = (TClonesArray*) fairRootMan->GetObject("TofTrack"); + if (fTofTracks == nullptr) { LOG(warning) << "LmvmTask::Init : no TofTrack array!"; } + InitHists(); fKFFitter.Init(); @@ -243,16 +317,14 @@ void LmvmTask::Exec(Option_t*) FillCands(); CalculateNofTopologyPairs("hNofTopoPairs_gamma", ELmvmSrc::Gamma); CalculateNofTopologyPairs("hNofTopoPairs_pi0", ELmvmSrc::Pi0); - DifferenceSignalAndBg(); - SignalAndBgReco(); - FillAccRecVsMomHist(); + AnalyseGlobalTracks(); + AnalyseCandidates(); fCandsTotal.insert(fCandsTotal.end(), fCands.begin(), fCands.end()); LOG(info) << "fCandsTotal.size = " << fCandsTotal.size(); //} } // Exec -//TODO: Move this functionality to RichUtil class void LmvmTask::FillRichRingNofHits() { fNofHitsInRingMap.clear(); @@ -268,6 +340,15 @@ void LmvmTask::FillRichRingNofHits() } } +double LmvmTask::MinvScale(double minv) +{ + double scale = -1.; + if (fParticle == "inmed") scale = LmvmUtils::GetMassScaleInmed(minv); + else if (fParticle == "qgp") scale = LmvmUtils::GetMassScaleQgp(minv); + else scale = 1.; + return scale; +} + void LmvmTask::FillMomHists(const CbmMCTrack* mct, const LmvmCand* cand, ELmvmSrc src, ELmvmAnaStep step) { @@ -277,15 +358,16 @@ void LmvmTask::FillMomHists(const CbmMCTrack* mct, const LmvmCand* cand, ELmvmSr } bool isMc = (mct != nullptr); string chargeStr = (isMc) ? LmvmUtils::GetChargeStr(mct) : LmvmUtils::GetChargeStr(cand); + double w = ((mct != nullptr && LmvmUtils::IsMcSignalEl(mct)) || (cand != nullptr && cand->IsMcSignal())) ? fW : 1.; for (const string& suff : {string(""), chargeStr}) { if (suff == "0") continue; - fH.FillH1("hMom" + suff, src, step, (isMc) ? mct->GetP() : cand->fMomentum.Mag(), fW); - fH.FillH1("hMomPx" + suff, src, step, (isMc) ? mct->GetPx() : cand->fMomentum.X(), fW); - fH.FillH1("hMomPy" + suff, src, step, (isMc) ? mct->GetPy() : cand->fMomentum.Y(), fW); - fH.FillH1("hMomPz" + suff, src, step, (isMc) ? mct->GetPz() : cand->fMomentum.Z(), fW); - fH.FillH1("hPt" + suff, src, step, (isMc) ? mct->GetPt() : cand->fMomentum.Perp(), fW); - fH.FillH1("hRapidity" + suff, src, step, (isMc) ? mct->GetRapidity() : cand->fRapidity, fW); + fH.FillH1("hMom" + suff, src, step, (isMc) ? mct->GetP() : cand->fMomentum.Mag(), w); + fH.FillH1("hMomPx" + suff, src, step, (isMc) ? mct->GetPx() : cand->fMomentum.X(), w); + fH.FillH1("hMomPy" + suff, src, step, (isMc) ? mct->GetPy() : cand->fMomentum.Y(), w); + fH.FillH1("hMomPz" + suff, src, step, (isMc) ? mct->GetPz() : cand->fMomentum.Z(), w); + fH.FillH1("hPt" + suff, src, step, (isMc) ? mct->GetPt() : cand->fMomentum.Perp(), w); + fH.FillH1("hRapidity" + suff, src, step, (isMc) ? mct->GetRapidity() : cand->fRapidity, w); } } @@ -301,15 +383,16 @@ void LmvmTask::DoMcTrack() string chargeStr = (mct->GetCharge() > 0) ? "+" : "-"; bool isAcc = IsMcTrackAccepted(i); double mom = mct->GetP(); + double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; FillMomHists(mct, nullptr, src, ELmvmAnaStep::Mc); if (isAcc) FillMomHists(mct, nullptr, src, ELmvmAnaStep::Acc); if (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4) - fH.FillH1("hMomAcc" + chargeStr + "_sts", src, mom, fW); - if (fNofHitsInRingMap[i] >= 7) fH.FillH1("hMomAcc" + chargeStr + "_rich", src, mom, fW); - 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); + fH.FillH1("hMomAcc" + chargeStr + "_sts", src, mom, w); + if (fNofHitsInRingMap[i] >= 7) fH.FillH1("hMomAcc" + chargeStr + "_rich", src, mom, w); + if (mct->GetNPoints(ECbmModuleId::kTrd) >= 2) fH.FillH1("hMomAcc" + chargeStr + "_trd", src, mom, w); + if (mct->GetNPoints(ECbmModuleId::kTof) >= 1) fH.FillH1("hMomAcc" + chargeStr + "_tof", src, mom, w); if (LmvmUtils::IsMcGammaEl(mct, fMCTracks)) { TVector3 v; @@ -336,17 +419,20 @@ void LmvmTask::DoMcTrack() if (isAcc) fH.FillH1("hMotherPdg_acc", mcMotherPdg); } } + + if (std::abs(mcPdg) == 211) fH.FillH1("hPionSupp_pi", ELmvmAnaStep::Mc, mom, w); + if (std::abs(mcPdg) == 211 && IsMcTrackAccepted(i)) fH.FillH1("hPionSupp_pi", ELmvmAnaStep::Acc, mom, w); } } void LmvmTask::DoMcPair() { int nMcTracks = fMCTracks->GetEntries(); - for (int iMc1 = 0; iMc1 < nMcTracks; iMc1++) { + for (int iMc1 = 0; iMc1 < nMcTracks; iMc1++) { // TODO: range until iMc1 < nMcTracks-1 ?? 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; bool isAcc1 = IsMcTrackAccepted(iMc1); @@ -362,14 +448,17 @@ void LmvmTask::DoMcPair() for (const auto step : {ELmvmAnaStep::Mc, ELmvmAnaStep::Acc}) { if (step == ELmvmAnaStep::Acc && !isAccPair) continue; - //fH.FillH1("hAnglePair", srcPair, step, pKin.fAngle, fW); + //fH.FillH1("hAnglePair", srcPair, step, pKin.fAngle, w); if (srcPair == ELmvmSrc::Signal) { fH.FillH2("hPtYPairSignal", step, pKin.fRapidity, pKin.fPt, fW); fH.FillH1("hMomPairSignal", step, pKin.fMomentumMag, fW); } + // MC and Acc minv only for signal, eta and pi0 if (srcPair == ELmvmSrc::Signal || srcPair == ELmvmSrc::Pi0 || srcPair == ELmvmSrc::Eta) { - fH.FillH1("hMinv", srcPair, step, pKin.fMinv, fW); + double w = 1; + if (srcPair == ELmvmSrc::Signal) w = fW * MinvScale(pKin.fMinv); + fH.FillH1("hMinv", srcPair, step, pKin.fMinv, w); } } } @@ -392,7 +481,8 @@ void LmvmTask::RichPmtXY() TVector3 v; mct->GetStartVertex(v); ELmvmSrc src = LmvmUtils::GetMcSrc(mct, fMCTracks); - if (v.Z() < 2.) { fH.FillH2("hPmtXY", src, richHit->GetX(), richHit->GetY(), fW); } + double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; + if (v.Z() < 2.) { fH.FillH2("hPmtXY", src, richHit->GetX(), richHit->GetY(), w); } } } @@ -405,80 +495,374 @@ bool LmvmTask::IsMcTrackAccepted(int mcTrackInd) && tr->GetNPoints(ECbmModuleId::kTrd) >= 2 && tr->GetNPoints(ECbmModuleId::kTof) > 1); } -void LmvmTask::FillAccRecVsMomHist() +void LmvmTask::AnalyseGlobalTracks() { - for (const int& pdg : {11, 211, 2212}) { // TODO: restructure this more efficiently (wihtout these pdg-loops) - int nofMcTracks = fMCTracks->GetEntriesFast(); - string hName = (pdg == 11) ? "hElMom" : (pdg == 211) ? "hPiMom" : "hProtonMom"; - - for (int i = 0; i < nofMcTracks; i++) { - CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i)); - bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4); - TVector3 vertex; - mct->GetStartVertex(vertex); - - if (std::abs(mct->GetPdgCode()) == pdg) { - fH.FillH1(hName + "_all_mc", mct->GetP()); - if (isAccSts) fH.FillH1(hName + "_all_acc", mct->GetP()); - - if (vertex.Mag() < 0.1) { - fH.FillH1(hName + "_prim_mc", mct->GetP()); - if (isAccSts) fH.FillH1(hName + "_prim_acc", mct->GetP()); - } + int nofMcTracks = fMCTracks->GetEntriesFast(); + + for (int i = 0; i < nofMcTracks; i++) { + CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i)); + int pdg = mct->GetPdgCode(); + bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4); + TVector3 vertex; + mct->GetStartVertex(vertex); + + if (std::abs(pdg) == 11 || std::abs(pdg) == 211 || pdg == 2212 || std::abs(pdg) == 321) { + string hName = (std::abs(pdg) == 11) ? "hEl" : (std::abs(pdg) == 211) ? "hPi" : pdg == 2212 ? "hProton" : "hKaon"; + fH.FillH1(hName + "Mom_all_mc", mct->GetP()); + if (isAccSts) fH.FillH1(hName + "Mom_all_acc", mct->GetP()); + + if (vertex.Mag() <= 0.1) { + fH.FillH1(hName + "Mom_prim_mc", mct->GetP()); + if (isAccSts) fH.FillH1(hName + "Mom_prim_acc", mct->GetP()); } } + } // MC tracks - int ngTracks = fGlobalTracks->GetEntriesFast(); - for (int i = 0; i < ngTracks; i++) { - CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(i)); - if (gTrack == nullptr) continue; - int stsInd = gTrack->GetStsTrackIndex(); - bool isRich = (gTrack->GetRichRingIndex() >= 0); - bool isTrd = (gTrack->GetTrdTrackIndex() >= 0); - bool isTof = (gTrack->GetTofHitIndex() >= 0); + int ngTracks = fGlobalTracks->GetEntriesFast(); + for (int i = 0; i < ngTracks; i++) { + CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(i)); + if (gTrack == nullptr) continue; + int stsInd = gTrack->GetStsTrackIndex(); + bool isRich = (gTrack->GetRichRingIndex() >= 0); + bool isTrd = (gTrack->GetTrdTrackIndex() >= 0); + bool isTof = (gTrack->GetTofHitIndex() >= 0); - if (stsInd < 0) continue; - CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd)); - if (stsTrack == nullptr) continue; - CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd)); - if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue; - int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex(); - if (stsMcTrackId < 0) continue; - CbmMCTrack* mct = (CbmMCTrack*) fMCTracks->At(stsMcTrackId); - if (mct == nullptr) continue; - TVector3 vertex; - mct->GetStartVertex(vertex); - - if (std::abs(mct->GetPdgCode()) == pdg) { - fH.FillH1(hName + "_all_recSts", mct->GetP()); - if (isRich) { fH.FillH1(hName + "_all_recStsRich", mct->GetP()); } - if (isRich && isTrd) { fH.FillH1(hName + "_all_recStsRichTrd", mct->GetP()); } - if (isRich && isTrd && isTof) { fH.FillH1(hName + "_all_recStsRichTrdTof", mct->GetP()); } - - if (vertex.Mag() < 0.1) { - fH.FillH1(hName + "_prim_recSts", mct->GetP()); - if (isRich) { fH.FillH1(hName + "_prim_recStsRich", mct->GetP()); } - if (isRich && isTrd) { fH.FillH1(hName + "_prim_recStsRichTrd", mct->GetP()); } - if (isRich && isTrd && isTof) { fH.FillH1(hName + "_prim_recStsRichTrdTof", mct->GetP()); } - } + if (stsInd < 0) continue; + CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd)); + if (stsTrack == nullptr) continue; + CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd)); + if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue; + int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex(); + if (stsMcTrackId < 0) continue; + CbmMCTrack* mct = (CbmMCTrack*) fMCTracks->At(stsMcTrackId); + if (mct == nullptr) continue; + + double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; + + int pdg = mct->GetPdgCode(); + double p = mct->GetP(); + double pt = mct->GetPt(); + double rap = mct->GetRapidity(); + double m2 = CbmLitGlobalElectronId::GetInstance().GetTofM2(i, p); + + TVector3 v; + mct->GetStartVertex(v); + bool isPrim = IsPrimary(v.Mag()); + if (std::abs(pdg) == 11 || std::abs(pdg) == 211 || pdg == 2212 || std::abs(pdg) == 321) { + string hName = (std::abs(pdg) == 11) ? "hEl" : (std::abs(pdg) == 211) ? "hPi" : pdg == 2212 ? "hProton" : "hKaon"; + fH.FillH1(hName + "Mom_all_recSts", p); + if (isRich) fH.FillH1(hName + "Mom_all_recStsRich", p); + if (isRich && isTrd) fH.FillH1(hName + "Mom_all_recStsRichTrd", p); + if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_all_recStsRichTrdTof", p); + + if (isPrim) { + fH.FillH1(hName + "Mom_prim_recSts", p); + if (isRich) fH.FillH1(hName + "Mom_prim_recStsRich", p); + if (isRich && isTrd) fH.FillH1(hName + "Mom_prim_recStsRichTrd", p); + if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_prim_recStsRichTrdTof", p); } } + + string ptcl = GetPidString(v.Mag(), pdg); + + fH.FillH2("hVertexGTrackRZ", v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y())); + /*fH.FillH2("hIndexStsRich_" + ptcl, stsInd, gTrack->GetRichRingIndex()); + fH.FillH2("hIndexStsTrd_" + ptcl, stsInd, gTrack->GetTrdTrackIndex()); + fH.FillH2("hIndexStsTof_" + ptcl, stsInd, gTrack->GetTofHitIndex());*/ + + bool isTofEl = (CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) ? true : false; + bool isElectron = (CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, p) && CbmLitGlobalElectronId::GetInstance().IsTrdElectron(i, p) && CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) ? true : false; + + if (gTrack != nullptr) { + CheckMismatches(gTrack, pdg, isElectron, ptcl, w); + } + if (std::abs(pdg) != 11) PidVsMom(gTrack, i, pdg, p); + + Double_t richDist = CbmRichUtil::GetRingTrackDistance(i); + fH.FillH2("hRichRingTrackDist_gTracks_" + ptcl, p, richDist); + + // investigate misidentifications: compare misidentified candidates with same not-misidentified particles + // first check if global track has entries in all detectors + if (!IsInAllDets(gTrack)) continue; // Mind: all following actions are done only for fully rec. gTracks + + CheckTofIdentification(gTrack, ptcl, p, m2, pdg, isTofEl); + + fH.FillH1("hMom_gTracks_" + ptcl, p); + fH.FillH2("hPtY_gTracks_" + ptcl, rap, pt); + fH.FillH2("hTofM2_gTracks_" + ptcl, p, m2); + + // check PIDs in "Tof pile" + if (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01) { + bool isSignal = (mct != nullptr && mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11); + double pdgBin = (pdg == 11 && isSignal) ? 0.5 : (pdg == -11 && isSignal) ? 1.5 + : (pdg == 11 && !isSignal && isPrim) ? 2.5 : (pdg == -11 && !isSignal && isPrim) ? 3.5 + : (pdg == 11 && !isSignal && !isPrim) ? 4.5 : (pdg == -11 && !isSignal && !isPrim) ? 5.5 + : (pdg == 211) ? 6.5 : (pdg == -211) ? 7.5 : (pdg == 2212) ? 8.5 : (pdg == 321) ? 9.5 : (pdg == -321) ? 10.5 : 11.5; + fH.FillH1("hTofPilePdgs_gTracks", pdgBin); + } + + string ptcl2 = GetPidString(mct, nullptr); + if (!isElectron && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) { + fH.FillH1("hMom_" + ptcl2 + "_true", p); + fH.FillH2("hPtY_" + ptcl2 + "_true", rap, pt); + fH.FillH2("hTofM2_" + ptcl2 + "_true", p, m2); + } + else if (isElectron && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) { + fH.FillH1("hMom_" + ptcl2 + "_misid", p); + fH.FillH2("hPtY_" + ptcl2 + "_misid", rap, pt); + fH.FillH2("hTofM2_" + ptcl2 + "_misid", p, m2); + } + + BetaMom(mct, gTrack, ptcl2); } +} - // Fill histo with misidentified pions, protons, kaons // TODO: move this to SignalAndBgReco()? - for (const auto& cand : fCands) { - int pdg = cand.fMcPdg; - if (!(pdg == 211 || pdg == 111 || pdg == 2212 || pdg == 321)) continue; - if (!cand.fIsElectron) continue; - cout << "LmvmTask::FillAccRecVsMomHist(): pdg = " << pdg << endl; - string hName = (pdg == 211) ? "hCandMisIdAsEl_pi" - : (pdg == 111) ? "hCandMisIdAsEl_pi0" - : (pdg == 2212) ? "hCandMisIdAsEl_proton" - : "hCandMisIdAsEl_kaon"; - fH.FillH1(hName, cand.fMomentum.Mag()); +bool LmvmTask::IsInTofPile(double p, double m2) +{ + return (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01); +} + +void LmvmTask::PidVsMom(const CbmGlobalTrack* gTrack, int iGTrack, int pdg, double p) +{ + bool isRichEl = CbmLitGlobalElectronId::GetInstance().IsRichElectron(iGTrack, p); + bool isTrdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iGTrack, p); + bool isTofEl = CbmLitGlobalElectronId::GetInstance().IsTofElectron(iGTrack, p); + + double pdgBin = (pdg == 211) ? 0.5 : (pdg == -211) ? 1.5 : (pdg == 2212) ? 2.5 : (pdg == 321) ? 3.5 : (pdg == -321) ? 4.5 : 5.5; + + if (isRichEl) fH.FillH2("hPdgVsMom_gTracks_rich_all", p, pdgBin); + if (isTrdEl) fH.FillH2("hPdgVsMom_gTracks_trd_all", p, pdgBin); + if (isTofEl) fH.FillH2("hPdgVsMom_gTracks_tof_all", p, pdgBin); + + if (isRichEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_rich_complete", p, pdgBin); + if (isTrdEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_trd_complete", p, pdgBin); + if (isTofEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_tof_complete", p, pdgBin); +} + +void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string& gtString, double pMc, double m2, int pdg, bool isTofEl) +{ + // Get STS ID + int stsInd = gTrack->GetStsTrackIndex(); + if (stsInd < 0) return; + CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd)); + if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) return; + int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex(); + if (stsMcTrackId < 0) return; + + // Get ToF Match + Int_t tofInd = gTrack->GetTofHitIndex(); + if (tofInd < 0) return; + CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd); + if (tofHit == nullptr) return; + CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(gTrack->GetStsTrackIndex())); + if (tofTrack == nullptr) return; + CbmMatch* tofHitMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofInd)); + if (tofHitMatch == nullptr) return; + int tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex(); + if (tofPointIndex < 0) return; + FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofPointIndex)); + if (tofPoint == nullptr) return; + int tofMcTrackId = tofPoint->GetTrackID(); + + bool isTofMismatch = (tofMcTrackId == stsMcTrackId) ? false : true; + bool isTofMisid = (std::abs(pdg) == 11 && isTofEl) ? false : (std::abs(pdg) != 11 && !isTofEl) ? false : true; + + double pointX = tofPoint->GetX(); + double pointY = tofPoint->GetY(); + double hitX = tofHit->GetX(); + double hitY = tofHit->GetY(); + double dX = pointX - hitX; + double dY = pointY - hitY; + double dR = sqrt(dX * dX + dY * dY); + double dist = tofTrack->GetDistance(); + + string hNameId = (!isTofMisid) ? "_trueid_" + gtString : "_misid_" + gtString; + fH.FillH2("hTofPointXY" + hNameId, pointX, pointY); + fH.FillH2("hTofHitXY" + hNameId, hitX, hitY); + fH.FillH1("hTofHitPointDist" + hNameId, dR); + fH.FillH2("hTofHitTrackDist_gTracks" + hNameId, pMc, dist); + + string hNameMatch = (!isTofMismatch) ? "_truematch_" + gtString : "_mismatch_" + gtString; + fH.FillH2("hTofPointXY" + hNameMatch, pointX, pointY); + fH.FillH2("hTofHitXY" + hNameMatch, hitX, hitY); + fH.FillH1("hTofHitPointDist" + hNameMatch, dR); + fH.FillH2("hTofHitTrackDist_gTracks" + hNameMatch, pMc, dist); + + // check m2 calculation in ToF + double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime(); + Double_t noOffsetTime = tofHit->GetTime() - eventTime; + Double_t t = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m + + fH.FillH2("hTofTimeVsMom_gTracks" + hNameId, pMc, t); + fH.FillH2("hTofTimeVsMom_gTracks" + hNameMatch, pMc, t); + + if (!isTofMisid) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 0.5); + if (isTofMisid) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 1.5); + if (!isTofMismatch) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 2.5); + if (isTofMismatch) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 3.5); + + // check ToF data for particles in ToF pile // TODO: needed or can be deleted? + if (isTofEl && std::abs(pdg) != 11 && IsInTofPile(pMc, m2)) { + string candString = (pdg == 211) ? fH.fCandNames[4] : (pdg == -211) ? fH.fCandNames[5] : (pdg == 2212) ? fH.fCandNames[6] : (pdg == 321) ? fH.fCandNames[7] : (pdg == -321) ? fH.fCandNames[8] : fH.fCandNames[9]; + fH.FillH2("hTofPileHitXY_" + candString, hitX, hitY); + fH.FillH2("hTofPilePointXY_" + candString, pointX, pointY); + fH.FillH1("hTofPileHitPointDist_" + candString, dR); + } +} + +void LmvmTask::BetaMom(const CbmMCTrack* mct, const CbmGlobalTrack* gTrack, const string& ptcl) +{ + Int_t tofInd = gTrack->GetTofHitIndex(); + if (tofInd < 0) return; + CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd); + if (NULL == tofHit) return; + double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime(); + Double_t noOffsetTime = tofHit->GetTime() - eventTime; + + double p = mct->GetP(); + double l = gTrack->GetLength(); + double t = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m + double q = (mct->GetCharge() > 0) ? 1 : (mct->GetCharge() < 0) ? -1. : 0.; + double beta = l / (t * 100); // cm to m + fH.FillH2("hBetaMom_gTracks_" + ptcl, q * p, beta); + fH.FillH2("hBetaMom_allGTracks", q * p, beta); +} + +void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isElectron, const string& ptcl, double w) +{ + // Chi2 of true-/mismatched; number of mismatches (general and seperated for detectors) + int stsInd = gTrack->GetStsTrackIndex(); + if (stsInd < 0) return; + CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd)); + if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) return; + int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex(); + if (stsMcTrackId < 0) return; + + bool isFull = IsInAllDets(gTrack); + + // first calculate Chi2 + CbmL1PFFitter fPFFitter; + CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd)); + vector<CbmStsTrack> stsTracks; + stsTracks.resize(1); + stsTracks[0] = *stsTrack; + vector<CbmL1PFFitter::PFFieldRegion> vField; + vector<float> chiPrim; + fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6); + double chi2 = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF(); + + CbmTrackMatchNew* richMatch = nullptr; + CbmTrackMatchNew* trdMatch = nullptr; + CbmMatch* tofMatch = nullptr; + + int richRingInd = gTrack->GetRichRingIndex(); + int trdTrackInd = gTrack->GetTrdTrackIndex(); + int tofHitInd = gTrack->GetTofHitIndex(); + + if (richRingInd >= 0) richMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(richRingInd)); + if (trdTrackInd >= 0) trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(trdTrackInd)); + if (tofHitInd >= 0) tofMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofHitInd)); + + int nofMismTrackSegmentsAll = 0; + int nofMismTrackSegmentsComp = 0; + + fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 0.5, w); + if (isFull) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 0.5, w); + + if (richMatch != nullptr) { + int richMcTrackId = richMatch->GetMatchedLink().GetIndex(); + if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 1.5, w); + if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_rich_" + ptcl, chi2, w); + if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) { + fH.FillH1("hChi2_mismatch_all_rich_" + ptcl, chi2, w); + fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 2.5, w); + nofMismTrackSegmentsAll++; + } + if (isFull) { + if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 1.5, w); + if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_complete_rich_" + ptcl, chi2, w); + if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) { + fH.FillH1("hChi2_mismatch_complete_rich_" + ptcl, chi2, w); + fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 2.5, w); + nofMismTrackSegmentsComp++; + } + } + } + + if (trdMatch != nullptr) { + int trdMcTrackId = trdMatch->GetMatchedLink().GetIndex(); + if (trdMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 3.5, w); + if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_trd_" + ptcl, chi2, w); + if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) { + fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 4.5, w); + fH.FillH1("hChi2_mismatch_all_trd_" + ptcl, chi2, w); + nofMismTrackSegmentsAll++; + } + if (isFull) { + if (trdMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 3.5, w); + if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_complete_trd_" + ptcl, chi2, w); + if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) { + fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 4.5, w); + fH.FillH1("hChi2_mismatch_complete_trd_" + ptcl, chi2, w); + nofMismTrackSegmentsComp++; + } + } + } + + if (tofMatch != nullptr && tofMatch->GetMatchedLink().GetIndex() >= 0) { + FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofMatch->GetMatchedLink().GetIndex())); + if (tofPoint != nullptr) { + int tofMcTrackId = tofPoint->GetTrackID(); + if (tofMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 5.5, w); + if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_tof_" + ptcl, chi2, w); + if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) { + fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 6.5, w); + fH.FillH1("hChi2_mismatch_all_tof_" + ptcl, chi2, w); + nofMismTrackSegmentsAll++; + } + if (isFull) { + if (tofMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 5.5, w); + if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_complete_tof_" + ptcl, chi2, w); + if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) { + fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 6.5, w); + fH.FillH1("hChi2_mismatch_complete_tof_" + ptcl, chi2, w); + nofMismTrackSegmentsComp++; + } + } + } } + + fH.FillH1("hNofMismatchedTrackSegments_all_" + ptcl, nofMismTrackSegmentsAll + 0.5, w); + if (isFull) fH.FillH1("hNofMismatchedTrackSegments_complete_" + ptcl, nofMismTrackSegmentsComp + 0.5, w); + + bool isTrueId = (std::abs(pdg) == 11 && isElectron) ? true : (std::abs(pdg) != 11 && !isElectron) ? true : false; + double binX = nofMismTrackSegmentsComp + 0.5; + double binY = (isTrueId == true) ? 0.5 : 1.5; + + fH.FillH2("hMatchId_gTracks_" + ptcl, binX, binY, w); } +bool LmvmTask::IsInAllDets(const CbmGlobalTrack* gTrack) +{ + if (gTrack == nullptr) return false; + + int stsInd = gTrack->GetStsTrackIndex(); + if (stsInd < 0) return false; + CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd)); + if (stsTrack == nullptr) return false; + + int richInd = gTrack->GetRichRingIndex(); + int trdInd = gTrack->GetTrdTrackIndex(); + int tofInd = gTrack->GetTofHitIndex(); + if (richInd < 0 || trdInd < 0 || tofInd < 0) return false; + + CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(richInd)); + CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(trdInd)); + CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(tofInd)); + if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr) return false; + + return true; +} void LmvmTask::FillTopologyCands() { @@ -497,9 +881,10 @@ void LmvmTask::FillTopologyCands() CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd)); if (stsTrack == nullptr) continue; - cand.fRichInd = gTrack->GetRichRingIndex(); - cand.fTrdInd = gTrack->GetTrdTrackIndex(); - cand.fTofInd = gTrack->GetTofHitIndex(); + cand.fRichInd = gTrack->GetRichRingIndex(); + cand.fTrdInd = gTrack->GetTrdTrackIndex(); + cand.fTofInd = gTrack->GetTofHitIndex(); // TODO: is TofHitIndex; change this name everywhere + cand.fTofTrackInd = gTrack->GetStsTrackIndex(); // TODO: is this right?? LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex); cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim); @@ -569,6 +954,24 @@ void LmvmTask::FillCands() CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd)); if (stsTrack == nullptr) continue; + // RICH - TRD - TOF + cand.fRichInd = gTrack->GetRichRingIndex(); + cand.fTrdInd = gTrack->GetTrdTrackIndex(); + cand.fTofInd = gTrack->GetTofHitIndex(); + cand.fTofTrackInd = gTrack->GetStsTrackIndex(); + + // Set length and time + cand.fLength = gTrack->GetLength() / 100.; + Int_t tofInd = gTrack->GetTofHitIndex(); + if (tofInd >= 0) { + CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd); + if (NULL != tofHit) { + double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime(); + Double_t noOffsetTime = tofHit->GetTime() - eventTime; + cand.fTime = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m + } + } + LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex); cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim); cand.fIsPtCut = fCuts.IsPtCutOk(cand.fMomentum.Perp()); @@ -593,21 +996,29 @@ void LmvmTask::FillCands() } } - // RICH - TRD - TOF - cand.fRichInd = gTrack->GetRichRingIndex(); - cand.fTrdInd = gTrack->GetTrdTrackIndex(); - cand.fTofInd = gTrack->GetTofHitIndex(); - if (cand.fRichInd < 0 || cand.fTrdInd < 0 || cand.fTofInd < 0) continue; + if (cand.fRichInd < 0 || cand.fTrdInd < 0 || cand.fTofInd < 0 || cand.fTofTrackInd < 0) continue; + + cand.fGTrackInd = iGTrack; CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(cand.fRichInd)); CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(cand.fTrdInd)); - CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd)); - if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr) continue; + CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd)); + CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(cand.fTofTrackInd)); + if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr || tofTrack == nullptr) continue; + + cand.fChi2Rich = (richRing->GetNDF() > 0.) ? richRing->GetChi2() / richRing->GetNDF() : 0.; + cand.fChi2Trd = (trdTrack->GetNDF() > 0.) ? trdTrack->GetChiSq() / trdTrack->GetNDF() : 0.; - LmvmUtils::IsElectron(iGTrack, cand.fMomentum.Mag(), fCuts.fMomentumCut, &cand); cand.fTrdLikeEl = trdTrack->GetPidLikeEL(); cand.fTrdLikePi = trdTrack->GetPidLikePI(); + cand.fTofDist = tofTrack->GetDistance(); + + LmvmUtils::IsElectron(iGTrack, cand.fMomentum.Mag(), fCuts.fMomentumCut, &cand); + LmvmUtils::IsRichElectron(iGTrack, cand.fMomentum.Mag(), &cand); + LmvmUtils::IsTrdElectron(iGTrack, cand.fMomentum.Mag(), &cand); + LmvmUtils::IsTofElectron(iGTrack, cand.fMomentum.Mag(), &cand); + fCands.push_back(cand); if (!cand.fIsElectron && cand.fIsChi2Prim) fTTCands.push_back(cand); @@ -624,76 +1035,75 @@ void LmvmTask::CombinatorialPairs() size_t nCand = fCandsTotal.size(); for (size_t iC1 = 0; iC1 < nCand - 1; iC1++) { const auto& cand1 = fCandsTotal[iC1]; - + for (size_t iC2 = iC1 + 1; iC2 < nCand; iC2++) { const auto& cand2 = fCandsTotal[iC2]; LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2); - double w = (cand1.IsMcSignal() && cand2.IsMcSignal()) ? fW : 1.; + double w = 1.; + if (cand1.IsMcSignal()) w = fW * MinvScale(cand1.fMassSig); + else if (cand2.IsMcSignal()) w = fW * MinvScale(cand2.fMassSig); bool isSameEvent = (cand1.fEventNumber == cand2.fEventNumber); + // PLUTO electrons from two different events have to be scaled with w^2! + if (!isSameEvent && cand1.IsMcSignal() && cand2.IsMcSignal()) w = fW * fW * MinvScale(cand1.fMassSig) * MinvScale(cand2.fMassSig); 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) { - // only PLUTO electrons - if (cand1.IsMcSignal() && cand2.IsMcSignal()) { + if (cand1.IsCutTill(step) && cand2.IsCutTill(step)) { + + // all particles + if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue; + if (cand1.fCharge * cand2.fCharge < 0) { 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); + fH.FillH1("hMinvCombPM_sameEv", step, pRec.fMinv, w); + cout << "LmvmTask::CombinatorialPairs(): e+e- same event (all). w = " << w << endl; //TODO: delete this line } - 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 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); + cout << "LmvmTask::CombinatorialPairs(): e-e- same event (all). w = " << w << endl; //TODO: delete this line } + else fH.FillH1("hMinvCombMM_mixedEv", step, pRec.fMinv, w); } - // only UrQMD electrons - else if (!cand1.IsMcSignal() && !cand2.IsMcSignal()) { + if (cand1.fCharge > 0 && cand2.fCharge > 0) { 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); + fH.FillH1("hMinvCombPP_sameEv", step, pRec.fMinv, w); + cout << "LmvmTask::CombinatorialPairs(): e+e+ same event (all). w = " << w << endl; //TODO: delete this line + } + else fH.FillH1("hMinvCombPP_mixedEv", step, pRec.fMinv, w); + } + + // only UrQMD particles (also misidentified) + if (!cand1.IsMcSignal() && !cand2.IsMcSignal()) { + if (isSameEvent) { + if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdAll_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdAll_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdAll_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); + if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdAll_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdAll_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdAll_mixedEv", step, pRec.fMinv, w); } } - } - // for common CB, don't use PLUTO signals - 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); - } - 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); - } - 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); + // only UrQMD electrons + if (!cand1.IsMcSignal() && !cand2.IsMcSignal() && std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11) { // TODO: delete double brackets + if (isSameEvent) { + if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdEl_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdEl_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdEl_sameEv", step, pRec.fMinv, w); + } + else { + if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdEl_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdEl_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdEl_mixedEv", step, pRec.fMinv, w); + } } - } - } // steps - } // cand2 - } // cand1 + } // isCutTillStep + } // steps + } // cand2 + } // cand1 } void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands) @@ -712,6 +1122,18 @@ void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands) cand.fMcPdg = mct->GetPdgCode(); cand.fMcSrc = LmvmUtils::GetMcSrc(mct, fMCTracks); + // Get mass of mother particle if it is signal (for mass-dependant scaling) + if (cand.IsMcSignal()) { + int nMcTracks = fMCTracks->GetEntries(); + for (int iMc2 = 0; iMc2 < nMcTracks; iMc2++) { + CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2)); + if (mct2->GetMotherId() == cand.fMcMotherId && iMc2 != cand.fStsMcTrackId) { + LmvmKinePar pKin = LmvmKinePar::Create(mct, mct2); + cand.fMassSig = pKin.fMinv; + } + } + } + if (std::abs(cand.fMcPdg) == 211 && fPionMisidLevel >= 0.) continue; // RICH @@ -759,7 +1181,8 @@ void LmvmTask::AssignMcToTopologyCands(vector<LmvmCand>& topoCands) void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step, const LmvmKinePar& parRec) { ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM); - fH.FillH1("hAnglePair", src, step, parRec.fAngle, fW); + double w = (candP.IsMcSignal() || candM.IsMcSignal()) ? fW : 1.; + fH.FillH1("hAnglePair", src, step, parRec.fAngle, w); if (src == ELmvmSrc::Bg) { // gamma=0.5, pi0=1.5, pions=2.5, other=3.5 @@ -768,73 +1191,73 @@ void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAna fH.FillH2("hSrcBgPairsEpEm", step, indP, indM); 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 (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); + } - if (std::abs(candP.fMcPdg) == 11) { + // 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 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 but not from Gamma or Pi0 + else { + if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "go", parRec.fMinv); + else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0()) + fH.FillH1(hName + "pi0o", parRec.fMinv); + else if (std::abs(candM.fMcPdg) == 211) + fH.FillH1(hName + "pio", parRec.fMinv); + else + fH.FillH1(hName + "oo", parRec.fMinv); + } } - // cand1 is El and from Pi0 - else if (candP.IsMcPi0()) { - if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gpi0", 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 + "pi0pi0", parRec.fMinv); - else if (std::abs(candM.fMcPdg) == 211) fH.FillH1(hName + "pipi0", parRec.fMinv); + else if (std::abs(candM.fMcPdg) == 211) + fH.FillH1(hName + "pipi", parRec.fMinv); else - fH.FillH1(hName + "pi0o", parRec.fMinv); + fH.FillH1(hName + "pio", parRec.fMinv); } - // cand1 is El but not from Gamma or Pi0 + // 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 + "pio", parRec.fMinv); + fH.FillH1(hName + "pipi0", 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); - } } - //} } } @@ -933,8 +1356,17 @@ void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const bool isMismatch = (LmvmUtils::IsMismatch(candP) || LmvmUtils::IsMismatch(candM)); ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM); - fH.FillH1("hMinv", src, step, parRec.fMinv, fW); - fH.FillH2("hMinvPt", src, step, parRec.fMinv, parMc.fPt, fW); + double w = 1.; + if (candP.IsMcSignal()) w = fW * MinvScale(candP.fMassSig); + if (candM.IsMcSignal()) w = fW * MinvScale(candM.fMassSig); + if (w < 0) LOG(warning) << "LmvmTask::FillPairHists(): Signal mass < 0!"; + + fH.FillH1("hMinv", src, step, parRec.fMinv, w); + if (!candP.IsMcSignal() && !candM.IsMcSignal()) fH.FillH1("hMinv_urqmdAll", src, step, parRec.fMinv, w); + if (!candP.IsMcSignal() && !candM.IsMcSignal() && std::abs(candP.fMcPdg) == 11 && std::abs(candM.fMcPdg) == 11) + fH.FillH1("hMinv_urqmdEl", src, step, parRec.fMinv, w); + + fH.FillH2("hMinvPt", src, step, parRec.fMinv, parMc.fPt, w); PairSource(candP, candM, step, parRec); @@ -955,7 +1387,58 @@ void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const } } -void LmvmTask::SignalAndBgReco() +string LmvmTask::GetPidString(const CbmMCTrack* mct, const LmvmCand* cand) +{ + if ((mct != nullptr && cand != nullptr) || (mct == nullptr && cand == nullptr)) { + LOG(error) << "LmvmTask::GetPidString: Both mct and cand are [not nullptr] or [nullptr]."; + return 0; + } + + int pdg = (mct != nullptr) ? mct->GetPdgCode() : cand->fMcPdg; + bool isMcSignal = (mct != nullptr) ? (mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11) : cand->IsMcSignal(); + + string pidString = fH.fCandNames[9]; + if (isMcSignal && pdg == -11) pidString = fH.fCandNames[0]; + else if (isMcSignal && pdg == 11) pidString = fH.fCandNames[1]; + else if (!isMcSignal && pdg == -11) pidString = fH.fCandNames[2]; + else if (!isMcSignal && pdg == 11) pidString = fH.fCandNames[3]; + else if (pdg == 211) pidString = fH.fCandNames[4]; + else if (pdg == -211) pidString = fH.fCandNames[5]; + else if (pdg == 2212) pidString = fH.fCandNames[6]; + else if (pdg == 321) pidString = fH.fCandNames[7]; + else if (pdg == -321) pidString = fH.fCandNames[8]; + + return pidString; +} + +string LmvmTask::GetPidString(double vertexMag, int pdg) +{ + string pidString = fH.fGTrackNames[14]; + bool isPrim = IsPrimary(vertexMag); + + if (isPrim) { + if (pdg == -11) pidString = fH.fGTrackNames[1]; + else if (pdg == 11) pidString = fH.fGTrackNames[3]; + else if (pdg == 211) pidString = fH.fGTrackNames[5]; + else if (pdg == -211) pidString = fH.fGTrackNames[7]; + else if (pdg == 2212) pidString = fH.fGTrackNames[9]; + else if (pdg == 321) pidString = fH.fGTrackNames[11]; + else if (pdg == -321) pidString = fH.fGTrackNames[13]; + } + else { + if (pdg == -11) pidString = fH.fGTrackNames[0]; + else if (pdg == 11) pidString = fH.fGTrackNames[2]; + else if (pdg == 211) pidString = fH.fGTrackNames[4]; + else if (pdg == -211) pidString = fH.fGTrackNames[6]; + else if (pdg == 2212) pidString = fH.fGTrackNames[8]; + else if (pdg == 321) pidString = fH.fGTrackNames[10]; + else if (pdg == -321) pidString = fH.fGTrackNames[12]; + } + + return pidString; +} + +void LmvmTask::AnalyseCandidates() { CheckGammaConvAndPi0(); CheckTopologyCut(ELmvmTopologyCut::ST, "hStCut"); @@ -966,20 +1449,58 @@ void LmvmTask::SignalAndBgReco() CheckClosestMvdHit(2, "hMvdCut_2", "hMvdCutQa_2"); } + // single candidates for (const auto& cand : fCands) { - int pdg = 0; - if (cand.fStsMcTrackId > 0) { - CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId)); - if (mcTrack != nullptr) pdg = mcTrack->GetPdgCode(); - } + double w = (cand.IsMcSignal()) ? fW : 1.; + CbmMCTrack* mcTrack = nullptr; + if (cand.fStsMcTrackId >= 0) mcTrack = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId)); + int pdg = mcTrack->GetPdgCode(); + double mom = mcTrack->GetP(); + string pidString = GetPidString(nullptr, &cand); + for (auto step : fH.fAnaSteps) { if (cand.IsCutTill(step)) { TrackSource(cand, step, pdg); - fH.FillH2("hPtYCandidate", step, cand.fRapidity, cand.fMomentum.Perp(), fW); + FillCandPidValues(mcTrack, cand, step); // if std::abs(pdg) != 11, this is misidentification + if (mcTrack != nullptr) CheckTofId(mcTrack, cand, step, pdg); + + if (step >= ELmvmAnaStep::ElId && std::abs(pdg) == 211) fH.FillH1("hPionSupp_idEl", step, mom, w); + + //double richDist = CbmRichUtil::GetRingTrackDistance(cand.fGTrackInd); + //fH.FillH2("hRichRingTrackDist_cands_" + pidString, step, mom, richDist, w); // TODO: causes 'Error in <TClonesArray::At>: index xxx out of bounds' error } } - } + if (cand.IsCutTill(ELmvmAnaStep::ElId)) { + fH.FillH2("hChi2VsMom_sts_" + pidString, cand.fMomentum.Mag(), cand.fChi2Sts, w); + fH.FillH2("hChi2VsMom_rich_" + pidString, cand.fMomentum.Mag(), cand.fChi2Rich, w); + fH.FillH2("hChi2VsMom_trd_" + pidString, cand.fMomentum.Mag(), cand.fChi2Trd, w); + + fH.FillH2("hTofTimeVsChi2_sts_" + pidString, cand.fChi2Sts, cand.fTime, w); + fH.FillH2("hTofTimeVsChi2_rich_" + pidString, cand.fChi2Rich, cand.fTime, w); + fH.FillH2("hTofTimeVsChi2_trd_" + pidString, cand.fChi2Trd, cand.fTime, w); + + fH.FillH2("hChi2Comb_StsRich_" + pidString, cand.fChi2Sts, cand.fChi2Rich, w); + fH.FillH2("hChi2Comb_StsTrd_" + pidString, cand.fChi2Sts, cand.fChi2Trd, w); + fH.FillH2("hChi2Comb_RichTrd_" + pidString, cand.fChi2Rich, cand.fChi2Trd, w); + + fH.FillH2("hTofTimeVsMom_cands_" + pidString, cand.fMomentum.Mag(), cand.fTime, w); + fH.FillH2("hTofHitTrackDist_cands_" + pidString, cand.fMomentum.Mag(), cand.fTofDist, w); + } + + DifferenceSignalAndBg(cand); + + // beta-momentum spectrum + string ptcl = GetPidString(nullptr, &cand); + double p = cand.fMomentum.Mag(); + double m = cand.fMass; + double r = p / m; + double beta = r / sqrt(1. + r * r); + double q = (cand.fCharge > 0) ? 1 : (cand.fCharge < 0) ? -1. : 0.; + fH.FillH2("hBetaMom_cands_" + ptcl, q * p, beta); + } // single candidates + + // candidate pairs for (const auto& candP : fCands) { if (candP.fCharge < 0) continue; CbmMCTrack* mctrackP = @@ -1000,6 +1521,50 @@ void LmvmTask::SignalAndBgReco() } } +void LmvmTask::CheckTofId(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step, int pdg) +{ + TVector3 v; + mct->GetStartVertex(v); + bool isPrim = IsPrimary(v.Mag()); + double pt = mct->GetPt(); + double rap = mct->GetRapidity(); + + // check PIDs in "Tof pile" + if (cand.fMomentum.Mag() > 0.3 && cand.fMomentum.Mag() < 1. && cand.fMass2 > -0.012 && cand.fMass2 < 0.01) { + string pidString = GetPidString(nullptr, &cand); + double pdgBin = (pdg == 11 && cand.IsMcSignal()) ? 0.5 : (pdg == -11 && cand.IsMcSignal()) ? 1.5 + : (pdg == 11 && !cand.IsMcSignal() && isPrim) ? 2.5 : (pdg == -11 && !cand.IsMcSignal() && isPrim) ? 3.5 + : (pdg == 11 && !cand.IsMcSignal() && !isPrim) ? 4.5 : (pdg == -11 && !cand.IsMcSignal() && !isPrim) ? 5.5 + : (pdg == 211) ? 6.5 : (pdg == -211) ? 7.5 : (pdg == 2212) ? 8.5 : (pdg == 321) ? 9.5 : (pdg == -321) ? 10.5 : 11.5; + fH.FillH1("hTofPilePdgs_cands", step, pdgBin); + if (step == ELmvmAnaStep::ElId && cand.IsCutTill(step)) fH.FillH2("hTofPilePty_cands_" + pidString, rap, pt); // TODO: check newly (25.10.22) added '&& IsTill(step)'!! + } + + // check vertex of misidentified particles in ToF after electron ID // TODO: split this up into single contributions? + if (std::abs(pdg) != 11 && cand.fIsTofElectron) { + fH.FillH2("hVertexXZ_misidTof", step, v.Z(), v.X()); + fH.FillH2("hVertexYZ_misidTof", step, v.Z(), v.Y()); + fH.FillH2("hVertexXY_misidTof", step, v.X(), v.Y()); + fH.FillH2("hVertexRZ_misidTof", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y())); + } +} + +void LmvmTask::FillCandPidValues(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step) { + double w = (cand.IsMcSignal()) ? fW : 1.; + string pidString = GetPidString(nullptr, &cand); + + fH.FillH1("hMom_cands_" + pidString, step, cand.fMomentum.Mag(), w); + fH.FillH2("hPtY_cands_" + pidString, step, cand.fRapidity, cand.fMomentum.Perp(), w); + fH.FillH2("hTofM2_cands_" + pidString, step, cand.fMomentum.Mag(), cand.fMass2, w); + + if (mct != nullptr) { + double rat = cand.fMomentum.Mag() / mct->GetP(); + string ptcl = GetPidString(nullptr, &cand); + fH.FillH1("hMomRatio_cands_" + pidString, step, rat); + fH.FillH2("hMomRatioVsMom_cands_" + ptcl, step, mct->GetP(), rat); + } +} + void LmvmTask::CheckGammaConvAndPi0() { for (auto& candP : fCands) { @@ -1145,63 +1710,58 @@ void LmvmTask::CalculateNofTopologyPairs(const string& name, ELmvmSrc src) } } -void LmvmTask::DifferenceSignalAndBg() +void LmvmTask::DifferenceSignalAndBg(const LmvmCand& cand) { - for (const auto& cand : fCands) { - fH.FillH1("hChi2PrimVertex", cand.fMcSrc, cand.fChi2Prim, fW); + fH.FillH1("hChi2PrimVertex", cand.fMcSrc, cand.fChi2Prim, fW); + + if (!cand.fIsChi2Prim) return; + fH.FillH1("hAnnRich", cand.fMcSrc, cand.fRichAnn, fW); + fH.FillH2("hAnnRichVsMom", cand.fMcSrc, cand.fMomentum.Mag(), cand.fRichAnn, fW); + //fH.FillH1("hAnnTrd", cand.fMcSrc, cand.fTrdAnn, fW); // TODO: uncomment when TRD ANN is working (CbmLitGlobalElectronId::GetTrdAnn() gives back El-Likelihood) + fH.FillH2("hTrdLike_El", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikeEl, fW); + fH.FillH2("hTrdLike_Pi", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikePi, fW); + fH.FillH2("hTofM2", cand.fMcSrc, cand.fMomentum.Mag(), cand.fMass2, fW); + + // electron purity + if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) == 11) { + fH.FillH2("hAnnRichVsMomPur_El", cand.fMomentum.Mag(), cand.fRichAnn, fW); + fH.FillH2("hTrdElLikePur_El", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW); + } + else if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) != 11) { + fH.FillH2("hAnnRichVsMomPur_Bg", cand.fMomentum.Mag(), cand.fRichAnn, fW); + fH.FillH2("hTrdElLikePur_Bg", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW); + } - if (!cand.fIsChi2Prim) continue; - fH.FillH1("hAnnRich", cand.fMcSrc, cand.fRichAnn, fW); - fH.FillH2("hAnnRichVsMom", cand.fMcSrc, cand.fMomentum.Mag(), cand.fRichAnn, fW); - //fH.FillH1("hAnnTrd", cand.fMcSrc, cand.fTrdAnn, fW); // TODO: uncomment when TRD ANN is working (CbmLitGlobalElectronId::GetTrdAnn() gives back El-Likelihood) - fH.FillH2("hTrdLike_El", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikeEl, fW); - fH.FillH2("hTrdLike_Pi", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikePi, fW); - fH.FillH2("hTofM2", cand.fMcSrc, cand.fMomentum.Mag(), cand.fMass2, fW); - - // electron purity - if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) == 11) { - fH.FillH2("hAnnRichVsMomPur_El", cand.fMomentum.Mag(), cand.fRichAnn, fW); - fH.FillH2("hTrdElLikePur_El", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW); - fH.FillH2("hTofM2Pur_El", cand.fMomentum.Mag(), cand.fMass2, fW); - } - else if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) != 11) { - fH.FillH2("hAnnRichVsMomPur_Bg", cand.fMomentum.Mag(), cand.fRichAnn, fW); - fH.FillH2("hTrdElLikePur_Bg", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW); - fH.FillH2("hTofM2Pur_Bg", cand.fMomentum.Mag(), cand.fMass2, fW); - } + if (!cand.IsCutTill(ELmvmAnaStep::ElId)) return; + //fH.FillSourceH1("hPt", cand.fMcSrc, cand.fMomentum.Perp(), fW); + //fH.FillSourceH1("hMom", cand.fMcSrc, cand.fMomentum.Mag(), fW); + fH.FillH1("hChi2Sts", cand.fMcSrc, cand.fChi2Sts, fW); - if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue; - //fH.FillSourceH1("hPt", cand.fMcSrc, cand.fMomentum.Perp(), fW); - //fH.FillSourceH1("hMom", cand.fMcSrc, cand.fMomentum.Mag(), fW); - fH.FillH1("hChi2Sts", cand.fMcSrc, cand.fChi2sts, fW); + CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd)); + if (stsTrack == nullptr) return; + fH.FillH1("hNofStsHits", cand.fMcSrc, stsTrack->GetNofHits(), fW); - CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd)); - if (stsTrack == nullptr) continue; - fH.FillH1("hNofStsHits", cand.fMcSrc, stsTrack->GetNofHits(), fW); - - if (fUseMvd) { - double mvd1x = 0., mvd1y = 0., mvd2x = 0., mvd2y = 0.; - for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) { - CbmMvdHit* mvdHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM))); - if (mvdHit == nullptr) continue; - if (mvdHit->GetStationNr() == 1) { - mvd1x = mvdHit->GetX(); - mvd1y = mvdHit->GetY(); - } - else if (mvdHit->GetStationNr() == 2) { - mvd2x = mvdHit->GetX(); - mvd2y = mvdHit->GetY(); - } + if (fUseMvd) { + double mvd1x = 0., mvd1y = 0., mvd2x = 0., mvd2y = 0.; + for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) { + CbmMvdHit* mvdHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM))); + if (mvdHit == nullptr) return; + if (mvdHit->GetStationNr() == 1) { + mvd1x = mvdHit->GetX(); + mvd1y = mvdHit->GetY(); + } + else if (mvdHit->GetStationNr() == 2) { + mvd2x = mvdHit->GetX(); + mvd2y = mvdHit->GetY(); } - double mvd1r = sqrt(mvd1x * mvd1x + mvd1y * mvd1y); - double mvd2r = sqrt(mvd2x * mvd2x + mvd2y * mvd2y); - - fH.FillH1("hNofMvdHits", cand.fMcSrc, stsTrack->GetNofMvdHits(), fW); - fH.FillH2("hMvdXY_1", cand.fMcSrc, mvd1x, mvd1y, fW); - fH.FillH1("hMvdR_1", cand.fMcSrc, mvd1r, fW); - fH.FillH2("hMvdXY_2", cand.fMcSrc, mvd2x, mvd2y, fW); - fH.FillH1("hMvdR_2", cand.fMcSrc, mvd2r, fW); } + double mvd1r = sqrt(mvd1x * mvd1x + mvd1y * mvd1y); + double mvd2r = sqrt(mvd2x * mvd2x + mvd2y * mvd2y); + fH.FillH1("hNofMvdHits", cand.fMcSrc, stsTrack->GetNofMvdHits(), fW); + fH.FillH2("hMvdXY_1", cand.fMcSrc, mvd1x, mvd1y, fW); + fH.FillH1("hMvdR_1", cand.fMcSrc, mvd1r, fW); + fH.FillH2("hMvdXY_2", cand.fMcSrc, mvd2x, mvd2y, fW); + fH.FillH1("hMvdR_2", cand.fMcSrc, mvd2r, fW); } } @@ -1221,7 +1781,6 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s if (fCands[iC].IsCutTill(ELmvmAnaStep::ElId)) { CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(fCands[iC].fStsInd)); if (stsTrack == nullptr) continue; - cout << "NofMvdHits = " << stsTrack->GetNofMvdHits() << endl; // TODO: delete for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) { CbmMvdHit* candHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM))); if (candHit != nullptr && candHit->GetStationNr() == mvdStationNum) { @@ -1231,9 +1790,6 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s } } - cout << "mvdV.size() = " << mvdV.size() << endl; // TODO: delete this line - cout << "candV.size() = " << candV.size() << endl; // TODO: delete this line - for (size_t iC = 0; iC < candV.size(); iC++) { LmvmCand& cand = fCands[candV[iC].fInd]; double minD = 9999999.; @@ -1277,9 +1833,6 @@ void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const s if (mvdStationNum == 1) cand.fIsMvd1Cut = isMvdCut; else if (mvdStationNum == 2) cand.fIsMvd2Cut = isMvdCut; - - cout << "cand.fIsMvd1Cut = " << cand.fIsMvd1Cut << endl; // TODO: delete these cout lines - cout << "cand.fIsMvd2Cut = " << cand.fIsMvd2Cut << endl; } } @@ -1311,6 +1864,12 @@ void LmvmTask::MvdCutMcDistance() } } +bool LmvmTask::IsPrimary(double vertexMag) +{ + if (vertexMag < fZ + 0.1 && vertexMag > fZ - 0.1) return true; + return false; +} + void LmvmTask::Finish() { CombinatorialPairs(); @@ -1326,4 +1885,5 @@ void LmvmTask::Finish() void LmvmTask::SetEnergyAndPlutoParticle(const string& energy, const string& particle) { this->SetWeight(LmvmSimParam::GetWeight(energy, particle)); + fParticle = particle; } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h old mode 100644 new mode 100755 index bb60cd1c90..8e8f8fcfbf --- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h @@ -12,6 +12,8 @@ #include "FairRootManager.h" #include "FairTask.h" +#include "CbmGlobalTrack.h" + #include <fstream> #include <map> #include <string> @@ -107,14 +109,22 @@ public: void AssignMcToTopologyCands(std::vector<LmvmCand>& topoCands); - void DifferenceSignalAndBg(); + void DifferenceSignalAndBg(const LmvmCand& cand); /* * \brief Initialize all histograms. */ void InitHists(); - void SignalAndBgReco(); + double MinvScale(double minv); + + void AnalyseCandidates(); + + void AnalyseGlobalTracks(); + void CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isElectron, const std::string& ptcl, double weight); + void BetaMom(const CbmMCTrack* mct, const CbmGlobalTrack* gTrack, const std::string& ptcl); + void CheckTofIdentification(const CbmGlobalTrack* gTrack, const std::string& pidString, double mom, double m2, int pdg, bool isTofEl); + void PidVsMom(const CbmGlobalTrack* gTrack, int iGTrack, int pdg, double mom); void CheckGammaConvAndPi0(); @@ -137,10 +147,22 @@ public: void CombinatorialPairs(); - virtual void Finish(); + void FillCandPidValues(const CbmMCTrack* mcTrack, const LmvmCand& cand, ELmvmAnaStep step); + + void CheckTofId(const CbmMCTrack* mcTrack, const LmvmCand& cand, ELmvmAnaStep step, int pdg); + bool IsInTofPile(double mom, double m2); + + std::string GetPidString(const CbmMCTrack* mct, const LmvmCand* cand); + std::string GetPidString(double vertexMag, int pdg); - void FillAccRecVsMomHist(); + /* + * \brief Check if global track has entries in all detectors. + */ + bool IsInAllDets(const CbmGlobalTrack* gTrack); + + bool IsPrimary(double vertexMag); + virtual void Finish(); ClassDef(LmvmTask, 1); @@ -168,6 +190,7 @@ private: TClonesArray* fTofHits = nullptr; TClonesArray* fTofHitsMatches = nullptr; TClonesArray* fTofPoints = nullptr; + TClonesArray* fTofTracks = nullptr; CbmVertex* fPrimVertex = nullptr; CbmKFVertex fKFVertex; CbmStsKFTrackFitter fKFFitter; @@ -194,6 +217,10 @@ private: std::map<int, int> fNofHitsInRingMap; // Number of hits in the MC RICH ring + double fZ = -44.; // z-position of target + + std::string fParticle = ""; + public: void SetUseMvd(bool use) { fUseMvd = use; } void SetWeight(double w) { fW = w; } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx old mode 100644 new mode 100755 index 620bc0b56b..149a5334bc --- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx @@ -4,6 +4,8 @@ #include "LmvmUtils.h" +#include "Logger.h" + #include "CbmKFVertex.h" #include "CbmL1PFFitter.h" #include "CbmMCTrack.h" @@ -17,6 +19,7 @@ #include <iostream> +//#include "L1Field.h" #include "LmvmCand.h" #include "LmvmDef.h" @@ -31,10 +34,11 @@ void LmvmUtils::CalculateAndSetTrackParams(LmvmCand* cand, CbmStsTrack* stsTrack vector<CbmStsTrack> stsTracks; stsTracks.resize(1); stsTracks[0] = *stsTrack; + //vector<L1FieldRegion> vField; // TODO: this line or next (think of #include "L1Field.h" in header!) vector<CbmL1PFFitter::PFFieldRegion> vField; vector<float> chiPrim; fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, kfVertex, 3e6); - cand->fChi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF(); + cand->fChi2Sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF(); cand->fChi2Prim = chiPrim[0]; const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst(); @@ -139,7 +143,8 @@ bool LmvmUtils::IsMcPairBg(const CbmMCTrack* mctP, const CbmMCTrack* mctM, TClon bool isGamma = IsMcPairGamma(mctP, mctM, mcTracks); bool isEta = IsMcPairEta(mctP, mctM, mcTracks); bool isPi0 = IsMcPairPi0(mctP, mctM, mcTracks); - return (!isEta) && (!isGamma) && (!isPi0) && (!(IsMcSignalEl(mctP) || IsMcSignalEl(mctM))); + //return (!isEta) && (!isGamma) && (!isPi0) && (!(IsMcSignalEl(mctP) && IsMcSignalEl(mctM))); // TODO: this line or next? + return (!isEta) && (!isPi0) && (!(IsMcSignalEl(mctP) && IsMcSignalEl(mctM))); } ELmvmSrc LmvmUtils::GetMcPairSrc(const CbmMCTrack* mctP, const CbmMCTrack* mctM, TClonesArray* mcTracks) @@ -179,7 +184,8 @@ bool LmvmUtils::IsMcPairBg(const LmvmCand& candP, const LmvmCand& candM) bool isGamma = IsMcPairGamma(candP, candM); bool isEta = IsMcPairEta(candP, candM); bool isPi0 = IsMcPairPi0(candP, candM); - return (!isEta) && (!isGamma) && (!isPi0) && (!(candP.IsMcSignal() || candM.IsMcSignal())); + //return (!isEta) && (!isGamma) && (!isPi0) && (!(candP.IsMcSignal() && candM.IsMcSignal())); // TODO: this line or next? + return (!isEta) && (!isPi0) && (!(candP.IsMcSignal() && candM.IsMcSignal())); } ELmvmSrc LmvmUtils::GetMcPairSrc(const LmvmCand& candP, const LmvmCand& candM) @@ -241,19 +247,35 @@ 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 = 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 trdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum); + // Additional TRD Cut + //if (cand->fChi2Trd > 6.) trdEl = false; bool tofEl = CbmLitGlobalElectronId::GetInstance().IsTofElectron(globalTrackIndex, momentum); cand->fMass2 = CbmLitGlobalElectronId::GetInstance().GetTofM2(globalTrackIndex, momentum); - bool isMomCut = (momentumCut > 0.) ? (momentum < momentumCut) : true; + // Additional ToF Cut + /*if (momentum >= 0.5 && momentum <= 2.) { + double slope = 4.; + double b = 0.; + if (cand->fTofDist > momentum*slope + b) tofEl = false; + }*/ + bool isMomCut = (momentumCut > 0.) ? (momentum < momentumCut) : true; cand->fIsElectron = (richEl && trdEl && tofEl && isMomCut); } +void LmvmUtils::IsRichElectron(int globalTrackIndex, double momentum, LmvmCand* cand) { + cand->fIsRichElectron = CbmLitGlobalElectronId::GetInstance().IsRichElectron(globalTrackIndex, momentum); +} + +void LmvmUtils::IsTrdElectron(int globalTrackIndex, double momentum, LmvmCand* cand) { + cand->fIsTrdElectron = (momentum < 1.) ? true : CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum); +} + +void LmvmUtils::IsTofElectron(int globalTrackIndex, double momentum, LmvmCand* cand) { + cand->fIsTofElectron = CbmLitGlobalElectronId::GetInstance().IsTofElectron(globalTrackIndex, momentum); +} + void LmvmUtils::IsElectronMc(LmvmCand* cand, TClonesArray* mcTracks, double pionMisidLevel) { // Use MC information for PID to set a required level of pion suppression @@ -278,3 +300,45 @@ string LmvmUtils::GetChargeStr(const CbmMCTrack* mct) if (mct->GetCharge() == 0) return "0"; return (mct->GetCharge() > 0) ? "+" : "-"; } + +double LmvmUtils::GetMassScaleInmed(double minv) // TODO: make these more elegant and delete cout's +{ + int nArray = sizeof(fMinvArray); + double weight = -1.; + + if (minv < fMinvArray[0]) return fScaleArrayInmed[0]; + else { + for (int i = 1; i < nArray; i++) { + if (fMinvArray[i] > minv) { + double dy = fScaleArrayInmed[i] - fScaleArrayInmed[i-1]; + double dx = fMinvArray[i] - fMinvArray[i-1]; + double slope = dy/dx; + double dLeft = minv - fMinvArray[i-1]; + weight = fScaleArrayInmed[i-1] + slope * dLeft; + return weight; + } + } + } + return weight; +} + +double LmvmUtils::GetMassScaleQgp(double minv) // TODO: make these more elegant and delete cout's +{ + int nArray = sizeof(fMinvArray); + double weight = -1.; + + if (minv < fMinvArray[0]) return fScaleArrayQgp[0]; + else { + for (int i = 1; i < nArray; i++) { + if (fMinvArray[i] > minv) { + double dy = fScaleArrayQgp[i] - fScaleArrayQgp[i-1]; + double dx = fMinvArray[i] - fMinvArray[i-1]; + double slope = dy/dx; + double dLeft = minv - fMinvArray[i-1]; + weight = fScaleArrayQgp[i-1] + slope * dLeft; + return weight; + } + } + } + return weight; +} diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h old mode 100644 new mode 100755 index 27e79ee656..cd04c23d52 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h @@ -88,13 +88,82 @@ public: static void IsElectron(int globalTrackIndex, double momentum, double momentumCut, LmvmCand* cand); + static void IsRichElectron(int globalTrackIndex, double momentum, LmvmCand* cand); + + static void IsTrdElectron(int globalTrackIndex, double momentum, LmvmCand* cand); + + static void IsTofElectron(int globalTrackIndex, double momentum, LmvmCand* cand); + static void IsElectronMc(LmvmCand* cand, TClonesArray* mcTracks, double pionMisidLevel); static std::string GetChargeStr(const LmvmCand* cand); static std::string GetChargeStr(const CbmMCTrack* mct); + + static double GetMassScaleInmed(double minv); + static double GetMassScaleQgp(double minv); + + ClassDef(LmvmUtils, 1); + +private: + + static constexpr double fMinvArray[170] = { + 0.0195, 0.0395, 0.0595, 0.0795, 0.0995, 0.1195, 0.1395, 0.1595, 0.1795, 0.1995, 0.2195, 0.2395, 0.2595, 0.2795, + 0.2995, 0.3195, 0.3395, 0.3595, 0.3795, 0.3995, 0.4195, 0.4395, 0.4595, 0.4795, 0.4995, 0.5195, 0.5395, 0.5595, + 0.5795, 0.5995, 0.6195, 0.6395, 0.6595, 0.6795, 0.6995, 0.7195, 0.7395, 0.7595, 0.7795, 0.7995, 0.8195, 0.8395, + 0.8595, 0.8795, 0.8995, 0.9195, 0.9395, 0.9595, 0.9795, 0.9995, 1.0195, 1.0395, 1.0595, 1.0795, 1.0995, 1.1195, + 1.1395, 1.1595, 1.1795, 1.1995, 1.2195, 1.2395, 1.2595, 1.2795, 1.2995, 1.3195, 1.3395, 1.3595, 1.3795, 1.3995, + 1.4195, 1.4395, 1.4595, 1.4795, 1.4995, 1.5195, 1.5395, 1.5595, 1.5795, 1.5995, 1.6195, 1.6395, 1.6595, 1.6795, + 1.6995, 1.7195, 1.7395, 1.7595, 1.7795, 1.7995, 1.8195, 1.8395, 1.8595, 1.8795, 1.8995, 1.9195, 1.9395, 1.9595, + 1.9795, 1.9995, 2.0195, 2.0395, 2.0595, 2.0795, 2.0995, 2.1195, 2.1395, 2.1595, 2.1795, 2.1995, 2.2195, 2.2395, + 2.2595, 2.2795, 2.2995, 2.3195, 2.3395, 2.3595, 2.3795, 2.3995, 2.4195, 2.4395, 2.4595, 2.4795, 2.4995, 2.5195, + 2.5395, 2.5595, 2.5795, 2.5995, 2.6195, 2.6395, 2.6595, 2.6795, 2.6995, 2.7195, 2.7395, 2.7595, 2.7795, 2.7995, + 2.8195, 2.8395, 2.8595, 2.8795, 2.8995, 2.9195, 2.9395, 2.9595, 2.9795, 2.9995, 3.0195, 3.0395, 3.0595, 3.0795, + 3.0995, 3.1195, 3.1395, 3.1595, 3.1795, 3.1995, 3.2195, 3.2395, 3.2595, 3.2795, 3.2995, 3.3195, 3.3395, 3.3595, + 3.3795, 3.3995}; + + static constexpr double fScaleArrayInmed[170] = { // for 12 AGeV + 41.706, 18.918, 11.465, 8.4388, 5.9176, 4.9025, 3.8087, 3.0387, 2.5856, + 2.1142, 1.7603, 1.5327, 1.28, 1.1579, 1.0367, 0.89355, 0.81317, 0.71582, + 0.65863, 0.59678, 0.53702, 0.45378, 0.41238, 0.37502, 0.33593, 0.28791, 0.26352, + 0.23939, 0.21167, 0.19479, 0.19204, 0.17492, 0.15811, 0.15479, 0.14935, 0.13803, + 0.1354, 0.11993, 0.1046, 0.08226, 0.073183, 0.055433, 0.043467, 0.033975, 0.028025, + 0.021504, 0.016863, 0.014108, 0.01094, 0.0088095, 0.007324, 0.0057162, 0.0046817, 0.0037459, + 0.0030017, 0.0024459, 0.0020671, 0.0016089, 0.0013754, 0.0011223, 0.00096256, 0.00081647, 0.00072656, + 0.00060776, 0.00051243, 0.00045705, 0.00039636, 0.00036259, 0.00033248, 0.0002953, 0.00027328, 0.00023776, + 0.00022163, 0.00019852, 0.000186, 0.00016846, 0.00015469, 0.00014169, 0.00013343, 0.00011594, 0.00010722, + 0.00010205, 9.1907e-05, 8.3718e-05, 7.5457e-05, 6.7192e-05, 6.2202e-05, 5.7372e-05, 4.8314e-05, 4.5502e-05, + 4.1334e-05, 3.7429e-05, 3.2131e-05, 3.0103e-05, 2.6125e-05, 2.3601e-05, 2.1167e-05, 1.94e-05, 1.7025e-05, + 1.5496e-05, 1.3704e-05, 1.1866e-05, 1.1135e-05, 9.8842e-06, 8.9101e-06, 7.9225e-06, 7.0706e-06, 6.3536e-06, + 5.3786e-06, 4.7179e-06, 4.2128e-06, 4.0015e-06, 3.4118e-06, 3.1864e-06, 2.734e-06, 2.3844e-06, 2.173e-06, + 1.8774e-06, 1.6468e-06, 1.501e-06, 1.3597e-06, 1.2113e-06, 1.0384e-06, 9.4105e-07, 8.4223e-07, 7.434e-07, + 6.5049e-07, 5.8824e-07, 5.3603e-07, 4.6756e-07, 4.1173e-07, 3.5872e-07, 3.2764e-07, 2.9889e-07, 2.5989e-07, + 2.219e-07, 1.9468e-07, 1.816e-07, 1.5707e-07, 1.3565e-07, 1.2619e-07, 1.0919e-07, 1.0071e-07, 8.4632e-08, + 7.6459e-08, 6.829e-08, 6.2046e-08, 5.5335e-08, 4.5937e-08, 4.2426e-08, 3.567e-08, 3.4051e-08, 2.9627e-08, + 2.5249e-08, 2.2767e-08, 2.1054e-08, 1.7873e-08, 1.574e-08, 1.3713e-08, 1.23e-08, 1.1045e-08, 9.5536e-09, + 8.5859e-09, 7.7217e-09, 6.9958e-09, 6.0992e-09, 5.3453e-09, 4.7659e-09, 4.3313e-09, 3.6575e-09}; + + static constexpr double fScaleArrayQgp[170] = { // for 12 AGeV + 39.496, 17.961, 11.024, 8.2093, 5.8331, 4.8995, 3.8612, 3.1258, 2.7006, + 2.2465, 1.908, 1.699, 1.4435, 1.3253, 1.2059, 1.049, 0.96753, 0.86685, + 0.81407, 0.75959, 0.70663, 0.61951, 0.58586, 0.55534, 0.51902, 0.46377, 0.4415, + 0.41412, 0.37414, 0.34883, 0.34494, 0.31141, 0.2762, 0.26331, 0.24693, 0.22286, + 0.21697, 0.1972, 0.1841, 0.16097, 0.16352, 0.14345, 0.13096, 0.11911, 0.11399, + 0.10111, 0.0913, 0.08764, 0.077745, 0.071417, 0.067561, 0.05987, 0.055543, 0.050193, + 0.045244, 0.04128, 0.03898, 0.03365, 0.031622, 0.028217, 0.026215, 0.023919, 0.022648, + 0.019915, 0.017524, 0.016145, 0.014357, 0.013362, 0.012368, 0.011036, 0.010198, 0.0088275, + 0.0081762, 0.0072697, 0.00675, 0.0060424, 0.0054788, 0.0049588, 0.0046174, 0.0039685, 0.00363, + 0.0034204, 0.0030534, 0.0027606, 0.0024723, 0.0021893, 0.0020174, 0.0018545, 0.0015584, 0.0014661, + 0.0013315, 0.0012065, 0.0010375, 0.00097456, 0.00084865, 0.00076982, 0.00069371, 0.00063931, 0.00056442, + 0.00051712, 0.00046054, 0.00040174, 0.00037996, 0.00034009, 0.00030921, 0.00027738, 0.00024981, 0.00022659, + 0.00019366, 0.00017153, 0.00015469, 0.00014841, 0.00012783, 0.00012061, 0.00010456, 9.2145e-05, 8.4856e-05, + 7.4087e-05, 6.5675e-05, 6.0496e-05, 5.5386e-05, 4.9865e-05, 4.3202e-05, 3.9571e-05, 3.5821e-05, 3.201e-05, + 2.8322e-05, 2.5886e-05, 2.384e-05, 2.1016e-05, 1.8703e-05, 1.6467e-05, 1.5199e-05, 1.4011e-05, 1.2311e-05, + 1.0621e-05, 9.4155e-06, 8.874e-06, 7.7548e-06, 6.7662e-06, 6.3589e-06, 5.5585e-06, 5.1791e-06, 4.3965e-06, + 4.012e-06, 3.6195e-06, 3.3215e-06, 2.9918e-06, 2.5084e-06, 2.3397e-06, 1.9865e-06, 1.915e-06, 1.6826e-06, + 1.448e-06, 1.3183e-06, 1.231e-06, 1.0551e-06, 9.3811e-07, 8.2511e-07, 7.4714e-07, 6.7735e-07, 5.9142e-07, + 5.3654e-07, 4.8709e-07, 4.4543e-07, 3.9199e-07, 3.4674e-07, 3.1203e-07, 2.862e-07, 2.4391e-07}; - ClassDef(LmvmUtils, 1); }; #endif diff --git a/macro/analysis/dielectron/batch_job.py b/macro/analysis/dielectron/batch_job.py index 7cf0390715..2dcfade19d 100755 --- a/macro/analysis/dielectron/batch_job.py +++ b/macro/analysis/dielectron/batch_job.py @@ -27,8 +27,8 @@ def main(): os.makedirs(workDir) os.chdir(workDir) - plutoFile = getPlutoPath(colSystem, colEnergy, plutoParticle, taskId) - #plutoFile = getPlutoPathTetyana(colSystem, colEnergy, plutoParticle, taskId) + #plutoFile = getPlutoPath(colSystem, colEnergy, plutoParticle, taskId) + plutoFile = getPlutoPathNew(colSystem, colEnergy, plutoParticle, taskId) urqmdFile = "/lustre/nyx/cbm/prod/gen/urqmd/auau/" + colEnergy + "/centr/urqmd.auau."+ colEnergy + ".centr." + str(taskId).zfill(5) + ".root" @@ -70,7 +70,7 @@ def getPlutoPath(colSystem, colEnergy, plutoParticle, taskId): elif plutoParticle == "urqmd": return "" -def getPlutoPathTetyana(colSystem, colEnergy, plutoParticle, taskId): +def getPlutoPathNew(colSystem, colEnergy, plutoParticle, taskId): if plutoParticle == "rho0": return "/lustre/nyx/cbm/users/galatyuk/pluto/epem/" + colEnergy + "/rho0_pluto/sum_rho0_1000_001.zip#rho0_" + str(taskId).zfill(5) + ".root" elif plutoParticle == "omegaepem": diff --git a/macro/analysis/dielectron/batch_send.py b/macro/analysis/dielectron/batch_send.py index df4835a343..91f0c91d94 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 = 1000 + nofJobs = 1000 timeLimit = "08:00:00" geoSetup = "sis100_electron" plutoParticles = ["inmed", "omegadalitz", "omegaepem", "phi", "qgp"] diff --git a/macro/analysis/dielectron/draw_all.py b/macro/analysis/dielectron/draw_all.py index 6ac7220aee..5325946929 100755 --- a/macro/analysis/dielectron/draw_all.py +++ b/macro/analysis/dielectron/draw_all.py @@ -26,7 +26,7 @@ def main(): resultDirAna = resultDirPP + "/lmvm/" inRootAna = dataDir + plutoParticle + "/analysis.all.root" - os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, inRootAna, resultDirAna, useMvd)) + #os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, inRootAna, resultDirAna, useMvd)) resultDirLitqa = resultDirPP + "/litqa/" inRootLitqa = dataDir + plutoParticle + "/litqa.all.root" @@ -40,7 +40,7 @@ def main(): allOmegaD = dataDir + "/omegadalitz/analysis.all.root" resultDirAll = resultDir + "/all/" - #os.system(('root -l -b -q {}/draw_analysis_all.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, allInmed, allQgp, allOmega, allPhi, allOmegaD, resultDirAll, useMvd)) + os.system(('root -l -b -q {}/draw_analysis_all.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\"\)').format(macroDir, allInmed, allQgp, allOmega, allPhi, allOmegaD, resultDirAll, useMvd)) if __name__ == '__main__': main() -- GitLab