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] */
4
5#include "LmvmTask.h"
6
7#include "CbmGlobalTrack.h"
8#include "CbmKF.h"
9#include "CbmL1PFFitter.h"
10#include "CbmMCTrack.h"
11#include "CbmMatch.h"
12#include "CbmMvdHit.h"
13#include "CbmRichHit.h"
14#include "CbmRichPoint.h"
15#include "CbmRichRing.h"
16#include "CbmRichUtil.h"
17#include "CbmStsHit.h"
18#include "CbmStsTrack.h"
19#include "CbmTofHit.h"
20#include "CbmTofPoint.h"
21#include "CbmTofTrack.h"
22#include "CbmTrackMatchNew.h"
23#include "CbmTrdHit.h"
24#include "CbmTrdTrack.h"
25#include "CbmVertex.h"
26#include "cbm/elid/CbmLitGlobalElectronId.h"
27#include "cbm/qa/mc/CbmLitMCTrackCreator.h"
28
29#include "FairEventHeader.h"
30#include "FairMCPoint.h"
31#include "FairRootManager.h"
32#include "FairRunAna.h"
33#include "FairTask.h"
34#include "FairTrackParam.h"
35
36#include "TClonesArray.h"
37#include "TDatabasePDG.h"
38#include "TFile.h"
39#include "TMCProcess.h"
40#include "TRandom3.h"
41#include "TVector3.h"
42
43#include <sstream>
44#include <vector>
45
46#include "LmvmHist.h"
47#include "LmvmSimParam.h"
48#include "LmvmUtils.h"
49
50using namespace std;
51
52ClassImp(LmvmTask);
53
54
55LmvmTask::LmvmTask() : FairTask("LmvmTask") {}
56
57
58LmvmTask::~LmvmTask() {}
59
60
61void LmvmTask::InitHists()
62{
63 string ax = "Yield";
64 string axMinv = "dN/dM_{ee}/N [GeV/c^{2}]^{-1}";
65
66 fH.CreateH2("hMomVsAnglePairSignalMc", "#sqrt{P_{e^{#pm}} P_{e^{#mp}}} [GeV/c]", "#theta_{e^{+},e^{-}} [deg]",
67 "Counter", 100, 0., 5., 1000, 0., 50.);
68
69 fH.CreateH1("hMotherPdg", {"mc", "acc"}, "Pdg code", "Particles / Event", 7000, -3500., 3500.);
70 fH.CreateH1("hCandPdg", fH.fAnaStepNames, "Pdg code", "Particles / Event", 7001, -3500., 3500.);
71 fH.CreateH2("hCandPdgVsMom", fH.fAnaStepNames, "P [GeV/c]", "Particle", "Yield/(Event * Bin)", 200, 0., 10., 6, 0.,
72 6.);
73 fH.CreateH2("hCandElSrc", "Analysis step", "Mother of Electron Candidate", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps,
74 8, 0., 8.);
75 fH.CreateH2("hBgPairPdg", fH.fAnaStepNames, "PDG of Candidate 1", "PDG of Candidate 2", ax, 8, 0., 8., 8, 0., 8.);
76
77 fH.CreateH2("hPmtXY", fH.fSrcNames, "X [cm]", "Y [cm]", "Counter", 110, -110, 110, 200, -200, 200);
78
79 fH.CreateH2("hVertexGammaXZ", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 200, -10., 190., 400, -130., 130.);
80 fH.CreateH2("hVertexGammaYZ", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 200, -10., 190., 400, -130., 130.);
81 fH.CreateH2("hVertexGammaXY", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 400, -130., 130., 400, -130., 130.);
82 fH.CreateH2("hVertexGammaRZ", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 50,
83 0., 50.);
84
85 // vertices of not-electron candidates that are misidentified as electrons in ToF
86 fH.CreateH2("hVertexXZ_misidTof", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.);
87 fH.CreateH2("hVertexYZ_misidTof", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.);
88 fH.CreateH2("hVertexXY_misidTof", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 100, -50., 50., 100, -50., 50.);
89 fH.CreateH2("hVertexRZ_misidTof", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100.,
90 100, 0., 50.);
91 fH.CreateH1("hNofBgTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
92 fH.CreateH1("hNofSignalTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
93 fH.CreateH2("hBgSrcTracks", "Analysis step", "Candidate Source", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.);
94
95 fH.CreateH1("hNofTopoPairs", {"gamma", "pi0"}, "Pair type", "Pairs/event", 8, 0., 8);
96 fH.CreateH1("hNofMismatches", {"all", "rich", "trd", "tof"}, "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0.,
97 fH.fNofAnaSteps);
98 fH.CreateH1("hNofMismatches_gTracks", {"all", "complete"}, fH.fGTrackNames, "Detector", ax, 7., 0., 7.);
99 fH.CreateH1("hNofMismatchedTrackSegments", {"all", "complete"}, fH.fGTrackNames, "nof mism. track segments", ax, 4.,
100 0., 4.); //"all": also partial tracks; "complete": only tracks with entries in all detectors
101 fH.CreateH2("hMatchId_gTracks", fH.fGTrackNames, "Nof mismatched Track Segments", "Identification", ax, 4., 0., 4.,
102 2., 0., 2.);
103
104 fH.CreateH1("hChi2_mismatch_all", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.);
105 fH.CreateH1("hChi2_truematch_all", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.);
106 fH.CreateH1("hChi2_mismatch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.);
107 fH.CreateH1("hChi2_truematch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.);
108
109 fH.CreateH1("hNofGhosts", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps);
110
111 fH.CreateH2("hSrcBgPairs", "Analysis step", "Pair", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, fH.fNofBgPairSrc, 0.,
112 fH.fNofBgPairSrc);
113 fH.CreateH2("hSrcBgPairsEpEm", fH.fAnaStepNames, "mother particle e+", "mother particle e-", ax, 4, 0., 4., 4, 0.,
114 4.);
115
116 fH.CreateH1("hEventNumber", "", "", 1, 0, 1.);
117 fH.CreateH1("hEventNumberMixed", "", "", 1, 0, 1.);
118
119 fH.CreateH1("hAnnRich", fH.fSrcNames, "RICH ANN output", ax, 100, -1.1, 1.1);
120 fH.CreateH2("hAnnRichVsMom", fH.fSrcNames, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 10., 100, -1.1, 1.1);
121 fH.CreateH1("hAnnTrd", fH.fSrcNames, "Likelihood output", ax, 100, -.1, 1.1); // TODO: change back to "TRD ANN"
122 fH.CreateH2("hTofM2", fH.fSrcNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 150, 0., 6., 500, -0.1, 1.0);
123
124 fH.CreateH1("hChi2Sts", fH.fSrcNames, "#chi^{2}", ax, 200, 0., 20.);
125 fH.CreateH1("hChi2PrimVertex", fH.fSrcNames, "#chi^{2}_{prim}", ax, 200, 0., 20.);
126 fH.CreateH1("hNofMvdHits", fH.fSrcNames, "Number of hits in MVD", ax, 5, -0.5, 4.5);
127 fH.CreateH1("hNofStsHits", fH.fSrcNames, "Number of hits in STS", ax, 9, -0.5, 8.5);
128 fH.CreateH2("hTrdLike", {"El", "Pi"}, fH.fSrcNames, "P [GeV/c]", "Likelihood output", ax, 100, 0., 6., 100, 0., 1.);
129
130 fH.CreateH2("hAnnRichVsMomPur", {"El", "Bg"}, "P [GeV/c]", "RICH ANN output", ax, 100, 0., 10., 100, -1.1, 1.1);
131 fH.CreateH2("hTrdElLikePur", {"El", "Bg"}, "P [GeV/c]", "Likelihood output", ax, 100, 0., 10., 100, 0., 1.);
132
133 fH.CreateH2("hTtCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
134 "#theta_{e^{+},e^{-}} [deg]", ax, 100, 0., 5., 100, 0., 5.);
135 fH.CreateH2("hStCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
136 "#theta_{e^{#pm},rec} [deg]", ax, 100, 0., 5., 100, 0., 5.);
137 fH.CreateH2("hRtCut", {"all", "pion", "truePair"}, fH.fSrcNames, "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
138 "#theta_{e^{#pm},rec} [deg]", ax, 100, 0., 5., 100, 0., 5.);
139
140 fH.CreateH2("hMvdCut", {"1", "2"}, fH.fSrcNames, "d_{MVD} [cm]", "P_{e} [GeV/c]", ax, 100, 0., 1., 100, 0., 5.);
141 fH.CreateH2("hMvdXY", {"1", "2"}, fH.fSrcNames, "X [cm]", "Y [cm]", ax, 60, -6., 6., 60, -6., 6.);
142 fH.CreateH1("hMvdR", {"1", "2"}, fH.fSrcNames, "#sqrt{X^{2}+Y^{2}} [cm]", ax, 60, 0., 6.);
143 fH.CreateH1("hMvdCutQa", {"1", "2"}, fH.fSrcNames, "MVD hit assignment", ax, 2, 0.,
144 2.); // [0.5]-correct, [1.5]-wrong
145 fH.CreateH1("hMvdMcDist", {"1", "2"}, fH.fSrcNames, "Track-Hit distance [cm]", ax, 100, 0., 10.);
146
147 fH.CreateH1("hMinv", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5);
148 fH.CreateH1("hMinv_urqmdAll", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0.,
149 2.5); // for UrQMD particles only
150 fH.CreateH1("hMinv_urqmdEl", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0.,
151 2.5); // for UrQMD electrons only
152
153 // Pair Yield histograms for combinatorial BG calculation
154 for (std::string comb : {"PM", "PP", "MM"}) {
155 for (std::string cat : {"", "_urqmdAll", "_urqmdEl"}) {
156 string hName = "hMinvComb" + comb + cat;
157 fH.CreateH1(hName, {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5);
158 }
159 }
160
161 fH.CreateH1("hMinvBgMatch", {"trueMatch", "trueMatchEl", "trueMatchNotEl", "mismatch"}, fH.fAnaStepNames,
162 "M_{ee} [GeV/c^{2}]", ax, 250, 0., 2.5);
163 fH.CreateH1("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5);
164 fH.CreateH1("hMinvBgSource2_elid", {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"},
165 "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); // "pi" are misid. charged pions
166
167 fH.CreateH2("hMinvPt", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", ax, 100, 0., 2., 25, 0.,
168 2.5);
169
170 fH.CreateH1("hMomPairSignal", fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 15.);
171 fH.CreateH2("hPtYPairSignal", fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
172 fH.CreateH1("hAnglePair", fH.fSrcNames, fH.fAnaStepNames, "#Theta_{1,2} [deg]", ax, 160, 0., 80.);
173
174 for (std::string suff : {"", "+", "-"}) {
175 fH.CreateH1("hMom" + suff, fH.fSrcNames, fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 10.);
176 fH.CreateH1("hMomPx" + suff, fH.fSrcNames, fH.fAnaStepNames, "Px [GeV/c]", ax, 100, -3., 3.);
177 fH.CreateH1("hMomPy" + suff, fH.fSrcNames, fH.fAnaStepNames, "Py [GeV/c]", ax, 100, -3., 3.);
178 fH.CreateH1("hMomPz" + suff, fH.fSrcNames, fH.fAnaStepNames, "Pz [GeV/c]", ax, 120, -1., 11.);
179 fH.CreateH1("hPt" + suff, fH.fSrcNames, fH.fAnaStepNames, "P_{t} [GeV/c]", ax, 100, 0., 4.);
180 fH.CreateH1("hRapidity" + suff, fH.fSrcNames, fH.fAnaStepNames, "Rapidity", ax, 100, 0., 5.);
181 }
182
183 fH.CreateH1("hMomAcc+", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.);
184 fH.CreateH1("hMomAcc-", {"sts", "rich", "trd", "tof"}, fH.fSrcNames, "P [GeV/c]", ax, 100, 0., 10.);
185
186 fH.CreateH1("hElMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
187 "P [GeV/c]", ax, 100, 0., 10.);
188 fH.CreateH1("hPiMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
189 "P [GeV/c]", ax, 100, 0., 10.);
190 fH.CreateH1("hProtonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
191 "P [GeV/c]", ax, 100, 0., 10.);
192 fH.CreateH1("hKaonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"},
193 "P [GeV/c]", ax, 100, 0., 10.);
194
195 fH.CreateH1("hPionSupp", {"pi", "idEl"}, fH.fAnaStepNames, "P [GeV/c]", ax, 20, 0., 10.);
196
197 fH.CreateH1("hMom_gTracks", fH.fGTrackNames, "P [GeV/c]", ax, 100, 0., 10.);
198 fH.CreateH1("hMom_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 10.);
199 fH.CreateH2("hPtY_gTracks", fH.fGTrackNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
200 fH.CreateH2("hPtY_cands", fH.fCandNames, fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
201 fH.CreateH2("hTofM2_gTracks", fH.fGTrackNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.);
202 fH.CreateH2("hTofM2_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10.,
203 500, -0.1, 1.);
204
205 fH.CreateH1("hTofPilePdgs_cands", fH.fAnaStepNames, "Particle", ax, 12, 0., 12.);
206 fH.CreateH1("hTofPilePdgs_gTracks", "Particle", ax, 12, 0., 12.);
207 fH.CreateH2("hTofPileHitXY", fH.fCandNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.);
208 fH.CreateH2("hTofPilePointXY", fH.fCandNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.);
209 fH.CreateH1("hTofPileHitPointDist", fH.fCandNames, "#sqrt{dX^{2} + dY^{2}} [cm]", ax, 1000., 0., 20.);
210 fH.CreateH2("hTofPilePty_cands", fH.fCandNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.);
211
212 // ToF Hits
213 fH.CreateH2("hTofPointXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax,
214 110., -550., 550., 90., -450.,
215 450.); // to see if maybe misid/mismatches stem from a certain region in ToF
216 fH.CreateH2("hTofHitXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, 110.,
217 -550., 550., 90., -450., 450.);
218 fH.CreateH1("hTofHitPointDist", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames,
219 "#sqrt{dX^{2} + dY^{2}} [cm]", ax, 400., 0., 20.);
220 fH.CreateH2("hTofHitTrackDist_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]",
221 "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.);
222 fH.CreateH2("hTofHitTrackDist_cands", fH.fCandNames, "P [GeV/c]", "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.);
223
224 fH.CreateH1("hMomRatio_cands", fH.fCandNames, fH.fAnaStepNames, "Ratio P_{rec}/P_{MC}", ax, 100, 0., 2.);
225 fH.CreateH2("hMomRatioVsMom_cands", fH.fCandNames, fH.fAnaStepNames, "P_{MC} [GeV/c]", "Ratio P_{rec}/P_{MC}", ax,
226 100, 0., 10., 200., 0., 2.);
227
228 // compare misid. candidates (pions, proton, ...) with not misid. ones
229 for (size_t iP = 4; iP < fH.fCandNames.size(); iP++) {
230 fH.CreateH1("hMom_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", ax, 100, 0., 10.);
231 fH.CreateH2("hPtY_" + fH.fCandNames[iP], {"true", "misid"}, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 40, 0.,
232 2.);
233 fH.CreateH2("hTofM2_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0.,
234 10., 500, -0.1, 1.);
235 }
236
237 // beta-momentum spectra
238 fH.CreateH2("hBetaMom", {"gTracks", "cands"}, fH.fCandNames, "P * Q ", "beta", ax, 200, -10., 10., 150., 0., 1.5);
239 fH.CreateH2("hBetaMom_allGTracks", "P * Q ", "beta", ax, 200, -10., 10., 150., 0., 1.5);
240
241 // z vertex of global tracks
242 fH.CreateH2("hVertexGTrackRZ", "Z [cm]", "R [cm]", ax, 1500, -50., 100., 1000, -50., 50.);
243
244 /*fH.CreateH2("hIndexStsRich", fH.fGTrackNames, "STS Index", "RICH Index", ax, 400., 0., 400., 100., 0., 100.); // TODO: can be deleted
245 fH.CreateH2("hIndexStsTrd", fH.fGTrackNames, "STS Index", "TRD Index", ax, 400., 0., 400., 400., 0., 400.);
246 fH.CreateH2("hIndexStsTof", fH.fGTrackNames, "STS Index", "ToF Index", ax, 400., 0., 400., 800., 0., 800.);*/
247
248 fH.CreateH2("hPdgVsMom_gTracks", {"rich", "trd", "tof"}, {"all", "complete"}, "P [GeV/c]", "Misid. Particle",
249 "Yield/(Event * Bin)", 200, 0., 10., 6, 0., 6.);
250
251 fH.CreateH2("hTofM2Calc_gTracks", fH.fGTrackNames, "Time", "", "Yield", 100., 5., 15., 4., 0., 4.);
252 fH.CreateH2("hTofTimeVsMom_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]",
253 "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., 15.);
254 fH.CreateH2("hTofTimeVsMom_cands", fH.fCandNames, "P [GeV/c]", "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5.,
255 15.);
256 fH.CreateH2("hRichRingTrackDist_gTracks", fH.fGTrackNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., 10., 200., 0.,
257 20.);
258 fH.CreateH2("hRichRingTrackDist_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0.,
259 10., 200., 0., 20.);
260
261 // check chi2 of candidates after El-ID step
262 fH.CreateH2("hChi2VsMom", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "P [GeV/c]", "#chi^{2}", ax, 100., 0., 10.,
263 200, 0., 20.);
264 fH.CreateH2("hTofTimeVsChi2", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "#chi^{2}", "Time (ct) [m]", ax, 200., 0.,
265 20., 200, 0., 20.);
266 fH.CreateH2("hChi2Comb", {"StsRich", "StsTrd", "StsTof", "RichTrd", "RichTof", "TrdTof"}, fH.fCandNames,
267 "#chi^{2}_{1}", "#chi^{2}_{2}", ax, 200., 0., 20., 200, 0., 20.);
268}
269
270InitStatus LmvmTask::Init()
271{
272 fMCEventHeader = InitOrFatal<FairMCEventHeader>("MCEventHeader.");
273 fMCTracks = InitOrFatal<TClonesArray>("MCTrack");
274 fRichHits = InitOrFatal<TClonesArray>("RichHit");
275 fRichRings = InitOrFatal<TClonesArray>("RichRing");
276 fRichPoints = InitOrFatal<TClonesArray>("RichPoint");
277 fRichRingMatches = InitOrFatal<TClonesArray>("RichRingMatch");
278 fRichProj = InitOrFatal<TClonesArray>("RichProjection");
279 fStsTrackMatches = InitOrFatal<TClonesArray>("StsTrackMatch");
280 fStsTracks = InitOrFatal<TClonesArray>("StsTrack");
281 fStsHits = InitOrFatal<TClonesArray>("StsHit");
282 if (fUseMvd) {
283 fMvdHits = InitOrFatal<TClonesArray>("MvdHit");
284 fMvdPoints = InitOrFatal<TClonesArray>("MvdPoint");
285 fMvdHitMatches = InitOrFatal<TClonesArray>("MvdHitMatch");
286 }
287 fGlobalTracks = InitOrFatal<TClonesArray>("GlobalTrack");
288 fTrdTracks = InitOrFatal<TClonesArray>("TrdTrack");
289 fTrdTrackMatches = InitOrFatal<TClonesArray>("TrdTrackMatch");
290 fTofPoints = InitOrFatal<TClonesArray>("TofPoint");
291 fTofHits = InitOrFatal<TClonesArray>("TofHit");
292 fTofHitsMatches = InitOrFatal<TClonesArray>("TofHitMatch");
293 fPrimVertex = InitOrFatal<CbmVertex>("PrimaryVertex.");
294
295 //fTofTracks = InitOrFatal<TClonesArray>("TofTracks"); // TODO: check this and next lines
296 FairRootManager* fairRootMan = FairRootManager::Instance();
297 fTofTracks = (TClonesArray*) fairRootMan->GetObject("TofTrack");
298 if (fTofTracks == nullptr) { LOG(warning) << "LmvmTask::Init : no TofTrack array!"; }
299
300 InitHists();
301
302 CbmLitMCTrackCreator::Instance();
303 CbmLitGlobalElectronId::GetInstance();
304
305 return kSUCCESS;
306}
307
308void LmvmTask::Exec(Option_t*)
309{
310 fH.FillH1("hEventNumber", 0.5);
311 fEventNumber++;
312 // bool useMbias = false; // false for 40% central agag collisions (b<7.7fm)
313 // bool isCentralCollision = false;
314
315 // if (!useMbias) {
316 // double impactPar = fMCEventHeader->GetB();
317 // if (impactPar <= 7.7) isCentralCollision = true;
318 // }
319
320 cout << "========================================================" << endl;
321 LOG(info) << "LmvmTask event number " << fEventNumber;
322 LOG(info) << "fPionMisidLevel = " << fPionMisidLevel;
323 LOG(info) << fCuts.ToString();
324 LOG(info) << "fW = " << fW;
325
326 if (fPrimVertex != nullptr) { fKFVertex = CbmKFVertex(*fPrimVertex); }
327 else {
328 Fatal("LmvmTask::Exec", "No PrimaryVertex array!");
329 }
330
331 //if (useMbias || (!useMbias && isCentralCollision)) {
332 FillRichRingNofHits();
333 DoMcTrack();
334 DoMcPair();
335 RichPmtXY();
336 FillTopologyCands();
337 FillCands();
338 CalculateNofTopologyPairs("hNofTopoPairs_gamma", ELmvmSrc::Gamma);
339 CalculateNofTopologyPairs("hNofTopoPairs_pi0", ELmvmSrc::Pi0);
340 AnalyseGlobalTracks();
341 AnalyseCandidates();
342
343 fCandsTotal.insert(fCandsTotal.end(), fCands.begin(), fCands.end());
344 LOG(info) << "fCandsTotal.size = " << fCandsTotal.size();
345 //}
346} // Exec
347
348void LmvmTask::FillRichRingNofHits()
349{
350 fNofHitsInRingMap.clear();
351 int nofRichHits = fRichHits->GetEntriesFast();
352 for (int iHit = 0; iHit < nofRichHits; iHit++) {
353 CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
354 if (hit == nullptr || hit->GetRefId() < 0) continue;
355 FairMCPoint* point = static_cast<FairMCPoint*>(fRichPoints->At(hit->GetRefId()));
356 if (point == nullptr) continue;
357 CbmMCTrack* track = static_cast<CbmMCTrack*>(fMCTracks->At(point->GetTrackID()));
358 if (track == nullptr || track->GetMotherId() < 0) continue;
359 fNofHitsInRingMap[track->GetMotherId()]++;
360 }
361}
362
363double LmvmTask::MinvScale(double minv)
364{
365 double scale = -1.;
366 if (fParticle == "inmed") scale = LmvmUtils::GetMassScaleInmed(minv);
367 else if (fParticle == "qgp")
368 scale = LmvmUtils::GetMassScaleQgp(minv);
369 else
370 scale = 1.;
371 return scale;
372}
373
374
375void LmvmTask::FillMomHists(const CbmMCTrack* mct, const LmvmCand* cand, ELmvmSrc src, ELmvmAnaStep step)
376{
377 if ((mct != nullptr && cand != nullptr) || (mct == nullptr && cand == nullptr)) {
378 LOG(error) << "LmvmTask::FillMomHists: Both mct and cand are [not nullptr] or [nullptr].";
379 return;
380 }
381 bool isMc = (mct != nullptr);
382 string chargeStr = (isMc) ? LmvmUtils::GetChargeStr(mct) : LmvmUtils::GetChargeStr(cand);
383 double w = ((mct != nullptr && LmvmUtils::IsMcSignalEl(mct)) || (cand != nullptr && cand->IsMcSignal())) ? fW : 1.;
384
385 for (const string& suff : {string(""), chargeStr}) {
386 if (suff == "0") continue;
387 fH.FillH1("hMom" + suff, src, step, (isMc) ? mct->GetP() : cand->fMomentum.Mag(), w);
388 fH.FillH1("hMomPx" + suff, src, step, (isMc) ? mct->GetPx() : cand->fMomentum.X(), w);
389 fH.FillH1("hMomPy" + suff, src, step, (isMc) ? mct->GetPy() : cand->fMomentum.Y(), w);
390 fH.FillH1("hMomPz" + suff, src, step, (isMc) ? mct->GetPz() : cand->fMomentum.Z(), w);
391 fH.FillH1("hPt" + suff, src, step, (isMc) ? mct->GetPt() : cand->fMomentum.Perp(), w);
392 fH.FillH1("hRapidity" + suff, src, step, (isMc) ? mct->GetRapidity() : cand->fRapidity, w);
393 }
394}
395
396void LmvmTask::DoMcTrack()
397{
398 int nMcTracks = fMCTracks->GetEntriesFast();
399 LOG(info) << "nMcTracks = " << nMcTracks;
400
401 for (int i = 0; i < nMcTracks; i++) {
402 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i));
403 if (mct == nullptr) continue;
404 ELmvmSrc src = LmvmUtils::GetMcSrc(mct, fMCTracks);
405 string chargeStr = (mct->GetCharge() > 0) ? "+" : "-";
406 bool isAcc = IsMcTrackAccepted(i);
407 double mom = mct->GetP();
408 double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.;
409
410 FillMomHists(mct, nullptr, src, ELmvmAnaStep::Mc);
411 if (isAcc) FillMomHists(mct, nullptr, src, ELmvmAnaStep::Acc);
412
413 if (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4)
414 fH.FillH1("hMomAcc" + chargeStr + "_sts", src, mom, w);
415 if (fNofHitsInRingMap[i] >= 7) fH.FillH1("hMomAcc" + chargeStr + "_rich", src, mom, w);
416 if (mct->GetNPoints(ECbmModuleId::kTrd) >= 2) fH.FillH1("hMomAcc" + chargeStr + "_trd", src, mom, w);
417 if (mct->GetNPoints(ECbmModuleId::kTof) >= 1) fH.FillH1("hMomAcc" + chargeStr + "_tof", src, mom, w);
418
419 if (LmvmUtils::IsMcGammaEl(mct, fMCTracks)) {
420 TVector3 v;
421 mct->GetStartVertex(v);
422 for (const auto step : {ELmvmAnaStep::Mc, ELmvmAnaStep::Acc}) {
423 if (step == ELmvmAnaStep::Acc && !isAcc) continue;
424 fH.FillH2("hVertexGammaXZ", step, v.Z(), v.X());
425 fH.FillH2("hVertexGammaYZ", step, v.Z(), v.Y());
426 fH.FillH2("hVertexGammaXY", step, v.X(), v.Y());
427 fH.FillH2("hVertexGammaRZ", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
428 }
429 }
430
431 // Fill PDG histos
432 int mcPdg = mct->GetPdgCode();
433 if (std::abs(mcPdg) == 11 || mcPdg == 99009911) {
434 int mcMotherPdg = 0;
435 if (mct->GetMotherId() != -1) {
436 CbmMCTrack* mother = static_cast<CbmMCTrack*>(fMCTracks->At(mct->GetMotherId()));
437 if (mother != nullptr) mcMotherPdg = mother->GetPdgCode();
438 }
439 if (std::abs(mcPdg) == 11) {
440 fH.FillH1("hMotherPdg_mc", mcMotherPdg);
441 if (isAcc) fH.FillH1("hMotherPdg_acc", mcMotherPdg);
442 }
443 }
444
445 if (std::abs(mcPdg) == 211) fH.FillH1("hPionSupp_pi", ELmvmAnaStep::Mc, mom, w);
446 if (std::abs(mcPdg) == 211 && IsMcTrackAccepted(i)) fH.FillH1("hPionSupp_pi", ELmvmAnaStep::Acc, mom, w);
447 }
448}
449
450void LmvmTask::DoMcPair()
451{
452 int nMcTracks = fMCTracks->GetEntries();
453 for (int iMc1 = 0; iMc1 < nMcTracks; iMc1++) { // TODO: range until iMc1 < nMcTracks-1 ??
454 CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc1));
455 ELmvmSrc src = LmvmUtils::GetMcSrc(mct1, fMCTracks);
456
457 // To speed up: select only signal, eta and pi0 electrons
458 if (!(src == ELmvmSrc::Signal || src == ELmvmSrc::Pi0 || src == ELmvmSrc::Eta)) continue;
459
460 bool isAcc1 = IsMcTrackAccepted(iMc1);
461 for (int iMc2 = iMc1 + 1; iMc2 < nMcTracks; iMc2++) {
462 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2));
463 bool isAccPair = isAcc1 && IsMcTrackAccepted(iMc2);
464 ELmvmSrc srcPair = LmvmUtils::GetMcPairSrc(mct1, mct2, fMCTracks);
465 LmvmKinePar pKin = LmvmKinePar::Create(mct1, mct2);
466
467 if (srcPair == ELmvmSrc::Signal) {
468 fH.FillH2("hMomVsAnglePairSignalMc", std::sqrt(mct1->GetP() * mct2->GetP()), pKin.fAngle);
469 }
470
471 for (const auto step : {ELmvmAnaStep::Mc, ELmvmAnaStep::Acc}) {
472 if (step == ELmvmAnaStep::Acc && !isAccPair) continue;
473 //fH.FillH1("hAnglePair", srcPair, step, pKin.fAngle, w);
474 if (srcPair == ELmvmSrc::Signal) {
475 fH.FillH2("hPtYPairSignal", step, pKin.fRapidity, pKin.fPt, fW);
476 fH.FillH1("hMomPairSignal", step, pKin.fMomentumMag, fW);
477 }
478
479 // MC and Acc minv only for signal, eta and pi0
480 if (srcPair == ELmvmSrc::Signal || srcPair == ELmvmSrc::Pi0 || srcPair == ELmvmSrc::Eta) {
481 double w = 1;
482 if (srcPair == ELmvmSrc::Signal) w = fW * MinvScale(pKin.fMinv);
483 fH.FillH1("hMinv", srcPair, step, pKin.fMinv, w);
484 }
485 }
486 }
487 }
488}
489
490void LmvmTask::RichPmtXY()
491{
492 int nofRichHits = fRichHits->GetEntriesFast();
493 for (int iH = 0; iH < nofRichHits; iH++) {
494 CbmRichHit* richHit = static_cast<CbmRichHit*>(fRichHits->At(iH));
495 if (richHit == nullptr || richHit->GetRefId() < 0) continue;
496 FairMCPoint* pointPhoton = static_cast<FairMCPoint*>(fRichPoints->At(richHit->GetRefId()));
497 if (pointPhoton == nullptr) continue;
498 CbmMCTrack* trackPhoton = static_cast<CbmMCTrack*>(fMCTracks->At(pointPhoton->GetTrackID()));
499 if (trackPhoton == nullptr || trackPhoton->GetMotherId() < 0) continue;
500 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(trackPhoton->GetMotherId()));
501 if (mct == nullptr) continue;
502
503 TVector3 v;
504 mct->GetStartVertex(v);
505 ELmvmSrc src = LmvmUtils::GetMcSrc(mct, fMCTracks);
506 double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.;
507 if (v.Z() < 2.) { fH.FillH2("hPmtXY", src, richHit->GetX(), richHit->GetY(), w); }
508 }
509}
510
511bool LmvmTask::IsMcTrackAccepted(int mcTrackInd)
512{
513 CbmMCTrack* tr = static_cast<CbmMCTrack*>(fMCTracks->At(mcTrackInd));
514 if (tr == nullptr) return false;
515 int nRichPoints = fNofHitsInRingMap[mcTrackInd];
516 return (tr->GetNPoints(ECbmModuleId::kMvd) + tr->GetNPoints(ECbmModuleId::kSts) >= 4 && nRichPoints >= 7
517 && tr->GetNPoints(ECbmModuleId::kTrd) >= 2 && tr->GetNPoints(ECbmModuleId::kTof) > 1);
518}
519
520void LmvmTask::AnalyseGlobalTracks()
521{
522 int nofMcTracks = fMCTracks->GetEntriesFast();
523
524 for (int i = 0; i < nofMcTracks; i++) {
525 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i));
526 int pdg = mct->GetPdgCode();
527 bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4);
528 TVector3 vertex;
529 mct->GetStartVertex(vertex);
530
531 if (std::abs(pdg) == 11 || std::abs(pdg) == 211 || pdg == 2212 || std::abs(pdg) == 321) {
532 string hName = (std::abs(pdg) == 11) ? "hEl" : (std::abs(pdg) == 211) ? "hPi" : pdg == 2212 ? "hProton" : "hKaon";
533 fH.FillH1(hName + "Mom_all_mc", mct->GetP());
534 if (isAccSts) fH.FillH1(hName + "Mom_all_acc", mct->GetP());
535
536 if (vertex.Mag() <= 0.1) {
537 fH.FillH1(hName + "Mom_prim_mc", mct->GetP());
538 if (isAccSts) fH.FillH1(hName + "Mom_prim_acc", mct->GetP());
539 }
540 }
541 } // MC tracks
542
543 int ngTracks = fGlobalTracks->GetEntriesFast();
544 for (int i = 0; i < ngTracks; i++) {
545 CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(i));
546 if (gTrack == nullptr) continue;
547 int stsInd = gTrack->GetStsTrackIndex();
548 bool isRich = (gTrack->GetRichRingIndex() >= 0);
549 bool isTrd = (gTrack->GetTrdTrackIndex() >= 0);
550 bool isTof = (gTrack->GetTofHitIndex() >= 0);
551
552 if (stsInd < 0) continue;
553 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
554 if (stsTrack == nullptr) continue;
555 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
556 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
557 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
558 if (stsMcTrackId < 0) continue;
559 CbmMCTrack* mct = (CbmMCTrack*) fMCTracks->At(stsMcTrackId);
560 if (mct == nullptr) continue;
561
562 double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.;
563
564 int pdg = mct->GetPdgCode();
565 double p = mct->GetP();
566 double pt = mct->GetPt();
567 double rap = mct->GetRapidity();
568 double m2 = CbmLitGlobalElectronId::GetInstance().GetTofM2(i, p);
569
570 TVector3 v;
571 mct->GetStartVertex(v);
572 bool isPrim = IsPrimary(v.Mag());
573 if (std::abs(pdg) == 11 || std::abs(pdg) == 211 || pdg == 2212 || std::abs(pdg) == 321) {
574 string hName = (std::abs(pdg) == 11) ? "hEl" : (std::abs(pdg) == 211) ? "hPi" : pdg == 2212 ? "hProton" : "hKaon";
575 fH.FillH1(hName + "Mom_all_recSts", p);
576 if (isRich) fH.FillH1(hName + "Mom_all_recStsRich", p);
577 if (isRich && isTrd) fH.FillH1(hName + "Mom_all_recStsRichTrd", p);
578 if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_all_recStsRichTrdTof", p);
579
580 if (isPrim) {
581 fH.FillH1(hName + "Mom_prim_recSts", p);
582 if (isRich) fH.FillH1(hName + "Mom_prim_recStsRich", p);
583 if (isRich && isTrd) fH.FillH1(hName + "Mom_prim_recStsRichTrd", p);
584 if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_prim_recStsRichTrdTof", p);
585 }
586 }
587
588 string ptcl = GetPidString(v.Mag(), pdg);
589
590 fH.FillH2("hVertexGTrackRZ", v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
591 /*fH.FillH2("hIndexStsRich_" + ptcl, stsInd, gTrack->GetRichRingIndex());
592 fH.FillH2("hIndexStsTrd_" + ptcl, stsInd, gTrack->GetTrdTrackIndex());
593 fH.FillH2("hIndexStsTof_" + ptcl, stsInd, gTrack->GetTofHitIndex());*/
594
595 bool isTofEl = (CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) ? true : false;
596 bool isElectron = (CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, p)
597 && CbmLitGlobalElectronId::GetInstance().IsTrdElectron(i, p)
598 && CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p))
599 ? true
600 : false;
601
602 if (gTrack != nullptr) { CheckMismatches(gTrack, pdg, isElectron, ptcl, w); }
603 if (std::abs(pdg) != 11) PidVsMom(gTrack, i, pdg, p);
604
605 Double_t richDist = CbmRichUtil::GetRingTrackDistance(i);
606 fH.FillH2("hRichRingTrackDist_gTracks_" + ptcl, p, richDist);
607
608 // investigate misidentifications: compare misidentified candidates with same not-misidentified particles
609 // first check if global track has entries in all detectors
610 if (!IsInAllDets(gTrack)) continue; // Mind: all following actions are done only for fully rec. gTracks
611
612 CheckTofIdentification(gTrack, ptcl, p, m2, pdg, isTofEl);
613
614 fH.FillH1("hMom_gTracks_" + ptcl, p);
615 fH.FillH2("hPtY_gTracks_" + ptcl, rap, pt);
616 fH.FillH2("hTofM2_gTracks_" + ptcl, p, m2);
617
618 // check PIDs in "Tof pile"
619 if (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01) {
620 bool isSignal = (mct != nullptr && mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11);
621 double pdgBin = (pdg == 11 && isSignal) ? 0.5
622 : (pdg == -11 && isSignal) ? 1.5
623 : (pdg == 11 && !isSignal && isPrim) ? 2.5
624 : (pdg == -11 && !isSignal && isPrim) ? 3.5
625 : (pdg == 11 && !isSignal && !isPrim) ? 4.5
626 : (pdg == -11 && !isSignal && !isPrim) ? 5.5
627 : (pdg == 211) ? 6.5
628 : (pdg == -211) ? 7.5
629 : (pdg == 2212) ? 8.5
630 : (pdg == 321) ? 9.5
631 : (pdg == -321) ? 10.5
632 : 11.5;
633 fH.FillH1("hTofPilePdgs_gTracks", pdgBin);
634 }
635
636 string ptcl2 = GetPidString(mct, nullptr);
637 if (!isElectron && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) {
638 fH.FillH1("hMom_" + ptcl2 + "_true", p);
639 fH.FillH2("hPtY_" + ptcl2 + "_true", rap, pt);
640 fH.FillH2("hTofM2_" + ptcl2 + "_true", p, m2);
641 }
642 else if (isElectron
643 && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) {
644 fH.FillH1("hMom_" + ptcl2 + "_misid", p);
645 fH.FillH2("hPtY_" + ptcl2 + "_misid", rap, pt);
646 fH.FillH2("hTofM2_" + ptcl2 + "_misid", p, m2);
647 }
648
649 BetaMom(mct, gTrack, ptcl2);
650 }
651}
652
653bool LmvmTask::IsInTofPile(double p, double m2) { return (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01); }
654
655void LmvmTask::PidVsMom(const CbmGlobalTrack* gTrack, int iGTrack, int pdg, double p)
656{
657 bool isRichEl = CbmLitGlobalElectronId::GetInstance().IsRichElectron(iGTrack, p);
658 bool isTrdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iGTrack, p);
659 bool isTofEl = CbmLitGlobalElectronId::GetInstance().IsTofElectron(iGTrack, p);
660
661 double pdgBin = (pdg == 211) ? 0.5
662 : (pdg == -211) ? 1.5
663 : (pdg == 2212) ? 2.5
664 : (pdg == 321) ? 3.5
665 : (pdg == -321) ? 4.5
666 : 5.5;
667
668 if (isRichEl) fH.FillH2("hPdgVsMom_gTracks_rich_all", p, pdgBin);
669 if (isTrdEl) fH.FillH2("hPdgVsMom_gTracks_trd_all", p, pdgBin);
670 if (isTofEl) fH.FillH2("hPdgVsMom_gTracks_tof_all", p, pdgBin);
671
672 if (isRichEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_rich_complete", p, pdgBin);
673 if (isTrdEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_trd_complete", p, pdgBin);
674 if (isTofEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_tof_complete", p, pdgBin);
675}
676
677void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string& gtString, double pMc, double m2,
678 int pdg, bool isTofEl)
679{
680 // Get STS ID
681 int stsInd = gTrack->GetStsTrackIndex();
682 if (stsInd < 0) return;
683 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
684 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) return;
685 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
686 if (stsMcTrackId < 0) return;
687
688 // Get ToF Match
689 Int_t tofInd = gTrack->GetTofHitIndex();
690 if (tofInd < 0) return;
691 CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd);
692 if (tofHit == nullptr) return;
693 CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(gTrack->GetStsTrackIndex()));
694 if (tofTrack == nullptr) return;
695 CbmMatch* tofHitMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofInd));
696 if (tofHitMatch == nullptr) return;
697 int tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
698 if (tofPointIndex < 0) return;
699 FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofPointIndex));
700 if (tofPoint == nullptr) return;
701 int tofMcTrackId = tofPoint->GetTrackID();
702
703 bool isTofMismatch = (tofMcTrackId == stsMcTrackId) ? false : true;
704 bool isTofMisid = (std::abs(pdg) == 11 && isTofEl) ? false : (std::abs(pdg) != 11 && !isTofEl) ? false : true;
705
706 double pointX = tofPoint->GetX();
707 double pointY = tofPoint->GetY();
708 double hitX = tofHit->GetX();
709 double hitY = tofHit->GetY();
710 double dX = pointX - hitX;
711 double dY = pointY - hitY;
712 double dR = sqrt(dX * dX + dY * dY);
713 double dist = tofTrack->GetDistance();
714
715 string hNameId = (!isTofMisid) ? "_trueid_" + gtString : "_misid_" + gtString;
716 fH.FillH2("hTofPointXY" + hNameId, pointX, pointY);
717 fH.FillH2("hTofHitXY" + hNameId, hitX, hitY);
718 fH.FillH1("hTofHitPointDist" + hNameId, dR);
719 fH.FillH2("hTofHitTrackDist_gTracks" + hNameId, pMc, dist);
720
721 string hNameMatch = (!isTofMismatch) ? "_truematch_" + gtString : "_mismatch_" + gtString;
722 fH.FillH2("hTofPointXY" + hNameMatch, pointX, pointY);
723 fH.FillH2("hTofHitXY" + hNameMatch, hitX, hitY);
724 fH.FillH1("hTofHitPointDist" + hNameMatch, dR);
725 fH.FillH2("hTofHitTrackDist_gTracks" + hNameMatch, pMc, dist);
726
727 // check m2 calculation in ToF
728 double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime();
729 Double_t noOffsetTime = tofHit->GetTime() - eventTime;
730 Double_t t = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m
731
732 fH.FillH2("hTofTimeVsMom_gTracks" + hNameId, pMc, t);
733 fH.FillH2("hTofTimeVsMom_gTracks" + hNameMatch, pMc, t);
734
735 if (!isTofMisid) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 0.5);
736 if (isTofMisid) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 1.5);
737 if (!isTofMismatch) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 2.5);
738 if (isTofMismatch) fH.FillH2("hTofM2Calc_gTracks_" + gtString, t, 3.5);
739
740 // check ToF data for particles in ToF pile // TODO: needed or can be deleted?
741 if (isTofEl && std::abs(pdg) != 11 && IsInTofPile(pMc, m2)) {
742 string candString = (pdg == 211) ? fH.fCandNames[4]
743 : (pdg == -211) ? fH.fCandNames[5]
744 : (pdg == 2212) ? fH.fCandNames[6]
745 : (pdg == 321) ? fH.fCandNames[7]
746 : (pdg == -321) ? fH.fCandNames[8]
747 : fH.fCandNames[9];
748 fH.FillH2("hTofPileHitXY_" + candString, hitX, hitY);
749 fH.FillH2("hTofPilePointXY_" + candString, pointX, pointY);
750 fH.FillH1("hTofPileHitPointDist_" + candString, dR);
751 }
752}
753
754void LmvmTask::BetaMom(const CbmMCTrack* mct, const CbmGlobalTrack* gTrack, const string& ptcl)
755{
756 Int_t tofInd = gTrack->GetTofHitIndex();
757 if (tofInd < 0) return;
758 CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd);
759 if (NULL == tofHit) return;
760 double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime();
761 Double_t noOffsetTime = tofHit->GetTime() - eventTime;
762
763 double p = mct->GetP();
764 double l = gTrack->GetLength();
765 double t = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m
766 double q = (mct->GetCharge() > 0) ? 1 : (mct->GetCharge() < 0) ? -1. : 0.;
767 double beta = l / (t * 100); // cm to m
768 fH.FillH2("hBetaMom_gTracks_" + ptcl, q * p, beta);
769 fH.FillH2("hBetaMom_allGTracks", q * p, beta);
770}
771
772void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isElectron, const string& ptcl, double w)
773{
774 // Chi2 of true-/mismatched; number of mismatches (general and seperated for detectors)
775 int stsInd = gTrack->GetStsTrackIndex();
776 if (stsInd < 0) return;
777 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(stsInd));
778 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) return;
779 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
780 if (stsMcTrackId < 0) return;
781
782 bool isFull = IsInAllDets(gTrack);
783
784 // first calculate Chi2
785 CbmL1PFFitter fPFFitter;
786 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
787 vector<CbmStsTrack> stsTracks;
788 stsTracks.resize(1);
789 stsTracks[0] = *stsTrack;
790 vector<CbmL1PFFitter::PFFieldRegion> vField;
791 vector<float> chiPrim;
792 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
793 double chi2 = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
794
795 CbmTrackMatchNew* richMatch = nullptr;
796 CbmTrackMatchNew* trdMatch = nullptr;
797 CbmMatch* tofMatch = nullptr;
798
799 int richRingInd = gTrack->GetRichRingIndex();
800 int trdTrackInd = gTrack->GetTrdTrackIndex();
801 int tofHitInd = gTrack->GetTofHitIndex();
802
803 if (richRingInd >= 0) richMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(richRingInd));
804 if (trdTrackInd >= 0) trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(trdTrackInd));
805 if (tofHitInd >= 0) tofMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofHitInd));
806
807 int nofMismTrackSegmentsAll = 0;
808 int nofMismTrackSegmentsComp = 0;
809
810 fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 0.5, w);
811 if (isFull) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 0.5, w);
812
813 if (richMatch != nullptr) {
814 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
815 if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 1.5, w);
816 if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_rich_" + ptcl, chi2, w);
817 if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) {
818 fH.FillH1("hChi2_mismatch_all_rich_" + ptcl, chi2, w);
819 fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 2.5, w);
820 nofMismTrackSegmentsAll++;
821 }
822 if (isFull) {
823 if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 1.5, w);
824 if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId)
825 fH.FillH1("hChi2_truematch_complete_rich_" + ptcl, chi2, w);
826 if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) {
827 fH.FillH1("hChi2_mismatch_complete_rich_" + ptcl, chi2, w);
828 fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 2.5, w);
829 nofMismTrackSegmentsComp++;
830 }
831 }
832 }
833
834 if (trdMatch != nullptr) {
835 int trdMcTrackId = trdMatch->GetMatchedLink().GetIndex();
836 if (trdMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 3.5, w);
837 if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_trd_" + ptcl, chi2, w);
838 if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) {
839 fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 4.5, w);
840 fH.FillH1("hChi2_mismatch_all_trd_" + ptcl, chi2, w);
841 nofMismTrackSegmentsAll++;
842 }
843 if (isFull) {
844 if (trdMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 3.5, w);
845 if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_complete_trd_" + ptcl, chi2, w);
846 if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) {
847 fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 4.5, w);
848 fH.FillH1("hChi2_mismatch_complete_trd_" + ptcl, chi2, w);
849 nofMismTrackSegmentsComp++;
850 }
851 }
852 }
853
854 if (tofMatch != nullptr && tofMatch->GetMatchedLink().GetIndex() >= 0) {
855 FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofMatch->GetMatchedLink().GetIndex()));
856 if (tofPoint != nullptr) {
857 int tofMcTrackId = tofPoint->GetTrackID();
858 if (tofMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 5.5, w);
859 if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_tof_" + ptcl, chi2, w);
860 if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) {
861 fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 6.5, w);
862 fH.FillH1("hChi2_mismatch_all_tof_" + ptcl, chi2, w);
863 nofMismTrackSegmentsAll++;
864 }
865 if (isFull) {
866 if (tofMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 5.5, w);
867 if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId)
868 fH.FillH1("hChi2_truematch_complete_tof_" + ptcl, chi2, w);
869 if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) {
870 fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 6.5, w);
871 fH.FillH1("hChi2_mismatch_complete_tof_" + ptcl, chi2, w);
872 nofMismTrackSegmentsComp++;
873 }
874 }
875 }
876 }
877
878 fH.FillH1("hNofMismatchedTrackSegments_all_" + ptcl, nofMismTrackSegmentsAll + 0.5, w);
879 if (isFull) fH.FillH1("hNofMismatchedTrackSegments_complete_" + ptcl, nofMismTrackSegmentsComp + 0.5, w);
880
881 bool isTrueId = (std::abs(pdg) == 11 && isElectron) ? true : (std::abs(pdg) != 11 && !isElectron) ? true : false;
882 double binX = nofMismTrackSegmentsComp + 0.5;
883 double binY = (isTrueId == true) ? 0.5 : 1.5;
884
885 fH.FillH2("hMatchId_gTracks_" + ptcl, binX, binY, w);
886}
887
888bool LmvmTask::IsInAllDets(const CbmGlobalTrack* gTrack)
889{
890 if (gTrack == nullptr) return false;
891
892 int stsInd = gTrack->GetStsTrackIndex();
893 if (stsInd < 0) return false;
894 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd));
895 if (stsTrack == nullptr) return false;
896
897 int richInd = gTrack->GetRichRingIndex();
898 int trdInd = gTrack->GetTrdTrackIndex();
899 int tofInd = gTrack->GetTofHitIndex();
900 if (richInd < 0 || trdInd < 0 || tofInd < 0) return false;
901
902 CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(richInd));
903 CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(trdInd));
904 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(tofInd));
905 if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr) return false;
906
907 return true;
908}
909
910void LmvmTask::FillTopologyCands()
911{
912 fSTCands.clear();
913 fRTCands.clear();
914 int ngTracks = fGlobalTracks->GetEntriesFast();
915
916 for (int iGTrack = 0; iGTrack < ngTracks; iGTrack++) {
917 LmvmCand cand;
918
919 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGTrack);
920 if (gTrack == nullptr) continue;
921
922 cand.fStsInd = gTrack->GetStsTrackIndex();
923 if (cand.fStsInd < 0) continue;
924 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
925 if (stsTrack == nullptr) continue;
926
927 cand.fRichInd = gTrack->GetRichRingIndex();
928 cand.fTrdInd = gTrack->GetTrdTrackIndex();
929 cand.fTofInd = gTrack->GetTofHitIndex(); // TODO: is TofHitIndex; change this name everywhere
930 cand.fTofTrackInd = gTrack->GetStsTrackIndex(); // TODO: is this right??
931
932 LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex);
933 cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim);
934 if (!cand.fIsChi2Prim) continue;
935
936 // ST cut candidates, only STS
937 if (cand.fRichInd < 0 && cand.fTrdInd < 0 && cand.fTofInd < 0) fSTCands.push_back(cand);
938
939 // RT cut candidates, STS + at least one detector (RICH, TRD, TOF) but not all
940 // Candidates must be identified as electron in registered detectors:
941 // if it is registered in RICH it must be identified in the RICH as electron
942 // RICH
943 bool isRichRT = (cand.fRichInd < 0) ? false : true;
944 if (isRichRT) {
945 CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(cand.fRichInd));
946 if (richRing == nullptr) isRichRT = false;
947 if (isRichRT) isRichRT = CbmLitGlobalElectronId::GetInstance().IsRichElectron(iGTrack, cand.fMomentum.Mag());
948 }
949
950 // TRD
951 bool isTrdRT = (cand.fTrdInd < 0) ? false : true;
952 if (isTrdRT) {
953 CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(cand.fTrdInd));
954 if (trdTrack == nullptr) isTrdRT = false;
955 if (isTrdRT) isTrdRT = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iGTrack, cand.fMomentum.Mag());
956 }
957
958 // ToF
959 bool isTofRT = (cand.fTofInd < 0) ? false : true;
960 if (isTofRT) {
961 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd));
962 if (tofHit == nullptr) isTofRT = false;
963 if (isTofRT) isTofRT = CbmLitGlobalElectronId::GetInstance().IsTofElectron(iGTrack, cand.fMomentum.Mag());
964 }
965
966 if (isRichRT || isTrdRT || isTofRT) {
967 if (!(cand.fRichInd >= 0 && cand.fTrdInd >= 0 && cand.fTofInd >= 0)) { fRTCands.push_back(cand); }
968 }
969 }
970 LOG(info) << "fSTCands.size = " << fSTCands.size();
971 LOG(info) << "fRTCands.size = " << fRTCands.size();
972
973 AssignMcToTopologyCands(fSTCands);
974 AssignMcToTopologyCands(fRTCands);
975}
976
977void LmvmTask::FillCands()
978{
979 fCands.clear();
980 fTTCands.clear();
981 int nGTracks = fGlobalTracks->GetEntriesFast();
982 fCands.reserve(nGTracks);
983
984 for (int iGTrack = 0; iGTrack < nGTracks; iGTrack++) {
985 LmvmCand cand;
986 // if MVD is not used set mvd cuts to true
987 cand.fIsMvd1Cut = !fUseMvd;
988 cand.fIsMvd2Cut = !fUseMvd;
989 cand.fEventNumber = fEventNumber;
990
991 CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iGTrack));
992 if (gTrack == nullptr) continue;
993
994 // STS
995 cand.fStsInd = gTrack->GetStsTrackIndex();
996 if (cand.fStsInd < 0) continue;
997 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
998 if (stsTrack == nullptr) continue;
999
1000 // RICH - TRD - TOF
1001 cand.fRichInd = gTrack->GetRichRingIndex();
1002 cand.fTrdInd = gTrack->GetTrdTrackIndex();
1003 cand.fTofInd = gTrack->GetTofHitIndex();
1004 cand.fTofTrackInd = gTrack->GetStsTrackIndex();
1005
1006 // Set length and time
1007 cand.fLength = gTrack->GetLength() / 100.;
1008 Int_t tofInd = gTrack->GetTofHitIndex();
1009 if (tofInd >= 0) {
1010 CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd);
1011 if (NULL != tofHit) {
1012 double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime();
1013 Double_t noOffsetTime = tofHit->GetTime() - eventTime;
1014 cand.fTime = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m
1015 }
1016 }
1017
1018 LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex);
1019 cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim);
1020 cand.fIsPtCut = fCuts.IsPtCutOk(cand.fMomentum.Perp());
1021
1022 // Add all pions from STS for pion misidentification level study
1023 if (fPionMisidLevel >= 0.0) {
1024 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(cand.fStsInd);
1025 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
1026 cand.fStsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1027 if (cand.fStsMcTrackId < 0) continue;
1028 CbmMCTrack* mcTrack1 = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1029 if (mcTrack1 == nullptr) continue;
1030
1031 //check that pion track has track projection onto the photodetector plane
1032 const FairTrackParam* richProj = static_cast<FairTrackParam*>(fRichProj->At(iGTrack));
1033 if (richProj == nullptr || richProj->GetX() == 0 || richProj->GetY() == 0) continue;
1034
1035 if (std::abs(mcTrack1->GetPdgCode()) == 211) {
1036 LmvmUtils::IsElectronMc(&cand, fMCTracks, fPionMisidLevel);
1037 fCands.push_back(cand);
1038 continue;
1039 }
1040 }
1041
1042 if (cand.fRichInd < 0 || cand.fTrdInd < 0 || cand.fTofInd < 0 || cand.fTofTrackInd < 0) continue;
1043
1044 cand.fGTrackInd = iGTrack;
1045
1046 CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(cand.fRichInd));
1047 CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(cand.fTrdInd));
1048 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd));
1049 CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(cand.fTofTrackInd));
1050 if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr || tofTrack == nullptr) continue;
1051
1052 cand.fChi2Rich = (richRing->GetNDF() > 0.) ? richRing->GetChi2() / richRing->GetNDF() : 0.;
1053 cand.fChi2Trd = (trdTrack->GetNDF() > 0.) ? trdTrack->GetChiSq() / trdTrack->GetNDF() : 0.;
1054
1055 cand.fTrdLikeEl = trdTrack->GetPidLikeEL();
1056 cand.fTrdLikePi = trdTrack->GetPidLikePI();
1057
1058 cand.fTofDist = tofTrack->GetDistance();
1059
1060 LmvmUtils::IsElectron(iGTrack, cand.fMomentum.Mag(), fCuts.fMomentumCut, &cand);
1061 LmvmUtils::IsRichElectron(iGTrack, cand.fMomentum.Mag(), &cand);
1062 LmvmUtils::IsTrdElectron(iGTrack, cand.fMomentum.Mag(), &cand);
1063 LmvmUtils::IsTofElectron(iGTrack, cand.fMomentum.Mag(), &cand);
1064
1065 fCands.push_back(cand);
1066
1067 if (!cand.fIsElectron && cand.fIsChi2Prim) fTTCands.push_back(cand);
1068 }
1069 LOG(info) << "fTTCands.size = " << fTTCands.size();
1070 LOG(info) << "fCands.size = " << fCands.size();
1071
1072 AssignMcToCands(fCands);
1073 AssignMcToTopologyCands(fTTCands);
1074}
1075
1076void LmvmTask::CombinatorialPairs()
1077{
1078 size_t nCand = fCandsTotal.size();
1079 for (size_t iC1 = 0; iC1 < nCand - 1; iC1++) {
1080 const auto& cand1 = fCandsTotal[iC1];
1081
1082 for (size_t iC2 = iC1 + 1; iC2 < nCand; iC2++) {
1083 const auto& cand2 = fCandsTotal[iC2];
1084 LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2);
1085 double w = 1.;
1086 if (cand1.IsMcSignal()) w = fW * MinvScale(cand1.fMassSig);
1087 else if (cand2.IsMcSignal())
1088 w = fW * MinvScale(cand2.fMassSig);
1089 bool isSameEvent = (cand1.fEventNumber == cand2.fEventNumber);
1090
1091 // PLUTO electrons from two different events have to be scaled with w^2!
1092 if (!isSameEvent && cand1.IsMcSignal() && cand2.IsMcSignal())
1093 w = fW * fW * MinvScale(cand1.fMassSig) * MinvScale(cand2.fMassSig);
1094 for (auto step : fH.fAnaSteps) {
1095 if (cand1.IsCutTill(step) && cand2.IsCutTill(step)) {
1096
1097 // all particles
1098 if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) continue;
1099 if (cand1.fCharge * cand2.fCharge < 0) {
1100 if (isSameEvent) {
1101 fH.FillH1("hMinvCombPM_sameEv", step, pRec.fMinv, w);
1102 cout << "LmvmTask::CombinatorialPairs(): e+e- same event (all). w = " << w
1103 << endl; //TODO: delete this line
1104 }
1105 else
1106 fH.FillH1("hMinvCombPM_mixedEv", step, pRec.fMinv, w);
1107 }
1108 if (cand1.fCharge < 0 && cand2.fCharge < 0) {
1109 if (isSameEvent) {
1110 fH.FillH1("hMinvCombMM_sameEv", step, pRec.fMinv, w);
1111 cout << "LmvmTask::CombinatorialPairs(): e-e- same event (all). w = " << w
1112 << endl; //TODO: delete this line
1113 }
1114 else
1115 fH.FillH1("hMinvCombMM_mixedEv", step, pRec.fMinv, w);
1116 }
1117 if (cand1.fCharge > 0 && cand2.fCharge > 0) {
1118 if (isSameEvent) {
1119 fH.FillH1("hMinvCombPP_sameEv", step, pRec.fMinv, w);
1120 cout << "LmvmTask::CombinatorialPairs(): e+e+ same event (all). w = " << w
1121 << endl; //TODO: delete this line
1122 }
1123 else
1124 fH.FillH1("hMinvCombPP_mixedEv", step, pRec.fMinv, w);
1125 }
1126
1127 // only UrQMD particles (also misidentified)
1128 if (!cand1.IsMcSignal() && !cand2.IsMcSignal()) {
1129 if (isSameEvent) {
1130 if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdAll_sameEv", step, pRec.fMinv, w);
1131 else if (cand1.fCharge < 0 && cand2.fCharge < 0)
1132 fH.FillH1("hMinvCombMM_urqmdAll_sameEv", step, pRec.fMinv, w);
1133 else if (cand1.fCharge > 0 && cand2.fCharge > 0)
1134 fH.FillH1("hMinvCombPP_urqmdAll_sameEv", step, pRec.fMinv, w);
1135 }
1136 else {
1137 if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdAll_mixedEv", step, pRec.fMinv, w);
1138 else if (cand1.fCharge < 0 && cand2.fCharge < 0)
1139 fH.FillH1("hMinvCombMM_urqmdAll_mixedEv", step, pRec.fMinv, w);
1140 else if (cand1.fCharge > 0 && cand2.fCharge > 0)
1141 fH.FillH1("hMinvCombPP_urqmdAll_mixedEv", step, pRec.fMinv, w);
1142 }
1143 }
1144
1145 // only UrQMD electrons
1146 if (!cand1.IsMcSignal() && !cand2.IsMcSignal() && std::abs(cand1.fMcPdg) == 11
1147 && std::abs(cand2.fMcPdg) == 11) { // TODO: delete double brackets
1148 if (isSameEvent) {
1149 if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdEl_sameEv", step, pRec.fMinv, w);
1150 else if (cand1.fCharge < 0 && cand2.fCharge < 0)
1151 fH.FillH1("hMinvCombMM_urqmdEl_sameEv", step, pRec.fMinv, w);
1152 else if (cand1.fCharge > 0 && cand2.fCharge > 0)
1153 fH.FillH1("hMinvCombPP_urqmdEl_sameEv", step, pRec.fMinv, w);
1154 }
1155 else {
1156 if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdEl_mixedEv", step, pRec.fMinv, w);
1157 else if (cand1.fCharge < 0 && cand2.fCharge < 0)
1158 fH.FillH1("hMinvCombMM_urqmdEl_mixedEv", step, pRec.fMinv, w);
1159 else if (cand1.fCharge > 0 && cand2.fCharge > 0)
1160 fH.FillH1("hMinvCombPP_urqmdEl_mixedEv", step, pRec.fMinv, w);
1161 }
1162 }
1163 } // isCutTillStep
1164 } // steps
1165 } // cand2
1166 } // cand1
1167}
1168
1169void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands)
1170{
1171 for (auto& cand : cands) {
1172 cand.ResetMcParams();
1173
1174 //STS. MCTrackId of the candidate is defined by STS track
1175 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(cand.fStsInd));
1176 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
1177 cand.fStsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1178 if (cand.fStsMcTrackId < 0) continue;
1179 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1180 if (mct == nullptr) continue;
1181 cand.fMcMotherId = mct->GetMotherId();
1182 cand.fMcPdg = mct->GetPdgCode();
1183 cand.fMcSrc = LmvmUtils::GetMcSrc(mct, fMCTracks);
1184
1185 // Get mass of mother particle if it is signal (for mass-dependant scaling)
1186 if (cand.IsMcSignal()) {
1187 int nMcTracks = fMCTracks->GetEntries();
1188 for (int iMc2 = 0; iMc2 < nMcTracks; iMc2++) {
1189 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2));
1190 if (mct2->GetMotherId() == cand.fMcMotherId && iMc2 != cand.fStsMcTrackId) {
1191 LmvmKinePar pKin = LmvmKinePar::Create(mct, mct2);
1192 cand.fMassSig = pKin.fMinv;
1193 }
1194 }
1195 }
1196
1197 if (std::abs(cand.fMcPdg) == 211 && fPionMisidLevel >= 0.) continue;
1198
1199 // RICH
1200 CbmTrackMatchNew* richMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(cand.fRichInd));
1201 if (richMatch == nullptr) continue;
1202 cand.fRichMcTrackId = richMatch->GetMatchedLink().GetIndex();
1203
1204 // TRD
1205 CbmTrackMatchNew* trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(cand.fTrdInd));
1206 if (trdMatch == nullptr) continue;
1207 cand.fTrdMcTrackId = trdMatch->GetMatchedLink().GetIndex();
1208
1209 // ToF
1210 if (cand.fTofInd < 0) continue;
1211 CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd));
1212 if (tofHit == nullptr) continue;
1213 CbmMatch* tofHitMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(cand.fTofInd));
1214 if (tofHitMatch == nullptr) continue;
1215 int tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex();
1216 if (tofPointIndex < 0) continue;
1217 FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofPointIndex));
1218 if (tofPoint == nullptr) continue;
1219 cand.fTofMcTrackId = tofPoint->GetTrackID();
1220 }
1221}
1222
1223void LmvmTask::AssignMcToTopologyCands(vector<LmvmCand>& topoCands)
1224{
1225 for (auto& cand : topoCands) {
1226 cand.ResetMcParams();
1227 if (cand.fStsInd < 0) continue;
1228 CbmTrackMatchNew* stsMatch = static_cast<CbmTrackMatchNew*>(fStsTrackMatches->At(cand.fStsInd));
1229 if (stsMatch == nullptr || stsMatch->GetNofLinks() == 0) continue;
1230 cand.fStsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1231 if (cand.fStsMcTrackId < 0) continue;
1232 CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1233 if (mct == nullptr) continue;
1234
1235 cand.fMcMotherId = mct->GetMotherId();
1236 cand.fMcPdg = mct->GetPdgCode();
1237 cand.fMcSrc = LmvmUtils::GetMcSrc(mct, fMCTracks);
1238 }
1239}
1240
1241void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step, const LmvmKinePar& parRec)
1242{
1243 ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM);
1244 double w = (candP.IsMcSignal() || candM.IsMcSignal()) ? fW : 1.;
1245 fH.FillH1("hAnglePair", src, step, parRec.fAngle, w);
1246
1247 if (src == ELmvmSrc::Bg) {
1248 // gamma=0.5, pi0=1.5, pions=2.5, other=3.5
1249 double indM = candM.IsMcGamma() ? 0.5 : (candM.IsMcPi0() ? 1.5 : (std::abs(candM.fMcPdg) == 211) ? 2.5 : 3.5);
1250 double indP = candP.IsMcGamma() ? 0.5 : (candP.IsMcPi0() ? 1.5 : (std::abs(candP.fMcPdg) == 211) ? 2.5 : 3.5);
1251 fH.FillH2("hSrcBgPairsEpEm", step, indP, indM);
1252
1253 ELmvmBgPairSrc bgSrc = LmvmUtils::GetBgPairSrc(candP, candM);
1254 if (bgSrc != ELmvmBgPairSrc::Undefined) {
1255 string name = fH.GetName("hMinvBgSource_" + fH.fBgPairSrcNames[static_cast<int>(bgSrc)], step);
1256 fH.FillH1(name, parRec.fMinv);
1257 fH.FillH2("hSrcBgPairs", static_cast<int>(step) + 0.5, static_cast<double>(bgSrc) + 0.5);
1258
1259 if (step == ELmvmAnaStep::ElId) {
1260 string hName = "hMinvBgSource2_elid_";
1261
1262 if (std::abs(candP.fMcPdg) == 11) {
1263
1264 // cand1 is El and from Gamma
1265 if (candP.IsMcGamma()) {
1266 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gg", parRec.fMinv);
1267 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1268 fH.FillH1(hName + "gpi0", parRec.fMinv);
1269 else if (std::abs(candM.fMcPdg) == 211)
1270 fH.FillH1(hName + "gpi", parRec.fMinv);
1271 else
1272 fH.FillH1(hName + "go", parRec.fMinv);
1273 }
1274
1275 // cand1 is El and from Pi0
1276 else if (candP.IsMcPi0()) {
1277 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gpi0", parRec.fMinv);
1278 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1279 fH.FillH1(hName + "pi0pi0", parRec.fMinv);
1280 else if (std::abs(candM.fMcPdg) == 211)
1281 fH.FillH1(hName + "pipi0", parRec.fMinv);
1282 else
1283 fH.FillH1(hName + "pi0o", parRec.fMinv);
1284 }
1285
1286 // cand1 is El but not from Gamma or Pi0
1287 else {
1288 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "go", parRec.fMinv);
1289 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1290 fH.FillH1(hName + "pi0o", parRec.fMinv);
1291 else if (std::abs(candM.fMcPdg) == 211)
1292 fH.FillH1(hName + "pio", parRec.fMinv);
1293 else
1294 fH.FillH1(hName + "oo", parRec.fMinv);
1295 }
1296 }
1297
1298 // cand1 is misid. charged pion
1299 else if (std::abs(candP.fMcPdg) == 211) {
1300 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "gpi", parRec.fMinv);
1301 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1302 fH.FillH1(hName + "pipi0", parRec.fMinv);
1303 else if (std::abs(candM.fMcPdg) == 211)
1304 fH.FillH1(hName + "pipi", parRec.fMinv);
1305 else
1306 fH.FillH1(hName + "pio", parRec.fMinv);
1307 }
1308
1309 // cand1 is neither electron nor misid. charged pion
1310 else {
1311 if (std::abs(candM.fMcPdg) == 11 && candM.IsMcGamma()) fH.FillH1(hName + "go", parRec.fMinv);
1312 else if (std::abs(candM.fMcPdg) == 11 && candM.IsMcPi0())
1313 fH.FillH1(hName + "pi0o", parRec.fMinv);
1314 else if (std::abs(candM.fMcPdg) == 211)
1315 fH.FillH1(hName + "pipi0", parRec.fMinv);
1316 else
1317 fH.FillH1(hName + "oo", parRec.fMinv);
1318 }
1319 }
1320 }
1321 }
1322}
1323
1324void LmvmTask::TrackSource(const LmvmCand& cand, ELmvmAnaStep step, int pdg)
1325{
1326 // no need to fill histograms for MC and Acc steps
1327 if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) return;
1328
1329 double stepBin = static_cast<double>(step) + 0.5;
1330
1331 FillMomHists(nullptr, &cand, cand.fMcSrc, step);
1332 fH.FillH1("hCandPdg", step, cand.fMcPdg);
1333
1334 int absPdg = std::abs(pdg);
1335 double pdgBin = (absPdg == 11 && cand.IsMcSignal()) ? 0.5
1336 : (absPdg == 11 && !cand.IsMcSignal()) ? 1.5
1337 : (absPdg == 211) ? 2.5
1338 : (absPdg == 2212) ? 3.5
1339 : (absPdg == 321) ? 4.5
1340 : 5.5;
1341 fH.FillH2("hCandPdgVsMom", step, cand.fMomentum.Mag(), pdgBin);
1342
1343 if (cand.IsMcSignal()) {
1344 fH.FillH1("hNofSignalTracks", stepBin, fW);
1345 fH.FillH2("hCandElSrc", stepBin, 7.5, fW);
1346 }
1347 else {
1348 if (LmvmUtils::IsMismatch(cand)) fH.FillH1("hNofMismatches_all", stepBin);
1349 if (cand.fStsMcTrackId != cand.fRichMcTrackId) fH.FillH1("hNofMismatches_rich", stepBin);
1350 if (cand.fStsMcTrackId != cand.fTrdMcTrackId) fH.FillH1("hNofMismatches_trd", stepBin);
1351 if (cand.fStsMcTrackId != cand.fTofMcTrackId) fH.FillH1("hNofMismatches_tof", stepBin);
1352 if (LmvmUtils::IsGhost(cand)) fH.FillH1("hNofGhosts", stepBin);
1353 fH.FillH1("hNofBgTracks", stepBin);
1354
1355 if (cand.IsMcGamma()) {
1356 CbmMCTrack* mctrack = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1357 if (mctrack != nullptr) {
1358 TVector3 v;
1359 mctrack->GetStartVertex(v);
1360 fH.FillH2("hVertexGammaXZ", step, v.Z(), v.X());
1361 fH.FillH2("hVertexGammaYZ", step, v.Z(), v.Y());
1362 fH.FillH2("hVertexGammaXY", step, v.X(), v.Y());
1363 fH.FillH2("hVertexGammaRZ", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
1364 }
1365 }
1366
1367 double srcBin = 0.0;
1368 if (cand.IsMcGamma()) srcBin = 0.5;
1369 else if (cand.IsMcPi0())
1370 srcBin = 1.5;
1371 else if (std::abs(pdg) == 211)
1372 srcBin = 2.5;
1373 else if (pdg == 2212)
1374 srcBin = 3.5;
1375 else if (std::abs(pdg) == 321)
1376 srcBin = 4.5;
1377 else if ((std::abs(pdg) == 11) && !cand.IsMcGamma() && !cand.IsMcPi0() && !cand.IsMcSignal())
1378 srcBin = 5.5;
1379 else
1380 srcBin = 6.5;
1381 fH.FillH2("hBgSrcTracks", stepBin, srcBin);
1382 if (std::abs(cand.fMcPdg) == 11) fH.FillH2("hCandElSrc", stepBin, srcBin);
1383 }
1384}
1385
1386void LmvmTask::BgPairPdg(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step)
1387{
1388 int pdgX = candP.fMcPdg;
1389 int pdgY = candM.fMcPdg;
1390
1391 double pdgBinX = (std::abs(pdgX) == 11 && candP.IsMcSignal()) ? 0.5
1392 : (std::abs(pdgX) == 11 && !candP.IsMcSignal()) ? 1.5
1393 : (std::abs(pdgX) == 211) ? 2.5
1394 : (pdgX == 2212) ? 3.5
1395 : (pdgX == 321) ? 4.5
1396 : (pdgX == 3112 or pdgX == 3222) ? 5.5
1397 : (std::abs(pdgX) == 13) ? 6.5
1398 : 7.5;
1399 double pdgBinY = (std::abs(pdgY) == 11 && candM.IsMcSignal()) ? 0.5
1400 : (std::abs(pdgY) == 11 && !candM.IsMcSignal()) ? 1.5
1401 : (std::abs(pdgY) == 211) ? 2.5
1402 : (pdgY == 2212) ? 3.5
1403 : (pdgY == 321) ? 4.5
1404 : (pdgY == 3112 or pdgY == 3222) ? 5.5
1405 : (std::abs(pdgY) == 13) ? 6.5
1406 : 7.5;
1407
1408 fH.FillH2("hBgPairPdg", step, pdgBinX, pdgBinY);
1409}
1410
1411void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const LmvmKinePar& parMc,
1412 const LmvmKinePar& parRec, ELmvmAnaStep step)
1413{
1414 // no need to fill histograms for MC and Acc steps
1415 if (step == ELmvmAnaStep::Mc || step == ELmvmAnaStep::Acc) return;
1416 bool isMismatch = (LmvmUtils::IsMismatch(candP) || LmvmUtils::IsMismatch(candM));
1417 ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM);
1418
1419 double w = 1.;
1420 if (candP.IsMcSignal()) w = fW * MinvScale(candP.fMassSig);
1421 if (candM.IsMcSignal()) w = fW * MinvScale(candM.fMassSig);
1422 if (w < 0) LOG(warning) << "LmvmTask::FillPairHists(): Signal mass < 0!";
1423
1424 fH.FillH1("hMinv", src, step, parRec.fMinv, w);
1425 if (!candP.IsMcSignal() && !candM.IsMcSignal()) fH.FillH1("hMinv_urqmdAll", src, step, parRec.fMinv, w);
1426 if (!candP.IsMcSignal() && !candM.IsMcSignal() && std::abs(candP.fMcPdg) == 11 && std::abs(candM.fMcPdg) == 11)
1427 fH.FillH1("hMinv_urqmdEl", src, step, parRec.fMinv, w);
1428
1429 fH.FillH2("hMinvPt", src, step, parRec.fMinv, parMc.fPt, w);
1430
1431 PairSource(candP, candM, step, parRec);
1432
1433 if (src == ELmvmSrc::Signal) {
1434 fH.FillH2("hPtYPairSignal", step, parMc.fRapidity, parMc.fPt, fW);
1435 fH.FillH1("hMomPairSignal", step, parMc.fMomentumMag, fW);
1436 }
1437 if (src == ELmvmSrc::Bg) {
1438 BgPairPdg(candP, candM, step);
1439 if (isMismatch) { fH.FillH1("hMinvBgMatch_mismatch", step, parRec.fMinv); }
1440 else {
1441 fH.FillH1("hMinvBgMatch_trueMatch", step, parRec.fMinv);
1442 if (std::abs(candP.fMcPdg) == 11 && std::abs(candM.fMcPdg) == 11)
1443 fH.FillH1("hMinvBgMatch_trueMatchEl", step, parRec.fMinv);
1444 if (std::abs(candP.fMcPdg) != 11 || std::abs(candM.fMcPdg) != 11)
1445 fH.FillH1("hMinvBgMatch_trueMatchNotEl", step, parRec.fMinv);
1446 }
1447 }
1448}
1449
1450string LmvmTask::GetPidString(const CbmMCTrack* mct, const LmvmCand* cand)
1451{
1452 if ((mct != nullptr && cand != nullptr) || (mct == nullptr && cand == nullptr)) {
1453 LOG(error) << "LmvmTask::GetPidString: Both mct and cand are [not nullptr] or [nullptr].";
1454 return 0;
1455 }
1456
1457 int pdg = (mct != nullptr) ? mct->GetPdgCode() : cand->fMcPdg;
1458 bool isMcSignal = (mct != nullptr) ? (mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11)
1459 : cand->IsMcSignal();
1460
1461 string pidString = fH.fCandNames[9];
1462 if (isMcSignal && pdg == -11) pidString = fH.fCandNames[0];
1463 else if (isMcSignal && pdg == 11)
1464 pidString = fH.fCandNames[1];
1465 else if (!isMcSignal && pdg == -11)
1466 pidString = fH.fCandNames[2];
1467 else if (!isMcSignal && pdg == 11)
1468 pidString = fH.fCandNames[3];
1469 else if (pdg == 211)
1470 pidString = fH.fCandNames[4];
1471 else if (pdg == -211)
1472 pidString = fH.fCandNames[5];
1473 else if (pdg == 2212)
1474 pidString = fH.fCandNames[6];
1475 else if (pdg == 321)
1476 pidString = fH.fCandNames[7];
1477 else if (pdg == -321)
1478 pidString = fH.fCandNames[8];
1479
1480 return pidString;
1481}
1482
1483string LmvmTask::GetPidString(double vertexMag, int pdg)
1484{
1485 string pidString = fH.fGTrackNames[14];
1486 bool isPrim = IsPrimary(vertexMag);
1487
1488 if (isPrim) {
1489 if (pdg == -11) pidString = fH.fGTrackNames[1];
1490 else if (pdg == 11)
1491 pidString = fH.fGTrackNames[3];
1492 else if (pdg == 211)
1493 pidString = fH.fGTrackNames[5];
1494 else if (pdg == -211)
1495 pidString = fH.fGTrackNames[7];
1496 else if (pdg == 2212)
1497 pidString = fH.fGTrackNames[9];
1498 else if (pdg == 321)
1499 pidString = fH.fGTrackNames[11];
1500 else if (pdg == -321)
1501 pidString = fH.fGTrackNames[13];
1502 }
1503 else {
1504 if (pdg == -11) pidString = fH.fGTrackNames[0];
1505 else if (pdg == 11)
1506 pidString = fH.fGTrackNames[2];
1507 else if (pdg == 211)
1508 pidString = fH.fGTrackNames[4];
1509 else if (pdg == -211)
1510 pidString = fH.fGTrackNames[6];
1511 else if (pdg == 2212)
1512 pidString = fH.fGTrackNames[8];
1513 else if (pdg == 321)
1514 pidString = fH.fGTrackNames[10];
1515 else if (pdg == -321)
1516 pidString = fH.fGTrackNames[12];
1517 }
1518
1519 return pidString;
1520}
1521
1522void LmvmTask::AnalyseCandidates()
1523{
1524 CheckGammaConvAndPi0();
1525 CheckTopologyCut(ELmvmTopologyCut::ST, "hStCut");
1526 CheckTopologyCut(ELmvmTopologyCut::TT, "hTtCut");
1527 CheckTopologyCut(ELmvmTopologyCut::RT, "hRtCut");
1528 if (fUseMvd) {
1529 CheckClosestMvdHit(1, "hMvdCut_1", "hMvdCutQa_1");
1530 CheckClosestMvdHit(2, "hMvdCut_2", "hMvdCutQa_2");
1531 }
1532
1533 // single candidates
1534 for (const auto& cand : fCands) {
1535 double w = (cand.IsMcSignal()) ? fW : 1.;
1536 CbmMCTrack* mcTrack = nullptr;
1537 if (cand.fStsMcTrackId >= 0) mcTrack = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1538 int pdg = mcTrack->GetPdgCode();
1539 double mom = mcTrack->GetP();
1540 string pidString = GetPidString(nullptr, &cand);
1541
1542 for (auto step : fH.fAnaSteps) {
1543 if (cand.IsCutTill(step)) {
1544 TrackSource(cand, step, pdg);
1545 FillCandPidValues(mcTrack, cand, step); // if std::abs(pdg) != 11, this is misidentification
1546 if (mcTrack != nullptr) CheckTofId(mcTrack, cand, step, pdg);
1547
1548 if (step >= ELmvmAnaStep::ElId && std::abs(pdg) == 211) fH.FillH1("hPionSupp_idEl", step, mom, w);
1549
1550 //double richDist = CbmRichUtil::GetRingTrackDistance(cand.fGTrackInd);
1551 //fH.FillH2("hRichRingTrackDist_cands_" + pidString, step, mom, richDist, w); // TODO: causes 'Error in <TClonesArray::At>: index xxx out of bounds' error
1552 }
1553 }
1554
1555 if (cand.IsCutTill(ELmvmAnaStep::ElId)) {
1556 fH.FillH2("hChi2VsMom_sts_" + pidString, cand.fMomentum.Mag(), cand.fChi2Sts, w);
1557 fH.FillH2("hChi2VsMom_rich_" + pidString, cand.fMomentum.Mag(), cand.fChi2Rich, w);
1558 fH.FillH2("hChi2VsMom_trd_" + pidString, cand.fMomentum.Mag(), cand.fChi2Trd, w);
1559
1560 fH.FillH2("hTofTimeVsChi2_sts_" + pidString, cand.fChi2Sts, cand.fTime, w);
1561 fH.FillH2("hTofTimeVsChi2_rich_" + pidString, cand.fChi2Rich, cand.fTime, w);
1562 fH.FillH2("hTofTimeVsChi2_trd_" + pidString, cand.fChi2Trd, cand.fTime, w);
1563
1564 fH.FillH2("hChi2Comb_StsRich_" + pidString, cand.fChi2Sts, cand.fChi2Rich, w);
1565 fH.FillH2("hChi2Comb_StsTrd_" + pidString, cand.fChi2Sts, cand.fChi2Trd, w);
1566 fH.FillH2("hChi2Comb_RichTrd_" + pidString, cand.fChi2Rich, cand.fChi2Trd, w);
1567
1568 fH.FillH2("hTofTimeVsMom_cands_" + pidString, cand.fMomentum.Mag(), cand.fTime, w);
1569 fH.FillH2("hTofHitTrackDist_cands_" + pidString, cand.fMomentum.Mag(), cand.fTofDist, w);
1570 }
1571
1572 DifferenceSignalAndBg(cand);
1573
1574 // beta-momentum spectrum
1575 string ptcl = GetPidString(nullptr, &cand);
1576 double p = cand.fMomentum.Mag();
1577 double m = cand.fMass;
1578 double r = p / m;
1579 double beta = r / sqrt(1. + r * r);
1580 double q = (cand.fCharge > 0) ? 1 : (cand.fCharge < 0) ? -1. : 0.;
1581 fH.FillH2("hBetaMom_cands_" + ptcl, q * p, beta);
1582 } // single candidates
1583
1584 // candidate pairs
1585 for (const auto& candP : fCands) {
1586 if (candP.fCharge < 0) continue;
1587 CbmMCTrack* mctrackP =
1588 (candP.fStsMcTrackId >= 0) ? static_cast<CbmMCTrack*>(fMCTracks->At(candP.fStsMcTrackId)) : nullptr;
1589 for (const auto& candM : fCands) {
1590 if (candM.fCharge > 0) continue;
1591
1592 CbmMCTrack* mctrackM =
1593 (candM.fStsMcTrackId >= 0) ? static_cast<CbmMCTrack*>(fMCTracks->At(candM.fStsMcTrackId)) : nullptr;
1594
1595 LmvmKinePar pMC = LmvmKinePar::Create(mctrackP, mctrackM);
1596 LmvmKinePar pRec = LmvmKinePar::Create(&candP, &candM);
1597
1598 for (auto step : fH.fAnaSteps) {
1599 if (candP.IsCutTill(step) && candM.IsCutTill(step)) FillPairHists(candP, candM, pMC, pRec, step);
1600 }
1601 }
1602 }
1603}
1604
1605void LmvmTask::CheckTofId(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step, int pdg)
1606{
1607 TVector3 v;
1608 mct->GetStartVertex(v);
1609 bool isPrim = IsPrimary(v.Mag());
1610 double pt = mct->GetPt();
1611 double rap = mct->GetRapidity();
1612
1613 // check PIDs in "Tof pile"
1614 if (cand.fMomentum.Mag() > 0.3 && cand.fMomentum.Mag() < 1. && cand.fMass2 > -0.012 && cand.fMass2 < 0.01) {
1615 string pidString = GetPidString(nullptr, &cand);
1616 double pdgBin = (pdg == 11 && cand.IsMcSignal()) ? 0.5
1617 : (pdg == -11 && cand.IsMcSignal()) ? 1.5
1618 : (pdg == 11 && !cand.IsMcSignal() && isPrim) ? 2.5
1619 : (pdg == -11 && !cand.IsMcSignal() && isPrim) ? 3.5
1620 : (pdg == 11 && !cand.IsMcSignal() && !isPrim) ? 4.5
1621 : (pdg == -11 && !cand.IsMcSignal() && !isPrim) ? 5.5
1622 : (pdg == 211) ? 6.5
1623 : (pdg == -211) ? 7.5
1624 : (pdg == 2212) ? 8.5
1625 : (pdg == 321) ? 9.5
1626 : (pdg == -321) ? 10.5
1627 : 11.5;
1628 fH.FillH1("hTofPilePdgs_cands", step, pdgBin);
1629 if (step == ELmvmAnaStep::ElId && cand.IsCutTill(step))
1630 fH.FillH2("hTofPilePty_cands_" + pidString, rap, pt); // TODO: check newly (25.10.22) added '&& IsTill(step)'!!
1631 }
1632
1633 // check vertex of misidentified particles in ToF after electron ID // TODO: split this up into single contributions?
1634 if (std::abs(pdg) != 11 && cand.fIsTofElectron) {
1635 fH.FillH2("hVertexXZ_misidTof", step, v.Z(), v.X());
1636 fH.FillH2("hVertexYZ_misidTof", step, v.Z(), v.Y());
1637 fH.FillH2("hVertexXY_misidTof", step, v.X(), v.Y());
1638 fH.FillH2("hVertexRZ_misidTof", step, v.Z(), sqrt(v.X() * v.X() + v.Y() * v.Y()));
1639 }
1640}
1641
1642void LmvmTask::FillCandPidValues(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step)
1643{
1644 double w = (cand.IsMcSignal()) ? fW : 1.;
1645 string pidString = GetPidString(nullptr, &cand);
1646
1647 fH.FillH1("hMom_cands_" + pidString, step, cand.fMomentum.Mag(), w);
1648 fH.FillH2("hPtY_cands_" + pidString, step, cand.fRapidity, cand.fMomentum.Perp(), w);
1649 fH.FillH2("hTofM2_cands_" + pidString, step, cand.fMomentum.Mag(), cand.fMass2, w);
1650
1651 if (mct != nullptr) {
1652 double rat = cand.fMomentum.Mag() / mct->GetP();
1653 string ptcl = GetPidString(nullptr, &cand);
1654 fH.FillH1("hMomRatio_cands_" + pidString, step, rat);
1655 fH.FillH2("hMomRatioVsMom_cands_" + ptcl, step, mct->GetP(), rat);
1656 }
1657}
1658
1659void LmvmTask::CheckGammaConvAndPi0()
1660{
1661 for (auto& candP : fCands) {
1662 if (candP.fCharge < 0) continue;
1663 for (auto& candM : fCands) {
1664 if (candM.fCharge > 0) continue;
1665 if (candP.IsCutTill(ELmvmAnaStep::ElId) && candM.IsCutTill(ELmvmAnaStep::ElId)) {
1666 LmvmKinePar pRec = LmvmKinePar::Create(&candP, &candM);
1667 if (!fCuts.IsGammaCutOk(pRec.fMinv)) {
1668 candM.fIsGammaCut = false;
1669 candP.fIsGammaCut = false;
1670 }
1671 }
1672 }
1673 }
1674}
1675
1676void LmvmTask::CheckTopologyCut(ELmvmTopologyCut cut, const string& name)
1677{
1678 string hcut = name + "_all";
1679 string hcutPion = name + "_pion";
1680 string hcutTruePair = name + "_truePair";
1681
1682 vector<LmvmDataAngMomInd> dataV;
1683
1684 vector<LmvmCand>& tpCands = fSTCands;
1685 if (cut == ELmvmTopologyCut::ST) { tpCands = fSTCands; }
1686 else if (cut == ELmvmTopologyCut::RT) {
1687 tpCands = fRTCands;
1688 }
1689 else if (cut == ELmvmTopologyCut::TT) {
1690 tpCands = fTTCands;
1691 }
1692 else {
1693 LOG(error) << "LmvmTask::CheckTopologyCut cut is not defined.";
1694 }
1695
1696 for (auto& cand : fCands) {
1697 if (cand.IsCutTill(ELmvmAnaStep::ElId)) {
1698 dataV.clear();
1699 for (size_t iM = 0; iM < tpCands.size(); iM++) {
1700 // different charges, charge iM != charge iP
1701 if (tpCands[iM].fCharge != cand.fCharge) {
1702 LmvmKinePar pRec = LmvmKinePar::Create(&cand, &tpCands[iM]);
1703 dataV.emplace_back(pRec.fAngle, tpCands[iM].fMomentum.Mag(), iM);
1704 }
1705 }
1706 //find min opening angle
1707 double minAng = 360.;
1708 int minInd = -1;
1709 for (size_t i = 0; i < dataV.size(); i++) {
1710 if (minAng > dataV[i].fAngle) {
1711 minAng = dataV[i].fAngle;
1712 minInd = i;
1713 }
1714 }
1715 if (minInd == -1) {
1716 cand.SetIsTopologyCutElectron(cut, true);
1717 continue;
1718 }
1719 bool isCut = fCuts.IsTopologyCutOk(cut, cand.fMomentum.Mag(), dataV[minInd].fMom, minAng);
1720 cand.SetIsTopologyCutElectron(cut, isCut);
1721
1722 // histogramms
1723 double sqrt_mom = TMath::Sqrt(cand.fMomentum.Mag() * dataV[minInd].fMom);
1724 int cutCandInd = dataV[minInd].fInd;
1725 int stsInd = tpCands[cutCandInd].fStsInd;
1726 if (stsInd < 0) continue;
1727 int pdgAbs = std::abs(tpCands[cutCandInd].fMcPdg);
1728 int motherId = tpCands[cutCandInd].fMcMotherId;
1729
1730 fH.FillH2(hcut, cand.fMcSrc, sqrt_mom, minAng, fW);
1731 if (pdgAbs == 211) fH.FillH2(hcutPion, cand.fMcSrc, sqrt_mom, minAng, fW);
1732 if (cand.IsMcSignal()) {
1733 if (motherId == cand.fMcMotherId) fH.FillH2(hcutTruePair, cand.fMcSrc, sqrt_mom, minAng, fW);
1734 }
1735 else {
1736 if (motherId != -1 && motherId == cand.fMcMotherId) fH.FillH2(hcutTruePair, cand.fMcSrc, sqrt_mom, minAng, fW);
1737 }
1738 }
1739 }
1740}
1741
1742void LmvmTask::CalculateNofTopologyPairs(const string& name, ELmvmSrc src)
1743{
1744 size_t nCand = fCands.size();
1745 for (size_t iP = 0; iP < nCand; iP++) {
1746 const LmvmCand& cand = fCands[iP];
1747 if (cand.fMcMotherId == -1) continue;
1748 if (src != cand.fMcSrc) continue;
1749 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue;
1750
1751 bool isAdded = false;
1752
1753 // 3 topology cuts: ST, RT, TT
1754 for (int i = 0; i < 3; i++) {
1755 if (isAdded) continue;
1756 vector<LmvmCand>& cands = fSTCands;
1757 double binNum = 4.5;
1758 if (i == 1) {
1759 cands = fRTCands;
1760 binNum = 5.5;
1761 }
1762 else if (i == 2) {
1763 cands = fTTCands;
1764 binNum = 6.5;
1765 }
1766 for (const auto& candT : cands) {
1767 if (candT.fMcMotherId == cand.fMcMotherId) {
1768 fH.FillH1(name, binNum);
1769 isAdded = true;
1770 break;
1771 }
1772 }
1773 }
1774 if (isAdded) continue;
1775
1776 for (size_t iM = 0; iM < fCands.size(); iM++) {
1777 if (iM != iP && fCands[iM].fMcMotherId == cand.fMcMotherId && fCands[iM].IsCutTill(ELmvmAnaStep::ElId)) {
1778 fH.FillH1(name, 7.5);
1779 isAdded = true;
1780 break;
1781 }
1782 }
1783 if (isAdded) continue;
1784
1785 int nofStsPoints = 0;
1786 int nofMcTracks = fMCTracks->GetEntriesFast();
1787 for (int iMc = 0; iMc < nofMcTracks; iMc++) {
1788 const CbmMCTrack* mcTrack = static_cast<const CbmMCTrack*>(fMCTracks->At(iMc));
1789 if (mcTrack == nullptr || mcTrack->GetMotherId() != cand.fMcMotherId || iMc == cand.fStsMcTrackId) continue;
1790
1791 int eventId = FairRun::Instance()->GetEventHeader()->GetMCEntryNumber();
1792 if (!CbmLitMCTrackCreator::Instance()->TrackExists(eventId, iMc)) continue;
1793 const CbmLitMCTrack& litMCTrack = CbmLitMCTrackCreator::Instance()->GetTrack(eventId, iMc);
1794 nofStsPoints = litMCTrack.GetNofPointsInDifferentStations(ECbmModuleId::kSts);
1795 break;
1796 }
1797 if (nofStsPoints == 0) fH.FillH1(name, 0.5);
1798 if (nofStsPoints >= 1 && nofStsPoints <= 3) fH.FillH1(name, 1.5);
1799 if (nofStsPoints >= 4 && nofStsPoints <= 5) fH.FillH1(name, 2.5);
1800 if (nofStsPoints >= 6) fH.FillH1(name, 3.5);
1801 }
1802}
1803
1804void LmvmTask::DifferenceSignalAndBg(const LmvmCand& cand)
1805{
1806 fH.FillH1("hChi2PrimVertex", cand.fMcSrc, cand.fChi2Prim, fW);
1807
1808 if (!cand.fIsChi2Prim) return;
1809 fH.FillH1("hAnnRich", cand.fMcSrc, cand.fRichAnn, fW);
1810 fH.FillH2("hAnnRichVsMom", cand.fMcSrc, cand.fMomentum.Mag(), cand.fRichAnn, fW);
1811 //fH.FillH1("hAnnTrd", cand.fMcSrc, cand.fTrdAnn, fW); // TODO: uncomment when TRD ANN is working (CbmLitGlobalElectronId::GetTrdAnn() gives back El-Likelihood)
1812 fH.FillH2("hTrdLike_El", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikeEl, fW);
1813 fH.FillH2("hTrdLike_Pi", cand.fMcSrc, cand.fMomentum.Mag(), cand.fTrdLikePi, fW);
1814 fH.FillH2("hTofM2", cand.fMcSrc, cand.fMomentum.Mag(), cand.fMass2, fW);
1815
1816 // electron purity
1817 if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) == 11) {
1818 fH.FillH2("hAnnRichVsMomPur_El", cand.fMomentum.Mag(), cand.fRichAnn, fW);
1819 fH.FillH2("hTrdElLikePur_El", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW);
1820 }
1821 else if (!cand.IsMcSignal() && std::abs(cand.fMcPdg) != 11) {
1822 fH.FillH2("hAnnRichVsMomPur_Bg", cand.fMomentum.Mag(), cand.fRichAnn, fW);
1823 fH.FillH2("hTrdElLikePur_Bg", cand.fMomentum.Mag(), cand.fTrdLikeEl, fW);
1824 }
1825
1826 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) return;
1827 //fH.FillSourceH1("hPt", cand.fMcSrc, cand.fMomentum.Perp(), fW);
1828 //fH.FillSourceH1("hMom", cand.fMcSrc, cand.fMomentum.Mag(), fW);
1829 fH.FillH1("hChi2Sts", cand.fMcSrc, cand.fChi2Sts, fW);
1830
1831 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
1832 if (stsTrack == nullptr) return;
1833 fH.FillH1("hNofStsHits", cand.fMcSrc, stsTrack->GetNofStsHits(), fW);
1834
1835 if (fUseMvd) {
1836 double mvd1x = 0., mvd1y = 0., mvd2x = 0., mvd2y = 0.;
1837 for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
1838 CbmMvdHit* mvdHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
1839 if (mvdHit == nullptr) return;
1840 if (mvdHit->GetStationNr() == 1) {
1841 mvd1x = mvdHit->GetX();
1842 mvd1y = mvdHit->GetY();
1843 }
1844 else if (mvdHit->GetStationNr() == 2) {
1845 mvd2x = mvdHit->GetX();
1846 mvd2y = mvdHit->GetY();
1847 }
1848 }
1849 double mvd1r = sqrt(mvd1x * mvd1x + mvd1y * mvd1y);
1850 double mvd2r = sqrt(mvd2x * mvd2x + mvd2y * mvd2y);
1851 fH.FillH1("hNofMvdHits", cand.fMcSrc, stsTrack->GetNofMvdHits(), fW);
1852 fH.FillH2("hMvdXY_1", cand.fMcSrc, mvd1x, mvd1y, fW);
1853 fH.FillH1("hMvdR_1", cand.fMcSrc, mvd1r, fW);
1854 fH.FillH2("hMvdXY_2", cand.fMcSrc, mvd2x, mvd2y, fW);
1855 fH.FillH1("hMvdR_2", cand.fMcSrc, mvd2r, fW);
1856 }
1857}
1858
1859void LmvmTask::CheckClosestMvdHit(int mvdStationNum, const string& hist, const string& histQa)
1860{
1861 vector<LmvmDataXYInd> mvdV;
1862 vector<LmvmDataXYInd> candV;
1863
1864 for (int iHit = 0; iHit < fMvdHits->GetEntriesFast(); iHit++) {
1865 CbmMvdHit* mvdHit = static_cast<CbmMvdHit*>(fMvdHits->At(iHit));
1866 if (mvdHit != nullptr && mvdHit->GetStationNr() == mvdStationNum) {
1867 mvdV.emplace_back(mvdHit->GetX(), mvdHit->GetY(), iHit);
1868 }
1869 }
1870
1871 for (size_t iC = 0; iC < fCands.size(); iC++) {
1872 if (fCands[iC].IsCutTill(ELmvmAnaStep::ElId)) {
1873 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(fCands[iC].fStsInd));
1874 if (stsTrack == nullptr) continue;
1875 for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
1876 CbmMvdHit* candHit = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
1877 if (candHit != nullptr && candHit->GetStationNr() == mvdStationNum) {
1878 candV.emplace_back(candHit->GetX(), candHit->GetY(), iC);
1879 }
1880 }
1881 }
1882 }
1883
1884 for (size_t iC = 0; iC < candV.size(); iC++) {
1885 LmvmCand& cand = fCands[candV[iC].fInd];
1886 double minD = 9999999.;
1887 int minMvdInd = -1;
1888 for (size_t iH = 0; iH < mvdV.size(); iH++) {
1889 double d2 = LmvmUtils::Distance2(mvdV[iH].fX, mvdV[iH].fY, candV[iC].fX, candV[iC].fY);
1890 if (d2 < 1.e-9) continue;
1891 if (d2 < minD) {
1892 minMvdInd = mvdV[iH].fInd;
1893 minD = d2;
1894 }
1895 }
1896 double dmvd = sqrt(minD);
1897
1898 // Check MVD cut quality
1899 double bin = -1.;
1900 const CbmMatch* hitMatch = static_cast<const CbmMatch*>(fMvdHitMatches->At(minMvdInd));
1901 if (hitMatch != nullptr) {
1902 CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(fMCTracks->At(hitMatch->GetMatchedLink().GetIndex()));
1903 int mcMvdHitPdg = TMath::Abs(mct1->GetPdgCode());
1904 int mvdMotherId = mct1->GetMotherId();
1905
1906 int stsMotherId = -2;
1907 if (cand.fStsMcTrackId >= 0) {
1908 CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId));
1909 stsMotherId = (mct2 != nullptr) ? mct2->GetMotherId() : -2;
1910 }
1911
1912 bin = (mvdMotherId != -1 && mvdMotherId == stsMotherId) ? 0.5 : 1.5; // correct or wrong assignment
1913 if (cand.IsMcSignal()) {
1914 bin = (mvdMotherId == stsMotherId && mcMvdHitPdg == 11) ? 0.5 : 1.5; // correct or wrong assignment
1915 }
1916 }
1917
1918 // Fill histograms
1919 fH.FillH1(histQa, cand.fMcSrc, bin, fW);
1920 fH.FillH2(hist, cand.fMcSrc, dmvd, cand.fMomentum.Mag(), fW);
1921
1922 // Apply MVD cut
1923 bool isMvdCut = fCuts.IsMvdCutOk(mvdStationNum, dmvd, cand.fMomentum.Mag());
1924 if (mvdStationNum == 1) cand.fIsMvd1Cut = isMvdCut;
1925 else if (mvdStationNum == 2)
1926 cand.fIsMvd2Cut = isMvdCut;
1927 }
1928}
1929
1930void LmvmTask::MvdCutMcDistance()
1931{
1932 if (!fUseMvd) return;
1933 for (const auto& cand : fCands) {
1934 if (!cand.IsCutTill(ELmvmAnaStep::ElId)) continue;
1935 CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(cand.fStsInd));
1936 if (stsTrack == nullptr) continue;
1937 for (int iM = 0; iM < stsTrack->GetNofMvdHits(); iM++) {
1938 CbmMvdHit* mvdHit1 = static_cast<CbmMvdHit*>(fMvdHits->At(stsTrack->GetMvdHitIndex(iM)));
1939 if (mvdHit1 == nullptr) continue;
1940
1941 int nofMvdHits = fMvdHitMatches->GetEntriesFast();
1942 for (int iMvd = 0; iMvd < nofMvdHits; iMvd++) {
1943 const CbmMatch* hitMatch = static_cast<const CbmMatch*>(fMvdHitMatches->At(iMvd));
1944 if (hitMatch == nullptr) continue;
1945 if (cand.fStsMcTrackId != hitMatch->GetMatchedLink().GetIndex()) continue;
1946 CbmMvdHit* mvdHit2 = static_cast<CbmMvdHit*>(fMvdHits->At(iMvd));
1947 if (mvdHit2 == nullptr || mvdHit2->GetStationNr() != mvdHit1->GetStationNr()) continue;
1948 double d = LmvmUtils::Distance(mvdHit1->GetX(), mvdHit1->GetY(), mvdHit2->GetX(), mvdHit2->GetY());
1949 if (mvdHit1->GetStationNr() == 1) { fH.FillH1("hMvdMcDist_1", cand.fMcSrc, d, fW); }
1950 else if (mvdHit1->GetStationNr() == 2) {
1951 fH.FillH1("hMvdMcDist_2", cand.fMcSrc, d, fW);
1952 }
1953 }
1954 }
1955 }
1956}
1957
1958bool LmvmTask::IsPrimary(double vertexMag)
1959{
1960 if (vertexMag < fZ + 0.1 && vertexMag > fZ - 0.1) return true;
1961 return false;
1962}
1963
1964void LmvmTask::Finish()
1965{
1966 CombinatorialPairs();
1967 TDirectory* oldir = gDirectory;
1968 TFile* outFile = FairRootManager::Instance()->GetOutFile();
1969 if (outFile != nullptr) {
1970 outFile->cd();
1971 fH.WriteToFile();
1972 }
1973 gDirectory->cd(oldir->GetPath());
1974}
1975
1976void LmvmTask::SetEnergyAndPlutoParticle(const string& energy, const string& particle)
1977{
1978 this->SetWeight(LmvmSimParam::GetWeight(energy, particle));
1979 fParticle = particle;
1980}
ECbmModuleId::kMvd
@ kMvd
Micro-Vertex Detector.
ECbmModuleId::kTrd
@ kTrd
Transition Radiation Detector.
ECbmModuleId::kTof
@ kTof
Time-of-flight Detector.
ECbmModuleId::kSts
@ kSts
Silicon Tracking System.
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
ELmvmBgPairSrc
ELmvmBgPairSrc
Definition LmvmDef.h:51
ELmvmAnaStep
ELmvmAnaStep
Definition LmvmDef.h:34
ELmvmTopologyCut
ELmvmTopologyCut
Definition LmvmDef.h:15
ClassImp
ClassImp(LmvmTask)
CbmGlobalTrack::GetStsTrackIndex
int32_t GetStsTrackIndex() const
Definition CbmGlobalTrack.h:44
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
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:76
CbmHit::GetRefId
int32_t GetRefId() const
Definition CbmHit.h:73
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition CbmL1PFFitter.cxx:501
CbmLitGlobalElectronId::IsTrdElectron
Bool_t IsTrdElectron(Int_t globalTrackindex, Double_t momentum)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:93
CbmLitGlobalElectronId::GetTofM2
Double_t GetTofM2(Int_t globalTrackIndex, Double_t momentum, Double_t eventTime=1000.)
Return TOF m2 value.
Definition CbmLitGlobalElectronId.cxx:161
CbmLitGlobalElectronId::IsTofElectron
Bool_t IsTofElectron(Int_t globalTrackIndex, Double_t momentum, Double_t eventTime=1000.)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:113
CbmLitGlobalElectronId::IsRichElectron
Bool_t IsRichElectron(Int_t globalTrackIndex, Double_t momentum)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:64
CbmLitGlobalElectronId::GetInstance
static CbmLitGlobalElectronId & GetInstance()
Definition CbmLitGlobalElectronId.h:30
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:80
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:72
CbmMCTrack::GetCharge
double GetCharge() const
Charge of the associated particle.
Definition CbmMCTrack.cxx:150
CbmMCTrack::GetPt
double GetPt() const
Definition CbmMCTrack.h:97
CbmMCTrack::GetP
double GetP() const
Definition CbmMCTrack.h:98
CbmMCTrack::GetPx
double GetPx() const
Definition CbmMCTrack.h:70
CbmMCTrack::GetGeantProcessId
uint32_t GetGeantProcessId() const
Definition CbmMCTrack.h:67
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetPy
double GetPy() const
Definition CbmMCTrack.h:71
CbmMCTrack::GetRapidity
double GetRapidity() const
Definition CbmMCTrack.cxx:187
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
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
CbmMvdHit::GetStationNr
virtual int32_t GetStationNr() const
Definition CbmMvdHit.h:61
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
CbmRichRing::GetChi2
double GetChi2() const
Definition CbmRichRing.h:92
CbmRichRing::GetNDF
double GetNDF() const
Definition CbmRichRing.h:93
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::GetNofStsHits
int32_t GetNofStsHits() const
Definition CbmStsTrack.h:93
CbmTofTrack::GetDistance
double GetDistance() const
Definition CbmTofTrack.h:79
CbmTrack::GetNDF
int32_t GetNDF() const
Definition CbmTrack.h:64
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:58
LmvmCand::IsMcPi0
bool IsMcPi0() const
Definition LmvmCand.h:103
LmvmCand::fMcPdg
int fMcPdg
Definition LmvmCand.h:81
LmvmCand::fRichInd
int fRichInd
Definition LmvmCand.h:77
LmvmCand::fTrdLikeEl
double fTrdLikeEl
Definition LmvmCand.h:86
LmvmCand::fChi2Sts
double fChi2Sts
Definition LmvmCand.h:57
LmvmCand::fGTrackInd
int fGTrackInd
Definition LmvmCand.h:82
LmvmCand::fStsMcTrackId
int fStsMcTrackId
Definition LmvmCand.h:72
LmvmCand::fTofMcTrackId
int fTofMcTrackId
Definition LmvmCand.h:75
LmvmCand::fEventNumber
int fEventNumber
Definition LmvmCand.h:71
LmvmCand::fMassSig
double fMassSig
Definition LmvmCand.h:52
LmvmCand::IsMcSignal
bool IsMcSignal() const
Definition LmvmCand.h:102
LmvmCand::fTrdLikePi
double fTrdLikePi
Definition LmvmCand.h:87
LmvmCand::fIsTofElectron
bool fIsTofElectron
Definition LmvmCand.h:68
LmvmCand::fTofTrackInd
int fTofTrackInd
Definition LmvmCand.h:80
LmvmCand::fRapidity
double fRapidity
Definition LmvmCand.h:54
LmvmCand::fChi2Trd
double fChi2Trd
Definition LmvmCand.h:59
LmvmCand::fIsElectron
bool fIsElectron
Definition LmvmCand.h:91
LmvmCand::fIsChi2Prim
bool fIsChi2Prim
Definition LmvmCand.h:90
LmvmCand::fTofDist
double fTofDist
Definition LmvmCand.h:63
LmvmCand::fMass2
double fMass2
Definition LmvmCand.h:85
LmvmCand::fRichAnn
double fRichAnn
Definition LmvmCand.h:83
LmvmCand::fTrdMcTrackId
int fTrdMcTrackId
Definition LmvmCand.h:74
LmvmCand::fTofInd
int fTofInd
Definition LmvmCand.h:79
LmvmCand::fIsMvd2Cut
bool fIsMvd2Cut
Definition LmvmCand.h:94
LmvmCand::fStsInd
int fStsInd
Definition LmvmCand.h:76
LmvmCand::fMcMotherId
int fMcMotherId
Definition LmvmCand.h:70
LmvmCand::fRichMcTrackId
int fRichMcTrackId
Definition LmvmCand.h:73
LmvmCand::fTime
double fTime
Definition LmvmCand.h:62
LmvmCand::IsMcGamma
bool IsMcGamma() const
Definition LmvmCand.h:104
LmvmCand::fMcSrc
ELmvmSrc fMcSrc
Definition LmvmCand.h:101
LmvmCand::fTrdInd
int fTrdInd
Definition LmvmCand.h:78
LmvmCand::fIsMvd1Cut
bool fIsMvd1Cut
Definition LmvmCand.h:93
LmvmCand::fChi2Prim
double fChi2Prim
Definition LmvmCand.h:56
LmvmCand::fMomentum
TVector3 fMomentum
Definition LmvmCand.h:50
LmvmCand::fIsPtCut
bool fIsPtCut
Definition LmvmCand.h:98
LmvmCand::fLength
double fLength
Definition LmvmCand.h:61
LmvmCand::IsCutTill
bool IsCutTill(ELmvmAnaStep step) const
Definition LmvmCand.h:33
LmvmCuts::IsChi2PrimaryOk
bool IsChi2PrimaryOk(double chi2Prim)
Definition LmvmCuts.h:43
LmvmCuts::IsMvdCutOk
bool IsMvdCutOk(int stationNum, double dmvd, double mom)
Definition LmvmCuts.h:49
LmvmCuts::IsTopologyCutOk
bool IsTopologyCutOk(ELmvmTopologyCut cut, double mom1, double mom2, double minAngle)
Definition LmvmCuts.h:19
LmvmCuts::IsGammaCutOk
bool IsGammaCutOk(double minv)
Definition LmvmCuts.h:45
LmvmCuts::ToString
std::string ToString()
Definition LmvmCuts.h:70
LmvmCuts::IsPtCutOk
bool IsPtCutOk(double pt)
Definition LmvmCuts.h:47
LmvmCuts::fMomentumCut
double fMomentumCut
Definition LmvmCuts.h:89
LmvmHist::CreateH1
void CreateH1(const std::string &name, const std::string &axisX, const std::string &axisY, double nBins, double min, double max)
LmvmHist::fCandNames
static const std::vector< std::string > fCandNames
Definition LmvmHist.h:57
LmvmHist::FillH1
void FillH1(const std::string &name, double x, double w=1.)
LmvmHist::fBgPairSrcNames
static const std::vector< std::string > fBgPairSrcNames
Definition LmvmHist.h:45
LmvmHist::CreateH2
void CreateH2(const std::string &name, const std::string &axisX, const std::string &axisY, const std::string &axisZ, double nBinsX, double minX, double maxX, double nBinsY, double minY, double maxY)
LmvmHist::fNofBgPairSrc
static const int fNofBgPairSrc
Definition LmvmHist.h:41
LmvmHist::fSrcNames
static const std::vector< std::string > fSrcNames
Definition LmvmHist.h:25
LmvmHist::GetName
std::string GetName(const std::string &name, ELmvmAnaStep step)
LmvmHist::FillH2
void FillH2(const std::string &name, double x, double y, double w=1.)
LmvmHist::fGTrackNames
static const std::vector< std::string > fGTrackNames
Definition LmvmHist.h:50
LmvmHist::fAnaSteps
static const std::vector< ELmvmAnaStep > fAnaSteps
Definition LmvmHist.h:29
LmvmHist::fNofAnaSteps
static const int fNofAnaSteps
Definition LmvmHist.h:31
LmvmHist::fAnaStepNames
static const std::vector< std::string > fAnaStepNames
Definition LmvmHist.h:33
LmvmHist::WriteToFile
void WriteToFile()
Definition LmvmHist.cxx:477
LmvmKinePar::fAngle
Double_t fAngle
Definition LmvmKinePar.h:23
LmvmKinePar::fMinv
Double_t fMinv
Definition LmvmKinePar.h:22
LmvmKinePar::fPt
Double_t fPt
Definition LmvmKinePar.h:20
LmvmKinePar::fMomentumMag
Double_t fMomentumMag
Definition LmvmKinePar.h:19
LmvmKinePar::fRapidity
Double_t fRapidity
Definition LmvmKinePar.h:21
LmvmKinePar::Create
static LmvmKinePar Create(const CbmMCTrack *mctrackP, const CbmMCTrack *mctrackM)
Definition LmvmKinePar.h:28
LmvmSimParam::GetWeight
static Double_t GetWeight(const std::string &energy, const std::string &particle)
Definition LmvmSimParam.h:16
LmvmTask::fRichHits
TClonesArray * fRichHits
Definition LmvmTask.h:178
LmvmTask::CombinatorialPairs
void CombinatorialPairs()
Definition LmvmTask.cxx:1076
LmvmTask::CheckGammaConvAndPi0
void CheckGammaConvAndPi0()
Definition LmvmTask.cxx:1659
LmvmTask::AnalyseCandidates
void AnalyseCandidates()
Definition LmvmTask.cxx:1522
LmvmTask::~LmvmTask
virtual ~LmvmTask()
Definition LmvmTask.cxx:58
LmvmTask::fNofHitsInRingMap
std::map< int, int > fNofHitsInRingMap
Definition LmvmTask.h:216
LmvmTask::InitHists
void InitHists()
Definition LmvmTask.cxx:61
LmvmTask::RichPmtXY
void RichPmtXY()
Definition LmvmTask.cxx:490
LmvmTask::fParticle
std::string fParticle
Definition LmvmTask.h:220
LmvmTask::fRichRingMatches
TClonesArray * fRichRingMatches
Definition LmvmTask.h:177
LmvmTask::fMCEventHeader
FairMCEventHeader * fMCEventHeader
Definition LmvmTask.h:172
LmvmTask::fSTCands
std::vector< LmvmCand > fSTCands
Definition LmvmTask.h:201
LmvmTask::fRichProj
TClonesArray * fRichProj
Definition LmvmTask.h:175
LmvmTask::AssignMcToTopologyCands
void AssignMcToTopologyCands(std::vector< LmvmCand > &topoCands)
Definition LmvmTask.cxx:1223
LmvmTask::fMCTracks
TClonesArray * fMCTracks
Definition LmvmTask.h:173
LmvmTask::fMvdHitMatches
TClonesArray * fMvdHitMatches
Definition LmvmTask.h:185
LmvmTask::fRichRings
TClonesArray * fRichRings
Definition LmvmTask.h:174
LmvmTask::DoMcTrack
void DoMcTrack()
Definition LmvmTask.cxx:396
LmvmTask::fPionMisidLevel
double fPionMisidLevel
Definition LmvmTask.h:209
LmvmTask::fCuts
LmvmCuts fCuts
Definition LmvmTask.h:212
LmvmTask::fTofPoints
TClonesArray * fTofPoints
Definition LmvmTask.h:191
LmvmTask::IsPrimary
bool IsPrimary(double vertexMag)
Definition LmvmTask.cxx:1958
LmvmTask::MinvScale
double MinvScale(double minv)
Definition LmvmTask.cxx:363
LmvmTask::fTofHitsMatches
TClonesArray * fTofHitsMatches
Definition LmvmTask.h:190
LmvmTask::SetWeight
void SetWeight(double w)
Definition LmvmTask.h:224
LmvmTask::fCandsTotal
std::vector< LmvmCand > fCandsTotal
Definition LmvmTask.h:199
LmvmTask::fMvdPoints
TClonesArray * fMvdPoints
Definition LmvmTask.h:184
LmvmTask::FillPairHists
void FillPairHists(const LmvmCand &candP, const LmvmCand &candM, const LmvmKinePar &parMc, const LmvmKinePar &parRec, ELmvmAnaStep step)
Definition LmvmTask.cxx:1411
LmvmTask::BgPairPdg
void BgPairPdg(const LmvmCand &candP, const LmvmCand &candM, ELmvmAnaStep step)
Definition LmvmTask.cxx:1386
LmvmTask::CalculateNofTopologyPairs
void CalculateNofTopologyPairs(const std::string &name, ELmvmSrc src)
Definition LmvmTask.cxx:1742
LmvmTask::DoMcPair
void DoMcPair()
Definition LmvmTask.cxx:450
LmvmTask::fKFVertex
CbmKFVertex fKFVertex
Definition LmvmTask.h:194
LmvmTask::SetEnergyAndPlutoParticle
void SetEnergyAndPlutoParticle(const std::string &energy, const std::string &particle)
Definition LmvmTask.cxx:1976
LmvmTask::fRichPoints
TClonesArray * fRichPoints
Definition LmvmTask.h:176
LmvmTask::fTTCands
std::vector< LmvmCand > fTTCands
Definition LmvmTask.h:203
LmvmTask::fCands
std::vector< LmvmCand > fCands
Definition LmvmTask.h:198
LmvmTask::fEventNumber
int fEventNumber
Definition LmvmTask.h:211
LmvmTask::FillCands
void FillCands()
Definition LmvmTask.cxx:977
LmvmTask::MvdCutMcDistance
void MvdCutMcDistance()
Definition LmvmTask.cxx:1930
LmvmTask::TrackSource
void TrackSource(const LmvmCand &cand, ELmvmAnaStep step, int pdg)
Definition LmvmTask.cxx:1324
LmvmTask::fMvdHits
TClonesArray * fMvdHits
Definition LmvmTask.h:183
LmvmTask::fUseMvd
bool fUseMvd
Definition LmvmTask.h:196
LmvmTask::IsInTofPile
bool IsInTofPile(double mom, double m2)
Definition LmvmTask.cxx:653
LmvmTask::AssignMcToCands
void AssignMcToCands(std::vector< LmvmCand > &cands)
Definition LmvmTask.cxx:1169
LmvmTask::AnalyseGlobalTracks
void AnalyseGlobalTracks()
Definition LmvmTask.cxx:520
LmvmTask::fStsHits
TClonesArray * fStsHits
Definition LmvmTask.h:182
LmvmTask::fH
LmvmHist fH
Definition LmvmTask.h:214
LmvmTask::CheckClosestMvdHit
void CheckClosestMvdHit(int mvdStationNum, const std::string &hist, const std::string &histQa)
Definition LmvmTask.cxx:1859
LmvmTask::fTofTracks
TClonesArray * fTofTracks
Definition LmvmTask.h:192
LmvmTask::IsMcTrackAccepted
bool IsMcTrackAccepted(int mcTrackInd)
Definition LmvmTask.cxx:511
LmvmTask::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition LmvmTask.h:181
LmvmTask::Init
virtual InitStatus Init()
Definition LmvmTask.cxx:270
LmvmTask::Exec
virtual void Exec(Option_t *option)
Definition LmvmTask.cxx:308
LmvmTask::fTofHits
TClonesArray * fTofHits
Definition LmvmTask.h:189
LmvmTask::BetaMom
void BetaMom(const CbmMCTrack *mct, const CbmGlobalTrack *gTrack, const std::string &ptcl)
Definition LmvmTask.cxx:754
LmvmTask::fTrdTracks
TClonesArray * fTrdTracks
Definition LmvmTask.h:186
LmvmTask::GetPidString
std::string GetPidString(const CbmMCTrack *mct, const LmvmCand *cand)
Definition LmvmTask.cxx:1450
LmvmTask::CheckTofId
void CheckTofId(const CbmMCTrack *mcTrack, const LmvmCand &cand, ELmvmAnaStep step, int pdg)
Definition LmvmTask.cxx:1605
LmvmTask::InitOrFatal
T * InitOrFatal(const std::string &name)
Definition LmvmTask.h:57
LmvmTask::DifferenceSignalAndBg
void DifferenceSignalAndBg(const LmvmCand &cand)
Definition LmvmTask.cxx:1804
LmvmTask::FillCandPidValues
void FillCandPidValues(const CbmMCTrack *mcTrack, const LmvmCand &cand, ELmvmAnaStep step)
Definition LmvmTask.cxx:1642
LmvmTask::fTrdTrackMatches
TClonesArray * fTrdTrackMatches
Definition LmvmTask.h:188
LmvmTask::CheckTopologyCut
void CheckTopologyCut(ELmvmTopologyCut cut, const std::string &name)
Definition LmvmTask.cxx:1676
LmvmTask::fGlobalTracks
TClonesArray * fGlobalTracks
Definition LmvmTask.h:179
LmvmTask::FillTopologyCands
void FillTopologyCands()
Definition LmvmTask.cxx:910
LmvmTask::CheckTofIdentification
void CheckTofIdentification(const CbmGlobalTrack *gTrack, const std::string &pidString, double mom, double m2, int pdg, bool isTofEl)
Definition LmvmTask.cxx:677
LmvmTask::PidVsMom
void PidVsMom(const CbmGlobalTrack *gTrack, int iGTrack, int pdg, double mom)
Definition LmvmTask.cxx:655
LmvmTask::fW
double fW
Definition LmvmTask.h:207
LmvmTask::fZ
double fZ
Definition LmvmTask.h:218
LmvmTask::fRTCands
std::vector< LmvmCand > fRTCands
Definition LmvmTask.h:205
LmvmTask::fPrimVertex
CbmVertex * fPrimVertex
Definition LmvmTask.h:193
LmvmTask::FillMomHists
void FillMomHists(const CbmMCTrack *mct, const LmvmCand *cand, ELmvmSrc src, ELmvmAnaStep step)
Definition LmvmTask.cxx:375
LmvmTask::Finish
virtual void Finish()
Definition LmvmTask.cxx:1964
LmvmTask::FillRichRingNofHits
void FillRichRingNofHits()
Definition LmvmTask.cxx:348
LmvmTask::PairSource
void PairSource(const LmvmCand &candP, const LmvmCand &candM, ELmvmAnaStep step, const LmvmKinePar &parRec)
Definition LmvmTask.cxx:1241
LmvmTask::IsInAllDets
bool IsInAllDets(const CbmGlobalTrack *gTrack)
Definition LmvmTask.cxx:888
LmvmTask::CheckMismatches
void CheckMismatches(const CbmGlobalTrack *gTrack, int pdg, bool isElectron, const std::string &ptcl, double weight)
Definition LmvmTask.cxx:772
LmvmTask::fStsTracks
TClonesArray * fStsTracks
Definition LmvmTask.h:180
LmvmUtils::IsTofElectron
static void IsTofElectron(int globalTrackIndex, double momentum, LmvmCand *cand)
Definition LmvmUtils.cxx:278
LmvmUtils::IsMcGammaEl
static bool IsMcGammaEl(const CbmMCTrack *mct, TClonesArray *mcTracks)
Definition LmvmUtils.cxx:97
LmvmUtils::Distance
static double Distance(double x1, double y1, double x2, double y2)
Definition LmvmUtils.cxx:238
LmvmUtils::GetMcPairSrc
static ELmvmSrc GetMcPairSrc(const CbmMCTrack *mctP, const CbmMCTrack *mctM, TClonesArray *mcTracks)
Definition LmvmUtils.cxx:150
LmvmUtils::IsElectronMc
static void IsElectronMc(LmvmCand *cand, TClonesArray *mcTracks, double pionMisidLevel)
Definition LmvmUtils.cxx:283
LmvmUtils::GetMassScaleQgp
static double GetMassScaleQgp(double minv)
Definition LmvmUtils.cxx:329
LmvmUtils::IsElectron
static void IsElectron(int globalTrackIndex, double momentum, double momentumCut, LmvmCand *cand)
Definition LmvmUtils.cxx:245
LmvmUtils::GetMcSrc
static ELmvmSrc GetMcSrc(CbmMCTrack *mctrack, TClonesArray *mcTracks)
Definition LmvmUtils.cxx:82
LmvmUtils::IsTrdElectron
static void IsTrdElectron(int globalTrackIndex, double momentum, LmvmCand *cand)
Definition LmvmUtils.cxx:272
LmvmUtils::IsRichElectron
static void IsRichElectron(int globalTrackIndex, double momentum, LmvmCand *cand)
Definition LmvmUtils.cxx:267
LmvmUtils::IsMismatch
static bool IsMismatch(const LmvmCand &cand)
Definition LmvmUtils.cxx:223
LmvmUtils::IsGhost
static bool IsGhost(const LmvmCand &cand)
Definition LmvmUtils.cxx:231
LmvmUtils::GetChargeStr
static std::string GetChargeStr(const LmvmCand *cand)
Definition LmvmUtils.cxx:296
LmvmUtils::GetBgPairSrc
static ELmvmBgPairSrc GetBgPairSrc(const LmvmCand &candP, const LmvmCand &candM)
Definition LmvmUtils.cxx:202
LmvmUtils::GetMassScaleInmed
static double GetMassScaleInmed(double minv)
Definition LmvmUtils.cxx:308
LmvmUtils::IsMcSignalEl
static bool IsMcSignalEl(const CbmMCTrack *mct)
Definition LmvmUtils.cxx:91
LmvmUtils::CalculateAndSetTrackParams
static void CalculateAndSetTrackParams(LmvmCand *cand, CbmStsTrack *stsTrack, CbmKFVertex &kfVertex)
Definition LmvmUtils.cxx:31
LmvmUtils::Distance2
static double Distance2(double x1, double y1, double x2, double y2)
Definition LmvmUtils.cxx:240
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21