CbmRoot
Loading...
Searching...
No Matches
CbmAnaConversionTest2.cxx
Go to the documentation of this file.
1/* Copyright (C) 2016-2019 Fakultaet fuer Mathematik und Naturwissenschaften, Bergische Universitaet Wuppertal, Wuppertal
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Sascha Reinecke [committer], Florian Uhlig */
4
12#include "CbmAnaConversionTest2.h"
13
14#include "CbmAnaConversionCutSettings.h"
15#include "CbmAnaConversionKinematicParams.h"
16#include "CbmGlobalTrack.h"
17#include "CbmL1PFFitter.h"
18#include "CbmLitGlobalElectronId.h"
19#include "CbmMCTrack.h"
20#include "CbmRichPoint.h"
21#include "CbmRichRing.h"
22#include "CbmStsTrack.h"
23#include "CbmTrackMatchNew.h"
24
25#include "FairRootManager.h"
26#include <Logger.h>
27
28#include <algorithm>
29#include <map>
30
31
32using namespace std;
33
34
35CbmAnaConversionTest2::CbmAnaConversionTest2()
36 : fRichPoints(NULL)
37 , fRichRings(NULL)
38 , fRichRingMatches(NULL)
39 , fMcTracks(NULL)
40 , fStsTracks(NULL)
41 , fStsTrackMatches(NULL)
42 , fGlobalTracks(NULL)
43 , fPrimVertex(NULL)
44 , fKFVertex()
45 , fHistoList_test2()
46 , fhTest2_invmass_RICHindex0(NULL)
47 , fhTest2_invmass_RICHindex1(NULL)
48 , fhTest2_invmass_RICHindex2(NULL)
49 , fhTest2_invmass_RICHindex3(NULL)
50 , fhTest2_invmass_RICHindex4(NULL)
51 , fVector_gt()
52 , fVector_momenta()
53 , fVector_mctrack()
54 , fVector_gtIndex()
55 , fVector_richIndex()
56 , fhTest2_invmass_gee_mc(NULL)
57 , fhTest2_invmass_gee_refitted(NULL)
58 , fhTest2_invmass_gg_mc(NULL)
59 , fhTest2_invmass_gg_refitted(NULL)
60 , fhTest2_invmass_all_mc(NULL)
61 , fhTest2_invmass_all_refitted(NULL)
62 , fhTest2_pt_vs_rap_gee(NULL)
63 , fhTest2_pt_vs_rap_gg(NULL)
64 , fhTest2_pt_vs_rap_all(NULL)
65 , fhTest2_startvertexElectrons_gee(NULL)
66 , fhTest2_startvertexElectrons_gg(NULL)
67 , fhTest2_startvertexElectrons_all(NULL)
68 , fhTest2_2rich_invmass_gee_mc(NULL)
69 , fhTest2_2rich_invmass_gee_refitted(NULL)
70 , fhTest2_2rich_invmass_gg_mc(NULL)
71 , fhTest2_2rich_invmass_gg_refitted(NULL)
72 , fhTest2_2rich_invmass_all_mc(NULL)
73 , fhTest2_2rich_invmass_all_refitted(NULL)
74 , fhTest2_2rich_pt_vs_rap_gee(NULL)
75 , fhTest2_2rich_pt_vs_rap_gg(NULL)
76 , fhTest2_2rich_pt_vs_rap_all(NULL)
77 , fhTest2_electrons_pt_vs_p(NULL)
78 , fhTest2_electrons_pt_vs_p_withRICH(NULL)
79 , fhTest2_electrons_pt_vs_p_noRICH(NULL)
80 , fhTest2_3rich_electrons_theta_included(NULL)
81 , fhTest2_3rich_electrons_theta_missing(NULL)
82 , fhTest2_3rich_electrons_thetaVSp_included(NULL)
83 , fhTest2_3rich_electrons_thetaVSp_missing(NULL)
84{
85}
86
87CbmAnaConversionTest2::~CbmAnaConversionTest2() {}
88
89
90void CbmAnaConversionTest2::Init()
91{
92 FairRootManager* ioman = FairRootManager::Instance();
93 if (NULL == ioman) { Fatal("CbmAnaConversion::Init", "RootManager not instantised!"); }
94
95 fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
96 if (NULL == fRichPoints) { Fatal("CbmAnaConversion::Init", "No RichPoint array!"); }
97
98 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
99 if (NULL == fMcTracks) { Fatal("CbmAnaConversion::Init", "No MCTrack array!"); }
100
101 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
102 if (NULL == fStsTracks) { Fatal("CbmAnaConversion::Init", "No StsTrack array!"); }
103
104 fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
105 if (NULL == fStsTrackMatches) { Fatal("CbmAnaConversion::Init", "No StsTrackMatch array!"); }
106
107 fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
108 if (NULL == fGlobalTracks) { Fatal("CbmAnaConversion::Init", "No GlobalTrack array!"); }
109
110 // Get pointer to PrimaryVertex object from IOManager if it exists
111 // The old name for the object is "PrimaryVertex" the new one
112 // "PrimaryVertex." Check first for the new name
113 fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex."));
114 if (nullptr == fPrimVertex) { fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex")); }
115 if (nullptr == fPrimVertex) { LOG(fatal) << "No PrimaryVertex array!"; }
116
117 fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
118 if (NULL == fRichRings) { Fatal("CbmAnaConversion::Init", "No RichRing array!"); }
119
120 fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
121 if (NULL == fRichRingMatches) { Fatal("CbmAnaConversion::Init", "No RichRingMatch array!"); }
122
123 InitHistos();
124}
125
126
127void CbmAnaConversionTest2::InitHistos()
128{
129 fHistoList_test2.clear();
130
131 Double_t invmassSpectra_nof = 800;
132 Double_t invmassSpectra_start = -0.00125;
133 Double_t invmassSpectra_end = 1.99875;
134
135
136 fhTest2_invmass_RICHindex0 =
137 new TH1D("fhTest2_invmass_RICHindex0", "fhTest2_invmass_RICHindex0; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
138 600, -0.0025, 2.9975);
139 fHistoList_test2.push_back(fhTest2_invmass_RICHindex0);
140 fhTest2_invmass_RICHindex1 =
141 new TH1D("fhTest2_invmass_RICHindex1", "fhTest2_invmass_RICHindex1; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
142 600, -0.0025, 2.9975);
143 fHistoList_test2.push_back(fhTest2_invmass_RICHindex1);
144 fhTest2_invmass_RICHindex2 =
145 new TH1D("fhTest2_invmass_RICHindex2", "fhTest2_invmass_RICHindex2; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
146 600, -0.0025, 2.9975);
147 fHistoList_test2.push_back(fhTest2_invmass_RICHindex2);
148 fhTest2_invmass_RICHindex3 =
149 new TH1D("fhTest2_invmass_RICHindex3", "fhTest2_invmass_RICHindex3; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
150 600, -0.0025, 2.9975);
151 fHistoList_test2.push_back(fhTest2_invmass_RICHindex3);
152 fhTest2_invmass_RICHindex4 =
153 new TH1D("fhTest2_invmass_RICHindex4", "fhTest2_invmass_RICHindex4; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
154 600, -0.0025, 2.9975);
155 fHistoList_test2.push_back(fhTest2_invmass_RICHindex4);
156
157
158 fhTest2_invmass_gee_mc =
159 new TH1D("fhTest2_invmass_gee_mc", "fhTest2_invmass_gee_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
160 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
161 fhTest2_invmass_gee_refitted =
162 new TH1D("fhTest2_invmass_gee_refitted", "fhTest2_invmass_gee_refitted;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
163 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
164 fhTest2_invmass_gg_mc =
165 new TH1D("fhTest2_invmass_gg_mc", "fhTest2_invmass_gg_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
166 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
167 fhTest2_invmass_gg_refitted =
168 new TH1D("fhTest2_invmass_gg_refitted", "fhTest2_invmass_gg_refitted;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
169 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
170 fhTest2_invmass_all_mc =
171 new TH1D("fhTest2_invmass_all_mc", "fhTest2_invmass_all_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
172 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
173 fhTest2_invmass_all_refitted =
174 new TH1D("fhTest2_invmass_all_refitted", "fhTest2_invmass_all_refitted;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
175 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
176 fHistoList_test2.push_back(fhTest2_invmass_gee_mc);
177 fHistoList_test2.push_back(fhTest2_invmass_gee_refitted);
178 fHistoList_test2.push_back(fhTest2_invmass_gg_mc);
179 fHistoList_test2.push_back(fhTest2_invmass_gg_refitted);
180 fHistoList_test2.push_back(fhTest2_invmass_all_mc);
181 fHistoList_test2.push_back(fhTest2_invmass_all_refitted);
182
183 fhTest2_pt_vs_rap_gee =
184 new TH2D("fhTest2_pt_vs_rap_gee", "fhTest2_pt_vs_rap_gee;p_{t} in GeV/c; rapidity y", 240, -2., 10., 270, -2., 7.);
185 fhTest2_pt_vs_rap_gg =
186 new TH2D("fhTest2_pt_vs_rap_gg", "fhTest2_pt_vs_rap_gg;p_{t} in GeV/c; rapidity y", 240, -2., 10., 270, -2., 7.);
187 fhTest2_pt_vs_rap_all =
188 new TH2D("fhTest2_pt_vs_rap_all", "fhTest2_pt_vs_rap_all;p_{t} in GeV/c; rapidity y", 240, -2., 10., 270, -2., 7.);
189 fHistoList_test2.push_back(fhTest2_pt_vs_rap_gee);
190 fHistoList_test2.push_back(fhTest2_pt_vs_rap_gg);
191 fHistoList_test2.push_back(fhTest2_pt_vs_rap_all);
192
193 fhTest2_startvertexElectrons_gee =
194 new TH1D("fhTest2_startvertexElectrons_gee", "fhTest2_startvertexElectrons_gee;z in cm;#", 411, -5.25, 200.25);
195 fhTest2_startvertexElectrons_gg =
196 new TH1D("fhTest2_startvertexElectrons_gg", "fhTest2_startvertexElectrons_gg;z in cm;#", 411, -5.25, 200.25);
197 fhTest2_startvertexElectrons_all =
198 new TH1D("fhTest2_startvertexElectrons_all", "fhTest2_startvertexElectrons_all;z in cm;#", 411, -5.25, 200.25);
199 fHistoList_test2.push_back(fhTest2_startvertexElectrons_gee);
200 fHistoList_test2.push_back(fhTest2_startvertexElectrons_gg);
201 fHistoList_test2.push_back(fhTest2_startvertexElectrons_all);
202
203
204 // 2 leptons in RICH
205 fhTest2_2rich_invmass_gee_mc =
206 new TH1D("fhTest2_2rich_invmass_gee_mc", "fhTest2_2rich_invmass_gee_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
207 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
208 fhTest2_2rich_invmass_gee_refitted = new TH1D("fhTest2_2rich_invmass_gee_refitted",
209 "fhTest2_2rich_invmass_gee_refitted;invariant mass of 4 e^{#pm} "
210 "in GeV/c^{2};#",
211 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
212 fhTest2_2rich_invmass_gg_mc =
213 new TH1D("fhTest2_2rich_invmass_gg_mc", "fhTest2_2rich_invmass_gg_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
214 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
215 fhTest2_2rich_invmass_gg_refitted = new TH1D("fhTest2_2rich_invmass_gg_refitted",
216 "fhTest2_2rich_invmass_gg_refitted;invariant mass of 4 e^{#pm} in "
217 "GeV/c^{2};#",
218 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
219 fhTest2_2rich_invmass_all_mc =
220 new TH1D("fhTest2_2rich_invmass_all_mc", "fhTest2_2rich_invmass_all_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
221 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
222 fhTest2_2rich_invmass_all_refitted = new TH1D("fhTest2_2rich_invmass_all_refitted",
223 "fhTest2_2rich_invmass_all_refitted;invariant mass of 4 e^{#pm} "
224 "in GeV/c^{2};#",
225 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
226 fHistoList_test2.push_back(fhTest2_2rich_invmass_gee_mc);
227 fHistoList_test2.push_back(fhTest2_2rich_invmass_gee_refitted);
228 fHistoList_test2.push_back(fhTest2_2rich_invmass_gg_mc);
229 fHistoList_test2.push_back(fhTest2_2rich_invmass_gg_refitted);
230 fHistoList_test2.push_back(fhTest2_2rich_invmass_all_mc);
231 fHistoList_test2.push_back(fhTest2_2rich_invmass_all_refitted);
232
233
234 fhTest2_2rich_pt_vs_rap_gee =
235 new TH2D("fhTest2_2rich_pt_vs_rap_gee", "fhTest2_2rich_pt_vs_rap_gee;p_{t} in GeV/c; rapidity y", 240, -2., 10.,
236 270, -2., 7.);
237 fhTest2_2rich_pt_vs_rap_gg = new TH2D(
238 "fhTest2_2rich_pt_vs_rap_gg", "fhTest2_2rich_pt_vs_rap_gg;p_{t} in GeV/c; rapidity y", 240, -2., 10., 270, -2., 7.);
239 fhTest2_2rich_pt_vs_rap_all =
240 new TH2D("fhTest2_2rich_pt_vs_rap_all", "fhTest2_2rich_pt_vs_rap_all;p_{t} in GeV/c; rapidity y", 240, -2., 10.,
241 270, -2., 7.);
242 fHistoList_test2.push_back(fhTest2_2rich_pt_vs_rap_gee);
243 fHistoList_test2.push_back(fhTest2_2rich_pt_vs_rap_gg);
244 fHistoList_test2.push_back(fhTest2_2rich_pt_vs_rap_all);
245
246
247 // further tests
248 fhTest2_electrons_pt_vs_p = new TH2D(
249 "fhTest2_electrons_pt_vs_p", "fhTest2_electrons_pt_vs_p;p_{t} in GeV/c; p in GeV/c", 240, -2., 10., 360, -2., 16.);
250 fHistoList_test2.push_back(fhTest2_electrons_pt_vs_p);
251 fhTest2_electrons_pt_vs_p_withRICH =
252 new TH2D("fhTest2_electrons_pt_vs_p_withRICH", "fhTest2_electrons_pt_vs_p_withRICH;p_{t} in GeV/c; p in GeV/c", 240,
253 -2., 10., 360, -2., 16.);
254 fHistoList_test2.push_back(fhTest2_electrons_pt_vs_p_withRICH);
255 fhTest2_electrons_pt_vs_p_noRICH =
256 new TH2D("fhTest2_electrons_pt_vs_p_noRICH", "fhTest2_electrons_pt_vs_p_noRICH;p_{t} in GeV/c; p in GeV/c", 240,
257 -2., 10., 360, -2., 16.);
258 fHistoList_test2.push_back(fhTest2_electrons_pt_vs_p_noRICH);
259
260 fhTest2_3rich_electrons_theta_included = new TH1D(
261 "fhTest2_3rich_electrons_theta_included", "fhTest2_3rich_electrons_theta_included;theta angle [deg];#", 90, 0, 90);
262 fHistoList_test2.push_back(fhTest2_3rich_electrons_theta_included);
263 fhTest2_3rich_electrons_theta_missing = new TH1D(
264 "fhTest2_3rich_electrons_theta_missing", "fhTest2_3rich_electrons_theta_missing;theta angle [deg];#", 90, 0, 90);
265 fHistoList_test2.push_back(fhTest2_3rich_electrons_theta_missing);
266 fhTest2_3rich_electrons_thetaVSp_included =
267 new TH2D("fhTest2_3rich_electrons_thetaVSp_included",
268 "fhTest2_3rich_electrons_thetaVSp_included;theta angle [deg];p in GeV/c", 90, 0, 90, 540, -2., 16.);
269 fHistoList_test2.push_back(fhTest2_3rich_electrons_thetaVSp_included);
270 fhTest2_3rich_electrons_thetaVSp_missing =
271 new TH2D("fhTest2_3rich_electrons_thetaVSp_missing",
272 "fhTest2_3rich_electrons_thetaVSp_missing;theta angle [deg];p in GeV/c", 90, 0, 90, 540, -2., 16.);
273 fHistoList_test2.push_back(fhTest2_3rich_electrons_thetaVSp_missing);
274}
275
276
277void CbmAnaConversionTest2::Finish()
278{
279 //gDirectory->cd("analysis-conversion");
280 gDirectory->mkdir("Test2");
281 gDirectory->cd("Test2");
282 for (UInt_t i = 0; i < fHistoList_test2.size(); i++) {
283 fHistoList_test2[i]->Write();
284 }
285 gDirectory->cd("..");
286}
287
288
289void CbmAnaConversionTest2::Exec()
290{
291 fVector_gt.clear();
292 fVector_momenta.clear();
293 fVector_mctrack.clear();
294 fVector_gtIndex.clear();
295 fVector_richIndex.clear();
296
297
298 Int_t ngTracks = fGlobalTracks->GetEntriesFast();
299 for (Int_t i = 0; i < ngTracks; i++) {
300 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
301 if (NULL == gTrack) continue;
302 int stsInd = gTrack->GetStsTrackIndex();
303 int richInd = gTrack->GetRichRingIndex();
304
305 if (stsInd < 0) continue;
306 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
307 if (stsTrack == NULL) continue;
308 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
309 if (stsMatch == NULL) continue;
310 if (stsMatch->GetNofLinks() <= 0) continue;
311 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
312 if (stsMcTrackId < 0) continue;
313 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
314 if (mcTrack1 == NULL) continue;
315
316
317 // calculate refitted momenta at primary vertex
318 TVector3 refittedMomentum;
319 CbmL1PFFitter fPFFitter;
320 vector<CbmStsTrack> stsTracks;
321 stsTracks.resize(1);
322 stsTracks[0] = *stsTrack;
323 vector<CbmL1PFFitter::PFFieldRegion> vField;
324 vector<float> chiPrim;
325 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
326 //cand.chi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
327 //cand.chi2Prim = chiPrim[0];
328 const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();
329 vtxTrack->Momentum(refittedMomentum);
330 //float result_chi = chiPrim[0];
331
332
333 // Doing refit of momenta with electron assumption
334 CbmL1PFFitter fPFFitter_electron;
335 vector<CbmStsTrack> stsTracks_electron;
336 stsTracks_electron.resize(1);
337 stsTracks_electron[0] = *stsTrack;
338 vector<CbmL1PFFitter::PFFieldRegion> vField_electron;
339 vector<float> chiPrim_electron;
340 vector<int> pidHypo_electron;
341 pidHypo_electron.push_back(11);
342 fPFFitter_electron.Fit(stsTracks_electron, pidHypo_electron);
343 fPFFitter_electron.GetChiToVertex(stsTracks_electron, vField_electron, chiPrim_electron, fKFVertex, 3e6);
344
345 TVector3 refittedMomentum_electron;
346 const FairTrackParam* vtxTrack_electron = stsTracks_electron[0].GetParamFirst();
347 vtxTrack_electron->Momentum(refittedMomentum_electron);
348 //float result_chi_electron = chiPrim_electron[0];
349 //float result_ndf_electron = stsTracks_electron[0].GetNDF();
350 //Double_t stsHits = stsTrack->GetNofHits();
351
352
353 Int_t pdg = mcTrack1->GetPdgCode();
354
355 if (TMath::Abs(pdg) == 11) {
356 fhTest2_electrons_pt_vs_p->Fill(refittedMomentum_electron.Perp(), refittedMomentum_electron.Mag());
357 if (richInd > 0) {
358 fhTest2_electrons_pt_vs_p_withRICH->Fill(refittedMomentum_electron.Perp(), refittedMomentum_electron.Mag());
359 }
360 if (richInd < 0) {
361 fhTest2_electrons_pt_vs_p_noRICH->Fill(refittedMomentum_electron.Perp(), refittedMomentum_electron.Mag());
362 }
363
364
365 fVector_momenta.push_back(refittedMomentum_electron);
366 fVector_mctrack.push_back(mcTrack1);
367 fVector_gtIndex.push_back(i);
368 fVector_gt.push_back(gTrack);
369
370 if (richInd < 0) {
371 fVector_richIndex.push_back(0);
372 continue;
373 }
374 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
375 if (richMatch == NULL) {
376 fVector_richIndex.push_back(0);
377 continue;
378 }
379 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
380 if (richMcTrackId < 0) {
381 fVector_richIndex.push_back(0);
382 continue;
383 }
384 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
385 if (mcTrack2 == NULL) {
386 fVector_richIndex.push_back(0);
387 continue;
388 }
389
390 fVector_richIndex.push_back(richInd);
391 }
392 }
393
394 InvariantMassTest_3RICH();
395 InvariantMassTest_2RICH();
396}
397
398
399void CbmAnaConversionTest2::InvariantMassTest_3RICH()
400// Calculating invariant mass of 4 ep/em, using MC data AND reconstructed momentum
401{
402 cout << "CbmAnaConversionTest2: InvariantMassTest_3RICH - Start..." << endl;
403 cout << "CbmAnaConversionTest2: InvariantMassTest_3RICH - " << fVector_mctrack.size() << endl;
404 int fill = 0;
405 if (fVector_mctrack.size() < 4) return;
406 for (unsigned int i = 0; i < fVector_mctrack.size() - 3; i++) {
407 if (i % 10 == 0) cout << "CbmAnaConversionTest2: InvariantMassTest_3RICH - iteration i = " << i << endl;
408 for (unsigned int j = i + 1; j < fVector_mctrack.size() - 2; j++) {
409 for (unsigned int k = j + 1; k < fVector_mctrack.size() - 1; k++) {
410 for (unsigned int l = k + 1; l < fVector_mctrack.size(); l++) {
411 if (fVector_mctrack[i]->GetPdgCode() + fVector_mctrack[j]->GetPdgCode() + fVector_mctrack[k]->GetPdgCode()
412 + fVector_mctrack[l]->GetPdgCode()
413 != 0)
414 continue;
415
416 if (fVector_mctrack.size() != fVector_richIndex.size()) {
417 cout << "CbmAnaConversionTest2: InvariantMassTest_4epem - not "
418 "matching number of entries!"
419 << endl;
420 continue;
421 }
422
423 Bool_t IsPi0 = false;
424
425 // starting points of each electron (-> i.e. conversion points of gamma OR decay points of pi0, depending on decay channel)
426 TVector3 pi0start_i;
427 fVector_mctrack[i]->GetStartVertex(pi0start_i);
428 TVector3 pi0start_j;
429 fVector_mctrack[j]->GetStartVertex(pi0start_j);
430 TVector3 pi0start_k;
431 fVector_mctrack[k]->GetStartVertex(pi0start_k);
432 TVector3 pi0start_l;
433 fVector_mctrack[l]->GetStartVertex(pi0start_l);
434
435
436 Int_t richIndex1 = (fVector_richIndex[i] > 0);
437 Int_t richIndex2 = (fVector_richIndex[j] > 0);
438 Int_t richIndex3 = (fVector_richIndex[k] > 0);
439 Int_t richIndex4 = (fVector_richIndex[l] > 0);
440
441 if (richIndex1 + richIndex2 + richIndex3 + richIndex4 != 3) continue;
442
443
444 int motherId1 = fVector_mctrack[i]->GetMotherId();
445 int motherId2 = fVector_mctrack[j]->GetMotherId();
446 int motherId3 = fVector_mctrack[k]->GetMotherId();
447 int motherId4 = fVector_mctrack[l]->GetMotherId();
448
449
450 if ((motherId1 == motherId2 && motherId3 == motherId4) || (motherId1 == motherId3 && motherId2 == motherId4)
451 || (motherId1 == motherId4 && motherId2 == motherId3)) { // start 1
452
453
454 int grandmotherId1 = -1;
455 int grandmotherId2 = -1;
456 int grandmotherId3 = -1;
457 int grandmotherId4 = -1;
458
459 int mcMotherPdg1 = -1;
460 // int mcMotherPdg2 = -1; (FU) not used
461 // int mcMotherPdg3 = -1; (FU) not used
462 // int mcMotherPdg4 = -1; (FU) not used
463 int mcGrandmotherPdg1 = -1;
464 // int mcGrandmotherPdg2 = -1;
465 // int mcGrandmotherPdg3 = -1;
466 // int mcGrandmotherPdg4 = -1;
467
468 CbmMCTrack* grandmother1 = nullptr;
469
470 if (motherId1 != -1) {
471 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
472 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
473 grandmotherId1 = mother1->GetMotherId();
474 if (grandmotherId1 != -1) {
475 grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
476 if (NULL != grandmother1) mcGrandmotherPdg1 = grandmother1->GetPdgCode();
477 }
478 }
479 if (motherId2 != -1) {
480 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
481 // if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode(); (FU) not used
482 grandmotherId2 = mother2->GetMotherId();
483 if (grandmotherId2 != -1) {
484 // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
485 // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
486 }
487 }
488 if (motherId3 != -1) {
489 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
490 // if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode(); (FU) not used
491 grandmotherId3 = mother3->GetMotherId();
492 if (grandmotherId3 != -1) {
493 // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
494 // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
495 }
496 }
497 if (motherId4 != -1) {
498 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
499 // if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode(); (FU) not used
500 grandmotherId4 = mother4->GetMotherId();
501 if (grandmotherId4 != -1) {
502 // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
503 // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
504 }
505 }
506
507 if (grandmotherId1 == -1) continue;
508
509 if (motherId1 == motherId2 && motherId3 == motherId4) {
510 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId3) != 2) continue;
511 if ((grandmotherId1 == motherId3 && mcGrandmotherPdg1 == 111)
512 || (motherId1 == grandmotherId3 && mcMotherPdg1 == 111)) {
513
514 fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
515 fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
516 fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
517 fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
518 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
519 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
520 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
521 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
522 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
523 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
524
525
526 Double_t invmass_mc = 0;
527 Double_t invmass_reco = 0;
528 CbmAnaConversionKinematicParams params_mc =
529 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(
530 fVector_mctrack[i], fVector_mctrack[j], fVector_mctrack[k], fVector_mctrack[l]);
531 CbmAnaConversionKinematicParams params_reco =
532 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
533 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
534 invmass_mc = params_mc.fMinv;
535 invmass_reco = params_reco.fMinv;
536
537 fhTest2_invmass_gee_mc->Fill(invmass_mc);
538 fhTest2_invmass_gee_refitted->Fill(invmass_reco);
539 fhTest2_invmass_all_mc->Fill(invmass_mc);
540 fhTest2_invmass_all_refitted->Fill(invmass_reco);
541
542 fhTest2_pt_vs_rap_gee->Fill(params_reco.fPt, params_reco.fRapidity);
543 fhTest2_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
544
545 IsPi0 = true;
546 }
547 }
548 if (motherId1 == motherId3 && motherId2 == motherId4) {
549 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2) continue;
550 if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
551 || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
552
553 fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
554 fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
555 fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
556 fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
557 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
558 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
559 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
560 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
561 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
562 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
563
564
565 Double_t invmass_mc = 0;
566 Double_t invmass_reco = 0;
567 CbmAnaConversionKinematicParams params_mc =
568 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(
569 fVector_mctrack[i], fVector_mctrack[j], fVector_mctrack[k], fVector_mctrack[l]);
570 CbmAnaConversionKinematicParams params_reco =
571 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
572 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
573 invmass_mc = params_mc.fMinv;
574 invmass_reco = params_reco.fMinv;
575
576 fhTest2_invmass_gee_mc->Fill(invmass_mc);
577 fhTest2_invmass_gee_refitted->Fill(invmass_reco);
578 fhTest2_invmass_all_mc->Fill(invmass_mc);
579 fhTest2_invmass_all_refitted->Fill(invmass_reco);
580
581 fhTest2_pt_vs_rap_gee->Fill(params_reco.fPt, params_reco.fRapidity);
582 fhTest2_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
583
584 IsPi0 = true;
585 }
586 }
587 if (motherId1 == motherId4 && motherId2 == motherId3) {
588 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2) continue;
589 if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
590 || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
591
592 fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
593 fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
594 fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
595 fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
596 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
597 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
598 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
599 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
600 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
601 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
602
603
604 Double_t invmass_mc = 0;
605 Double_t invmass_reco = 0;
606 CbmAnaConversionKinematicParams params_mc =
607 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(
608 fVector_mctrack[i], fVector_mctrack[j], fVector_mctrack[k], fVector_mctrack[l]);
609 CbmAnaConversionKinematicParams params_reco =
610 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
611 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
612 invmass_mc = params_mc.fMinv;
613 invmass_reco = params_reco.fMinv;
614
615 fhTest2_invmass_gee_mc->Fill(invmass_mc);
616 fhTest2_invmass_gee_refitted->Fill(invmass_reco);
617 fhTest2_invmass_all_mc->Fill(invmass_mc);
618 fhTest2_invmass_all_refitted->Fill(invmass_reco);
619
620 fhTest2_pt_vs_rap_gee->Fill(params_reco.fPt, params_reco.fRapidity);
621 fhTest2_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
622
623 IsPi0 = true;
624 }
625 }
626
627
628 // ===================================================================================================
629 // HERE DECAY pi0 -> gamma gamma -> e+e- e+e-
630 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
631 && grandmotherId1 == grandmotherId4) {
632 if (mcGrandmotherPdg1 != 111) continue; // 111 = pi0, 221 = eta
633
634 TVector3 pi0start;
635 grandmother1->GetStartVertex(pi0start);
636
637 fhTest2_startvertexElectrons_gg->Fill(pi0start_i.Z());
638 fhTest2_startvertexElectrons_gg->Fill(pi0start_j.Z());
639 fhTest2_startvertexElectrons_gg->Fill(pi0start_k.Z());
640 fhTest2_startvertexElectrons_gg->Fill(pi0start_l.Z());
641 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
642 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
643 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
644 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
645
646
647 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
648 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
649
650 Double_t invmass_mc = 0;
651 Double_t invmass_reco = 0;
652 CbmAnaConversionKinematicParams params_mc =
653 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(fVector_mctrack[i], fVector_mctrack[j],
654 fVector_mctrack[k], fVector_mctrack[l]);
655 CbmAnaConversionKinematicParams params_reco =
656 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
657 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
658 invmass_mc = params_mc.fMinv;
659 invmass_reco = params_reco.fMinv;
660
661 fhTest2_invmass_gg_mc->Fill(invmass_mc);
662 fhTest2_invmass_gg_refitted->Fill(invmass_reco);
663 fhTest2_invmass_all_mc->Fill(invmass_mc);
664 fhTest2_invmass_all_refitted->Fill(invmass_reco);
665
666 fhTest2_pt_vs_rap_gg->Fill(params_reco.fPt, params_reco.fRapidity);
667 fhTest2_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
668
669 IsPi0 = true;
670
671 cout << "########################################################"
672 "##############"
673 << endl;
674
675 fill++;
676 }
677
678
679 if (richIndex1 <= 0 && IsPi0) {
680 fhTest2_3rich_electrons_theta_missing->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi());
681 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi());
682 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi());
683 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi());
684 fhTest2_3rich_electrons_thetaVSp_missing->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi(),
685 fVector_momenta[i].Mag());
686 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi(),
687 fVector_momenta[j].Mag());
688 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi(),
689 fVector_momenta[k].Mag());
690 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi(),
691 fVector_momenta[l].Mag());
692 }
693 if (richIndex2 <= 0 && IsPi0) {
694 fhTest2_3rich_electrons_theta_missing->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi());
695 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi());
696 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi());
697 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi());
698 fhTest2_3rich_electrons_thetaVSp_missing->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi(),
699 fVector_momenta[j].Mag());
700 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi(),
701 fVector_momenta[i].Mag());
702 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi(),
703 fVector_momenta[k].Mag());
704 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi(),
705 fVector_momenta[l].Mag());
706 }
707 if (richIndex3 <= 0 && IsPi0) {
708 fhTest2_3rich_electrons_theta_missing->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi());
709 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi());
710 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi());
711 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi());
712 fhTest2_3rich_electrons_thetaVSp_missing->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi(),
713 fVector_momenta[k].Mag());
714 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi(),
715 fVector_momenta[i].Mag());
716 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi(),
717 fVector_momenta[j].Mag());
718 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi(),
719 fVector_momenta[l].Mag());
720 }
721 if (richIndex4 <= 0 && IsPi0) {
722 fhTest2_3rich_electrons_theta_missing->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi());
723 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi());
724 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi());
725 fhTest2_3rich_electrons_theta_included->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi());
726 fhTest2_3rich_electrons_thetaVSp_missing->Fill(fVector_momenta[l].Theta() * 180 / TMath::Pi(),
727 fVector_momenta[l].Mag());
728 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[i].Theta() * 180 / TMath::Pi(),
729 fVector_momenta[i].Mag());
730 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[j].Theta() * 180 / TMath::Pi(),
731 fVector_momenta[j].Mag());
732 fhTest2_3rich_electrons_thetaVSp_included->Fill(fVector_momenta[k].Theta() * 180 / TMath::Pi(),
733 fVector_momenta[k].Mag());
734 }
735
736 } // end 1
737 }
738 }
739 }
740 }
741}
742
743
744void CbmAnaConversionTest2::InvariantMassTest_2RICH()
745// Calculating invariant mass of 4 ep/em, using MC data AND reconstructed momentum
746{
747 cout << "CbmAnaConversionTest2: InvariantMassTest_2RICH - Start..." << endl;
748 cout << "CbmAnaConversionTest2: InvariantMassTest_2RICH - " << fVector_mctrack.size() << endl;
749 int fill = 0;
750 if (fVector_mctrack.size() < 4) return;
751 for (unsigned int i = 0; i < fVector_mctrack.size() - 3; i++) {
752 if (i % 10 == 0) cout << "CbmAnaConversionTest2: InvariantMassTest_2RICH - iteration i = " << i << endl;
753 for (unsigned int j = i + 1; j < fVector_mctrack.size() - 2; j++) {
754 for (unsigned int k = j + 1; k < fVector_mctrack.size() - 1; k++) {
755 for (unsigned int l = k + 1; l < fVector_mctrack.size(); l++) {
756 if (fVector_mctrack[i]->GetPdgCode() + fVector_mctrack[j]->GetPdgCode() + fVector_mctrack[k]->GetPdgCode()
757 + fVector_mctrack[l]->GetPdgCode()
758 != 0)
759 continue;
760
761 if (fVector_mctrack.size() != fVector_richIndex.size()) {
762 cout << "CbmAnaConversionTest2: InvariantMassTest_4epem - not "
763 "matching number of entries!"
764 << endl;
765 continue;
766 }
767
768
769 // starting points of each electron (-> i.e. conversion points of gamma OR decay points of pi0, depending on decay channel)
770 TVector3 pi0start_i;
771 fVector_mctrack[i]->GetStartVertex(pi0start_i);
772 TVector3 pi0start_j;
773 fVector_mctrack[j]->GetStartVertex(pi0start_j);
774 TVector3 pi0start_k;
775 fVector_mctrack[k]->GetStartVertex(pi0start_k);
776 TVector3 pi0start_l;
777 fVector_mctrack[l]->GetStartVertex(pi0start_l);
778
779
780 Int_t richIndex1 = (fVector_richIndex[i] > 0);
781 Int_t richIndex2 = (fVector_richIndex[j] > 0);
782 Int_t richIndex3 = (fVector_richIndex[k] > 0);
783 Int_t richIndex4 = (fVector_richIndex[l] > 0);
784
785 if (richIndex1 + richIndex2 + richIndex3 + richIndex4 != 2) continue;
786
787 int motherId1 = fVector_mctrack[i]->GetMotherId();
788 int motherId2 = fVector_mctrack[j]->GetMotherId();
789 int motherId3 = fVector_mctrack[k]->GetMotherId();
790 int motherId4 = fVector_mctrack[l]->GetMotherId();
791
792
793 if ((motherId1 == motherId2 && motherId3 == motherId4) || (motherId1 == motherId3 && motherId2 == motherId4)
794 || (motherId1 == motherId4 && motherId2 == motherId3)) {
795
796
797 int grandmotherId1 = -1;
798 int grandmotherId2 = -1;
799 int grandmotherId3 = -1;
800 int grandmotherId4 = -1;
801
802 int mcMotherPdg1 = -1;
803 // int mcMotherPdg2 = -1; (FU) not used
804 // int mcMotherPdg3 = -1; (FU) not used
805 // int mcMotherPdg4 = -1; (FU) not used
806 int mcGrandmotherPdg1 = -1;
807 // int mcGrandmotherPdg2 = -1;
808 // int mcGrandmotherPdg3 = -1;
809 // int mcGrandmotherPdg4 = -1;
810
811 CbmMCTrack* grandmother1 = nullptr;
812
813 if (motherId1 != -1) {
814 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
815 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
816 grandmotherId1 = mother1->GetMotherId();
817 if (grandmotherId1 != -1) {
818 grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
819 if (NULL != grandmother1) mcGrandmotherPdg1 = grandmother1->GetPdgCode();
820 }
821 }
822 if (motherId2 != -1) {
823 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
824 // if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode(); (FU) not used
825 grandmotherId2 = mother2->GetMotherId();
826 if (grandmotherId2 != -1) {
827 // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
828 // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
829 }
830 }
831 if (motherId3 != -1) {
832 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
833 // if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode(); (FU) not used
834 grandmotherId3 = mother3->GetMotherId();
835 if (grandmotherId3 != -1) {
836 // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
837 // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
838 }
839 }
840 if (motherId4 != -1) {
841 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
842 // if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode(); (FU) not used
843 grandmotherId4 = mother4->GetMotherId();
844 if (grandmotherId4 != -1) {
845 // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
846 // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
847 }
848 }
849
850 if (grandmotherId1 == -1) continue;
851
852 if (motherId1 == motherId2 && motherId3 == motherId4) {
853 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId3) != 2) continue;
854 if ((grandmotherId1 == motherId3 && mcGrandmotherPdg1 == 111)
855 || (motherId1 == grandmotherId3 && mcMotherPdg1 == 111)) {
856 /*
857 fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
858 fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
859 fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
860 fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
861 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
862 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
863 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
864 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
865 */
866 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
867 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
868
869
870 Double_t invmass_mc = 0;
871 Double_t invmass_reco = 0;
872 CbmAnaConversionKinematicParams params_mc =
873 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(
874 fVector_mctrack[i], fVector_mctrack[j], fVector_mctrack[k], fVector_mctrack[l]);
875 CbmAnaConversionKinematicParams params_reco =
876 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
877 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
878 invmass_mc = params_mc.fMinv;
879 invmass_reco = params_reco.fMinv;
880
881 fhTest2_2rich_invmass_gee_mc->Fill(invmass_mc);
882 fhTest2_2rich_invmass_gee_refitted->Fill(invmass_reco);
883 fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
884 fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
885
886 fhTest2_2rich_pt_vs_rap_gee->Fill(params_reco.fPt, params_reco.fRapidity);
887 fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
888 }
889 }
890 if (motherId1 == motherId3 && motherId2 == motherId4) {
891 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2) continue;
892 if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
893 || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
894 /*
895 fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
896 fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
897 fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
898 fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
899 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
900 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
901 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
902 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
903 */
904 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
905 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
906
907
908 Double_t invmass_mc = 0;
909 Double_t invmass_reco = 0;
910 CbmAnaConversionKinematicParams params_mc =
911 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(
912 fVector_mctrack[i], fVector_mctrack[j], fVector_mctrack[k], fVector_mctrack[l]);
913 CbmAnaConversionKinematicParams params_reco =
914 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
915 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
916 invmass_mc = params_mc.fMinv;
917 invmass_reco = params_reco.fMinv;
918
919 fhTest2_2rich_invmass_gee_mc->Fill(invmass_mc);
920 fhTest2_2rich_invmass_gee_refitted->Fill(invmass_reco);
921 fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
922 fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
923
924 fhTest2_2rich_pt_vs_rap_gee->Fill(params_reco.fPt, params_reco.fRapidity);
925 fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
926 }
927 }
928 if (motherId1 == motherId4 && motherId2 == motherId3) {
929 if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2) continue;
930 if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
931 || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
932 /*
933 fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
934 fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
935 fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
936 fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
937 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
938 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
939 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
940 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
941 */
942 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
943 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
944
945
946 Double_t invmass_mc = 0;
947 Double_t invmass_reco = 0;
948 CbmAnaConversionKinematicParams params_mc =
949 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(
950 fVector_mctrack[i], fVector_mctrack[j], fVector_mctrack[k], fVector_mctrack[l]);
951 CbmAnaConversionKinematicParams params_reco =
952 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
953 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
954 invmass_mc = params_mc.fMinv;
955 invmass_reco = params_reco.fMinv;
956
957 fhTest2_2rich_invmass_gee_mc->Fill(invmass_mc);
958 fhTest2_2rich_invmass_gee_refitted->Fill(invmass_reco);
959 fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
960 fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
961
962 fhTest2_2rich_pt_vs_rap_gee->Fill(params_reco.fPt, params_reco.fRapidity);
963 fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
964 }
965 }
966
967
968 // ===================================================================================================
969 // HERE DECAY pi0 -> gamma gamma -> e+e- e+e-
970 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
971 && grandmotherId1 == grandmotherId4) {
972 if (mcGrandmotherPdg1 != 111) continue; // 111 = pi0, 221 = eta
973
974 TVector3 pi0start;
975 grandmother1->GetStartVertex(pi0start);
976 /*
977 fhTest2_startvertexElectrons_gg->Fill(pi0start_i.Z());
978 fhTest2_startvertexElectrons_gg->Fill(pi0start_j.Z());
979 fhTest2_startvertexElectrons_gg->Fill(pi0start_k.Z());
980 fhTest2_startvertexElectrons_gg->Fill(pi0start_l.Z());
981 fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
982 fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
983 fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
984 fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
985 */
986
987 // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
988 if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1 || pi0start_l.Z() > 1) continue;
989
990 Double_t invmass_mc = 0;
991 Double_t invmass_reco = 0;
992 CbmAnaConversionKinematicParams params_mc =
993 CbmAnaConversionKinematicParams::KinematicParams_4particles_MC(fVector_mctrack[i], fVector_mctrack[j],
994 fVector_mctrack[k], fVector_mctrack[l]);
995 CbmAnaConversionKinematicParams params_reco =
996 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
997 fVector_momenta[i], fVector_momenta[j], fVector_momenta[k], fVector_momenta[l]);
998 invmass_mc = params_mc.fMinv;
999 invmass_reco = params_reco.fMinv;
1000
1001 fhTest2_2rich_invmass_gg_mc->Fill(invmass_mc);
1002 fhTest2_2rich_invmass_gg_refitted->Fill(invmass_reco);
1003 fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
1004 fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
1005
1006 fhTest2_2rich_pt_vs_rap_gg->Fill(params_reco.fPt, params_reco.fRapidity);
1007 fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt, params_reco.fRapidity);
1008
1009 cout << "########################################################"
1010 "##############"
1011 << endl;
1012
1013 fill++;
1014 }
1015 }
1016 }
1017 }
1018 }
1019 }
1020}
1021
1022
1023Int_t CbmAnaConversionTest2::NofDaughters(Int_t motherId)
1024{
1025 Int_t nofDaughters = 0;
1026 for (unsigned int i = 0; i < fVector_mctrack.size(); i++) {
1027 Int_t motherId_temp = fVector_mctrack[i]->GetMotherId();
1028 if (motherId == motherId_temp) nofDaughters++;
1029 }
1030
1031
1032 return nofDaughters;
1033}
CbmStsTrack.h
Data class for STS tracks.
CbmAnaConversionKinematicParams::KinematicParams_4particles_MC
static CbmAnaConversionKinematicParams KinematicParams_4particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
Definition CbmAnaConversionKinematicParams.h:34
CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco
static CbmAnaConversionKinematicParams KinematicParams_4particles_Reco(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
Definition CbmAnaConversionKinematicParams.h:87
CbmAnaConversionTest2::fVector_gt
std::vector< CbmGlobalTrack * > fVector_gt
Definition CbmAnaConversionTest2.h:74
CbmAnaConversionTest2::fVector_momenta
std::vector< TVector3 > fVector_momenta
Definition CbmAnaConversionTest2.h:75
CbmAnaConversionTest2::fVector_richIndex
std::vector< int > fVector_richIndex
Definition CbmAnaConversionTest2.h:78
CbmAnaConversionTest2::fHistoList_test2
std::vector< TH1 * > fHistoList_test2
Definition CbmAnaConversionTest2.h:62
CbmAnaConversionTest2::fVector_gtIndex
std::vector< int > fVector_gtIndex
Definition CbmAnaConversionTest2.h:77
CbmAnaConversionTest2::fVector_mctrack
std::vector< CbmMCTrack * > fVector_mctrack
Definition CbmAnaConversionTest2.h:76
CbmAnaConversionTest2::NofDaughters
Int_t NofDaughters(Int_t motherId)
Definition CbmAnaConversionTest2.cxx:1023
CbmGlobalTrack::GetStsTrackIndex
int32_t GetStsTrackIndex() const
Definition CbmGlobalTrack.h:44
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition CbmL1PFFitter.cxx:501
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, const std::vector< CbmMvdHit > &vMvdHits, const std::vector< CbmStsHit > &vStsHits, const std::vector< int > &pidHypo)
Definition CbmL1PFFitter.cxx:172
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21