CbmRoot
Loading...
Searching...
No Matches
CbmAnaConversionReco.cxx
Go to the documentation of this file.
1/* Copyright (C) 2014-2020 Fakultaet fuer Mathematik und Naturwissenschaften, Bergische Universitaet Wuppertal, Wuppertal
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Sascha Reinecke [committer], Florian Uhlig, Andrey Lebedev */
4
15#include "CbmAnaConversionReco.h"
16
17// standard includes
18#include <iostream>
19
20// includes from ROOT
21#include "TRandom3.h"
22
23// included from CbmRoot
24#include "CbmAnaConversionCutSettings.h"
25#include "CbmAnaConversionGlobalFunctions.h"
26#include "CbmGlobalTrack.h"
27#include "CbmMCTrack.h"
28#include "CbmRichElectronIdAnn.h"
29#include "CbmRichHit.h"
30#include "CbmRichRing.h"
31#include "CbmRichRingLight.h"
32#include "CbmRichUtil.h"
33
34#include "FairRootManager.h"
35
36#define M2E 2.6112004954086e-7
37using namespace std;
38
39
40CbmAnaConversionReco::CbmAnaConversionReco()
41 : fMcTracks(NULL)
42 , fGlobalTracks(NULL)
43 , fRichRings(NULL)
44 , fRichHits(NULL)
45 , fMCTracklist_all()
46 , fRecoTracklistEPEM()
47 , fRecoTracklistEPEM_ids()
48 , fRecoTracklistEPEM_chi()
49 , fRecoTracklistEPEM_gtid()
50 , fRecoMomentum()
51 , fRecoRefittedMomentum()
52 , fHistoList_MC()
53 , fHistoList_reco()
54 , fHistoList_reco_mom()
55 , fHistoList_gg()
56 , fHistoList_gee()
57 , fHistoList_eeee()
58 , fHistoList_all()
59 , fHistoList_eta()
60 , fhInvariantMass_MC_all(NULL)
61 , fhInvariantMass_MC_pi0(NULL)
62 , fhInvariantMass_MC_pi0_epem(NULL)
63 , fhInvariantMass_MC_pi0_gepem(NULL)
64 , fhInvariantMass_MC_pi0_gg(NULL)
65 , fhInvariantMass_MC_eta(NULL)
66 , fhInvariantMass_MC_etaPrime(NULL)
67 , fhMC_electrons_theta(NULL)
68 , fhMC_electrons_p(NULL)
69 , fhMC_electrons_theta_vs_p(NULL)
70 , fhEta_openingAngleGG(NULL)
71 , fhMC_grandmotherPDGs(NULL)
72 , fhInvariantMassReco_pi0(NULL)
73 , fhMCtest(NULL)
74 , fhEPEM_invmass_gg_mc(NULL)
75 , fhEPEM_invmass_gg_refitted(NULL)
76 , fhEPEM_invmass_gee_mc(NULL)
77 , fhEPEM_invmass_gee_refitted(NULL)
78 , fhEPEM_invmass_eeee_mc(NULL)
79 , fhEPEM_invmass_eeee_refitted(NULL)
80 , fhEPEM_invmass_all_mc(NULL)
81 , fhEPEM_invmass_all_refitted(NULL)
82 , fhEPEM_openingAngle_gg_mc(NULL)
83 , fhEPEM_openingAngle_gg_refitted(NULL)
84 , fhEPEM_openingAngle_gee_mc(NULL)
85 , fhEPEM_openingAngle_gee_refitted(NULL)
86 , fhEPEM_openingAngle_gee_mc_dalitz(NULL)
87 , fhEPEM_openingAngle_gee_refitted_dalitz(NULL)
88 , fhEPEM_openingAngle_vs_pt_gg_mc(NULL)
89 , fhEPEM_openingAngle_vs_pt_gg_reco(NULL)
90 , fhEPEM_openingAngle_betweenGammas_mc(NULL)
91 , fhEPEM_openingAngle_betweenGammas_reco(NULL)
92 , fhPi0_pt_vs_rap_gg(NULL)
93 , fhPi0_pt_vs_rap_gee(NULL)
94 , fhPi0_pt_vs_rap_all(NULL)
95 , fhPi0_pt_gg(NULL)
96 , fhPi0_pt_gee(NULL)
97 , fhPi0_pt_all(NULL)
98 , fhEPEM_efficiencyCuts(NULL)
99 , fhEPEM_efficiencyCuts2(NULL)
100 , fhEPEM_pi0_nofLeptons_ann(NULL)
101 , fhEPEM_pi0_ANNvalues_noCuts(NULL)
102 , fhEPEM_pi0_ANNvalues_angleCut(NULL)
103 , fhEPEM_pi0_ANNefficiencies(NULL)
104 , fhEPEM_rap_vs_chi(NULL)
105 , fhEPEM_rap_vs_invmass(NULL)
106 , fhInvMass_EPEM_mc(NULL)
107 , fhInvMass_EPEM_stsMomVec(NULL)
108 , fhInvMass_EPEM_refitted(NULL)
109 , fhInvMass_EPEM_error_stsMomVec(NULL)
110 , fhInvMass_EPEM_error_refitted(NULL)
111 , fhInvMass_EPEM_openingAngleRef(NULL)
112 , fhUsedMomenta_stsMomVec(NULL)
113 , fhUsedMomenta_mc(NULL)
114 , fhUsedMomenta_error_stsMomVec(NULL)
115 , fhUsedMomenta_error_refitted(NULL)
116 , fhUsedMomenta_errorX_stsMomVec(NULL)
117 , fhUsedMomenta_vsX_stsMomVec(NULL)
118 , fhUsedMomenta_errorY_stsMomVec(NULL)
119 , fhUsedMomenta_vsY_stsMomVec(NULL)
120 , fhUsedMomenta_errorZ_stsMomVec(NULL)
121 , fhUsedMomenta_vsZ_stsMomVec(NULL)
122 , fhUsedMomenta_errorX_refitted(NULL)
123 , fhUsedMomenta_vsX_refitted(NULL)
124 , fhUsedMomenta_errorY_refitted(NULL)
125 , fhUsedMomenta_vsY_refitted(NULL)
126 , fhUsedMomenta_errorZ_refitted(NULL)
127 , fhUsedMomenta_vsZ_refitted(NULL)
128 , fhInvariantMass_pi0epem(NULL)
129 , fhPi0_startvertex(NULL)
130 , fhPi0_startvertexElectrons_all(NULL)
131 , fhPi0_startvertexElectrons_gg(NULL)
132 , fhPi0_startvertexElectrons_gee(NULL)
133 , fhPi0_startvertex_vs_chi(NULL)
134 , fhPi0_startvertex_vs_momentum(NULL)
135 , fhInvMassWithFullRecoCuts(NULL)
136 , fhEPEM_InDetector_invmass_gg_mc(NULL)
137 , fhEPEM_InDetector_invmass_gg_refitted(NULL)
138 , fhEPEM_InDetector_invmass_gee_mc(NULL)
139 , fhEPEM_InDetector_invmass_gee_refitted(NULL)
140 , fhEPEM_InDetector_invmass_all_mc(NULL)
141 , fhEPEM_InDetector_invmass_all_refitted(NULL)
142 , fhEPEM_pt_vs_p_all_mc(NULL)
143 , fhEPEM_pt_vs_p_all_refitted(NULL)
144 , fhEPEM_missingLepton_nofRingHits(NULL)
145 , fhEPEM_missingLepton_ringMid(NULL)
146 , fhEPEM_missingLepton_ringRadius(NULL)
147 , fhEPEM_missingLepton_distance(NULL)
148 , fhEPEM_missingLepton_selectionNN(NULL)
149 , fhEPEM_missingLepton_rings(NULL)
150 , fhEPEM_missingLepton_radius_vs_p(NULL)
151 , fhEPEM_missingLepton_ANNvalue(NULL)
152 , fhEPEM_identifiedLepton_nofRingHits(NULL)
153 , fhEPEM_identifiedLepton_ringMid(NULL)
154 , fhEPEM_identifiedLepton_ringRadius(NULL)
155 , fhEPEM_identifiedLepton_distance(NULL)
156 , fhEPEM_identifiedLepton_selectionNN(NULL)
157 , fhEPEM_identifiedLepton_rings(NULL)
158 , fhEPEM_identifiedLepton_radius_vs_p(NULL)
159 , fhEPEM_identifiedLepton_ANNvalue(NULL)
160 , fhEPEM_invmass_eta_mc(NULL)
161 , fhEPEM_invmass_eta_refitted(NULL)
162 , fhEPEM_efficiencyCuts_eta(NULL)
163 , timer()
164 , fTime(0.)
165{
166}
167
168CbmAnaConversionReco::~CbmAnaConversionReco() {}
169
170
171void CbmAnaConversionReco::Init()
172{
173 FairRootManager* ioman = FairRootManager::Instance();
174 if (NULL == ioman) { Fatal("CbmAnaConversion::Init", "RootManager not instantised!"); }
175
176 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
177 if (NULL == fMcTracks) { Fatal("CbmAnaConversion::Init", "No MCTrack array!"); }
178
179 fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
180 if (NULL == fGlobalTracks) { Fatal("CbmAnaConversion::Init", "No GlobalTrack array!"); }
181
182 fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
183 if (NULL == fRichRings) { Fatal("CbmAnaConversion::Init", "No RichRing array!"); }
184
185 fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
186 if (NULL == fRichHits) { Fatal("CbmAnaConversion::Init", "No RichHit array!"); }
187
188
189 InitHistos();
190}
191
192
193void CbmAnaConversionReco::InitHistos()
194{
195 fHistoList_MC.clear();
196 fHistoList_reco.clear();
197 fHistoList_reco_mom.clear();
198
199 fHistoList_gg.clear();
200 fHistoList_gee.clear();
201 fHistoList_eeee.clear();
202 fHistoList_all.clear();
203
204 fHistoList_eta.clear();
205
206
207 Double_t invmassSpectra_nof = 800;
208 Double_t invmassSpectra_start = -0.00125;
209 Double_t invmassSpectra_end = 1.99875;
210
211
212 fhInvariantMass_MC_all =
213 new TH1D("fhInvariantMass_MC_all", "fhInvariantMass_MC_all;invariant mass in GeV/c^{2};#", 2001, -0.0005, 2.0005);
214 fhInvariantMass_MC_pi0 =
215 new TH1D("fhInvariantMass_MC_pi0", "fhInvariantMass_MC_pi0;invariant mass in GeV/c^{2};#", 2001, -0.0005, 2.0005);
216 fhInvariantMass_MC_pi0_epem = new TH1D(
217 "fhInvariantMass_MC_pi0_epem", "fhInvariantMass_MC_pi0_epem;invariant mass in GeV/c^{2};#", 2001, -0.0005, 2.0005);
218 fhInvariantMass_MC_pi0_gepem =
219 new TH1D("fhInvariantMass_MC_pi0_gepem", "fhInvariantMass_MC_pi0_gepem;invariant mass in GeV/c^{2};#", 2001,
220 -0.0005, 2.0005);
221 fhInvariantMass_MC_pi0_gg = new TH1D(
222 "fhInvariantMass_MC_pi0_gg", "fhInvariantMass_MC_pi0_gg;invariant mass in GeV/c^{2};#", 2001, -0.0005, 2.0005);
223 fhInvariantMass_MC_eta =
224 new TH1D("fhInvariantMass_MC_eta", "fhInvariantMass_MC_eta;invariant mass in GeV/c^{2};#", 2001, -0.0005, 2.0005);
225 fhInvariantMass_MC_etaPrime = new TH1D(
226 "fhInvariantMass_MC_etaPrime", "fhInvariantMass_MC_etaPrime;invariant mass in GeV/c^{2};#", 2001, -0.0005, 2.0005);
227 fHistoList_MC.push_back(fhInvariantMass_MC_all);
228 fHistoList_MC.push_back(fhInvariantMass_MC_pi0);
229 fHistoList_MC.push_back(fhInvariantMass_MC_pi0_epem);
230 fHistoList_MC.push_back(fhInvariantMass_MC_pi0_gepem);
231 fHistoList_MC.push_back(fhInvariantMass_MC_pi0_gg);
232 fHistoList_MC.push_back(fhInvariantMass_MC_eta);
233 fHistoList_MC.push_back(fhInvariantMass_MC_etaPrime);
234
235
236 fhMC_electrons_theta = new TH1D("fhMC_electrons_theta", "fhMC_electrons_theta;theta in deg;#", 90, 0., 90.);
237 fhMC_electrons_p = new TH1D("fhMC_electrons_p", "fhMC_electrons_p;momentum p in GeV/c;#", 100, 0., 10.);
238 fhMC_electrons_theta_vs_p =
239 new TH2D("fhMC_electrons_theta_vs_p", "fhMC_electrons_theta_vs_p;theta in deg;momentum p in GeV/c", 90, 0., 90.,
240 100, 0., 10.);
241 fHistoList_MC.push_back(fhMC_electrons_theta);
242 fHistoList_MC.push_back(fhMC_electrons_p);
243 fHistoList_MC.push_back(fhMC_electrons_theta_vs_p);
244
245 fhEta_openingAngleGG = new TH1D("fhEta_openingAngleGG", "fhEta_openingAngleGG;theta in deg;#", 900, 0., 90.);
246 fHistoList_MC.push_back(fhEta_openingAngleGG);
247
248 fhMC_grandmotherPDGs = new TH1D("fhMC_grandmotherPDGs", "fhMC_grandmotherPDGs;particle pdg;#", 1000, 0., 1000.);
249 fHistoList_MC.push_back(fhMC_grandmotherPDGs);
250
251 fhMCtest = new TH1D("fhMCtest", "fhMCtest;invariant mass in GeV/c^{2};#", 2000, 0., 2.);
252 fHistoList_MC.push_back(fhMCtest);
253
254
255 fhEPEM_invmass_gg_mc = new TH1D("fhEPEM_invmass_gg_mc", "fhEPEM_invmass_gg_mc;invariant mass in GeV/c^{2};#",
256 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
257 fhEPEM_invmass_gg_refitted =
258 new TH1D("fhEPEM_invmass_gg_refitted", "fhEPEM_invmass_gg_refitted;invariant mass in GeV/c^{2};#",
259 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
260 fhEPEM_invmass_gee_mc = new TH1D("fhEPEM_invmass_gee_mc", "fhEPEM_invmass_gee_mc;invariant mass in GeV/c^{2};#",
261 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
262 fhEPEM_invmass_gee_refitted =
263 new TH1D("fhEPEM_invmass_gee_refitted", "fhEPEM_invmass_gee_refitted;invariant mass in GeV/c^{2};#",
264 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
265 fhEPEM_invmass_eeee_mc = new TH1D("fhEPEM_invmass_eeee_mc", "fhEPEM_invmass_eeee_mc;invariant mass in GeV/c^{2};#",
266 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
267 fhEPEM_invmass_eeee_refitted =
268 new TH1D("fhEPEM_invmass_eeee_refitted", "fhEPEM_invmass_eeee_refitted;invariant mass in GeV/c^{2};#",
269 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
270 fhEPEM_invmass_all_mc = new TH1D("fhEPEM_invmass_all_mc", "fhEPEM_invmass_all_mc;invariant mass in GeV/c^{2};#",
271 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
272 fhEPEM_invmass_all_refitted =
273 new TH1D("fhEPEM_invmass_all_refitted", "fhEPEM_invmass_all_refitted;invariant mass in GeV/c^{2};#",
274 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
275 fHistoList_gg.push_back(fhEPEM_invmass_gg_mc);
276 fHistoList_gg.push_back(fhEPEM_invmass_gg_refitted);
277 fHistoList_gee.push_back(fhEPEM_invmass_gee_mc);
278 fHistoList_gee.push_back(fhEPEM_invmass_gee_refitted);
279 fHistoList_eeee.push_back(fhEPEM_invmass_eeee_mc);
280 fHistoList_eeee.push_back(fhEPEM_invmass_eeee_refitted);
281 fHistoList_all.push_back(fhEPEM_invmass_all_mc);
282 fHistoList_all.push_back(fhEPEM_invmass_all_refitted);
283
284 fhEPEM_openingAngle_gg_mc = new TH1D(
285 "fhEPEM_openingAngle_gg_mc", "fhEPEM_openingAngle_gg_mc (between e+e- from g);angle [deg];#", 1010, -0.1, 100.9);
286 fhEPEM_openingAngle_gg_refitted =
287 new TH1D("fhEPEM_openingAngle_gg_refitted", "fhEPEM_openingAngle_gg_refitted (between e+e- from g);angle [deg];#",
288 1010, -0.1, 100.9);
289 fhEPEM_openingAngle_gee_mc = new TH1D(
290 "fhEPEM_openingAngle_gee_mc", "fhEPEM_openingAngle_gee_mc (between e+e- from g);angle [deg];#", 1010, -0.1, 100.9);
291 fhEPEM_openingAngle_gee_refitted =
292 new TH1D("fhEPEM_openingAngle_gee_refitted", "fhEPEM_openingAngle_gee_refitted (between e+e- from g);angle [deg];#",
293 1010, -0.1, 100.9);
294 fHistoList_gg.push_back(fhEPEM_openingAngle_gg_mc);
295 fHistoList_gg.push_back(fhEPEM_openingAngle_gg_refitted);
296 fHistoList_gee.push_back(fhEPEM_openingAngle_gee_mc);
297 fHistoList_gee.push_back(fhEPEM_openingAngle_gee_refitted);
298 fhEPEM_openingAngle_gee_mc_dalitz =
299 new TH1D("fhEPEM_openingAngle_gee_mc_dalitz",
300 "fhEPEM_openingAngle_gee_mc_dalitz (between e+e- from pi0);angle [deg];#", 1010, -0.1, 100.9);
301 fhEPEM_openingAngle_gee_refitted_dalitz = new TH1D("fhEPEM_openingAngle_gee_refitted_dalitz",
302 "fhEPEM_openingAngle_gee_refitted_dalitz (between e+e- from "
303 "pi0);angle [deg];#",
304 1010, -0.1, 100.9);
305 fHistoList_gee.push_back(fhEPEM_openingAngle_gee_mc_dalitz);
306 fHistoList_gee.push_back(fhEPEM_openingAngle_gee_refitted_dalitz);
307
308 fhEPEM_openingAngle_vs_pt_gg_mc =
309 new TH2D("fhEPEM_openingAngle_vs_pt_gg_mc", "fhEPEM_openingAngle_vs_pt_gg_mc;pt [GeV]; opening angle [deg]", 220,
310 -1., 10., 100, 0., 10.);
311 fHistoList_gg.push_back(fhEPEM_openingAngle_vs_pt_gg_mc);
312 fhEPEM_openingAngle_vs_pt_gg_reco =
313 new TH2D("fhEPEM_openingAngle_vs_pt_gg_reco", "fhEPEM_openingAngle_vs_pt_gg_reco;pt [GeV]; opening angle [deg]",
314 220, -1., 10., 100, 0., 10.);
315 fHistoList_gg.push_back(fhEPEM_openingAngle_vs_pt_gg_reco);
316
317 fhEPEM_openingAngle_betweenGammas_mc = new TH1D(
318 "fhEPEM_openingAngle_betweenGammas_mc", "fhEPEM_openingAngle_betweenGammas_mc;angle [deg];#", 1010, -0.1, 100.9);
319 fhEPEM_openingAngle_betweenGammas_reco =
320 new TH1D("fhEPEM_openingAngle_betweenGammas_reco", "fhEPEM_openingAngle_betweenGammas_reco;angle [deg];#", 1010,
321 -0.1, 100.9);
322 fHistoList_gg.push_back(fhEPEM_openingAngle_betweenGammas_mc);
323 fHistoList_gg.push_back(fhEPEM_openingAngle_betweenGammas_reco);
324
325
326 fhPi0_pt_vs_rap_gg =
327 new TH2D("fhPi0_pt_vs_rap_gg", "fhPi0_pt_vs_rap_gg;pt [GeV]; rap [GeV]", 240, -2., 10., 270, -2., 7.);
328 fhPi0_pt_vs_rap_gee =
329 new TH2D("fhPi0_pt_vs_rap_gee", "fhPi0_pt_vs_rap_gee;pt [GeV]; rap [GeV]", 240, -2., 10., 270, -2., 7.);
330 fhPi0_pt_vs_rap_all =
331 new TH2D("fhPi0_pt_vs_rap_all", "fhPi0_pt_vs_rap_all;pt [GeV]; rap [GeV]", 240, -2., 10., 270, -2., 7.);
332 fHistoList_gg.push_back(fhPi0_pt_vs_rap_gg);
333 fHistoList_gee.push_back(fhPi0_pt_vs_rap_gee);
334 fHistoList_all.push_back(fhPi0_pt_vs_rap_all);
335
336 fhPi0_pt_gg = new TH1D("fhPi0_pt_gg", "fhPi0_pt_gg;pt [GeV];#", 200, 0., 10.);
337 fhPi0_pt_gee = new TH1D("fhPi0_pt_gee", "fhPi0_pt_gee;pt [GeV];#", 200, 0., 10.);
338 fhPi0_pt_all = new TH1D("fhPi0_pt_all", "fhPi0_pt_all;pt [GeV];#", 200, 0., 10.);
339 fHistoList_gg.push_back(fhPi0_pt_gg);
340 fHistoList_gee.push_back(fhPi0_pt_gee);
341 fHistoList_all.push_back(fhPi0_pt_all);
342
343 fhEPEM_efficiencyCuts = new TH1D("fhEPEM_efficiencyCuts", "fhEPEM_efficiencyCuts;;#", 13, 0., 13.);
344 fHistoList_all.push_back(fhEPEM_efficiencyCuts);
345 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(1, "no cuts");
346 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(2, "ANN: 4 rich electrons");
347 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(3, "ANN: #chi^{2}-cut");
348 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(4, "ANN: #theta of e^{+}e^{-} pairs");
349 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(5, "ANN: m_{inv} of e^{+}e^{-} pairs");
350 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(6, "Normal: 4 rich electrons");
351 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(7, "Normal: #chi^{2}-cut");
352 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(8, "Normal: #theta of e^{+}e^{-} pairs");
353 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(9, "Normal: m_{inv} of e^{+}e^{-} pairs");
354 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(10, "MC: 4 rich electrons");
355 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(11, "MC: #chi^{2}-cut");
356 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(12, "MC: #theta of e^{+}e^{-} pairs");
357 fhEPEM_efficiencyCuts->GetXaxis()->SetBinLabel(13, "MC: m_{inv} of e^{+}e^{-} pairs");
358
359 fhEPEM_efficiencyCuts2 = new TH1D("fhEPEM_efficiencyCuts2", "fhEPEM_efficiencyCuts2;;#", 10, 0., 10.);
360 fHistoList_all.push_back(fhEPEM_efficiencyCuts2);
361
362 fhEPEM_pi0_nofLeptons_ann = new TH1D("fhEPEM_pi0_nofLeptons_ann", "fhEPEM_pi0_nofLeptons_ann;;#", 5, -0.5, 4.5);
363 fhEPEM_pi0_ANNvalues_noCuts = new TH1D("fhEPEM_pi0_ANNvalues_noCuts", "fhEPEM_pi0_ANNvalues_noCuts;;#", 400, -2, 2);
364 fhEPEM_pi0_ANNvalues_angleCut =
365 new TH1D("fhEPEM_pi0_ANNvalues_angleCut", "fhEPEM_pi0_ANNvalues_angleCut;;#", 400, -2, 2);
366 fhEPEM_pi0_ANNefficiencies = new TH1D("fhEPEM_pi0_ANNefficiencies", "fhEPEM_pi0_ANNefficiencies;;#", 15, 0, 15);
367 fHistoList_all.push_back(fhEPEM_pi0_ANNvalues_noCuts);
368 fHistoList_all.push_back(fhEPEM_pi0_ANNvalues_angleCut);
369 fHistoList_all.push_back(fhEPEM_pi0_nofLeptons_ann);
370 fHistoList_all.push_back(fhEPEM_pi0_ANNefficiencies);
371
372 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(1, "no cuts");
373 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(2, "-1.0");
374 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(3, "-0.9");
375 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(4, "-0.8");
376 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(5, "-0.7");
377 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(6, "-0.6");
378 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(7, "-0.5");
379 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(8, "-0.4");
380 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(9, "-0.3");
381 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(10, "-0.2");
382 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(11, "-0.1");
383 fhEPEM_pi0_ANNefficiencies->GetXaxis()->SetBinLabel(12, "0.0");
384
385 fhEPEM_rap_vs_chi =
386 new TH2D("fhEPEM_rap_vs_chi", "fhEPEM_rap_vs_chi; rap [GeV]; chi of electrons", 300, 0., 10., 100, 0., 100.);
387 fHistoList_all.push_back(fhEPEM_rap_vs_chi);
388 fhEPEM_rap_vs_invmass =
389 new TH2D("fhEPEM_rap_vs_invmass", "fhEPEM_rap_vs_invmass; rap [GeV]; invmass", 300, 0., 10., 100, 0., 10.);
390 fHistoList_all.push_back(fhEPEM_rap_vs_invmass);
391
392
393 fhInvMass_EPEM_mc = new TH1D("fhInvMass_EPEM_mc", "fhInvariantMass_recoMomentum1 (mc);mass [GeV/c^2];#",
394 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
395 fhInvMass_EPEM_stsMomVec =
396 new TH1D("fhInvMass_EPEM_stsMomVec", "fhInvariantMass_recoMomentum2 (stsMomentumVec);mass [GeV/c^2];#",
397 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
398 fhInvMass_EPEM_refitted =
399 new TH1D("fhInvMass_EPEM_refitted", "fhInvariantMass_recoMomentum3 (refitted at primary);mass [GeV/c^2];#",
400 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
401 fhInvMass_EPEM_error_stsMomVec =
402 new TH1D("fhInvMass_EPEM_error_stsMomVec", "fhInvariantMass_recoMomentum4 (error, stsMomentumVec);(mc-reco)/mc;#",
403 500, -0.005, 4.995);
404 fhInvMass_EPEM_error_refitted =
405 new TH1D("fhInvMass_EPEM_error_refitted", "fhInvariantMass_recoMomentum5 (error, refitted);(mc-reco)/mc;#", 500,
406 -0.005, 4.995);
407 fhInvMass_EPEM_openingAngleRef =
408 new TH1D("fhInvMass_EPEM_openingAngleRef", "fhInvariantMass_openingAngleRef;angle [deg];#", 1010, -0.1, 100.9);
409 fhUsedMomenta_stsMomVec = new TH1D("fhUsedMomenta_stsMomVec", "fhMomentumtest1;momentum;#", 100, 0., 2.);
410 fhUsedMomenta_mc = new TH1D("fhUsedMomenta_mc", "fhMomentumtest2;momentum;#", 100, 0., 2.);
411 fhUsedMomenta_error_stsMomVec =
412 new TH1D("fhUsedMomenta_error_stsMomVec", "fhMomentumtest3 (error);(mc-reco)/mc;#", 400, -2.005, 1.995);
413 fhUsedMomenta_error_refitted =
414 new TH1D("fhUsedMomenta_error_refitted", "fhMomentumtest4 (error);(mc-reco_refitted)/mc;#", 400, -2.005, 1.995);
415
416 fhUsedMomenta_errorX_stsMomVec = new TH1D(
417 "fhUsedMomenta_errorX_stsMomVec", "fhMomentumtest5 (error of x-momentum);(mc-reco_reco)/mc;#", 401, -4.01, 4.01);
418 fhUsedMomenta_vsX_stsMomVec =
419 new TH2D("fhUsedMomenta_vsX_stsMomVec", "fhMomentumtest5vs (error of x-momentum);mc;reco", 101, -1.01, 1.01, 101,
420 -1.01, 1.01);
421 fhUsedMomenta_errorY_stsMomVec = new TH1D(
422 "fhUsedMomenta_errorY_stsMomVec", "fhMomentumtest6 (error of y-momentum);(mc-reco_reco)/mc;#", 401, -4.01, 4.01);
423 fhUsedMomenta_vsY_stsMomVec =
424 new TH2D("fhUsedMomenta_vsY_stsMomVec", "fhMomentumtest6vs (error of y-momentum);mc;reco", 101, -1.01, 1.01, 101,
425 -1.01, 1.01);
426 fhUsedMomenta_errorZ_stsMomVec = new TH1D(
427 "fhUsedMomenta_errorZ_stsMomVec", "fhMomentumtest7 (error of z-momentum);(mc-reco_reco)/mc;#", 401, -4.01, 4.01);
428 fhUsedMomenta_vsZ_stsMomVec =
429 new TH2D("fhUsedMomenta_vsZ_stsMomVec", "fhMomentumtest7vs (error of z-momentum);mc;reco", 201, -0.01, 4.01, 201,
430 -0.01, 4.01);
431 fhUsedMomenta_errorX_refitted = new TH1D(
432 "fhUsedMomenta_errorX_refitted", "fhMomentumtest5 (error of x-momentum);(mc-reco_reco)/mc;#", 401, -4.01, 4.01);
433 fhUsedMomenta_vsX_refitted = new TH2D("fhUsedMomenta_vsX_refitted", "fhMomentumtest5vs (error of x-momentum);mc;reco",
434 101, -1.01, 1.01, 101, -1.01, 1.01);
435 fhUsedMomenta_errorY_refitted = new TH1D(
436 "fhUsedMomenta_errorY_refitted", "fhMomentumtest6 (error of y-momentum);(mc-reco_reco)/mc;#", 401, -4.01, 4.01);
437 fhUsedMomenta_vsY_refitted = new TH2D("fhUsedMomenta_vsY_refitted", "fhMomentumtest6vs (error of y-momentum);mc;reco",
438 101, -1.01, 1.01, 101, -1.01, 1.01);
439 fhUsedMomenta_errorZ_refitted = new TH1D(
440 "fhUsedMomenta_errorZ_refitted", "fhMomentumtest7 (error of z-momentum);(mc-reco_reco)/mc;#", 401, -4.01, 4.01);
441 fhUsedMomenta_vsZ_refitted = new TH2D("fhUsedMomenta_vsZ_refitted", "fhMomentumtest7vs (error of z-momentum);mc;reco",
442 201, -0.01, 4.01, 201, -0.01, 4.01);
443 fHistoList_reco.push_back(fhInvMass_EPEM_mc);
444 fHistoList_reco.push_back(fhInvMass_EPEM_stsMomVec);
445 fHistoList_reco.push_back(fhInvMass_EPEM_refitted);
446 fHistoList_reco.push_back(fhInvMass_EPEM_error_stsMomVec);
447 fHistoList_reco.push_back(fhInvMass_EPEM_error_refitted);
448 fHistoList_reco.push_back(fhInvMass_EPEM_openingAngleRef);
449 fHistoList_reco_mom.push_back(fhUsedMomenta_stsMomVec);
450 fHistoList_reco_mom.push_back(fhUsedMomenta_mc);
451 fHistoList_reco_mom.push_back(fhUsedMomenta_error_stsMomVec);
452 fHistoList_reco_mom.push_back(fhUsedMomenta_error_refitted);
453 fHistoList_reco_mom.push_back(fhUsedMomenta_errorX_stsMomVec); // error of x-component of reconstructed momentum
454 fHistoList_reco_mom.push_back(fhUsedMomenta_errorY_stsMomVec); // error of y-component of reconstructed momentum
455 fHistoList_reco_mom.push_back(fhUsedMomenta_errorZ_stsMomVec); // error of z-component of reconstructed momentum
456 fHistoList_reco_mom.push_back(fhUsedMomenta_vsX_stsMomVec); // x-component of reconstructed momentum vs mc-momentum
457 fHistoList_reco_mom.push_back(fhUsedMomenta_vsY_stsMomVec); // y-component of reconstructed momentum vs mc-momentum
458 fHistoList_reco_mom.push_back(fhUsedMomenta_vsZ_stsMomVec); // z-component of reconstructed momentum vs mc-momentum
459 fHistoList_reco_mom.push_back(fhUsedMomenta_errorX_refitted); // error of x-component of reconstructed momentum
460 fHistoList_reco_mom.push_back(fhUsedMomenta_errorY_refitted); // error of y-component of reconstructed momentum
461 fHistoList_reco_mom.push_back(fhUsedMomenta_errorZ_refitted); // error of z-component of reconstructed momentum
462 fHistoList_reco_mom.push_back(fhUsedMomenta_vsX_refitted); // x-component of reconstructed momentum vs mc-momentum
463 fHistoList_reco_mom.push_back(fhUsedMomenta_vsY_refitted); // y-component of reconstructed momentum vs mc-momentum
464 fHistoList_reco_mom.push_back(fhUsedMomenta_vsZ_refitted); // z-component of reconstructed momentum vs mc-momentum
465
466 fhInvariantMass_pi0epem =
467 new TH1D("fhInvariantMass_pi0epem", "fhInvariantMass_pi0epem;mass [GeV/c^2];#", 400, 0., 2.);
468 fHistoList_reco.push_back(fhInvariantMass_pi0epem);
469
470 fhPi0_startvertex = new TH1D("fhPi0_startvertex", "fhPi0_startvertex;z in cm;#", 210, -5., 100.);
471 fHistoList_reco.push_back(fhPi0_startvertex);
472
473 fhPi0_startvertexElectrons_all =
474 new TH1D("fhPi0_startvertexElectrons_all", "fhPi0_startvertexElectrons_all;z in cm;#", 411, -5.25, 200.25);
475 fHistoList_reco.push_back(fhPi0_startvertexElectrons_all);
476
477 fhPi0_startvertexElectrons_gg =
478 new TH1D("fhPi0_startvertexElectrons_gg", "fhPi0_startvertexElectrons_gg;z in cm;#", 411, -5.25, 200.25);
479 fHistoList_reco.push_back(fhPi0_startvertexElectrons_gg);
480
481 fhPi0_startvertexElectrons_gee =
482 new TH1D("fhPi0_startvertexElectrons_gee", "fhPi0_startvertexElectrons_gee;z in cm;#", 411, -5.25, 200.25);
483 fHistoList_reco.push_back(fhPi0_startvertexElectrons_gee);
484
485 fhPi0_startvertex_vs_chi =
486 new TH2D("fhPi0_startvertex_vs_chi", "fhPi0_startvertex_vs_chi;z[cm];chi", 210, -5., 100., 1000, 0., 100.);
487 fHistoList_reco.push_back(fhPi0_startvertex_vs_chi);
488
489 fhPi0_startvertex_vs_momentum =
490 new TH2D("fhPi0_startvertex_vs_momentum", "fhPi0_startvertex_vs_momentum;z in cm;momentum (MC-true)", 210, -5.,
491 100., 1000, 0., 100.);
492 fHistoList_reco.push_back(fhPi0_startvertex_vs_momentum);
493
494 fhInvMassWithFullRecoCuts =
495 new TH1D("fhInvMassWithFullRecoCuts", "fhInvMassWithFullRecoCuts;mass [GeV/c^2];#", 800, 0., 2.);
496 fHistoList_reco.push_back(fhInvMassWithFullRecoCuts);
497
498
499 fhEPEM_InDetector_invmass_gg_mc =
500 new TH1D("fhEPEM_InDetector_invmass_gg_mc", "fhEPEM_InDetector_invmass_gg_mc;mass [GeV/c^2];#", invmassSpectra_nof,
501 invmassSpectra_start, invmassSpectra_end);
502 fhEPEM_InDetector_invmass_gg_refitted =
503 new TH1D("fhEPEM_InDetector_invmass_gg_refitted", "fhEPEM_InDetector_invmass_gg_refitted;mass [GeV/c^2];#",
504 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
505 fhEPEM_InDetector_invmass_gee_mc =
506 new TH1D("fhEPEM_InDetector_invmass_gee_mc", "fhEPEM_InDetector_invmass_gee_mc;mass [GeV/c^2];#",
507 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
508 fhEPEM_InDetector_invmass_gee_refitted =
509 new TH1D("fhEPEM_InDetector_invmass_gee_refitted", "fhEPEM_InDetector_invmass_gee_refitted;mass [GeV/c^2];#",
510 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
511 fhEPEM_InDetector_invmass_all_mc =
512 new TH1D("fhEPEM_InDetector_invmass_all_mc", "fhEPEM_InDetector_invmass_all_mc;mass [GeV/c^2];#",
513 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
514 fhEPEM_InDetector_invmass_all_refitted =
515 new TH1D("fhEPEM_InDetector_invmass_all_refitted", "fhEPEM_InDetector_invmass_all_refitted;mass [GeV/c^2];#",
516 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
517 fHistoList_gg.push_back(fhEPEM_InDetector_invmass_gg_mc);
518 fHistoList_gg.push_back(fhEPEM_InDetector_invmass_gg_refitted);
519 fHistoList_gee.push_back(fhEPEM_InDetector_invmass_gee_mc);
520 fHistoList_gee.push_back(fhEPEM_InDetector_invmass_gee_refitted);
521 fHistoList_all.push_back(fhEPEM_InDetector_invmass_all_mc);
522 fHistoList_all.push_back(fhEPEM_InDetector_invmass_all_refitted);
523
524
525 fhEPEM_pt_vs_p_all_mc =
526 new TH2D("fhEPEM_pt_vs_p_all_mc", "fhEPEM_pt_vs_p_all_mc;p_{t} in GeV/c; p in GeV/c", 240, -2., 10., 360, -2., 16.);
527 fHistoList_all.push_back(fhEPEM_pt_vs_p_all_mc);
528 fhEPEM_pt_vs_p_all_refitted =
529 new TH2D("fhEPEM_pt_vs_p_all_refitted", "fhTest2_electrons_pt_vs_p;p_{t} in GeV/c; p in GeV/c", 240, -2., 10., 360,
530 -2., 16.);
531 fHistoList_all.push_back(fhEPEM_pt_vs_p_all_refitted);
532
533 fhEPEM_missingLepton_nofRingHits =
534 new TH1D("fhEPEM_missingLepton_nofRingHits", "fhEPEM_missingLepton_nofRingHits;nofringhits;#", 41, -0.5, 40.5);
535 fhEPEM_missingLepton_ringMid =
536 new TH2D("fhEPEM_missingLepton_ringMid", "fhEPEM_missingLepton_ringMid;X;Y", 100, -150, 150, 150, -250, 250);
537 fhEPEM_missingLepton_ringRadius =
538 new TH1D("fhEPEM_missingLepton_ringRadius", "fhEPEM_missingLepton_ringRadius;ringRadius;#", 100, 0, 10);
539 fhEPEM_missingLepton_distance =
540 new TH1D("fhEPEM_missingLepton_distance", "fhEPEM_missingLepton_distance;distance;#", 100, 0, 20);
541 fhEPEM_missingLepton_selectionNN =
542 new TH1D("fhEPEM_missingLepton_selectionNN", "fhEPEM_missingLepton_selectionNN;selectionNN;#", 40, -2, 2);
543 fhEPEM_missingLepton_rings =
544 new TH2D("fhEPEM_missingLepton_rings", "fhEPEM_missingLepton_rings;selectionNN;#", 400, -200, 200, 600, -300, 300);
545 fhEPEM_missingLepton_radius_vs_p =
546 new TH2D("fhEPEM_missingLepton_radius_vs_p",
547 "fhEPEM_missingLepton_radius_vs_p;momentum p in GeV/c;ring radius in cm", 120, 0, 12, 100, 0, 10);
548 fhEPEM_missingLepton_ANNvalue =
549 new TH1D("fhEPEM_missingLepton_ANNvalue", "fhEPEM_missingLepton_ANNvalue;ANNvalue;#", 40, -2, 2);
550 fHistoList_all.push_back(fhEPEM_missingLepton_nofRingHits);
551 fHistoList_all.push_back(fhEPEM_missingLepton_ringMid);
552 fHistoList_all.push_back(fhEPEM_missingLepton_ringRadius);
553 fHistoList_all.push_back(fhEPEM_missingLepton_distance);
554 fHistoList_all.push_back(fhEPEM_missingLepton_selectionNN);
555 fHistoList_all.push_back(fhEPEM_missingLepton_rings);
556 fHistoList_all.push_back(fhEPEM_missingLepton_radius_vs_p);
557 fHistoList_all.push_back(fhEPEM_missingLepton_ANNvalue);
558
559
560 fhEPEM_identifiedLepton_nofRingHits = new TH1D("fhEPEM_identifiedLepton_nofRingHits",
561 "fhEPEM_identifiedLepton_nofRingHits;nofringhits;#", 41, -0.5, 40.5);
562 fhEPEM_identifiedLepton_ringMid =
563 new TH2D("fhEPEM_identifiedLepton_ringMid", "fhEPEM_identifiedLepton_ringMid;X;Y", 100, -150, 150, 150, -250, 250);
564 fhEPEM_identifiedLepton_ringRadius =
565 new TH1D("fhEPEM_identifiedLepton_ringRadius", "fhEPEM_identifiedLepton_ringRadius;ringRadius;#", 100, 0, 10);
566 fhEPEM_identifiedLepton_distance =
567 new TH1D("fhEPEM_identifiedLepton_distance", "fhEPEM_identifiedLepton_distance;distance;#", 100, 0, 20);
568 fhEPEM_identifiedLepton_selectionNN =
569 new TH1D("fhEPEM_identifiedLepton_selectionNN", "fhEPEM_identifiedLepton_selectionNN;selectionNN;#", 40, -2, 2);
570 fhEPEM_identifiedLepton_rings = new TH2D(
571 "fhEPEM_identifiedLepton_rings", "fhEPEM_identifiedLepton_rings;selectionNN;#", 400, -200, 200, 600, -300, 300);
572 fhEPEM_identifiedLepton_radius_vs_p =
573 new TH2D("fhEPEM_identifiedLepton_radius_vs_p",
574 "fhEPEM_identifiedLepton_radius_vs_p;momentum p in GeV/c;ring radius in cm", 120, 0, 12, 100, 0, 10);
575 fhEPEM_identifiedLepton_ANNvalue =
576 new TH1D("fhEPEM_identifiedLepton_ANNvalue", "fhEPEM_identifiedLepton_ANNvalue;ANNvalue;#", 40, -2, 2);
577 fHistoList_all.push_back(fhEPEM_identifiedLepton_nofRingHits);
578 fHistoList_all.push_back(fhEPEM_identifiedLepton_ringMid);
579 fHistoList_all.push_back(fhEPEM_identifiedLepton_ringRadius);
580 fHistoList_all.push_back(fhEPEM_identifiedLepton_distance);
581 fHistoList_all.push_back(fhEPEM_identifiedLepton_selectionNN);
582 fHistoList_all.push_back(fhEPEM_identifiedLepton_rings);
583 fHistoList_all.push_back(fhEPEM_identifiedLepton_radius_vs_p);
584 fHistoList_all.push_back(fhEPEM_identifiedLepton_ANNvalue);
585
586
587 // histograms for eta analysis
588 fhEPEM_invmass_eta_mc = new TH1D("fhEPEM_invmass_eta_mc", "fhEPEM_Iinvmass_eta_mc;mass [GeV/c^2];#",
589 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
590 fhEPEM_invmass_eta_refitted = new TH1D("fhEPEM_invmass_eta_refitted", "fhEPEM_Iinvmass_eta_refitted;mass [GeV/c^2];#",
591 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
592 fHistoList_eta.push_back(fhEPEM_invmass_eta_mc);
593 fHistoList_eta.push_back(fhEPEM_invmass_eta_refitted);
594
595 fhEPEM_efficiencyCuts_eta = new TH1D("fhEPEM_efficiencyCuts_eta", "fhEPEM_efficiencyCuts_eta;;#", 13, 0., 13.);
596 fHistoList_eta.push_back(fhEPEM_efficiencyCuts_eta);
597 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(1, "no cuts");
598 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(2, "ANN: 4 rich electrons");
599 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(3, "ANN: #chi^{2}-cut");
600 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(4, "ANN: #theta of e^{+}e^{-} pairs");
601 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(5, "ANN: m_{inv} of e^{+}e^{-} pairs");
602 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(6, "Normal: 4 rich electrons");
603 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(7, "Normal: #chi^{2}-cut");
604 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(8, "Normal: #theta of e^{+}e^{-} pairs");
605 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(9, "Normal: m_{inv} of e^{+}e^{-} pairs");
606 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(10, "MC: 4 rich electrons");
607 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(11, "MC: #chi^{2}-cut");
608 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(12, "MC: #theta of e^{+}e^{-} pairs");
609 fhEPEM_efficiencyCuts_eta->GetXaxis()->SetBinLabel(13, "MC: m_{inv} of e^{+}e^{-} pairs");
610}
611
612
613void CbmAnaConversionReco::Finish()
614{
615 //gDirectory->cd("analysis-conversion");
616 gDirectory->mkdir("MCreco");
617 gDirectory->cd("MCreco");
618 for (UInt_t i = 0; i < fHistoList_MC.size(); i++) {
619 fHistoList_MC[i]->Write();
620 }
621 gDirectory->cd("..");
622
623 gDirectory->mkdir("Reconstruction2");
624 gDirectory->cd("Reconstruction2");
625
626 gDirectory->mkdir("pi0 -> gg");
627 gDirectory->cd("pi0 -> gg");
628 for (UInt_t i = 0; i < fHistoList_gg.size(); i++) {
629 fHistoList_gg[i]->Write();
630 }
631 gDirectory->cd("..");
632
633 gDirectory->mkdir("pi0 -> ge+e-");
634 gDirectory->cd("pi0 -> ge+e-");
635 for (UInt_t i = 0; i < fHistoList_gee.size(); i++) {
636 fHistoList_gee[i]->Write();
637 }
638 gDirectory->cd("..");
639
640 gDirectory->mkdir("pi0 -> e+e-e+e-");
641 gDirectory->cd("pi0 -> e+e-e+e-");
642 for (UInt_t i = 0; i < fHistoList_eeee.size(); i++) {
643 fHistoList_eeee[i]->Write();
644 }
645 gDirectory->cd("..");
646
647 gDirectory->mkdir("pi0 -> all");
648 gDirectory->cd("pi0 -> all");
649 for (UInt_t i = 0; i < fHistoList_all.size(); i++) {
650 fHistoList_all[i]->Write();
651 }
652 gDirectory->cd("..");
653
654 gDirectory->mkdir("eta");
655 gDirectory->cd("eta");
656 for (UInt_t i = 0; i < fHistoList_eta.size(); i++) {
657 fHistoList_eta[i]->Write();
658 }
659 gDirectory->cd("..");
660
661 for (UInt_t i = 0; i < fHistoList_reco.size(); i++) {
662 fHistoList_reco[i]->Write();
663 }
664 gDirectory->mkdir("Momenta2");
665 gDirectory->cd("Momenta2");
666 for (UInt_t i = 0; i < fHistoList_reco_mom.size(); i++) {
667 fHistoList_reco_mom[i]->Write();
668 }
669 gDirectory->cd("../..");
670
671 cout << "CbmAnaConversionReco: Realtime - " << fTime << endl;
672 //timer.Print();
673}
674
675
676void CbmAnaConversionReco::SetTracklistMC(vector<CbmMCTrack*> MCTracklist) { fMCTracklist_all = MCTracklist; }
677
678
679void CbmAnaConversionReco::SetTracklistReco(vector<CbmMCTrack*> MCTracklist, vector<TVector3> RecoTracklist1,
680 vector<TVector3> RecoTracklist2, vector<int> ids, vector<Double_t> chi,
681 vector<Int_t> GlobalTrackId)
682{
683 fRecoTracklistEPEM = MCTracklist;
684 fRecoTracklistEPEM_ids = ids;
685 fRecoTracklistEPEM_chi = chi;
686 fRecoTracklistEPEM_gtid = GlobalTrackId;
687 fRecoMomentum = RecoTracklist1;
688 fRecoRefittedMomentum = RecoTracklist2;
689}
690
691
692void CbmAnaConversionReco::InvariantMassMC_all()
693// calculation of invariant mass via X -> gamma gamma -> e+ e- e+ e-, only MC data with several cuts (see FillMCTrackslists())
694{
695 timer.Start();
696
697
698 cout << "CbmAnaConversionReco: InvariantMassTestMC - Start..." << endl;
699 cout << "CbmAnaConversionReco: InvariantMassTestMC - Size of "
700 "fTracklistMC_all:\t "
701 << fMCTracklist_all.size() << endl;
702 if (fMCTracklist_all.size() >= 4) {
703 for (unsigned int i = 0; i < fMCTracklist_all.size() - 3; i++) {
704 if (i % 10 == 0) cout << "CbmAnaConversionReco: InvariantMassTestMC - iteration i = " << i << endl;
705 for (unsigned int j = i + 1; j < fMCTracklist_all.size() - 2; j++) {
706 for (unsigned int k = j + 1; k < fMCTracklist_all.size() - 1; k++) {
707 for (unsigned int l = k + 1; l < fMCTracklist_all.size(); l++) {
708
709 if (fMCTracklist_all[i]->GetPdgCode() + fMCTracklist_all[j]->GetPdgCode()
710 + fMCTracklist_all[k]->GetPdgCode() + fMCTracklist_all[l]->GetPdgCode()
711 != 0)
712 continue;
713
714 if (i == j || i == k || i == l || j == k || j == l || k == l) continue;
715
716 int motherId1 = fMCTracklist_all[i]->GetMotherId();
717 int motherId2 = fMCTracklist_all[j]->GetMotherId();
718 int motherId3 = fMCTracklist_all[k]->GetMotherId();
719 int motherId4 = fMCTracklist_all[l]->GetMotherId();
720
721 TVector3 momentum1, momentum2, momentum3, momentum4;
722 fMCTracklist_all[i]->GetMomentum(momentum1);
723 fMCTracklist_all[j]->GetMomentum(momentum2);
724 fMCTracklist_all[k]->GetMomentum(momentum3);
725 fMCTracklist_all[l]->GetMomentum(momentum4);
726
727
728 // decay pi0 -> e+ e- e+ e-
729 if (motherId1 == motherId2 && motherId1 == motherId3 && motherId1 == motherId4) {
730 cout << "testxyz" << endl;
731 Double_t invmass =
732 Invmass_4particles(fMCTracklist_all[i], fMCTracklist_all[j], fMCTracklist_all[k], fMCTracklist_all[l]);
733 fhMCtest->Fill(invmass);
734
735 if (motherId1 != -1) {
736 int mcMotherPdg1 = -1;
737 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
738 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
739 if (mcMotherPdg1 == 111) {
740 Double_t invmass2 = 0;
741 invmass2 = Invmass_4particles(fMCTracklist_all[i], fMCTracklist_all[j], fMCTracklist_all[k],
742 fMCTracklist_all[l]);
743 fhMCtest->Fill(invmass2);
744 fhInvariantMass_MC_pi0->Fill(invmass2);
745 fhInvariantMass_MC_pi0_epem->Fill(invmass2);
746 cout << "####################################### Decay pi0 "
747 "-> e+e-e+e- detected!\t\t"
748 << invmass2 << endl;
749 }
750 }
751 else {
752 continue;
753 }
754 }
755
756
757 int grandmotherId1 = -1;
758 int grandmotherId2 = -1;
759 int grandmotherId3 = -1;
760 int grandmotherId4 = -1;
761
762 int mcMotherPdg1 = -1;
763 int mcMotherPdg2 = -1;
764 int mcMotherPdg3 = -1;
765 // int mcMotherPdg4 = -1;
766 int mcGrandmotherPdg1 = -1;
767 // int mcGrandmotherPdg2 = -1;
768 // int mcGrandmotherPdg3 = -1;
769 // int mcGrandmotherPdg4 = -1;
770
771
772 if (motherId1 != -1) {
773 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
774 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
775 grandmotherId1 = mother1->GetMotherId();
776 if (grandmotherId1 != -1) {
777 CbmMCTrack* grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
778 if (NULL != grandmother1) mcGrandmotherPdg1 = grandmother1->GetPdgCode();
779 }
780 }
781 if (motherId2 != -1) {
782 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
783 if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode();
784 grandmotherId2 = mother2->GetMotherId();
785 if (grandmotherId2 != -1) {
786 // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
787 // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
788 }
789 }
790 if (motherId3 != -1) {
791 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
792 if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode();
793 grandmotherId3 = mother3->GetMotherId();
794 if (grandmotherId3 != -1) {
795 // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
796 // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
797 }
798 }
799 if (motherId4 != -1) {
800 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
801 // if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode();
802 grandmotherId4 = mother4->GetMotherId();
803 if (grandmotherId4 != -1) {
804 // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
805 // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
806 }
807 }
808
809
810 if (motherId1 == motherId2 && motherId1 == motherId3) {}
811 if (motherId1 == motherId2 && motherId1 == motherId4) {}
812 if (motherId1 == motherId3 && motherId1 == motherId4) {}
813 if (motherId2 == motherId3 && motherId2 == motherId4) {}
814
815 // if( ((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 22 || mcMotherPdg1 == 111) && (mcMotherPdg3 == 22 || mcMotherPdg3 == 111))
816 // || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 22 || mcMotherPdg1 == 111) && (mcMotherPdg2 == 22 || mcMotherPdg2 == 111))
817 // || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 22 || mcMotherPdg1 == 111) && (mcMotherPdg2 == 22 || mcMotherPdg2 == 111))) {
818
819 // decay X -> gg -> e+ e- e+ e-
820 if (((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg3 == 22)
821 && grandmotherId1 == grandmotherId3)
822 || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 22)
823 && grandmotherId1 == grandmotherId2)
824 || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 22)
825 && grandmotherId1 == grandmotherId2)) {
826 Double_t invmass =
827 Invmass_4particles(fMCTracklist_all[i], fMCTracklist_all[j], fMCTracklist_all[k], fMCTracklist_all[l]);
828 //fhMCtest->Fill(invmass);
829 fhInvariantMass_MC_all->Fill(invmass);
830 if (mcGrandmotherPdg1 == 111) {
831 fhInvariantMass_MC_pi0->Fill(invmass);
832 fhInvariantMass_MC_pi0_gg->Fill(invmass);
833 fhMC_electrons_theta->Fill(momentum1.Theta() * 180. / TMath::Pi());
834 fhMC_electrons_theta->Fill(momentum2.Theta() * 180. / TMath::Pi());
835 fhMC_electrons_theta->Fill(momentum3.Theta() * 180. / TMath::Pi());
836 fhMC_electrons_theta->Fill(momentum4.Theta() * 180. / TMath::Pi());
837 fhMC_electrons_p->Fill(momentum1.Mag());
838 fhMC_electrons_p->Fill(momentum2.Mag());
839 fhMC_electrons_p->Fill(momentum3.Mag());
840 fhMC_electrons_p->Fill(momentum4.Mag());
841 fhMC_electrons_theta_vs_p->Fill(momentum1.Theta() * 180. / TMath::Pi(), momentum1.Mag());
842 fhMC_electrons_theta_vs_p->Fill(momentum2.Theta() * 180. / TMath::Pi(), momentum2.Mag());
843 fhMC_electrons_theta_vs_p->Fill(momentum3.Theta() * 180. / TMath::Pi(), momentum3.Mag());
844 fhMC_electrons_theta_vs_p->Fill(momentum4.Theta() * 180. / TMath::Pi(), momentum4.Mag());
845 }
846 if (mcGrandmotherPdg1 == 221) {
847 fhInvariantMass_MC_eta->Fill(invmass);
848
849 Double_t opening_angle_gg = 0;
850 if (motherId1 == motherId2) {
851 opening_angle_gg = CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma(momentum1, momentum2,
852 momentum3, momentum4);
853 }
854 if (motherId1 == motherId3) {
855 opening_angle_gg = CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma(momentum1, momentum3,
856 momentum2, momentum4);
857 }
858 if (motherId1 == motherId4) {
859 opening_angle_gg = CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma(momentum1, momentum4,
860 momentum2, momentum3);
861 }
862
863 fhEta_openingAngleGG->Fill(opening_angle_gg);
864 }
865 if (TMath::Abs(mcGrandmotherPdg1) == 331) { // eta prime (958)
866 fhInvariantMass_MC_etaPrime->Fill(invmass);
867 }
868 fhMC_grandmotherPDGs->Fill(TMath::Abs(mcGrandmotherPdg1));
869 }
870
871
872 // decay pi0 -> g e+ e- -> e+ e- e+ e-
873 if (((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg3 == 111)
874 && grandmotherId1 == motherId3)
875 || ((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 111) && (mcMotherPdg3 == 22)
876 && grandmotherId3 == motherId1)
877 || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 111)
878 && grandmotherId1 == motherId2)
879 || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 111) && (mcMotherPdg2 == 22)
880 && grandmotherId2 == motherId1)
881 || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 111)
882 && grandmotherId1 == motherId2)
883 || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 111) && (mcMotherPdg2 == 22)
884 && grandmotherId2 == motherId1)) {
885 Double_t invmass =
886 Invmass_4particles(fMCTracklist_all[i], fMCTracklist_all[j], fMCTracklist_all[k], fMCTracklist_all[l]);
887 fhInvariantMass_MC_pi0_gepem->Fill(invmass);
888 fhInvariantMass_MC_pi0->Fill(invmass);
889
890 fhMC_electrons_theta->Fill(momentum1.Theta() * 180. / TMath::Pi());
891 fhMC_electrons_theta->Fill(momentum2.Theta() * 180. / TMath::Pi());
892 fhMC_electrons_theta->Fill(momentum3.Theta() * 180. / TMath::Pi());
893 fhMC_electrons_theta->Fill(momentum4.Theta() * 180. / TMath::Pi());
894 fhMC_electrons_p->Fill(momentum1.Mag());
895 fhMC_electrons_p->Fill(momentum2.Mag());
896 fhMC_electrons_p->Fill(momentum3.Mag());
897 fhMC_electrons_p->Fill(momentum4.Mag());
898 fhMC_electrons_theta_vs_p->Fill(momentum1.Theta() * 180. / TMath::Pi(), momentum1.Mag());
899 fhMC_electrons_theta_vs_p->Fill(momentum2.Theta() * 180. / TMath::Pi(), momentum2.Mag());
900 fhMC_electrons_theta_vs_p->Fill(momentum3.Theta() * 180. / TMath::Pi(), momentum3.Mag());
901 fhMC_electrons_theta_vs_p->Fill(momentum4.Theta() * 180. / TMath::Pi(), momentum4.Mag());
902 }
903 // decay eta -> g e+ e- -> e+ e- e+ e-
904 if (((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg3 == 221)
905 && grandmotherId1 == motherId3)
906 || ((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 221) && (mcMotherPdg3 == 22)
907 && grandmotherId3 == motherId1)
908 || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 221)
909 && grandmotherId1 == motherId2)
910 || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 221) && (mcMotherPdg2 == 22)
911 && grandmotherId2 == motherId1)
912 || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 221)
913 && grandmotherId1 == motherId2)
914 || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 221) && (mcMotherPdg2 == 22)
915 && grandmotherId2 == motherId1)) {
916 Double_t invmass =
917 Invmass_4particles(fMCTracklist_all[i], fMCTracklist_all[j], fMCTracklist_all[k], fMCTracklist_all[l]);
918 fhInvariantMass_MC_eta->Fill(invmass);
919 }
920 }
921 }
922 }
923 }
924 }
925 cout << "CbmAnaConversionReco: InvariantMassTestMC - End!" << endl;
926
927 timer.Stop();
928 fTime += timer.RealTime();
929}
930
931
932Double_t CbmAnaConversionReco::Invmass_4particles(const CbmMCTrack* mctrack1, const CbmMCTrack* mctrack2,
933 const CbmMCTrack* mctrack3, const CbmMCTrack* mctrack4)
934// calculation of invariant mass from four electrons/positrons
935{
936 /* TVector3 mom1;
937 mctrack1->GetMomentum(mom1);
938 TVector3 tempmom1;
939 tempmom1.SetX(SmearValue(mom1.X()));
940 tempmom1.SetY(SmearValue(mom1.Y()));
941 tempmom1.SetZ(SmearValue(mom1.Z()));
942 Double_t energy1 = TMath::Sqrt(tempmom1.Mag2() + M2E);
943 TLorentzVector lorVec1(tempmom1, energy1);
944
945 TVector3 mom2;
946 mctrack2->GetMomentum(mom2);
947 TVector3 tempmom2;
948 tempmom2.SetX(SmearValue(mom2.X()));
949 tempmom2.SetY(SmearValue(mom2.Y()));
950 tempmom2.SetZ(SmearValue(mom2.Z()));
951 Double_t energy2 = TMath::Sqrt(tempmom2.Mag2() + M2E);
952 TLorentzVector lorVec2(tempmom2, energy2);
953
954 TVector3 mom3;
955 mctrack3->GetMomentum(mom3);
956 TVector3 tempmom3;
957 tempmom3.SetX(SmearValue(mom3.X()));
958 tempmom3.SetY(SmearValue(mom3.Y()));
959 tempmom3.SetZ(SmearValue(mom3.Z()));
960 Double_t energy3 = TMath::Sqrt(tempmom3.Mag2() + M2E);
961 TLorentzVector lorVec3(tempmom3, energy3);
962
963 TVector3 mom4;
964 mctrack4->GetMomentum(mom4);
965 TVector3 tempmom4;
966 tempmom4.SetX(SmearValue(mom4.X()));
967 tempmom4.SetY(SmearValue(mom4.Y()));
968 tempmom4.SetZ(SmearValue(mom4.Z()));
969 Double_t energy4 = TMath::Sqrt(tempmom4.Mag2() + M2E);
970 TLorentzVector lorVec4(tempmom4, energy4);
971*/
972 TLorentzVector lorVec1;
973 mctrack1->Get4Momentum(lorVec1);
974
975 TLorentzVector lorVec2;
976 mctrack2->Get4Momentum(lorVec2);
977
978 TLorentzVector lorVec3;
979 mctrack3->Get4Momentum(lorVec3);
980
981 TLorentzVector lorVec4;
982 mctrack4->Get4Momentum(lorVec4);
983
984
985 TLorentzVector sum;
986 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
987 //cout << "mc: \t" << sum.Px() << " / " << sum.Py() << " / " << sum.Pz() << " / " << sum.E() << "\t => mag = " << sum.Mag() << endl;
988
989
990 return sum.Mag();
991}
992
993
994Double_t CbmAnaConversionReco::SmearValue(Double_t value)
995{
996 TRandom3 generator(0);
997 Double_t result = 0;
998 // Double_t smear = 0;
999 Int_t plusminus = 0;
1000 while (plusminus == 0) { // should be either 1 or -1, not 0
1001 plusminus = generator.Uniform(-2, 2);
1002 }
1003 // Double_t gaus = generator.Gaus(1,1);
1004 // smear = gaus * plusminus;
1005 // result = value * (1. + 1.0*smear/100); //smearing as wished
1006
1007 result = value; // -> no smearing
1008
1009 return result;
1010}
1011
1012
1013Double_t CbmAnaConversionReco::Invmass_4particlesRECO(const TVector3 part1, const TVector3 part2, const TVector3 part3,
1014 const TVector3 part4)
1015// calculation of invariant mass from four electrons/positrons
1016{
1017 Double_t energy1 = TMath::Sqrt(part1.Mag2() + M2E);
1018 TLorentzVector lorVec1(part1, energy1);
1019
1020 Double_t energy2 = TMath::Sqrt(part2.Mag2() + M2E);
1021 TLorentzVector lorVec2(part2, energy2);
1022
1023 Double_t energy3 = TMath::Sqrt(part3.Mag2() + M2E);
1024 TLorentzVector lorVec3(part3, energy3);
1025
1026 Double_t energy4 = TMath::Sqrt(part4.Mag2() + M2E);
1027 TLorentzVector lorVec4(part4, energy4);
1028
1029 TLorentzVector sum;
1030 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
1031 //cout << "reco: \t" << sum.Px() << " / " << sum.Py() << " / " << sum.Pz() << " / " << sum.E() << "\t => mag = " << sum.Mag() << endl;
1032
1033
1034 return sum.Mag();
1035}
1036
1037
1038void CbmAnaConversionReco::InvariantMassTest_4epem()
1039// Calculating invariant mass of 4 ep/em, using MC data AND reconstructed momentum
1040{
1041 timer.Start();
1042
1043
1044 cout << "CbmAnaConversionReco: InvariantMassTest_4epem - Start..." << endl;
1045 cout << "CbmAnaConversionReco: InvariantMassTest_4epem - " << fRecoTracklistEPEM.size() << "\t"
1046 << fRecoMomentum.size() << endl;
1047 int fill = 0;
1048 if (fRecoTracklistEPEM.size() < 4) return;
1049 for (unsigned int i = 0; i < fRecoTracklistEPEM.size(); i++) {
1050 if (i % 10 == 0) cout << "CbmAnaConversionReco: InvariantMassTest_4epem - iteration i = " << i << endl;
1051 for (unsigned int j = i + 1; j < fRecoTracklistEPEM.size(); j++) {
1052 for (unsigned int k = j + 1; k < fRecoTracklistEPEM.size(); k++) {
1053 for (unsigned int l = k + 1; l < fRecoTracklistEPEM.size(); l++) {
1054 if (fRecoTracklistEPEM[i]->GetPdgCode() + fRecoTracklistEPEM[j]->GetPdgCode()
1055 + fRecoTracklistEPEM[k]->GetPdgCode() + fRecoTracklistEPEM[l]->GetPdgCode()
1056 != 0)
1057 continue;
1058
1059 if (fRecoTracklistEPEM.size() != fRecoMomentum.size()
1060 || fRecoTracklistEPEM.size() != fRecoRefittedMomentum.size()) {
1061 cout << "CbmAnaConversionReco: InvariantMassTest_4epem - not "
1062 "matching number of entries!"
1063 << endl;
1064 continue;
1065 }
1066
1067
1068 // starting points of each electron (-> i.e. conversion points of gamma OR decay points of pi0, depending on decay channel)
1069 TVector3 pi0start_i;
1070 fRecoTracklistEPEM[i]->GetStartVertex(pi0start_i);
1071 TVector3 pi0start_j;
1072 fRecoTracklistEPEM[j]->GetStartVertex(pi0start_j);
1073 TVector3 pi0start_k;
1074 fRecoTracklistEPEM[k]->GetStartVertex(pi0start_k);
1075 TVector3 pi0start_l;
1076 fRecoTracklistEPEM[l]->GetStartVertex(pi0start_l);
1077
1078
1079 int motherId1 = fRecoTracklistEPEM[i]->GetMotherId();
1080 int motherId2 = fRecoTracklistEPEM[j]->GetMotherId();
1081 int motherId3 = fRecoTracklistEPEM[k]->GetMotherId();
1082 int motherId4 = fRecoTracklistEPEM[l]->GetMotherId();
1083
1084 // decay pi0 -> e+ e- e+ e-
1085 if (motherId1 == motherId2 && motherId1 == motherId3 && motherId1 == motherId4) {
1086 if (motherId1 != -1) {
1087 int mcMotherPdg1 = -1;
1088 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
1089 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
1090 if ((mcMotherPdg1 == 111 || mcMotherPdg1 == 221)) { // && NofDaughters(motherId1) == 4) {
1091 Double_t invmass1 = 0;
1092 // Double_t invmass2 = 0; // momenta from stsMomentumVec
1093 Double_t invmass3 = 0;
1094 // invmass2 = Invmass_4particlesRECO(fRecoMomentum[i], fRecoMomentum[j], fRecoMomentum[k], fRecoMomentum[l]);
1095 invmass1 = Invmass_4particles(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j], fRecoTracklistEPEM[k],
1096 fRecoTracklistEPEM[l]); // true MC values
1097 invmass3 = Invmass_4particlesRECO(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1098 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1099 // fhInvariantMass_pi0epem->Fill(invmass2);
1100 cout << "Decay pi0 -> e+e-e+e- detected!\t\t mc mass: " << invmass1 << "\t, reco mass: " << invmass3
1101 << endl;
1102 cout << "motherids: " << motherId1 << "/" << motherId2 << "/" << motherId3 << "/" << motherId4
1103 << "\t motherpdg: " << mcMotherPdg1 << "\t mctrack mass: " << mother1->GetMass()
1104 << "\t nofdaughters: " << NofDaughters(motherId1) << endl;
1105 cout << "pdgs " << fRecoTracklistEPEM[i]->GetPdgCode() << "/" << fRecoTracklistEPEM[j]->GetPdgCode()
1106 << "/" << fRecoTracklistEPEM[k]->GetPdgCode() << "/" << fRecoTracklistEPEM[l]->GetPdgCode()
1107 << endl;
1108 TVector3 start1;
1109 fRecoTracklistEPEM[i]->GetStartVertex(start1);
1110 TVector3 start2;
1111 fRecoTracklistEPEM[j]->GetStartVertex(start2);
1112 TVector3 start3;
1113 fRecoTracklistEPEM[k]->GetStartVertex(start3);
1114 TVector3 start4;
1115 fRecoTracklistEPEM[l]->GetStartVertex(start4);
1116 cout << "start: " << start1.Z() << "/" << start2.Z() << "/" << start3.Z() << "/" << start4.Z() << endl;
1117
1118 fhEPEM_invmass_eeee_mc->Fill(invmass1);
1119 fhEPEM_invmass_eeee_refitted->Fill(invmass3);
1120 //fhEPEM_invmass_all_mc->Fill(invmass1);
1121 //fhEPEM_invmass_all_refitted->Fill(invmass3);
1122 }
1123 }
1124 else { // all 4 particles come directly from the vertex
1125 continue;
1126 }
1127 }
1128
1129 if ((motherId1 == motherId2 && motherId3 == motherId4) || (motherId1 == motherId3 && motherId2 == motherId4)
1130 || (motherId1 == motherId4 && motherId2 == motherId3)) {
1131
1132
1133 int grandmotherId1 = -1;
1134 int grandmotherId2 = -1;
1135 int grandmotherId3 = -1;
1136 int grandmotherId4 = -1;
1137
1138 int mcMotherPdg1 = -1;
1139 int mcMotherPdg2 = -1;
1140 int mcMotherPdg3 = -1;
1141 int mcMotherPdg4 = -1;
1142 int mcGrandmotherPdg1 = -1;
1143 // int mcGrandmotherPdg2 = -1;
1144 // int mcGrandmotherPdg3 = -1;
1145 // int mcGrandmotherPdg4 = -1;
1146
1147 CbmMCTrack* grandmother1 = NULL;
1148
1149 if (motherId1 != -1) {
1150 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
1151 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
1152 grandmotherId1 = mother1->GetMotherId();
1153 if (grandmotherId1 != -1) {
1154 grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
1155 if (NULL != grandmother1) mcGrandmotherPdg1 = grandmother1->GetPdgCode();
1156 }
1157 }
1158 if (motherId2 != -1) {
1159 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
1160 if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode();
1161 grandmotherId2 = mother2->GetMotherId();
1162 if (grandmotherId2 != -1) {
1163 // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
1164 // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
1165 }
1166 }
1167 if (motherId3 != -1) {
1168 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
1169 if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode();
1170 grandmotherId3 = mother3->GetMotherId();
1171 if (grandmotherId3 != -1) {
1172 // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
1173 // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
1174 }
1175 }
1176 if (motherId4 != -1) {
1177 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
1178 if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode();
1179 grandmotherId4 = mother4->GetMotherId();
1180 if (grandmotherId4 != -1) {
1181 // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
1182 // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
1183 }
1184 }
1185
1186 //if(grandmotherId1 == -1) continue;
1187
1188 if (motherId1 == motherId2 && motherId3 == motherId4) {
1189 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId3) != 2) continue;
1190 if ((grandmotherId1 == motherId3 && mcMotherPdg3 == 111)
1191 || (motherId1 == grandmotherId3 && mcMotherPdg1 == 111)) {
1192
1193 fhPi0_startvertexElectrons_gee->Fill(pi0start_i.Z());
1194 fhPi0_startvertexElectrons_gee->Fill(pi0start_j.Z());
1195 fhPi0_startvertexElectrons_gee->Fill(pi0start_k.Z());
1196 fhPi0_startvertexElectrons_gee->Fill(pi0start_l.Z());
1197 fhPi0_startvertexElectrons_all->Fill(pi0start_i.Z());
1198 fhPi0_startvertexElectrons_all->Fill(pi0start_j.Z());
1199 fhPi0_startvertexElectrons_all->Fill(pi0start_k.Z());
1200 fhPi0_startvertexElectrons_all->Fill(pi0start_l.Z());
1201
1202 Double_t invmass1 = 0;
1203 Double_t invmass3 = 0;
1204 invmass1 = Invmass_4particles(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j], fRecoTracklistEPEM[k],
1205 fRecoTracklistEPEM[l]); // true MC values
1206 invmass3 = Invmass_4particlesRECO(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1207 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1208
1209 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1210 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) {
1211 fhEPEM_InDetector_invmass_gee_mc->Fill(invmass1);
1212 fhEPEM_InDetector_invmass_gee_refitted->Fill(invmass3);
1213 fhEPEM_InDetector_invmass_all_mc->Fill(invmass1);
1214 fhEPEM_InDetector_invmass_all_refitted->Fill(invmass3);
1215 continue;
1216 }
1217
1218
1219 cout << "HURRAY! .-.-.-.-.-.-..-.-.-.-.-.-.-.-.-.-.-.-." << endl;
1220 CutEfficiencyStudies(i, j, k, l, motherId1, motherId2, motherId3, motherId4);
1221 fhInvariantMass_pi0epem->Fill(invmass1);
1222
1223 fhEPEM_invmass_gee_mc->Fill(invmass1);
1224 fhEPEM_invmass_gee_refitted->Fill(invmass3);
1225 fhEPEM_invmass_all_mc->Fill(invmass1);
1226 fhEPEM_invmass_all_refitted->Fill(invmass3);
1227
1228 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[i]->GetPt(), fRecoTracklistEPEM[i]->GetP());
1229 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[j]->GetPt(), fRecoTracklistEPEM[j]->GetP());
1230 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[k]->GetPt(), fRecoTracklistEPEM[k]->GetP());
1231 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[l]->GetPt(), fRecoTracklistEPEM[l]->GetP());
1232 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[i].Perp(), fRecoRefittedMomentum[i].Mag());
1233 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[j].Perp(), fRecoRefittedMomentum[j].Mag());
1234 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[k].Perp(), fRecoRefittedMomentum[k].Mag());
1235 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[l].Perp(), fRecoRefittedMomentum[l].Mag());
1236
1237 LmvmKinePar params1 =
1238 CalculateKinematicParams_4particles(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1239 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1240 fhPi0_pt_vs_rap_gee->Fill(params1.fPt, params1.fRapidity);
1241 fhPi0_pt_vs_rap_all->Fill(params1.fPt, params1.fRapidity);
1242 fhPi0_pt_gee->Fill(params1.fPt);
1243 fhPi0_pt_all->Fill(params1.fPt);
1244
1245 if (mcGrandmotherPdg1 == 111) { // case: i,j = electrons from gamma, k,l = electrons from pi0
1246 Double_t opening_angle1_mc = 0;
1247 Double_t opening_angle1_refitted = 0;
1248 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j]);
1249 opening_angle1_refitted =
1250 CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j]);
1251 fhEPEM_openingAngle_gee_mc->Fill(opening_angle1_mc);
1252 fhEPEM_openingAngle_gee_refitted->Fill(opening_angle1_refitted);
1253
1254 Double_t opening_angle1_mc_dalitz = 0;
1255 Double_t opening_angle1_refitted_dalitz = 0;
1256 opening_angle1_mc_dalitz = CalculateOpeningAngleMC(fRecoTracklistEPEM[k], fRecoTracklistEPEM[l]);
1257 opening_angle1_refitted_dalitz =
1258 CalculateOpeningAngleReco(fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1259 fhEPEM_openingAngle_gee_mc_dalitz->Fill(opening_angle1_mc_dalitz);
1260 fhEPEM_openingAngle_gee_refitted_dalitz->Fill(opening_angle1_refitted_dalitz);
1261 }
1262
1263 if (mcMotherPdg1 == 111) { // case: i,j = electrons from pi0, k,l = electrons from gamma
1264 Double_t opening_angle1_mc = 0;
1265 Double_t opening_angle1_refitted = 0;
1266 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[k], fRecoTracklistEPEM[l]);
1267 opening_angle1_refitted =
1268 CalculateOpeningAngleReco(fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1269 fhEPEM_openingAngle_gee_mc->Fill(opening_angle1_mc);
1270 fhEPEM_openingAngle_gee_refitted->Fill(opening_angle1_refitted);
1271
1272 Double_t opening_angle1_mc_dalitz = 0;
1273 Double_t opening_angle1_refitted_dalitz = 0;
1274 opening_angle1_mc_dalitz = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j]);
1275 opening_angle1_refitted_dalitz =
1276 CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j]);
1277 fhEPEM_openingAngle_gee_mc_dalitz->Fill(opening_angle1_mc_dalitz);
1278 fhEPEM_openingAngle_gee_refitted_dalitz->Fill(opening_angle1_refitted_dalitz);
1279 }
1280 }
1281 }
1282 if (motherId1 == motherId3 && motherId2 == motherId4) {
1283 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2) continue;
1284 if ((grandmotherId1 == motherId2 && mcMotherPdg2 == 111)
1285 || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
1286
1287 fhPi0_startvertexElectrons_gee->Fill(pi0start_i.Z());
1288 fhPi0_startvertexElectrons_gee->Fill(pi0start_j.Z());
1289 fhPi0_startvertexElectrons_gee->Fill(pi0start_k.Z());
1290 fhPi0_startvertexElectrons_gee->Fill(pi0start_l.Z());
1291 fhPi0_startvertexElectrons_all->Fill(pi0start_i.Z());
1292 fhPi0_startvertexElectrons_all->Fill(pi0start_j.Z());
1293 fhPi0_startvertexElectrons_all->Fill(pi0start_k.Z());
1294 fhPi0_startvertexElectrons_all->Fill(pi0start_l.Z());
1295
1296 Double_t invmass1 = 0;
1297 Double_t invmass3 = 0;
1298 invmass1 = Invmass_4particles(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j], fRecoTracklistEPEM[k],
1299 fRecoTracklistEPEM[l]); // true MC values
1300 invmass3 = Invmass_4particlesRECO(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1301 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1302
1303 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1304 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) {
1305 fhEPEM_InDetector_invmass_gee_mc->Fill(invmass1);
1306 fhEPEM_InDetector_invmass_gee_refitted->Fill(invmass3);
1307 fhEPEM_InDetector_invmass_all_mc->Fill(invmass1);
1308 fhEPEM_InDetector_invmass_all_refitted->Fill(invmass3);
1309 continue;
1310 }
1311
1312
1313 cout << "HURRAY! .-.-.-.-.-.-..-.-.-.-.-.-.-.-.-.-.-.-." << endl;
1314 CutEfficiencyStudies(i, j, k, l, motherId1, motherId2, motherId3, motherId4);
1315 fhInvariantMass_pi0epem->Fill(invmass1);
1316
1317 fhEPEM_invmass_gee_mc->Fill(invmass1);
1318 fhEPEM_invmass_gee_refitted->Fill(invmass3);
1319 fhEPEM_invmass_all_mc->Fill(invmass1);
1320 fhEPEM_invmass_all_refitted->Fill(invmass3);
1321
1322 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[i]->GetPt(), fRecoTracklistEPEM[i]->GetP());
1323 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[j]->GetPt(), fRecoTracklistEPEM[j]->GetP());
1324 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[k]->GetPt(), fRecoTracklistEPEM[k]->GetP());
1325 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[l]->GetPt(), fRecoTracklistEPEM[l]->GetP());
1326 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[i].Perp(), fRecoRefittedMomentum[i].Mag());
1327 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[j].Perp(), fRecoRefittedMomentum[j].Mag());
1328 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[k].Perp(), fRecoRefittedMomentum[k].Mag());
1329 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[l].Perp(), fRecoRefittedMomentum[l].Mag());
1330
1331 LmvmKinePar params1 =
1332 CalculateKinematicParams_4particles(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1333 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1334 fhPi0_pt_vs_rap_gee->Fill(params1.fPt, params1.fRapidity);
1335 fhPi0_pt_vs_rap_all->Fill(params1.fPt, params1.fRapidity);
1336 fhPi0_pt_gee->Fill(params1.fPt);
1337 fhPi0_pt_all->Fill(params1.fPt);
1338
1339 if (mcGrandmotherPdg1 == 111) {
1340 Double_t opening_angle1_mc = 0;
1341 Double_t opening_angle1_refitted = 0;
1342 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[k]);
1343 opening_angle1_refitted =
1344 CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[k]);
1345 fhEPEM_openingAngle_gee_mc->Fill(opening_angle1_mc);
1346 fhEPEM_openingAngle_gee_refitted->Fill(opening_angle1_refitted);
1347
1348 Double_t opening_angle1_mc_dalitz = 0;
1349 Double_t opening_angle1_refitted_dalitz = 0;
1350 opening_angle1_mc_dalitz = CalculateOpeningAngleMC(fRecoTracklistEPEM[j], fRecoTracklistEPEM[l]);
1351 opening_angle1_refitted_dalitz =
1352 CalculateOpeningAngleReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[l]);
1353 fhEPEM_openingAngle_gee_mc_dalitz->Fill(opening_angle1_mc_dalitz);
1354 fhEPEM_openingAngle_gee_refitted_dalitz->Fill(opening_angle1_refitted_dalitz);
1355 }
1356
1357 if (mcMotherPdg1 == 111) {
1358 Double_t opening_angle1_mc = 0;
1359 Double_t opening_angle1_refitted = 0;
1360 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[j], fRecoTracklistEPEM[l]);
1361 opening_angle1_refitted =
1362 CalculateOpeningAngleReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[l]);
1363 fhEPEM_openingAngle_gee_mc->Fill(opening_angle1_mc);
1364 fhEPEM_openingAngle_gee_refitted->Fill(opening_angle1_refitted);
1365
1366 Double_t opening_angle1_mc_dalitz = 0;
1367 Double_t opening_angle1_refitted_dalitz = 0;
1368 opening_angle1_mc_dalitz = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[k]);
1369 opening_angle1_refitted_dalitz =
1370 CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[k]);
1371 fhEPEM_openingAngle_gee_mc_dalitz->Fill(opening_angle1_mc_dalitz);
1372 fhEPEM_openingAngle_gee_refitted_dalitz->Fill(opening_angle1_refitted_dalitz);
1373 }
1374 }
1375 }
1376 if (motherId1 == motherId4 && motherId2 == motherId3) {
1377 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2) continue;
1378 if ((grandmotherId1 == motherId2 && mcMotherPdg2 == 111)
1379 || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
1380
1381 fhPi0_startvertexElectrons_gee->Fill(pi0start_i.Z());
1382 fhPi0_startvertexElectrons_gee->Fill(pi0start_j.Z());
1383 fhPi0_startvertexElectrons_gee->Fill(pi0start_k.Z());
1384 fhPi0_startvertexElectrons_gee->Fill(pi0start_l.Z());
1385 fhPi0_startvertexElectrons_all->Fill(pi0start_i.Z());
1386 fhPi0_startvertexElectrons_all->Fill(pi0start_j.Z());
1387 fhPi0_startvertexElectrons_all->Fill(pi0start_k.Z());
1388 fhPi0_startvertexElectrons_all->Fill(pi0start_l.Z());
1389
1390 Double_t invmass1 = 0;
1391 Double_t invmass3 = 0;
1392 invmass1 = Invmass_4particles(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j], fRecoTracklistEPEM[k],
1393 fRecoTracklistEPEM[l]); // true MC values
1394 invmass3 = Invmass_4particlesRECO(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1395 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1396
1397 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1398 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) {
1399 fhEPEM_InDetector_invmass_gee_mc->Fill(invmass1);
1400 fhEPEM_InDetector_invmass_gee_refitted->Fill(invmass3);
1401 fhEPEM_InDetector_invmass_all_mc->Fill(invmass1);
1402 fhEPEM_InDetector_invmass_all_refitted->Fill(invmass3);
1403 continue;
1404 }
1405
1406
1407 cout << "HURRAY! .-.-.-.-.-.-..-.-.-.-.-.-.-.-.-.-.-.-." << endl;
1408 CutEfficiencyStudies(i, j, k, l, motherId1, motherId2, motherId3, motherId4);
1409 fhInvariantMass_pi0epem->Fill(invmass1);
1410
1411 fhEPEM_invmass_gee_mc->Fill(invmass1);
1412 fhEPEM_invmass_gee_refitted->Fill(invmass3);
1413 fhEPEM_invmass_all_mc->Fill(invmass1);
1414 fhEPEM_invmass_all_refitted->Fill(invmass3);
1415
1416 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[i]->GetPt(), fRecoTracklistEPEM[i]->GetP());
1417 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[j]->GetPt(), fRecoTracklistEPEM[j]->GetP());
1418 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[k]->GetPt(), fRecoTracklistEPEM[k]->GetP());
1419 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[l]->GetPt(), fRecoTracklistEPEM[l]->GetP());
1420 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[i].Perp(), fRecoRefittedMomentum[i].Mag());
1421 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[j].Perp(), fRecoRefittedMomentum[j].Mag());
1422 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[k].Perp(), fRecoRefittedMomentum[k].Mag());
1423 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[l].Perp(), fRecoRefittedMomentum[l].Mag());
1424
1425 LmvmKinePar params1 =
1426 CalculateKinematicParams_4particles(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1427 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1428 fhPi0_pt_vs_rap_gee->Fill(params1.fPt, params1.fRapidity);
1429 fhPi0_pt_vs_rap_all->Fill(params1.fPt, params1.fRapidity);
1430 fhPi0_pt_gee->Fill(params1.fPt);
1431 fhPi0_pt_all->Fill(params1.fPt);
1432
1433 if (mcGrandmotherPdg1 == 111) {
1434 Double_t opening_angle1_mc = 0;
1435 Double_t opening_angle1_refitted = 0;
1436 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[l]);
1437 opening_angle1_refitted =
1438 CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[l]);
1439 fhEPEM_openingAngle_gee_mc->Fill(opening_angle1_mc);
1440 fhEPEM_openingAngle_gee_refitted->Fill(opening_angle1_refitted);
1441
1442 Double_t opening_angle1_mc_dalitz = 0;
1443 Double_t opening_angle1_refitted_dalitz = 0;
1444 opening_angle1_mc_dalitz = CalculateOpeningAngleMC(fRecoTracklistEPEM[j], fRecoTracklistEPEM[k]);
1445 opening_angle1_refitted_dalitz =
1446 CalculateOpeningAngleReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[k]);
1447 fhEPEM_openingAngle_gee_mc_dalitz->Fill(opening_angle1_mc_dalitz);
1448 fhEPEM_openingAngle_gee_refitted_dalitz->Fill(opening_angle1_refitted_dalitz);
1449 }
1450
1451 if (mcMotherPdg1 == 111) {
1452 Double_t opening_angle1_mc = 0;
1453 Double_t opening_angle1_refitted = 0;
1454 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[j], fRecoTracklistEPEM[k]);
1455 opening_angle1_refitted =
1456 CalculateOpeningAngleReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[k]);
1457 fhEPEM_openingAngle_gee_mc->Fill(opening_angle1_mc);
1458 fhEPEM_openingAngle_gee_refitted->Fill(opening_angle1_refitted);
1459
1460 Double_t opening_angle1_mc_dalitz = 0;
1461 Double_t opening_angle1_refitted_dalitz = 0;
1462 opening_angle1_mc_dalitz = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[l]);
1463 opening_angle1_refitted_dalitz =
1464 CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[l]);
1465 fhEPEM_openingAngle_gee_mc_dalitz->Fill(opening_angle1_mc_dalitz);
1466 fhEPEM_openingAngle_gee_refitted_dalitz->Fill(opening_angle1_refitted_dalitz);
1467 }
1468 }
1469 }
1470
1471
1472 // ===================================================================================================
1473 // HERE DECAY pi0 -> gamma gamma -> e+e- e+e-
1474 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
1475 && grandmotherId1 == grandmotherId4) {
1476 // if(mcGrandmotherPdg1 != 111 && mcGrandmotherPdg1 != 221) continue; // 111 = pi0, 221 = eta
1477
1478
1479 if (pi0start_i.Z() <= 1 && pi0start_j.Z() <= 1 && pi0start_k.Z() <= 1 && pi0start_l.Z() <= 1
1480 && mcGrandmotherPdg1 == 221) {
1481 Double_t invmass_eta_mc = Invmass_4particles(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j],
1482 fRecoTracklistEPEM[k], fRecoTracklistEPEM[l]);
1483 Double_t invmass_eta_reco = Invmass_4particlesRECO(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1484 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1485
1486 fhEPEM_invmass_eta_mc->Fill(invmass_eta_mc);
1487 fhEPEM_invmass_eta_refitted->Fill(invmass_eta_reco);
1488
1489 CutEfficiencyStudies(i, j, k, l, motherId1, motherId2, motherId3, motherId4, 1);
1490 }
1491
1492
1493 if (mcGrandmotherPdg1 != 111) continue; // 111 = pi0, 221 = eta
1494
1495 TVector3 pi0start;
1496 grandmother1->GetStartVertex(pi0start);
1497 fhPi0_startvertex->Fill(pi0start.Z());
1498
1499 fhPi0_startvertexElectrons_gg->Fill(pi0start_i.Z());
1500 fhPi0_startvertexElectrons_gg->Fill(pi0start_j.Z());
1501 fhPi0_startvertexElectrons_gg->Fill(pi0start_k.Z());
1502 fhPi0_startvertexElectrons_gg->Fill(pi0start_l.Z());
1503 fhPi0_startvertexElectrons_all->Fill(pi0start_i.Z());
1504 fhPi0_startvertexElectrons_all->Fill(pi0start_j.Z());
1505 fhPi0_startvertexElectrons_all->Fill(pi0start_k.Z());
1506 fhPi0_startvertexElectrons_all->Fill(pi0start_l.Z());
1507
1508 fhPi0_startvertex_vs_chi->Fill(pi0start_i.Z(), fRecoTracklistEPEM_chi[i]);
1509 fhPi0_startvertex_vs_chi->Fill(pi0start_j.Z(), fRecoTracklistEPEM_chi[j]);
1510 fhPi0_startvertex_vs_chi->Fill(pi0start_k.Z(), fRecoTracklistEPEM_chi[l]);
1511 fhPi0_startvertex_vs_chi->Fill(pi0start_l.Z(), fRecoTracklistEPEM_chi[l]);
1512
1513 fhPi0_startvertex_vs_momentum->Fill(pi0start_i.Z(), fRecoTracklistEPEM[i]->GetP());
1514 fhPi0_startvertex_vs_momentum->Fill(pi0start_j.Z(), fRecoTracklistEPEM[j]->GetP());
1515 fhPi0_startvertex_vs_momentum->Fill(pi0start_k.Z(), fRecoTracklistEPEM[k]->GetP());
1516 fhPi0_startvertex_vs_momentum->Fill(pi0start_l.Z(), fRecoTracklistEPEM[l]->GetP());
1517
1518 Double_t invmass1 = 0; // true MC values
1519 invmass1 = Invmass_4particles(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j], fRecoTracklistEPEM[k],
1520 fRecoTracklistEPEM[l]);
1521 Double_t invmass2 = 0; // momenta from stsMomentumVec
1522 invmass2 = Invmass_4particlesRECO(fRecoMomentum[i], fRecoMomentum[j], fRecoMomentum[k], fRecoMomentum[l]);
1523 Double_t invmass3 = 0; // momenta from refitted at primary vertex
1524 invmass3 = Invmass_4particlesRECO(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j],
1525 fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1526
1527 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1528 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) {
1529 fhEPEM_InDetector_invmass_gg_mc->Fill(invmass1);
1530 fhEPEM_InDetector_invmass_gg_refitted->Fill(invmass3);
1531 fhEPEM_InDetector_invmass_all_mc->Fill(invmass1);
1532 fhEPEM_InDetector_invmass_all_refitted->Fill(invmass3);
1533 continue;
1534 }
1535
1536 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[i]->GetPt(), fRecoTracklistEPEM[i]->GetP());
1537 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[j]->GetPt(), fRecoTracklistEPEM[j]->GetP());
1538 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[k]->GetPt(), fRecoTracklistEPEM[k]->GetP());
1539 fhEPEM_pt_vs_p_all_mc->Fill(fRecoTracklistEPEM[l]->GetPt(), fRecoTracklistEPEM[l]->GetP());
1540 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[i].Perp(), fRecoRefittedMomentum[i].Mag());
1541 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[j].Perp(), fRecoRefittedMomentum[j].Mag());
1542 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[k].Perp(), fRecoRefittedMomentum[k].Mag());
1543 fhEPEM_pt_vs_p_all_refitted->Fill(fRecoRefittedMomentum[l].Perp(), fRecoRefittedMomentum[l].Mag());
1544
1545 cout << "########################################################"
1546 "##############"
1547 << endl;
1548 cout << fRecoMomentum[i].X() << "\t" << fRecoRefittedMomentum[i].X() << endl;
1549 cout << fRecoMomentum[i].Y() << "\t" << fRecoRefittedMomentum[i].Y() << endl;
1550 cout << fRecoMomentum[i].Z() << "\t" << fRecoRefittedMomentum[i].Z() << endl;
1551
1552 if (fRecoTracklistEPEM_ids.size() != fRecoTracklistEPEM.size()) {
1553 cout << "size mismatch! " << fRecoTracklistEPEM_ids.size() << "/" << fRecoTracklistEPEM.size() << endl;
1554 }
1555
1556 cout << "########################################################"
1557 "##############"
1558 << endl;
1559 cout << "index: " << i << "\t" << j << "\t" << k << "\t" << l << endl;
1560 cout << "mc id: " << fRecoTracklistEPEM_ids[i] << "\t" << fRecoTracklistEPEM_ids[j] << "\t"
1561 << fRecoTracklistEPEM_ids[k] << "\t" << fRecoTracklistEPEM_ids[l] << endl;
1562 cout << "motherid: " << motherId1 << "\t" << motherId2 << "\t" << motherId3 << "\t" << motherId4 << endl;
1563 cout << "motherpdg: " << mcMotherPdg1 << "\t" << mcMotherPdg2 << "\t" << mcMotherPdg3 << "\t"
1564 << mcMotherPdg4 << endl;
1565 cout << "grandmotherid: " << grandmotherId1 << "\t" << grandmotherId2 << "\t" << grandmotherId3 << "\t"
1566 << grandmotherId4 << endl;
1567 cout << "pdg: " << fRecoTracklistEPEM[i]->GetPdgCode() << "\t" << fRecoTracklistEPEM[j]->GetPdgCode()
1568 << "\t" << fRecoTracklistEPEM[k]->GetPdgCode() << "\t" << fRecoTracklistEPEM[l]->GetPdgCode()
1569 << endl;
1570 cout << "invmass reco: " << invmass2 << "\t invmass mc: " << invmass1 << endl;
1571
1572 fhInvMass_EPEM_mc->Fill(invmass1);
1573 fhInvMass_EPEM_stsMomVec->Fill(invmass2);
1574 fhInvMass_EPEM_refitted->Fill(invmass3);
1575 fhInvMass_EPEM_error_stsMomVec->Fill(1.0 * TMath::Abs(invmass1 - invmass2) / invmass1);
1576 fhInvMass_EPEM_error_refitted->Fill(1.0 * TMath::Abs(invmass1 - invmass3) / invmass1);
1577
1578 fhUsedMomenta_stsMomVec->Fill(fRecoMomentum[i].Mag());
1579 fhUsedMomenta_stsMomVec->Fill(fRecoMomentum[j].Mag());
1580 fhUsedMomenta_stsMomVec->Fill(fRecoMomentum[k].Mag());
1581 fhUsedMomenta_stsMomVec->Fill(fRecoMomentum[l].Mag());
1582
1583 fhUsedMomenta_mc->Fill(fRecoTracklistEPEM[i]->GetP());
1584 fhUsedMomenta_mc->Fill(fRecoTracklistEPEM[j]->GetP());
1585 fhUsedMomenta_mc->Fill(fRecoTracklistEPEM[k]->GetP());
1586 fhUsedMomenta_mc->Fill(fRecoTracklistEPEM[l]->GetP());
1587
1588 fhUsedMomenta_error_stsMomVec->Fill(TMath::Abs(fRecoTracklistEPEM[i]->GetP() - fRecoMomentum[i].Mag())
1589 / fRecoTracklistEPEM[i]->GetP());
1590 fhUsedMomenta_error_stsMomVec->Fill(TMath::Abs(fRecoTracklistEPEM[j]->GetP() - fRecoMomentum[j].Mag())
1591 / fRecoTracklistEPEM[j]->GetP());
1592 fhUsedMomenta_error_stsMomVec->Fill(TMath::Abs(fRecoTracklistEPEM[k]->GetP() - fRecoMomentum[k].Mag())
1593 / fRecoTracklistEPEM[k]->GetP());
1594 fhUsedMomenta_error_stsMomVec->Fill(TMath::Abs(fRecoTracklistEPEM[l]->GetP() - fRecoMomentum[l].Mag())
1595 / fRecoTracklistEPEM[l]->GetP());
1596
1597 // fhUsedMomenta_error_refitted->Fill(TMath::Abs(fRecoRefittedMomentum[i].Mag() - fRecoMomentum[i].Mag())/fRecoRefittedMomentum[i].Mag());
1598 // fhUsedMomenta_error_refitted->Fill(TMath::Abs(fRecoRefittedMomentum[j].Mag() - fRecoMomentum[j].Mag())/fRecoRefittedMomentum[j].Mag());
1599 // fhUsedMomenta_error_refitted->Fill(TMath::Abs(fRecoRefittedMomentum[k].Mag() - fRecoMomentum[k].Mag())/fRecoRefittedMomentum[k].Mag());
1600 // fhUsedMomenta_error_refitted->Fill(TMath::Abs(fRecoRefittedMomentum[l].Mag() - fRecoMomentum[l].Mag())/fRecoRefittedMomentum[l].Mag());
1601
1602 fhUsedMomenta_error_refitted->Fill(
1603 TMath::Abs(fRecoTracklistEPEM[i]->GetP() - fRecoRefittedMomentum[i].Mag())
1604 / fRecoTracklistEPEM[i]->GetP());
1605 fhUsedMomenta_error_refitted->Fill(
1606 TMath::Abs(fRecoTracklistEPEM[j]->GetP() - fRecoRefittedMomentum[j].Mag())
1607 / fRecoTracklistEPEM[j]->GetP());
1608 fhUsedMomenta_error_refitted->Fill(
1609 TMath::Abs(fRecoTracklistEPEM[k]->GetP() - fRecoRefittedMomentum[k].Mag())
1610 / fRecoTracklistEPEM[k]->GetP());
1611 fhUsedMomenta_error_refitted->Fill(
1612 TMath::Abs(fRecoTracklistEPEM[l]->GetP() - fRecoRefittedMomentum[l].Mag())
1613 / fRecoTracklistEPEM[l]->GetP());
1614
1615
1616 TVector3 momentumtest5a;
1617 fRecoTracklistEPEM[i]->GetMomentum(momentumtest5a);
1618 fhUsedMomenta_errorX_stsMomVec->Fill(TMath::Abs(momentumtest5a.X() - fRecoMomentum[i].X())
1619 / momentumtest5a.X());
1620 fhUsedMomenta_errorY_stsMomVec->Fill(TMath::Abs(momentumtest5a.Y() - fRecoMomentum[i].Y())
1621 / momentumtest5a.Y());
1622 fhUsedMomenta_errorZ_stsMomVec->Fill(TMath::Abs(momentumtest5a.Z() - fRecoMomentum[i].Z())
1623 / momentumtest5a.Z());
1624 fhUsedMomenta_vsX_stsMomVec->Fill(momentumtest5a.X(), fRecoMomentum[i].X());
1625 fhUsedMomenta_vsY_stsMomVec->Fill(momentumtest5a.Y(), fRecoMomentum[i].Y());
1626 fhUsedMomenta_vsZ_stsMomVec->Fill(momentumtest5a.Z(), fRecoMomentum[i].Z());
1627 fhUsedMomenta_errorX_refitted->Fill(TMath::Abs(momentumtest5a.X() - fRecoRefittedMomentum[i].X())
1628 / momentumtest5a.X());
1629 fhUsedMomenta_errorY_refitted->Fill(TMath::Abs(momentumtest5a.Y() - fRecoRefittedMomentum[i].Y())
1630 / momentumtest5a.Y());
1631 fhUsedMomenta_errorZ_refitted->Fill(TMath::Abs(momentumtest5a.Z() - fRecoRefittedMomentum[i].Z())
1632 / momentumtest5a.Z());
1633 fhUsedMomenta_vsX_refitted->Fill(momentumtest5a.X(), fRecoRefittedMomentum[i].X());
1634 fhUsedMomenta_vsY_refitted->Fill(momentumtest5a.Y(), fRecoRefittedMomentum[i].Y());
1635 fhUsedMomenta_vsZ_refitted->Fill(momentumtest5a.Z(), fRecoRefittedMomentum[i].Z());
1636
1637 TVector3 momentumtest5b;
1638 fRecoTracklistEPEM[j]->GetMomentum(momentumtest5b);
1639 fhUsedMomenta_errorX_stsMomVec->Fill(TMath::Abs(momentumtest5b.X() - fRecoMomentum[j].X())
1640 / momentumtest5b.X());
1641 fhUsedMomenta_errorY_stsMomVec->Fill(TMath::Abs(momentumtest5b.Y() - fRecoMomentum[j].Y())
1642 / momentumtest5b.Y());
1643 fhUsedMomenta_errorZ_stsMomVec->Fill(TMath::Abs(momentumtest5b.Z() - fRecoMomentum[j].Z())
1644 / momentumtest5b.Z());
1645 fhUsedMomenta_vsX_stsMomVec->Fill(momentumtest5b.X(), fRecoMomentum[j].X());
1646 fhUsedMomenta_vsY_stsMomVec->Fill(momentumtest5b.Y(), fRecoMomentum[j].Y());
1647 fhUsedMomenta_vsZ_stsMomVec->Fill(momentumtest5b.Z(), fRecoMomentum[j].Z());
1648 fhUsedMomenta_errorX_refitted->Fill(TMath::Abs(momentumtest5b.X() - fRecoRefittedMomentum[j].X())
1649 / momentumtest5b.X());
1650 fhUsedMomenta_errorY_refitted->Fill(TMath::Abs(momentumtest5b.Y() - fRecoRefittedMomentum[j].Y())
1651 / momentumtest5b.Y());
1652 fhUsedMomenta_errorZ_refitted->Fill(TMath::Abs(momentumtest5b.Z() - fRecoRefittedMomentum[j].Z())
1653 / momentumtest5b.Z());
1654 fhUsedMomenta_vsX_refitted->Fill(momentumtest5b.X(), fRecoRefittedMomentum[j].X());
1655 fhUsedMomenta_vsY_refitted->Fill(momentumtest5b.Y(), fRecoRefittedMomentum[j].Y());
1656 fhUsedMomenta_vsZ_refitted->Fill(momentumtest5b.Z(), fRecoRefittedMomentum[j].Z());
1657
1658 TVector3 momentumtest5c;
1659 fRecoTracklistEPEM[k]->GetMomentum(momentumtest5c);
1660 fhUsedMomenta_errorX_stsMomVec->Fill(TMath::Abs(momentumtest5c.X() - fRecoMomentum[k].X())
1661 / momentumtest5c.X());
1662 fhUsedMomenta_errorY_stsMomVec->Fill(TMath::Abs(momentumtest5c.Y() - fRecoMomentum[k].Y())
1663 / momentumtest5c.Y());
1664 fhUsedMomenta_errorZ_stsMomVec->Fill(TMath::Abs(momentumtest5c.Z() - fRecoMomentum[k].Z())
1665 / momentumtest5c.Z());
1666 fhUsedMomenta_vsX_stsMomVec->Fill(momentumtest5c.X(), fRecoMomentum[k].X());
1667 fhUsedMomenta_vsY_stsMomVec->Fill(momentumtest5c.Y(), fRecoMomentum[k].Y());
1668 fhUsedMomenta_vsZ_stsMomVec->Fill(momentumtest5c.Z(), fRecoMomentum[k].Z());
1669 fhUsedMomenta_errorX_refitted->Fill(TMath::Abs(momentumtest5c.X() - fRecoRefittedMomentum[k].X())
1670 / momentumtest5c.X());
1671 fhUsedMomenta_errorY_refitted->Fill(TMath::Abs(momentumtest5c.Y() - fRecoRefittedMomentum[k].Y())
1672 / momentumtest5c.Y());
1673 fhUsedMomenta_errorZ_refitted->Fill(TMath::Abs(momentumtest5c.Z() - fRecoRefittedMomentum[k].Z())
1674 / momentumtest5c.Z());
1675 fhUsedMomenta_vsX_refitted->Fill(momentumtest5c.X(), fRecoRefittedMomentum[k].X());
1676 fhUsedMomenta_vsY_refitted->Fill(momentumtest5c.Y(), fRecoRefittedMomentum[k].Y());
1677 fhUsedMomenta_vsZ_refitted->Fill(momentumtest5c.Z(), fRecoRefittedMomentum[k].Z());
1678
1679 TVector3 momentumtest5d;
1680 fRecoTracklistEPEM[l]->GetMomentum(momentumtest5d);
1681 fhUsedMomenta_errorX_stsMomVec->Fill(TMath::Abs(momentumtest5d.X() - fRecoMomentum[l].X())
1682 / momentumtest5d.X());
1683 fhUsedMomenta_errorY_stsMomVec->Fill(TMath::Abs(momentumtest5d.Y() - fRecoMomentum[l].Y())
1684 / momentumtest5d.Y());
1685 fhUsedMomenta_errorZ_stsMomVec->Fill(TMath::Abs(momentumtest5d.Z() - fRecoMomentum[l].Z())
1686 / momentumtest5d.Z());
1687 fhUsedMomenta_vsX_stsMomVec->Fill(momentumtest5d.X(), fRecoMomentum[l].X());
1688 fhUsedMomenta_vsY_stsMomVec->Fill(momentumtest5d.Y(), fRecoMomentum[l].Y());
1689 fhUsedMomenta_vsZ_stsMomVec->Fill(momentumtest5d.Z(), fRecoMomentum[l].Z());
1690 fhUsedMomenta_errorX_refitted->Fill(TMath::Abs(momentumtest5d.X() - fRecoRefittedMomentum[l].X())
1691 / momentumtest5d.X());
1692 fhUsedMomenta_errorY_refitted->Fill(TMath::Abs(momentumtest5d.Y() - fRecoRefittedMomentum[l].Y())
1693 / momentumtest5d.Y());
1694 fhUsedMomenta_errorZ_refitted->Fill(TMath::Abs(momentumtest5d.Z() - fRecoRefittedMomentum[l].Z())
1695 / momentumtest5d.Z());
1696 fhUsedMomenta_vsX_refitted->Fill(momentumtest5d.X(), fRecoRefittedMomentum[l].X());
1697 fhUsedMomenta_vsY_refitted->Fill(momentumtest5d.Y(), fRecoRefittedMomentum[l].Y());
1698 fhUsedMomenta_vsZ_refitted->Fill(momentumtest5d.Z(), fRecoRefittedMomentum[l].Z());
1699
1700
1701 Double_t opening_angle1_mc = 0;
1702 Double_t opening_angle2_mc = 0;
1703 Double_t opening_angle1_refitted = 0;
1704 Double_t opening_angle2_refitted = 0;
1705
1706 if (motherId1 == motherId2) {
1707 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[j]);
1708 opening_angle2_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[k], fRecoTracklistEPEM[l]);
1709 opening_angle1_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j]);
1710 opening_angle2_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1711 }
1712 if (motherId1 == motherId3) {
1713 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[k]);
1714 opening_angle2_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[j], fRecoTracklistEPEM[l]);
1715 opening_angle1_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[k]);
1716 opening_angle2_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[l]);
1717 }
1718 if (motherId1 == motherId4) {
1719 opening_angle1_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[i], fRecoTracklistEPEM[l]);
1720 opening_angle2_mc = CalculateOpeningAngleMC(fRecoTracklistEPEM[j], fRecoTracklistEPEM[k]);
1721 opening_angle1_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[l]);
1722 opening_angle2_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[k]);
1723 }
1724
1725
1726 fhInvMass_EPEM_openingAngleRef->Fill(opening_angle1_refitted);
1727 fhInvMass_EPEM_openingAngleRef->Fill(opening_angle2_refitted);
1728
1729
1730 fhEPEM_invmass_gg_mc->Fill(invmass1);
1731 fhEPEM_invmass_gg_refitted->Fill(invmass3);
1732 fhEPEM_invmass_all_mc->Fill(invmass1);
1733 fhEPEM_invmass_all_refitted->Fill(invmass3);
1734 fhEPEM_openingAngle_gg_mc->Fill(opening_angle1_mc);
1735 fhEPEM_openingAngle_gg_mc->Fill(opening_angle2_mc);
1736 fhEPEM_openingAngle_gg_refitted->Fill(opening_angle1_refitted);
1737 fhEPEM_openingAngle_gg_refitted->Fill(opening_angle2_refitted);
1738
1739 Double_t openingAngleBetweenGammas = CalculateOpeningAngleBetweenGammasMC(
1740 fRecoTracklistEPEM[i], fRecoTracklistEPEM[j], fRecoTracklistEPEM[k], fRecoTracklistEPEM[l]);
1741 fhEPEM_openingAngle_betweenGammas_mc->Fill(openingAngleBetweenGammas);
1742 Double_t openingAngleBetweenGammasReco = CalculateOpeningAngleBetweenGammasReco(
1743 fRecoRefittedMomentum[i], fRecoRefittedMomentum[j], fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1744 fhEPEM_openingAngle_betweenGammas_reco->Fill(openingAngleBetweenGammasReco);
1745
1746
1747 LmvmKinePar params1 = CalculateKinematicParams_4particles(
1748 fRecoRefittedMomentum[i], fRecoRefittedMomentum[j], fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1749 fhPi0_pt_vs_rap_gg->Fill(params1.fPt, params1.fRapidity);
1750 fhPi0_pt_vs_rap_all->Fill(params1.fPt, params1.fRapidity);
1751 fhPi0_pt_gg->Fill(params1.fPt);
1752 fhPi0_pt_all->Fill(params1.fPt);
1753
1754
1755 fhEPEM_openingAngle_vs_pt_gg_reco->Fill(params1.fPt, opening_angle1_refitted);
1756 fhEPEM_openingAngle_vs_pt_gg_reco->Fill(params1.fPt, opening_angle2_refitted);
1757
1758 fhEPEM_rap_vs_chi->Fill(params1.fRapidity, fRecoTracklistEPEM_chi[i]);
1759 fhEPEM_rap_vs_chi->Fill(params1.fRapidity, fRecoTracklistEPEM_chi[j]);
1760 fhEPEM_rap_vs_chi->Fill(params1.fRapidity, fRecoTracklistEPEM_chi[k]);
1761 fhEPEM_rap_vs_chi->Fill(params1.fRapidity, fRecoTracklistEPEM_chi[l]);
1762
1763 fhEPEM_rap_vs_invmass->Fill(params1.fRapidity, invmass3);
1764
1765
1766 // ####################################
1767 // STUDIES: efficiency of cuts
1768 // ####################################
1769 //fRichElIdAnn = new CbmRichElectronIdAnn();
1770 //fRichElIdAnn->Init();
1771 //Double_t ann = fRichElIdAnn->DoSelect(ring, momentum);
1772 //if (ann > fCuts.fRichAnnCut) return true; // cut = 0.0
1773
1774
1775 CutEfficiencyStudies(i, j, k, l, motherId1, motherId2, motherId3, motherId4);
1776
1777 //Bool_t IsRichElectron1 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[i], fRecoRefittedMomentum[i].Mag());
1778 //Bool_t IsRichElectron2 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[j], fRecoRefittedMomentum[j].Mag());
1779 //Bool_t IsRichElectron3 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[k], fRecoRefittedMomentum[k].Mag());
1780 //Bool_t IsRichElectron4 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[l], fRecoRefittedMomentum[l].Mag());
1781
1782 Bool_t IsRichElectron1ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[i], fRecoRefittedMomentum[i].Mag());
1783 Bool_t IsRichElectron2ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[j], fRecoRefittedMomentum[j].Mag());
1784 Bool_t IsRichElectron3ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[k], fRecoRefittedMomentum[k].Mag());
1785 Bool_t IsRichElectron4ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[l], fRecoRefittedMomentum[l].Mag());
1786
1787 Bool_t IsRichElectron1normal =
1788 IsRichElectronNormal(fRecoTracklistEPEM_gtid[i], fRecoRefittedMomentum[i].Mag());
1789 Bool_t IsRichElectron2normal =
1790 IsRichElectronNormal(fRecoTracklistEPEM_gtid[j], fRecoRefittedMomentum[j].Mag());
1791 Bool_t IsRichElectron3normal =
1792 IsRichElectronNormal(fRecoTracklistEPEM_gtid[k], fRecoRefittedMomentum[k].Mag());
1793 Bool_t IsRichElectron4normal =
1794 IsRichElectronNormal(fRecoTracklistEPEM_gtid[l], fRecoRefittedMomentum[l].Mag());
1795
1796 Bool_t IsRichElectron1MC = (TMath::Abs(fRecoTracklistEPEM[i]->GetPdgCode()) == 11);
1797 Bool_t IsRichElectron2MC = (TMath::Abs(fRecoTracklistEPEM[j]->GetPdgCode()) == 11);
1798 Bool_t IsRichElectron3MC = (TMath::Abs(fRecoTracklistEPEM[k]->GetPdgCode()) == 11);
1799 Bool_t IsRichElectron4MC = (TMath::Abs(fRecoTracklistEPEM[l]->GetPdgCode()) == 11);
1800
1801 LmvmKinePar paramsCut1;
1802 LmvmKinePar paramsCut2;
1803
1804 if (motherId1 == motherId2) {
1805 paramsCut1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[j]);
1806 paramsCut2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[k], fRecoRefittedMomentum[l]);
1807 }
1808 if (motherId1 == motherId3) {
1809 paramsCut1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[k]);
1810 paramsCut2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[l]);
1811 }
1812 if (motherId1 == motherId4) {
1813 paramsCut1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[i], fRecoRefittedMomentum[l]);
1814 paramsCut2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[j], fRecoRefittedMomentum[k]);
1815 }
1816
1817
1818 Double_t Value_invariantMassCut = 0.03;
1819 Double_t Value_openingAngleCut1 = 1.8 - 0.6 * paramsCut1.fPt;
1820 Double_t Value_openingAngleCut2 = 1.8 - 0.6 * paramsCut2.fPt;
1821
1822 Bool_t OpeningAngleCut1 = (opening_angle1_refitted < Value_openingAngleCut1);
1823 Bool_t OpeningAngleCut2 = (opening_angle2_refitted < Value_openingAngleCut2);
1824 Bool_t InvariantMassCut1 = (paramsCut1.fMinv < Value_invariantMassCut);
1825 Bool_t InvariantMassCut2 = (paramsCut2.fMinv < Value_invariantMassCut);
1826
1827
1828 fhEPEM_efficiencyCuts2->Fill(0); // no further cuts applied
1829 // first ANN usage for electron identification
1830 if (IsRichElectron1ann && IsRichElectron2ann && IsRichElectron3ann
1831 && IsRichElectron4ann) { // all 4 electrons correctly identified with the RICH via ANN
1832 fhEPEM_efficiencyCuts2->Fill(1);
1833 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
1834 fhEPEM_efficiencyCuts2->Fill(2);
1835 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts2->Fill(3); }
1836 }
1837 }
1838 // then standard method for electron identification
1839 if (IsRichElectron1normal && IsRichElectron2normal && IsRichElectron3normal
1840 && IsRichElectron4normal) { // all 4 electrons correctly identified with the RICH via "normal way"
1841 fhEPEM_efficiencyCuts2->Fill(4);
1842 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
1843 fhEPEM_efficiencyCuts2->Fill(5);
1844 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts2->Fill(6); }
1845 }
1846 }
1847 // MC-true data for electron identification
1848 if (IsRichElectron1MC && IsRichElectron2MC && IsRichElectron3MC
1849 && IsRichElectron4MC) { // all 4 electrons correctly identified with the RICH via MC-true data
1850 fhEPEM_efficiencyCuts2->Fill(7);
1851 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
1852 fhEPEM_efficiencyCuts2->Fill(8);
1853 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts2->Fill(9); }
1854 }
1855 }
1856
1857
1858 cout << "reco/mc: " << fRecoMomentum[i].Mag() << " / " << fRecoTracklistEPEM[i]->GetP() << " ### "
1859 << fRecoMomentum[j].Mag() << " / " << fRecoTracklistEPEM[j]->GetP() << " ### "
1860 << fRecoMomentum[k].Mag() << " / " << fRecoTracklistEPEM[k]->GetP() << " ### "
1861 << fRecoMomentum[l].Mag() << " / " << fRecoTracklistEPEM[l]->GetP() << endl;
1862
1863 fill++;
1864 }
1865 }
1866 }
1867 }
1868 }
1869 }
1870 cout << "CbmAnaConversionReco: InvariantMassTest_4epem - Filled events: " << fill << endl;
1871 cout << "CbmAnaConversionReco: InvariantMassTest_4epem - End!" << endl;
1872
1873
1874 timer.Stop();
1875 fTime += timer.RealTime();
1876}
1877
1878
1879void CbmAnaConversionReco::CutEfficiencyStudies(int e1, int e2, int e3, int e4, int motherE1, int motherE2,
1880 int motherE3, int motherE4, int IsEta)
1881{
1882 // ####################################
1883 // STUDIES: efficiency of cuts
1884 // ####################################
1885 //fRichElIdAnn = new CbmRichElectronIdAnn();
1886 //fRichElIdAnn->Init();
1887 //Double_t ann = fRichElIdAnn->DoSelect(ring, momentum);
1888 //if (ann > fCuts.fRichAnnCut) return true; // cut = 0.0
1889
1890 Double_t opening_angle1_refitted = 0;
1891 Double_t opening_angle2_refitted = 0;
1892
1893 if (motherE1 == motherE2) {
1894 opening_angle1_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[e1], fRecoRefittedMomentum[e2]);
1895 opening_angle2_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[e3], fRecoRefittedMomentum[e4]);
1896 }
1897 if (motherE1 == motherE3) {
1898 opening_angle1_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[e1], fRecoRefittedMomentum[e3]);
1899 opening_angle2_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[e2], fRecoRefittedMomentum[e4]);
1900 }
1901 if (motherE1 == motherE4) {
1902 opening_angle1_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[e1], fRecoRefittedMomentum[e4]);
1903 opening_angle2_refitted = CalculateOpeningAngleReco(fRecoRefittedMomentum[e2], fRecoRefittedMomentum[e3]);
1904 }
1905
1906
1907 Bool_t IsWithinChiCut1 =
1908 (fRecoTracklistEPEM_chi[e1] < CbmAnaConversionCutSettings::CalcChiCut(fRecoRefittedMomentum[e1].Perp()));
1909 Bool_t IsWithinChiCut2 =
1910 (fRecoTracklistEPEM_chi[e2] < CbmAnaConversionCutSettings::CalcChiCut(fRecoRefittedMomentum[e2].Perp()));
1911 Bool_t IsWithinChiCut3 =
1912 (fRecoTracklistEPEM_chi[e3] < CbmAnaConversionCutSettings::CalcChiCut(fRecoRefittedMomentum[e3].Perp()));
1913 Bool_t IsWithinChiCut4 =
1914 (fRecoTracklistEPEM_chi[e4] < CbmAnaConversionCutSettings::CalcChiCut(fRecoRefittedMomentum[e4].Perp()));
1915
1916 Bool_t AllWithinChiCut = (IsWithinChiCut1 && IsWithinChiCut2 && IsWithinChiCut3 && IsWithinChiCut4);
1917
1918
1919 //Bool_t IsRichElectron1 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[i], fRecoRefittedMomentum[i].Mag());
1920 //Bool_t IsRichElectron2 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[j], fRecoRefittedMomentum[j].Mag());
1921 //Bool_t IsRichElectron3 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[k], fRecoRefittedMomentum[k].Mag());
1922 //Bool_t IsRichElectron4 = CbmLitGlobalElectronId::GetInstance().IsRichElectron(fRecoTracklistEPEM_gtid[l], fRecoRefittedMomentum[l].Mag());
1923
1924 Bool_t IsRichElectron1ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[e1], fRecoRefittedMomentum[e1].Mag());
1925 Bool_t IsRichElectron2ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[e2], fRecoRefittedMomentum[e2].Mag());
1926 Bool_t IsRichElectron3ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[e3], fRecoRefittedMomentum[e3].Mag());
1927 Bool_t IsRichElectron4ann = IsRichElectronANN(fRecoTracklistEPEM_gtid[e4], fRecoRefittedMomentum[e4].Mag());
1928
1929 Bool_t IsRichElectron1normal = IsRichElectronNormal(fRecoTracklistEPEM_gtid[e1], fRecoRefittedMomentum[e1].Mag());
1930 Bool_t IsRichElectron2normal = IsRichElectronNormal(fRecoTracklistEPEM_gtid[e2], fRecoRefittedMomentum[e2].Mag());
1931 Bool_t IsRichElectron3normal = IsRichElectronNormal(fRecoTracklistEPEM_gtid[e3], fRecoRefittedMomentum[e3].Mag());
1932 Bool_t IsRichElectron4normal = IsRichElectronNormal(fRecoTracklistEPEM_gtid[e4], fRecoRefittedMomentum[e4].Mag());
1933
1934 Bool_t IsRichElectron1MC = (TMath::Abs(fRecoTracklistEPEM[e1]->GetPdgCode()) == 11);
1935 Bool_t IsRichElectron2MC = (TMath::Abs(fRecoTracklistEPEM[e2]->GetPdgCode()) == 11);
1936 Bool_t IsRichElectron3MC = (TMath::Abs(fRecoTracklistEPEM[e3]->GetPdgCode()) == 11);
1937 Bool_t IsRichElectron4MC = (TMath::Abs(fRecoTracklistEPEM[e4]->GetPdgCode()) == 11);
1938
1939 Double_t ANNvalueE1 = ElectronANNvalue(fRecoTracklistEPEM_gtid[e1], fRecoRefittedMomentum[e1].Mag());
1940 Double_t ANNvalueE2 = ElectronANNvalue(fRecoTracklistEPEM_gtid[e2], fRecoRefittedMomentum[e2].Mag());
1941 Double_t ANNvalueE3 = ElectronANNvalue(fRecoTracklistEPEM_gtid[e3], fRecoRefittedMomentum[e3].Mag());
1942 Double_t ANNvalueE4 = ElectronANNvalue(fRecoTracklistEPEM_gtid[e4], fRecoRefittedMomentum[e4].Mag());
1943
1944 LmvmKinePar paramsCut1;
1945 LmvmKinePar paramsCut2;
1946
1947 if (motherE1 == motherE2) {
1948 paramsCut1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[e1], fRecoRefittedMomentum[e2]);
1949 paramsCut2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[e3], fRecoRefittedMomentum[e4]);
1950 }
1951 if (motherE1 == motherE3) {
1952 paramsCut1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[e1], fRecoRefittedMomentum[e3]);
1953 paramsCut2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[e2], fRecoRefittedMomentum[e4]);
1954 }
1955 if (motherE1 == motherE4) {
1956 paramsCut1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[e1], fRecoRefittedMomentum[e4]);
1957 paramsCut2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[e2], fRecoRefittedMomentum[e3]);
1958 }
1959
1960
1961 Double_t Value_invariantMassCut = 0.03;
1962 Double_t Value_openingAngleCut1 = 1.8 - 0.6 * paramsCut1.fPt;
1963 Double_t Value_openingAngleCut2 = 1.8 - 0.6 * paramsCut2.fPt;
1964
1965 Bool_t OpeningAngleCut1 = (opening_angle1_refitted < Value_openingAngleCut1);
1966 Bool_t OpeningAngleCut2 = (opening_angle2_refitted < Value_openingAngleCut2);
1967 Bool_t InvariantMassCut1 = (paramsCut1.fMinv < Value_invariantMassCut);
1968 Bool_t InvariantMassCut2 = (paramsCut2.fMinv < Value_invariantMassCut);
1969
1970
1971 if (IsEta == 0) {
1972 fhEPEM_pi0_nofLeptons_ann->Fill(IsRichElectron1ann + IsRichElectron2ann + IsRichElectron3ann + IsRichElectron4ann);
1973
1974 if (IsRichElectron1ann + IsRichElectron2ann + IsRichElectron3ann + IsRichElectron4ann == 3) {
1975 CbmGlobalTrack* gTrack = NULL;
1976 if (IsRichElectron1ann == 0) gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e1]);
1977 if (IsRichElectron2ann == 0) gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e2]);
1978 if (IsRichElectron3ann == 0) gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e3]);
1979 if (IsRichElectron4ann == 0) gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e4]);
1980 int richInd = gTrack->GetRichRingIndex();
1981 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(richInd);
1982 Int_t nofringhits = richRing->GetNofHits();
1983 fhEPEM_missingLepton_nofRingHits->Fill(nofringhits);
1984 fhEPEM_missingLepton_ringMid->Fill(richRing->GetCenterX(), richRing->GetCenterY());
1985 fhEPEM_missingLepton_ringRadius->Fill(richRing->GetRadius());
1986 // CbmRichRing::GetDistance() method is no longer supported
1987 // If you wan to use cuts update code using CbmRichUtil::GetRingTrackDistance()
1988 fhEPEM_missingLepton_distance->Fill(1. /*richRing->GetDistance()*/);
1989 fhEPEM_missingLepton_selectionNN->Fill(richRing->GetSelectionNN());
1990 if (IsRichElectron1ann == 0) {
1991 fhEPEM_missingLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e1].Mag(), richRing->GetRadius());
1992 fhEPEM_missingLepton_ANNvalue->Fill(ANNvalueE1);
1993 }
1994 if (IsRichElectron2ann == 0) {
1995 fhEPEM_missingLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e2].Mag(), richRing->GetRadius());
1996 fhEPEM_missingLepton_ANNvalue->Fill(ANNvalueE2);
1997 }
1998 if (IsRichElectron3ann == 0) {
1999 fhEPEM_missingLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e3].Mag(), richRing->GetRadius());
2000 fhEPEM_missingLepton_ANNvalue->Fill(ANNvalueE3);
2001 }
2002 if (IsRichElectron4ann == 0) {
2003 fhEPEM_missingLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e4].Mag(), richRing->GetRadius());
2004 fhEPEM_missingLepton_ANNvalue->Fill(ANNvalueE4);
2005 }
2006
2007 for (int i = 0; i < nofringhits; i++) {
2008 Int_t hitInd = richRing->GetHit(i);
2009 CbmRichHit* hit = (CbmRichHit*) fRichHits->At(hitInd);
2010 if (NULL == hit) continue;
2011 //CbmRichHitLight hl(hit->GetX(), hit->GetY());
2012 fhEPEM_missingLepton_rings->Fill(hit->GetX(), hit->GetY());
2013 }
2014 }
2015
2016
2017 if (IsRichElectron1ann > 0) {
2018 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e1]);
2019 int richInd = gTrack->GetRichRingIndex();
2020 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(richInd);
2021 Int_t nofringhits = richRing->GetNofHits();
2022 fhEPEM_identifiedLepton_nofRingHits->Fill(nofringhits);
2023 fhEPEM_identifiedLepton_ringMid->Fill(richRing->GetCenterX(), richRing->GetCenterY());
2024 fhEPEM_identifiedLepton_ringRadius->Fill(richRing->GetRadius());
2025 // CbmRichRing::GetDistance() method is no longer supported
2026 // If you wan to use cuts update code using CbmRichUtil::GetRingTrackDistance()
2027 fhEPEM_identifiedLepton_distance->Fill(1. /*richRing->GetDistance()*/);
2028 fhEPEM_identifiedLepton_selectionNN->Fill(richRing->GetSelectionNN());
2029 fhEPEM_identifiedLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e1].Mag(), richRing->GetRadius());
2030 fhEPEM_identifiedLepton_ANNvalue->Fill(ANNvalueE1);
2031 }
2032 if (IsRichElectron2ann > 0) {
2033 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e2]);
2034 int richInd = gTrack->GetRichRingIndex();
2035 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(richInd);
2036 Int_t nofringhits = richRing->GetNofHits();
2037 fhEPEM_identifiedLepton_nofRingHits->Fill(nofringhits);
2038 fhEPEM_identifiedLepton_ringMid->Fill(richRing->GetCenterX(), richRing->GetCenterY());
2039 fhEPEM_identifiedLepton_ringRadius->Fill(richRing->GetRadius());
2040 // CbmRichRing::GetDistance() method is no longer supported
2041 // If you wan to use cuts update code using CbmRichUtil::GetRingTrackDistance()
2042 fhEPEM_identifiedLepton_distance->Fill(1. /*richRing->GetDistance()*/);
2043 fhEPEM_identifiedLepton_selectionNN->Fill(richRing->GetSelectionNN());
2044 fhEPEM_identifiedLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e2].Mag(), richRing->GetRadius());
2045 fhEPEM_identifiedLepton_ANNvalue->Fill(ANNvalueE2);
2046 }
2047 if (IsRichElectron3ann > 0) {
2048 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e3]);
2049 int richInd = gTrack->GetRichRingIndex();
2050 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(richInd);
2051 Int_t nofringhits = richRing->GetNofHits();
2052 fhEPEM_identifiedLepton_nofRingHits->Fill(nofringhits);
2053 fhEPEM_identifiedLepton_ringMid->Fill(richRing->GetCenterX(), richRing->GetCenterY());
2054 fhEPEM_identifiedLepton_ringRadius->Fill(richRing->GetRadius());
2055 // CbmRichRing::GetDistance() method is no longer supported
2056 // If you wan to use cuts update code using CbmRichUtil::GetRingTrackDistance()
2057 fhEPEM_identifiedLepton_distance->Fill(1. /*richRing->GetDistance() */);
2058 fhEPEM_identifiedLepton_selectionNN->Fill(richRing->GetSelectionNN());
2059 fhEPEM_identifiedLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e3].Mag(), richRing->GetRadius());
2060 fhEPEM_identifiedLepton_ANNvalue->Fill(ANNvalueE3);
2061 }
2062 if (IsRichElectron4ann > 0) {
2063 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(fRecoTracklistEPEM_gtid[e4]);
2064 int richInd = gTrack->GetRichRingIndex();
2065 CbmRichRing* richRing = (CbmRichRing*) fRichRings->At(richInd);
2066 Int_t nofringhits = richRing->GetNofHits();
2067 fhEPEM_identifiedLepton_nofRingHits->Fill(nofringhits);
2068 fhEPEM_identifiedLepton_ringMid->Fill(richRing->GetCenterX(), richRing->GetCenterY());
2069 fhEPEM_identifiedLepton_ringRadius->Fill(richRing->GetRadius());
2070 // CbmRichRing::GetDistance() method is no longer supported
2071 // If you wan to use cuts update code using CbmRichUtil::GetRingTrackDistance()
2072 fhEPEM_identifiedLepton_distance->Fill(1. /*richRing->GetDistance()*/);
2073 fhEPEM_identifiedLepton_selectionNN->Fill(richRing->GetSelectionNN());
2074 fhEPEM_identifiedLepton_radius_vs_p->Fill(fRecoRefittedMomentum[e4].Mag(), richRing->GetRadius());
2075 fhEPEM_identifiedLepton_ANNvalue->Fill(ANNvalueE4);
2076 }
2077
2078
2079 fhEPEM_pi0_ANNvalues_noCuts->Fill(ANNvalueE1);
2080 fhEPEM_pi0_ANNvalues_noCuts->Fill(ANNvalueE2);
2081 fhEPEM_pi0_ANNvalues_noCuts->Fill(ANNvalueE3);
2082 fhEPEM_pi0_ANNvalues_noCuts->Fill(ANNvalueE4);
2083 if (OpeningAngleCut1 && OpeningAngleCut2) {
2084 fhEPEM_pi0_ANNvalues_angleCut->Fill(ANNvalueE1);
2085 fhEPEM_pi0_ANNvalues_angleCut->Fill(ANNvalueE2);
2086 fhEPEM_pi0_ANNvalues_angleCut->Fill(ANNvalueE3);
2087 fhEPEM_pi0_ANNvalues_angleCut->Fill(ANNvalueE4);
2088 }
2089
2090 fhEPEM_pi0_ANNefficiencies->Fill(0);
2091 if (ANNvalueE1 > -1 && ANNvalueE2 > -1 && ANNvalueE3 > -1 && ANNvalueE4 > -1) {
2092 fhEPEM_pi0_ANNefficiencies->Fill(1);
2093 }
2094 if (ANNvalueE1 > -0.9 && ANNvalueE2 > -0.9 && ANNvalueE3 > -0.9 && ANNvalueE4 > -0.9) {
2095 fhEPEM_pi0_ANNefficiencies->Fill(2);
2096 }
2097 if (ANNvalueE1 > -0.8 && ANNvalueE2 > -0.8 && ANNvalueE3 > -0.8 && ANNvalueE4 > -0.8) {
2098 fhEPEM_pi0_ANNefficiencies->Fill(3);
2099 }
2100 if (ANNvalueE1 > -0.7 && ANNvalueE2 > -0.7 && ANNvalueE3 > -0.7 && ANNvalueE4 > -0.7) {
2101 fhEPEM_pi0_ANNefficiencies->Fill(4);
2102 }
2103 if (ANNvalueE1 > -0.6 && ANNvalueE2 > -0.6 && ANNvalueE3 > -0.6 && ANNvalueE4 > -0.6) {
2104 fhEPEM_pi0_ANNefficiencies->Fill(5);
2105 }
2106 if (ANNvalueE1 > -0.5 && ANNvalueE2 > -0.5 && ANNvalueE3 > -0.5 && ANNvalueE4 > -0.5) {
2107 fhEPEM_pi0_ANNefficiencies->Fill(6);
2108 }
2109 if (ANNvalueE1 > -0.4 && ANNvalueE2 > -0.4 && ANNvalueE3 > -0.4 && ANNvalueE4 > -0.4) {
2110 fhEPEM_pi0_ANNefficiencies->Fill(7);
2111 }
2112 if (ANNvalueE1 > -0.3 && ANNvalueE2 > -0.3 && ANNvalueE3 > -0.3 && ANNvalueE4 > -0.3) {
2113 fhEPEM_pi0_ANNefficiencies->Fill(8);
2114 }
2115 if (ANNvalueE1 > -0.2 && ANNvalueE2 > -0.2 && ANNvalueE3 > -0.2 && ANNvalueE4 > -0.2) {
2116 fhEPEM_pi0_ANNefficiencies->Fill(9);
2117 }
2118 if (ANNvalueE1 > -0.1 && ANNvalueE2 > -0.1 && ANNvalueE3 > -0.1 && ANNvalueE4 > -0.1) {
2119 fhEPEM_pi0_ANNefficiencies->Fill(10);
2120 }
2121 if (ANNvalueE1 > -0.0 && ANNvalueE2 > -0.0 && ANNvalueE3 > -0.0 && ANNvalueE4 > -0.0) {
2122 fhEPEM_pi0_ANNefficiencies->Fill(11);
2123 }
2124
2125 fhEPEM_efficiencyCuts->Fill(0); // no further cuts applied
2126 // first ANN usage for electron identification
2127 if (IsRichElectron1ann && IsRichElectron2ann && IsRichElectron3ann
2128 && IsRichElectron4ann) { // all 4 electrons correctly identified with the RICH via ANN
2129 fhEPEM_efficiencyCuts->Fill(1);
2130 if (AllWithinChiCut) { // refitted momenta are within chi cut
2131 fhEPEM_efficiencyCuts->Fill(2);
2132 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
2133 fhEPEM_efficiencyCuts->Fill(3);
2134 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts->Fill(4); }
2135 }
2136 }
2137 }
2138 // then standard method for electron identification
2139 if (IsRichElectron1normal && IsRichElectron2normal && IsRichElectron3normal
2140 && IsRichElectron4normal) { // all 4 electrons correctly identified with the RICH via "normal way"
2141 fhEPEM_efficiencyCuts->Fill(5);
2142 if (AllWithinChiCut) { // refitted momenta are within chi cut
2143 fhEPEM_efficiencyCuts->Fill(6);
2144 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
2145 fhEPEM_efficiencyCuts->Fill(7);
2146 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts->Fill(8); }
2147 }
2148 }
2149 }
2150 // MC-true data for electron identification
2151 if (IsRichElectron1MC && IsRichElectron2MC && IsRichElectron3MC
2152 && IsRichElectron4MC) { // all 4 electrons correctly identified with the RICH via MC-true data
2153 fhEPEM_efficiencyCuts->Fill(9);
2154 if (AllWithinChiCut) { // refitted momenta are within chi cut
2155 fhEPEM_efficiencyCuts->Fill(10);
2156 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
2157 fhEPEM_efficiencyCuts->Fill(11);
2158 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts->Fill(12); }
2159 }
2160 }
2161 }
2162 }
2163
2164 if (IsEta == 1) {
2165 fhEPEM_efficiencyCuts_eta->Fill(0); // no further cuts applied
2166 // first ANN usage for electron identification
2167 if (IsRichElectron1ann && IsRichElectron2ann && IsRichElectron3ann
2168 && IsRichElectron4ann) { // all 4 electrons correctly identified with the RICH via ANN
2169 fhEPEM_efficiencyCuts_eta->Fill(1);
2170 if (AllWithinChiCut) { // refitted momenta are within chi cut
2171 fhEPEM_efficiencyCuts_eta->Fill(2);
2172 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
2173 fhEPEM_efficiencyCuts_eta->Fill(3);
2174 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts_eta->Fill(4); }
2175 }
2176 }
2177 }
2178 // then standard method for electron identification
2179 if (IsRichElectron1normal && IsRichElectron2normal && IsRichElectron3normal
2180 && IsRichElectron4normal) { // all 4 electrons correctly identified with the RICH via "normal way"
2181 fhEPEM_efficiencyCuts_eta->Fill(5);
2182 if (AllWithinChiCut) { // refitted momenta are within chi cut
2183 fhEPEM_efficiencyCuts_eta->Fill(6);
2184 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
2185 fhEPEM_efficiencyCuts_eta->Fill(7);
2186 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts_eta->Fill(8); }
2187 }
2188 }
2189 }
2190 // MC-true data for electron identification
2191 if (IsRichElectron1MC && IsRichElectron2MC && IsRichElectron3MC
2192 && IsRichElectron4MC) { // all 4 electrons correctly identified with the RICH via MC-true data
2193 fhEPEM_efficiencyCuts_eta->Fill(9);
2194 if (AllWithinChiCut) { // refitted momenta are within chi cut
2195 fhEPEM_efficiencyCuts_eta->Fill(10);
2196 if (OpeningAngleCut1 && OpeningAngleCut2) { // opening angle of e+e- pairs below x
2197 fhEPEM_efficiencyCuts_eta->Fill(11);
2198 if (InvariantMassCut1 && InvariantMassCut2) { fhEPEM_efficiencyCuts_eta->Fill(12); }
2199 }
2200 }
2201 }
2202 }
2203}
2204
2205
2206Int_t CbmAnaConversionReco::NofDaughters(Int_t motherId)
2207{
2208 Int_t nofDaughters = 0;
2209 for (unsigned int i = 0; i < fRecoTracklistEPEM.size(); i++) {
2210 Int_t motherId_temp = fRecoTracklistEPEM[i]->GetMotherId();
2211 if (motherId == motherId_temp) nofDaughters++;
2212 }
2213
2214
2215 return nofDaughters;
2216}
2217
2218
2219Double_t CbmAnaConversionReco::CalculateOpeningAngleReco(TVector3 electron1, TVector3 electron2)
2220{
2221 Double_t energyP = TMath::Sqrt(electron1.Mag2() + M2E);
2222 TLorentzVector lorVecP(electron1, energyP);
2223
2224 Double_t energyM = TMath::Sqrt(electron2.Mag2() + M2E);
2225 TLorentzVector lorVecM(electron2, energyM);
2226
2227 Double_t anglePair = lorVecM.Angle(lorVecP.Vect());
2228 Double_t theta = 180. * anglePair / TMath::Pi();
2229
2230 return theta;
2231}
2232
2233
2234Double_t CbmAnaConversionReco::CalculateOpeningAngleMC(CbmMCTrack* mctrack1, CbmMCTrack* mctrack2)
2235{
2236 TVector3 electron1;
2237 mctrack1->GetMomentum(electron1);
2238 Double_t energyP = TMath::Sqrt(electron1.Mag2() + M2E);
2239 TLorentzVector lorVecP(electron1, energyP);
2240
2241 TVector3 electron2;
2242 mctrack2->GetMomentum(electron2);
2243 Double_t energyM = TMath::Sqrt(electron2.Mag2() + M2E);
2244 TLorentzVector lorVecM(electron2, energyM);
2245
2246 Double_t anglePair = lorVecM.Angle(lorVecP.Vect());
2247 Double_t theta = 180. * anglePair / TMath::Pi();
2248
2249 return theta;
2250}
2251
2252
2253Double_t CbmAnaConversionReco::CalculateOpeningAngleBetweenGammasMC(CbmMCTrack* mctrack1, CbmMCTrack* mctrack2,
2254 CbmMCTrack* mctrack3, CbmMCTrack* mctrack4)
2255{
2256 Double_t openingAngle;
2257 TLorentzVector gamma1;
2258 TLorentzVector gamma2;
2259
2260 if (mctrack1->GetMotherId() == mctrack2->GetMotherId() && mctrack3->GetMotherId() == mctrack4->GetMotherId()) {
2261 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(mctrack1->GetMotherId());
2262 mother1->Get4Momentum(gamma1);
2263 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(mctrack3->GetMotherId());
2264 mother2->Get4Momentum(gamma2);
2265 }
2266 if (mctrack1->GetMotherId() == mctrack3->GetMotherId() && mctrack2->GetMotherId() == mctrack4->GetMotherId()) {
2267 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(mctrack1->GetMotherId());
2268 mother1->Get4Momentum(gamma1);
2269 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(mctrack2->GetMotherId());
2270 mother2->Get4Momentum(gamma2);
2271 }
2272 if (mctrack1->GetMotherId() == mctrack4->GetMotherId() && mctrack2->GetMotherId() == mctrack3->GetMotherId()) {
2273 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(mctrack1->GetMotherId());
2274 mother1->Get4Momentum(gamma1);
2275 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(mctrack2->GetMotherId());
2276 mother2->Get4Momentum(gamma2);
2277 }
2278
2279 Double_t angle = gamma1.Angle(gamma2.Vect());
2280 openingAngle = 180. * angle / TMath::Pi();
2281
2282 return openingAngle;
2283}
2284
2285
2286Double_t CbmAnaConversionReco::CalculateOpeningAngleBetweenGammasReco(TVector3 electron1, TVector3 electron2,
2287 TVector3 electron3, TVector3 electron4)
2288{
2289 Double_t energy1 = TMath::Sqrt(electron1.Mag2() + M2E);
2290 TLorentzVector lorVec1(electron1, energy1);
2291
2292 Double_t energy2 = TMath::Sqrt(electron2.Mag2() + M2E);
2293 TLorentzVector lorVec2(electron2, energy2);
2294
2295 Double_t energy3 = TMath::Sqrt(electron3.Mag2() + M2E);
2296 TLorentzVector lorVec3(electron3, energy3);
2297
2298 Double_t energy4 = TMath::Sqrt(electron4.Mag2() + M2E);
2299 TLorentzVector lorVec4(electron4, energy4);
2300
2301 TLorentzVector lorPhoton1 = lorVec1 + lorVec2;
2302 TLorentzVector lorPhoton2 = lorVec3 + lorVec4;
2303
2304 Double_t angleBetweenPhotons = lorPhoton1.Angle(lorPhoton2.Vect());
2305 Double_t theta = 180. * angleBetweenPhotons / TMath::Pi();
2306
2307 return theta;
2308}
2309
2310
2311void CbmAnaConversionReco::CalculateInvMassWithFullRecoCuts()
2312{
2313 Int_t nof = fRecoTracklistEPEM.size();
2314 cout << "CbmAnaConversionReco::CalculateInvMassWithFullRecoCuts: nof entries - " << nof << endl;
2315 if (nof >= 4) {
2316 for (int a = 0; a < nof - 3; a++) {
2317 for (int b = a + 1; b < nof - 2; b++) {
2318 for (int c = b + 1; c < nof - 1; c++) {
2319 for (int d = c + 1; d < nof; d++) {
2320 Int_t check1 = (fRecoTracklistEPEM[a]->GetPdgCode() > 0);
2321 Int_t check2 = (fRecoTracklistEPEM[b]->GetPdgCode() > 0);
2322 Int_t check3 = (fRecoTracklistEPEM[c]->GetPdgCode() > 0);
2323 Int_t check4 = (fRecoTracklistEPEM[d]->GetPdgCode() > 0);
2324 Int_t test = check1 + check2 + check3 + check4;
2325 if (test != 2) continue; // need two electrons and two positrons
2326
2327
2328 Double_t invmass = Invmass_4particlesRECO(fRecoRefittedMomentum[a], fRecoRefittedMomentum[b],
2329 fRecoRefittedMomentum[c], fRecoRefittedMomentum[d]);
2330 //fhElectrons_invmass->Fill(invmass);
2331
2332
2333 LmvmKinePar params1 = CalculateKinematicParamsReco(fRecoRefittedMomentum[a], fRecoRefittedMomentum[b]);
2334 LmvmKinePar params2 = CalculateKinematicParamsReco(fRecoRefittedMomentum[a], fRecoRefittedMomentum[c]);
2335 LmvmKinePar params3 = CalculateKinematicParamsReco(fRecoRefittedMomentum[a], fRecoRefittedMomentum[d]);
2336 LmvmKinePar params4 = CalculateKinematicParamsReco(fRecoRefittedMomentum[b], fRecoRefittedMomentum[c]);
2337 LmvmKinePar params5 = CalculateKinematicParamsReco(fRecoRefittedMomentum[b], fRecoRefittedMomentum[d]);
2338 LmvmKinePar params6 = CalculateKinematicParamsReco(fRecoRefittedMomentum[c], fRecoRefittedMomentum[d]);
2339
2340 Double_t openingAngleCut = 1;
2341 Int_t IsPhoton_openingAngle1 = (params1.fAngle < openingAngleCut);
2342 Int_t IsPhoton_openingAngle2 = (params2.fAngle < openingAngleCut);
2343 Int_t IsPhoton_openingAngle3 = (params3.fAngle < openingAngleCut);
2344 Int_t IsPhoton_openingAngle4 = (params4.fAngle < openingAngleCut);
2345 Int_t IsPhoton_openingAngle5 = (params5.fAngle < openingAngleCut);
2346 Int_t IsPhoton_openingAngle6 = (params6.fAngle < openingAngleCut);
2347
2348 Double_t invMassCut = 0.03;
2349 Int_t IsPhoton_invMass1 = (params1.fMinv < invMassCut);
2350 Int_t IsPhoton_invMass2 = (params2.fMinv < invMassCut);
2351 Int_t IsPhoton_invMass3 = (params3.fMinv < invMassCut);
2352 Int_t IsPhoton_invMass4 = (params4.fMinv < invMassCut);
2353 Int_t IsPhoton_invMass5 = (params5.fMinv < invMassCut);
2354 Int_t IsPhoton_invMass6 = (params6.fMinv < invMassCut);
2355
2356 if (IsPhoton_openingAngle1 && IsPhoton_openingAngle6 && IsPhoton_invMass1 && IsPhoton_invMass6
2357 && (check1 + check2 == 1) && (check3 + check4 == 1)) {
2358 fhInvMassWithFullRecoCuts->Fill(invmass);
2359 }
2360 if (IsPhoton_openingAngle2 && IsPhoton_openingAngle5 && IsPhoton_invMass2 && IsPhoton_invMass5
2361 && (check1 + check3 == 1) && (check2 + check4 == 1)) {
2362 fhInvMassWithFullRecoCuts->Fill(invmass);
2363 }
2364 if (IsPhoton_openingAngle3 && IsPhoton_openingAngle4 && IsPhoton_invMass3 && IsPhoton_invMass4
2365 && (check1 + check4 == 1) && (check2 + check3 == 1)) {
2366 fhInvMassWithFullRecoCuts->Fill(invmass);
2367 }
2368 }
2369 }
2370 }
2371 }
2372 }
2373}
2374
2375
2376LmvmKinePar CbmAnaConversionReco::CalculateKinematicParamsReco(const TVector3 electron1, const TVector3 electron2)
2377{
2378 LmvmKinePar params;
2379
2380 Double_t energyP = TMath::Sqrt(electron1.Mag2() + M2E);
2381 TLorentzVector lorVecP(electron1, energyP);
2382
2383 Double_t energyM = TMath::Sqrt(electron2.Mag2() + M2E);
2384 TLorentzVector lorVecM(electron2, energyM);
2385
2386 TVector3 momPair = electron1 + electron2;
2387 Double_t energyPair = energyP + energyM;
2388 Double_t ptPair = momPair.Perp();
2389 Double_t pzPair = momPair.Pz();
2390 Double_t yPair = 0.5 * TMath::Log((energyPair + pzPair) / (energyPair - pzPair));
2391 Double_t anglePair = lorVecM.Angle(lorVecP.Vect());
2392 Double_t theta = 180. * anglePair / TMath::Pi();
2393 Double_t minv = 2. * TMath::Sin(anglePair / 2.) * TMath::Sqrt(electron1.Mag() * electron2.Mag());
2394
2395 params.fMomentumMag = momPair.Mag();
2396 params.fPt = ptPair;
2397 params.fRapidity = yPair;
2398 params.fMinv = minv;
2399 params.fAngle = theta;
2400 return params;
2401}
2402
2403
2404LmvmKinePar CbmAnaConversionReco::CalculateKinematicParams_4particles(const TVector3 part1, const TVector3 part2,
2405 const TVector3 part3, const TVector3 part4)
2406{
2407 LmvmKinePar params;
2408
2409 Double_t energy1 = TMath::Sqrt(part1.Mag2() + M2E);
2410 TLorentzVector lorVec1(part1, energy1);
2411
2412 Double_t energy2 = TMath::Sqrt(part2.Mag2() + M2E);
2413 TLorentzVector lorVec2(part2, energy2);
2414
2415 Double_t energy3 = TMath::Sqrt(part3.Mag2() + M2E);
2416 TLorentzVector lorVec3(part3, energy3);
2417
2418 Double_t energy4 = TMath::Sqrt(part4.Mag2() + M2E);
2419 TLorentzVector lorVec4(part4, energy4);
2420
2421 TLorentzVector sum;
2422 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
2423
2424 TVector3 momPair = part1 + part2 + part3 + part4;
2425 Double_t energyPair = energy1 + energy2 + energy3 + energy4;
2426 Double_t ptPair = momPair.Perp();
2427 Double_t pzPair = momPair.Pz();
2428 Double_t yPair = 0.5 * TMath::Log((energyPair + pzPair) / (energyPair - pzPair));
2429 // Double_t anglePair = 0;
2430 // Double_t theta = 180.*anglePair/TMath::Pi();
2431 Double_t minv = sum.Mag();
2432
2433 params.fMomentumMag = momPair.Mag();
2434 params.fPt = ptPair;
2435 params.fRapidity = yPair;
2436 params.fMinv = minv;
2437 params.fAngle = 0;
2438 return params;
2439}
2440
2441
2442Bool_t CbmAnaConversionReco::IsRichElectronANN(Int_t globalTrackIndex, Double_t momentum)
2443{
2444 if (NULL == fGlobalTracks || NULL == fRichRings) return false;
2445 //CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(globalTrackIndex);
2446 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(globalTrackIndex));
2447 Int_t richId = globalTrack->GetRichRingIndex();
2448 if (richId < 0) return false;
2449 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richId));
2450 if (NULL == ring) return false;
2451
2452 Double_t fRichAnnCut = -0.8;
2453 Double_t ann = CbmRichElectronIdAnn::GetInstance().CalculateAnnValue(globalTrackIndex, momentum);
2454 if (ann > fRichAnnCut) return true;
2455 else
2456 return false;
2457}
2458
2459
2460Double_t CbmAnaConversionReco::ElectronANNvalue(Int_t globalTrackIndex, Double_t momentum)
2461{
2462 if (NULL == fGlobalTracks || NULL == fRichRings) return -2;
2463 //CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(globalTrackIndex);
2464 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(globalTrackIndex));
2465 Int_t richId = globalTrack->GetRichRingIndex();
2466 if (richId < 0) return -2;
2467 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richId));
2468 if (NULL == ring) return -2;
2469
2470 Double_t ann = CbmRichElectronIdAnn::GetInstance().CalculateAnnValue(globalTrackIndex, momentum);
2471 return ann;
2472}
2473
2474
2475Bool_t CbmAnaConversionReco::IsRichElectronNormal(Int_t globalTrackIndex, Double_t momentum)
2476{
2477 if (NULL == fGlobalTracks || NULL == fRichRings) return false;
2478 CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(globalTrackIndex);
2479 Int_t richId = globalTrack->GetRichRingIndex();
2480 if (richId < 0) return false;
2481 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richId));
2482 if (NULL == ring) return false;
2483
2484 Double_t axisA = ring->GetAaxis();
2485 Double_t axisB = ring->GetBaxis();
2486 Double_t dist = CbmRichUtil::GetRingTrackDistance(globalTrackIndex);
2487
2488 Double_t fMeanA = 4.95;
2489 Double_t fMeanB = 4.54;
2490 Double_t fRmsA = 0.30;
2491 Double_t fRmsB = 0.22;
2492 Double_t fRmsCoeff = 4.5;
2493 Double_t fDistCut = 1.;
2494
2495 Bool_t isElectronRICH = 0;
2496
2497 if (momentum < 5.5) {
2498 if (fabs(axisA - fMeanA) < fRmsCoeff * fRmsA && fabs(axisB - fMeanB) < fRmsCoeff * fRmsB && dist < fDistCut)
2499 isElectronRICH = 1;
2500 }
2501 else {
2503 // Double_t polAaxis = 5.80008 - 4.10118 / (momentum - 3.67402);
2504 // Double_t polBaxis = 5.58839 - 4.75980 / (momentum - 3.31648);
2505 // Double_t polRaxis = 5.87252 - 7.64641/(momentum - 1.62255);
2507 Double_t polAaxis = 5.64791 - 4.24077 / (momentum - 3.65494);
2508 Double_t polBaxis = 5.41106 - 4.49902 / (momentum - 3.52450);
2509 //Double_t polRaxis = 5.66516 - 6.62229/(momentum - 2.25304);
2510 if (axisA < (fMeanA + fRmsCoeff * fRmsA) && axisA > polAaxis && axisB < (fMeanB + fRmsCoeff * fRmsB)
2511 && axisB > polBaxis && dist < fDistCut)
2512 isElectronRICH = 1;
2513 }
2514 return isElectronRICH;
2515}
M2E
#define M2E
Definition CbmAnaJpsiKinematicParams.h:14
CbmRichElectronIdAnn.h
Implementation of the electron identification algorithm in the RICH detector using Artificial Neural ...
CbmAnaConversionCutSettings::CalcChiCut
static Double_t CalcChiCut(Double_t pt)
Definition CbmAnaConversionCutSettings.h:80
CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma
static Double_t OpeningAngleBetweenGamma(const TVector3 part11, const TVector3 part12, const TVector3 part21, const TVector3 part22)
Definition CbmAnaConversionGlobalFunctions.h:80
CbmAnaConversionReco::fHistoList_gee
std::vector< TH1 * > fHistoList_gee
Definition CbmAnaConversionReco.h:94
CbmAnaConversionReco::CalculateKinematicParams_4particles
LmvmKinePar CalculateKinematicParams_4particles(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
Definition CbmAnaConversionReco.cxx:2404
CbmAnaConversionReco::fHistoList_eta
std::vector< TH1 * > fHistoList_eta
Definition CbmAnaConversionReco.h:97
CbmAnaConversionReco::ElectronANNvalue
Double_t ElectronANNvalue(Int_t globalTrackIndex, Double_t momentum)
Definition CbmAnaConversionReco.cxx:2460
CbmAnaConversionReco::IsRichElectronNormal
Bool_t IsRichElectronNormal(Int_t globalTrackIndex, Double_t momentum)
Definition CbmAnaConversionReco.cxx:2475
CbmAnaConversionReco::SetTracklistMC
void SetTracklistMC(std::vector< CbmMCTrack * > MCTracklist)
Definition CbmAnaConversionReco.cxx:676
CbmAnaConversionReco::fRecoRefittedMomentum
std::vector< TVector3 > fRecoRefittedMomentum
Definition CbmAnaConversionReco.h:87
CbmAnaConversionReco::CutEfficiencyStudies
void CutEfficiencyStudies(int e1, int e2, int e3, int e4, int motherE1, int motherE2, int motherE3, int motherE4, int IsEta=0)
Definition CbmAnaConversionReco.cxx:1879
CbmAnaConversionReco::CalculateOpeningAngleReco
Double_t CalculateOpeningAngleReco(TVector3 electron1, TVector3 electron2)
Definition CbmAnaConversionReco.cxx:2219
CbmAnaConversionReco::CalculateOpeningAngleBetweenGammasMC
Double_t CalculateOpeningAngleBetweenGammasMC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2, CbmMCTrack *mctrack3, CbmMCTrack *mctrack4)
Definition CbmAnaConversionReco.cxx:2253
CbmAnaConversionReco::fRecoTracklistEPEM
std::vector< CbmMCTrack * > fRecoTracklistEPEM
Definition CbmAnaConversionReco.h:82
CbmAnaConversionReco::fRecoTracklistEPEM_gtid
std::vector< Int_t > fRecoTracklistEPEM_gtid
Definition CbmAnaConversionReco.h:85
CbmAnaConversionReco::CalculateOpeningAngleMC
Double_t CalculateOpeningAngleMC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2)
Definition CbmAnaConversionReco.cxx:2234
CbmAnaConversionReco::Invmass_4particles
Double_t Invmass_4particles(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
Definition CbmAnaConversionReco.cxx:932
CbmAnaConversionReco::SmearValue
Double_t SmearValue(Double_t value)
Definition CbmAnaConversionReco.cxx:994
CbmAnaConversionReco::fRecoMomentum
std::vector< TVector3 > fRecoMomentum
Definition CbmAnaConversionReco.h:86
CbmAnaConversionReco::NofDaughters
Int_t NofDaughters(Int_t motherId)
Definition CbmAnaConversionReco.cxx:2206
CbmAnaConversionReco::fRecoTracklistEPEM_chi
std::vector< Double_t > fRecoTracklistEPEM_chi
Definition CbmAnaConversionReco.h:84
CbmAnaConversionReco::fHistoList_reco_mom
std::vector< TH1 * > fHistoList_reco_mom
Definition CbmAnaConversionReco.h:92
CbmAnaConversionReco::CalculateOpeningAngleBetweenGammasReco
Double_t CalculateOpeningAngleBetweenGammasReco(TVector3 electron1, TVector3 electron2, TVector3 electron3, TVector3 electron4)
Definition CbmAnaConversionReco.cxx:2286
CbmAnaConversionReco::CalculateKinematicParamsReco
LmvmKinePar CalculateKinematicParamsReco(const TVector3 electron1, const TVector3 electron2)
Definition CbmAnaConversionReco.cxx:2376
CbmAnaConversionReco::fHistoList_eeee
std::vector< TH1 * > fHistoList_eeee
Definition CbmAnaConversionReco.h:95
CbmAnaConversionReco::IsRichElectronANN
Bool_t IsRichElectronANN(Int_t globalTrackIndex, Double_t momentum)
Definition CbmAnaConversionReco.cxx:2442
CbmAnaConversionReco::fHistoList_reco
std::vector< TH1 * > fHistoList_reco
Definition CbmAnaConversionReco.h:91
CbmAnaConversionReco::fRecoTracklistEPEM_ids
std::vector< int > fRecoTracklistEPEM_ids
Definition CbmAnaConversionReco.h:83
CbmAnaConversionReco::fMCTracklist_all
std::vector< CbmMCTrack * > fMCTracklist_all
Definition CbmAnaConversionReco.h:81
CbmAnaConversionReco::SetTracklistReco
void SetTracklistReco(std::vector< CbmMCTrack * > MCTracklist, std::vector< TVector3 > RecoTracklist1, std::vector< TVector3 > RecoTracklist2, std::vector< int > ids, std::vector< Double_t > chi, std::vector< Int_t > GlobalTrackId)
Definition CbmAnaConversionReco.cxx:679
CbmAnaConversionReco::fHistoList_MC
std::vector< TH1 * > fHistoList_MC
Definition CbmAnaConversionReco.h:90
CbmAnaConversionReco::fHistoList_gg
std::vector< TH1 * > fHistoList_gg
Definition CbmAnaConversionReco.h:93
CbmAnaConversionReco::Invmass_4particlesRECO
Double_t Invmass_4particlesRECO(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
Definition CbmAnaConversionReco.cxx:1013
CbmAnaConversionReco::fHistoList_all
std::vector< TH1 * > fHistoList_all
Definition CbmAnaConversionReco.h:96
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:170
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetMass
double GetMass() const
Mass of the associated particle.
Definition CbmMCTrack.cxx:115
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMCTrack::Get4Momentum
void Get4Momentum(TLorentzVector &momentum) const
Definition CbmMCTrack.h:173
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
CbmRichElectronIdAnn::CalculateAnnValue
double CalculateAnnValue(int globalTrackIndex, double momentum)
Calculate output value of the ANN.
Definition CbmRichElectronIdAnn.cxx:85
CbmRichElectronIdAnn::GetInstance
static CbmRichElectronIdAnn & GetInstance()
Definition CbmRichElectronIdAnn.h:47
CbmRichRing::GetHit
uint32_t GetHit(int32_t i) const
Definition CbmRichRing.h:39
CbmRichRing::GetRadius
float GetRadius() const
Definition CbmRichRing.h:79
CbmRichRing::GetAaxis
double GetAaxis() const
Definition CbmRichRing.h:80
CbmRichRing::GetNofHits
int32_t GetNofHits() const
Definition CbmRichRing.h:37
CbmRichRing::GetCenterX
float GetCenterX() const
Definition CbmRichRing.h:77
CbmRichRing::GetCenterY
float GetCenterY() const
Definition CbmRichRing.h:78
CbmRichRing::GetSelectionNN
double GetSelectionNN() const
Definition CbmRichRing.h:91
CbmRichRing::GetBaxis
double GetBaxis() const
Definition CbmRichRing.h:81
CbmRichUtil::GetRingTrackDistance
static double GetRingTrackDistance(int globalTrackId)
Definition CbmRichUtil.cxx:24
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
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21