CbmRoot
Loading...
Searching...
No Matches
LmvmTask.cxx
Go to the documentation of this file.
1/* Copyright (C) 2010-2021 UGiessen, JINR-LIT
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Andrey Lebedev, Elena Lebedeva, Semen Lebedev [committer], Cornelius Feier-Riesen */
4
5#include "LmvmTask.h"
6
7#include "CbmEvent.h"
8#include "CbmGlobalTrack.h"
9#include "CbmKF.h"
10#include "CbmL1PFFitter.h"
11#include "CbmMCTrack.h"
12#include "CbmMatch.h"
13#include "CbmMvdHit.h"
14#include "CbmRichDigi.h"
15#include "CbmRichHit.h"
16#include "CbmRichPoint.h"
17#include "CbmRichRing.h"
18#include "CbmRichUtil.h"
19#include "CbmStsHit.h"
20#include "CbmStsKFTrackFitter.h"
21#include "CbmStsTrack.h"
22#include "CbmTofHit.h"
23#include "CbmTofPoint.h"
24#include "CbmTofTrack.h"
25#include "CbmTrackMatchNew.h"
26#include "CbmTrdHit.h"
27#include "CbmTrdTrack.h"
28#include "CbmUtils.h"
29#include "CbmVertex.h"
30#include "FairEventHeader.h"
31#include "FairMCPoint.h"
32#include "FairRootManager.h"
33#include "FairRunAna.h"
34#include "FairTask.h"
35#include "FairTrackParam.h"
36#include "LmvmHist.h"
37#include "LmvmSimParam.h"
38#include "LmvmUtils.h"
39#include "TClonesArray.h"
40#include "TDatabasePDG.h"
41#include "TFile.h"
42#include "THnSparse.h"
43#include "TMCProcess.h"
44#include "TRandom3.h"
45#include "TVector3.h"
46#include "cbm/elid/CbmLitGlobalElectronId.h"
47#include "cbm/qa/mc/CbmLitMCTrackCreator.h"
48
49#include <sstream>
50#include <vector>
51
52using namespace std;
53
54ClassImp(LmvmTask);
55
56
57LmvmTask::LmvmTask() : FairTask("LmvmTask") {}
58
59
60LmvmTask::~LmvmTask() {}
61
62
63void LmvmTask::InitHists()
64{
65 string ax = "Yield/Event";
66 string axMinvX = "M_{ee} [GeV/c^{2}]";
67 string axMinvY = "dN/dM_{ee} [GeV/c^{2}]^{-1}";
68
69 // Numbers
70 fH.CreateH1("hEventNumber", "", "", 1, 0, 1.);
71 fH.CreateH1("hJobNumber", "", "", 1, 0, 1.);
72
73
74 // Momentum
75 fH.CreateH1("hMom", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 10.);
76 fH.CreateH2("hMomPt", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "P_{t} [GeV/c]", ax, 100, 0., 10., 100, 0., 4.);
77 fH.CreateH2("hRecoPrec_Mom", fH.fCandNames, fH.fAnaStepNames, "P_{MC} [GeV/c]", "Ratio P_{rec}/P_{MC}", ax, 100, 0.,
78 10., 100., 0., 2.);
79 fH.CreateH2("hRecoPrec_Mom_ID", {"true", "misid"}, "P_{MC} [GeV/c]", "Ratio P_{rec}/P_{MC}", ax, 100, 0., 10., 100.,
80 0., 2.);
81 fH.CreateH2("hMomVsAnglePairSignalMc", "#sqrt{P_{e^{#pm}} P_{e^{#mp}}} [GeV/c]", "#theta_{e^{+},e^{-}} [deg]",
82 "Counter", 100, 0., 5., 100, 0., 50.);
83 fH.CreateH1("hElMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
84 "P [GeV/c]", ax, 100, 0., 10.);
85 fH.CreateH1("hPiMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
86 "P [GeV/c]", ax, 100, 0., 10.);
87 fH.CreateH1("hProtonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
88 "P [GeV/c]", ax, 100, 0., 10.);
89 fH.CreateH1("hKaonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
90 "P [GeV/c]", ax, 100, 0., 10.);
91 fH.CreateH1("hMomPairSignal", fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 15.);
92 for (const string suff : {"", "+", "-"}) {
93 fH.CreateH1("hMom" + suff, fH.fSrcNames, fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 10.);
94 fH.CreateH1("hMomPx" + suff, fH.fSrcNames, fH.fAnaStepNames, "Px [GeV/c]", ax, 120, -3., 3.);
95 fH.CreateH1("hMomPy" + suff, fH.fSrcNames, fH.fAnaStepNames, "Py [GeV/c]", ax, 120, -3., 3.);
96 fH.CreateH1("hMomPz" + suff, fH.fSrcNames, fH.fAnaStepNames, "Pz [GeV/c]", ax, 120, -1., 11.);
97 fH.CreateH1("hPt" + suff, fH.fSrcNames, fH.fAnaStepNames, "P_{t} [GeV/c]", ax, 100, 0., 4.);
98 fH.CreateH1("hRapidity" + suff, fH.fSrcNames, fH.fAnaStepNames, "Rapidity", ax, 100, 0., 5.);
99 if (suff == "") continue;
100 fH.CreateH1("hMom" + suff + "_trueEl_bg_elid", "P [GeV/c]", ax, 100, 0., 10.);
101 }
102 fH.CreateH1("hMomAcc+", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.);
103 fH.CreateH1("hMomAcc-", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.);
104
105
106 // PDG
107 fH.CreateH1("hMotherPdg_mc", "Pdg code", "Particles / Event", 7000, -3500., 3500.);
108 fH.CreateH1("hCandPdg", fH.fAnaStepNames, "Pdg code", "Particles / Event", 7001, -3500., 3500.);
109 fH.CreateH2("hCandPdgVsMom", fH.fAnaStepNames, "P [GeV/c]", "Particle", "Yield/(Event * Bin)", 100, 0., 10., 7, 0.,
110 7.);
111 fH.CreateH2("hCandPdgVsMom+", "P [GeV/c]", "Particle", "Yield/(Event * Bin)", 100, 0., 10., 7, 0., 7.); // only El-ID
112 fH.CreateH2("hCandPdgVsMom-", "P [GeV/c]", "Particle", "Yield/(Event * Bin)", 100, 0., 10., 7, 0., 7.); // only El-ID
113
114 fH.CreateH2("hCandElSrc", "Analysis step", "Mother of Electron Candidate", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps,
115 8, 0., 8.);
116 fH.CreateH2("hPdgVsMom", {"rich", "trd", "tof"}, {"rec", "acc", "chi2prim"}, "P [GeV/c]", "Misid. Particle",
117 "Yield/(Event * Bin)", 100, 0., 10., 6, 0., 6.);
118
119
120 // Vertex
121 fH.CreateH2("hVertexGammaXZ", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 100, -10., 190., 100, -130., 130.);
122 fH.CreateH2("hVertexGammaYZ", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 100, -10., 190., 100, -130., 130.);
123 fH.CreateH2("hVertexGammaXY", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 100, -130., 130., 100, -130., 130.);
124 fH.CreateH2("hVertexGammaRZ", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 50,
125 0., 50.);
126 fH.CreateH2("hVertexGTrackRZ", "Z [cm]", "R [cm]", ax, 100, -50., 100., 100, -50., 50.);
127
128 // Vertices of not-electron candidates that are misidentified as electrons in ToF
129 fH.CreateH2("hVertexXZ_misidTof", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.);
130 fH.CreateH2("hVertexYZ_misidTof", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.);
131 fH.CreateH2("hVertexXY_misidTof", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 100, -50., 50., 100, -50., 50.);
132 fH.CreateH2("hVertexRZ_misidTof", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100.,
133 100, 0., 50.);
134
135
136 // Numbers
137 fH.CreateH1("hNofMcTracks", "Particle", "Particles / Event", 7, 0., 7.);
138 fH.CreateH1("hNofGTracks", "Particle", "Particles / Event", 6, 0., 6.);
139 fH.CreateH1("hNofBgTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
140 fH.CreateH1("hNofSignalTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
141 fH.CreateH2("hBgSrcTracks", "Analysis step", "Candidate Source", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.);
142 fH.CreateH1("hNofTopoPairs", {"gamma", "pi0"}, "Pair type", "Pairs/event", 8, 0., 8);
143 fH.CreateH1("hNofMismatches", {"all", "rich", "trd", "tof"}, "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0.,
144 fH.fNofAnaSteps);
145 fH.CreateH1("hNofMismatches_gTracks", fH.fGTrackNames, "Detector", ax, 7., 0., 7.);
146 fH.CreateH1("hNofMismatchedTrackSegments", fH.fGTrackNames, "nof mism. track segments", ax, 4., 0., 4.);
147 fH.CreateH2("hMatchId_gTracks", fH.fGTrackNames, "Nof mismatched Track Segments", "Identification", ax, 4., 0., 4.,
148 2., 0., 2.);
149 fH.CreateH1("hNofGhosts", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
150 fH.CreateH1("hNofMvdHits", fH.fSrcNames, "Number of hits in MVD", ax, 5, -0.5, 4.5);
151 fH.CreateH1("hNofStsHits", fH.fSrcNames, "Number of hits in STS", ax, 9, -0.5, 8.5);
152 fH.CreateH2("hNofHitsVsMom", {"sts", "rich", "trd", "tof"}, fH.fGTrackNames, "P [GeV/c]", "Number of Hits", ax, 100.,
153 0., 10., 50, -0.5, 49.5);
154
155
156 // Background
157 fH.CreateH2("hSrcBgPairs", "Analysis step", "Pair", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, fH.fNofBgPairSrc, 0.,
158 fH.fNofBgPairSrc);
159 fH.CreateH2("hSrcBgPairsEpEm", fH.fAnaStepNames, "mother particle e+", "mother particle e-", ax, 4, 0., 4., 4, 0.,
160 4.);
161 for (int iMinv = 0; iMinv <= 24; iMinv++) {
162 for (const string comb : {"PM", "PP", "MM"}) {
163 for (const string type : {"pMc", "pRec"}) { // TODO: delete one of these histos
164 string hName = "hBgPairPdg" + comb + "_" + type + "_" + Cbm::NumberToString(iMinv, 0);
165 fH.CreateH2(hName, "PDG e^{+}", "PDG e^{-}", "Yield", 8, 0., 8., 8, 0., 8.);
166 }
167 }
168 }
169 fH.CreateH1("hMinvBgMatch", {"trueMatch", "trueMatchEl", "trueMatchNotEl", "mismatch"}, fH.fAnaStepNames, axMinvX, ax,
170 250, 0., 2.5);
171 fH.CreateH1("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, axMinvX, axMinvY, 250, 0., 2.5);
172 fH.CreateH1("hMinvBgSource2", fH.fAnaStepNames,
173 {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"}, axMinvX, axMinvY, 250, 0.,
174 2.5); // "pi" are misid. charged pions
175
176
177 // Chi2
178 fH.CreateH1("hChi2_truematch", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.);
179 fH.CreateH1("hChi2_mismatch", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.);
180 fH.CreateH1("hChi2Sts", fH.fSrcNames, "#chi^{2}", ax, 200, 0., 20.);
181 fH.CreateH1("hChi2PrimVertex", fH.fSrcNames, "#chi^{2}_{prim}", ax, 200, 0., 20.);
182 for (const string det : {"sts", "rich", "trd", "tof"}) {
183 fH.CreateH2("hChi2Dets_" + det, fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "#chi^{2}", ax, 100., 0., 10., 100,
184 0., 20.);
185 }
186 fH.CreateH2("hChi2Comb", {"StsRich", "StsTrd", "StsTof", "RichTrd", "RichTof", "TrdTof"}, fH.fCandNames,
187 "#chi^{2}_{1}", "#chi^{2}_{2}", ax, 100., 0., 20., 100, 0., 20.);
188
189
190 // Minv for PLUTO + UrQMD / UrQMD / UrQMD electrons
191 for (const string cat : {"", "_trueEl", "_urqmdAll", "_urqmdEl"}) {
192 fH.CreateH1("hMinv" + cat, fH.fSrcNames, fH.fAnaStepNames, axMinvX, axMinvY, 250, 0., 2.5, true);
193 }
194 for (const string comb : {"PM", "PP", "MM"}) { // to compare LmvmEventMix yields with LmvmTask yields
195 fH.CreateH1("hMinv" + comb, {"sameEv", "mixedEv"}, axMinvX, axMinvY, 250, 0., 2.5, true);
196 }
197
198 fH.CreateH2("hRecoPrec_Minv", "M_{ee, MC} [GeV/c^{2}]", "Ratio M_{ee, rec}/M_{ee, MC}", ax, 250, 0., 2.5, 200., 0.,
199 2.);
200 fH.CreateH1("hMinv_unscaled", axMinvX, axMinvY, 250, 0., 2.5); // for signals only
201 fH.CreateH1("hMinv", {"tt", "tm", "mt", "mm"}, axMinvX, axMinvY, 250, 0., 2.5, true); // xx = +-
202
203
204 // Topology Cuts
205 fH.CreateH2("hTtCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
206 "#theta_{e^{+},e^{-}} [deg]", ax, 100, 0., 5., 100, 0., 5.);
207 fH.CreateH2("hStCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
208 "#theta_{e^{#pm},rec} [deg]", ax, 100, 0., 5., 100, 0., 5.);
209 fH.CreateH2("hRtCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
210 "#theta_{e^{#pm},rec} [deg]", ax, 100, 0., 5., 100, 0., 5.);
211
212
213 // Combinatorial Background
214 for (const string comb : {"PM", "PP", "MM"}) {
215 string hName = "hMinvComb" + comb;
216 for (const string cat : {"", "_trueEl", "_urqmdAll", "_urqmdEl"}) {
217 string hNameCat = hName + cat;
218 fH.CreateH1(hNameCat, {"sameEv", "mixedEv"}, fH.fAnaStepNames, axMinvX, axMinvY, 250, 0., 2.5, true);
219 }
220 fH.CreateH2(hName, {"px", "py", "px_mother", "py_mother"}, "P [GeV/c]", "P [GeV/c]", ax, 30, -3., 3., 30, -3., 3.);
221 fH.CreateH2(hName, {"R", "R_mother"}, "Start Vertex #sqrt{X_{1}^{2}+Y_{1}^{2}} [cm]",
222 "Start Vertex #sqrt{X_{2}^{2}+Y_{2}^{2}} [cm]", ax, 50., 0., 50., 50., 0., 50.);
223 fH.CreateH2(hName + "_theta-minv", axMinvX, "#theta [#circ]", ax, 100, 0., 2.5, 20, 0., 80.);
224 }
225 fH.CreateH1("hMinvCombPM_same", {"tt", "tm", "mt", "mm"}, axMinvX, axMinvY, 250, 0., 2.5); // tm = true e+, misid -
226
227
228 // Histograms for Fast Simulations (see manual for naming conventions)
229 const Int_t nDim = 3; // Dimension of momentum histograms
230 int nBinsMult = 6; // Number of bins in multiplicity histogram. Maximum multiplicity is nBinsMult-1
231 int nBinsMom = 5000; // Number of bins for momentum histogram
232 std::array<Int_t, nDim> FSnBinsMom = {nBinsMom, nBinsMom, nBinsMom};
233 std::array<Double_t, nDim> FSMomMin = {-3., -3., 0.}; // Minimum values for momentum in X, Y, Z
234 std::array<Double_t, nDim> FSMomMax = {3., 3., 10.}; // Maximum values for momentum in X, Y, Z
235 for (ELmvmAnaStep step : fH.fAnaStepsFS) {
236 string stepname = fH.fAnaStepNames[static_cast<int>(step)];
237
238 // charge-based and Plutos
239 for (const string charge : {"plus", "minus", "plusEl", "minusEl"}) {
240 // UrQMD
241 string hNameMultUrqmd = "hfsc_mult_urqmd_urqmd_" + charge + "_" + stepname;
242 string hNameMomUrqmd = "hfsc_mom_urqmd_urqmd_" + charge + "_" + stepname;
243 string titleMomUrqmd = hNameMomUrqmd + "; P_{x} [GeV/c]; P_{y} [GeV/c]; P_{z} [GeV/c]";
244 fH.CreateH1(hNameMultUrqmd, "N / Event", ax, nBinsMult, 0., nBinsMult);
245 fH.fHM.CreateSparse<THnSparseD, nDim>(hNameMomUrqmd, titleMomUrqmd, FSnBinsMom, FSMomMin, FSMomMax);
246 // Plutos
247 if (charge == "plusEl" || charge == "minusEl") continue;
248 string hNameMultPluto = "hfsc_mult_pluto_" + fParticle + "_" + charge + "_" + stepname;
249 string hNameMomPluto = "hfsc_mom_pluto_" + fParticle + "_" + charge + "_" + stepname;
250 string titleMomPluto = hNameMomPluto + "; P_{x} [GeV/c]; P_{y} [GeV/c]; P_{z} [GeV/c]";
251 fH.CreateH1(hNameMultPluto, "N / Event", ax, nBinsMult, 0., nBinsMult);
252 fH.fHM.CreateSparse<THnSparseD, nDim>(hNameMomPluto, titleMomPluto, FSnBinsMom, FSMomMin, FSMomMax);
253 }
254
255 // particle-based (Plutos will be taken from charge-based)
256 for (const string& ptcl : fH.fFSCandNames) {
257 for (const string charge : {"plus", "minus"}) {
258 string hNameMultUrqmd = "hfsp_mult_urqmd_" + ptcl + "_" + charge + "_" + stepname;
259 string hNameMomUrqmd = "hfsp_mom_urqmd_" + ptcl + "_" + charge + "_" + stepname;
260 string titleMomUrQmd = hNameMomUrqmd + "; P_{x} [GeV/c]; P_{y} [GeV/c]; P_{z} [GeV/c]";
261 fH.CreateH1(hNameMultUrqmd, "N / Event", ax, nBinsMult, 0., nBinsMult);
262 fH.fHM.CreateSparse<THnSparseD, nDim>(hNameMomUrqmd, titleMomUrQmd, FSnBinsMom, FSMomMin, FSMomMax);
263 }
264 }
265 }
266
267
268 // Like-sign correlations
269 fH.CreateH1("hLikeSignCorr", fH.fAnaStepNames, "Decay Mode", "N / Event", 6, 0., 6.);
270
271
272 // STS
273 fH.CreateH2("hELossSts", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "E-Loss [e/300 #mum]", ax, 100., 0., 10., 100,
274 1e4, 1e5);
275
276
277 // MVD
278 fH.CreateH2("hMvdCut", {"1", "2"}, fH.fSrcNames, "d_{MVD} [cm]", "P_{e} [GeV/c]", ax, 100, 0., 1., 100, 0., 5.);
279 fH.CreateH2("hMvdXY", {"1", "2"}, fH.fSrcNames, "X [cm]", "Y [cm]", ax, 60, -6., 6., 60, -6., 6.);
280 fH.CreateH1("hMvdR", {"1", "2"}, fH.fSrcNames, "#sqrt{X^{2}+Y^{2}} [cm]", ax, 60, 0., 6.);
281 fH.CreateH1("hMvdCutQa", {"1", "2"}, fH.fSrcNames, "MVD hit assignment", ax, 2, 0., 2.); //[0.5]-correct, [1.5]-wrong
282 fH.CreateH1("hMvdMcDist", {"1", "2"}, fH.fSrcNames, "Track-Hit distance [cm]", ax, 100, 0., 10.);
283
284
285 // RICH
286 fH.CreateH2("hPmtXY", fH.fSrcNames, "X [cm]", "Y [cm]", "Counter", 100, -110, 110, 100, -200, 200);
287 fH.CreateH1("hNofRichHits", "NofRichHits / Event", ax, 30, 0., 3000.);
288 fH.CreateH1("hNofRichRings", "NofRichRings / Event", ax, 30, 0., 60.);
289 fH.CreateH2("hRichAnnSrc", fH.fSrcNames, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 10., 100, -1.1, 1.1);
290 fH.CreateH2("hRichAnn", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 10., 100, -1.1,
291 1.1);
292 fH.CreateH1("hRichRingTrackDist", fH.fCandNames, fH.fAnaStepNames, "D [cm]", ax, 100, 0., 20.);
293 fH.CreateH2("hRichRingTrackDistVsMom", fH.fCandNames, fH.fAnaStepNames, "P_{rec} [GeV/c]", "D [cm]", "Yield", 100, 0.,
294 10., 100, 0., 20.);
295 fH.CreateH1("hRichRingBoA", fH.fAnaStepNames, "B/A", ax, 24, 0., 1.2);
296 fH.CreateH2("RichRingNeighbour_all", {"next", "secondnext"}, fH.fCandNames, "r [cm]", "Ring-Track-Dist [cm]", ax, 150,
297 100., 250., 100, 0., 20.);
298 fH.CreateH2("RichRingNeighbour_elid", {"next", "secondnext"}, fH.fCandNames, "r [cm]", "Ring-Track-Dist [cm]", ax,
299 150, 100., 250., 100, 0., 20.);
300
301
302 // TRD
303 fH.CreateH2("hTrdElLike", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "Likelihood output", ax, 100, 0., 10., 100,
304 0., 1.);
305
306
307 // TOF
308 fH.CreateH2("hTofM2Src", fH.fSrcNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 6., 100, -0.3, 1.0);
309 fH.CreateH2("hTofM2", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 400,
310 -0.5, 1.);
311 fH.CreateH2("hTofM2_woDoubleIndex", {"1", "2", "3+"}, fH.fCandNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100,
312 0., 10., 300, -0.5, 1.);
313 fH.CreateH2("hTofM2", {"truematch", "mismatch"}, fH.fCandNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0.,
314 10., 100, -0.5, 1.);
315 fH.CreateH2("hTofDist", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "Distance (norm.)", ax, 100, 0., 10., 100, 0.,
316 100.);
317 fH.CreateH2("hTofDist", {"truematch", "mismatch"}, fH.fCandNames, "P [GeV/c]", "Distance (norm.)", ax, 100, 0., 10.,
318 100, 0., 100.);
319 fH.CreateH2("hTofM2_chi2prim", {"50", "100", "150", "200", "300", "400", "500"}, fH.fCandNames, "P [GeV/c]",
320 "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 100, -0.5, 1.);
321 fH.CreateH2("hTofM2_elid", {"50", "100", "150", "200", "300", "400", "500"}, fH.fCandNames, "P [GeV/c]",
322 "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 100, -0.5, 1.);
323 fH.CreateH2("hTofChi2", {"all", "elid"}, fH.fCandNames, "P [GeV/c]", "#chi^2_{TOF}", ax, 100, 0., 10., 100., 0., 20.);
324 fH.CreateH1("hMatch", {"rich", "trd", "tof"}, "Match", ax, 2, 0., 2.);
325
326
327 // compare misid. candidates (pions, proton, ...) with not misid. ones
328 for (size_t iP = 4; iP < fH.fCandNames.size(); iP++) {
329 fH.CreateH1("hMom_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", ax, 100, 0., 10.);
330 fH.CreateH2("hPtY_" + fH.fCandNames[iP], {"true", "misid"}, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 40, 0.,
331 2.);
332 fH.CreateH2("hTofM2_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0.,
333 10., 100, -0.1, 1.);
334 }
335
336
337 // Misc.
338 fH.CreateH2("hElossTime_chi2prim", {"Seg1Mom1", "Seg1Mom2", "Seg2Mom1", "Seg2Mom2"}, fH.fCandNames, "Time (ct) [m]",
339 "E-Loss [e/300 #mum]", ax, 80, 8., 16., 100, 1e4, 1e5);
340 fH.CreateH2("hElossTime_elid", {"Seg1Mom1", "Seg1Mom2", "Seg2Mom1", "Seg2Mom2"}, fH.fCandNames, "Time (ct) [m]",
341 "E-Loss [e/300 #mum]", ax, 80, 8., 16., 100, 1e4, 1e5);
342 fH.CreateH2("hMinvPt", fH.fSrcNames, fH.fAnaStepNames, axMinvX, "P_{t} [GeV/c]", ax, 100, 0., 2., 25, 0., 2.5);
343 fH.CreateH2("hPtYPairSignal", fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
344 fH.CreateH1("hAnglePair", fH.fSrcNames, fH.fAnaStepNames, "#Theta_{1,2} [deg]", ax, 80, 0., 80.);
345 fH.CreateH1("hSuppression", {"pion", "proton"}, {"mc", "reco", "acc", "elid", "ttcut", "ptcut"}, "P_{MC} [GeV/c]", ax,
346 20, 0., 10.);
347 fH.CreateH2("hPtY", fH.fCandNames, fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
348 fH.CreateH2("hBetaMom", fH.fCandNames, "P * Q ", "#beta", ax, 200, -10., 10., 150., 0., 1.5);
349 fH.CreateH2("hCandProperties", {"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11"}, fH.fCandNames, "Property",
350 "Value", ax, 16, 0., 16., 100, 0., 10.);
351}
352
353
354InitStatus LmvmTask::Init()
355{
356 fMCEventHeader = InitOrFatal<FairMCEventHeader>("MCEventHeader.");
357 fCbmEvents = InitOrFatal<TClonesArray>("CbmEvent");
358 fMCTracks = InitOrFatal<TClonesArray>("MCTrack");
359 fRichHits = InitOrFatal<TClonesArray>("RichHit");
360 fRichRings = InitOrFatal<TClonesArray>("RichRing");
361 fRichPoints = InitOrFatal<TClonesArray>("RichPoint");
362 fRichRingMatches = InitOrFatal<TClonesArray>("RichRingMatch");
363 fRichProj = InitOrFatal<TClonesArray>("RichProjection");
364 fStsTrackMatches = InitOrFatal<TClonesArray>("StsTrackMatch");
365 fStsTracks = InitOrFatal<TClonesArray>("StsTrack");
366 fStsHits = InitOrFatal<TClonesArray>("StsHit");
367 if (fUseMvd) {
368 fMvdHits = InitOrFatal<TClonesArray>("MvdHit");
369 fMvdPoints = InitOrFatal<TClonesArray>("MvdPoint");
370 fMvdHitMatches = InitOrFatal<TClonesArray>("MvdHitMatch");
371 }
372 fGlobalTracks = InitOrFatal<TClonesArray>("GlobalTrack");
373 fTrdTracks = InitOrFatal<TClonesArray>("TrdTrack");
374 fTrdTrackMatches = InitOrFatal<TClonesArray>("TrdTrackMatch");
375 fTofPoints = InitOrFatal<TClonesArray>("TofPoint");
376 fTofHits = InitOrFatal<TClonesArray>("TofHit");
377 fTofHitsMatches = InitOrFatal<TClonesArray>("TofHitMatch");
378 fTofTracks = InitOrFatal<TClonesArray>("TofTrack");
379 fPrimVertex = InitOrFatal<CbmVertex>("PrimaryVertex.");
380 fDigiManager = CbmDigiManager::Instance();
381 fDigiManager->Init();
382
383 InitHists();
384
385 fKFFitter.Init();
386 CbmLitMCTrackCreator::Instance();
387 CbmLitGlobalElectronId::GetInstance();
388
389 return kSUCCESS;
390}
391
392void LmvmTask::Exec(Option_t*)
393{
394 fH.FillH1("hEventNumber", 0.5);
395 fEventNumber++;
396 // bool useMbias = false; // false for 40% central agag collisions (b<7.7fm)
397 // bool isCentralCollision = false;
398
399 // if (!useMbias) {
400 // double impactPar = fMCEventHeader->GetB();
401 // if (impactPar <= 7.7) isCentralCollision = true;
402 // }
403
404 cout << "========================================================" << endl;
405 LOG(info) << "LmvmTask event number " << fEventNumber;
406 LOG(info) << "nMcTracks = " << fMCTracks->GetEntriesFast();
407 LOG(info) << "nGlobalTracks = " << fGlobalTracks->GetEntriesFast();
408 LOG(info) << "nRichRings = " << fRichRings->GetEntriesFast();
409 LOG(info) << "nTrdTacks = " << fTrdTracks->GetEntriesFast();
410 LOG(info) << "nTofTracks = " << fTofTracks->GetEntriesFast();
411 LOG(info) << "nTofHits = " << fTofHits->GetEntriesFast();
412 LOG(info) << "fPionMisidLevel = " << fPionMisidLevel << endl;
413 //LOG(info) << fCuts.ToString();
414 LOG(info) << "fW = " << fW;
415
416 if (fPrimVertex != nullptr) {
417 fKFVertex = CbmKFVertex(*fPrimVertex);
418 }
419 else {
420 Fatal("LmvmTask::Exec", "No PrimaryVertex array!");
421 }
422
423 //if (useMbias || (!useMbias && isCentralCollision)) {
424 DoMcTrack();
425 DoMcPair();
426 RichPmtXY();
427 FillTopologyCands();
428 FillCands();
429 CalculateNofTopologyPairs("hNofTopoPairs_gamma", ELmvmSrc::Gamma);
430 CalculateNofTopologyPairs("hNofTopoPairs_pi0", ELmvmSrc::Pi0);
431 AnalyseGlobalTracks();
432 AnalyseCandidates();
433 AnalyseProperties();
434 AnalyseRichRings();
435 FillRichRingNofHits();
436 CheckLikeSignCorrelations();
437 FillHistosForFastSim();
438 FillCandsForEventMix();
439 //}
440} // Exec
441
442
443void LmvmTask::FillCandsForEventMix()
444{
445 for (size_t iC = 0; iC < fCands.size(); iC++) {
446 const LmvmCand& cand = fCands[iC];
447 if (cand.IsCutTill(ELmvmAnaStep::ElId)) fCandsTotal.push_back(cand);
448 }
449 LOG(info) << "fCandsTotal.size = " << fCandsTotal.size();
450}
451
452void LmvmTask::FillRichRingNofHits()
453{
454 fNofHitsInRingMap.clear();
455 int nofRichHits = fRichHits->GetEntriesFast();
456 fH.FillH1("hNofRichHits", nofRichHits);
457 for (int iHit = 0; iHit < nofRichHits; iHit++) {
458 CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
459 if (hit == nullptr || hit->GetRefId() < 0) continue;
460 Int_t digiIndex = hit->GetRefId();
461 if (digiIndex < 0) continue;
462 const CbmRichDigi* digi = fDigiManager->Get<CbmRichDigi>(digiIndex);
463 if (!digi) continue;
464 const CbmMatch* digiMatch = fDigiManager->GetMatch(ECbmModuleId::kRich, digiIndex);
465 if (!digiMatch) continue;
466 vector<CbmLink> links = digiMatch->GetLinks();
467 for (const auto& link : links) {
468 if (link.IsNoise()) continue;
469 const CbmRichPoint* point = static_cast<const CbmRichPoint*>(fRichPoints->At(link.GetIndex()));
470 if (!point) continue;
471 Int_t mcTrackId = point->GetTrackID();
472 if (mcTrackId < 0) continue;
473 const CbmMCTrack* mcTrack = static_cast<const CbmMCTrack*>(fMCTracks->At(mcTrackId));
474 if (!mcTrack) continue;
475 if (mcTrack->GetPdgCode() != 50000050) continue; // select only Cherenkov photons
476 Int_t motherId = mcTrack->GetMotherId();
477 fNofHitsInRingMap[motherId]++;
478 }
479 }
480}
481
482double LmvmTask::MinvScale(const CbmMCTrack* mct, const string& signal)
483{
484 if (!LmvmUtils::IsMcSignalEl(mct)) {
485 LOG(warn) << "LmvmTask::MinvScale(CbmMCTrack*): MCTrack is not signal! Will return -1.";
486 return -1.;
487 }
488 int nMcTracks = fMCTracks->GetEntries();
489 for (int iMc2 = 0; iMc2 < nMcTracks; iMc2++) {
490 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2));
491 if (!LmvmUtils::IsMcSignalEl(mct2)) continue;
492 if (mct->GetCharge() * mct2->GetCharge() < 0.) {
493 LmvmKinePar pKin = LmvmKinePar::Create(mct, mct2);
494 return LmvmUtils::MinvScale(pKin.fMinv, signal);
495 }
496 }
497 return -1.;
498}
499
500void LmvmTask::CheckLikeSignCorrelations()
501{
502 /*size_t nCand = fCands.size();
503 int nCharge1 = 0;
504 for (auto step : fH.fAnaSteps) {
505 for (size_t iC1 = 0; iC1 < nCand-1; iC1++) {
506 const auto& cand1 = fCands[iC1];
507 if (!cand1.IsCutTill(step)) continue;
508 if (std::abs(cand1.fMcPdg) != 11) continue;
509 CbmMCTrack* mct1 = GetMcTrackSts(cand1.fStsInd);
510 if (mct1 == nullptr) continue;
511 CbmMCTrack* mother1 = static_cast<CbmMCTrack*>(fMCTracks->At(mct1->GetMotherId()));
512 if (mother1 == nullptr) continue;
513 CbmMCTrack* grmother1 = static_cast<CbmMCTrack*>(fMCTracks->At(mother1->GetMotherId()));
514 int pdgGrMother1 = (grmother1 == nullptr) ? -9999999 : grmother1->GetPdgCode();
515
516 for (size_t iC2 = iC1+1; iC2 < nCand; iC2++) {
517 const auto& cand2 = fCands[iC2];
518 if (!cand2.IsCutTill(step)) continue;
519 if (std::abs(cand2.fMcPdg) != 11) continue;
520 if (cand1.fCharge != cand2.fCharge) continue;
521 CbmMCTrack* mct2 = GetMcTrackSts(cand2.fStsInd);
522 if (mct2 == nullptr) continue;
523 CbmMCTrack* mother2 = static_cast<CbmMCTrack*>(fMCTracks->At(mct2->GetMotherId()));
524 if (mother2 == nullptr) continue;
525 CbmMCTrack* grmother2 = static_cast<CbmMCTrack*>(fMCTracks->At(mother2->GetMotherId()));
526
527 bool isPi0To2Gamma = (cand1.IsMcGamma() && cand2.IsMcGamma() && mother1->GetMotherId() == mother2->GetMotherId() && pdgGrMother1 == 111);
528 bool isPi0Dalitz = ((cand1.IsMcPi0() && cand2.IsMcGamma() && (mct1->GetMotherId() == mother2->GetMotherId())) ||
529 (cand2.IsMcPi0() && cand1.IsMcGamma() && (mct2->GetMotherId() == mother1->GetMotherId())));
530 bool isEtaDalitz = ((cand1.IsMcEta() && cand2.IsMcGamma() && (mct1->GetMotherId() == mother2->GetMotherId())) ||
531 (cand2.IsMcEta() && cand1.IsMcGamma() && (mct2->GetMotherId() == mother1->GetMotherId())));
532
533 if (isPi0To2Gamma) fH.FillH1("hLikeSignCorr", step, 0.5);
534 if (isPi0Dalitz) fH.FillH1("hLikeSignCorr", step, 1.5);
535 if (isEtaDalitz) fH.FillH1("hLikeSignCorr", step, 2.5);
536 if (mct1->GetMotherId() == mct2->GetMotherId()) fH.FillH1("hLikeSignCorr", step, 3.5);
537 if ((mct1->GetMotherId() == mother2->GetMotherId() || mct2->GetMotherId() == mother1->GetMotherId()) && !isPi0Dalitz && !isEtaDalitz) fH.FillH1("hLikeSignCorr", step, 4.5);
538 if (mother1->GetMotherId() == mother2->GetMotherId() && !isPi0To2Gamma && grmother1 != nullptr && grmother2 != nullptr) {
539 fH.FillH1("hLikeSignCorr", step, 5.5);
540 cout << "CheckLikeSignCorrelations(): PDG McT1 = " << mct1->GetPdgCode() << "; PDG Mother1 = " << mother1->GetPdgCode() << "; PDG GrMother1 = " << grmother1->GetPdgCode() << endl;
541 cout << "CheckLikeSignCorrelations(): PDG McT2 = " << mct2->GetPdgCode() << "; PDG Mother2 = " << mother2->GetPdgCode() << "; PDG GrMother2 = " << grmother2->GetPdgCode() << endl;
542 }
543 }
544 }
545 }
546 if (nCharge1 > 0) LOG(warn) << "LmvmTask::CheckLikeSignCorrelations(): " << nCharge1 << " Candidates with Charge=0 are in sample. CHECK!";*/
547}
548
549void LmvmTask::FillHistosForFastSim()
550{
551 int mPlus, mMinus, mPlusEl, mMinusEl; // charge-based
552 int mPi0ElP, mGammaElP, mEtaElP, mOtherElP, mPionP, mKaonP, mOtherP, mProton; // positive charged UrQMD particles
553 int mPi0ElM, mGammaElM, mEtaElM, mOtherElM, mPionM, mKaonM, mOtherM; // negative charged UrQMD particles
554
555 for (ELmvmAnaStep step : fH.fAnaStepsFS) {
556 string stepname = fH.fAnaStepNames[static_cast<int>(step)];
557
558 // charge-based multiplicity counters
559 mPlus = 0; // Multiplicity of positive UrQMD particles
560 mMinus = 0; // Multiplicity of negative UrQMD particles
561 mPlusEl = 0; // Multiplicity of true UrQMD positrons (pdg=-11)
562 mMinusEl = 0; // Multiplicity of true UrQMD electrons (pdg=11)
563
564 // particle-based multiplicity counters
565 mPi0ElP = 0;
566 mPi0ElM = 0;
567 mGammaElP = 0;
568 mGammaElM = 0;
569 mEtaElP = 0;
570 mEtaElM = 0;
571 mOtherElP = 0;
572 mOtherElM = 0;
573 mPionP = 0;
574 mPionM = 0;
575 mProton = 0;
576 mKaonP = 0;
577 mKaonM = 0;
578 mOtherP = 0;
579 mOtherM = 0;
580
581 for (const auto& cand : fCands) {
582 if (!cand.IsCutTill(step)) continue;
583 int pdg = cand.fMcPdg;
584 double w = cand.fWeight;
585 double mom[3] = {cand.fMomentum.X(), cand.fMomentum.Y(), cand.fMomentum.Z()};
586
587 // Pluto particles
588 if (cand.IsMcSignal()) {
589 if (cand.fCharge > 0) {
590 fH.FillH1("hfsc_mult_pluto_" + fParticle + "_plus_" + stepname, 1.5, w);
591 fH.fHM.HnSparse("hfsc_mom_pluto_" + fParticle + "_plus_" + stepname)->Fill(mom, w);
592 }
593 else if (cand.fCharge < 0) {
594 fH.FillH1("hfsc_mult_pluto_" + fParticle + "_minus_" + stepname, 1.5, w);
595 fH.fHM.HnSparse("hfsc_mom_pluto_" + fParticle + "_minus_" + stepname)->Fill(mom, w);
596 }
597 }
598
599 // UrQMD particles
600 else {
601 // Charge-based
602 if (cand.fCharge > 0) {
603 fH.fHM.HnSparse("hfsc_mom_urqmd_urqmd_plus_" + stepname)->Fill(mom, w);
604 mPlus++;
605 if (std::abs(pdg) == 11) {
606 fH.fHM.HnSparse("hfsc_mom_urqmd_urqmd_plusEl_" + stepname)->Fill(mom, w);
607 mPlusEl++;
608 }
609 }
610 else if (cand.fCharge < 0) {
611 fH.fHM.HnSparse("hfsc_mom_urqmd_urqmd_minus_" + stepname)->Fill(mom, w);
612 mMinus++;
613 if (std::abs(pdg) == 11) {
614 fH.fHM.HnSparse("hfsc_mom_urqmd_urqmd_minusEl_" + stepname)->Fill(mom, w);
615 mMinusEl++;
616 }
617 }
618
619 // Particle-based
620 if (std::abs(pdg) == 11) {
621 // positrons
622 if (cand.fCharge > 0) {
623 if (cand.IsMcPi0()) {
624 fH.fHM.HnSparse("hfsp_mom_urqmd_pi0El_plus_" + stepname)->Fill(mom, w);
625 mPi0ElP++;
626 }
627 else if (cand.IsMcGamma()) {
628 fH.fHM.HnSparse("hfsp_mom_urqmd_gammaEl_plus_" + stepname)->Fill(mom, w);
629 mGammaElP++;
630 }
631 else if (cand.IsMcEta()) {
632 fH.fHM.HnSparse("hfsp_mom_urqmd_etaEl_plus_" + stepname)->Fill(mom, w);
633 mEtaElP++;
634 }
635 else {
636 fH.fHM.HnSparse("hfsp_mom_urqmd_otherEl_plus_" + stepname)->Fill(mom, w);
637 mOtherElP++;
638 }
639 }
640 // electrons
641 else if (cand.fCharge < 0) {
642 if (cand.IsMcPi0()) {
643 fH.fHM.HnSparse("hfsp_mom_urqmd_pi0El_minus_" + stepname)->Fill(mom, w);
644 mPi0ElM++;
645 }
646 else if (cand.IsMcGamma()) {
647 fH.fHM.HnSparse("hfsp_mom_urqmd_gammaEl_minus_" + stepname)->Fill(mom, w);
648 mGammaElM++;
649 }
650 else if (cand.IsMcEta()) {
651 fH.fHM.HnSparse("hfsp_mom_urqmd_etaEl_minus_" + stepname)->Fill(mom, w);
652 mEtaElM++;
653 }
654 else {
655 fH.fHM.HnSparse("hfsp_mom_urqmd_otherEl_minus_" + stepname)->Fill(mom, w);
656 mOtherElM++;
657 }
658 }
659 } // ...if (pdg == 11)
660
661 // misidentified other particles
662 else {
663 if (pdg == 211) {
664 fH.fHM.HnSparse("hfsp_mom_urqmd_pion_plus_" + stepname)->Fill(mom, w);
665 mPionP++;
666 }
667 else if (pdg == -211) {
668 fH.fHM.HnSparse("hfsp_mom_urqmd_pion_minus_" + stepname)->Fill(mom, w);
669 mPionM++;
670 }
671 else if (pdg == 2212) {
672 fH.fHM.HnSparse("hfsp_mom_urqmd_proton_plus_" + stepname)->Fill(mom, w);
673 mProton++;
674 }
675 else if (pdg == 321) {
676 fH.fHM.HnSparse("hfsp_mom_urqmd_kaon_plus_" + stepname)->Fill(mom, w);
677 mKaonP++;
678 }
679 else if (pdg == -321) {
680 fH.fHM.HnSparse("hfsp_mom_urqmd_kaon_minus_" + stepname)->Fill(mom, w);
681 mKaonM++;
682 }
683 else {
684 if (cand.fCharge > 0) {
685 fH.fHM.HnSparse("hfsp_mom_urqmd_other_plus_" + stepname)->Fill(mom, w);
686 mOtherP++;
687 }
688 else {
689 fH.fHM.HnSparse("hfsp_mom_urqmd_other_minus_" + stepname)->Fill(mom, w);
690 mOtherM++;
691 }
692 }
693 } // ... if (pdg != 11)
694 } // ... if !signal
695 } // ... for: cands
696
697 // Filling mulitplicity values for UrQMD
698 // For Pluto, no multiplicity has to be filled here, since only max. one per event exists and the first
699 // bin ('0/event') will be calculated in FastSim procedure (see there).
700
701 // charge-based
702 fH.FillH1("hfsc_mult_urqmd_urqmd_plus_" + stepname, mPlus);
703 fH.FillH1("hfsc_mult_urqmd_urqmd_minus_" + stepname, mMinus);
704 fH.FillH1("hfsc_mult_urqmd_urqmd_plusEl_" + stepname, mPlusEl);
705 fH.FillH1("hfsc_mult_urqmd_urqmd_minusEl_" + stepname, mMinusEl);
706
707 // particle-based
708 fH.FillH1("hfsp_mult_urqmd_pi0El_plus_" + stepname, mPi0ElP);
709 fH.FillH1("hfsp_mult_urqmd_pi0El_minus_" + stepname, mPi0ElM);
710 fH.FillH1("hfsp_mult_urqmd_gammaEl_plus_" + stepname, mGammaElP);
711 fH.FillH1("hfsp_mult_urqmd_gammaEl_minus_" + stepname, mGammaElM);
712 fH.FillH1("hfsp_mult_urqmd_etaEl_plus_" + stepname, mEtaElP);
713 fH.FillH1("hfsp_mult_urqmd_etaEl_minus_" + stepname, mEtaElM);
714 fH.FillH1("hfsp_mult_urqmd_otherEl_plus_" + stepname, mOtherElP);
715 fH.FillH1("hfsp_mult_urqmd_otherEl_minus_" + stepname, mOtherElM);
716 fH.FillH1("hfsp_mult_urqmd_pion_plus_" + stepname, mPionP);
717 fH.FillH1("hfsp_mult_urqmd_pion_minus_" + stepname, mPionM);
718 fH.FillH1("hfsp_mult_urqmd_proton_plus_" + stepname, mProton);
719 fH.FillH1("hfsp_mult_urqmd_kaon_plus_" + stepname, mKaonP);
720 fH.FillH1("hfsp_mult_urqmd_kaon_minus_" + stepname, mKaonM);
721 fH.FillH1("hfsp_mult_urqmd_other_plus_" + stepname, mOtherP);
722 fH.FillH1("hfsp_mult_urqmd_other_minus_" + stepname, mOtherM);
723 }
724}
725
726void LmvmTask::AnalyseRichRings()
727{
728 int nofRichRings = fRichRings->GetEntriesFast();
729 fH.FillH1("hNofRichRings", nofRichRings);
730 for (int iR = 0; iR < nofRichRings; iR++) {
731 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(iR);
732 if (richRing == nullptr) continue;
733 double rX = richRing->GetCenterX();
734 double rY = richRing->GetCenterY();
735 double dR = sqrt(rX * rX + rY * rY);
736 double aAxis = richRing->GetAaxis();
737 double bAxis = richRing->GetBaxis();
738 double BoA = bAxis / aAxis;
739 double dNext = 999.;
740 double dSecondNext = 999.;
741 string pidString = "";
742 bool isAcc = false;
743 bool isAccMin = false;
744 bool isChi2Prim = false;
745 bool isChi2PrimMin = false;
746 bool isElectron = false;
747 bool isElectronMin = false;
748 for (int iT = 0; iT < fGlobalTracks->GetEntriesFast(); iT++) {
749 CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iT));
750 if (gTrack == nullptr) continue;
751 int stsInd = gTrack->GetStsTrackIndex();
752 if (stsInd < 0) continue;
753 CbmMCTrack* mct = GetMcTrackSts(stsInd);
754
755 if (mct == nullptr) continue;
756 double p = mct->GetP();
757 const FairTrackParam* pTrack = static_cast<const FairTrackParam*>(fRichProj->At(stsInd));
758 if (pTrack == NULL) continue;
759 double dX = rX - pTrack->GetX();
760 double dY = rY - pTrack->GetY();
761 double rtDist = sqrt(dX * dX + dY * dY);
762
763 //Calculate Chi2
764 CbmL1PFFitter fPFFitter;
765 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
766 vector<CbmStsTrack> stsTracks;
767 stsTracks.resize(1);
768 stsTracks[0] = *stsTrack;
769 vector<CbmL1PFFitter::PFFieldRegion> vField;
770 vector<float> chiPrim;
771 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
772 double chi2 = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
773
774 isAcc = IsRecoTrackAccepted(gTrack);
775 isChi2Prim = (chi2 <= fCuts.fChi2PrimCut);
776 isElectron = (CbmLitGlobalElectronId::GetInstance().IsRichElectron(iT, p)
777 && CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iT, p)
778 && CbmLitGlobalElectronId::GetInstance().IsTofElectron(iT, p));
779
780 fH.FillH1("hRichRingBoA_mc", BoA); // TODO: scale Plutos?
781 if (isAcc) fH.FillH1("hRichRingBoA_acc", BoA);
782 if (isAcc && isChi2Prim) fH.FillH1("hRichRingBoA_chi2prim", BoA);
783 if (isAcc && isChi2Prim && isElectron) fH.FillH1("hRichRingBoA_elid", BoA);
784
785 if (rtDist < dNext) {
786 dSecondNext = dNext;
787 dNext = rtDist;
788 pidString = GetPidString(mct, nullptr);
789 isAccMin = isAcc;
790 isChi2PrimMin = isChi2Prim;
791 isElectronMin = isElectron;
792 }
793 }
794 if (pidString == "") continue;
795 fH.FillH2("RichRingNeighbour_all_next_" + pidString, dR, dNext);
796 fH.FillH2("RichRingNeighbour_all_secondnext_" + pidString, dR, dSecondNext);
797 if (isAccMin && isChi2PrimMin && isElectronMin) {
798 fH.FillH2("RichRingNeighbour_elid_next_" + pidString, dR, dNext);
799 fH.FillH2("RichRingNeighbour_elid_secondnext_" + pidString, dR, dSecondNext);
800 }
801 } // iR
802}
803
804void LmvmTask::FillMomHists(const CbmMCTrack* mct, const LmvmCand* cand, ELmvmSrc src, ELmvmAnaStep step)
805{
806 if ((mct != nullptr && cand != nullptr) || (mct == nullptr && cand == nullptr)) {
807 LOG(error) << "LmvmTask::FillMomHists: Both mct and cand are [not nullptr] or [nullptr].";
808 return;
809 }
810 bool isMc = (mct != nullptr);
811 string chargeStr = (isMc) ? LmvmUtils::GetChargeStr(mct) : LmvmUtils::GetChargeStr(cand);
812
813 int pdg = (isMc) ? mct->GetPdgCode() : cand->fMcPdg;
814
815 double w = 1.;
816 if (!isMc)
817 w = cand->fWeight; // changed on 10.03.25
818 else if (isMc && LmvmUtils::IsMcSignalEl(mct))
819 w = fW * MinvScale(mct, fParticle);
820
821 for (const string& suff : {string(""), chargeStr}) {
822 if (suff == "0") continue;
823 fH.FillH1("hMom" + suff, src, step, (isMc) ? mct->GetP() : cand->fMomentum.Mag(), w);
824 fH.FillH1("hMomPx" + suff, src, step, (isMc) ? mct->GetPx() : cand->fMomentum.X(), w);
825 fH.FillH1("hMomPy" + suff, src, step, (isMc) ? mct->GetPy() : cand->fMomentum.Y(), w);
826 fH.FillH1("hMomPz" + suff, src, step, (isMc) ? mct->GetPz() : cand->fMomentum.Z(), w);
827 fH.FillH1("hPt" + suff, src, step, (isMc) ? mct->GetPt() : cand->fMomentum.Perp(), w);
828 fH.FillH1("hRapidity" + suff, src, step, (isMc) ? mct->GetRapidity() : cand->fRapidity, w);
829 if ((suff == "+" || suff == "-") && src == ELmvmSrc::Bg && std::abs(pdg) == 11 && step == ELmvmAnaStep::ElId)
830 fH.FillH1("hMom" + suff + "_trueEl_bg_elid", (isMc) ? mct->GetP() : cand->fMomentum.Mag(), w);
831 }
832}
833
834void LmvmTask::DoMcTrack()
835{
836 int nMcTracks = fMCTracks->GetEntriesFast();
837
838 for (int i = 0; i < nMcTracks; i++) {
839 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i));
840 if (mct == nullptr) continue;
841 ELmvmSrc src = LmvmUtils::GetMcSrc(mct, fMCTracks);
842 string chargeStr = (mct->GetCharge() > 0) ? "+" : "-";
843 bool isSignal = LmvmUtils::IsMcSignalEl(mct);
844 double mom = mct->GetP();
845 double w = (isSignal) ? fW * MinvScale(mct, fParticle) : 1.;
846 int pdg = mct->GetPdgCode();
847
848 double binX = (std::abs(pdg) == 11) ? 0.5
849 : (std::abs(pdg) == 211) ? 1.5
850 : (pdg == 111) ? 2.5
851 : (pdg == 2212) ? 3.5
852 : (pdg == 2112) ? 4.5
853 : (pdg == 50000050) ? 5.5
854 : 6.5;
855 fH.FillH1("hNofMcTracks", binX, w);
856
857 FillMomHists(mct, nullptr, src, ELmvmAnaStep::Mc);
858
859 if (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4)
860 fH.FillH1("hMomAcc" + chargeStr + "_sts", src, mom, w);
861 if (fNofHitsInRingMap[i] >= 7) fH.FillH1("hMomAcc" + chargeStr + "_rich", src, mom, w);
862 if (mct->GetNPoints(ECbmModuleId::kTrd) >= 2) fH.FillH1("hMomAcc" + chargeStr + "_trd", src, mom, w);
863 if (mct->GetNPoints(ECbmModuleId::kTof) >= 1) fH.FillH1("hMomAcc" + chargeStr + "_tof", src, mom, w);
864
865 if (LmvmUtils::IsMcGammaEl(mct, fMCTracks)) {
866 TVector3 v;
867 mct->GetStartVertex(v);
868 fH.FillH2("hVertexGammaXZ", ELmvmAnaStep::Mc, v.Z(), v.X());
869 fH.FillH2("hVertexGammaYZ", ELmvmAnaStep::Mc, v.Z(), v.Y());
870 fH.FillH2("hVertexGammaXY", ELmvmAnaStep::Mc, v.X(), v.Y());
871 fH.FillH2("hVertexGammaRZ", ELmvmAnaStep::Mc, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
872 }
873
874 // Fill PDG histos
875 if (std::abs(pdg) == 11 || pdg == 99009911) {
876 int mcMotherPdg = 0;
877 if (mct->GetMotherId() != -1) {
878 CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMCTracks->At(mct->GetMotherId()));
879 if (mother != nullptr) mcMotherPdg = mother->GetPdgCode();
880 }
881 if (std::abs(pdg) == 11) fH.FillH1("hMotherPdg_mc", mcMotherPdg);
882 }
883
884 if (std::abs(pdg) == 211)
885 fH.FillH1("hSuppression_pion_mc", mom, w);
886 else if (pdg == 2212)
887 fH.FillH1("hSuppression_proton_mc", mom, w);
888
889 string pidString = GetPidString(mct, nullptr);
890 fH.FillH1("hMom_" + pidString, ELmvmAnaStep::Mc, mom, w);
891 fH.FillH1("hCandPdg", ELmvmAnaStep::Mc, pdg, w);
892 /*fH.FillH2("hMomPt_" + pidString, ELmvmAnaStep::Mc, mom, cand.fMomentum.Perp(), w);
893 fH.FillH2("hTrdElLike_" + pidString, ELmvmAnaStep::Mc, mom, cand.fTrdLikeEl, w); use here mcTrack, not cand!
894 fH.FillH2("hELossSts_" + pidString, ELmvmAnaStep::Mc, mom, cand.fELossSts, w);
895 fH.FillH2("hRichAnn_" + pidString, ELmvmAnaStep::Mc, mom, cand.fRichAnn, w);
896 fH.FillH2("hChi2Dets_sts_" + pidString, ELmvmAnaStep::Mc, mom, cand.fChi2Sts, w);
897 fH.FillH2("hChi2Dets_rich_" + pidString, ELmvmAnaStep::Mc, mom, cand.fChi2Rich, w);
898 fH.FillH2("hChi2Dets_trd_" + pidString, ELmvmAnaStep::Mc, pRec, cand.fChi2Trd, w);
899 fH.FillH2("hTofM2_" + pidString, ELmvmAnaStep::Mc, pRec, cand.fMass2, w);*/
900 }
901}
902
903void LmvmTask::DoMcPair()
904{
905 int nMcTracks = fMCTracks->GetEntries();
906 for (int iMc1 = 0; iMc1 < nMcTracks - 1; iMc1++) {
907 CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc1));
908 ELmvmSrc src1 = LmvmUtils::GetMcSrc(mct1, fMCTracks);
909
910 // To speed up: select only signal, eta and pi0 electrons
911 if (!(src1 == ELmvmSrc::Signal || src1 == ELmvmSrc::Pi0 || src1 == ELmvmSrc::Eta)) continue;
912
913
914 for (int iMc2 = iMc1 + 1; iMc2 < nMcTracks; iMc2++) {
915 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2));
916 ELmvmSrc src2 = LmvmUtils::GetMcSrc(mct2, fMCTracks);
917 ELmvmSrc srcPair = LmvmUtils::GetMcPairSrc(mct1, mct2, fMCTracks);
918 LmvmKinePar pKin = LmvmKinePar::Create(mct1, mct2);
919 double w = (src1 == ELmvmSrc::Signal || src2 == ELmvmSrc::Signal)
920 ? fW * LmvmUtils::MinvScale(pKin.fMinv, fParticle)
921 : 1.; // changed on 10.03.25
922 if (srcPair == ELmvmSrc::Signal) {
923 fH.FillH2("hMomVsAnglePairSignalMc", std::sqrt(mct1->GetP() * mct2->GetP()), pKin.fAngle, w);
924 fH.FillH2("hPtYPairSignal", ELmvmAnaStep::Mc, pKin.fRapidity, pKin.fPt, w);
925 fH.FillH1("hMomPairSignal", ELmvmAnaStep::Mc, pKin.fMomentumMag, w);
926 fH.FillH1("hMinv_unscaled", pKin.fMinv);
927 }
928 fH.FillH1("hMinv", srcPair, ELmvmAnaStep::Mc, pKin.fMinv, w);
929 }
930 }
931}
932
933void LmvmTask::RichPmtXY()
934{
935 int nofRichHits = fRichHits->GetEntriesFast();
936 for (int iH = 0; iH < nofRichHits; iH++) {
937 CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
938 if (richHit == nullptr || richHit->GetRefId() < 0) continue;
939
940 Int_t digiIndex = richHit->GetRefId();
941 if (digiIndex < 0) continue;
942 const CbmRichDigi* digi = fDigiManager->Get<CbmRichDigi>(digiIndex);
943 if (!digi) continue;
944 const CbmMatch* digiMatch = fDigiManager->GetMatch(ECbmModuleId::kRich, digiIndex);
945 if (!digiMatch) continue;
946 vector<CbmLink> links = digiMatch->GetLinks();
947 for (const auto& link : links) {
948 if (link.IsNoise()) continue;
949 FairMCPoint* pointPhoton = static_cast<FairMCPoint*>(fRichPoints->At(link.GetIndex()));
950 if (!pointPhoton) continue;
951 CbmMCTrack* trackPhoton = static_cast<CbmMCTrack*>(fMCTracks->At(pointPhoton->GetTrackID()));
952 if (!trackPhoton || trackPhoton->GetMotherId() < 0) continue;
953 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(trackPhoton->GetMotherId()));
954 if (mct == nullptr) continue;
955 TVector3 v;
956 mct->GetStartVertex(v);
957 ELmvmSrc src = LmvmUtils::GetMcSrc(mct, fMCTracks);
958 double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW * MinvScale(mct, fParticle) : 1.;
959 if (IsPrimary(v.Z())) {
960 fH.FillH2("hPmtXY", src, richHit->GetX(), richHit->GetY(), w);
961 }
962 }
963 }
964}
965
966bool LmvmTask::IsRecoTrackAccepted(const CbmGlobalTrack* gTrack)
967{
968 if (gTrack->GetStsTrackIndex() < 0 || gTrack->GetRichRingIndex() < 0 || gTrack->GetTrdTrackIndex() < 0
969 || gTrack->GetTofHitIndex() < 0)
970 return false;
971 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(gTrack->GetStsTrackIndex());
972 if (stsTrack == nullptr) return false;
973 int nStsHits = stsTrack->GetNofStsHits();
974 int nMvdHits = stsTrack->GetNofMvdHits();
975
976 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(gTrack->GetRichRingIndex());
977 if (richRing == nullptr) return false;
978 int nRichHits = richRing->GetNofHits();
979
980 CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTracks->At(gTrack->GetTrdTrackIndex());
981 if (trdTrack == nullptr) return false;
982 int nTrdHits = trdTrack->GetNofHits();
983
984 CbmTofTrack* tofTrack = (CbmTofTrack*) fTofTracks->At(gTrack->GetTofTrackIndex());
985 if (tofTrack == nullptr) return false;
986 int nTofHits = tofTrack->GetNofTofHits();
987
988 if (nStsHits + nMvdHits >= fNofMinHitsStsMvd && nRichHits >= fNofMinHitsRich && nTrdHits >= fNofMinHitsTrd
989 && nTofHits >= fNofMinHitsTof)
990 return true;
991 return false;
992}
993
994void LmvmTask::AnalyseGlobalTracks()
995{
996 int ngTracks = fGlobalTracks->GetEntriesFast();
997
998 for (int i = 0; i < ngTracks; i++) {
999 CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(i));
1000 if (gTrack == nullptr) continue;
1001 int stsInd = gTrack->GetStsTrackIndex();
1002 bool isRich = (gTrack->GetRichRingIndex() >= 0);
1003 bool isTrd = (gTrack->GetTrdTrackIndex() >= 0);
1004 bool isTof = (gTrack->GetTofHitIndex() >= 0);
1005
1006 if (stsInd < 0) continue;
1007 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
1008 if (stsTrack == nullptr) continue;
1009 CbmMCTrack* mct = GetMcTrackSts(stsInd);
1010 if (mct == nullptr) continue;
1011
1012 CbmRichRing* richRing = nullptr;
1013 CbmTrdTrack* trdTrack = nullptr;
1014 CbmTofHit* tofHit = nullptr;
1015
1016 if (gTrack->GetRichRingIndex() >= 0) richRing = (CbmRichRing*) fRichRings->At(gTrack->GetRichRingIndex());
1017 if (gTrack->GetTrdTrackIndex() >= 0) trdTrack = (CbmTrdTrack*) fTrdTracks->At(gTrack->GetTrdTrackIndex());
1018 if (gTrack->GetTofHitIndex() >= 0) tofHit = (CbmTofHit*) fTofHits->At(gTrack->GetTofHitIndex());
1019
1020 double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW * MinvScale(mct, fParticle) : 1.;
1021
1022 bool isSignalEl = LmvmUtils::IsMcSignalEl(mct);
1023 bool isGammaEl = LmvmUtils::IsMcGammaEl(mct, fMCTracks);
1024 int pdg = mct->GetPdgCode();
1025 double p = mct->GetP();
1026 double pt = mct->GetPt();
1027 double rap = mct->GetRapidity();
1028 double m2 = CbmLitGlobalElectronId::GetInstance().GetTofM2(i, p);
1029 bool isAcc = IsRecoTrackAccepted(gTrack);
1030
1031 double binX = (std::abs(pdg) == 11 && isSignalEl) ? 0.5
1032 : (std::abs(pdg) == 11 && isGammaEl) ? 1.5
1033 : (std::abs(pdg) == 11 && !isSignalEl && !isGammaEl) ? 2.5
1034 : (std::abs(pdg) == 211) ? 3.5
1035 : (pdg == 2212) ? 4.5
1036 : 5.5;
1037 fH.FillH1("hNofGTracks", binX, w);
1038
1039 TVector3 v;
1040 mct->GetStartVertex(v);
1041 string ptcl = GetPidString(v.Z(), pdg);
1042 bool isPrim = IsPrimary(v.Z());
1043
1044 if (std::abs(pdg) == 11 || std::abs(pdg) == 211 || pdg == 2212 || std::abs(pdg) == 321) {
1045 string hName = (std::abs(pdg) == 11) ? "hEl" : (std::abs(pdg) == 211) ? "hPi" : pdg == 2212 ? "hProton" : "hKaon";
1046 fH.FillH1(hName + "Mom_all_recSts", p);
1047 if (isRich) fH.FillH1(hName + "Mom_all_recStsRich", p);
1048 if (isRich && isTrd) fH.FillH1(hName + "Mom_all_recStsRichTrd", p);
1049 if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_all_recStsRichTrdTof", p);
1050
1051 if (isPrim) {
1052 fH.FillH1(hName + "Mom_prim_recSts", p);
1053 if (isRich) fH.FillH1(hName + "Mom_prim_recStsRich", p);
1054 if (isRich && isTrd) fH.FillH1(hName + "Mom_prim_recStsRichTrd", p);
1055 if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_prim_recStsRichTrdTof", p);
1056 }
1057 }
1058
1059 // Fill NofHits in diff. detectors
1060 int nRichHits = -1;
1061 int nTrdHits = -1;
1062 int nTofHits = -1;
1063
1064 int nStsHits = stsTrack->GetNofStsHits();
1065 if (richRing != nullptr) nRichHits = richRing->GetNofHits();
1066 if (trdTrack != nullptr) nTrdHits = trdTrack->GetNofHits();
1067 if (tofHit != nullptr) nTofHits = 1;
1068
1069 for (const string det : {"sts", "rich", "trd", "tof"}) {
1070 int nHits = (det == "sts") ? nStsHits : (det == "rich") ? nRichHits : (det == "trd") ? nTrdHits : nTofHits;
1071 fH.FillH2("hNofHitsVsMom_" + det + "_" + ptcl, p, nHits, w);
1072 }
1073
1074 fH.FillH2("hVertexGTrackRZ", v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
1075
1076 bool isElectron = (CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, p)
1077 && CbmLitGlobalElectronId::GetInstance().IsTrdElectron(i, p)
1078 && CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p))
1079 ? true
1080 : false;
1081
1082 CheckMismatches(gTrack, pdg, isElectron, ptcl, w);
1083 if (std::abs(pdg) != 11) PidVsMom(gTrack, i, pdg, p, isAcc);
1084
1085 // investigate misidentifications: compare misidentified candidates with same not-misidentified particles
1086 // first check if global track has entries in all detectors
1087 if (!IsInAllDets(gTrack)) continue; // Mind: all following actions are done only for fully rec. gTracks
1088
1089 string ptcl2 = GetPidString(mct, nullptr);
1090 if (!isElectron && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) {
1091 fH.FillH1("hMom_" + ptcl2 + "_true", p);
1092 fH.FillH2("hPtY_" + ptcl2 + "_true", rap, pt);
1093 fH.FillH2("hTofM2_" + ptcl2 + "_true", p, m2);
1094 }
1095 else if (isElectron
1096 && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) {
1097 fH.FillH1("hMom_" + ptcl2 + "_misid", p);
1098 fH.FillH2("hPtY_" + ptcl2 + "_misid", rap, pt);
1099 fH.FillH2("hTofM2_" + ptcl2 + "_misid", p, m2);
1100 }
1101 }
1102}
1103
1104void LmvmTask::PidVsMom(const CbmGlobalTrack* gTrack, int iGTrack, int pdg, double p, bool isAcc)
1105{
1106 // Calculate Chi2
1107 CbmL1PFFitter fPFFitter;
1108 int stsInd = gTrack->GetStsTrackIndex();
1109 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
1110 vector<CbmStsTrack> stsTracks;
1111 stsTracks.resize(1);
1112 stsTracks[0] = *stsTrack;
1113 vector<CbmL1PFFitter::PFFieldRegion> vField;
1114 vector<float> chiPrim;
1115 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
1116 double chi2 = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
1117 bool isChi2 = fCuts.IsChi2PrimaryOk(chi2);
1118
1119 bool isRichEl = CbmLitGlobalElectronId::GetInstance().IsRichElectron(iGTrack, p);
1120 bool isTrdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iGTrack, p);
1121 bool isTofEl = CbmLitGlobalElectronId::GetInstance().IsTofElectron(iGTrack, p);
1122
1123 double pdgBin = (pdg == 211) ? 0.5
1124 : (pdg == -211) ? 1.5
1125 : (pdg == 2212) ? 2.5
1126 : (pdg == 321) ? 3.5
1127 : (pdg == -321) ? 4.5
1128 : 5.5;
1129
1130 if (isRichEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_rich_rec", p, pdgBin);
1131 if (isTrdEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_trd_rec", p, pdgBin);
1132 if (isTofEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_tof_rec", p, pdgBin);
1133
1134 if (isRichEl && isAcc) fH.FillH2("hPdgVsMom_rich_acc", p, pdgBin);
1135 if (isTrdEl && isAcc) fH.FillH2("hPdgVsMom_trd_acc", p, pdgBin);
1136 if (isTofEl && isAcc) fH.FillH2("hPdgVsMom_tof_acc", p, pdgBin);
1137
1138 if (isRichEl && isAcc && isChi2) fH.FillH2("hPdgVsMom_rich_chi2prim", p, pdgBin);
1139 if (isTrdEl && isAcc && isChi2) fH.FillH2("hPdgVsMom_trd_chi2prim", p, pdgBin);
1140 if (isTofEl && isAcc && isChi2) fH.FillH2("hPdgVsMom_tof_chi2prim", p, pdgBin);
1141}
1142
1143void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isElectron, const string& ptcl, double w)
1144{
1145 // Chi2 of true-/mismatched; number of mismatches (general and seperated for detectors)
1146 int stsInd = gTrack->GetStsTrackIndex();
1147 if (stsInd < 0) return;
1148 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
1149 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) return;
1150 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1151 if (stsMcTrackId < 0) return;
1152
1153 bool isFull = IsInAllDets(gTrack);
1154
1155 // first calculate Chi2
1156 CbmL1PFFitter fPFFitter;
1157 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
1158 vector<CbmStsTrack> stsTracks;
1159 stsTracks.resize(1);
1160 stsTracks[0] = *stsTrack;
1161 vector<CbmL1PFFitter::PFFieldRegion> vField;
1162 vector<float> chiPrim;
1163 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
1164 double chi2 = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
1165
1166 CbmTrackMatchNew* richMatch = nullptr;
1167 CbmTrackMatchNew* trdMatch = nullptr;
1168 CbmMatch* tofMatch = nullptr;
1169
1170 int richRingInd = gTrack->GetRichRingIndex();
1171 int trdTrackInd = gTrack->GetTrdTrackIndex();
1172 int tofHitInd = gTrack->GetTofHitIndex();
1173
1174 if (richRingInd >= 0) richMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(richRingInd));
1175 if (trdTrackInd >= 0) trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(trdTrackInd));
1176 if (tofHitInd >= 0) tofMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofHitInd));
1177
1178 int nofMismTrackSegmentsAll = 0;
1179 int nofMismTrackSegmentsComp = 0;
1180
1181 if (isFull) fH.FillH1("hNofMismatches_gTracks_" + ptcl, 0.5, w);
1182
1183 if (richMatch != nullptr && richMatch->GetNofLinks() > 0) {
1184 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
1185 if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) {
1186 fH.FillH1("hMatch_rich", 1.5, w);
1187 nofMismTrackSegmentsAll++;
1188 }
1189 else if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId)
1190 fH.FillH1("hMatch_rich", 0.5, w);
1191 if (isFull) {
1192 if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_" + ptcl, 1.5, w);
1193 if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_rich_" + ptcl, chi2, w);
1194 if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) {
1195 fH.FillH1("hChi2_mismatch_rich_" + ptcl, chi2, w);
1196 fH.FillH1("hNofMismatches_gTracks_" + ptcl, 2.5, w);
1197 nofMismTrackSegmentsComp++;
1198 }
1199 }
1200 }
1201
1202 if (trdMatch != nullptr && trdMatch->GetNofLinks() > 0) {
1203 int trdMcTrackId = trdMatch->GetMatchedLink().GetIndex();
1204 if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) {
1205 fH.FillH1("hMatch_trd", 1.5, w);
1206 nofMismTrackSegmentsAll++;
1207 }
1208 else if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId)
1209 fH.FillH1("hMatch_trd", 0.5, w);
1210 if (isFull) {
1211 if (trdMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_" + ptcl, 3.5, w);
1212 if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_trd_" + ptcl, chi2, w);
1213 if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) {
1214 fH.FillH1("hNofMismatches_gTracks_" + ptcl, 4.5, w);
1215 fH.FillH1("hChi2_mismatch_trd_" + ptcl, chi2, w);
1216 nofMismTrackSegmentsComp++;
1217 }
1218 }
1219 }
1220
1221 if (tofMatch != nullptr && tofMatch->GetMatchedLink().GetIndex() >= 0) {
1222 FairMCPoint* tofPoint = nullptr;
1223 tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofMatch->GetMatchedLink().GetIndex()));
1224 if (tofPoint != nullptr) {
1225 int tofMcTrackId = tofPoint->GetTrackID();
1226 if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) {
1227 fH.FillH1("hMatch_tof", 1.5, w);
1228 nofMismTrackSegmentsAll++;
1229 }
1230 else if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId)
1231 fH.FillH1("hMatch_tof", 0.5, w);
1232 if (isFull) {
1233 if (tofMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_" + ptcl, 5.5, w);
1234 if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_tof_" + ptcl, chi2, w);
1235 if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) {
1236 fH.FillH1("hNofMismatches_gTracks_" + ptcl, 6.5, w);
1237 fH.FillH1("hChi2_mismatch_tof_" + ptcl, chi2, w);
1238 nofMismTrackSegmentsComp++;
1239 }
1240 }
1241 }
1242 }
1243
1244 if (isFull) fH.FillH1("hNofMismatchedTrackSegments_" + ptcl, nofMismTrackSegmentsComp + 0.5, w);
1245
1246 bool isTrueId = (std::abs(pdg) == 11 && isElectron) ? true : (std::abs(pdg) != 11 && !isElectron) ? true : false;
1247 double binX = nofMismTrackSegmentsComp + 0.5;
1248 double binY = (isTrueId == true) ? 0.5 : 1.5;
1249
1250 fH.FillH2("hMatchId_gTracks_" + ptcl, binX, binY, w);
1251}
1252
1253bool LmvmTask::IsInAllDets(const CbmGlobalTrack* gTrack)
1254{
1255 if (gTrack == nullptr) return false;
1256
1257 int stsInd = gTrack->GetStsTrackIndex();
1258 if (stsInd < 0) return false;
1259 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
1260 if (stsTrack == nullptr) return false;
1261
1262 int richInd = gTrack->GetRichRingIndex();
1263 int trdInd = gTrack->GetTrdTrackIndex();
1264 int tofInd = gTrack->GetTofHitIndex();
1265 if (richInd < 0 || trdInd < 0 || tofInd < 0) return false;
1266
1267 CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(richInd));
1268 CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(trdInd));
1269 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(tofInd));
1270 if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr) return false;
1271
1272 return true;
1273}
1274
1275void LmvmTask::FillTopologyCands()
1276{
1277 fSTCands.clear();
1278 fRTCands.clear();
1279 int ngTracks = fGlobalTracks->GetEntriesFast();
1280
1281 for (int iGTrack = 0; iGTrack < ngTracks; iGTrack++) {
1282 LmvmCand cand;
1283
1284 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGTrack);
1285 if (gTrack == nullptr) continue;
1286
1287 cand.fStsInd = gTrack->GetStsTrackIndex();
1288 if (cand.fStsInd < 0) continue;
1289 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
1290 if (stsTrack == nullptr) continue;
1291
1292 cand.fRichInd = gTrack->GetRichRingIndex();
1293 cand.fTrdInd = gTrack->GetTrdTrackIndex();
1294 cand.fTofHitInd = gTrack->GetTofHitIndex();
1295
1296 LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex);
1297 if (cand.fCharge == 0) {
1298 LOG(fatal) << "LmvmTask::FillTopologyCands: Charge of Global Track is 0. Check!";
1299 continue;
1300 }
1301 cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim);
1302 if (!cand.fIsChi2Prim) continue;
1303
1304 // ST cut candidates, only STS
1305 if (cand.fRichInd < 0 && cand.fTrdInd < 0 && cand.fTofHitInd < 0) fSTCands.push_back(cand);
1306
1307 // RT cut candidates, STS + at least one detector (RICH, TRD, TOF) but not all
1308 // Candidates must be identified as electron in registered detectors:
1309 // if it is registered in RICH it must be identified in the RICH as electron
1310 // RICH
1311 bool isRichRT = (cand.fRichInd < 0) ? false : true;
1312 if (isRichRT) {
1313 CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(cand.fRichInd));
1314 if (richRing == nullptr) isRichRT = false;
1315 if (isRichRT) isRichRT = CbmLitGlobalElectronId::GetInstance().IsRichElectron(iGTrack, cand.fMomentum.Mag());
1316 }
1317
1318 // TRD
1319 bool isTrdRT = (cand.fTrdInd < 0) ? false : true;
1320 if (isTrdRT) {
1321 CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(cand.fTrdInd));
1322 if (trdTrack == nullptr) isTrdRT = false;
1323 if (isTrdRT) isTrdRT = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iGTrack, cand.fMomentum.Mag());
1324 }
1325
1326
1327 // ToF
1328 bool isTofRT = (cand.fTofHitInd < 0) ? false : true;
1329 if (isTofRT) {
1330 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofHitInd));
1331 if (tofHit == nullptr) isTofRT = false;
1332 if (isTofRT) isTofRT = CbmLitGlobalElectronId::GetInstance().IsTofElectron(iGTrack, cand.fMomentum.Mag());
1333 }
1334
1335 if (isRichRT || isTrdRT || isTofRT) {
1336 if (!(cand.fRichInd >= 0 && cand.fTrdInd >= 0 && cand.fTofHitInd >= 0)) {
1337 fRTCands.push_back(cand);
1338 }
1339 }
1340 }
1341 LOG(info) << "fSTCands.size = " << fSTCands.size();
1342 LOG(info) << "fRTCands.size = " << fRTCands.size();
1343
1344 AssignMcToTopologyCands(fSTCands);
1345 AssignMcToTopologyCands(fRTCands);
1346}
1347
1348void LmvmTask::FillCands()
1349{
1350 fCands.clear();
1351 fTTCands.clear();
1352 int nGTracks = fGlobalTracks->GetEntriesFast();
1353 fCands.reserve(nGTracks);
1354
1355 for (int iGTrack = 0; iGTrack < nGTracks; iGTrack++) {
1356 LmvmCand cand;
1357 cand.fEventNumber = fEventNumber;
1358 cand.fTaskId = fTaskId;
1359
1360 // if MVD is not used set mvd cuts to true
1361 cand.fIsMvd1Cut = !fUseMvd;
1362 cand.fIsMvd2Cut = !fUseMvd;
1363
1364 CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iGTrack));
1365 if (gTrack == nullptr) continue;
1366
1367 // STS
1368 cand.fStsInd = gTrack->GetStsTrackIndex();
1369 if (cand.fStsInd < 0) continue;
1370 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
1371 if (stsTrack == nullptr) continue;
1372
1373 CbmMCTrack* mct = GetMcTrackSts(cand.fStsInd);
1374 if (mct == nullptr) continue;
1375
1376 // RICH - TRD - TOF
1377 cand.fRichInd = gTrack->GetRichRingIndex();
1378 cand.fTrdInd = gTrack->GetTrdTrackIndex();
1379 cand.fTofHitInd = gTrack->GetTofHitIndex();
1380 cand.fTofTrackInd = gTrack->GetTofTrackIndex();
1381 cand.fIsRec = true;
1382 if (cand.fRichInd < 0 || cand.fTrdInd < 0 || cand.fTofHitInd < 0) continue;
1383
1384 CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(cand.fRichInd));
1385 CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(cand.fTrdInd));
1386 CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(cand.fTofTrackInd));
1387 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofHitInd));
1388
1389 if (richRing == nullptr || trdTrack == nullptr || tofTrack == nullptr) continue;
1390 LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex);
1391 if (cand.fCharge == 0) {
1392 LOG(fatal) << "LmvmTask::FillCands: Charge of Global Track is 0. Check!";
1393 continue;
1394 }
1395 cand.fGTrackInd = iGTrack;
1396 cand.fIsAcc = IsRecoTrackAccepted(gTrack);
1397 cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim);
1398 cand.fIsPtCut = fCuts.IsPtCutOk(cand.fMomentum.Mag(), cand.fMomentum.Perp());
1399 cand.fNofHitsSts = stsTrack->GetNofStsHits();
1400 cand.fNofHitsMvd = stsTrack->GetNofMvdHits();
1401 cand.fELossSts = stsTrack->GetELoss();
1402
1403 if (richRing != nullptr) cand.fNofHitsRich = richRing->GetNofHits(); // TODO: if-conditions redundant
1404 if (trdTrack != nullptr) cand.fNofHitsTrd = trdTrack->GetNofHits();
1405 if (tofTrack != nullptr) cand.fNofHitsTof = tofTrack->GetNofTofHits();
1406
1407 cand.fChi2Rich = (richRing->GetNDF() > 0.) ? richRing->GetChi2() / richRing->GetNDF() : 0.;
1408 cand.fChi2Trd = (trdTrack->GetNDF() > 0.) ? trdTrack->GetChiSq() / trdTrack->GetNDF() : 0.;
1409 cand.fChi2Tof = tofTrack->GetChiSq();
1410 cand.fTofDist = tofTrack->GetDistance();
1411 cand.fTrdLikeEl = trdTrack->GetPidLikeEL();
1412 cand.fTrdLikePi = trdTrack->GetPidLikePI();
1413
1414 LmvmUtils::IsElectron(iGTrack, cand.fMomentum.Mag(), fCuts.fMomentumCut, &cand);
1415 LmvmUtils::IsRichElectron(iGTrack, cand.fMomentum.Mag(), &cand);
1416 LmvmUtils::IsTrdElectron(iGTrack, cand.fMomentum.Mag(), &cand);
1417 LmvmUtils::IsTofElectron(iGTrack, cand.fMomentum.Mag(), &cand);
1418
1419 // Add all pions from STS for pion misidentification level study
1420 if (fPionMisidLevel >= 0.0 && mct != nullptr) {
1421 //check that pion track has track projection onto the photodetector plane
1422 const FairTrackParam* richProj = static_cast<FairTrackParam*>(fRichProj->At(iGTrack));
1423 if (richProj == nullptr || richProj->GetX() == 0 || richProj->GetY() == 0) continue;
1424
1425 if (std::abs(mct->GetPdgCode()) == 211) {
1426 LmvmUtils::IsElectronMc(&cand, fMCTracks, fPionMisidLevel);
1427 fCands.push_back(cand);
1428 continue;
1429 }
1430 }
1431
1432 // Set length and time
1433 cand.fLength = gTrack->GetLength() / 100.;
1434 // TODO: get time via (CbmEvent*) cbmevent->GetTzero();
1435 double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime();
1436 Double_t noOffsetTime = tofHit->GetTime() - eventTime;
1437 cand.fTime = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m
1438
1439 fCands.push_back(cand);
1440 if (!cand.fIsElectron && cand.fIsChi2Prim) fTTCands.push_back(cand);
1441 }
1442 LOG(info) << "fTTCands.size = " << fTTCands.size();
1443 LOG(info) << "fCands.size = " << fCands.size();
1444
1445 AssignMcToCands(fCands);
1446 AssignMcToTopologyCands(fTTCands);
1447}
1448
1449void LmvmTask::CombinatorialPairs()
1450{
1451 /************************************************************
1452* This function is implemented for comparing yields from *
1453* event mixing procedure LmvmEventMix with yields from this *
1454* analysis. See Histogram 'CheckWithSim' in LmvmDrawAll. *
1455************************************************************/
1456
1457 size_t nCands = fCandsTotal.size();
1458 for (size_t iC1 = 0; iC1 < nCands - 1; iC1++) {
1459 const auto& cand1 = fCandsTotal[iC1];
1460 for (size_t iC2 = iC1 + 1; iC2 < nCands; iC2++) {
1461 const auto& cand2 = fCandsTotal[iC2];
1462 string comb = "";
1463 if (cand1.fCharge * cand2.fCharge < 0)
1464 comb = "PM";
1465 else if (cand1.fCharge > 0 && cand2.fCharge > 0)
1466 comb = "PP";
1467 else if (cand1.fCharge < 0 && cand2.fCharge < 0)
1468 comb = "MM";
1469
1470 int pdg1 = std::abs(cand1.fMcPdg);
1471 int pdg2 = std::abs(cand2.fMcPdg);
1472 LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2);
1473 bool isSameEvent = (cand1.fEventNumber == cand2.fEventNumber);
1474 double w = (isSameEvent) ? LmvmUtils::GetWeightPair(cand1, cand2) : cand1.fWeight * cand2.fWeight;
1475
1476 for (auto step : fH.fAnaSteps) {
1477 if (step <= ELmvmAnaStep::Reco) continue;
1478 if (!(cand1.IsCutTill(step) && cand2.IsCutTill(step))) continue;
1479
1480 // Fill minv histos for all combinations of e+e- to compare LmvmEventMix yields with MC
1481 if (step == ELmvmAnaStep::ElId) {
1482 if (isSameEvent)
1483 fH.FillH1("hMinv" + comb + "_sameEv", pRec.fMinv, w);
1484 else
1485 fH.FillH1("hMinv" + comb + "_mixedEv", pRec.fMinv, w);
1486 }
1487
1488 // Check symmetries of misidentifications
1489 if (step == ELmvmAnaStep::ElId && isSameEvent && cand1.fCharge * cand2.fCharge < 0) {
1490 bool isC1El = (pdg1 == 11);
1491 bool isC2El = (pdg2 == 11);
1492 bool isC1Plus = (cand1.fCharge > 0);
1493 bool isC2Plus = (cand2.fCharge > 0);
1494
1495 bool cand1PlusTrue = (isC1Plus && isC1El);
1496 bool cand1MinusTrue = (!isC1Plus && isC1El);
1497 bool cand1PlusMis = (isC1Plus && !isC1El);
1498 bool cand1MinusMis = (!isC1Plus && !isC1El);
1499
1500 bool cand2PlusTrue = (isC2Plus && isC2El);
1501 bool cand2MinusTrue = (!isC2Plus && isC2El);
1502 bool cand2PlusMis = (isC2Plus && !isC2El);
1503 bool cand2MinusMis = (!isC2Plus && !isC2El);
1504
1505 if (isC1El && isC2El)
1506 fH.FillH1("hMinvCombPM_same_tt", pRec.fMinv, w);
1507 else if (!isC1El && !isC2El)
1508 fH.FillH1("hMinvCombPM_same_mm", pRec.fMinv, w);
1509 else if ((cand1PlusTrue && cand2MinusMis) || (cand1MinusMis && cand2PlusTrue))
1510 fH.FillH1("hMinvCombPM_same_tm", pRec.fMinv, w);
1511 else if ((cand1PlusMis && cand2MinusTrue) || (cand1MinusTrue && cand2PlusMis))
1512 fH.FillH1("hMinvCombPM_same_mt", pRec.fMinv, w);
1513 }
1514
1515 // Histograms for k factor in dependance on various quantities
1516 if (step == ELmvmAnaStep::ElId && !isSameEvent) {
1517 CbmMCTrack* mct1 = GetMcTrackSts(cand1.fStsInd);
1518 CbmMCTrack* mct2 = GetMcTrackSts(cand2.fStsInd);
1519 if (mct1 == nullptr || mct2 == nullptr || mct1->GetMotherId() < 0 || mct2->GetMotherId() < 0) continue;
1520 CbmMCTrack* mother1 = static_cast<CbmMCTrack*>(fMCTracks->At(mct1->GetMotherId()));
1521 CbmMCTrack* mother2 = static_cast<CbmMCTrack*>(fMCTracks->At(mct2->GetMotherId()));
1522 TVector3 v1, v2, v1Mother, v2Mother;
1523 mct1->GetStartVertex(v1);
1524 mct2->GetStartVertex(v2);
1525 double theta = pRec.fAngle;
1526 double R1 = sqrt(v1.X() * v1.X() + v1.Y() * v1.Y());
1527 double R2 = sqrt(v2.X() * v2.X() + v2.Y() * v2.Y());
1528 double R1mother = 0.;
1529 double R2mother = 0.;
1530 if (mother1 != nullptr && mother2 != nullptr) {
1531 mother1->GetStartVertex(v1Mother);
1532 mother2->GetStartVertex(v2Mother);
1533 R1mother = sqrt(v1Mother.X() * v1Mother.X() + v1Mother.Y() * v1Mother.Y());
1534 R2mother = sqrt(v2Mother.X() * v2Mother.X() + v2Mother.Y() * v2Mother.Y());
1535 }
1536
1537 if (cand1.fCharge * cand2.fCharge < 0) {
1538 fH.FillH2("hMinvCombPM_px", cand1.fMomentum.X(), cand2.fMomentum.X(), w);
1539 fH.FillH2("hMinvCombPM_py", cand1.fMomentum.Y(), cand2.fMomentum.Y(), w);
1540 fH.FillH2("hMinvCombPM_R", R1, R2, w);
1541 fH.FillH2("hMinvCombPM_theta-minv", pRec.fMinv, theta, w);
1542 if (mother1 != nullptr && mother2 != nullptr) {
1543 fH.FillH2("hMinvCombPM_px_mother", mother1->GetPx(), mother2->GetPx(), w);
1544 fH.FillH2("hMinvCombPM_py_mother", mother1->GetPy(), mother2->GetPy(), w);
1545 fH.FillH2("hMinvCombPM_R_mother", R1mother, R2mother, w);
1546 }
1547 }
1548 else if (cand1.fCharge > 0 && cand2.fCharge > 0) {
1549 fH.FillH2("hMinvCombPP_px", cand1.fMomentum.X(), cand2.fMomentum.X(), w);
1550 fH.FillH2("hMinvCombPP_py", cand1.fMomentum.Y(), cand2.fMomentum.Y(), w);
1551 fH.FillH2("hMinvCombPP_R", R1, R2, w);
1552 fH.FillH2("hMinvCombPP_theta-minv", pRec.fMinv, theta, w);
1553 if (mother1 != nullptr && mother2 != nullptr) {
1554 fH.FillH2("hMinvCombPP_px_mother", mother1->GetPx(), mother2->GetPx(), w);
1555 fH.FillH2("hMinvCombPP_py_mother", mother1->GetPy(), mother2->GetPy(), w);
1556 fH.FillH2("hMinvCombPP_R_mother", R1mother, R2mother, w);
1557 }
1558 }
1559 else if (cand1.fCharge < 0 && cand2.fCharge < 0) {
1560 fH.FillH2("hMinvCombMM_px", cand1.fMomentum.X(), cand2.fMomentum.X(), w);
1561 fH.FillH2("hMinvCombMM_py", cand1.fMomentum.Y(), cand2.fMomentum.Y(), w);
1562 fH.FillH2("hMinvCombMM_R", R1, R2, w);
1563 fH.FillH2("hMinvCombMM_theta-minv", pRec.fMinv, theta, w);
1564 if (mother1 != nullptr && mother2 != nullptr) {
1565 fH.FillH2("hMinvCombMM_px_mother", mother1->GetPx(), mother2->GetPx(), w);
1566 fH.FillH2("hMinvCombMM_py_mother", mother1->GetPy(), mother2->GetPy(), w);
1567 fH.FillH2("hMinvCombMM_R_mother", R1mother, R2mother, w);
1568 }
1569 }
1570 } // k Factor dependencies
1571 } // step
1572 } // cand2
1573 } // cand1
1574}
1575
1576void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands)
1577{
1578 for (auto& cand : cands) {
1579 cand.ResetMcParams();
1580
1581 //STS. MCTrackId of the candidate is defined by STS track
1582 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(cand.fStsInd));
1583 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
1584 cand.fStsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1585 if (cand.fStsMcTrackId < 0) continue;
1586 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1587 if (mct == nullptr) continue;
1588
1589 cand.fMcMotherId = mct->GetMotherId();
1590 cand.fMcPdg = mct->GetPdgCode();
1591 cand.fMcSrc = LmvmUtils::GetMcSrc(mct, fMCTracks);
1592
1593 // Get mass of mother particle if it is signal (for mass-dependant scaling)
1594 if (cand.IsMcSignal()) {
1595 int nMcTracks = fMCTracks->GetEntries();
1596 for (int iMc2 = 0; iMc2 < nMcTracks; iMc2++) {
1597 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2));
1598 if (mct2->GetMotherId() == cand.fMcMotherId && iMc2 != cand.fStsMcTrackId) {
1599 LmvmKinePar pKin = LmvmKinePar::Create(mct, mct2);
1600 cand.fMassSig = pKin.fMinv;
1601 cand.fWeight = fW * LmvmUtils::MinvScale(pKin.fMinv, fParticle);
1602 }
1603 }
1604 }
1605 else
1606 cand.fWeight = 1.;
1607 cand.fName = fParticle;
1608
1609 if (std::abs(cand.fMcPdg) == 211 && fPionMisidLevel >= 0.) continue;
1610
1611 // RICH
1612 CbmTrackMatchNew* richMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(cand.fRichInd));
1613 if (richMatch == nullptr && richMatch->GetNofLinks() > 0) continue;
1614 cand.fRichMcTrackId = richMatch->GetMatchedLink().GetIndex();
1615
1616 // TRD
1617 CbmTrackMatchNew* trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(cand.fTrdInd));
1618 if (trdMatch == nullptr) continue;
1619 cand.fTrdMcTrackId = trdMatch->GetMatchedLink().GetIndex();
1620
1621 // ToF
1622 if (cand.fTofHitInd < 0) continue;
1623 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofHitInd));
1624 if (tofHit == nullptr) continue;
1625 CbmMatch* tofHitMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(cand.fTofHitInd));
1626 if (tofHitMatch == nullptr) continue;
1627 int tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
1628 if (tofPointIndex < 0) continue;
1629 FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofPointIndex));
1630 if (tofPoint == nullptr) continue;
1631 cand.fTofMcTrackId = tofPoint->GetTrackID();
1632 }
1633}
1634
1635void LmvmTask::AssignMcToTopologyCands(vector<LmvmCand>& topoCands)
1636{
1637 for (auto& cand : topoCands) {
1638 cand.ResetMcParams();
1639 if (cand.fStsInd < 0) continue;
1640 CbmMCTrack* mct = GetMcTrackSts(cand.fStsInd);
1641 if (mct == nullptr) continue;
1642 cand.fMcMotherId = mct->GetMotherId();
1643 cand.fMcPdg = mct->GetPdgCode();
1644 cand.fMcSrc = LmvmUtils::GetMcSrc(mct, fMCTracks);
1645 }
1646}
1647
1648void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step, const LmvmKinePar& parRec)
1649{
1650 ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM);
1651 double w = LmvmUtils::GetWeightPair(candP, candM);
1652
1653 fH.FillH1("hAnglePair", src, step, parRec.fAngle, w);
1654
1655 if (src == ELmvmSrc::Bg) {
1656 // gamma=0.5, pi0=1.5, pions=2.5, other=3.5
1657 double indM = candM.IsMcGamma() ? 0.5 : (candM.IsMcPi0() ? 1.5 : (std::abs(candM.fMcPdg) == 211) ? 2.5 : 3.5);
1658 double indP = candP.IsMcGamma() ? 0.5 : (candP.IsMcPi0() ? 1.5 : (std::abs(candP.fMcPdg) == 211) ? 2.5 : 3.5);
1659 fH.FillH2("hSrcBgPairsEpEm", step, indP, indM);
1660
1661 // Background
1662 ELmvmBgPairSrc bgSrc = LmvmUtils::GetBgPairSrc(candP, candM);
1663 if (bgSrc != ELmvmBgPairSrc::Undefined) {
1664 string name = fH.GetName("hMinvBgSource_" + fH.fBgPairSrcNames[static_cast<int>(bgSrc)], step);
1665 fH.FillH1(name, parRec.fMinv, w);
1666 fH.FillH2("hSrcBgPairs", static_cast<int>(step) + 0.5, static_cast<double>(bgSrc) + 0.5, w);
1667
1668 string hName = "hMinvBgSource2_" + fH.fAnaStepNames[static_cast<int>(step)] + "_";
1669
1670 if (std::abs(candP.fMcPdg) == 11) {
1671
1672 // cand1 is El and from Gamma
1673 if (candP.IsMcGamma()) {
1674 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma())
1675 fH.FillH1(hName + "gg", parRec.fMinv, w);
1676 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1677 fH.FillH1(hName + "gpi0", parRec.fMinv, w);
1678 else if (std::abs(candM.fMcPdg) == 211)
1679 fH.FillH1(hName + "gpi", parRec.fMinv, w);
1680 else
1681 fH.FillH1(hName + "go", parRec.fMinv, w);
1682 }
1683
1684 // cand1 is El and from Pi0
1685 else if (candP.IsMcPi0()) {
1686 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma())
1687 fH.FillH1(hName + "gpi0", parRec.fMinv, w);
1688 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1689 fH.FillH1(hName + "pi0pi0", parRec.fMinv, w);
1690 else if (std::abs(candM.fMcPdg) == 211)
1691 fH.FillH1(hName + "pipi0", parRec.fMinv, w);
1692 else
1693 fH.FillH1(hName + "pi0o", parRec.fMinv, w);
1694 }
1695
1696 // cand1 is El but not from Gamma or Pi0
1697 else {
1698 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma())
1699 fH.FillH1(hName + "go", parRec.fMinv, w);
1700 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1701 fH.FillH1(hName + "pi0o", parRec.fMinv, w);
1702 else if (std::abs(candM.fMcPdg) == 211)
1703 fH.FillH1(hName + "pio", parRec.fMinv, w);
1704 else
1705 fH.FillH1(hName + "oo", parRec.fMinv, w);
1706 }
1707 }
1708
1709 // cand1 is misid. charged pion
1710 else if (std::abs(candP.fMcPdg) == 211) {
1711 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma())
1712 fH.FillH1(hName + "gpi", parRec.fMinv, w);
1713 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1714 fH.FillH1(hName + "pipi0", parRec.fMinv, w);
1715 else if (std::abs(candM.fMcPdg) == 211)
1716 fH.FillH1(hName + "pipi", parRec.fMinv, w);
1717 else
1718 fH.FillH1(hName + "pio", parRec.fMinv, w);
1719 }
1720
1721 // cand1 is neither electron nor misid. charged pion
1722 else {
1723 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma())
1724 fH.FillH1(hName + "go", parRec.fMinv, w);
1725 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1726 fH.FillH1(hName + "pi0o", parRec.fMinv, w);
1727 else if (std::abs(candM.fMcPdg) == 211)
1728 fH.FillH1(hName + "pipi0", parRec.fMinv, w);
1729 else
1730 fH.FillH1(hName + "oo", parRec.fMinv, w);
1731 }
1732 }
1733 }
1734}
1735
1736void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg)
1737{
1738 // Histograms for MC step is filled in DoMcTrack()
1739 if (step == ELmvmAnaStep::Mc) return;
1740 FillMomHists(nullptr, &cand, cand.fMcSrc, step);
1741
1742 double w = cand.fWeight; // changed on 10.03.25
1743 double pdgBin = (std::abs(pdg) == 11 && cand.IsMcSignal()) ? 0.5
1744 : (std::abs(pdg) == 11 && !cand.IsMcGamma()) ? 1.5
1745 : (std::abs(pdg) == 11 && !cand.IsMcSignal()) ? 2.5
1746 : (std::abs(pdg) == 211) ? 3.5
1747 : (pdg == 2212) ? 4.5
1748 : (std::abs(pdg) == 321) ? 5.5
1749 : 6.5;
1750 fH.FillH2("hCandPdgVsMom", step, cand.fMomentum.Mag(), pdgBin, w);
1751 if (step == ELmvmAnaStep::ElId) {
1752 if (cand.fCharge > 0) fH.FillH2("hCandPdgVsMom+", cand.fMomentum.Mag(), pdgBin, w);
1753 if (cand.fCharge < 0) fH.FillH2("hCandPdgVsMom-", cand.fMomentum.Mag(), pdgBin, w);
1754 }
1755
1756 double stepBin = static_cast<double>(step) + 0.5;
1757 if (cand.IsMcSignal()) {
1758 fH.FillH1("hNofSignalTracks", stepBin, w);
1759 fH.FillH2("hCandElSrc", stepBin, 7.5, w);
1760 }
1761 else {
1762 if (LmvmUtils::IsMismatch(cand)) fH.FillH1("hNofMismatches_all", stepBin, w);
1763 if (cand.fStsMcTrackId != cand.fRichMcTrackId) fH.FillH1("hNofMismatches_rich", stepBin, w);
1764 if (cand.fStsMcTrackId != cand.fTrdMcTrackId) fH.FillH1("hNofMismatches_trd", stepBin, w);
1765 if (cand.fStsMcTrackId != cand.fTofMcTrackId) fH.FillH1("hNofMismatches_tof", stepBin, w);
1766 if (LmvmUtils::IsGhost(cand)) fH.FillH1("hNofGhosts", stepBin, w);
1767 fH.FillH1("hNofBgTracks", stepBin, w);
1768
1769 if (cand.IsMcGamma()) {
1770 CbmMCTrack* mctrack = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1771 if (mctrack != nullptr) {
1772 TVector3 v;
1773 mctrack->GetStartVertex(v);
1774 fH.FillH2("hVertexGammaXZ", step, v.Z(), v.X());
1775 fH.FillH2("hVertexGammaYZ", step, v.Z(), v.Y());
1776 fH.FillH2("hVertexGammaXY", step, v.X(), v.Y());
1777 fH.FillH2("hVertexGammaRZ", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
1778 }
1779 }
1780
1781 double srcBin = 0.0;
1782 if (cand.IsMcGamma())
1783 srcBin = 0.5;
1784 else if (cand.IsMcPi0())
1785 srcBin = 1.5;
1786 else if (std::abs(pdg) == 211)
1787 srcBin = 2.5;
1788 else if (pdg == 2212)
1789 srcBin = 3.5;
1790 else if (std::abs(pdg) == 321)
1791 srcBin = 4.5;
1792 else if ((std::abs(pdg) == 11) && !cand.IsMcGamma() && !cand.IsMcPi0() && !cand.IsMcSignal())
1793 srcBin = 5.5;
1794 else
1795 srcBin = 6.5;
1796 fH.FillH2("hBgSrcTracks", stepBin, srcBin, w);
1797 if (std::abs(cand.fMcPdg) == 11) fH.FillH2("hCandElSrc", stepBin, srcBin, w);
1798 }
1799}
1800
1801void LmvmTask::FillPairHists(const LmvmCand& cand1, const LmvmCand& cand2, const LmvmKinePar& pMc,
1802 const LmvmKinePar& pRec, ELmvmAnaStep step)
1803{
1804 // no need to fill histograms for MC step
1805 if (step == ELmvmAnaStep::Mc) return;
1806 ELmvmSrc src = LmvmUtils::GetMcPairSrc(cand1, cand2);
1807 int pdg1 = cand1.fMcPdg;
1808 int pdg2 = cand2.fMcPdg;
1809 double w = LmvmUtils::GetWeightPair(cand1, cand2);
1810 string comb = "";
1811 if (cand1.fCharge * cand2.fCharge < 0)
1812 comb = "PM";
1813 else if (cand1.fCharge > 0 && cand2.fCharge > 0)
1814 comb = "PP";
1815 else if (cand1.fCharge < 0 && cand2.fCharge < 0)
1816 comb = "MM";
1817
1818 if (step == ELmvmAnaStep::ElId) {
1819 double bin1 = -1.;
1820 double bin2 = -1.;
1821
1822 if (cand1.IsMcSignal())
1823 bin1 = 0.5;
1824 else if (cand1.IsMcPi0())
1825 bin1 = 1.5;
1826 else if (cand1.IsMcGamma())
1827 bin1 = 2.5;
1828 else if (std::abs(pdg1) == 11 && !cand1.IsMcSignal() && !cand1.IsMcGamma() && !cand1.IsMcPi0())
1829 bin1 = 3.5;
1830 else if (std::abs(pdg1) == 211)
1831 bin1 = 4.5;
1832 else if (pdg1 == 2212)
1833 bin1 = 5.5;
1834 else if (std::abs(pdg1) == 321)
1835 bin1 = 6.5;
1836 else
1837 bin1 = 7.5;
1838
1839 if (cand2.IsMcSignal())
1840 bin2 = 0.5;
1841 else if (cand2.IsMcPi0())
1842 bin2 = 1.5;
1843 else if (cand2.IsMcGamma())
1844 bin2 = 2.5;
1845 else if (std::abs(pdg2) == 11 && !cand2.IsMcSignal() && !cand2.IsMcGamma() && !cand2.IsMcPi0())
1846 bin2 = 3.5;
1847 else if (std::abs(pdg2) == 211)
1848 bin2 = 4.5;
1849 else if (pdg2 == 2212)
1850 bin2 = 5.5;
1851 else if (std::abs(pdg2) == 321)
1852 bin2 = 6.5;
1853 else
1854 bin2 = 7.5;
1855
1856 int mBinMc = (int) (pMc.fMinv * 10);
1857 int mBinRec = (int) (pRec.fMinv * 10);
1858 //if (mBinMc >= 25) mBinMc = 24; // consider overflow bins
1859 //if (mBinRec >= 25) mBinRec = 24;
1860 string hNameMc = "hBgPairPdg" + comb + "_pMc_" + Cbm::NumberToString(mBinMc, 0);
1861 string hNameRec = "hBgPairPdg" + comb + "_pRec_" + Cbm::NumberToString(mBinRec, 0);
1862 if (cand1.fCharge < 0 && cand2.fCharge > 0) {
1863 if (mBinMc <= 24) fH.FillH2(hNameMc, bin2, bin1, w);
1864 if (mBinRec <= 24) fH.FillH2(hNameRec, bin2, bin1, w);
1865 }
1866 else {
1867 if (mBinMc <= 24) fH.FillH2(hNameMc, bin1, bin2, w);
1868 if (mBinRec <= 24) fH.FillH2(hNameRec, bin1, bin2, w);
1869 }
1870 }
1871
1872 // From here on only unlike-sign pairs!
1873 if (comb != "PM") return;
1874 bool isMismatch = (LmvmUtils::IsMismatch(cand1) || LmvmUtils::IsMismatch(cand2));
1875 fH.FillH1("hMinv", src, step, pRec.fMinv, w);
1876 if (std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11) fH.FillH1("hMinv_trueEl", src, step, pRec.fMinv, w);
1877 if (!cand1.IsMcSignal() && !cand2.IsMcSignal()) fH.FillH1("hMinv_urqmdAll", src, step, pRec.fMinv, w);
1878 if (!cand1.IsMcSignal() && !cand2.IsMcSignal() && std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11)
1879 fH.FillH1("hMinv_urqmdEl", src, step, pRec.fMinv, w);
1880 if (step == ELmvmAnaStep::ElId && src == ELmvmSrc::Bg) {
1881 if (std::abs(pdg1) == 11 && std::abs(pdg2) == 11) fH.FillH1("hMinv_tt", pRec.fMinv, w);
1882 if (std::abs(pdg1) == 11 && std::abs(pdg2) != 11) fH.FillH1("hMinv_tm", pRec.fMinv, w);
1883 if (std::abs(pdg1) != 11 && std::abs(pdg2) == 11) fH.FillH1("hMinv_mt", pRec.fMinv, w);
1884 if (std::abs(pdg1) != 11 && std::abs(pdg2) != 11) fH.FillH1("hMinv_mm", pRec.fMinv, w);
1885 }
1886
1887 if (step == ELmvmAnaStep::ElId) {
1888 double rat = pRec.fMinv / pMc.fMinv;
1889 fH.FillH2("hRecoPrec_Minv", pMc.fMinv, rat, w);
1890 }
1891
1892 fH.FillH2("hMinvPt", src, step, pRec.fMinv, pMc.fPt, w);
1893
1894 PairSource(cand1, cand2, step, pRec);
1895 if (src == ELmvmSrc::Signal) {
1896 fH.FillH2("hPtYPairSignal", step, pMc.fRapidity, pMc.fPt, w);
1897 fH.FillH1("hMomPairSignal", step, pMc.fMomentumMag, w);
1898 }
1899 if (src == ELmvmSrc::Bg) {
1900 if (isMismatch) {
1901 fH.FillH1("hMinvBgMatch_mismatch", step, pRec.fMinv, w);
1902 }
1903 else {
1904 fH.FillH1("hMinvBgMatch_trueMatch", step, pRec.fMinv, w);
1905 if (std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11)
1906 fH.FillH1("hMinvBgMatch_trueMatchEl", step, pRec.fMinv, w);
1907 if (std::abs(cand1.fMcPdg) != 11 || std::abs(cand2.fMcPdg) != 11)
1908 fH.FillH1("hMinvBgMatch_trueMatchNotEl", step, pRec.fMinv, w);
1909 }
1910 }
1911}
1912
1913string LmvmTask::GetPidString(const CbmMCTrack* mct, const LmvmCand* cand)
1914{
1915 if ((mct != nullptr && cand != nullptr) || (mct == nullptr && cand == nullptr)) {
1916 LOG(error) << "LmvmTask::GetPidString: Both mct and cand are [not nullptr] or [nullptr].";
1917 return 0;
1918 }
1919
1920 int pdg = (mct != nullptr) ? mct->GetPdgCode() : cand->fMcPdg;
1921 bool isMcSignal = (mct != nullptr) ? (mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11)
1922 : cand->IsMcSignal();
1923
1924 if (isMcSignal && pdg == -11)
1925 return fH.fCandNames[0];
1926 else if (isMcSignal && pdg == 11)
1927 return fH.fCandNames[1];
1928 else if (!isMcSignal && pdg == -11)
1929 return fH.fCandNames[2];
1930 else if (!isMcSignal && pdg == 11)
1931 return fH.fCandNames[3];
1932 else if (pdg == 211)
1933 return fH.fCandNames[4];
1934 else if (pdg == -211)
1935 return fH.fCandNames[5];
1936 else if (pdg == 2212)
1937 return fH.fCandNames[6];
1938 else if (pdg == 321)
1939 return fH.fCandNames[7];
1940 else if (pdg == -321)
1941 return fH.fCandNames[8];
1942 return fH.fCandNames[9];
1943}
1944
1945string LmvmTask::GetPidString(double vertexZ, int pdg)
1946{
1947 string pidString = fH.fGTrackNames[14];
1948 bool isPrim = IsPrimary(vertexZ);
1949
1950 if (isPrim) {
1951 if (pdg == -11)
1952 pidString = fH.fGTrackNames[1];
1953 else if (pdg == 11)
1954 pidString = fH.fGTrackNames[3];
1955 else if (pdg == 211)
1956 pidString = fH.fGTrackNames[5];
1957 else if (pdg == -211)
1958 pidString = fH.fGTrackNames[7];
1959 else if (pdg == 2212)
1960 pidString = fH.fGTrackNames[9];
1961 else if (pdg == 321)
1962 pidString = fH.fGTrackNames[11];
1963 else if (pdg == -321)
1964 pidString = fH.fGTrackNames[13];
1965 }
1966 else {
1967 if (pdg == -11)
1968 pidString = fH.fGTrackNames[0];
1969 else if (pdg == 11)
1970 pidString = fH.fGTrackNames[2];
1971 else if (pdg == 211)
1972 pidString = fH.fGTrackNames[4];
1973 else if (pdg == -211)
1974 pidString = fH.fGTrackNames[6];
1975 else if (pdg == 2212)
1976 pidString = fH.fGTrackNames[8];
1977 else if (pdg == 321)
1978 pidString = fH.fGTrackNames[10];
1979 else if (pdg == -321)
1980 pidString = fH.fGTrackNames[12];
1981 }
1982
1983 return pidString;
1984}
1985
1986CbmMCTrack* LmvmTask::GetMcTrackSts(int stsIndex)
1987{
1988 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsIndex));
1989 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) return nullptr;
1990 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1991 if (stsMcTrackId < 0) return nullptr;
1992 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(stsMcTrackId));
1993 return mct;
1994}
1995
1996void LmvmTask::AnalyseCandidates()
1997{
1998 CheckGammaConvAndPi0();
1999 CheckTopologyCut(ELmvmTopologyCut::ST, "hStCut");
2000 CheckTopologyCut(ELmvmTopologyCut::TT, "hTtCut");
2001 CheckTopologyCut(ELmvmTopologyCut::RT, "hRtCut");
2002 if (fUseMvd) {
2003 CheckClosestMvdHit(1, "hMvdCut_1", "hMvdCutQa_1");
2004 CheckClosestMvdHit(2, "hMvdCut_2", "hMvdCutQa_2");
2005 }
2006 // single candidates
2007 for (const auto& cand : fCands) {
2008 double w = cand.fWeight;
2009 CbmMCTrack* mcTrack = GetMcTrackSts(cand.fStsInd);
2010 if (mcTrack == nullptr) continue;
2011
2012 int pdg = mcTrack->GetPdgCode();
2013 double pRec = cand.fMomentum.Mag();
2014 double pMc = mcTrack->GetP();
2015 string pidString = GetPidString(nullptr, &cand);
2016
2017 bool isMismatch = LmvmUtils::IsMismatch(cand);
2018 string matchString = (isMismatch) ? "mismatch" : "truematch";
2019 fH.FillH2("hTofDist_" + matchString + "_" + pidString, pRec, cand.fTofDist, w);
2020
2021 for (auto step : fH.fAnaSteps) {
2022 fH.FillH2("hPtY_" + pidString, step, cand.fRapidity, cand.fMomentum.Perp(), w);
2023 if (step == ELmvmAnaStep::Mc) continue; // Mc has to be filled in DoMcTracks()
2024 if (cand.IsCutTill(step)) {
2025 TrackSource(cand, step, pdg);
2026 fH.FillH1("hCandPdg", step, cand.fMcPdg, w);
2027 fH.FillH1("hMom_" + pidString, step, pRec, w); // changed pMc to pRec on 2.7.24
2028 fH.FillH2("hMomPt_" + pidString, step, pRec, cand.fMomentum.Perp(), w);
2029 fH.FillH2("hELossSts_" + pidString, step, pRec, cand.fELossSts, w);
2030 fH.FillH2("hRichAnn_" + pidString, step, pRec, cand.fRichAnn, w);
2031 fH.FillH2("hTrdElLike_" + pidString, step, pRec, cand.fTrdLikeEl, w);
2032 fH.FillH2("hChi2Dets_sts_" + pidString, step, pRec, cand.fChi2Sts, w);
2033 fH.FillH2("hChi2Dets_rich_" + pidString, step, pRec, cand.fChi2Rich, w);
2034 fH.FillH2("hChi2Dets_trd_" + pidString, step, pRec, cand.fChi2Trd, w);
2035 fH.FillH2("hTofM2_" + pidString, step, pRec, cand.fMass2, w);
2036 fH.FillH2("hTofDist_" + pidString, step, pRec, cand.fTofDist, w);
2037 if (step == ELmvmAnaStep::ElId) fH.FillH2("hTofChi2_elid_" + pidString, pRec, cand.fChi2Tof, w);
2038
2039 if (mcTrack != nullptr) {
2040 RatioMomentum(mcTrack, cand, step, pdg);
2041 CheckTofId(mcTrack, cand, step, pdg);
2042 }
2043
2044 // Pion and proton suppression (filled with pMc, not pRec, because compared to MC track later)
2045 if (step == ELmvmAnaStep::Reco) {
2046 if (std::abs(pdg) == 211) fH.FillH1("hSuppression_pion_reco", pMc, w);
2047 if (pdg == 2212) fH.FillH1("hSuppression_proton_reco", pMc, w);
2048 }
2049 if (step == ELmvmAnaStep::Acc) {
2050 if (std::abs(pdg) == 211) fH.FillH1("hSuppression_pion_acc", pMc, w);
2051 if (pdg == 2212) fH.FillH1("hSuppression_proton_acc", pMc, w);
2052 }
2053 if (step == ELmvmAnaStep::ElId) {
2054 if (std::abs(pdg) == 211) fH.FillH1("hSuppression_pion_elid", pMc, w);
2055 if (pdg == 2212) fH.FillH1("hSuppression_proton_elid", pMc, w);
2056 }
2057 if (step == ELmvmAnaStep::TtCut) {
2058 if (std::abs(pdg) == 211) fH.FillH1("hSuppression_pion_ttcut", pMc, w);
2059 if (pdg == 2212) fH.FillH1("hSuppression_proton_ttcut", pMc, w);
2060 }
2061 if (step == ELmvmAnaStep::PtCut) {
2062 if (std::abs(pdg) == 211) fH.FillH1("hSuppression_pion_ptcut", pMc, w);
2063 if (pdg == 2212) fH.FillH1("hSuppression_proton_ptcut", pMc, w);
2064 }
2065
2066 double richDist = CbmRichUtil::GetRingTrackDistance(cand.fGTrackInd);
2067 fH.FillH1("hRichRingTrackDist_" + pidString, step, richDist);
2068 fH.FillH1("hRichRingTrackDistVsMom_" + pidString, step, pRec, richDist);
2069 }
2070 } // steps
2071
2072 FillSourceHistos(cand);
2073
2074 if (cand.IsCutTill(ELmvmAnaStep::ElId)) {
2075 fH.FillH2("hChi2Comb_StsRich_" + pidString, cand.fChi2Sts, cand.fChi2Rich, w);
2076 fH.FillH2("hChi2Comb_StsTrd_" + pidString, cand.fChi2Sts, cand.fChi2Trd, w);
2077 fH.FillH2("hChi2Comb_RichTrd_" + pidString, cand.fChi2Rich, cand.fChi2Trd, w);
2078
2079 // Double Hits in TOF
2080 int nInd = 0;
2081 for (int iT = 0; iT < fTofTracks->GetEntriesFast(); iT++) {
2082 CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(iT));
2083 if (tofTrack->GetTofHitIndex() == cand.fTofHitInd) nInd++;
2084 }
2085 if (nInd == 1)
2086 fH.FillH2("hTofM2_woDoubleIndex_1_" + pidString, pRec, cand.fMass2, w);
2087 else if (nInd == 2)
2088 fH.FillH2("hTofM2_woDoubleIndex_2_" + pidString, pRec, cand.fMass2, w);
2089 else if (nInd > 2)
2090 fH.FillH2("hTofM2_woDoubleIndex_3+_" + pidString, pRec, cand.fMass2, w);
2091 }
2092
2093 //TOF
2094 fH.FillH2("hTofChi2_all_" + pidString, pRec, cand.fChi2Tof, w);
2095 CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(cand.fTofHitInd);
2096 double hitX = tofHit->GetX();
2097 double hitY = tofHit->GetY();
2098
2099 // TOF-M2 vs. position in TOF
2100 double width = 20.;
2101 for (int pos : {50, 100, 150, 200, 300, 400, 500}) {
2102 if (cand.fIsChi2Prim) {
2103 string hName = "hTofM2_chi2prim_" + Cbm::NumberToString(pos, 0) + "_" + pidString;
2104 if (IsInBox(pos, -pos, hitX, hitY, width)) fH.FillH2(hName, pMc, cand.fMass2, w);
2105 }
2106 if (cand.fIsElectron) {
2107 string hName = "hTofM2_elid_" + Cbm::NumberToString(pos, 0) + "_" + pidString;
2108 if (IsInBox(pos, -pos, hitX, hitY, width)) fH.FillH2(hName, pMc, cand.fMass2, w);
2109 }
2110 }
2111
2112 // TOF-m2 vs. true-/mismatch
2113 if (cand.IsCutTill(ELmvmAnaStep::Chi2Prim)) {
2114 if (cand.fStsMcTrackId == cand.fTofMcTrackId)
2115 fH.FillH2("hTofM2_truematch_" + pidString, pMc, cand.fMass2, w);
2116 else
2117 fH.FillH2("hTofM2_mismatch_" + pidString, pMc, cand.fMass2, w);
2118 }
2119
2120 // E-Loss vs. M2
2121 if (tofHit != nullptr) {
2122 int x1 = -300;
2123 int y1 = -300;
2124 int x2 = 80;
2125 int y2 = 80;
2126 int wi = 20;
2127 if (cand.fIsChi2Prim) {
2128 if (hitX > (x1 - wi) && hitX < (x1 + wi) && hitY > (y1 - wi) && hitY < (y1 + wi) && pMc > 0.5 && pMc <= 1.0)
2129 fH.FillH2("hElossTime_chi2prim_Seg1Mom1_" + pidString, cand.fTime, cand.fELossSts, w);
2130 if (hitX > (x1 - wi) && hitX < (x1 + wi) && hitY > (y1 - wi) && hitY < (y1 + wi) && pMc > 1.0 && pMc <= 2.0)
2131 fH.FillH2("hElossTime_chi2prim_Seg1Mom2_" + pidString, cand.fTime, cand.fELossSts, w);
2132 if (hitX > (x2 - wi) && hitX < (x2 + wi) && hitY > (y2 - wi) && hitY < (y2 + wi) && pMc > 0.5 && pMc <= 1.0)
2133 fH.FillH2("hElossTime_chi2prim_Seg2Mom1_" + pidString, cand.fTime, cand.fELossSts, w);
2134 if (hitX > (x2 - wi) && hitX < (x2 + wi) && hitY > (y2 - wi) && hitY < (y2 + wi) && pMc > 1.0 && pMc <= 2.0)
2135 fH.FillH2("hElossTime_chi2prim_Seg2Mom2_" + pidString, cand.fTime, cand.fELossSts, w);
2136 }
2137 if (cand.fIsElectron) {
2138 if (hitX > (x1 - wi) && hitX < (x1 + wi) && hitY > (y1 - wi) && hitY < (y1 + wi) && pMc > 0.5 && pMc <= 1.0)
2139 fH.FillH2("hElossTime_elid_Seg1Mom1_" + pidString, cand.fTime, cand.fELossSts, w);
2140 if (hitX > (x1 - wi) && hitX < (x1 + wi) && hitY > (y1 - wi) && hitY < (y1 + wi) && pMc > 1.0 && pMc <= 2.0)
2141 fH.FillH2("hElossTime_elid_Seg1Mom2_" + pidString, cand.fTime, cand.fELossSts, w);
2142 if (hitX > (x2 - wi) && hitX < (x2 + wi) && hitY > (y2 - wi) && hitY < (y2 + wi) && pMc > 0.5 && pMc <= 1.0)
2143 fH.FillH2("hElossTime_elid_Seg2Mom1_" + pidString, cand.fTime, cand.fELossSts, w);
2144 if (hitX > (x2 - wi) && hitX < (x2 + wi) && hitY > (y2 - wi) && hitY < (y2 + wi) && pMc > 1.0 && pMc <= 2.0)
2145 fH.FillH2("hElossTime_elid_Seg2Mom2_" + pidString, cand.fTime, cand.fELossSts, w);
2146 }
2147 }
2148
2149 // beta-momentum spectrum
2150 string ptcl = GetPidString(nullptr, &cand);
2151 double m = cand.fMass;
2152 double r = pRec / m;
2153 double beta = r / sqrt(1. + r * r);
2154 double q = (cand.fCharge > 0) ? 1 : (cand.fCharge < 0) ? -1. : 0.;
2155 fH.FillH2("hBetaMom_" + ptcl, q * pRec, beta);
2156 } // single candidates
2157
2158 // candidate pairs
2159 for (size_t iC1 = 0; iC1 < fCands.size() - 1; iC1++) {
2160 const LmvmCand& cand1 = fCands[iC1];
2161 CbmMCTrack* mctrack1 = GetMcTrackSts(cand1.fStsInd);
2162 if (mctrack1 == nullptr) continue;
2163 for (size_t iC2 = iC1 + 1; iC2 < fCands.size(); iC2++) {
2164 const LmvmCand& cand2 = fCands[iC2];
2165 CbmMCTrack* mctrack2 = GetMcTrackSts(cand2.fStsInd);
2166 if (mctrack2 == nullptr) continue;
2167 LmvmKinePar pMC = LmvmKinePar::Create(mctrack1, mctrack2);
2168 LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2);
2169 for (auto step : fH.fAnaSteps) {
2170 if (cand1.IsCutTill(step) && cand2.IsCutTill(step)) FillPairHists(cand1, cand2, pMC, pRec, step);
2171 }
2172 }
2173 }
2174}
2175
2176bool LmvmTask::IsInBox(double boxX, double boxY, double xVal, double yVal, double width)
2177{
2178 return (xVal >= (boxX - width) && xVal <= (boxX + width) && yVal >= (boxY - width) && yVal <= (boxY + width));
2179}
2180
2181void LmvmTask::AnalyseProperties()
2182{
2183 for (const auto& cand : fCands) {
2184 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue;
2185 double w = cand.fWeight;
2186 string hName = "hCandProperties";
2187 string ptclString = GetPidString(nullptr, &cand);
2188 string pRange = "";
2189 if (cand.fMomentum.Mag() >= 0 && cand.fMomentum.Mag() < 1.)
2190 pRange = "1";
2191 else if (cand.fMomentum.Mag() >= 1. && cand.fMomentum.Mag() < 2.)
2192 pRange = "2";
2193 else if (cand.fMomentum.Mag() >= 2. && cand.fMomentum.Mag() < 3.)
2194 pRange = "3";
2195 else if (cand.fMomentum.Mag() >= 3. && cand.fMomentum.Mag() < 4.)
2196 pRange = "4";
2197 else if (cand.fMomentum.Mag() >= 4. && cand.fMomentum.Mag() < 5.)
2198 pRange = "5";
2199 else if (cand.fMomentum.Mag() >= 5. && cand.fMomentum.Mag() < 6.)
2200 pRange = "6";
2201 else if (cand.fMomentum.Mag() >= 6. && cand.fMomentum.Mag() < 7.)
2202 pRange = "7";
2203 else if (cand.fMomentum.Mag() >= 7. && cand.fMomentum.Mag() < 8.)
2204 pRange = "8";
2205 else if (cand.fMomentum.Mag() >= 8. && cand.fMomentum.Mag() < 9.)
2206 pRange = "9";
2207 else if (cand.fMomentum.Mag() >= 9. && cand.fMomentum.Mag() <= 10.)
2208 pRange = "10";
2209 else if (cand.fMomentum.Mag() > 10.)
2210 pRange = "11";
2211 else
2212 continue;
2213 string hFullName = hName + "_" + pRange + "_" + ptclString;
2214
2215 double px = cand.fMomentum.X();
2216 double py = cand.fMomentum.Y();
2217 double chi2Prim = cand.fChi2Prim;
2218 double chi2Sts = cand.fChi2Sts;
2219 double chi2Rich = cand.fChi2Rich;
2220 double fChi2Trd = cand.fChi2Trd / 2.;
2221 double stsELoss = cand.fELossSts / 1e4;
2222 double richAnn = cand.fRichAnn + 2.;
2223 double trdLikeEl = cand.fTrdLikeEl;
2224 double trdLikePi = cand.fTrdLikePi;
2225 double tofM2 = cand.fMass2 + 1.;
2226 double toftime = cand.fTime;
2227 double nHitsSts = cand.fNofHitsSts;
2228 double nHitsRich = cand.fNofHitsRich / 5.;
2229 double nHitsTrd = cand.fNofHitsTrd;
2230 double chi2Tof = cand.fChi2Tof;
2231
2232 fH.FillH2(hFullName, 0.5, px, w);
2233 fH.FillH2(hFullName, 1.5, py, w);
2234 fH.FillH2(hFullName, 2.5, chi2Prim, w);
2235 fH.FillH2(hFullName, 3.5, chi2Sts, w);
2236 fH.FillH2(hFullName, 4.5, chi2Rich, w);
2237 fH.FillH2(hFullName, 5.5, fChi2Trd, w);
2238 fH.FillH2(hFullName, 6.5, stsELoss, w);
2239 fH.FillH2(hFullName, 7.5, richAnn, w);
2240 fH.FillH2(hFullName, 8.5, trdLikeEl, w);
2241 fH.FillH2(hFullName, 9.5, trdLikePi, w);
2242 fH.FillH2(hFullName, 10.5, tofM2, w);
2243 fH.FillH2(hFullName, 11.5, toftime, w);
2244 fH.FillH2(hFullName, 12.5, chi2Tof, w);
2245 fH.FillH2(hFullName, 13.5, nHitsSts, w);
2246 fH.FillH2(hFullName, 14.5, nHitsRich, w);
2247 fH.FillH2(hFullName, 15.5, nHitsTrd, w);
2248 }
2249}
2250
2251void LmvmTask::CheckTofId(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step, int pdg)
2252{
2253 TVector3 v;
2254 mct->GetStartVertex(v);
2255 // check vertex of misidentified particles in ToF after electron ID
2256 if (std::abs(pdg) != 11 && cand.fIsTofElectron) {
2257 fH.FillH2("hVertexXZ_misidTof", step, v.Z(), v.X());
2258 fH.FillH2("hVertexYZ_misidTof", step, v.Z(), v.Y());
2259 fH.FillH2("hVertexXY_misidTof", step, v.X(), v.Y());
2260 fH.FillH2("hVertexRZ_misidTof", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
2261 }
2262}
2263
2264void LmvmTask::RatioMomentum(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step, int pdg)
2265{
2266 double w = cand.fWeight; // changed on 10.03.25
2267 double rat = cand.fMomentum.Mag() / mct->GetP();
2268 string pidString = GetPidString(nullptr, &cand);
2269 fH.FillH2("hRecoPrec_Mom_" + pidString, step, mct->GetP(), rat, w);
2270 if (step == ELmvmAnaStep::ElId && cand.IsCutTill(step) && std::abs(pdg) == 11)
2271 fH.FillH2("hRecoPrec_Mom_ID_true", mct->GetP(), rat, w);
2272 if (step == ELmvmAnaStep::ElId && cand.IsCutTill(step) && std::abs(pdg) != 11)
2273 fH.FillH2("hRecoPrec_Mom_ID_misid", mct->GetP(), rat, w);
2274}
2275
2276void LmvmTask::CheckGammaConvAndPi0()
2277{
2278 for (auto& candP : fCands) {
2279 if (candP.fCharge < 0) continue;
2280 for (auto& candM : fCands) {
2281 if (candM.fCharge > 0) continue;
2282 if (candP.IsCutTill(ELmvmAnaStep::ElId) && candM.IsCutTill(ELmvmAnaStep::ElId)) {
2283 LmvmKinePar pRec = LmvmKinePar::Create(&candP, &candM);
2284 if (!fCuts.IsGammaCutOk(pRec.fMinv)) {
2285 candM.fIsGammaCut = false;
2286 candP.fIsGammaCut = false;
2287 }
2288 }
2289 }
2290 }
2291}
2292
2293void LmvmTask::CheckTopologyCut(ELmvmTopologyCut cut, const string& name)
2294{
2295 string hcut = name + "_all";
2296 string hcutPion = name + "_pion";
2297 string hcutTruePair = name + "_truePair";
2298
2299 vector<LmvmDataAngMomInd> dataV;
2300
2301 vector<LmvmCand>& tpCands = fSTCands;
2302 if (cut == ELmvmTopologyCut::ST) {
2303 tpCands = fSTCands;
2304 }
2305 else if (cut == ELmvmTopologyCut::RT) {
2306 tpCands = fRTCands;
2307 }
2308 else if (cut == ELmvmTopologyCut::TT) {
2309 tpCands = fTTCands;
2310 }
2311 else {
2312 LOG(error) << "LmvmTask::CheckTopologyCut(): cut is not defined.";
2313 }
2314
2315 for (auto& cand : fCands) {
2316 if (cand.IsCutTill(ELmvmAnaStep::ElId)) {
2317 dataV.clear();
2318 for (size_t iM = 0; iM < tpCands.size(); iM++) {
2319 // different charges, charge iM != charge iP
2320 if (tpCands[iM].fCharge != cand.fCharge) {
2321 LmvmKinePar pRec = LmvmKinePar::Create(&cand, &tpCands[iM]);
2322 dataV.emplace_back(pRec.fAngle, tpCands[iM].fMomentum.Mag(), iM);
2323 }
2324 }
2325 //find min opening angle
2326 double minAng = 360.;
2327 int minInd = -1;
2328 for (size_t i = 0; i < dataV.size(); i++) {
2329 if (minAng > dataV[i].fAngle) {
2330 minAng = dataV[i].fAngle;
2331 minInd = i;
2332 }
2333 }
2334 if (minInd == -1) {
2335 cand.SetIsTopologyCutElectron(cut, true);
2336 continue;
2337 }
2338 bool isCut = fCuts.IsTopologyCutOk(cut, cand.fMomentum.Mag(), dataV[minInd].fMom, minAng);
2339 cand.SetIsTopologyCutElectron(cut, isCut);
2340
2341 // histogramms
2342 double sqrt_mom = TMath::Sqrt(cand.fMomentum.Mag() * dataV[minInd].fMom);
2343 int cutCandInd = dataV[minInd].fInd;
2344 int stsInd = tpCands[cutCandInd].fStsInd;
2345 if (stsInd < 0) continue;
2346 int pdgAbs = std::abs(tpCands[cutCandInd].fMcPdg);
2347 int motherId = tpCands[cutCandInd].fMcMotherId;
2348 double w = cand.fWeight; // changed on 10.03.25
2349
2350 fH.FillH2(hcut, cand.fMcSrc, sqrt_mom, minAng, w);
2351 if (pdgAbs == 211) fH.FillH2(hcutPion, cand.fMcSrc, sqrt_mom, minAng, w);
2352 if (cand.IsMcSignal()) {
2353 if (motherId == cand.fMcMotherId) fH.FillH2(hcutTruePair, cand.fMcSrc, sqrt_mom, minAng, w);
2354 }
2355 else {
2356 if (motherId != -1 && motherId == cand.fMcMotherId) fH.FillH2(hcutTruePair, cand.fMcSrc, sqrt_mom, minAng, w);
2357 }
2358 }
2359 }
2360}
2361
2362void LmvmTask::CalculateNofTopologyPairs(const string& name, ELmvmSrc src)
2363{
2364 size_t nCand = fCands.size();
2365 for (size_t iP = 0; iP < nCand; iP++) {
2366 const LmvmCand& cand = fCands[iP];
2367 if (cand.fMcMotherId == -1) continue;
2368 if (src != cand.fMcSrc) continue;
2369 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue;
2370
2371 bool isAdded = false;
2372
2373 // 3 topology cuts: ST, RT, TT
2374 for (int i = 0; i < 3; i++) {
2375 if (isAdded) continue;
2376 vector<LmvmCand>& cands = fSTCands;
2377 double binNum = 4.5;
2378 if (i == 1) {
2379 cands = fRTCands;
2380 binNum = 5.5;
2381 }
2382 else if (i == 2) {
2383 cands = fTTCands;
2384 binNum = 6.5;
2385 }
2386 for (const auto& candT : cands) {
2387 if (candT.fMcMotherId == cand.fMcMotherId) {
2388 fH.FillH1(name, binNum);
2389 isAdded = true;
2390 break;
2391 }
2392 }
2393 }
2394 if (isAdded) continue;
2395
2396 for (size_t iM = 0; iM < fCands.size(); iM++) {
2397 if (iM != iP && fCands[iM].fMcMotherId == cand.fMcMotherId && fCands[iM].IsCutTill(ELmvmAnaStep::ElId)) {
2398 fH.FillH1(name, 7.5);
2399 isAdded = true;
2400 break;
2401 }
2402 }
2403 if (isAdded) continue;
2404
2405 int nofStsPoints = 0;
2406 int nofMcTracks = fMCTracks->GetEntriesFast();
2407 for (int iMc = 0; iMc < nofMcTracks; iMc++) {
2408 const CbmMCTrack* mcTrack = static_cast<const CbmMCTrack*>(fMCTracks->At(iMc));
2409 if (mcTrack == nullptr || mcTrack->GetMotherId() != cand.fMcMotherId || iMc == cand.fStsMcTrackId) continue;
2410
2411 int eventId = FairRun::Instance()->GetEventHeader()->GetMCEntryNumber();
2412 if (!CbmLitMCTrackCreator::Instance()->TrackExists(eventId, iMc)) continue;
2413 const CbmLitMCTrack& litMCTrack = CbmLitMCTrackCreator::Instance()->GetTrack(eventId, iMc);
2414 nofStsPoints = litMCTrack.GetNofPointsInDifferentStations(ECbmModuleId::kSts);
2415 break;
2416 }
2417 if (nofStsPoints == 0) fH.FillH1(name, 0.5);
2418 if (nofStsPoints >= 1 && nofStsPoints <= 3) fH.FillH1(name, 1.5);
2419 if (nofStsPoints >= 4 && nofStsPoints <= 5) fH.FillH1(name, 2.5);
2420 if (nofStsPoints >= 6) fH.FillH1(name, 3.5);
2421 }
2422}
2423
2424void LmvmTask::FillSourceHistos(const LmvmCand& cand)
2425{
2426 double w = cand.fWeight; // changed on 10.03.25
2427 fH.FillH1("hChi2PrimVertex", cand.fMcSrc, cand.fChi2Prim, w);
2428
2429 if (!cand.fIsChi2Prim) return;
2430 fH.FillH2("hRichAnnSrc", cand.fMcSrc, cand.fMomentum.Mag(), cand.fRichAnn, w);
2431 fH.FillH2("hTofM2Src", cand.fMcSrc, cand.fMomentum.Mag(), cand.fMass2, w);
2432
2433 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) return;
2434 //fH.FillSourceH1("hPt", cand.fMcSrc, cand.fMomentum.Perp(), w);
2435 //fH.FillSourceH1("hMom", cand.fMcSrc, cand.fMomentum.Mag(), w);
2436 fH.FillH1("hChi2Sts", cand.fMcSrc, cand.fChi2Sts, w);
2437
2438 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
2439 if (stsTrack == nullptr) return;
2440 fH.FillH1("hNofStsHits", cand.fMcSrc, stsTrack->GetNofStsHits(), w);
2441
2442 if (fUseMvd) {
2443 double mvd1x = 0., mvd1y = 0., mvd2x = 0., mvd2y = 0.;
2444 for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
2445 CbmMvdHit* mvdHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
2446 if (mvdHit == nullptr) return;
2447 if (mvdHit->GetStationNr() == 1) {
2448 mvd1x = mvdHit->GetX();
2449 mvd1y = mvdHit->GetY();
2450 }
2451 else if (mvdHit->GetStationNr() == 2) {
2452 mvd2x = mvdHit->GetX();
2453 mvd2y = mvdHit->GetY();
2454 }
2455 }
2456 double mvd1r = sqrt(mvd1x * mvd1x + mvd1y * mvd1y);
2457 double mvd2r = sqrt(mvd2x * mvd2x + mvd2y * mvd2y);
2458 fH.FillH1("hNofMvdHits", cand.fMcSrc, stsTrack->GetNofMvdHits(), w);
2459 fH.FillH2("hMvdXY_1", cand.fMcSrc, mvd1x, mvd1y, w);
2460 fH.FillH1("hMvdR_1", cand.fMcSrc, mvd1r, w);
2461 fH.FillH2("hMvdXY_2", cand.fMcSrc, mvd2x, mvd2y, w);
2462 fH.FillH1("hMvdR_2", cand.fMcSrc, mvd2r, w);
2463 }
2464}
2465
2466void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const string& histQa)
2467{
2468 vector<LmvmDataXYInd> mvdV;
2469 vector<LmvmDataXYInd> candV;
2470
2471 for (int iHit = 0; iHit < fMvdHits->GetEntriesFast(); iHit++) {
2472 CbmMvdHit* mvdHit = static_cast<CbmMvdHit*>(fMvdHits->At(iHit));
2473 if (mvdHit != nullptr && mvdHit->GetStationNr() == mvdStationNum) {
2474 mvdV.emplace_back(mvdHit->GetX(), mvdHit->GetY(), iHit);
2475 }
2476 }
2477
2478 for (size_t iC = 0; iC < fCands.size(); iC++) {
2479 if (fCands[iC].IsCutTill(ELmvmAnaStep::ElId)) {
2480 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(fCands[iC].fStsInd));
2481 if (stsTrack == nullptr) continue;
2482 for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
2483 CbmMvdHit* candHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
2484 if (candHit != nullptr && candHit->GetStationNr() == mvdStationNum) {
2485 candV.emplace_back(candHit->GetX(), candHit->GetY(), iC);
2486 }
2487 }
2488 }
2489 }
2490
2491 for (size_t iC = 0; iC < candV.size(); iC++) {
2492 LmvmCand& cand = fCands[candV[iC].fInd];
2493 double minD = 9999999.;
2494 int minMvdInd = -1;
2495 for (size_t iH = 0; iH < mvdV.size(); iH++) {
2496 double d2 = LmvmUtils::Distance2(mvdV[iH].fX, mvdV[iH].fY, candV[iC].fX, candV[iC].fY);
2497 if (d2 < 1.e-9) continue;
2498 if (d2 < minD) {
2499 minMvdInd = mvdV[iH].fInd;
2500 minD = d2;
2501 }
2502 }
2503 double dmvd = sqrt(minD);
2504
2505 // Check MVD cut quality
2506 double bin = -1.;
2507 const CbmMatch* hitMatch = static_cast<const CbmMatch*>(fMvdHitMatches->At(minMvdInd));
2508 if (hitMatch != nullptr) {
2509 CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(fMCTracks->At(hitMatch->GetMatchedLink().GetIndex()));
2510 int mcMvdHitPdg = TMath::Abs(mct1->GetPdgCode());
2511 int mvdMotherId = mct1->GetMotherId();
2512
2513 int stsMotherId = -2;
2514 if (cand.fStsMcTrackId >= 0) {
2515 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
2516 stsMotherId = (mct2 != nullptr) ? mct2->GetMotherId() : -2;
2517 }
2518
2519 bin = (mvdMotherId != -1 && mvdMotherId == stsMotherId) ? 0.5 : 1.5; // correct or wrong assignment
2520 if (cand.IsMcSignal()) {
2521 bin = (mvdMotherId == stsMotherId && mcMvdHitPdg == 11) ? 0.5 : 1.5; // correct or wrong assignment
2522 }
2523 }
2524
2525 // Fill histograms
2526 double w = cand.fWeight; // changed on 10.03.25
2527 fH.FillH1(histQa, cand.fMcSrc, bin, w);
2528 fH.FillH2(hist, cand.fMcSrc, dmvd, cand.fMomentum.Mag(), w);
2529
2530 // Apply MVD cut
2531 bool isMvdCut = fCuts.IsMvdCutOk(mvdStationNum, dmvd, cand.fMomentum.Mag());
2532 if (mvdStationNum == 1)
2533 cand.fIsMvd1Cut = isMvdCut;
2534 else if (mvdStationNum == 2)
2535 cand.fIsMvd2Cut = isMvdCut;
2536 }
2537}
2538
2539void LmvmTask::MvdCutMcDistance()
2540{
2541 if (!fUseMvd) return;
2542 for (const auto& cand : fCands) {
2543 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue;
2544 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
2545 if (stsTrack == nullptr) continue;
2546 for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
2547 CbmMvdHit* mvdHit1 = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
2548 if (mvdHit1 == nullptr) continue;
2549
2550 int nofMvdHits = fMvdHitMatches->GetEntriesFast();
2551 for (int iMvd = 0; iMvd < nofMvdHits; iMvd++) {
2552 const CbmMatch* hitMatch = static_cast<const CbmMatch*>(fMvdHitMatches->At(iMvd));
2553 if (hitMatch == nullptr) continue;
2554 if (cand.fStsMcTrackId != hitMatch->GetMatchedLink().GetIndex()) continue;
2555 CbmMvdHit* mvdHit2 = static_cast<CbmMvdHit*>(fMvdHits->At(iMvd));
2556 if (mvdHit2 == nullptr || mvdHit2->GetStationNr() != mvdHit1->GetStationNr()) continue;
2557 double d = LmvmUtils::Distance(mvdHit1->GetX(), mvdHit1->GetY(), mvdHit2->GetX(), mvdHit2->GetY());
2558 double w = cand.fWeight; // changed on 10.03.25
2559 if (mvdHit1->GetStationNr() == 1) {
2560 fH.FillH1("hMvdMcDist_1", cand.fMcSrc, d, w);
2561 }
2562 else if (mvdHit1->GetStationNr() == 2) {
2563 fH.FillH1("hMvdMcDist_2", cand.fMcSrc, d, w);
2564 }
2565 }
2566 }
2567 }
2568}
2569
2570void LmvmTask::Finish()
2571{
2572 if (fCandsTotal.size() >= 1)
2573 CombinatorialPairs(); // 'if..' to prevent analysis macro not finishing in case of empty vector
2574 else
2575 LOG(warn) << "LmvmTask::CombinatorialPairs: fCandsTotal.size() = 0!";
2576
2577 fH.FillH1("hJobNumber", 0.5);
2578 TDirectory* oldir = gDirectory;
2579 TFile* outFile = FairRootManager::Instance()->GetOutFile();
2580 if (outFile != nullptr) {
2581 outFile->cd();
2582 fH.WriteToFile();
2583 }
2584
2585 // store 'vector<LmvmCand>' for large event mixing depths in ROOT file
2586 string lmvmcandfilename =
2587 "/lustre/cbm/users/criesen/data/lmvm/" + fParticle + "/lmvmcand." + fTaskId + ".root"; // path on lustre
2588 //string lmvmcandfilename = "/home/cornelius/soft/data/output/" + fParticle + "/lmvmcand." + fTaskId + ".root"; // local path
2589 std::unique_ptr<TFile> lmvmCandFile(TFile::Open(lmvmcandfilename.c_str(), "RECREATE"));
2590 lmvmCandFile->WriteObject(&fCandsTotal, "fCandsTotal");
2591
2592 gDirectory->cd(oldir->GetPath());
2593 LOG(info) << "Finished LmvmTask.";
2594}
2595
2596void LmvmTask::SetEnergyAndPlutoParticle(const string& energy, const string& particle)
2597{
2598 this->SetWeight(LmvmSimParam::GetWeight(energy, particle));
2599 fParticle = particle;
2600}
ECbmModuleId::kMvd
@ kMvd
Micro-Vertex Detector.
Definition CbmDefs.h:47
ECbmModuleId::kTrd
@ kTrd
Transition Radiation Detector.
Definition CbmDefs.h:51
ECbmModuleId::kTof
@ kTof
Time-of-flight Detector.
Definition CbmDefs.h:52
ECbmModuleId::kSts
@ kSts
Silicon Tracking System.
Definition CbmDefs.h:48
ECbmModuleId::kRich
@ kRich
Ring-Imaging Cherenkov Detector.
Definition CbmDefs.h:49
nTofHits
Int_t nTofHits
Definition CbmHadronAnalysis.cxx:72
nStsHits
Int_t nStsHits
Definition CbmHadronAnalysis.cxx:74
CbmLitMCTrackCreator.h
Creates CbmLitMCTrack objects.
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsTrack.h
Data class for STS tracks.
CbmTrdHit.h
Class for hits in TRD detector.
richProj
TClonesArray * richProj
Definition Compute_distance.h:22
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
v
fscal v[fmask::Size]
Definition KfSimdPseudo.h:4
ELmvmSrc
ELmvmSrc
Definition LmvmDef.h:23
ELmvmSrc::Signal
@ Signal
Definition LmvmDef.h:24
ELmvmSrc::Gamma
@ Gamma
Definition LmvmDef.h:27
ELmvmBgPairSrc
ELmvmBgPairSrc
Definition LmvmDef.h:51
ELmvmAnaStep
ELmvmAnaStep
Definition LmvmDef.h:34
ELmvmTopologyCut
ELmvmTopologyCut
Definition LmvmDef.h:15
ClassImp
ClassImp(LmvmTask)
Int_t
int Int_t
Definition RootTypesDef.h:16
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition CbmDigiManager.h:126
CbmGlobalTrack::GetStsTrackIndex
int32_t GetStsTrackIndex() const
Definition CbmGlobalTrack.h:44
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmGlobalTrack::GetTofTrackIndex
int32_t GetTofTrackIndex() const
Definition CbmGlobalTrack.h:49
CbmGlobalTrack::GetTofHitIndex
int32_t GetTofHitIndex() const
Definition CbmGlobalTrack.h:48
CbmGlobalTrack::GetLength
double GetLength() const
Definition CbmGlobalTrack.h:59
CbmGlobalTrack::GetTrdTrackIndex
int32_t GetTrdTrackIndex() const
Definition CbmGlobalTrack.h:45
CbmHit::GetTime
double GetTime() const
Definition CbmHit.h:79
CbmHit::GetRefId
int32_t GetRefId() const
Definition CbmHit.h:76
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition CbmL1PFFitter.cxx:506
CbmLitGlobalElectronId::IsTrdElectron
Bool_t IsTrdElectron(Int_t globalTrackindex, Double_t momentum)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:109
CbmLitGlobalElectronId::GetTofM2
Double_t GetTofM2(Int_t globalTrackIndex, Double_t momentum, Double_t eventTime=1000.)
Return TOF m2 value.
Definition CbmLitGlobalElectronId.cxx:152
CbmLitGlobalElectronId::IsTofElectron
Bool_t IsTofElectron(Int_t globalTrackIndex, Double_t momentum, Double_t eventTime=1000.)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:122
CbmLitGlobalElectronId::IsRichElectron
Bool_t IsRichElectron(Int_t globalTrackIndex, Double_t momentum)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:81
CbmLitGlobalElectronId::GetInstance
static CbmLitGlobalElectronId & GetInstance()
Definition CbmLitGlobalElectronId.h:33
CbmLitMCTrackCreator::Instance
static CbmLitMCTrackCreator * Instance()
Singleton instance.
Definition CbmLitMCTrackCreator.cxx:68
CbmLitMCTrackCreator::GetTrack
const CbmLitMCTrack & GetTrack(int mcEventId, int mcId) const
Return CbmLitMCTrack by its index.
Definition CbmLitMCTrackCreator.h:79
CbmLitMCTrack
Monte-Carlo track.
Definition CbmLitMCTrack.h:35
CbmLitMCTrack::GetNofPointsInDifferentStations
UInt_t GetNofPointsInDifferentStations(ECbmModuleId detId) const
Return number of MC points in different stations for specified detector id.
Definition CbmLitMCTrack.h:121
CbmMCTrack::GetPz
double GetPz() const
Definition CbmMCTrack.h:71
CbmMCTrack::GetCharge
double GetCharge() const
Charge of the associated particle.
Definition CbmMCTrack.cxx:150
CbmMCTrack::GetPt
double GetPt() const
Definition CbmMCTrack.h:96
CbmMCTrack::GetP
double GetP() const
Definition CbmMCTrack.h:97
CbmMCTrack::GetPx
double GetPx() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetGeantProcessId
uint32_t GetGeantProcessId() const
Definition CbmMCTrack.h:66
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:67
CbmMCTrack::GetPy
double GetPy() const
Definition CbmMCTrack.h:70
CbmMCTrack::GetRapidity
double GetRapidity() const
Definition CbmMCTrack.cxx:187
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:175
CbmMCTrack::GetNPoints
int32_t GetNPoints(ECbmModuleId detId) const
Definition CbmMCTrack.cxx:197
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
CbmMatch::GetLinks
const std::vector< CbmLink > & GetLinks() const
Definition CbmMatch.h:40
CbmMvdHit::GetStationNr
virtual int32_t GetStationNr() const
Definition CbmMvdHit.h:64
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
CbmRichRing::GetAaxis
double GetAaxis() const
Definition CbmRichRing.h:80
CbmRichRing::GetNofHits
int32_t GetNofHits() const
Definition CbmRichRing.h:37
CbmRichRing::GetCenterX
float GetCenterX() const
Definition CbmRichRing.h:77
CbmRichRing::GetChi2
double GetChi2() const
Definition CbmRichRing.h:92
CbmRichRing::GetNDF
double GetNDF() const
Definition CbmRichRing.h:93
CbmRichRing::GetCenterY
float GetCenterY() const
Definition CbmRichRing.h:78
CbmRichRing::GetBaxis
double GetBaxis() const
Definition CbmRichRing.h:81
CbmRichUtil::GetRingTrackDistance
static double GetRingTrackDistance(int globalTrackId)
Definition CbmRichUtil.cxx:24
CbmStsTrack::GetNofMvdHits
int32_t GetNofMvdHits() const
Definition CbmStsTrack.h:87
CbmStsTrack::GetMvdHitIndex
int32_t GetMvdHitIndex(int32_t iHit) const
Definition CbmStsTrack.h:75
CbmStsTrack::GetELoss
float GetELoss() const
Definition CbmStsTrack.h:121
CbmStsTrack::GetNofStsHits
int32_t GetNofStsHits() const
Definition CbmStsTrack.h:93
CbmTofTrack::GetTofHitIndex
int32_t GetTofHitIndex() const
Definition CbmTofTrack.h:55
CbmTofTrack::GetDistance
double GetDistance() const
Definition CbmTofTrack.h:79
CbmTofTrack::GetNofTofHits
int32_t GetNofTofHits() const
Definition CbmTofTrack.h:111
CbmTrack::GetNDF
int32_t GetNDF() const
Definition CbmTrack.h:64
CbmTrack::GetNofHits
virtual int32_t GetNofHits() const
Definition CbmTrack.h:58
CbmTrack::GetChiSq
double GetChiSq() const
Definition CbmTrack.h:63
CbmTrdTrack::GetPidLikeEL
double GetPidLikeEL() const
Definition CbmTrdTrack.h:43
CbmTrdTrack::GetPidLikePI
double GetPidLikePI() const
Definition CbmTrdTrack.h:44
LmvmCand::fChi2Rich
double fChi2Rich
Definition LmvmCand.h:66
LmvmCand::fChi2Tof
double fChi2Tof
Definition LmvmCand.h:68
LmvmCand::IsMcPi0
bool IsMcPi0() const
Definition LmvmCand.h:123
LmvmCand::fNofHitsMvd
int fNofHitsMvd
Definition LmvmCand.h:92
LmvmCand::fMcPdg
int fMcPdg
Definition LmvmCand.h:78
LmvmCand::fELossSts
double fELossSts
Definition LmvmCand.h:101
LmvmCand::fNofHitsRich
int fNofHitsRich
Definition LmvmCand.h:93
LmvmCand::fRichInd
int fRichInd
Definition LmvmCand.h:86
LmvmCand::fTrdLikeEl
double fTrdLikeEl
Definition LmvmCand.h:99
LmvmCand::fChi2Sts
double fChi2Sts
Definition LmvmCand.h:65
LmvmCand::fGTrackInd
int fGTrackInd
Definition LmvmCand.h:90
LmvmCand::fStsMcTrackId
int fStsMcTrackId
Definition LmvmCand.h:81
LmvmCand::fCharge
int fCharge
Definition LmvmCand.h:63
LmvmCand::fTofMcTrackId
int fTofMcTrackId
Definition LmvmCand.h:84
LmvmCand::fEventNumber
int fEventNumber
Definition LmvmCand.h:80
LmvmCand::IsMcSignal
bool IsMcSignal() const
Definition LmvmCand.h:122
LmvmCand::fTofHitInd
int fTofHitInd
Definition LmvmCand.h:88
LmvmCand::fTrdLikePi
double fTrdLikePi
Definition LmvmCand.h:100
LmvmCand::fNofHitsSts
int fNofHitsSts
Definition LmvmCand.h:91
LmvmCand::fIsTofElectron
bool fIsTofElectron
Definition LmvmCand.h:76
LmvmCand::fTofTrackInd
int fTofTrackInd
Definition LmvmCand.h:89
LmvmCand::fRapidity
double fRapidity
Definition LmvmCand.h:62
LmvmCand::fTaskId
std::string fTaskId
Definition LmvmCand.h:117
LmvmCand::fChi2Trd
double fChi2Trd
Definition LmvmCand.h:67
LmvmCand::fIsRec
bool fIsRec
Definition LmvmCand.h:104
LmvmCand::fIsElectron
bool fIsElectron
Definition LmvmCand.h:107
LmvmCand::fIsChi2Prim
bool fIsChi2Prim
Definition LmvmCand.h:106
LmvmCand::fTofDist
double fTofDist
Definition LmvmCand.h:71
LmvmCand::fMass2
double fMass2
Definition LmvmCand.h:98
LmvmCand::fRichAnn
double fRichAnn
Definition LmvmCand.h:96
LmvmCand::fTrdMcTrackId
int fTrdMcTrackId
Definition LmvmCand.h:83
LmvmCand::fIsAcc
bool fIsAcc
Definition LmvmCand.h:105
LmvmCand::fIsMvd2Cut
bool fIsMvd2Cut
Definition LmvmCand.h:110
LmvmCand::fStsInd
int fStsInd
Definition LmvmCand.h:85
LmvmCand::fMcMotherId
int fMcMotherId
Definition LmvmCand.h:79
LmvmCand::fRichMcTrackId
int fRichMcTrackId
Definition LmvmCand.h:82
LmvmCand::fTime
double fTime
Definition LmvmCand.h:70
LmvmCand::IsMcGamma
bool IsMcGamma() const
Definition LmvmCand.h:124
LmvmCand::fNofHitsTof
int fNofHitsTof
Definition LmvmCand.h:95
LmvmCand::fMcSrc
ELmvmSrc fMcSrc
Definition LmvmCand.h:121
LmvmCand::fTrdInd
int fTrdInd
Definition LmvmCand.h:87
LmvmCand::fIsMvd1Cut
bool fIsMvd1Cut
Definition LmvmCand.h:109
LmvmCand::fWeight
double fWeight
Definition LmvmCand.h:60
LmvmCand::fChi2Prim
double fChi2Prim
Definition LmvmCand.h:64
LmvmCand::fMomentum
TVector3 fMomentum
Definition LmvmCand.h:57
LmvmCand::fIsPtCut
bool fIsPtCut
Definition LmvmCand.h:114
LmvmCand::fNofHitsTrd
int fNofHitsTrd
Definition LmvmCand.h:94
LmvmCand::fLength
double fLength
Definition LmvmCand.h:69
LmvmCand::IsCutTill
bool IsCutTill(ELmvmAnaStep step) const
Definition LmvmCand.h:38
LmvmKinePar::fAngle
Double_t fAngle
Definition LmvmKinePar.h:21
LmvmKinePar::fMinv
Double_t fMinv
Definition LmvmKinePar.h:20
LmvmKinePar::fPt
Double_t fPt
Definition LmvmKinePar.h:18
LmvmKinePar::fMomentumMag
Double_t fMomentumMag
Definition LmvmKinePar.h:17
LmvmKinePar::fRapidity
Double_t fRapidity
Definition LmvmKinePar.h:19
LmvmKinePar::Create
static LmvmKinePar Create(const CbmMCTrack *mctrackP, const CbmMCTrack *mctrackM)
Definition LmvmKinePar.h:26
LmvmSimParam::GetWeight
static Double_t GetWeight(const std::string &energy, const std::string &particle)
Definition LmvmSimParam.h:15
LmvmTask::fRichHits
TClonesArray * fRichHits
Definition LmvmTask.h:180
LmvmTask::CombinatorialPairs
void CombinatorialPairs()
Definition LmvmTask.cxx:1449
LmvmTask::CheckGammaConvAndPi0
void CheckGammaConvAndPi0()
Definition LmvmTask.cxx:2276
LmvmTask::AnalyseCandidates
void AnalyseCandidates()
Definition LmvmTask.cxx:1996
LmvmTask::~LmvmTask
virtual ~LmvmTask()
Definition LmvmTask.cxx:60
LmvmTask::fNofHitsInRingMap
std::map< int, int > fNofHitsInRingMap
Definition LmvmTask.h:222
LmvmTask::fDigiManager
CbmDigiManager * fDigiManager
Definition LmvmTask.h:196
LmvmTask::InitHists
void InitHists()
Definition LmvmTask.cxx:63
LmvmTask::RichPmtXY
void RichPmtXY()
Definition LmvmTask.cxx:933
LmvmTask::fParticle
std::string fParticle
Definition LmvmTask.h:231
LmvmTask::fRichRingMatches
TClonesArray * fRichRingMatches
Definition LmvmTask.h:179
LmvmTask::fCbmEvents
TClonesArray * fCbmEvents
Definition LmvmTask.h:174
LmvmTask::fNofMinHitsRich
int fNofMinHitsRich
Definition LmvmTask.h:227
LmvmTask::fMCEventHeader
FairMCEventHeader * fMCEventHeader
Definition LmvmTask.h:173
LmvmTask::fTaskId
std::string fTaskId
Definition LmvmTask.h:232
LmvmTask::fSTCands
std::vector< LmvmCand > fSTCands
Definition LmvmTask.h:207
LmvmTask::fRichProj
TClonesArray * fRichProj
Definition LmvmTask.h:177
LmvmTask::AssignMcToTopologyCands
void AssignMcToTopologyCands(std::vector< LmvmCand > &topoCands)
Definition LmvmTask.cxx:1635
LmvmTask::RatioMomentum
void RatioMomentum(const CbmMCTrack *mct, const LmvmCand &cand, ELmvmAnaStep step, int pdg)
Definition LmvmTask.cxx:2264
LmvmTask::fMCTracks
TClonesArray * fMCTracks
Definition LmvmTask.h:175
LmvmTask::fMvdHitMatches
TClonesArray * fMvdHitMatches
Definition LmvmTask.h:187
LmvmTask::fRichRings
TClonesArray * fRichRings
Definition LmvmTask.h:176
LmvmTask::DoMcTrack
void DoMcTrack()
Definition LmvmTask.cxx:834
LmvmTask::fPionMisidLevel
double fPionMisidLevel
Definition LmvmTask.h:215
LmvmTask::fCuts
LmvmCuts fCuts
Definition LmvmTask.h:218
LmvmTask::fTofPoints
TClonesArray * fTofPoints
Definition LmvmTask.h:193
LmvmTask::fTofHitsMatches
TClonesArray * fTofHitsMatches
Definition LmvmTask.h:192
LmvmTask::CheckLikeSignCorrelations
void CheckLikeSignCorrelations()
Definition LmvmTask.cxx:500
LmvmTask::SetWeight
void SetWeight(double w)
Definition LmvmTask.h:236
LmvmTask::fCandsTotal
std::vector< LmvmCand > fCandsTotal
Definition LmvmTask.h:203
LmvmTask::fMvdPoints
TClonesArray * fMvdPoints
Definition LmvmTask.h:186
LmvmTask::FillPairHists
void FillPairHists(const LmvmCand &candP, const LmvmCand &candM, const LmvmKinePar &parMc, const LmvmKinePar &parRec, ELmvmAnaStep step)
Definition LmvmTask.cxx:1801
LmvmTask::fEventNumber
Int_t fEventNumber
Definition LmvmTask.h:217
LmvmTask::AnalyseRichRings
void AnalyseRichRings()
Definition LmvmTask.cxx:726
LmvmTask::CalculateNofTopologyPairs
void CalculateNofTopologyPairs(const std::string &name, ELmvmSrc src)
Definition LmvmTask.cxx:2362
LmvmTask::FillCandsForEventMix
void FillCandsForEventMix()
Definition LmvmTask.cxx:443
LmvmTask::MinvScale
double MinvScale(const CbmMCTrack *mct, const std::string &signal)
Definition LmvmTask.cxx:482
LmvmTask::DoMcPair
void DoMcPair()
Definition LmvmTask.cxx:903
LmvmTask::fNofMinHitsStsMvd
int fNofMinHitsStsMvd
Definition LmvmTask.h:226
LmvmTask::fKFVertex
CbmKFVertex fKFVertex
Definition LmvmTask.h:197
LmvmTask::SetEnergyAndPlutoParticle
void SetEnergyAndPlutoParticle(const std::string &energy, const std::string &particle)
Definition LmvmTask.cxx:2596
LmvmTask::fRichPoints
TClonesArray * fRichPoints
Definition LmvmTask.h:178
LmvmTask::IsPrimary
bool IsPrimary(double vertexZ)
Definition LmvmTask.h:163
LmvmTask::fTTCands
std::vector< LmvmCand > fTTCands
Definition LmvmTask.h:209
LmvmTask::fCands
std::vector< LmvmCand > fCands
Definition LmvmTask.h:202
LmvmTask::FillCands
void FillCands()
Definition LmvmTask.cxx:1348
LmvmTask::MvdCutMcDistance
void MvdCutMcDistance()
Definition LmvmTask.cxx:2539
LmvmTask::TrackSource
void TrackSource(const LmvmCand &cand, ELmvmAnaStep step, int pdg)
Definition LmvmTask.cxx:1736
LmvmTask::fMvdHits
TClonesArray * fMvdHits
Definition LmvmTask.h:185
LmvmTask::fUseMvd
bool fUseMvd
Definition LmvmTask.h:200
LmvmTask::GetMcTrackSts
CbmMCTrack * GetMcTrackSts(int stsIndex)
Definition LmvmTask.cxx:1986
LmvmTask::AssignMcToCands
void AssignMcToCands(std::vector< LmvmCand > &cands)
Definition LmvmTask.cxx:1576
LmvmTask::AnalyseGlobalTracks
void AnalyseGlobalTracks()
Definition LmvmTask.cxx:994
LmvmTask::fStsHits
TClonesArray * fStsHits
Definition LmvmTask.h:184
LmvmTask::fH
LmvmHist fH
Definition LmvmTask.h:220
LmvmTask::CheckClosestMvdHit
void CheckClosestMvdHit(int mvdStationNum, const std::string &hist, const std::string &histQa)
Definition LmvmTask.cxx:2466
LmvmTask::FillSourceHistos
void FillSourceHistos(const LmvmCand &cand)
Definition LmvmTask.cxx:2424
LmvmTask::fTofTracks
TClonesArray * fTofTracks
Definition LmvmTask.h:194
LmvmTask::AnalyseProperties
void AnalyseProperties()
Definition LmvmTask.cxx:2181
LmvmTask::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition LmvmTask.h:183
LmvmTask::Init
virtual InitStatus Init()
Definition LmvmTask.cxx:354
LmvmTask::fNofMinHitsTof
int fNofMinHitsTof
Definition LmvmTask.h:229
LmvmTask::IsRecoTrackAccepted
bool IsRecoTrackAccepted(const CbmGlobalTrack *gTrack)
Definition LmvmTask.cxx:966
LmvmTask::Exec
virtual void Exec(Option_t *option)
Definition LmvmTask.cxx:392
LmvmTask::fTofHits
TClonesArray * fTofHits
Definition LmvmTask.h:191
LmvmTask::fTrdTracks
TClonesArray * fTrdTracks
Definition LmvmTask.h:188
LmvmTask::GetPidString
std::string GetPidString(const CbmMCTrack *mct, const LmvmCand *cand)
Definition LmvmTask.cxx:1913
LmvmTask::CheckTofId
void CheckTofId(const CbmMCTrack *mcTrack, const LmvmCand &cand, ELmvmAnaStep step, int pdg)
Definition LmvmTask.cxx:2251
LmvmTask::InitOrFatal
T * InitOrFatal(const std::string &name)
Definition LmvmTask.h:57
LmvmTask::fKFFitter
CbmStsKFTrackFitter fKFFitter
Definition LmvmTask.h:198
LmvmTask::fNofMinHitsTrd
int fNofMinHitsTrd
Definition LmvmTask.h:228
LmvmTask::fTrdTrackMatches
TClonesArray * fTrdTrackMatches
Definition LmvmTask.h:190
LmvmTask::IsInBox
bool IsInBox(double boxX, double boxY, double xVal, double yVal, double width)
Definition LmvmTask.cxx:2176
LmvmTask::CheckTopologyCut
void CheckTopologyCut(ELmvmTopologyCut cut, const std::string &name)
Definition LmvmTask.cxx:2293
LmvmTask::fGlobalTracks
TClonesArray * fGlobalTracks
Definition LmvmTask.h:181
LmvmTask::FillTopologyCands
void FillTopologyCands()
Definition LmvmTask.cxx:1275
LmvmTask::fW
double fW
Definition LmvmTask.h:213
LmvmTask::fZ
double fZ
Definition LmvmTask.h:224
LmvmTask::fRTCands
std::vector< LmvmCand > fRTCands
Definition LmvmTask.h:211
LmvmTask::fPrimVertex
CbmVertex * fPrimVertex
Definition LmvmTask.h:195
LmvmTask::FillMomHists
void FillMomHists(const CbmMCTrack *mct, const LmvmCand *cand, ELmvmSrc src, ELmvmAnaStep step)
Definition LmvmTask.cxx:804
LmvmTask::Finish
virtual void Finish()
Definition LmvmTask.cxx:2570
LmvmTask::PidVsMom
void PidVsMom(const CbmGlobalTrack *gTrack, int iGTrack, int pdg, double mom, bool isAcc)
Definition LmvmTask.cxx:1104
LmvmTask::FillHistosForFastSim
void FillHistosForFastSim()
Definition LmvmTask.cxx:549
LmvmTask::FillRichRingNofHits
void FillRichRingNofHits()
Definition LmvmTask.cxx:452
LmvmTask::PairSource
void PairSource(const LmvmCand &candP, const LmvmCand &candM, ELmvmAnaStep step, const LmvmKinePar &parRec)
Definition LmvmTask.cxx:1648
LmvmTask::IsInAllDets
bool IsInAllDets(const CbmGlobalTrack *gTrack)
Definition LmvmTask.cxx:1253
LmvmTask::CheckMismatches
void CheckMismatches(const CbmGlobalTrack *gTrack, int pdg, bool isElectron, const std::string &ptcl, double weight)
Definition LmvmTask.cxx:1143
LmvmTask::fStsTracks
TClonesArray * fStsTracks
Definition LmvmTask.h:182
LmvmUtils::IsTofElectron
static void IsTofElectron(int globalTrackIndex, double momentum, LmvmCand *cand)
Definition LmvmUtils.cxx:266
LmvmUtils::IsMcGammaEl
static bool IsMcGammaEl(const CbmMCTrack *mct, TClonesArray *mcTracks)
Definition LmvmUtils.cxx:93
LmvmUtils::Distance
static double Distance(double x1, double y1, double x2, double y2)
Definition LmvmUtils.cxx:233
LmvmUtils::GetMcPairSrc
static ELmvmSrc GetMcPairSrc(const CbmMCTrack *mctP, const CbmMCTrack *mctM, TClonesArray *mcTracks)
Definition LmvmUtils.cxx:145
LmvmUtils::IsElectronMc
static void IsElectronMc(LmvmCand *cand, TClonesArray *mcTracks, double pionMisidLevel)
Definition LmvmUtils.cxx:271
LmvmUtils::GetWeightPair
static double GetWeightPair(const LmvmCand &cand1, const LmvmCand &cand2)
Definition LmvmUtils.cxx:300
LmvmUtils::IsElectron
static void IsElectron(int globalTrackIndex, double momentum, double momentumCut, LmvmCand *cand)
Definition LmvmUtils.cxx:240
LmvmUtils::GetMcSrc
static ELmvmSrc GetMcSrc(CbmMCTrack *mctrack, TClonesArray *mcTracks)
Definition LmvmUtils.cxx:78
LmvmUtils::IsTrdElectron
static void IsTrdElectron(int globalTrackIndex, double momentum, LmvmCand *cand)
Definition LmvmUtils.cxx:261
LmvmUtils::MinvScale
static double MinvScale(double minv, const std::string &particle)
Definition LmvmUtils.cxx:312
LmvmUtils::IsRichElectron
static void IsRichElectron(int globalTrackIndex, double momentum, LmvmCand *cand)
Definition LmvmUtils.cxx:256
LmvmUtils::IsMismatch
static bool IsMismatch(const LmvmCand &cand)
Definition LmvmUtils.cxx:218
LmvmUtils::IsGhost
static bool IsGhost(const LmvmCand &cand)
Definition LmvmUtils.cxx:226
LmvmUtils::GetChargeStr
static std::string GetChargeStr(const LmvmCand *cand)
Definition LmvmUtils.cxx:288
LmvmUtils::GetBgPairSrc
static ELmvmBgPairSrc GetBgPairSrc(const LmvmCand &candP, const LmvmCand &candM)
Definition LmvmUtils.cxx:197
LmvmUtils::IsMcSignalEl
static bool IsMcSignalEl(const CbmMCTrack *mct)
Definition LmvmUtils.cxx:87
LmvmUtils::CalculateAndSetTrackParams
static void CalculateAndSetTrackParams(LmvmCand *cand, CbmStsTrack *stsTrack, CbmKFVertex &kfVertex)
Definition LmvmUtils.cxx:28
LmvmUtils::Distance2
static double Distance2(double x1, double y1, double x2, double y2)
Definition LmvmUtils.cxx:235
Cbm::NumberToString
std::string NumberToString(const T &value, int precision=1)
Definition CbmUtils.h:34
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:22