CbmRoot
Loading...
Searching...
No Matches
CbmAnaConversion.cxx
Go to the documentation of this file.
1/* Copyright (C) 2014-2019 GSI/JINR-LIT, Darmstadt/Dubna
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Sascha Reinecke, Florian Uhlig, Andrey Lebedev [committer] */
4
14#include "CbmAnaConversion.h"
15
16// includes of CBMROOT classes
17#include "CbmDrawHist.h"
18#include "CbmGlobalTrack.h"
19#include "CbmKFParticleFinder.h"
21#include "CbmL1PFFitter.h"
22#include "CbmMCTrack.h"
23#include "CbmRichHit.h"
24#include "CbmRichPoint.h"
25#include "CbmRichRing.h"
26#include "CbmStsTrack.h"
27#include "CbmTrackMatchNew.h"
28#include "CbmUtils.h"
29
30#include "FairMCPoint.h"
31#include "FairTrackParam.h"
32#include <Logger.h>
33
34#include "KFParticleTopoReconstructor.h"
35
36
37// includes of standard c++ classes or ROOT classes
38#include "TCanvas.h"
39#include "TClonesArray.h"
40#include "TH1.h"
41#include "TH1D.h"
42#include "TH2D.h"
43#include "TH3.h"
44#include "TRandom3.h"
45
46#include <boost/assign/list_of.hpp>
47
48#include <iomanip>
49#include <iostream>
50#include <string>
51
52
53// includes of further conversion classes
55#include "CbmAnaConversionKF.h"
64
65
66#define M2E 2.6112004954086e-7
67
68using namespace std;
69using boost::assign::list_of;
70
72 : FairTask("CbmAnaConversion")
73 , DoTomography(0)
74 , DoRichAnalysis(0)
75 , DoKFAnalysis(0)
76 , DoReconstruction(0)
77 , DoPhotons(0)
78 , DoPhotons2(0)
79 , DoRecoFull(0)
80 , DoTest(0)
81 , fhPdgCodes(NULL)
82 , fhNofElPrim(NULL)
83 , fhNofElSec(NULL)
84 , fhNofElAll(NULL)
85 , fhElectronSources(NULL)
86 , fhNofPi0_perEvent(NULL)
87 , fhNofPi0_perEvent_cut(NULL)
88 , fhNofPi0_perEvent_cut2(NULL)
89 , fhNofPi0_perEvent_cut3(NULL)
90 , fhNofEta_perEvent(NULL)
91 , fhNofEta_perEvent_cut(NULL)
92 , fhNofEta_perEvent_cut2(NULL)
93 , fhPi0_z(NULL)
94 , fhPi0_z_cut(NULL)
95 , fhPi0_pt(NULL)
96 , fhPi0_pt_vs_rap(NULL)
97 , fhPi0_theta(NULL)
98 , fhPi0_theta_vs_rap(NULL)
99 , fhEta_pt(NULL)
100 , fhEta_pt_vs_rap(NULL)
101 , fhEta_theta(NULL)
102 , fhEta_theta_vs_rap(NULL)
103 , fhRho_pt(NULL)
104 , fhRho_pt_vs_rap(NULL)
105 , fhRho_theta(NULL)
106 , fhRho_theta_vs_rap(NULL)
107 , fhOmega_pt(NULL)
108 , fhOmega_pt_vs_rap(NULL)
109 , fhOmega_theta(NULL)
110 , fhOmega_theta_vs_rap(NULL)
111 , fhElectronsFromPi0_z(NULL)
112 , fhNofTracks_mctrack(NULL)
113 , fhNofTracks_globaltrack(NULL)
114 , fhInvariantMass_test(NULL)
115 , fhInvariantMass_test2(NULL)
116 , fhInvariantMass_test3(NULL)
117 , fhInvariantMassReco_test(NULL)
118 , fhInvariantMassReco_test2(NULL)
119 , fhInvariantMassReco_test3(NULL)
120 , fhInvariantMassReco_pi0(NULL)
121 , fhMomentum_MCvsReco(NULL)
122 , fhMomentum_MCvsReco_diff(NULL)
123 , fhSearchGammas(NULL)
124 , fPrimVertex(NULL)
125 , fKFVertex()
126 , fRichPoints(NULL)
127 , fRichRings(NULL)
128 , fRichRingMatches(NULL)
129 , fMcTracks(NULL)
130 , fStsTracks(NULL)
131 , fStsTrackMatches(NULL)
132 , fGlobalTracks(NULL)
133 , fhANN_output_electrons(NULL)
134 , fhANN_output_electrons2(NULL)
135 , fhANN_output_electrons_chiCut(NULL)
136 , fhANN_output_electrons_vs_p(NULL)
137 , fhANN_output_else(NULL)
138 , fhANN_output_else2(NULL)
139 , fhANN_output_else_chiCut(NULL)
140 , fhANN_output_else_vs_p(NULL)
141 , fEventNum(0)
142 , test(0)
143 , testint(0)
144 , fAnalyseMode(0)
145 , fKFparticle(NULL)
146 , fKFparticleFinderQA(NULL)
147 , fKFtopo(NULL)
148 , trackindexarray()
149 , particlecounter(0)
150 , particlecounter_2daughters(0)
151 , particlecounter_3daughters(0)
152 , particlecounter_4daughters(0)
153 , particlecounter_all(0)
154 , fNofGeneratedPi0_allEvents(0)
155 , fNofPi0_kfparticle_allEvents(0)
156 , fNofGeneratedPi0(0)
157 , fNofPi0_kfparticle(0)
158 , fhPi0Ratio(NULL)
159 , fhPi0_mass(NULL)
160 , fhPi0_NDaughters(NULL)
161 , fHistoList()
162 , fHistoList_MC()
163 , fHistoList_tomography()
164 , fHistoList_reco()
165 , fHistoList_reco_mom()
166 , fHistoList_kfparticle()
167 , fHistoList_richrings()
168 , fHistoList_furtherAnalyses()
169 , fMCTracklist()
170 , fMCTracklist_all()
171 , fRecoTracklist()
172 , fRecoTracklistEPEM()
173 , fRecoTracklistEPEM_id()
174 , fRecoTracklistEPEM_chi()
175 , fRecoTracklistEPEM_gtid()
176 , fTestTracklist()
177 , fTestTracklist_momentum()
178 , fTestTracklist_chi()
179 , fTestTracklist_richInd()
180 , fTestTracklist_ndf()
181 , fTestTracklist_nofhits()
182 , fTestTracklist_noRichInd()
183 , fTestTracklist_noRichInd_MCindex()
184 , fTestTracklist_noRichInd_momentum()
185 , fTestTracklist_noRichInd_chi()
186 , fTestTracklist_noRichInd_richInd()
187 , fTestTracklist_noRichInd_gTrackId()
188 , fTestTracklist_noRichInd_ndf()
189 , fTestTracklist_noRichInd_nofhits()
190 , fRecoMomentum()
191 , fRecoRefittedMomentum()
192 , fhNofElectrons_4epem(NULL)
193 , fhPi0_MC_occurence(NULL)
194 , fhPi0_MC_occurence2(NULL)
195 , fhPi0_Reco_occurence(NULL)
196 , fhPi0_Reco_occurence2(NULL)
197 , fhPi0_Reco_angles(NULL)
198 , fhPi0_Reco_chi(NULL)
199 , fhPi0_Reco_ndf(NULL)
200 , fhPi0_Reco_ndf_vs_chi(NULL)
201 , fhPi0_Reco_ndf_vs_startvertex(NULL)
202 , fhPi0_Reco_startvertex_vs_chi(NULL)
203 , fhPi0_Reco_startvertex_vs_nofhits(NULL)
204 , fhPi0_noRichInd_diffPhi(NULL)
205 , fhPi0_noRichInd_diffTheta(NULL)
206 , fhPi0_Reco_invmass_cases(NULL)
207 , fhPi0_Reco_noRichInd_invmass_cases(NULL)
208 , fhPi0_Reco_invmass(NULL)
209 , fhPi0_Reco_invmass_mc(NULL)
210 , fhPi0_Reco_noRichInd_invmass(NULL)
211 , fhPi0_Reco_noRichInd_invmass_mc(NULL)
212 , fhPi0_Reco_noRichInd_ndf_vs_nofhits(NULL)
213 , fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi(NULL)
214 , fhPi0_Reco_ndf_vs_nofhits(NULL)
215 , fhPi0_Reco_ndf_vs_nofhits_withChi(NULL)
216 , fhPi0_Reco_noRichInd_chi_vs_momentum(NULL)
217 , fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0(NULL)
218 , fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target(NULL)
219 , fhPi0_Reco_noRichInd_chi_vs_momentum_eRest(NULL)
220 , fhPi0_Reco_noRichInd_chi_vs_pt(NULL)
221 , fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0(NULL)
222 , fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target(NULL)
223 , fhPi0_Reco_noRichInd_chi_vs_pt_eRest(NULL)
224 , fhPi0_Reco_chi_vs_momentum(NULL)
225 , fhPi0_Reco_chi_vs_momentum_eFromPi0(NULL)
226 , fhPi0_Reco_chi_vs_momentum_eFromPi0_Target(NULL)
227 , fhPi0_Reco_chi_vs_pt(NULL)
228 , fhPi0_Reco_chi_vs_pt_eFromPi0(NULL)
229 , fhPi0_Reco_chi_vs_pt_eFromPi0_Target(NULL)
230 , timer_all()
231 , fTime_all(0.)
232 , timer_exec()
233 , fTime_exec(0.)
234 , timer_mc()
235 , fTime_mc(0.)
236 , timer_rec()
237 , fTime_rec(0.)
238 , fAnaTomography(NULL)
239 , fAnaRich(NULL)
240 , fAnaKF(NULL)
241 , fAnaReco(NULL)
242 , fAnaPhotons(NULL)
243 , fAnaPhotons2(NULL)
244 , fAnaRecoFull(NULL)
245 , fAnaTest(NULL)
246 , fAnaTest2(NULL)
247{
248}
249
251
253{
254 //timer_all.Reset();
255 timer_all.Start();
256
257 cout << "CbmAnaConversion::Init" << endl;
258 FairRootManager* ioman = FairRootManager::Instance();
259 if (NULL == ioman) { Fatal("CbmAnaConversion::Init", "RootManager not instantised!"); }
260
261 fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
262 if (NULL == fRichPoints) { Fatal("CbmAnaConversion::Init", "No RichPoint array!"); }
263
264 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
265 if (NULL == fMcTracks) { Fatal("CbmAnaConversion::Init", "No MCTrack array!"); }
266
267 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
268 if (NULL == fStsTracks) { Fatal("CbmAnaConversion::Init", "No StsTrack array!"); }
269
270 fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
271 if (NULL == fStsTrackMatches) { Fatal("CbmAnaConversion::Init", "No StsTrackMatch array!"); }
272
273 fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
274 if (NULL == fGlobalTracks) { Fatal("CbmAnaConversion::Init", "No GlobalTrack array!"); }
275
276 // Get pointer to PrimaryVertex object from IOManager if it exists
277 // The old name for the object is "PrimaryVertex" the new one
278 // "PrimaryVertex." Check first for the new name
279 fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex."));
280 if (nullptr == fPrimVertex) { fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex")); }
281 if (nullptr == fPrimVertex) { LOG(fatal) << "No PrimaryVertex array!"; }
282
283 fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
284 if (NULL == fRichRings) { Fatal("CbmAnaConversion::Init", "No RichRing array!"); }
285
286 fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
287 if (NULL == fRichRingMatches) { Fatal("CbmAnaConversion::Init", "No RichRingMatch array!"); }
288
290
291 particlecounter = 0;
292
293 testint = 0;
294 test = 0;
295
298
299
300 DoTomography = 1;
301 DoRichAnalysis = 1;
302 DoKFAnalysis = 0;
304 DoPhotons = 1;
305 DoPhotons2 = 1;
306 DoRecoFull = 1;
307 DoTest = 1;
308
309 if (DoTomography) {
312 }
313 if (DoRichAnalysis) {
315 fAnaRich->Init();
316 }
317 if (DoKFAnalysis) {
320 fAnaKF->Init();
321 }
322 if (DoReconstruction) {
324 fAnaReco->Init();
325 }
326 if (DoPhotons) {
328 fAnaPhotons->Init();
329 }
330 if (DoPhotons2) {
333 }
334 if (DoRecoFull) {
337 }
338 if (DoTest) {
340 fAnaTest->Init();
342 fAnaTest2->Init();
343 }
344
345 timer_all.Stop();
346 fTime_all += timer_all.RealTime();
347
348 return kSUCCESS;
349}
350
352{
353 fHistoList.clear();
354 fHistoList_kfparticle.clear();
356
357
358 fhPdgCodes = new TH1D("fhPdfCodes", "fhPdgCodes;pdg code;#", 1000, 0, 1000);
359 fhNofElPrim = new TH1D("fhNofElPrim", "fhNofElPrim;Nof prim El;Entries", 10., -0.5, 9.5);
360 fhNofElSec = new TH1D("fhNofElSec", "fhNofElSec;Nof Sec El;Entries", 20., -0.5, 19.5);
361 fhNofElAll = new TH1D("fhNofElAll", "fhNofElAll;Nof All El;Entries", 30., -0.5, 29.5);
362 fhNofPi0_perEvent = new TH1D("fhNofPi0_perEvent", "fhNofPi0_perEvent;Nof pi0;Entries", 1000., -0.5, 999.5);
364 new TH1D("fhNofPi0_perEvent_cut", "fhNofPi0_perEvent_cut (Z<10cm);Nof pi0;Entries", 800., -0.5, 799.5);
366 new TH1D("fhNofPi0_perEvent_cut2", "fhNofPi0_perEvent_cut2 (motherId = -1);Nof pi0;Entries", 800., -0.5, 799.5);
367 fhNofPi0_perEvent_cut3 = new TH1D(
368 "fhNofPi0_perEvent_cut3", "fhNofPi0_perEvent_cut3 (conversion before 70cm);Nof pi0;Entries", 100., -0.5, 99.5);
369 fhNofEta_perEvent = new TH1D("fhNofEta_perEvent", "fhNofEta_perEvent;Nof eta;Entries", 100., -0.5, 99.5);
371 new TH1D("fhNofEta_perEvent_cut", "fhNofEta_perEvent_cut (Z<4cm);Nof eta;Entries", 100., -0.5, 99.5);
373 new TH1D("fhNofEta_perEvent_cut2", "fhNofEta_perEvent_cut2 (motherId = -1);Nof eta;Entries", 100., -0.5, 99.5);
374 fhPi0_z = new TH1D("fhPi0_z", "fhPi0_z;z in cm;Entries", 600., -0.5, 599.5);
375 fhPi0_z_cut = new TH1D("fhPi0_z_cut", "fhPi0_z_cut;z in cm;Entries", 600., -0.5, 599.5);
376 fhPi0_pt = new TH1D("fhPi0_pt", "fhPi0_pt;p_{t} in GeV/c;Entries", 200., 0., 10.);
378 new TH2D("fhPi0_pt_vs_rap", "fhPi0_pt_vs_rap;p_{t} in GeV/c; rapidity y", 240, -2., 10., 270, -2., 7.);
379 fhPi0_theta = new TH1D("fhPi0_theta", "fhPi0_theta;emission angle #theta in deg;Entries", 180., 0., 180.);
381 new TH2D("fhPi0_theta_vs_rap", "fhPi0_theta_vs_rap;theta in deg;rapidity y", 180., 0., 180., 270, -2., 7.);
382 fhEta_pt = new TH1D("fhEta_pt", "fhEta_pt;p_{t} in GeV;Entries", 200., 0., 10.);
384 new TH2D("fhEta_pt_vs_rap", "fhEta_pt_vs_rap;p_{t} in GeV; rapidity y", 240, -2., 10., 270, -2., 7.);
385 fhEta_theta = new TH1D("fhEta_theta", "fhEta_theta;emission angle #theta in deg;Entries", 180., 0., 180.);
386 fhEta_theta_vs_rap = new TH2D("fhEta_theta_vs_rap", "fhEta_theta_vs_rap;emission angle #theta in deg;rapidity y",
387 180., 0., 180., 270, -2., 7.);
388 fhRho_pt = new TH1D("fhRho_pt", "fhRho_pt;p_{t} in GeV/c;Entries", 200., 0., 10.);
390 new TH2D("fhRho_pt_vs_rap", "fhRho_pt_vs_rap;p_{t} in GeV/c; rapidity y", 240, -2., 10., 270, -2., 7.);
391 fhRho_theta = new TH1D("fhRho_theta", "fhRho_theta;emission angle #theta in deg;Entries", 180., 0., 180.);
392 fhRho_theta_vs_rap = new TH2D("fhRho_theta_vs_rap", "fhRho_theta_vs_rap;emission angle #theta in deg;rapidity y",
393 180., 0., 180., 270, -2., 7.);
394 fhOmega_pt = new TH1D("fhOmega_pt", "fhOmega_pt;p_{t} in GeV/c;Entries", 200., 0., 10.);
396 new TH2D("fhOmega_pt_vs_rap", "fhOmega_pt_vs_rap;p_{t} in GeV; rapidity y", 240, -2., 10., 270, -2., 7.);
397 fhOmega_theta = new TH1D("fhOmega_theta", "fhOmega_theta;emission angle #theta in deg;Entries", 180., 0., 180.);
399 new TH2D("fhOmega_theta_vs_rap", "fhOmega_theta_vs_rap;emission angle #theta in deg;rapidity y", 180., 0., 180.,
400 270, -2., 7.);
401
402
403 fhElectronSources = new TH1D("fhElectronSources", "fhElectronSources;Source;Entries", 6., 0., 6.);
404 fhElectronsFromPi0_z = new TH1D(
405 "fhElectronsFromPi0_z", "fhElectronsFromPi0_z (= pos. of gamma conversion);z [cm];Entries", 600., -0.5, 599.5);
406 fHistoList.push_back(fhPdgCodes);
407 fHistoList.push_back(fhNofPi0_perEvent);
411 fHistoList.push_back(fhNofEta_perEvent);
414 fHistoList.push_back(fhPi0_z);
415 fHistoList.push_back(fhPi0_z_cut);
416 fHistoList.push_back(fhPi0_pt);
417 fHistoList.push_back(fhPi0_pt_vs_rap);
418 fHistoList.push_back(fhPi0_theta);
420 fHistoList.push_back(fhEta_pt);
421 fHistoList.push_back(fhEta_pt_vs_rap);
422 fHistoList.push_back(fhEta_theta);
424 fHistoList.push_back(fhRho_pt);
425 fHistoList.push_back(fhRho_pt_vs_rap);
426 fHistoList.push_back(fhRho_theta);
428 fHistoList.push_back(fhOmega_pt);
429 fHistoList.push_back(fhOmega_pt_vs_rap);
430 fHistoList.push_back(fhOmega_theta);
432 fHistoList.push_back(fhElectronSources);
434
435 fhElectronSources->GetXaxis()->SetBinLabel(1, "gamma");
436 fhElectronSources->GetXaxis()->SetBinLabel(2, "pi0");
437 fhElectronSources->GetXaxis()->SetBinLabel(3, "eta");
438 fhElectronSources->GetXaxis()->SetBinLabel(4, "else");
439 fhElectronSources->GetXaxis()->SetBinLabel(5, "gamma from pi0");
440 fhElectronSources->GetXaxis()->SetBinLabel(6, "gamma from eta");
441
442
443 fhNofTracks_mctrack = new TH1D("fhNofTracks_mctrack", "fhNofTracks_mctrack;nof;#", 1000., 0., 1000.);
445 fhNofTracks_globaltrack = new TH1D("fhNofTracks_globaltrack", "fhNofTracks_globaltrack;nof;#", 1000., 0., 1000.);
447
448
449 // for UrQMD events (invariant mass from pi0 -> gamma + gamma
450 fhInvariantMass_test = new TH1D("fhInvariant", "fhInvariant;mass in GeV/c^{2}];#", 2000, 0., 2.);
451 fhInvariantMass_test2 = new TH1D("fhInvariant2", "fhInvariant2;mass in GeV/c^{2}];#", 2000, 0., 2.);
452 fhInvariantMass_test3 = new TH1D("fhInvariant3", "fhInvariant3;mass in GeV/c^{2}];#", 2000, 0., 2.);
456
457 fhInvariantMassReco_test = new TH1D("fhInvariantReco", "fhInvariantReco;mass [GeV/c^2];#", 2000, 0., 2.);
458 fhInvariantMassReco_test2 = new TH1D("fhInvariantReco2", "fhInvariantReco2;mass [GeV/c^2];#", 2000, 0., 2.);
459 fhInvariantMassReco_test3 = new TH1D("fhInvariantReco3", "fhInvariantReco3;mass [GeV/c^2];#", 2000, 0., 2.);
463
464 fhInvariantMassReco_pi0 = new TH1D("fhInvariantReco_pi0", "fhInvariantReco_pi0;mass [GeV/c^2];#", 2000, 0., 2.);
466
467
468 // for reconstructed tracks
469 fhMomentum_MCvsReco = new TH2D("fhMomentum_MCvsReco", "fhMomentum_MCvsReco;MC;Reco", 400, 0., 40., 400, 0., 40.);
471 new TH1D("fhMomentum_MCvsReco_diff", "fhMomentum_MCvsReco_diff;(MC-Reco)/MC", 500, -0.01, 4.);
474
475
476 fhSearchGammas = new TH1D("fhSearchGammas", "fhSearchGammas;mass;#", 100, -0.005, 0.995);
477 fHistoList.push_back(fhSearchGammas);
478
479
480 fhANN_output_electrons = new TH1D("fhANN_output_electrons", "fhANN_output_electrons;ann output", 400, -2, 2.);
481 fhANN_output_electrons2 = new TH1D("fhANN_output_electrons2", "fhANN_output_electrons2;ann output", 400, -2, 2.);
483 new TH1D("fhANN_output_electrons_chiCut", "fhANN_output_electrons_chiCut;ann output", 400, -2, 2.);
485 new TH2D("fhANN_output_electrons_vs_p", "fhANN_output_electrons_vs_p;momentum in GeV/c; ann output", 100, 0., 10.,
486 400, -2, 2.);
487 fhANN_output_else = new TH1D("fhANN_output_else", "fhANN_output_else;ann output", 400, -2, 2.);
488 fhANN_output_else2 = new TH1D("fhANN_output_else2", "fhANN_output_else2;ann output", 400, -2, 2.);
489 fhANN_output_else_chiCut = new TH1D("fhANN_output_else_chiCut", "fhANN_output_else_chiCut;ann output", 400, -2, 2.);
490 fhANN_output_else_vs_p = new TH2D("fhANN_output_else_vs_p", "fhANN_output_else_vs_p;momentum in GeV/c; ann output",
491 100, 0., 10., 400, -2, 2.);
496 fHistoList.push_back(fhANN_output_else);
500
501
502 // #############################################
503 // Histograms related to KFParticle results
504 fhPi0_NDaughters = new TH1D("fhPi0_NDaughters", "fhPi0_NDaughters;number of daughers;#", 4, 0.5, 4.5);
505 fhPi0Ratio = new TH1D("fhPi0Ratio", "fhPi0Ratio; ratio;#", 1000, 0., 0.02);
506 fhPi0_mass = new TH1D("fhPi0_mass", "fhPi0_mass;mass;#", 500, 0., 0.5);
510
511
513 new TH1D("fhNofElectrons_4epem", "fhNofElectrons_4epem;number of electrons per event;#", 101, -0.5, 100.5);
515
516 fhPi0_MC_occurence = new TH1D("fhPi0_MC_occurence", "fhPi0_MC_occurence;;#", 20, 0, 20);
518 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(1, "== -1: all pi0 from target");
519 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(2, "all pi0 -> gg");
520 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(3, "all g -> e+e-");
521 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(4, "both conv before 70cm");
522 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(5, "both conv in target");
523 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(6, "...");
524 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(7, "!= -1: all pi0 not from target");
525 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(8, "all pi0 -> gg");
526 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(9, "all g -> e+e-");
527 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(10, "both conv before 70cm");
528 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(11, "both conv in target");
529 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(12, "...");
530 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(13, "daughters == 0");
531 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(14, "daughters == 1");
532 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(15, "daughters == 2");
533 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(16, "daughters == 3");
534 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(17, "daughters == 4");
535 fhPi0_MC_occurence->GetXaxis()->SetBinLabel(18, "daughters > 4");
536
537 fhPi0_MC_occurence2 = new TH1D("fhPi0_MC_occurence2", "fhPi0_MC_occurence2;;#", 20, 0, 20);
539 fhPi0_MC_occurence2->GetXaxis()->SetBinLabel(1, "!= -1: all pi0 not from target");
540 fhPi0_MC_occurence2->GetXaxis()->SetBinLabel(2, "all pi0 -> gg");
541 fhPi0_MC_occurence2->GetXaxis()->SetBinLabel(3, "all g -> e+e-");
542 fhPi0_MC_occurence2->GetXaxis()->SetBinLabel(4, "both conv before 70cm");
543 fhPi0_MC_occurence2->GetXaxis()->SetBinLabel(5, "both conv in target");
544 fhPi0_MC_occurence2->GetXaxis()->SetBinLabel(6, "...");
545
546 fhPi0_Reco_occurence = new TH1D("fhPi0_Reco_occurence", "fhPi0_Reco_occurence;;#", 16, 0, 16);
548 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(1, "4 e from pi0 (not same)");
549 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(2, "test");
550 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(3, "x 1 e from same pi0");
551 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(4, "x 2 e from same pi0");
552 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(5, "x 3 e from same pi0");
553 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(6, "x 4 e from same pi0");
554 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(7, "x >4 e from same pi0");
555 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(8, "...");
556 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(9, "==4, within cuts");
557 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(10, "==4, chi-check");
558 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(11, "==4, chi+cuts");
559 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(12, "richInd: ==1");
560 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(13, "richInd: ==2");
561 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(14, "richInd: ==3");
562 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(15, "richInd: ==4");
563 fhPi0_Reco_occurence->GetXaxis()->SetBinLabel(16, "richInd: >4");
564
565 fhPi0_Reco_occurence2 = new TH1D("fhPi0_Reco_occurence2", "fhPi0_Reco_occurence2;;#", 16, 0, 16);
567 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(1, "4 e from pi0 (not same)");
568 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(2, "test");
569 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(3, "x 1 e from same pi0");
570 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(4, "x 2 e from same pi0");
571 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(5, "x 3 e from same pi0");
572 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(6, "x 4 e from same pi0");
573 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(7, "x >4 e from same pi0");
574 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(8, "...");
575 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(9, "==4, within cuts");
576 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(10, "==4, chi-check");
577 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(11, "==4, chi+cuts");
578 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(12, "...");
579 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(13, "...");
580 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(14, "...");
581 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(15, "...");
582 fhPi0_Reco_occurence2->GetXaxis()->SetBinLabel(16, "...");
583
584
585 fhPi0_Reco_angles = new TH1D("fhPi0_Reco_angles", "fhPi0_Reco_angles;angle;#", 500, 0, 50);
587 fhPi0_Reco_chi = new TH1D("fhPi0_Reco_chi", "fhPi0_Reco_chi;chi;#", 500, 0, 500);
589 fhPi0_Reco_ndf = new TH1D("fhPi0_Reco_ndf", "fhPi0_Reco_ndf;ndf;#", 500, 0, 50);
592 new TH2D("fhPi0_Reco_ndf_vs_chi", "fhPi0_Reco_ndf_vs_chi;ndf;chi", 51, -0.5, 50.5, 500, 0, 50);
595 "fhPi0_Reco_ndf_vs_startvertex", "fhPi0_Reco_ndf_vs_startvertex;ndf;startvertex", 51, -0.5, 50.5, 101, -0.5, 100.5);
598 "fhPi0_Reco_startvertex_vs_chi", "fhPi0_Reco_startvertex_vs_chi;startvertex;chi", 101, -0.5, 100.5, 100, 0, 100);
601 new TH2D("fhPi0_Reco_startvertex_vs_nofhits", "fhPi0_Reco_startvertex_vs_nofhits;startvertex;nofhits", 101, -0.5,
602 100.5, 21, -0.5, 20.5);
605 new TH1D("fhPi0_noRichInd_diffPhi", "fhPi0_noRichInd_diffPhi;phi difference;#", 150, 0, 150);
608 new TH1D("fhPi0_noRichInd_diffTheta", "fhPi0_noRichInd_diffTheta;theta difference;#", 150, 0, 150);
610
612 new TH2D("fhPi0_Reco_invmass_cases", "fhPi0_Reco_invmass_cases;cases;invmass [GeV]", 6, 0, 6, 300, 0, 3);
614 fhPi0_Reco_invmass_cases->GetXaxis()->SetBinLabel(1, "..");
615 fhPi0_Reco_invmass_cases->GetXaxis()->SetBinLabel(2, "no cuts");
616 fhPi0_Reco_invmass_cases->GetXaxis()->SetBinLabel(3, "only chi");
617 fhPi0_Reco_invmass_cases->GetXaxis()->SetBinLabel(4, "only cuts");
618 fhPi0_Reco_invmass_cases->GetXaxis()->SetBinLabel(5, "both");
619 fhPi0_Reco_invmass_cases->GetXaxis()->SetBinLabel(6, "mc-true");
621 new TH2D("fhPi0_Reco_noRichInd_invmass_cases", "fhPi0_Reco_noRichInd_invmass_cases;cases;invmass in GeV^{2}", 11, 0,
622 11, 300, 0, 3);
624 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(1, "..");
625 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(2, "no cuts");
626 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(3, "only chi");
627 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(4, "only cuts");
628 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(5, "both");
629 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(6, "mc-true");
630 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(7, "richInd-no cuts");
631 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(8, "richInd-chi");
632 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(9, "richInd-cuts");
633 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(10, "richInd-both");
634 fhPi0_Reco_noRichInd_invmass_cases->GetXaxis()->SetBinLabel(11, "richInd-mctrue");
635
636 fhPi0_Reco_invmass = new TH1D("fhPi0_Reco_invmass", "fhPi0_Reco_invmass;invmass [GeV];#", 300, 0, 3);
638 fhPi0_Reco_invmass_mc = new TH1D("fhPi0_Reco_invmass_mc", "fhPi0_Reco_invmass_mc;invmass_mc [GeV];#", 300, 0, 3);
641 new TH1D("fhPi0_Reco_noRichInd_invmass", "fhPi0_Reco_noRichInd_invmass;invmass_noRichInd [GeV];#", 300, 0, 3);
643 fhPi0_Reco_noRichInd_invmass_mc = new TH1D("fhPi0_Reco_noRichInd_invmass_mc",
644 "fhPi0_Reco_noRichInd_invmass_mc;invmass_noRichInd_mc [GeV];#", 300, 0, 3);
646
648 new TH2D("fhPi0_Reco_noRichInd_ndf_vs_nofhits", "fhPi0_Reco_noRichInd_ndf_vs_nofhits;ndf;nofhits", 31, -0.5, 30.5,
649 21, -0.5, 20.5);
652 new TH2D("fhPi0_Reco_ndf_vs_nofhits", "fhPi0_Reco_ndf_vs_nofhits;ndf;nofhits", 31, -0.5, 30.5, 21, -0.5, 20.5);
655 new TH2D("fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi", "fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi;ndf;nofhits",
656 31, -0.5, 30.5, 21, -0.5, 20.5);
659 new TH2D("fhPi0_Reco_ndf_vs_nofhits_withChi", "fhPi0_Reco_ndf_vs_nofhits_withChi;ndf;nofhits", 31, -0.5, 30.5, 21,
660 -0.5, 20.5);
662
663
664 // ################################################
665
666 fhPi0_Reco_noRichInd_chi_vs_momentum = new TH2D("fhPi0_Reco_noRichInd_chi_vs_momentum",
667 "fhPi0_Reco_noRichInd_chi_vs_momentum;momentum [GeV];#chi^{2} of "
668 "momentum fit",
669 200, 0., 10., 1000, 0., 100.);
671 fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0 = new TH2D("fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0",
672 "fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0;momentum "
673 "[GeV];#chi^{2} of momentum fit",
674 200, 0., 10., 1000, 0., 100.);
677 new TH2D("fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target",
678 "fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target;momentum "
679 "[GeV];#chi^{2} of momentum fit",
680 200, 0., 10., 1000, 0., 100.);
682 fhPi0_Reco_noRichInd_chi_vs_momentum_eRest = new TH2D("fhPi0_Reco_noRichInd_chi_vs_momentum_eRest",
683 "fhPi0_Reco_noRichInd_chi_vs_momentum_eRest;momentum "
684 "[GeV];#chi^{2} of momentum fit",
685 200, 0., 10., 1000, 0., 100.);
687
689 new TH2D("fhPi0_Reco_noRichInd_chi_vs_pt", "fhPi0_Reco_noRichInd_chi_vs_pt;pt [GeV];#chi^{2} of momentum fit", 200,
690 0., 10., 1000, 0., 100.);
693 new TH2D("fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0",
694 "fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0;pt [GeV];#chi^{2} of momentum fit", 200, 0., 10., 1000, 0., 100.);
697 new TH2D("fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target",
698 "fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target;pt [GeV];#chi^{2} "
699 "of momentum fit",
700 200, 0., 10., 1000, 0., 100.);
703 new TH2D("fhPi0_Reco_noRichInd_chi_vs_pt_eRest",
704 "fhPi0_Reco_noRichInd_chi_vs_pt_eRest;pt [GeV];#chi^{2} of momentum fit", 200, 0., 10., 1000, 0., 100.);
706
707
708 // ################################################
709
711 new TH2D("fhPi0_Reco_chi_vs_momentum", "fhPi0_Reco_chi_vs_momentum;momentum [GeV];#chi^{2} of momentum fit", 200,
712 0., 10., 1000, 0., 100.);
714 fhPi0_Reco_chi_vs_momentum_eFromPi0 = new TH2D("fhPi0_Reco_chi_vs_momentum_eFromPi0",
715 "fhPi0_Reco_chi_vs_momentum_eFromPi0;momentum [GeV];#chi^{2} of "
716 "momentum fit",
717 200, 0., 10., 1000, 0., 100.);
719 fhPi0_Reco_chi_vs_momentum_eFromPi0_Target = new TH2D("fhPi0_Reco_chi_vs_momentum_eFromPi0_Target",
720 "fhPi0_Reco_chi_vs_momentum_eFromPi0_Target;momentum "
721 "[GeV];#chi^{2} of momentum fit",
722 200, 0., 10., 1000, 0., 100.);
724
725 fhPi0_Reco_chi_vs_pt = new TH2D("fhPi0_Reco_chi_vs_pt", "fhPi0_Reco_chi_vs_pt;pt [GeV];#chi^{2} of momentum fit", 200,
726 0., 10., 1000, 0., 100.);
729 new TH2D("fhPi0_Reco_chi_vs_pt_eFromPi0", "fhPi0_Reco_chi_vs_pt_eFromPi0;pt [GeV];#chi^{2} of momentum fit", 200,
730 0., 10., 1000, 0., 100.);
733 new TH2D("fhPi0_Reco_chi_vs_pt_eFromPi0_Target",
734 "fhPi0_Reco_chi_vs_pt_eFromPi0_Target;pt [GeV];#chi^{2} of momentum fit", 200, 0., 10., 1000, 0., 100.);
736}
737
738
740{
741 timer_exec.Start();
742 timer_all.Start();
743
744
745 cout << "=======================================================================" << endl;
746 cout << "========== CbmAnaConversion, event No. " << fEventNum << endl;
747 cout << "=======================================================================" << endl;
748
749 fEventNum++;
750
753
754
755 // arrays of tracks
756 fMCTracklist.clear();
757 fMCTracklist_all.clear();
758 fRecoTracklist.clear();
759 fRecoTracklistEPEM.clear();
760 fRecoTracklistEPEM_id.clear();
763 fRecoMomentum.clear();
764 fRecoRefittedMomentum.clear();
765
766 fTestTracklist.clear();
768 fTestTracklist_chi.clear();
770 fTestTracklist_ndf.clear();
772
781
782 // several counters
783 int countPrimEl = 0;
784 int countSecEl = 0;
785 int countAllEl = 0;
786 // int countGammaEl = 0;
787 // int countMothers = 0;
788 int countPrimPart = 0;
789 int countPi0MC = 0;
790 int countPi0MC_cut = 0;
791 int countPi0MC_fromPrimary = 0;
792 int countPi0MC_reconstructible = 0;
793 int countEtaMC = 0;
794 int countEtaMC_cut = 0;
795 int countEtaMC_fromPrimary = 0;
796
797 if (fPrimVertex != NULL) { fKFVertex = CbmKFVertex(*fPrimVertex); }
798 else {
799 Fatal("CbmAnaConversion::Exec", "No PrimaryVertex array!");
800 }
801
802 if (DoKFAnalysis) {
803 // fAnaKF->SetSignalIds(fKFparticleFinderQA->GetSignalIds());
804 // fAnaKF->SetGhostIds(fKFparticleFinderQA->GetGhostIds());
805 fAnaKF->Exec();
806 }
807
809
810 if (DoPhotons) { fAnaPhotons->Exec(); }
811
812 if (DoPhotons2) { fAnaPhotons2->Exec(); }
813
814 if (DoRecoFull) { fAnaRecoFull->Exec(); }
815
816 if (DoTomography) {
817 fAnaTomography->Exec(); // analyse gamma-conversions with MC data
818 }
819
820 if (DoTest) {
821 fAnaTest->Exec();
822 fAnaTest2->Exec();
823 }
824
825 // ========================================================================================
826 // START - Analyse MC tracks
827 timer_mc.Start();
828
829 Int_t nofMcTracks = fMcTracks->GetEntriesFast();
830 fhNofTracks_mctrack->Fill(nofMcTracks);
831 for (int i = 0; i < nofMcTracks; i++) {
832 CbmMCTrack* mctrack = (CbmMCTrack*) fMcTracks->At(i);
833 if (mctrack == NULL) continue;
834
835 FillMCTracklists(mctrack, i); // fill tracklists for further analyses
836
837 fhPdgCodes->Fill(TMath::Abs(mctrack->GetPdgCode()));
838
839 if (mctrack->GetMotherId() == -1) { countPrimPart++; }
840 if (mctrack->GetMotherId() == -1 && TMath::Abs(mctrack->GetPdgCode()) == 11) { countPrimEl++; }
841 if (mctrack->GetMotherId() != -1 && TMath::Abs(mctrack->GetPdgCode()) == 11) { countSecEl++; }
842 if (TMath::Abs(mctrack->GetPdgCode()) == 11) { countAllEl++; }
843
844
845 if (mctrack->GetPdgCode() == 111) { // particle is pi0
846 countPi0MC++;
847 TVector3 v;
848 mctrack->GetStartVertex(v);
849 if (v.Z() <= 4) { countPi0MC_cut++; }
850 fhPi0_z->Fill(v.Z());
851 //Double_t r2 = v.Z()*v.Z() * tan(25./180*TMath::Pi()) * tan(25./180*TMath::Pi());
852 //if( (v.X()*v.X() + v.Y()*v.Y()) <= r2) {
853 // fhPi0_z_cut->Fill(v.Z());
854 //}
855
856
857 TVector3 momentum;
858 mctrack->GetMomentum(momentum);
859
860 int motherId = mctrack->GetMotherId();
861 if (motherId == -1) {
862 countPi0MC_fromPrimary++;
863 fhPi0_pt->Fill(mctrack->GetPt());
864 fhPi0_pt_vs_rap->Fill(mctrack->GetPt(), mctrack->GetRapidity());
865 fhPi0_theta->Fill(momentum.Theta() * 180. / TMath::Pi());
866 fhPi0_theta_vs_rap->Fill(momentum.Theta() * 180. / TMath::Pi(), mctrack->GetRapidity());
867 fhPi0_z_cut->Fill(v.Z());
868 }
869
870 Bool_t reconstructible = AnalysePi0_MC(mctrack, i);
871 if (reconstructible) countPi0MC_reconstructible++;
872 }
873
874
875 if (TMath::Abs(mctrack->GetPdgCode()) == 221) { // particle is eta
876 countEtaMC++;
877 TVector3 v, momentum;
878 mctrack->GetStartVertex(v);
879 mctrack->GetMomentum(momentum);
880
881 if (v.Z() <= 4) { countEtaMC_cut++; }
882 int motherId = mctrack->GetMotherId();
883 if (motherId == -1) {
884 countEtaMC_fromPrimary++;
885 fhEta_pt->Fill(mctrack->GetPt());
886 fhEta_pt_vs_rap->Fill(mctrack->GetPt(), mctrack->GetRapidity());
887 fhEta_theta->Fill(momentum.Theta() * 180. / TMath::Pi());
888 fhEta_theta_vs_rap->Fill(momentum.Theta() * 180. / TMath::Pi(), mctrack->GetRapidity());
889 }
890 }
891
892
893 if (TMath::Abs(mctrack->GetPdgCode()) == 113
894 || TMath::Abs(mctrack->GetPdgCode()) == 213) { // particle is rho(770)^0 or rho^+/-
895 TVector3 v, momentum;
896 mctrack->GetStartVertex(v);
897 mctrack->GetMomentum(momentum);
898
899 //if (v.Z() <= 1) {
900 fhRho_pt->Fill(mctrack->GetPt());
901 fhRho_pt_vs_rap->Fill(mctrack->GetPt(), mctrack->GetRapidity());
902 fhRho_theta->Fill(momentum.Theta() * 180. / TMath::Pi());
903 fhRho_theta_vs_rap->Fill(momentum.Theta() * 180. / TMath::Pi(), mctrack->GetRapidity());
904 //}
905 }
906
907
908 if (TMath::Abs(mctrack->GetPdgCode()) == 223
909 || TMath::Abs(mctrack->GetPdgCode()) == 3334) { // particle is omega(782)^0 or omega^+/-
910 TVector3 v, momentum;
911 mctrack->GetStartVertex(v);
912 mctrack->GetMomentum(momentum);
913
914 //if (v.Z() <= 1) {
915 fhOmega_pt->Fill(mctrack->GetPt());
916 fhOmega_pt_vs_rap->Fill(mctrack->GetPt(), mctrack->GetRapidity());
917 fhOmega_theta->Fill(momentum.Theta() * 180. / TMath::Pi());
918 fhOmega_theta_vs_rap->Fill(momentum.Theta() * 180. / TMath::Pi(), mctrack->GetRapidity());
919 //}
920 }
921
922 if (TMath::Abs(mctrack->GetPdgCode()) == 11) { // particle is electron
923 AnalyseElectrons(mctrack);
924 }
925 }
926
927
928 cout << "CbmAnaConversion::Exec - Number of pi0 in MC sample: " << countPi0MC << endl;
929 cout << "CbmAnaConversion::Exec - Number of pi0 from primary: " << countPi0MC_fromPrimary << endl;
930 fhNofPi0_perEvent->Fill(countPi0MC);
931 fhNofPi0_perEvent_cut->Fill(countPi0MC_cut);
932 fhNofPi0_perEvent_cut2->Fill(countPi0MC_fromPrimary);
933 fhNofPi0_perEvent_cut3->Fill(countPi0MC_reconstructible);
934 fhNofEta_perEvent->Fill(countEtaMC);
935 fhNofEta_perEvent_cut->Fill(countEtaMC_cut);
936 fhNofEta_perEvent_cut2->Fill(countEtaMC_fromPrimary);
937
938 fNofGeneratedPi0 = countPi0MC_fromPrimary;
942
944
945 if (DoReconstruction) {
947 //fAnaReco->InvariantMassMC_all();
948 }
949
951
952
953 fhNofElPrim->Fill(countPrimEl);
954 fhNofElSec->Fill(countSecEl);
955 fhNofElAll->Fill(countAllEl);
956
957
958 // END - Analyse MC tracks
959 // ========================================================================================
960 timer_mc.Stop();
961 fTime_mc += timer_mc.RealTime();
962
963
964 // ========================================================================================
965 // START - Analyse reconstructed tracks
966 timer_rec.Start();
967
968 Int_t nofElectrons4epem = 0;
969
970 Int_t ngTracks = fGlobalTracks->GetEntriesFast();
971 fhNofTracks_globaltrack->Fill(ngTracks);
972 for (Int_t iGTrack = 0; iGTrack < ngTracks; iGTrack++) {
973 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGTrack);
974 if (NULL == gTrack) continue;
975 int stsInd = gTrack->GetStsTrackIndex();
976 int richInd = gTrack->GetRichRingIndex();
977 // int trdInd = gTrack->GetTrdTrackIndex();
978 // int tofInd = gTrack->GetTofHitIndex();
979
980 if (stsInd < 0) continue;
981 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
982 if (stsTrack == NULL) continue;
983 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
984 if (stsMatch == NULL) continue;
985 if (stsMatch->GetNofLinks() <= 0) continue;
986 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
987 if (stsMcTrackId < 0) continue;
988 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
989 if (mcTrack1 == NULL) continue;
990
991
992 // calculate refitted momenta at primary vertex
993 TVector3 refittedMomentum;
994 CbmL1PFFitter fPFFitter;
995 vector<CbmStsTrack> stsTracks;
996 stsTracks.resize(1);
997 stsTracks[0] = *stsTrack;
998 vector<CbmL1PFFitter::PFFieldRegion> vField;
999 vector<float> chiPrim;
1000 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
1001 //cand.chi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
1002 //cand.chi2Prim = chiPrim[0];
1003 const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();
1004 vtxTrack->Momentum(refittedMomentum);
1005 float result_chi = chiPrim[0];
1006
1007
1008 // Doing refit of momenta with electron assumption
1009 CbmL1PFFitter fPFFitter_electron;
1010 vector<CbmStsTrack> stsTracks_electron;
1011 stsTracks_electron.resize(1);
1012 stsTracks_electron[0] = *stsTrack;
1013 vector<CbmL1PFFitter::PFFieldRegion> vField_electron;
1014 vector<float> chiPrim_electron;
1015 vector<int> pidHypo_electron;
1016 pidHypo_electron.push_back(11);
1017 fPFFitter_electron.Fit(stsTracks_electron, pidHypo_electron);
1018 fPFFitter_electron.GetChiToVertex(stsTracks_electron, vField_electron, chiPrim_electron, fKFVertex, 3e6);
1019
1020 TVector3 refittedMomentum_electron;
1021 const FairTrackParam* vtxTrack_electron = stsTracks_electron[0].GetParamFirst();
1022 vtxTrack_electron->Momentum(refittedMomentum_electron);
1023 float result_chi_electron = chiPrim_electron[0];
1024 float result_ndf_electron = stsTracks_electron[0].GetNDF();
1025
1026
1027 Double_t startvertexZ = vtxTrack_electron->GetZ();
1028 fhPi0_Reco_ndf->Fill(result_ndf_electron);
1029 fhPi0_Reco_chi->Fill(result_chi_electron);
1030 fhPi0_Reco_ndf_vs_chi->Fill(result_ndf_electron, result_chi_electron);
1031 fhPi0_Reco_ndf_vs_startvertex->Fill(result_ndf_electron, startvertexZ);
1032 fhPi0_Reco_startvertex_vs_chi->Fill(startvertexZ, result_chi_electron);
1033
1034 Double_t nofhits_sts = stsTrack->GetTotalNofHits();
1035 fhPi0_Reco_startvertex_vs_nofhits->Fill(startvertexZ, nofhits_sts);
1036
1037 fhPi0_Reco_noRichInd_chi_vs_momentum->Fill(refittedMomentum_electron.Mag(), result_chi_electron);
1038 fhPi0_Reco_noRichInd_chi_vs_pt->Fill(refittedMomentum_electron.Perp(), result_chi_electron);
1039
1040 fTestTracklist_noRichInd.push_back(mcTrack1);
1041 fTestTracklist_noRichInd_MCindex.push_back(stsMcTrackId);
1042 fTestTracklist_noRichInd_momentum.push_back(refittedMomentum_electron);
1043 fTestTracklist_noRichInd_chi.push_back(result_chi_electron);
1044 fTestTracklist_noRichInd_richInd.push_back(richInd);
1045 fTestTracklist_noRichInd_gTrackId.push_back(iGTrack);
1046 fTestTracklist_noRichInd_ndf.push_back(result_ndf_electron);
1047 fTestTracklist_noRichInd_nofhits.push_back(nofhits_sts);
1048
1049
1050 if (richInd < 0) continue;
1051
1052 // ANN output for rings
1053 Int_t pdg = mcTrack1->GetPdgCode();
1054 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richInd));
1055 if (NULL != ring) {
1056 Double_t ann = CbmRichElectronIdAnn::GetInstance().CalculateAnnValue(iGTrack, refittedMomentum_electron.Mag());
1057
1058 if (TMath::Abs(pdg) == 11) {
1059 fhANN_output_electrons->Fill(ann);
1060 if (result_chi_electron <= CbmAnaConversionCutSettings::CalcChiCut(refittedMomentum_electron.Perp())) {
1062 }
1063 }
1064 if (TMath::Abs(pdg) != 11) {
1065 fhANN_output_else->Fill(ann);
1066 if (result_chi_electron <= CbmAnaConversionCutSettings::CalcChiCut(refittedMomentum_electron.Perp())) {
1067 fhANN_output_else_chiCut->Fill(ann);
1068 }
1069 }
1070 }
1071
1072
1073 TVector3 stsMomentumVec; // momenta as measured by STS
1074 stsTrack->GetParamFirst()->Momentum(stsMomentumVec);
1075 Double_t stsMomentum = stsMomentumVec.Mag();
1076
1077 TVector3 mcMomentumVec; // momenta from true-MC data
1078 mcTrack1->GetMomentum(mcMomentumVec);
1079 Double_t mcMomentum = mcMomentumVec.Mag();
1080 fhMomentum_MCvsReco->Fill(mcMomentum, stsMomentum);
1081 fhMomentum_MCvsReco_diff->Fill(TMath::Abs(mcMomentum - stsMomentum) / mcMomentum);
1082
1083 TVector3 bothtogether; // combination of measured (STS) momenta and MC momenta
1084 bothtogether.SetX(mcMomentumVec.X());
1085 bothtogether.SetY(stsMomentumVec.Y());
1086 bothtogether.SetZ(stsMomentumVec.Z());
1087
1088
1089 // Fill tracklists containing momenta from mc-true, measured in sts, refitted at primary
1090 Bool_t isFilled =
1091 FillRecoTracklistEPEM(mcTrack1, stsMomentumVec, refittedMomentum, stsMcTrackId, result_chi, iGTrack);
1092 if (isFilled) nofElectrons4epem++;
1093
1094
1095 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
1096 if (richMatch == NULL) continue;
1097 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
1098 if (richMcTrackId < 0) continue;
1099 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
1100 if (mcTrack2 == NULL) continue;
1101
1102 //if(stsMcTrackId != richMcTrackId) continue;
1103
1104
1105 if (stsMcTrackId == richMcTrackId) {
1106 //CbmRichRing* ring = static_cast<CbmRichRing*> (fRichRings->At(richInd));
1107 if (NULL != ring) {
1108 Double_t ann = CbmRichElectronIdAnn::GetInstance().CalculateAnnValue(iGTrack, refittedMomentum_electron.Mag());
1109 if (TMath::Abs(pdg) == 11) {
1110 fhANN_output_electrons2->Fill(ann);
1111 fhANN_output_electrons_vs_p->Fill(refittedMomentum_electron.Mag(), ann);
1112 }
1113 if (TMath::Abs(pdg) != 11) {
1114 fhANN_output_else2->Fill(ann);
1115 fhANN_output_else_vs_p->Fill(refittedMomentum_electron.Mag(), ann);
1116 }
1117 }
1118 }
1119
1120
1121 //int pdg = TMath::Abs(mcTrack1->GetPdgCode());
1122 //int motherId = mcTrack1->GetMotherId();
1123 //double momentum = mcTrack1->GetP();
1124 stsMatch->GetTrueOverAllHitsRatio();
1125
1126
1127 if (DoTomography) { fAnaTomography->TomographyReco(mcTrack1); }
1128 FillRecoTracklist(mcTrack1);
1129
1130
1131 fTestTracklist.push_back(mcTrack1);
1132 fTestTracklist_momentum.push_back(refittedMomentum_electron);
1133 fTestTracklist_chi.push_back(result_chi_electron);
1134 fTestTracklist_richInd.push_back(richInd);
1135 fTestTracklist_ndf.push_back(result_ndf_electron);
1136 fTestTracklist_nofhits.push_back(nofhits_sts);
1137
1138 fhPi0_Reco_chi_vs_pt->Fill(refittedMomentum_electron.Perp(), result_chi_electron);
1139 fhPi0_Reco_chi_vs_momentum->Fill(refittedMomentum_electron.Mag(), result_chi_electron);
1140 }
1141
1144
1145
1146 fhNofElectrons_4epem->Fill(nofElectrons4epem);
1147
1148
1149 // InvariantMassTestReco();
1150
1151 if (DoReconstruction) {
1156 }
1157
1158 // END - analyse reconstructed tracks
1159 // ========================================================================================
1160 timer_rec.Stop();
1161 fTime_rec += timer_rec.RealTime();
1162
1163
1164 // =========================================================================================================================
1165 // ============================================== END - EXEC function ======================================================
1166 // =========================================================================================================================
1167 timer_exec.Stop();
1168 fTime_exec += timer_exec.RealTime();
1169 timer_all.Stop();
1170 fTime_all += timer_all.RealTime();
1171}
1172
1173
1175{
1176 timer_all.Start();
1177
1178 cout << "\n\n############### CALLING FINISH ROUTINES... ############" << endl;
1179
1180
1181 // Write histograms to a file
1182 gDirectory->mkdir("analysis-conversion");
1183 gDirectory->cd("analysis-conversion");
1184
1185
1186 gDirectory->mkdir("KFParticle");
1187 gDirectory->cd("KFParticle");
1188 for (UInt_t i = 0; i < fHistoList_kfparticle.size(); i++) {
1189 fHistoList_kfparticle[i]->Write();
1190 }
1191 gDirectory->cd("..");
1192
1193
1195 if (DoRichAnalysis) { fAnaRich->Finish(); }
1196 if (DoKFAnalysis) { fAnaKF->Finish(); }
1197 if (DoReconstruction) { fAnaReco->Finish(); }
1198 if (DoRecoFull) { fAnaRecoFull->Finish(); }
1199 if (DoPhotons) { fAnaPhotons->Finish(); }
1200 if (DoPhotons2) { fAnaPhotons2->Finish(); }
1201 if (DoTest) { fAnaTest->Finish(); }
1202 if (DoTest) { fAnaTest2->Finish(); }
1203
1204
1205 gDirectory->mkdir("further analyses");
1206 gDirectory->cd("further analyses");
1207 for (UInt_t i = 0; i < fHistoList_furtherAnalyses.size(); i++) {
1208 fHistoList_furtherAnalyses[i]->Write();
1209 }
1210 gDirectory->cd("..");
1211
1212 for (UInt_t i = 0; i < fHistoList.size(); i++) {
1213 fHistoList[i]->Write();
1214 }
1215 gDirectory->cd("..");
1216
1217
1218 timer_all.Stop();
1219 fTime_all += timer_all.RealTime();
1220
1221 cout << endl;
1222 cout << "############### FINISHED MAIN TASK ##############" << endl;
1223 cout << "Particlecounter: " << particlecounter << endl;
1224 cout << "Particlecounter (2 daughters): " << particlecounter_2daughters << endl;
1225 cout << "Particlecounter (3 daughters): " << particlecounter_3daughters << endl;
1226 cout << "Particlecounter (4 daughters): " << particlecounter_4daughters << endl;
1227 cout << "Particlecounter_all: " << particlecounter_all << endl;
1228 cout << "#####################################" << endl;
1229 cout << "Number of generated pi0 (all events): " << fNofGeneratedPi0_allEvents << endl;
1230 cout << "Number of reconstructed pi0 (all events): " << fNofPi0_kfparticle_allEvents
1231 << "\t - fraction: " << 1.0 * fNofPi0_kfparticle_allEvents / fNofGeneratedPi0_allEvents << endl;
1232 cout << "#####################################" << endl;
1233 cout << "############### OVERALL TIMERS ###############" << endl;
1234 cout << std::fixed;
1235 cout << std::setprecision(1);
1236 cout << "Complete time: " << fTime_all << endl;
1237 cout << "Exec time: " << fTime_exec << endl;
1238 cout << "MC time: " << fTime_mc << "\t RECO time: " << fTime_rec << endl;
1239 cout << "############### ############## ###############" << endl;
1240 cout << "Number of events in fhNofPi0_perEvent histogram: " << fhNofPi0_perEvent->GetEntries() << endl;
1241 cout << "############### ############## ###############" << endl;
1242 // =========================================================================================================================
1243 // ============================================== END - FINISH function ====================================================
1244 // =========================================================================================================================
1245}
1246
1247
1249{
1250 int motherId = mctrack->GetMotherId();
1251 if (motherId == -1) return;
1252 CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
1253 int mcMotherPdg = -1;
1254 if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
1255 //cout << mcMotherPdg << endl;
1256
1257 if (mcMotherPdg == 22) {
1258 fhElectronSources->Fill(0);
1259 int grandmotherId = mother->GetMotherId();
1260 if (grandmotherId == -1) return;
1261 CbmMCTrack* grandmother = (CbmMCTrack*) fMcTracks->At(grandmotherId);
1262 int mcGrandmotherPdg = -1;
1263 if (NULL != grandmother) mcGrandmotherPdg = grandmother->GetPdgCode();
1264 if (mcGrandmotherPdg == 111) fhElectronSources->Fill(4);
1265 if (mcGrandmotherPdg == 221) fhElectronSources->Fill(5);
1266
1267 if (mcGrandmotherPdg == 111) {
1268 TVector3 v;
1269 mctrack->GetStartVertex(v);
1270 fhElectronsFromPi0_z->Fill(v.Z());
1271 }
1272 }
1273 if (mcMotherPdg == 111) fhElectronSources->Fill(1);
1274 if (mcMotherPdg == 221) fhElectronSources->Fill(2);
1275 if (mcMotherPdg != 22 && mcMotherPdg != 111 && mcMotherPdg != 221) fhElectronSources->Fill(3);
1276
1277 if (mcMotherPdg == 22) {
1278 TVector3 v;
1279 mctrack->GetStartVertex(v);
1280 // fhGammaZ->Fill(v.Z());
1281 // countGammaEl++;
1282 }
1283}
1284
1285
1287{
1288 LmvmKinePar params;
1289
1290 TVector3 momP; //momentum e+
1291 mctrackP->GetMomentum(momP);
1292 Double_t energyP = TMath::Sqrt(momP.Mag2() + M2E);
1293 TLorentzVector lorVecP(momP, energyP);
1294
1295 TVector3 momM; //momentum e-
1296 mctrackM->GetMomentum(momM);
1297 Double_t energyM = TMath::Sqrt(momM.Mag2() + M2E);
1298 TLorentzVector lorVecM(momM, energyM);
1299
1300 TVector3 momPair = momP + momM;
1301 Double_t energyPair = energyP + energyM;
1302 Double_t ptPair = momPair.Perp();
1303 Double_t pzPair = momPair.Pz();
1304 Double_t yPair = 0.5 * TMath::Log((energyPair + pzPair) / (energyPair - pzPair));
1305 Double_t anglePair = lorVecM.Angle(lorVecP.Vect());
1306 Double_t theta = 180. * anglePair / TMath::Pi();
1307 Double_t minv = 2. * TMath::Sin(anglePair / 2.) * TMath::Sqrt(momM.Mag() * momP.Mag());
1308
1309 params.fMomentumMag = momPair.Mag();
1310 params.fPt = ptPair;
1311 params.fRapidity = yPair;
1312 params.fMinv = minv;
1313 params.fAngle = theta;
1314 return params;
1315}
1316
1317Double_t CbmAnaConversion::Invmass_2gammas(const CbmMCTrack* mctrack1, const CbmMCTrack* mctrack2)
1318// calculation of invariant mass from two gammas (m=0)
1319{
1320 TVector3 mom1;
1321 mctrack1->GetMomentum(mom1);
1322 Double_t energy1 = TMath::Sqrt(mom1.Mag2());
1323 TLorentzVector lorVec1(mom1, energy1);
1324
1325 TVector3 mom2;
1326 mctrack2->GetMomentum(mom2);
1327 Double_t energy2 = TMath::Sqrt(mom2.Mag2());
1328 TLorentzVector lorVec2(mom2, energy2);
1329
1330 TLorentzVector sum = lorVec1 + lorVec2;
1331 return sum.Mag();
1332}
1333
1334
1335Double_t CbmAnaConversion::Invmass_2particles(const CbmMCTrack* mctrack1, const CbmMCTrack* mctrack2)
1336// calculation of invariant mass from two electrons/positrons
1337{
1338 TVector3 mom1;
1339 mctrack1->GetMomentum(mom1);
1340 Double_t energy1 = TMath::Sqrt(mom1.Mag2() + M2E);
1341 TLorentzVector lorVec1(mom1, energy1);
1342
1343 TVector3 mom2;
1344 mctrack2->GetMomentum(mom2);
1345 Double_t energy2 = TMath::Sqrt(mom2.Mag2() + M2E);
1346 TLorentzVector lorVec2(mom2, energy2);
1347
1348 TLorentzVector sum = lorVec1 + lorVec2;
1349 return sum.Mag();
1350}
1351
1352
1353Double_t CbmAnaConversion::Invmass_4particles(const CbmMCTrack* mctrack1, const CbmMCTrack* mctrack2,
1354 const CbmMCTrack* mctrack3, const CbmMCTrack* mctrack4)
1355// calculation of invariant mass from four electrons/positrons
1356{
1357 /* TVector3 mom1;
1358 mctrack1->GetMomentum(mom1);
1359 TVector3 tempmom1;
1360 tempmom1.SetX(SmearValue(mom1.X()));
1361 tempmom1.SetY(SmearValue(mom1.Y()));
1362 tempmom1.SetZ(SmearValue(mom1.Z()));
1363 Double_t energy1 = TMath::Sqrt(tempmom1.Mag2() + M2E);
1364 TLorentzVector lorVec1(tempmom1, energy1);
1365
1366 TVector3 mom2;
1367 mctrack2->GetMomentum(mom2);
1368 TVector3 tempmom2;
1369 tempmom2.SetX(SmearValue(mom2.X()));
1370 tempmom2.SetY(SmearValue(mom2.Y()));
1371 tempmom2.SetZ(SmearValue(mom2.Z()));
1372 Double_t energy2 = TMath::Sqrt(tempmom2.Mag2() + M2E);
1373 TLorentzVector lorVec2(tempmom2, energy2);
1374
1375 TVector3 mom3;
1376 mctrack3->GetMomentum(mom3);
1377 TVector3 tempmom3;
1378 tempmom3.SetX(SmearValue(mom3.X()));
1379 tempmom3.SetY(SmearValue(mom3.Y()));
1380 tempmom3.SetZ(SmearValue(mom3.Z()));
1381 Double_t energy3 = TMath::Sqrt(tempmom3.Mag2() + M2E);
1382 TLorentzVector lorVec3(tempmom3, energy3);
1383
1384 TVector3 mom4;
1385 mctrack4->GetMomentum(mom4);
1386 TVector3 tempmom4;
1387 tempmom4.SetX(SmearValue(mom4.X()));
1388 tempmom4.SetY(SmearValue(mom4.Y()));
1389 tempmom4.SetZ(SmearValue(mom4.Z()));
1390 Double_t energy4 = TMath::Sqrt(tempmom4.Mag2() + M2E);
1391 TLorentzVector lorVec4(tempmom4, energy4);
1392*/
1393 TLorentzVector lorVec1;
1394 mctrack1->Get4Momentum(lorVec1);
1395
1396 TLorentzVector lorVec2;
1397 mctrack2->Get4Momentum(lorVec2);
1398
1399 TLorentzVector lorVec3;
1400 mctrack3->Get4Momentum(lorVec3);
1401
1402 TLorentzVector lorVec4;
1403 mctrack4->Get4Momentum(lorVec4);
1404
1405
1406 TLorentzVector sum;
1407 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
1408 cout << "mc: \t" << sum.Px() << " / " << sum.Py() << " / " << sum.Pz() << " / " << sum.E()
1409 << "\t => mag = " << sum.Mag() << endl;
1410 return sum.Mag();
1411}
1412
1413
1414Double_t CbmAnaConversion::Invmass_4particlesRECO(const TVector3 part1, const TVector3 part2, const TVector3 part3,
1415 const TVector3 part4)
1416// calculation of invariant mass from four electrons/positrons
1417{
1418 Double_t energy1 = TMath::Sqrt(part1.Mag2() + M2E);
1419 TLorentzVector lorVec1(part1, energy1);
1420
1421 Double_t energy2 = TMath::Sqrt(part2.Mag2() + M2E);
1422 TLorentzVector lorVec2(part2, energy2);
1423
1424 Double_t energy3 = TMath::Sqrt(part3.Mag2() + M2E);
1425 TLorentzVector lorVec3(part3, energy3);
1426
1427 Double_t energy4 = TMath::Sqrt(part4.Mag2() + M2E);
1428 TLorentzVector lorVec4(part4, energy4);
1429
1430 TLorentzVector sum;
1431 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
1432
1433 return sum.Mag();
1434}
1435
1436
1437// =========================================================================================================================
1438// ============================================== UP TO HERE: Invariant mass calculation ===================================
1439// =========================================================================================================================
1440
1441
1443// fill all relevant tracklists containing MC tracks
1444{
1445 Bool_t electrons = true;
1446 Bool_t gammas = true;
1447
1448 // fill gamma tracklist
1449 if (TMath::Abs(mctrack->GetPdgCode()) == 22 && gammas) {
1450 int motherId = mctrack->GetMotherId();
1451 if (motherId != -1) {
1452 CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
1453 int mcMotherPdg = -1;
1454 if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
1455 if (mcMotherPdg == 111) { // pdg code 111 = pi0
1456 fMCTracklist.push_back(mctrack);
1457 }
1458 }
1459 }
1460
1461 // fill electron tracklists
1462 if (TMath::Abs(mctrack->GetPdgCode()) == 11 && electrons) {
1463 TVector3 v;
1464 mctrack->GetStartVertex(v);
1465 if (v.Z() <= 70) {
1466 int motherId = mctrack->GetMotherId();
1467 if (motherId >= 0) {
1468 CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
1469 int mcMotherPdg = -1;
1470 if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
1471 if (mcMotherPdg == 22 || mcMotherPdg == 111 || mcMotherPdg == 221) {
1472 int grandmotherId = mother->GetMotherId();
1473 if (grandmotherId >= 0 || true) {
1474 // CbmMCTrack* grandmother = (CbmMCTrack*) fMcTracks->At(grandmotherId);
1475 // int mcGrandmotherPdg = -1;
1476 // if (NULL != grandmother) mcGrandmotherPdg = grandmother->GetPdgCode();
1477 // if(mcGrandmotherPdg == 111) {
1478 fMCTracklist_all.push_back(mctrack);
1479 // }
1480 }
1481 }
1482 }
1483 }
1484 }
1485}
1486
1487
1489{
1490 if (TMath::Abs(mctrack->GetPdgCode()) == 11) {
1491 int motherId = mctrack->GetMotherId();
1492 if (motherId != -1) {
1493 CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(motherId);
1494 int mcMotherPdg = -1;
1495 if (NULL != mother) mcMotherPdg = mother->GetPdgCode();
1496 if (mcMotherPdg == 111 || mcMotherPdg == 22) { // pdg code 111 = pi0, 22 = gamma
1497 fRecoTracklist.push_back(mctrack);
1498 // cout << "pdg " << mctrack->GetPdgCode() << "\t motherid " << motherId << endl;
1499 test++;
1500 }
1501 }
1502 }
1503}
1504
1505
1506Bool_t CbmAnaConversion::FillRecoTracklistEPEM(CbmMCTrack* mctrack, TVector3 stsMomentum, TVector3 refittedMom, int i,
1507 Double_t chi, Int_t GlobalTrackId)
1508{
1509 Bool_t isFilled = false;
1510 if (TMath::Abs(mctrack->GetPdgCode()) == 11) {
1511 int motherId = mctrack->GetMotherId();
1512 if (motherId != -1) {
1513 fRecoTracklistEPEM.push_back(mctrack);
1514 fRecoTracklistEPEM_id.push_back(i);
1515 fRecoTracklistEPEM_chi.push_back(chi);
1516 fRecoTracklistEPEM_gtid.push_back(GlobalTrackId);
1517 fRecoMomentum.push_back(stsMomentum);
1518 fRecoRefittedMomentum.push_back(refittedMom);
1519 isFilled = true;
1520 }
1521 }
1522 return isFilled;
1523}
1524
1525
1527{
1528 for (unsigned int i = 0; i < fMCTracklist_all.size(); i++) {
1529 for (unsigned int j = i + 1; j < fMCTracklist_all.size(); j++) {
1530 // Double_t invmass = Invmass_2particles(fMCTracklist_all[i],fMCTracklist_all[j]);
1531
1532 int motherId_i = fMCTracklist_all[i]->GetMotherId();
1533 int motherId_j = fMCTracklist_all[j]->GetMotherId();
1534
1535 CbmMCTrack* mother_i = (CbmMCTrack*) fMcTracks->At(motherId_i);
1536 int mcMotherPdg_i = -1;
1537 if (NULL != mother_i) mcMotherPdg_i = mother_i->GetPdgCode();
1538
1539 // CbmMCTrack* mother_j = (CbmMCTrack*) fMcTracks->At(motherId_j);
1540 // int mcMotherPdg_j = -1;
1541 // if (NULL != mother_j) mcMotherPdg_j = mother_j->GetPdgCode();
1542
1543 if (motherId_i == motherId_j
1544 && ((fMCTracklist_all[i]->GetPdgCode() == 11 && fMCTracklist_all[j]->GetPdgCode() == -11)
1545 || (fMCTracklist_all[i]->GetPdgCode() == -11 && fMCTracklist_all[j]->GetPdgCode() == 11))) {
1546 // fhInvariantMass_MC_all->Fill(invmass);
1547 //cout << "e+e- decay detected! MotherId " << motherId_i << "\t invariant mass: " << invmass << endl;
1548
1549 if (mcMotherPdg_i == 111) {
1550 //fhInvariantMass_MC_pi0->Fill(invmass);
1551 cout << "-- mother particle of decay: pi0" << endl;
1552 }
1553 }
1554 }
1555 }
1556}
1557
1558
1560// calculation of invariant mass via pi0-> gamma gamma, ONLY FROM MC DATA!
1561{
1562 for (unsigned int i = 0; i < fMCTracklist.size(); i++) {
1563 for (unsigned int j = i + 1; j < fMCTracklist.size(); j++) {
1564
1565 //if(fMCTracklist[i]->GetPx() != fMCTracklist[j]->GetPx() && fMCTracklist[i]->GetPy() != fMCTracklist[j]->GetPy()) {
1566 Double_t invmass = Invmass_2gammas(fMCTracklist[i], fMCTracklist[j]);
1567 fhInvariantMass_test->Fill(invmass);
1568 TVector3 vi;
1569 fMCTracklist[i]->GetStartVertex(vi);
1570 TVector3 vj;
1571 fMCTracklist[j]->GetStartVertex(vj);
1572
1573 int motherId_i = fMCTracklist[i]->GetMotherId();
1574 int motherId_j = fMCTracklist[j]->GetMotherId();
1575
1576 if (motherId_i == motherId_j) {
1577 fhInvariantMass_test2->Fill(invmass);
1578
1579 if (invmass < 0.001) {
1580 // cout << "gamma1 " << fMCTracklist[i]->GetPx() << "\t" << fMCTracklist[i]->GetPy() << "\t" << fMCTracklist[i]->GetPz() << endl;
1581 // cout << "gamma2 " << fMCTracklist[j]->GetPx() << "\t" << fMCTracklist[j]->GetPy() << "\t" << fMCTracklist[j]->GetPz() << endl;
1582 }
1583 }
1584 if (vi.Z() == vj.Z() && vi.Y() == vj.Y() && vi.X() == vj.X()) { fhInvariantMass_test3->Fill(invmass); }
1585 //}
1586 }
1587 }
1588}
1589
1590
1592// calculation of invariant mass via pi0 -> .. -> e+ e- e+ e-, ONLY FROM RECONSTRUCTED TRACKS!
1593{
1594 // cout << "InvariantMassTestReco - Start..." << endl;
1595 // cout << "InvariantMassTestReco - Size of fRecoTracklist:\t " << fRecoTracklist.size() << endl;
1596 if (fRecoTracklist.size() >= 4) {
1597 for (unsigned int i = 0; i < fRecoTracklist.size(); i++) {
1598 for (unsigned int j = i + 1; j < fRecoTracklist.size(); j++) {
1599 for (unsigned int k = j + 1; k < fRecoTracklist.size(); k++) {
1600 for (unsigned int l = k + 1; l < fRecoTracklist.size(); l++) {
1601
1602 if (fRecoTracklist[i]->GetPdgCode() + fRecoTracklist[j]->GetPdgCode() + fRecoTracklist[k]->GetPdgCode()
1603 + fRecoTracklist[l]->GetPdgCode()
1604 != 0)
1605 continue;
1606
1607 int motherId1 = fRecoTracklist[i]->GetMotherId();
1608 int motherId2 = fRecoTracklist[j]->GetMotherId();
1609 int motherId3 = fRecoTracklist[k]->GetMotherId();
1610 int motherId4 = fRecoTracklist[l]->GetMotherId();
1611 int grandmotherId1 = -1;
1612 int grandmotherId2 = -1;
1613 int grandmotherId3 = -1;
1614 int grandmotherId4 = -1;
1615
1616 int mcMotherPdg1 = -1;
1617 int mcMotherPdg2 = -1;
1618 int mcMotherPdg3 = -1;
1619 /*
1620 int mcMotherPdg4 = -1;
1621 int mcGrandmotherPdg1 = -1;
1622 int mcGrandmotherPdg2 = -1;
1623 int mcGrandmotherPdg3 = -1;
1624 int mcGrandmotherPdg4 = -1;
1625*/
1626
1627 if (motherId1 != -1) {
1628 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
1629 if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
1630 grandmotherId1 = mother1->GetMotherId();
1631 if (grandmotherId1 != -1) {
1632 // CbmMCTrack* grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
1633 // if (NULL != grandmother1) mcGrandmotherPdg1 = grandmother1->GetPdgCode();
1634 }
1635 }
1636 if (motherId2 != -1) {
1637 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
1638 if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode();
1639 grandmotherId2 = mother2->GetMotherId();
1640 if (grandmotherId2 != -1) {
1641 // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
1642 // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
1643 }
1644 }
1645 if (motherId3 != -1) {
1646 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
1647 if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode();
1648 grandmotherId3 = mother3->GetMotherId();
1649 if (grandmotherId3 != -1) {
1650 // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
1651 // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
1652 }
1653 }
1654 if (motherId4 != -1) {
1655 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
1656 // if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode();
1657 grandmotherId4 = mother4->GetMotherId();
1658 if (grandmotherId4 != -1) {
1659 // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
1660 // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
1661 }
1662 }
1663
1664
1665 if (motherId1 == motherId2 && motherId2 == motherId3 && motherId3 == motherId4) {
1666 Double_t invmass =
1668 fhInvariantMassReco_pi0->Fill(invmass);
1669 // cout << "Decay pi0 -> e+ e- e+ e- detected!" << endl;
1670 }
1671
1672 if (motherId1 == motherId2 && motherId1 == motherId3) {}
1673 if (motherId1 == motherId2 && motherId1 == motherId4) {}
1674 if (motherId1 == motherId3 && motherId1 == motherId4) {}
1675 if (motherId2 == motherId3 && motherId2 == motherId4) {}
1676
1677 // if( ((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 22 || mcMotherPdg1 == 111) && (mcMotherPdg3 == 22 || mcMotherPdg3 == 111))
1678 // || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 22 || mcMotherPdg1 == 111) && (mcMotherPdg2 == 22 || mcMotherPdg2 == 111))
1679 // || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 22 || mcMotherPdg1 == 111) && (mcMotherPdg2 == 22 || mcMotherPdg2 == 111))) {
1680 if (((motherId1 == motherId2 && motherId3 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg3 == 22)
1681 && grandmotherId1 == grandmotherId3)
1682 || ((motherId1 == motherId3 && motherId2 == motherId4) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 22)
1683 && grandmotherId1 == grandmotherId2)
1684 || ((motherId1 == motherId4 && motherId2 == motherId3) && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 22)
1685 && grandmotherId1 == grandmotherId2)) {
1686 Double_t invmass =
1688 fhInvariantMassReco_pi0->Fill(invmass);
1689 // cout << "Decay pi0 -> 2gamma -> e+ e-! \t MotherId " << motherId1 << "\t" << motherId2 << "\t" << motherId3 << "\t" << motherId4 <<
1690 // "\t GrandmotherId " << grandmotherId1 << "\t" << grandmotherId2 << "\t" << grandmotherId3 << "\t" << grandmotherId4 << endl;
1691 }
1692 }
1693 }
1694 }
1695 }
1696 }
1697 // cout << "InvariantMassTestReco - End!" << endl;
1698}
1699
1700
1702
1703
1705
1706
1708// mode 1 = tomography
1709// mode 2 = urqmd
1710// mode 3 = pluto
1711{
1712 fAnalyseMode = mode;
1713}
1714
1715
1717{
1718 if (fMCTracklist_all.size() >= 2) {
1719 for (unsigned int i = 0; i < fMCTracklist_all.size(); i++) {
1720 for (unsigned int j = i + 1; j < fMCTracklist_all.size(); j++) {
1721 if (fMCTracklist_all[i]->GetPdgCode() + fMCTracklist_all[j]->GetPdgCode() != 0) continue;
1722
1723 int motherId1 = fMCTracklist_all[i]->GetMotherId();
1724 int motherId2 = fMCTracklist_all[j]->GetMotherId();
1725
1726 if (motherId1 == motherId2) {
1727 Double_t invmass = Invmass_2particles(fMCTracklist_all[i], fMCTracklist_all[j]);
1728 fhSearchGammas->Fill(invmass);
1729 }
1730 }
1731 }
1732 }
1733}
1734
1735
1737{
1738
1739 fKFparticle = kfparticle;
1740 fKFparticleFinderQA = kfparticleQA;
1741 if (fKFparticle) { cout << "kf works" << endl; }
1742 else {
1743 cout << "kf works not" << endl;
1744 }
1745}
1746
1747
1748Bool_t CbmAnaConversion::AnalysePi0_MC(CbmMCTrack* mctrack, int trackid)
1749{
1750 Bool_t reconstructible = false;
1751 if (mctrack->GetMotherId() == -1) {
1752 fhPi0_MC_occurence->Fill(0);
1753
1754 Int_t daughters = 0;
1755 Int_t gammaDaughters = 0;
1756 Int_t electronDaughters = 0;
1757 Int_t electronDaughtersBeforeRICH = 0;
1758 Int_t electronDaughtersInTarget = 0;
1759 vector<int> gammaDaughterIDs;
1760 gammaDaughterIDs.clear();
1761 vector<int> electronDaughterIDs;
1762 electronDaughterIDs.clear();
1763
1764 Int_t nofMcTracks = fMcTracks->GetEntriesFast();
1765 for (int i = 0; i < nofMcTracks; i++) {
1766 CbmMCTrack* mctrack_switch = (CbmMCTrack*) fMcTracks->At(i);
1767 if (mctrack_switch == NULL) continue;
1768 Int_t motherID = mctrack_switch->GetMotherId();
1769 Int_t pdg = mctrack_switch->GetPdgCode();
1770 if (motherID == trackid && pdg == 22) {
1771 gammaDaughters++;
1772 gammaDaughterIDs.push_back(i);
1773 }
1774 if (motherID == trackid) { daughters++; }
1775 }
1776
1777
1778 if (daughters == 0) fhPi0_MC_occurence->Fill(12);
1779 if (daughters == 1) fhPi0_MC_occurence->Fill(13);
1780 if (daughters == 2) fhPi0_MC_occurence->Fill(14);
1781 if (daughters == 3) fhPi0_MC_occurence->Fill(15);
1782 if (daughters == 4) fhPi0_MC_occurence->Fill(16);
1783 if (daughters > 4) fhPi0_MC_occurence->Fill(17);
1784
1785
1786 if (gammaDaughters == 2) {
1787 fhPi0_MC_occurence->Fill(1);
1788
1789 for (int i = 0; i < nofMcTracks; i++) {
1790 CbmMCTrack* mctrack_switch = (CbmMCTrack*) fMcTracks->At(i);
1791 if (mctrack_switch == NULL) continue;
1792 Int_t motherID = mctrack_switch->GetMotherId();
1793 Int_t pdg = mctrack_switch->GetPdgCode();
1794 if (TMath::Abs(pdg) == 11 && (motherID == gammaDaughterIDs[0] || motherID == gammaDaughterIDs[1])) {
1795 electronDaughters++;
1796 electronDaughterIDs.push_back(i);
1797
1798 TVector3 startvertex;
1799 mctrack_switch->GetStartVertex(startvertex);
1800 if (startvertex.Z() <= 70) { electronDaughtersBeforeRICH++; }
1801 if (startvertex.Z() <= 4) { electronDaughtersInTarget++; }
1802 }
1803 }
1804
1805 if (electronDaughters == 4) { fhPi0_MC_occurence->Fill(2); }
1806 if (electronDaughtersBeforeRICH == 4) {
1807 fhPi0_MC_occurence->Fill(3);
1808 reconstructible = true;
1809 }
1810 if (electronDaughtersInTarget == 4) { fhPi0_MC_occurence->Fill(4); }
1811 }
1812 }
1813 else {
1814 fhPi0_MC_occurence2->Fill(0);
1815 Int_t daughters = 0;
1816 Int_t gammaDaughters = 0;
1817 Int_t electronDaughters = 0;
1818 Int_t electronDaughtersBeforeRICH = 0;
1819 Int_t electronDaughtersInTarget = 0;
1820 vector<int> gammaDaughterIDs;
1821 gammaDaughterIDs.clear();
1822 vector<int> electronDaughterIDs;
1823 electronDaughterIDs.clear();
1824
1825 Int_t nofMcTracks = fMcTracks->GetEntriesFast();
1826 for (int i = 0; i < nofMcTracks; i++) {
1827 CbmMCTrack* mctrack_switch = (CbmMCTrack*) fMcTracks->At(i);
1828 if (mctrack_switch == NULL) continue;
1829 Int_t motherID = mctrack_switch->GetMotherId();
1830 Int_t pdg = mctrack_switch->GetPdgCode();
1831 if (motherID == trackid && pdg == 22) {
1832 gammaDaughters++;
1833 gammaDaughterIDs.push_back(i);
1834 }
1835 if (motherID == trackid) { daughters++; }
1836 }
1837
1838 if (gammaDaughters == 2) {
1839 fhPi0_MC_occurence2->Fill(1);
1840 for (int i = 0; i < nofMcTracks; i++) {
1841 CbmMCTrack* mctrack_switch = (CbmMCTrack*) fMcTracks->At(i);
1842 if (mctrack_switch == NULL) continue;
1843 Int_t motherID = mctrack_switch->GetMotherId();
1844 Int_t pdg = mctrack_switch->GetPdgCode();
1845 if (TMath::Abs(pdg) == 11 && (motherID == gammaDaughterIDs[0] || motherID == gammaDaughterIDs[1])) {
1846 electronDaughters++;
1847 electronDaughterIDs.push_back(i);
1848
1849 TVector3 startvertex;
1850 mctrack_switch->GetStartVertex(startvertex);
1851 if (startvertex.Z() <= 70) { electronDaughtersBeforeRICH++; }
1852 if (startvertex.Z() <= 4) { electronDaughtersInTarget++; }
1853 }
1854 }
1855 if (electronDaughters == 4) { fhPi0_MC_occurence2->Fill(2); }
1856 if (electronDaughtersBeforeRICH == 4) { fhPi0_MC_occurence2->Fill(3); }
1857 if (electronDaughtersInTarget == 4) { fhPi0_MC_occurence2->Fill(4); }
1858 }
1859 }
1860
1861 return reconstructible;
1862}
1863
1864
1866{
1867 Int_t electrons = 0;
1868 std::multimap<int, int> electronPi0ID_map; // contains all IDs of pi0, from which an electron has been detected
1869 std::multimap<int, int> electronPi0ID_map_richInd;
1870
1871 Int_t nof = fTestTracklist.size();
1872 for (int i = 0; i < nof; i++) {
1873 Int_t motherID = fTestTracklist[i]->GetMotherId();
1874 Int_t pdg = fTestTracklist[i]->GetPdgCode();
1875
1876 if (TMath::Abs(pdg) == 11 && motherID != -1) {
1877 CbmMCTrack* mctrack_mother = (CbmMCTrack*) fMcTracks->At(motherID);
1878 if (mctrack_mother == NULL) continue;
1879 Int_t grandmotherID = mctrack_mother->GetMotherId();
1880 Int_t motherPdg = mctrack_mother->GetPdgCode();
1881
1882 if (motherPdg == 22 && grandmotherID != -1) {
1883 CbmMCTrack* mctrack_grandmother = (CbmMCTrack*) fMcTracks->At(grandmotherID);
1884 if (mctrack_grandmother == NULL) continue;
1885 Int_t grandmotherPdg = mctrack_grandmother->GetPdgCode();
1886
1887 if (grandmotherPdg == 111) {
1888 electrons++;
1889 electronPi0ID_map.insert(std::pair<int, int>(grandmotherID, i));
1890 if (!(fTestTracklist_richInd[i] < 0)) {
1891 electronPi0ID_map_richInd.insert(std::pair<int, int>(grandmotherID, i));
1892 }
1895
1896 TVector3 startvertex;
1897 fTestTracklist[i]->GetStartVertex(startvertex);
1898 if (startvertex.Z() < 1) {
1901 }
1902 }
1903 }
1904 }
1905 }
1906
1907 if (electrons >= 4) { // at least 4 electrons from pi0 (NOT from the same one) have been detected
1908 fhPi0_Reco_occurence->Fill(0);
1909 }
1910
1911
1912 int samePi0counter = 0;
1913 int check = 0;
1914 for (std::map<int, int>::iterator it = electronPi0ID_map.begin(); it != electronPi0ID_map.end(); ++it) {
1915 if (it == electronPi0ID_map.begin()) check = 1;
1916 if (it != electronPi0ID_map.begin()) {
1917 std::map<int, int>::iterator zwischen = it;
1918 zwischen--;
1919 int id = it->first;
1920 int id_old = zwischen->first;
1921 if (id == id_old) {
1922 check++;
1923 if (check > 3) {
1924 fhPi0_Reco_occurence->Fill(1);
1925 samePi0counter++;
1926 }
1927 }
1928 else {
1929 if (check == 1) fhPi0_Reco_occurence->Fill(2);
1930 if (check == 2) fhPi0_Reco_occurence->Fill(3);
1931 if (check == 3) fhPi0_Reco_occurence->Fill(4);
1932 if (check == 4) fhPi0_Reco_occurence->Fill(5);
1933 if (check > 4) fhPi0_Reco_occurence->Fill(6);
1934
1935 if (check == 4) {
1936 std::map<int, int>::iterator alt3 = zwischen;
1937 alt3--;
1938 std::map<int, int>::iterator alt4 = alt3;
1939 alt4--;
1940 std::map<int, int>::iterator alt5 = alt4;
1941 alt5--;
1942 cout << "CbmAnaConversion: AnalysePi0_Reco: " << alt5->first << "/" << zwischen->first << "/" << alt3->first
1943 << "/" << alt4->first << endl;
1944 cout << "CbmAnaConversion: AnalysePi0_Reco: " << alt5->second << "/" << zwischen->second << "/"
1945 << alt3->second << "/" << alt4->second << endl;
1946 Bool_t IsWithinCuts = AnalysePi0_Reco_calc(alt5->second, zwischen->second, alt3->second, alt4->second);
1947
1948 Double_t chi_e1 = fTestTracklist_chi[alt5->second];
1949 Double_t chi_e2 = fTestTracklist_chi[zwischen->second];
1950 Double_t chi_e3 = fTestTracklist_chi[alt3->second];
1951 Double_t chi_e4 = fTestTracklist_chi[alt4->second];
1952
1953
1954 Double_t invmass =
1956 fTestTracklist_momentum[alt3->second], fTestTracklist_momentum[alt4->second]);
1957 fhPi0_Reco_invmass->Fill(invmass);
1958 fhPi0_Reco_invmass_cases->Fill(1, invmass);
1959 Double_t invmass_mc = Invmass_4particles(fTestTracklist[alt5->second], fTestTracklist[zwischen->second],
1960 fTestTracklist[alt3->second], fTestTracklist[alt4->second]);
1961 fhPi0_Reco_invmass_mc->Fill(invmass_mc);
1962 fhPi0_Reco_invmass_cases->Fill(5, invmass_mc);
1963
1964
1965 std::map<int, int>::iterator temp = zwischen;
1966 for (int i = 0; i < 4; i++) {
1967 int ndf = fTestTracklist_ndf[temp->second];
1968 int nofhits = fTestTracklist_nofhits[temp->second];
1969 fhPi0_Reco_ndf_vs_nofhits->Fill(ndf, nofhits);
1970 temp--;
1971 }
1972
1973
1974 if (chi_e1 <= 3 && chi_e2 <= 3 && chi_e3 <= 3 && chi_e4 <= 3) {
1975 fhPi0_Reco_occurence->Fill(9);
1976 fhPi0_Reco_invmass_cases->Fill(2, invmass);
1977 if (IsWithinCuts) {
1978 fhPi0_Reco_occurence->Fill(10);
1979 fhPi0_Reco_invmass_cases->Fill(4, invmass);
1980 }
1981
1982 std::map<int, int>::iterator temp2 = zwischen;
1983 for (int i = 0; i < 4; i++) {
1984 int ndf = fTestTracklist_ndf[temp2->second];
1985 int nofhits = fTestTracklist_nofhits[temp2->second];
1986 fhPi0_Reco_ndf_vs_nofhits_withChi->Fill(ndf, nofhits);
1987 temp2--;
1988 }
1989 }
1990 if (IsWithinCuts) fhPi0_Reco_invmass_cases->Fill(3, invmass);
1991
1992 fhPi0_Reco_chi->Fill(chi_e1);
1993 fhPi0_Reco_chi->Fill(chi_e2);
1994 fhPi0_Reco_chi->Fill(chi_e3);
1995 fhPi0_Reco_chi->Fill(chi_e4);
1996 }
1997
1998 check = 1;
1999 }
2000 }
2001 }
2002
2003 if (samePi0counter >= 4) {
2004 //fhPi0_MC_occurence->Fill(13);
2005 }
2006
2007
2008 int samePi0counter2 = 0;
2009 int check2 = 0;
2010 cout << "CbmAnaConversion: RecoPi0: electronmapsize: " << electronPi0ID_map_richInd.size() << endl;
2011 for (std::map<int, int>::iterator it = electronPi0ID_map_richInd.begin(); it != electronPi0ID_map_richInd.end();
2012 ++it) {
2013 if (it == electronPi0ID_map_richInd.begin()) check = 1;
2014 if (it != electronPi0ID_map_richInd.begin()) {
2015 std::map<int, int>::iterator zwischen = it;
2016 zwischen--;
2017 int id = it->first;
2018 int id_old = zwischen->first;
2019 if (id == id_old) {
2020 check2++;
2021 if (check2 > 3) {
2022 //fhPi0_Reco_occurence->Fill(1);
2023 samePi0counter2++;
2024 }
2025 }
2026 else {
2027 if (check2 == 1) fhPi0_Reco_occurence->Fill(11);
2028 if (check2 == 2) fhPi0_Reco_occurence->Fill(12);
2029 if (check2 == 3) fhPi0_Reco_occurence->Fill(13);
2030 if (check2 == 4) fhPi0_Reco_occurence->Fill(14);
2031 if (check2 > 4) fhPi0_Reco_occurence->Fill(15);
2032 /*
2033 if(check == 4) {
2034 std::map<int,int>::iterator alt3 = zwischen;
2035 alt3--;
2036 std::map<int,int>::iterator alt4 = alt3;
2037 alt4--;
2038 std::map<int,int>::iterator alt5 = alt4;
2039 alt5--;
2040 cout << "CbmAnaConversion: AnalysePi0_Reco: " << alt5->first << "/" << zwischen->first << "/" << alt3->first << "/" << alt4->first << endl;
2041 cout << "CbmAnaConversion: AnalysePi0_Reco: " << alt5->second << "/" << zwischen->second << "/" << alt3->second << "/" << alt4->second << endl;
2042 AnalysePi0_Reco_calc(alt5->second, zwischen->second, alt3->second, alt4->second);
2043
2044 }
2045 */
2046 check2 = 1;
2047 }
2048 }
2049 }
2050}
2051
2052
2053Bool_t CbmAnaConversion::AnalysePi0_Reco_calc(int e1, int e2, int e3, int e4)
2054{
2055 CbmMCTrack* mctrack_e1 = fTestTracklist[e1];
2056 CbmMCTrack* mctrack_e2 = fTestTracklist[e2];
2057 CbmMCTrack* mctrack_e3 = fTestTracklist[e3];
2058 CbmMCTrack* mctrack_e4 = fTestTracklist[e4];
2059
2060 Int_t motherID_e1 = mctrack_e1->GetMotherId();
2061 Int_t motherID_e2 = mctrack_e2->GetMotherId();
2062 Int_t motherID_e3 = mctrack_e3->GetMotherId();
2063 Int_t motherID_e4 = mctrack_e4->GetMotherId();
2064
2065
2066 Double_t energy1 = TMath::Sqrt(fTestTracklist_momentum[e1].Mag2() + M2E);
2067 TLorentzVector lorVec1(fTestTracklist_momentum[e1], energy1);
2068
2069 Double_t energy2 = TMath::Sqrt(fTestTracklist_momentum[e2].Mag2() + M2E);
2070 TLorentzVector lorVec2(fTestTracklist_momentum[e2], energy2);
2071
2072 Double_t energy3 = TMath::Sqrt(fTestTracklist_momentum[e3].Mag2() + M2E);
2073 TLorentzVector lorVec3(fTestTracklist_momentum[e3], energy3);
2074
2075 Double_t energy4 = TMath::Sqrt(fTestTracklist_momentum[e4].Mag2() + M2E);
2076 TLorentzVector lorVec4(fTestTracklist_momentum[e4], energy4);
2077
2078
2079 Bool_t IsWithinCuts = false;
2080 Double_t OpeningAngleCut = 1.5;
2081
2082 if (motherID_e1 == motherID_e2 && motherID_e3 == motherID_e4) {
2083 Double_t anglePair1 = lorVec1.Angle(lorVec2.Vect());
2084 Double_t theta1 = 180. * anglePair1 / TMath::Pi();
2085
2086 Double_t anglePair2 = lorVec3.Angle(lorVec4.Vect());
2087 Double_t theta2 = 180. * anglePair2 / TMath::Pi();
2088
2089 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut) IsWithinCuts = true;
2090 fhPi0_Reco_angles->Fill(theta1);
2091 fhPi0_Reco_angles->Fill(theta2);
2092 cout << "CbmAnaConversion: AnalysePi0_Reco_calc: " << theta1 << "/" << theta2 << endl;
2093 }
2094
2095 if (motherID_e1 == motherID_e3 && motherID_e2 == motherID_e4) {
2096 Double_t anglePair1 = lorVec1.Angle(lorVec3.Vect());
2097 Double_t theta1 = 180. * anglePair1 / TMath::Pi();
2098
2099 Double_t anglePair2 = lorVec2.Angle(lorVec4.Vect());
2100 Double_t theta2 = 180. * anglePair2 / TMath::Pi();
2101
2102 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut) IsWithinCuts = true;
2103 fhPi0_Reco_angles->Fill(theta1);
2104 fhPi0_Reco_angles->Fill(theta2);
2105 cout << "CbmAnaConversion: AnalysePi0_Reco_calc: " << theta1 << "/" << theta2 << endl;
2106 }
2107
2108 if (motherID_e1 == motherID_e4 && motherID_e2 == motherID_e3) {
2109 Double_t anglePair1 = lorVec1.Angle(lorVec4.Vect());
2110 Double_t theta1 = 180. * anglePair1 / TMath::Pi();
2111
2112 Double_t anglePair2 = lorVec2.Angle(lorVec3.Vect());
2113 Double_t theta2 = 180. * anglePair2 / TMath::Pi();
2114
2115 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut) IsWithinCuts = true;
2116 fhPi0_Reco_angles->Fill(theta1);
2117 fhPi0_Reco_angles->Fill(theta2);
2118 cout << "CbmAnaConversion: AnalysePi0_Reco_calc: " << theta1 << "/" << theta2 << endl;
2119 }
2120
2121 if (IsWithinCuts == true) { fhPi0_Reco_occurence->Fill(8); }
2122
2123 return IsWithinCuts;
2124}
2125
2126
2128{
2129 Int_t electrons = 0;
2130 std::multimap<int, int> electronPi0ID_map; // contains all IDs of pi0, from which an electron has been detected
2131 std::multimap<int, int>
2132 electronPi0ID_map_richInd; // contains all IDs of pi0, from which an electron has been detected
2133
2134 Int_t nof = fTestTracklist_noRichInd.size();
2135 for (int i = 0; i < nof; i++) {
2136 Int_t motherID = fTestTracklist_noRichInd[i]->GetMotherId();
2137 Int_t pdg = fTestTracklist_noRichInd[i]->GetPdgCode();
2138
2139 if (TMath::Abs(pdg) == 11 && motherID != -1) {
2140 CbmMCTrack* mctrack_mother = (CbmMCTrack*) fMcTracks->At(motherID);
2141 if (mctrack_mother == NULL) continue;
2142 Int_t grandmotherID = mctrack_mother->GetMotherId();
2143 Int_t motherPdg = mctrack_mother->GetPdgCode();
2144
2145 if (motherPdg == 22 && grandmotherID != -1) {
2146 CbmMCTrack* mctrack_grandmother = (CbmMCTrack*) fMcTracks->At(grandmotherID);
2147 if (mctrack_grandmother == NULL) continue;
2148 Int_t grandmotherPdg = mctrack_grandmother->GetPdgCode();
2149
2150 if (grandmotherPdg == 111) {
2151 electrons++;
2152 electronPi0ID_map.insert(std::pair<int, int>(grandmotherID, i));
2153 if (!(fTestTracklist_noRichInd_richInd[i] < 0)) {
2154 electronPi0ID_map_richInd.insert(std::pair<int, int>(grandmotherID, i));
2155 }
2160
2161 TVector3 startvertex;
2162 fTestTracklist_noRichInd[i]->GetStartVertex(startvertex);
2163 if (startvertex.Z() < 1) {
2168 }
2169 }
2170 else {
2175 }
2176 }
2177 else {
2182 }
2183 }
2184 else {
2185 if (TMath::Abs(pdg) != 11) continue;
2190 }
2191 }
2192
2193 if (electrons >= 4) { // at least 4 electrons from pi0 (NOT from the same one) have been detected
2194 fhPi0_Reco_occurence2->Fill(0);
2195 }
2196
2197
2198 int samePi0counter = 0;
2199 int check = 0;
2200 int check_withRichInd = 0;
2201 for (std::map<int, int>::iterator it = electronPi0ID_map.begin(); it != electronPi0ID_map.end(); ++it) {
2202 if (it == electronPi0ID_map.begin()) check = 1;
2203 if (it != electronPi0ID_map.begin()) {
2204 std::map<int, int>::iterator zwischen = it;
2205 zwischen--;
2206 int id = it->first;
2207 int id_old = zwischen->first;
2208 if (id == id_old) {
2209 check++;
2210 if (check > 3) {
2211 //fhPi0_Reco_occurence->Fill(1);
2212 samePi0counter++;
2213 }
2214 if (fTestTracklist_noRichInd_richInd[it->second] >= 0) { check_withRichInd++; }
2215 }
2216 else {
2217 if (check == 1) fhPi0_Reco_occurence2->Fill(2);
2218 if (check == 2) fhPi0_Reco_occurence2->Fill(3);
2219 if (check == 3) fhPi0_Reco_occurence2->Fill(4);
2220 if (check == 4) fhPi0_Reco_occurence2->Fill(5);
2221 if (check > 4) fhPi0_Reco_occurence2->Fill(6);
2222
2223 if (check == 4) {
2224 std::map<int, int>::iterator alt3 = zwischen;
2225 alt3--;
2226 std::map<int, int>::iterator alt4 = alt3;
2227 alt4--;
2228 std::map<int, int>::iterator alt5 = alt4;
2229 alt5--;
2230 cout << "CbmAnaConversion: AnalysePi0_Reco_noRichInd: " << alt5->first << "/" << zwischen->first << "/"
2231 << alt3->first << "/" << alt4->first << endl;
2232 cout << "CbmAnaConversion: AnalysePi0_Reco_noRichInd: " << alt5->second << "/" << zwischen->second << "/"
2233 << alt3->second << "/" << alt4->second << endl;
2234 Bool_t IsWithinCuts =
2235 AnalysePi0_Reco_noRichInd_calc(alt5->second, zwischen->second, alt3->second, alt4->second);
2236
2237 Double_t chi_e1 = fTestTracklist_noRichInd_chi[alt5->second];
2238 Double_t chi_e2 = fTestTracklist_noRichInd_chi[zwischen->second];
2239 Double_t chi_e3 = fTestTracklist_noRichInd_chi[alt3->second];
2240 Double_t chi_e4 = fTestTracklist_noRichInd_chi[alt4->second];
2241
2242 Int_t richInd_e1 = fTestTracklist_noRichInd_richInd[alt5->second];
2243 Int_t richInd_e2 = fTestTracklist_noRichInd_richInd[zwischen->second];
2244 Int_t richInd_e3 = fTestTracklist_noRichInd_richInd[alt3->second];
2245 Int_t richInd_e4 = fTestTracklist_noRichInd_richInd[alt4->second];
2246
2247
2248 Double_t invmass = Invmass_4particlesRECO(
2251 fhPi0_Reco_noRichInd_invmass->Fill(invmass);
2252 fhPi0_Reco_noRichInd_invmass_cases->Fill(1, invmass);
2253 Double_t invmass_mc =
2255 fTestTracklist_noRichInd[alt3->second], fTestTracklist_noRichInd[alt4->second]);
2256 fhPi0_Reco_noRichInd_invmass_mc->Fill(invmass_mc);
2257 fhPi0_Reco_noRichInd_invmass_cases->Fill(5, invmass_mc);
2258
2259
2260 std::map<int, int>::iterator temp = zwischen;
2261 for (int i = 0; i < 4; i++) {
2262 int ndf = fTestTracklist_noRichInd_ndf[temp->second];
2263 int nofhits = fTestTracklist_noRichInd_nofhits[temp->second];
2264 fhPi0_Reco_noRichInd_ndf_vs_nofhits->Fill(ndf, nofhits);
2265 temp--;
2266 }
2267
2268 if (chi_e1 <= 3 && chi_e2 <= 3 && chi_e3 <= 3 && chi_e4 <= 3) {
2269 fhPi0_Reco_occurence2->Fill(9);
2270 fhPi0_Reco_noRichInd_invmass_cases->Fill(2, invmass);
2271 if (IsWithinCuts) {
2272 fhPi0_Reco_occurence2->Fill(10);
2273 fhPi0_Reco_noRichInd_invmass_cases->Fill(4, invmass);
2274 }
2275
2276 std::map<int, int>::iterator temp2 = zwischen;
2277 for (int i = 0; i < 4; i++) {
2278 int ndf = fTestTracklist_noRichInd_ndf[temp2->second];
2279 int nofhits = fTestTracklist_noRichInd_nofhits[temp2->second];
2281 temp2--;
2282 }
2283 }
2284
2285 if (IsWithinCuts) { fhPi0_Reco_noRichInd_invmass_cases->Fill(3, invmass); }
2286
2287 if (richInd_e1 >= 0 && richInd_e2 >= 0 && richInd_e3 >= 0 && richInd_e4 >= 0) {
2288 fhPi0_Reco_noRichInd_invmass_cases->Fill(6, invmass);
2289 fhPi0_Reco_noRichInd_invmass_cases->Fill(10, invmass_mc);
2290 if (IsWithinCuts) { fhPi0_Reco_noRichInd_invmass_cases->Fill(8, invmass); }
2291 if (chi_e1 <= 3 && chi_e2 <= 3 && chi_e3 <= 3 && chi_e4 <= 3) {
2292 fhPi0_Reco_noRichInd_invmass_cases->Fill(7, invmass);
2293 if (IsWithinCuts) { fhPi0_Reco_noRichInd_invmass_cases->Fill(9, invmass); }
2294 }
2295 }
2296 }
2297
2298 if (check > 4) {
2299 std::map<int, int>::iterator temp = zwischen;
2300 cout << "CbmAnaConversion: AnalysePi0_Reco_noRichInd, check>4: ";
2301 for (int i = 0; i < check; i++) {
2302 TVector3 momentum_mc;
2303 fTestTracklist_noRichInd[temp->second]->GetMomentum(momentum_mc);
2304 Double_t theta_mc = 180.0 * momentum_mc.Theta() / TMath::Pi();
2305 Double_t phi_mc = 180.0 * momentum_mc.Phi() / TMath::Pi();
2306 TVector3 momentum_reco = fTestTracklist_noRichInd_momentum[temp->second];
2307 Double_t theta_reco = 180.0 * momentum_reco.Theta() / TMath::Pi();
2308 Double_t phi_reco = 180.0 * momentum_reco.Phi() / TMath::Pi();
2309 cout << "(" << temp->first << "/" << temp->second << "/" << fTestTracklist_noRichInd_MCindex[temp->second]
2310 << "/" << momentum_reco.Mag() << "/" << theta_mc << "-" << phi_mc << "/" << theta_reco << "-"
2311 << phi_reco << ") \t";
2312 temp--;
2313 fhPi0_noRichInd_diffPhi->Fill(TMath::Abs(phi_mc - phi_reco));
2314 fhPi0_noRichInd_diffTheta->Fill(TMath::Abs(theta_mc - theta_reco));
2315 }
2316 cout << endl;
2317 }
2318
2319
2320 if (check_withRichInd == 1) fhPi0_Reco_occurence2->Fill(11);
2321 if (check_withRichInd == 2) fhPi0_Reco_occurence2->Fill(12);
2322 if (check_withRichInd == 3) fhPi0_Reco_occurence2->Fill(13);
2323 if (check_withRichInd == 4) fhPi0_Reco_occurence2->Fill(14);
2324 if (check_withRichInd > 4) fhPi0_Reco_occurence2->Fill(15);
2325
2326
2327 // reset of check-parameters for next loop
2328 check = 1;
2329 if (fTestTracklist_noRichInd_richInd[it->second] >= 0) { check_withRichInd = 1; }
2330 else {
2331 check_withRichInd = 0;
2332 }
2333 }
2334 }
2335 }
2336
2337 if (samePi0counter >= 4) {
2338 //fhPi0_MC_occurence->Fill(13);
2339 }
2340}
2341
2342
2343Bool_t CbmAnaConversion::AnalysePi0_Reco_noRichInd_calc(int e1, int e2, int e3, int e4)
2344{
2345 CbmMCTrack* mctrack_e1 = fTestTracklist_noRichInd[e1];
2346 CbmMCTrack* mctrack_e2 = fTestTracklist_noRichInd[e2];
2347 CbmMCTrack* mctrack_e3 = fTestTracklist_noRichInd[e3];
2348 CbmMCTrack* mctrack_e4 = fTestTracklist_noRichInd[e4];
2349
2350 Int_t motherID_e1 = mctrack_e1->GetMotherId();
2351 Int_t motherID_e2 = mctrack_e2->GetMotherId();
2352 Int_t motherID_e3 = mctrack_e3->GetMotherId();
2353 Int_t motherID_e4 = mctrack_e4->GetMotherId();
2354
2355
2356 Double_t energy1 = TMath::Sqrt(fTestTracklist_noRichInd_momentum[e1].Mag2() + M2E);
2357 TLorentzVector lorVec1(fTestTracklist_noRichInd_momentum[e1], energy1);
2358
2359 Double_t energy2 = TMath::Sqrt(fTestTracklist_noRichInd_momentum[e2].Mag2() + M2E);
2360 TLorentzVector lorVec2(fTestTracklist_noRichInd_momentum[e2], energy2);
2361
2362 Double_t energy3 = TMath::Sqrt(fTestTracklist_noRichInd_momentum[e3].Mag2() + M2E);
2363 TLorentzVector lorVec3(fTestTracklist_noRichInd_momentum[e3], energy3);
2364
2365 Double_t energy4 = TMath::Sqrt(fTestTracklist_noRichInd_momentum[e4].Mag2() + M2E);
2366 TLorentzVector lorVec4(fTestTracklist_noRichInd_momentum[e4], energy4);
2367
2368
2369 Bool_t IsWithinCuts = false;
2370 Double_t OpeningAngleCut = 1.5;
2371
2372 if (motherID_e1 == motherID_e2 && motherID_e3 == motherID_e4) {
2373 Double_t anglePair1 = lorVec1.Angle(lorVec2.Vect());
2374 Double_t theta1 = 180. * anglePair1 / TMath::Pi();
2375
2376 Double_t anglePair2 = lorVec3.Angle(lorVec4.Vect());
2377 Double_t theta2 = 180. * anglePair2 / TMath::Pi();
2378
2379 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut) IsWithinCuts = true;
2380 cout << "CbmAnaConversion: AnalysePi0_Reco_noRichInd_calc: " << theta1 << "/" << theta2 << endl;
2381 }
2382
2383 if (motherID_e1 == motherID_e3 && motherID_e2 == motherID_e4) {
2384 Double_t anglePair1 = lorVec1.Angle(lorVec3.Vect());
2385 Double_t theta1 = 180. * anglePair1 / TMath::Pi();
2386
2387 Double_t anglePair2 = lorVec2.Angle(lorVec4.Vect());
2388 Double_t theta2 = 180. * anglePair2 / TMath::Pi();
2389
2390 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut) IsWithinCuts = true;
2391 cout << "CbmAnaConversion: AnalysePi0_Reco_noRichInd_calc: " << theta1 << "/" << theta2 << endl;
2392 }
2393
2394 if (motherID_e1 == motherID_e4 && motherID_e2 == motherID_e3) {
2395 Double_t anglePair1 = lorVec1.Angle(lorVec4.Vect());
2396 Double_t theta1 = 180. * anglePair1 / TMath::Pi();
2397
2398 Double_t anglePair2 = lorVec2.Angle(lorVec3.Vect());
2399 Double_t theta2 = 180. * anglePair2 / TMath::Pi();
2400
2401 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut) IsWithinCuts = true;
2402 cout << "CbmAnaConversion: AnalysePi0_Reco_noRichInd_calc: " << theta1 << "/" << theta2 << endl;
2403 }
2404
2405 if (IsWithinCuts == true) { fhPi0_Reco_occurence2->Fill(8); }
2406
2407 return IsWithinCuts;
2408}
2409
2410
#define M2E
ClassImp(CbmConverterManager)
Helper functions for drawing 1D and 2D histograms and graphs.
Data class for STS tracks.
fscal v[fmask::Size]
Definition KfSimdPseudo.h:4
static Double_t CalcChiCut(Double_t pt)
void SetKF(CbmKFParticleFinder *kfparticle, CbmKFParticleFinderQa *kfparticleQA)
void SetTracklistMC(std::vector< CbmMCTrack * > MCTracklist)
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)
void TomographyReco(CbmMCTrack *mctrack)
Bool_t FillRecoTracklistEPEM(CbmMCTrack *mctrack, TVector3 stsMomentum, TVector3 refittedMom, int i, Double_t chi, Int_t GlobalTrackId)
std::vector< CbmMCTrack * > fTestTracklist
std::vector< TVector3 > fTestTracklist_momentum
TH2D * fhPi0_Reco_noRichInd_chi_vs_pt
TH2D * fhPi0_Reco_ndf_vs_startvertex
TH2D * fhPi0_Reco_startvertex_vs_nofhits
TH2D * fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0
std::vector< TH1 * > fHistoList_kfparticle
std::vector< CbmMCTrack * > fMCTracklist_all
TH2D * fhPi0_Reco_noRichInd_chi_vs_momentum
TH2D * fhPi0_Reco_chi_vs_momentum
std::vector< int > fTestTracklist_nofhits
std::vector< CbmMCTrack * > fRecoTracklist
std::vector< Double_t > fRecoTracklistEPEM_chi
std::vector< CbmMCTrack * > fRecoTracklistEPEM
TClonesArray * fGlobalTracks
TH2D * fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target
std::vector< TVector3 > fTestTracklist_noRichInd_momentum
TClonesArray * fRichRings
std::vector< TH1 * > fHistoList_furtherAnalyses
Double_t Invmass_4particles(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
TH2D * fhPi0_Reco_chi_vs_pt_eFromPi0
CbmAnaConversion()
Standard constructor.
std::vector< int > fTestTracklist_noRichInd_ndf
void AnalyseElectrons(CbmMCTrack *mctrack)
std::vector< float > fTestTracklist_noRichInd_chi
TH2D * fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0
CbmKFParticleFinder * fKFparticle
TH2D * fhPi0_Reco_chi_vs_momentum_eFromPi0
CbmAnaConversionPhotons2 * fAnaPhotons2
std::vector< int > fTestTracklist_richInd
Bool_t AnalysePi0_Reco_noRichInd_calc(int e1, int e2, int e3, int e4)
TH2D * fhPi0_Reco_noRichInd_chi_vs_pt_eRest
CbmKFParticleFinderQa * fKFparticleFinderQA
Double_t Invmass_2particles(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2)
std::vector< int > fTestTracklist_ndf
virtual void Finish()
Inherited from FairTask.
std::vector< TVector3 > fRecoMomentum
CbmAnaConversionReco * fAnaReco
std::vector< CbmMCTrack * > fMCTracklist
virtual InitStatus Init()
Inherited from FairTask.
TH2D * fhPi0_Reco_ndf_vs_nofhits_withChi
std::vector< TH1 * > fHistoList
Double_t Invmass_4particlesRECO(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
LmvmKinePar CalculateKinematicParams(const CbmMCTrack *mctrackP, const CbmMCTrack *mctrackM)
Bool_t AnalysePi0_MC(CbmMCTrack *mctrack, int i)
TClonesArray * fStsTrackMatches
std::vector< TVector3 > fRecoRefittedMomentum
void FillMCTracklists(CbmMCTrack *mctrack, int i)
std::vector< Int_t > fRecoTracklistEPEM_gtid
TH2D * fhANN_output_electrons_vs_p
CbmAnaConversionKF * fAnaKF
void FillRecoTracklist(CbmMCTrack *mtrack)
TH2D * fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target
TH2D * fhPi0_Reco_chi_vs_pt_eFromPi0_Target
void InitHistograms()
Initialize histograms.
virtual void Exec(Option_t *option)
Inherited from FairTask.
void SetKF(CbmKFParticleFinder *kfparticle, CbmKFParticleFinderQa *kfparticleQA)
CbmAnaConversionRecoFull * fAnaRecoFull
TClonesArray * fRichRingMatches
std::vector< int > fTestTracklist_noRichInd_gTrackId
std::vector< int > fRecoTracklistEPEM_id
TH1D * fhANN_output_electrons_chiCut
TH2D * fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi
Bool_t AnalysePi0_Reco_calc(int e1, int e2, int e3, int e4)
CbmAnaConversionTest2 * fAnaTest2
TH1D * fhPi0_Reco_noRichInd_invmass
TH2D * fhPi0_Reco_noRichInd_chi_vs_momentum_eRest
CbmAnaConversionTomography * fAnaTomography
CbmAnaConversionRich * fAnaRich
TClonesArray * fMcTracks
TClonesArray * fRichPoints
TH1D * fhPi0_Reco_noRichInd_invmass_mc
TH2D * fhPi0_Reco_chi_vs_momentum_eFromPi0_Target
CbmVertex * fPrimVertex
virtual ~CbmAnaConversion()
Standard destructor.
Double_t Invmass_2gammas(const CbmMCTrack *gamma1, const CbmMCTrack *gamma2)
TH2D * fhPi0_Reco_noRichInd_invmass_cases
CbmAnaConversionPhotons * fAnaPhotons
TClonesArray * fStsTracks
CbmAnaConversionTest * fAnaTest
std::vector< CbmMCTrack * > fTestTracklist_noRichInd
TH2D * fhPi0_Reco_startvertex_vs_chi
std::vector< int > fTestTracklist_noRichInd_MCindex
std::vector< float > fTestTracklist_chi
std::vector< int > fTestTracklist_noRichInd_richInd
TH2D * fhPi0_Reco_noRichInd_ndf_vs_nofhits
std::vector< int > fTestTracklist_noRichInd_nofhits
int32_t GetStsTrackIndex() const
int32_t GetRichRingIndex() const
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
void Fit(std::vector< CbmStsTrack > &Tracks, const std::vector< CbmMvdHit > &vMvdHits, const std::vector< CbmStsHit > &vStsHits, const std::vector< int > &pidHypo)
double GetPt() const
Definition CbmMCTrack.h:97
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:170
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
void Get4Momentum(TLorentzVector &momentum) const
Definition CbmMCTrack.h:173
double GetRapidity() const
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:176
int32_t GetNofLinks() const
Definition CbmMatch.h:42
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
double CalculateAnnValue(int globalTrackIndex, double momentum)
Calculate output value of the ANN.
static CbmRichElectronIdAnn & GetInstance()
virtual int32_t GetTotalNofHits() const
Definition CbmStsTrack.h:81
double GetTrueOverAllHitsRatio() const
const FairTrackParam * GetParamFirst() const
Definition CbmTrack.h:68
Double_t fAngle
Definition LmvmKinePar.h:23
Double_t fMinv
Definition LmvmKinePar.h:22
Double_t fPt
Definition LmvmKinePar.h:20
Double_t fMomentumMag
Definition LmvmKinePar.h:19
Double_t fRapidity
Definition LmvmKinePar.h:21
Hash for CbmL1LinkKey.