CbmRoot
Loading...
Searching...
No Matches
CbmAnaConversionTest.cxx
Go to the documentation of this file.
1/* Copyright (C) 2015-2019 Fakultaet fuer Mathematik und Naturwissenschaften, Bergische Universitaet Wuppertal, Wuppertal
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Sascha Reinecke [committer], Florian Uhlig, Andrey Lebedev */
4
12#include "CbmAnaConversionTest.h"
13
14#include "CbmAnaConversionCutSettings.h"
15#include "CbmAnaConversionGlobalFunctions.h"
16#include "CbmAnaConversionKinematicParams.h"
17#include "CbmGlobalTrack.h"
18#include "CbmL1PFFitter.h"
19#include "CbmLitGlobalElectronId.h"
20#include "CbmMCTrack.h"
21#include "CbmRichPoint.h"
22#include "CbmRichRing.h"
23#include "CbmStsTrack.h"
24#include "CbmTrackMatchNew.h"
25
26#include "FairRootManager.h"
27#include <Logger.h>
28
29#include <algorithm>
30#include <map>
31
32
33using namespace std;
34
35
36CbmAnaConversionTest::CbmAnaConversionTest()
37 : fRichPoints(NULL)
38 , fRichRings(NULL)
39 , fRichRingMatches(NULL)
40 , fMcTracks(NULL)
41 , fStsTracks(NULL)
42 , fStsTrackMatches(NULL)
43 , fGlobalTracks(NULL)
44 , fPrimVertex(NULL)
45 , fKFVertex()
46 , fHistoList_test()
47 , fElectrons_gtid()
48 , fElectrons_mcid()
49 , fElectrons_richInd()
50 , fElectrons_pi0mcid()
51 , fElectrons_same()
52 , fElectrons_nofPerPi0(NULL)
53 , fElectrons_nofPerPi0_withRichInd(NULL)
54 , fhElectronsTest_invmass(NULL)
55 , fhTest_ParticlesPerEvent(NULL)
56 , fhTest_RICHelectronsPerEvent(NULL)
57 , fhTest_PhotonsPerEvent_RICHonly(NULL)
58 , fhTest_PhotonsPerEvent_STSandRICH(NULL)
59 , fhTest_ReconstructedPi0PerEvent(NULL)
60 , fhTest_invmass(NULL)
61 , fhTest_invmass_pCut(NULL)
62 , fhTest_invmass_GGcut(NULL)
63 , fhTest_invmass_RICHindex0(NULL)
64 , fhTest_invmass_RICHindex1(NULL)
65 , fhTest_invmass_RICHindex2(NULL)
66 , fhTest_invmass_RICHindex3(NULL)
67 , fhTest_invmass_RICHindex4(NULL)
68 , fhTest_invmass_MCcutAll(NULL)
69 , fhTest_peakCheck1(NULL)
70 , fhTest_peakCheck2(NULL)
71 , fhTest_peakCheck3(NULL)
72 , fhTest_invmass_ANNcuts(NULL)
73 , fhTest_phaseSpace_pi0(NULL)
74 , fhTest_phaseSpace_eta(NULL)
75 , fVector_AllMomenta()
76 , fVector_gt()
77 , fVector_momenta()
78 , fVector_chi()
79 , fVector_gtIndex()
80 , fVector_richIndex()
81 , fVector_mcIndex()
82 , fVector_withRichSignal()
83 , fVector_reconstructedPhotons_FromSTSandRICH()
84 , fVector_electronRICH_gt()
85 , fVector_electronRICH_gtIndex()
86 , fVector_electronRICH_mcIndex()
87 , fVector_electronRICH_momenta()
88 , fVector_electronRICH_reconstructedPhotons()
89 , fVector_reconstructedPhotons_STSonly()
90 , globalEventNo(0)
91 , fMixedTest_STSonly_photons()
92 , fMixedTest_STSonly_eventno()
93 , fMixedTest_STSonly_hasRichInd()
94 , fhTest_eventMixing_STSonly_2p2(NULL)
95 , fhTest_eventMixing_STSonly_3p1(NULL)
96 , fhTest_eventMixing_STSonly_4p0(NULL)
97 , fMixedTest_3p1_photons()
98 , fMixedTest_3p1_eventno()
99 , fMixedTest_3p1_combined()
100 , fMixedTest_3p1_ann()
101 , fhTest_eventMixing_3p1(NULL)
102 , fhTest_eventMixing_3p1_pCut(NULL)
103 , fhTest_eventMixing_3p1_GGcut(NULL)
104 , fhTest_eventMixing_3p1_ANNcuts(NULL)
105{
106}
107
108CbmAnaConversionTest::~CbmAnaConversionTest() {}
109
110
111void CbmAnaConversionTest::Init()
112{
113 FairRootManager* ioman = FairRootManager::Instance();
114 if (NULL == ioman) { Fatal("CbmAnaConversion::Init", "RootManager not instantised!"); }
115
116 fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
117 if (NULL == fRichPoints) { Fatal("CbmAnaConversion::Init", "No RichPoint array!"); }
118
119 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
120 if (NULL == fMcTracks) { Fatal("CbmAnaConversion::Init", "No MCTrack array!"); }
121
122 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
123 if (NULL == fStsTracks) { Fatal("CbmAnaConversion::Init", "No StsTrack array!"); }
124
125 fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
126 if (NULL == fStsTrackMatches) { Fatal("CbmAnaConversion::Init", "No StsTrackMatch array!"); }
127
128 fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
129 if (NULL == fGlobalTracks) { Fatal("CbmAnaConversion::Init", "No GlobalTrack array!"); }
130
131 // Get pointer to PrimaryVertex object from IOManager if it exists
132 // The old name for the object is "PrimaryVertex" the new one
133 // "PrimaryVertex." Check first for the new name
134 fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex."));
135 if (nullptr == fPrimVertex) { fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex")); }
136 if (nullptr == fPrimVertex) { LOG(fatal) << "No PrimaryVertex array!"; }
137
138 fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
139 if (NULL == fRichRings) { Fatal("CbmAnaConversion::Init", "No RichRing array!"); }
140
141 fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
142 if (NULL == fRichRingMatches) { Fatal("CbmAnaConversion::Init", "No RichRingMatch array!"); }
143
144 InitHistos();
145 CbmLitGlobalElectronId::GetInstance().SetRichAnnCut(-0.8);
146
147
148 globalEventNo = 0;
149}
150
151
152void CbmAnaConversionTest::InitHistos()
153{
154 fHistoList_test.clear();
155
156 Double_t invmassSpectra_nof = 800;
157 Double_t invmassSpectra_start = -0.00125;
158 Double_t invmassSpectra_end = 1.99875;
159
160
161 fElectrons_nofPerPi0 = new TH1I("fElectrons_nofPerPi0", "fElectrons_nofPerPi0; nof; #", 7, -0.5, 6.5);
162 fHistoList_test.push_back(fElectrons_nofPerPi0);
163 fElectrons_nofPerPi0_withRichInd =
164 new TH1I("fElectrons_nofPerPi0_withRichInd", "fElectrons_nofPerPi0_withRichInd; nof; #", 7, -0.5, 6.5);
165 fHistoList_test.push_back(fElectrons_nofPerPi0_withRichInd);
166
167 fhElectronsTest_invmass =
168 new TH1D("fhElectronsTest_invmass", "fhElectronsTest_invmass; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
169 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
170 fHistoList_test.push_back(fhElectronsTest_invmass);
171
172
173 fhTest_PhotonsPerEvent_RICHonly =
174 new TH1I("fhTest_PhotonsPerEvent_RICHonly", "fhTest_PhotonsPerEvent_RICHonly; nof; #", 501, -0.5, 500.5);
175 fHistoList_test.push_back(fhTest_PhotonsPerEvent_RICHonly);
176 fhTest_PhotonsPerEvent_STSandRICH =
177 new TH1I("fhTest_PhotonsPerEvent_STSandRICH", "fhTest_PhotonsPerEvent_STSandRICH; nof; #", 501, -0.5, 500.5);
178 fHistoList_test.push_back(fhTest_PhotonsPerEvent_STSandRICH);
179 fhTest_ParticlesPerEvent =
180 new TH1I("fhTest_ParticlesPerEvent", "fhTest_ParticlesPerEvent; nof; #", 1001, -0.5, 1000.5);
181 fHistoList_test.push_back(fhTest_ParticlesPerEvent);
182 fhTest_ReconstructedPi0PerEvent =
183 new TH1I("fhTest_ReconstructedPi0PerEvent", "fhTest_ReconstructedPi0PerEvent; nof; #", 101, -0.5, 100.5);
184 fHistoList_test.push_back(fhTest_ReconstructedPi0PerEvent);
185 fhTest_RICHelectronsPerEvent =
186 new TH1I("fhTest_RICHelectronsPerEvent", "fhTest_RICHelectronsPerEvent; nof; #", 501, -0.5, 500.5);
187 fHistoList_test.push_back(fhTest_RICHelectronsPerEvent);
188
189 fhTest_invmass =
190 new TH1D("fhTest_invmass", "fhTest_invmass; invariant mass of 3 e^{#pm} + 1 particle in GeV/c^{2}; #",
191 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
192 fHistoList_test.push_back(fhTest_invmass);
193 fhTest_invmass_pCut = new TH1D("fhTest_invmass_pCut",
194 "fhTest_invmass_pCut; invariant mass of 3 "
195 "e^{#pm} + 1 particle in GeV/c^{2}; #",
196 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
197 fHistoList_test.push_back(fhTest_invmass_pCut);
198 fhTest_invmass_GGcut = new TH1D("fhTest_invmass_GGcut",
199 "fhTest_invmass_GGcut; invariant mass of 3 "
200 "e^{#pm} + 1 particle in GeV/c^{2}; #",
201 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
202 fHistoList_test.push_back(fhTest_invmass_GGcut);
203
204
205 fhTest_invmass_RICHindex0 =
206 new TH1D("fhTest_invmass_RICHindex0", "fhTest_invmass_RICHindex0; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
207 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
208 fHistoList_test.push_back(fhTest_invmass_RICHindex0);
209 fhTest_invmass_RICHindex1 =
210 new TH1D("fhTest_invmass_RICHindex1", "fhTest_invmass_RICHindex1; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
211 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
212 fHistoList_test.push_back(fhTest_invmass_RICHindex1);
213 fhTest_invmass_RICHindex2 =
214 new TH1D("fhTest_invmass_RICHindex2", "fhTest_invmass_RICHindex2; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
215 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
216 fHistoList_test.push_back(fhTest_invmass_RICHindex2);
217 fhTest_invmass_RICHindex3 =
218 new TH1D("fhTest_invmass_RICHindex3", "fhTest_invmass_RICHindex3; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
219 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
220 fHistoList_test.push_back(fhTest_invmass_RICHindex3);
221 fhTest_invmass_RICHindex4 =
222 new TH1D("fhTest_invmass_RICHindex4", "fhTest_invmass_RICHindex4; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
223 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
224 fHistoList_test.push_back(fhTest_invmass_RICHindex4);
225
226 fhTest_invmass_MCcutAll = new TH2D("fhTest_invmass_MCcutAll", "fhTest_invmass_MCcutAll; case; invariant mass", 25, 0,
227 25, invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
228 fhTest_peakCheck1 = new TH1D("fhTest_peakCheck1", "fhTest_peakCheck1; sum; #", 15, -0.5, 14.5);
229 fhTest_peakCheck2 = new TH1D("fhTest_peakCheck2", "fhTest_peakCheck2; sum; #", 15, -0.5, 14.5);
230 fhTest_peakCheck3 = new TH1D("fhTest_peakCheck3", "fhTest_peakCheck3; sum; #", 20, -5.5, 14.5);
231 fHistoList_test.push_back(fhTest_invmass_MCcutAll);
232 fHistoList_test.push_back(fhTest_peakCheck1);
233 fHistoList_test.push_back(fhTest_peakCheck2);
234 fHistoList_test.push_back(fhTest_peakCheck3);
235
236 fhTest_invmass_ANNcuts = new TH2D("fhTest_invmass_ANNcuts",
237 "fhTest_invmass_ANNcuts;ann;invariant mass "
238 "of 3 e^{#pm} + 1 particle in GeV/c^{2}",
239 10, 0, 10, invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
240 fHistoList_test.push_back(fhTest_invmass_ANNcuts);
241
242 fhTest_phaseSpace_pi0 =
243 new TH2D("fhTest_phaseSpace_pi0", "fhTest_phaseSpace_pi0;p_{t} in GeV/c;rapidity y", 240, -2., 10., 270, -2., 7.);
244 fhTest_phaseSpace_eta =
245 new TH2D("fhTest_phaseSpace_eta", "fhTest_phaseSpace_eta;p_{t} in GeV/c;rapidity y", 240, -2., 10., 270, -2., 7.);
246 fHistoList_test.push_back(fhTest_phaseSpace_pi0);
247 fHistoList_test.push_back(fhTest_phaseSpace_eta);
248
249 fhTest_eventMixing_STSonly_2p2 = new TH1D("fhTest_eventMixing_STSonly_2p2",
250 "fhTest_eventMixing_STSonly_2p2; invariant mass of 4 e^{#pm} in "
251 "GeV/c^{2}; #",
252 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
253 fHistoList_test.push_back(fhTest_eventMixing_STSonly_2p2);
254 fhTest_eventMixing_STSonly_3p1 = new TH1D("fhTest_eventMixing_STSonly_3p1",
255 "fhTest_eventMixing_STSonly_3p1; invariant mass of 4 e^{#pm} in "
256 "GeV/c^{2}; #",
257 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
258 fHistoList_test.push_back(fhTest_eventMixing_STSonly_3p1);
259 fhTest_eventMixing_STSonly_4p0 = new TH1D("fhTest_eventMixing_STSonly_4p0",
260 "fhTest_eventMixing_STSonly_4p0; invariant mass of 4 e^{#pm} in "
261 "GeV/c^{2}; #",
262 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
263 fHistoList_test.push_back(fhTest_eventMixing_STSonly_4p0);
264
265 fhTest_eventMixing_3p1 =
266 new TH1D("fhTest_eventMixing_3p1", "fhTest_eventMixing_3p1; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
267 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
268 fHistoList_test.push_back(fhTest_eventMixing_3p1);
269 fhTest_eventMixing_3p1_pCut =
270 new TH1D("fhTest_eventMixing_3p1_pCut", "fhTest_eventMixing_3p1_pCut; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
271 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
272 fHistoList_test.push_back(fhTest_eventMixing_3p1_pCut);
273 fhTest_eventMixing_3p1_GGcut = new TH1D("fhTest_eventMixing_3p1_GGcut",
274 "fhTest_eventMixing_3p1_GGcut; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
275 invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
276 fHistoList_test.push_back(fhTest_eventMixing_3p1_GGcut);
277
278 fhTest_eventMixing_3p1_ANNcuts = new TH2D("fhTest_eventMixing_3p1_ANNcuts",
279 "fhTest_eventMixing_3p1_ANNcuts;ann;invariant mass of 4 e^{#pm} "
280 "in GeV/c^{2}",
281 10, 0, 10, invmassSpectra_nof, invmassSpectra_start, invmassSpectra_end);
282 fHistoList_test.push_back(fhTest_eventMixing_3p1_ANNcuts);
283}
284
285
286void CbmAnaConversionTest::Finish()
287{
288 //gDirectory->cd("analysis-conversion");
289 gDirectory->mkdir("Test");
290 gDirectory->cd("Test");
291 for (UInt_t i = 0; i < fHistoList_test.size(); i++) {
292 fHistoList_test[i]->Write();
293 }
294 gDirectory->cd("..");
295}
296
297
298void CbmAnaConversionTest::Exec()
299{
300 cout << "CbmAnaConversionTest: Exec()" << endl;
301
302 globalEventNo++;
303
304 if (fPrimVertex != NULL) { fKFVertex = CbmKFVertex(*fPrimVertex); }
305 else {
306 Fatal("CbmAnaConversionTest::Exec", "No PrimaryVertex array!");
307 }
308
309 fElectrons_mcid.clear();
310 fElectrons_gtid.clear();
311 fElectrons_richInd.clear();
312 fElectrons_pi0mcid.clear();
313 fElectrons_same.clear();
314
315
316 Int_t nGTracks = fGlobalTracks->GetEntriesFast();
317 for (Int_t i = 0; i < nGTracks; i++) {
318 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
319 if (NULL == gTrack) continue;
320 int stsInd = gTrack->GetStsTrackIndex();
321 int richInd = gTrack->GetRichRingIndex();
322 if (stsInd < 0) continue;
323 if (richInd < 0) continue;
324
325 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
326 if (stsTrack == NULL) continue;
327
328 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
329 if (stsMatch == NULL) continue;
330 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
331 if (stsMcTrackId < 0) continue;
332 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
333 if (mcTrack1 == NULL) continue;
334
335 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
336 if (richMatch == NULL) continue;
337 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
338 if (richMcTrackId < 0) continue;
339 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
340 if (mcTrack2 == NULL) continue;
341
342 int pdg = TMath::Abs(mcTrack1->GetPdgCode());
343 int motherId = mcTrack1->GetMotherId();
344
345 if (motherId < 0) continue;
346 CbmMCTrack* mcTrack1mother = (CbmMCTrack*) fMcTracks->At(motherId);
347 if (mcTrack1mother == NULL) continue;
348 int grandmotherId = mcTrack1mother->GetMotherId();
349 if (grandmotherId < 0) continue;
350 CbmMCTrack* mcTrack1grandmother = (CbmMCTrack*) fMcTracks->At(grandmotherId);
351 if (mcTrack1grandmother == NULL) continue;
352
353 int pdg_mother = TMath::Abs(mcTrack1mother->GetPdgCode());
354 int pdg_grandmother = TMath::Abs(mcTrack1grandmother->GetPdgCode());
355
356 cout << "CbmAnaConversionTest: pdg particle/mother/grandmother: " << pdg << " / " << pdg_mother << " / "
357 << pdg_grandmother << endl;
358
359 Int_t sameRichSts = (stsMcTrackId == richMcTrackId);
360
361 if (TMath::Abs(pdg) == 11 && (pdg_mother == 111 || pdg_grandmother == 111)) {
362 fElectrons_richInd.push_back(richInd);
363 fElectrons_gtid.push_back(i);
364 fElectrons_mcid.push_back(stsMcTrackId);
365 if (pdg_mother == 111) fElectrons_pi0mcid.push_back(motherId);
366 if (pdg_grandmother == 111) fElectrons_pi0mcid.push_back(grandmotherId);
367 fElectrons_same.push_back(sameRichSts);
368 }
369 }
370
371 GetNofRichElectrons();
372
373
374 DoSTSonlyAnalysis();
375}
376
377
378void CbmAnaConversionTest::GetNofRichElectrons()
379{
380
381 std::multimap<int, int> electronMap;
382 for (unsigned int i = 0; i < fElectrons_pi0mcid.size(); i++) {
383 electronMap.insert(std::pair<int, int>(fElectrons_pi0mcid[i], i));
384 }
385
386 int check = 0;
387 int nofRich = 0;
388 for (std::map<int, int>::iterator it = electronMap.begin(); it != electronMap.end(); ++it) {
389 if (it == electronMap.begin()) {
390 check = 1;
391 if (fElectrons_same[it->second] > 0) nofRich++;
392 }
393 if (it != electronMap.begin()) {
394 std::map<int, int>::iterator zwischen = it;
395 zwischen--;
396 int id = it->first;
397 int id_old = zwischen->first;
398 if (id == id_old) {
399 check++;
400 if (fElectrons_same[it->second] > 0) nofRich++;
401 if (check > 3) {
402 std::map<int, int>::iterator zwischen2 = zwischen--;
403 std::map<int, int>::iterator zwischen3 = zwischen2--;
404 Double_t invmass = CalcInvMass(it->second, zwischen->second, zwischen2->second, zwischen3->second);
405 if (invmass > 0) fhElectronsTest_invmass->Fill(invmass);
406 if (nofRich == 4) {}
407 }
408 }
409 else {
410 fElectrons_nofPerPi0->Fill(check);
411 fElectrons_nofPerPi0_withRichInd->Fill(nofRich);
412 check = 1;
413 nofRich = 0;
414 if (fElectrons_same[it->second] > 0) nofRich++;
415 }
416 }
417 }
418}
419
420
421Double_t CbmAnaConversionTest::CalcInvMass(Int_t e1, Int_t e2, Int_t e3, Int_t e4)
422{
423
424 CbmMCTrack* mcTrack_e1 = (CbmMCTrack*) fMcTracks->At(fElectrons_mcid[e1]);
425 CbmMCTrack* mcTrack_e2 = (CbmMCTrack*) fMcTracks->At(fElectrons_mcid[e2]);
426 CbmMCTrack* mcTrack_e3 = (CbmMCTrack*) fMcTracks->At(fElectrons_mcid[e3]);
427 CbmMCTrack* mcTrack_e4 = (CbmMCTrack*) fMcTracks->At(fElectrons_mcid[e4]);
428
429 if (mcTrack_e1->GetPdgCode() + mcTrack_e2->GetPdgCode() + mcTrack_e3->GetPdgCode() + mcTrack_e4->GetPdgCode())
430 return -1;
431
432 TLorentzVector lorentz1;
433 mcTrack_e1->Get4Momentum(lorentz1);
434 TLorentzVector lorentz2;
435 mcTrack_e2->Get4Momentum(lorentz2);
436 TLorentzVector lorentz3;
437 mcTrack_e3->Get4Momentum(lorentz3);
438 TLorentzVector lorentz4;
439 mcTrack_e4->Get4Momentum(lorentz4);
440
441 TLorentzVector sum;
442 sum = lorentz1 + lorentz2 + lorentz3 + lorentz4;
443
444 return sum.Mag();
445}
446
447
448void CbmAnaConversionTest::DoSTSonlyAnalysis()
449{
450 fVector_AllMomenta.clear();
451
452 fVector_gt.clear();
453 fVector_momenta.clear();
454 fVector_chi.clear();
455 fVector_gtIndex.clear();
456 fVector_richIndex.clear();
457 fVector_mcIndex.clear();
458 fVector_withRichSignal.clear();
459 fVector_reconstructedPhotons_FromSTSandRICH.clear();
460
461 fVector_electronRICH_momenta.clear();
462 fVector_electronRICH_gt.clear();
463 fVector_electronRICH_gtIndex.clear();
464 fVector_electronRICH_mcIndex.clear();
465 fVector_electronRICH_reconstructedPhotons.clear();
466
467 fVector_reconstructedPhotons_STSonly.clear();
468
469
470 if (globalEventNo % 200 == 0) {
471 MixedEventTest_3p1();
472 fMixedTest_3p1_photons.clear();
473 fMixedTest_3p1_eventno.clear();
474 fMixedTest_3p1_combined.clear();
475 fMixedTest_3p1_ann.clear();
476 }
477
478 if (globalEventNo % 20 == 0) {
479 MixedEventTest_STSonly();
480 fMixedTest_STSonly_photons.clear();
481 fMixedTest_STSonly_eventno.clear();
482 fMixedTest_STSonly_hasRichInd.clear();
483 }
484
485
486 Int_t nofRICHelectrons = 0;
487
488 Int_t ngTracks = fGlobalTracks->GetEntriesFast();
489 fhTest_ParticlesPerEvent->Fill(ngTracks);
490 for (Int_t i = 0; i < ngTracks; i++) {
491 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
492 if (NULL == gTrack) continue;
493 int stsInd = gTrack->GetStsTrackIndex();
494 int richInd = gTrack->GetRichRingIndex();
495 // int trdInd = gTrack->GetTrdTrackIndex();
496 // int tofInd = gTrack->GetTofHitIndex();
497
498 if (stsInd < 0) continue;
499 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
500 if (stsTrack == NULL) continue;
501 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
502 if (stsMatch == NULL) continue;
503 if (stsMatch->GetNofLinks() <= 0) continue;
504 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
505 if (stsMcTrackId < 0) continue;
506 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
507 if (mcTrack1 == NULL) continue;
508
509
510 // calculate refitted momenta at primary vertex
511 TVector3 refittedMomentum;
512 CbmL1PFFitter fPFFitter;
513 vector<CbmStsTrack> stsTracks;
514 stsTracks.resize(1);
515 stsTracks[0] = *stsTrack;
516 vector<CbmL1PFFitter::PFFieldRegion> vField;
517 vector<float> chiPrim;
518 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
519 //cand.chi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
520 //cand.chi2Prim = chiPrim[0];
521 const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();
522 vtxTrack->Momentum(refittedMomentum);
523 //float result_chi = chiPrim[0];
524
525
526 // Doing refit of momenta with electron assumption
527 CbmL1PFFitter fPFFitter_electron;
528 vector<CbmStsTrack> stsTracks_electron;
529 stsTracks_electron.resize(1);
530 stsTracks_electron[0] = *stsTrack;
531 vector<CbmL1PFFitter::PFFieldRegion> vField_electron;
532 vector<float> chiPrim_electron;
533 vector<int> pidHypo_electron;
534 pidHypo_electron.push_back(11);
535 fPFFitter_electron.Fit(stsTracks_electron, pidHypo_electron);
536 fPFFitter_electron.GetChiToVertex(stsTracks_electron, vField_electron, chiPrim_electron, fKFVertex, 3e6);
537
538 TVector3 refittedMomentum_electron;
539 const FairTrackParam* vtxTrack_electron = stsTracks_electron[0].GetParamFirst();
540 vtxTrack_electron->Momentum(refittedMomentum_electron);
541 float result_chi_electron = chiPrim_electron[0];
542 //float result_ndf_electron = stsTracks_electron[0].GetNDF();
543 // Double_t stsHits = stsTrack->GetNofHits();
544
545
546 fVector_AllMomenta.push_back(refittedMomentum_electron);
547
548 Double_t chiCut = 0;
549 chiCut = CbmAnaConversionCutSettings::CalcChiCut(refittedMomentum_electron.Perp());
550
551 if (result_chi_electron > chiCut) continue;
552
553 //if(stsHits < 10) continue;
554 //if(stsHits > 9) {
555 fVector_momenta.push_back(refittedMomentum_electron);
556 fVector_chi.push_back(result_chi_electron);
557 fVector_gtIndex.push_back(i);
558 fVector_gt.push_back(gTrack);
559 fVector_richIndex.push_back(richInd);
560 fVector_mcIndex.push_back(stsMcTrackId);
561 //}
562
563 Bool_t WithRichSignal = false;
564
565 /* if (richInd < 0) continue;
566 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*)fRichRingMatches->At(richInd);
567 if (richMatch == NULL) continue;
568 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
569 if (richMcTrackId < 0) continue;
570 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
571 if (mcTrack2 == NULL) continue;
572
573 Bool_t electron_rich = CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, refittedMomentum_electron.Mag());
574 if(electron_rich) {
575 fVector_electronRICH_momenta.push_back(refittedMomentum_electron);
576 fVector_electronRICH_gt.push_back(gTrack);
577 fVector_electronRICH_gtIndex.push_back(i);
578 nofRICHelectrons++;
579 WithRichSignal = true;
580 }
581*/
582 if (richInd >= 0) {
583 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
584 if (richMatch != NULL) {
585 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
586 if (richMcTrackId >= 0) {
587 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
588 if (mcTrack2 != NULL) {
589
590
591 Bool_t electron_rich =
592 CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, refittedMomentum_electron.Mag());
593 if (electron_rich) {
594 fVector_electronRICH_momenta.push_back(refittedMomentum_electron);
595 fVector_electronRICH_gt.push_back(gTrack);
596 fVector_electronRICH_gtIndex.push_back(i);
597 fVector_electronRICH_mcIndex.push_back(stsMcTrackId);
598 nofRICHelectrons++;
599 WithRichSignal = true;
600 }
601 }
602 }
603 }
604 }
605 fVector_withRichSignal.push_back(WithRichSignal);
606 }
607
608 fhTest_RICHelectronsPerEvent->Fill(nofRICHelectrons);
609
610 // combination of 3 identified electrons in RICH and 1 track from STS with no RICH signal
611 CombineElectrons_FromRICH();
612 CombineElectrons_FromSTSandRICH();
613 CombinePhotons();
614
615 // starting point: track with STS signal only
616 CombineElectrons_STSonly();
617 CombinePhotons_STSonly();
618}
619
620
621void CbmAnaConversionTest::CombineElectrons_FromSTSandRICH()
622{
623 Int_t nof_sts = fVector_momenta.size();
624 Int_t nof_rich = fVector_electronRICH_momenta.size();
625 cout << "CbmAnaConversionTest: CombineElectrons_FromSTSandRICH, nof sts/rich - " << nof_sts << " / " << nof_rich
626 << endl;
627 Int_t nofPhotons = 0;
628 if (nof_sts + nof_rich >= 2) {
629 for (int a = 0; a < nof_sts; a++) {
630 for (int b = 0; b < nof_rich; b++) {
631 Int_t check1 = (fVector_gt[a]->GetParamLast()->GetQp()
632 > 0); // positive or negative charge (qp = charge over momentum ratio)
633 Int_t check2 = (fVector_electronRICH_gt[b]->GetParamLast()->GetQp() > 0);
634 Int_t test = check1 + check2;
635 if (test != 1) continue; // need one electron and one positron
636 if (fVector_gtIndex[a] == fVector_electronRICH_gtIndex[b]) continue;
637 if (fVector_richIndex[a] >= 0) continue; // 4th lepton should have no signal in RICH, only in STS
638
639 //LmvmKinePar params1 = CalculateKinematicParamsReco(fVector_momenta[a], fVector_electronRICH_momenta[b]);
640 CbmAnaConversionKinematicParams paramsTest = CbmAnaConversionKinematicParams::KinematicParams_2particles_Reco(
641 fVector_momenta[a], fVector_electronRICH_momenta[b]);
642
643 // opening angle cut depending on pt of e+e- pair
644 Double_t openingAngleCut = 1.8 - 0.6 * paramsTest.fPt;
645 //Double_t openingAngleCut = 1.5 - 0.5 * paramsTest.fPt;
646
647 Double_t invMassCut = 0.03;
648
649 Int_t IsPhoton_openingAngle1 = (paramsTest.fAngle < openingAngleCut);
650 Int_t IsPhoton_invMass1 = (paramsTest.fMinv < invMassCut);
651
652
653 if (IsPhoton_openingAngle1 && IsPhoton_invMass1) {
654 nofPhotons++;
655 vector<int> pair; // = {a, b};
656 pair.push_back(a);
657 pair.push_back(b);
658 fVector_reconstructedPhotons_FromSTSandRICH.push_back(pair);
659 //fhElectrons_invmass_cut->Fill(params1.fMinv);
660
661 vector<TVector3> pairmomenta;
662 pairmomenta.push_back(fVector_momenta[a]);
663 pairmomenta.push_back(fVector_electronRICH_momenta[b]);
664
665 vector<Double_t> ann;
666 ann.push_back(ElectronANNvalue(fVector_gtIndex[a], fVector_momenta[a].Mag()));
667 ann.push_back(ElectronANNvalue(fVector_electronRICH_gtIndex[b], fVector_electronRICH_momenta[b].Mag()));
668
669 fMixedTest_3p1_photons.push_back(pairmomenta);
670 fMixedTest_3p1_eventno.push_back(globalEventNo);
671 fMixedTest_3p1_combined.push_back(1);
672 fMixedTest_3p1_ann.push_back(ann);
673 }
674 }
675 }
676 }
677 fhTest_PhotonsPerEvent_STSandRICH->Fill(nofPhotons);
678 cout << "CbmAnaConversionTest: CombineElectronsFromSTSandRICH: Crosscheck - "
679 "nof reconstructed photons: "
680 << nofPhotons << endl;
681}
682
683
684void CbmAnaConversionTest::CombineElectrons_FromRICH()
685{
686 Int_t nof = fVector_electronRICH_momenta.size();
687 cout << "CbmAnaConversionTest: CombineElectronsFromRICH, nof - " << nof << endl;
688 Int_t nofPhotons = 0;
689 if (nof >= 2) {
690 for (int a = 0; a < nof - 1; a++) {
691 for (int b = a + 1; b < nof; b++) {
692 Int_t check1 = (fVector_electronRICH_gt[a]->GetParamLast()->GetQp()
693 > 0); // positive or negative charge (qp = charge over momentum ratio)
694 Int_t check2 = (fVector_electronRICH_gt[b]->GetParamLast()->GetQp() > 0);
695 Int_t test = check1 + check2;
696 if (test != 1) continue; // need one electron and one positron
697 if (fVector_electronRICH_gtIndex[a] == fVector_electronRICH_gtIndex[b]) continue;
698
699 //LmvmKinePar params1 = CalculateKinematicParamsReco(fVector_electronRICH_momenta[a], fVector_electronRICH_momenta[b]);
700 CbmAnaConversionKinematicParams paramsTest = CbmAnaConversionKinematicParams::KinematicParams_2particles_Reco(
701 fVector_electronRICH_momenta[a], fVector_electronRICH_momenta[b]);
702
703 // opening angle cut depending on pt of e+e- pair
704 Double_t openingAngleCut = 1.8 - 0.6 * paramsTest.fPt;
705 //Double_t openingAngleCut = 1.5 - 0.5 * paramsTest.fPt;
706
707 Double_t invMassCut = 0.03;
708
709 Int_t IsPhoton_openingAngle1 = (paramsTest.fAngle < openingAngleCut);
710 Int_t IsPhoton_invMass1 = (paramsTest.fMinv < invMassCut);
711
712
713 if (IsPhoton_openingAngle1 && IsPhoton_invMass1) {
714 nofPhotons++;
715 vector<int> pair; // = {a, b};
716 pair.push_back(a);
717 pair.push_back(b);
718 fVector_electronRICH_reconstructedPhotons.push_back(pair);
719 //fhElectrons_invmass_cut->Fill(params1.fMinv);
720
721 vector<TVector3> pairmomenta;
722 pairmomenta.push_back(fVector_electronRICH_momenta[a]);
723 pairmomenta.push_back(fVector_electronRICH_momenta[b]);
724
725 vector<Double_t> ann;
726 ann.push_back(ElectronANNvalue(fVector_electronRICH_gtIndex[a], fVector_electronRICH_momenta[a].Mag()));
727 ann.push_back(ElectronANNvalue(fVector_electronRICH_gtIndex[b], fVector_electronRICH_momenta[b].Mag()));
728
729 fMixedTest_3p1_photons.push_back(pairmomenta);
730 fMixedTest_3p1_eventno.push_back(globalEventNo);
731 fMixedTest_3p1_combined.push_back(0);
732 fMixedTest_3p1_ann.push_back(ann);
733 }
734 }
735 }
736 }
737 fhTest_PhotonsPerEvent_RICHonly->Fill(nofPhotons);
738 cout << "CbmAnaConversionTest: CombineElectronsFromRICH: Crosscheck - nof "
739 "reconstructed photons: "
740 << nofPhotons << endl;
741}
742
743
744void CbmAnaConversionTest::CombinePhotons()
745{
746 Int_t nof_STSandRICH = fVector_reconstructedPhotons_FromSTSandRICH.size();
747 Int_t nof_RICH = fVector_electronRICH_reconstructedPhotons.size();
748 cout << "CbmAnaConversionTest: CombinePhotons, nof - " << nof_STSandRICH << "/" << nof_RICH << endl;
749 Int_t nofPi0 = 0;
750 if (nof_STSandRICH + nof_RICH >= 2) {
751 for (int a = 0; a < nof_STSandRICH; a++) {
752 for (int b = 0; b < nof_RICH; b++) {
753 Int_t electron11 = fVector_reconstructedPhotons_FromSTSandRICH[a][0]; // track with STS signal only
754 Int_t electron12 = fVector_reconstructedPhotons_FromSTSandRICH[a][1]; // track with STS + RICH signal
755 Int_t electron21 = fVector_electronRICH_reconstructedPhotons[b][0]; // track with STS + RICH signal
756 Int_t electron22 = fVector_electronRICH_reconstructedPhotons[b][1]; // track with STS + RICH signal
757
758 Int_t gtIndex11 = fVector_gtIndex[electron11];
759 Int_t gtIndex12 = fVector_electronRICH_gtIndex[electron12];
760 Int_t gtIndex21 = fVector_electronRICH_gtIndex[electron21];
761 Int_t gtIndex22 = fVector_electronRICH_gtIndex[electron22];
762
763 if (gtIndex11 == gtIndex12 || gtIndex11 == gtIndex21 || gtIndex11 == gtIndex22 || gtIndex12 == gtIndex21
764 || gtIndex12 == gtIndex22 || gtIndex21 == gtIndex22) {
765 //if(electron12 == electron21 || electron12 == electron22) {
766 cout << "CbmAnaConversionTest: Test_DoubleIndex." << endl;
767 continue;
768 }
769
770 //Double_t invmass = Invmass_4particlesRECO(fVector_momenta[electron11], fVector_electronRICH_momenta[electron12], fVector_electronRICH_momenta[electron21], fVector_electronRICH_momenta[electron22]);
771
772 CbmAnaConversionKinematicParams paramsTest = CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
773 fVector_momenta[electron11], fVector_electronRICH_momenta[electron12],
774 fVector_electronRICH_momenta[electron21], fVector_electronRICH_momenta[electron22]);
775
776 Double_t invmass = paramsTest.fMinv;
777
778 fhTest_invmass->Fill(invmass);
779
780 if (fVector_momenta[electron11].Mag() < 0.6) { fhTest_invmass_pCut->Fill(invmass); }
781
782 Double_t OpeningAngleGG = CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma(
783 fVector_momenta[electron11], fVector_electronRICH_momenta[electron12],
784 fVector_electronRICH_momenta[electron21], fVector_electronRICH_momenta[electron22]);
785 if (OpeningAngleGG > 15 && OpeningAngleGG < 45) fhTest_invmass_GGcut->Fill(invmass);
786
787 //Double_t pt = Pt_4particlesRECO(momenta[electron11], momenta[electron12], momenta[electron21], momenta[electron22]);
788 //Double_t rap = Rap_4particlesRECO(momenta[electron11], momenta[electron12], momenta[electron21], momenta[electron22]);
789
790 //Double_t opening_angle = OpeningAngleBetweenPhotons(momenta, reconstructedPhotons[a], reconstructedPhotons[b]);
791
792
793 //LmvmKinePar params1 = CalculateKinematicParams_4particles(momenta[electron11], momenta[electron12], momenta[electron21], momenta[electron22]);
794
795
796 Double_t ANNe11 = ElectronANNvalue(gtIndex11, fVector_momenta[electron11].Mag());
797 Double_t ANNe12 = ElectronANNvalue(gtIndex12, fVector_electronRICH_momenta[electron12].Mag());
798 Double_t ANNe21 = ElectronANNvalue(gtIndex21, fVector_electronRICH_momenta[electron21].Mag());
799 Double_t ANNe22 = ElectronANNvalue(gtIndex22, fVector_electronRICH_momenta[electron22].Mag());
800
801 cout << "CbmAnaConversionTest: CombinePhotons, anns: " << ANNe11 << "/" << ANNe12 << "/" << ANNe21 << "/"
802 << ANNe22 << endl;
803
804 // ann-check only for those tracks, which have a signal in the RICH, i.e. not the first one
805 if (ANNe12 > -1 && ANNe21 > -1 && ANNe22 > -1) fhTest_invmass_ANNcuts->Fill(1, invmass);
806 if (ANNe12 > -0.9 && ANNe21 > -0.9 && ANNe22 > -0.9) fhTest_invmass_ANNcuts->Fill(2, invmass);
807 if (ANNe12 > -0.8 && ANNe21 > -0.8 && ANNe22 > -0.8) fhTest_invmass_ANNcuts->Fill(3, invmass);
808 if (ANNe12 > -0.7 && ANNe21 > -0.7 && ANNe22 > -0.7) fhTest_invmass_ANNcuts->Fill(4, invmass);
809 if (ANNe12 > -0.6 && ANNe21 > -0.6 && ANNe22 > -0.6) fhTest_invmass_ANNcuts->Fill(5, invmass);
810 if (ANNe12 > -0.5 && ANNe21 > -0.5 && ANNe22 > -0.5) fhTest_invmass_ANNcuts->Fill(6, invmass);
811 if (ANNe12 > 0.0 && ANNe21 > 0.0 && ANNe22 > 0.0) fhTest_invmass_ANNcuts->Fill(7, invmass);
812
813
814 cout << "CbmAnaConversionTest: CombinePhotons, photon combined! " << invmass << "/" << gtIndex11 << "/"
815 << gtIndex12 << "/" << gtIndex21 << "/" << gtIndex22 << endl;
816 nofPi0++;
817
818
819 // MC-TRUE crosschecks
820 CbmMCTrack* mctrack11 = (CbmMCTrack*) fMcTracks->At(fVector_mcIndex[electron11]); // mctracks of four leptons
821 CbmMCTrack* mctrack12 = (CbmMCTrack*) fMcTracks->At(fVector_electronRICH_mcIndex[electron12]);
822 CbmMCTrack* mctrack21 = (CbmMCTrack*) fMcTracks->At(fVector_electronRICH_mcIndex[electron21]);
823 CbmMCTrack* mctrack22 = (CbmMCTrack*) fMcTracks->At(fVector_electronRICH_mcIndex[electron22]);
824
825 // Int_t pdg11 = mctrack11->GetPdgCode(); // pdg codes of four leptons
826 // Int_t pdg12 = mctrack12->GetPdgCode();
827 // Int_t pdg21 = mctrack21->GetPdgCode();
828 // Int_t pdg22 = mctrack22->GetPdgCode();
829
830 Int_t motherId11 = mctrack11->GetMotherId(); // motherIDs of four leptons
831 Int_t motherId12 = mctrack12->GetMotherId();
832 Int_t motherId21 = mctrack21->GetMotherId();
833 Int_t motherId22 = mctrack22->GetMotherId();
834
835 CbmMCTrack* mothermctrack11 = NULL; // mctracks of mother particles of the four leptons
836 CbmMCTrack* mothermctrack12 = NULL;
837 CbmMCTrack* mothermctrack21 = NULL;
838 CbmMCTrack* mothermctrack22 = NULL;
839 if (motherId11 > 0) mothermctrack11 = (CbmMCTrack*) fMcTracks->At(motherId11);
840 if (motherId11 > 0) mothermctrack12 = (CbmMCTrack*) fMcTracks->At(motherId12);
841 if (motherId11 > 0) mothermctrack21 = (CbmMCTrack*) fMcTracks->At(motherId21);
842 if (motherId11 > 0) mothermctrack22 = (CbmMCTrack*) fMcTracks->At(motherId22);
843
844 Int_t motherpdg11 = -2; // pdg codes of the mother particles
845 // Int_t motherpdg12 = -2;
846 Int_t motherpdg21 = -2;
847 // Int_t motherpdg22 = -2;
848 if (mothermctrack11 != NULL) motherpdg11 = mothermctrack11->GetPdgCode();
849 // if(mothermctrack12 != NULL) motherpdg12 = mothermctrack12->GetPdgCode();
850 if (mothermctrack21 != NULL) motherpdg21 = mothermctrack21->GetPdgCode();
851 // if(mothermctrack22 != NULL) motherpdg22 = mothermctrack22->GetPdgCode();
852
853 Int_t grandmotherId11 = -2; // grandmotherIDs of four leptons
854 Int_t grandmotherId12 = -2;
855 Int_t grandmotherId21 = -2;
856 Int_t grandmotherId22 = -2;
857 if (mothermctrack11 != NULL) grandmotherId11 = mothermctrack11->GetMotherId();
858 if (mothermctrack12 != NULL) grandmotherId12 = mothermctrack12->GetMotherId();
859 if (mothermctrack21 != NULL) grandmotherId21 = mothermctrack21->GetMotherId();
860 if (mothermctrack22 != NULL) grandmotherId22 = mothermctrack22->GetMotherId();
861
862
863 Int_t sameGrandmothers1 = 0;
864 Int_t sameGrandmothers2 = 0;
865 Int_t sameGrandmothers3 = 0;
866 Int_t sameGrandmothers4 = 0;
867 if (grandmotherId11 == grandmotherId12) sameGrandmothers1++;
868 if (grandmotherId11 == grandmotherId21) sameGrandmothers1++;
869 if (grandmotherId11 == grandmotherId22) sameGrandmothers1++;
870 if (grandmotherId12 == grandmotherId11) sameGrandmothers2++;
871 if (grandmotherId12 == grandmotherId21) sameGrandmothers2++;
872 if (grandmotherId12 == grandmotherId22) sameGrandmothers2++;
873 if (grandmotherId21 == grandmotherId11) sameGrandmothers3++;
874 if (grandmotherId21 == grandmotherId12) sameGrandmothers3++;
875 if (grandmotherId21 == grandmotherId22) sameGrandmothers3++;
876 if (grandmotherId22 == grandmotherId11) sameGrandmothers4++;
877 if (grandmotherId22 == grandmotherId12) sameGrandmothers4++;
878 if (grandmotherId22 == grandmotherId21) sameGrandmothers4++;
879 Int_t sameGrandmothersSum = sameGrandmothers1 + sameGrandmothers2 + sameGrandmothers3 + sameGrandmothers4;
880
881 Int_t sameMothers1 = 0;
882 Int_t sameMothers2 = 0;
883 Int_t sameMothers3 = 0;
884 Int_t sameMothers4 = 0;
885 if (motherId11 == motherId12) sameMothers1++;
886 if (motherId11 == motherId21) sameMothers1++;
887 if (motherId11 == motherId22) sameMothers1++;
888 if (motherId12 == motherId11) sameMothers2++;
889 if (motherId12 == motherId21) sameMothers2++;
890 if (motherId12 == motherId22) sameMothers2++;
891 if (motherId21 == motherId11) sameMothers3++;
892 if (motherId21 == motherId12) sameMothers3++;
893 if (motherId21 == motherId22) sameMothers3++;
894 if (motherId22 == motherId11) sameMothers4++;
895 if (motherId22 == motherId12) sameMothers4++;
896 if (motherId22 == motherId21) sameMothers4++;
897 Int_t sameMothersSum = sameMothers1 + sameMothers2 + sameMothers3 + sameMothers4;
898
899
900 if (
901 motherId11 == motherId12
902 && motherId21
903 == motherId22) { // both combined e+e- pairs come from the same mother (which can be gamma, pi0, or whatever)
904 fhTest_invmass_MCcutAll->Fill(1, invmass);
905 if (TMath::Abs(motherpdg11) == 22 && TMath::Abs(motherpdg21) == 22) {
906 fhTest_invmass_MCcutAll->Fill(2, invmass);
907 }
908 if (TMath::Abs(motherpdg11) == 22 && TMath::Abs(motherpdg21) == 22 && grandmotherId11 == grandmotherId21
909 && grandmotherId11 > 0) { // decay in to gg of pi0 and eta
910 fhTest_invmass_MCcutAll->Fill(3, invmass);
911 if (invmass < 0.3) fhTest_phaseSpace_pi0->Fill(paramsTest.fPt, paramsTest.fRapidity);
912 if (invmass > 0.3) fhTest_phaseSpace_eta->Fill(paramsTest.fPt, paramsTest.fRapidity);
913 }
914 if (TMath::Abs(motherpdg11) == 22 && TMath::Abs(motherpdg21) == 22 && grandmotherId11 != grandmotherId21) {
915 fhTest_invmass_MCcutAll->Fill(4, invmass);
916 }
917 if ((TMath::Abs(motherpdg11) == 22 && TMath::Abs(motherpdg21) == 111)
918 || (TMath::Abs(motherpdg11) == 111 && TMath::Abs(motherpdg21) == 22)) {
919 fhTest_invmass_MCcutAll->Fill(5, invmass);
920 if (grandmotherId11 == motherId21 || motherId11 == grandmotherId21) { // Dalitz decay of pi0
921 fhTest_invmass_MCcutAll->Fill(12, invmass);
922 fhTest_phaseSpace_pi0->Fill(paramsTest.fPt, paramsTest.fRapidity);
923 }
924 }
925 if (TMath::Abs(motherpdg11) == 111 && TMath::Abs(motherpdg21) == 111) {
926 fhTest_invmass_MCcutAll->Fill(6, invmass);
927 }
928 if ((TMath::Abs(motherpdg11) != 22 && TMath::Abs(motherpdg11) != 111)
929 || (TMath::Abs(motherpdg21) != 22 && TMath::Abs(motherpdg21) != 111)) {
930 fhTest_invmass_MCcutAll->Fill(7, invmass);
931 }
932 if (TMath::Abs(motherpdg11) != 22 && TMath::Abs(motherpdg11) != 111 && TMath::Abs(motherpdg21) != 22
933 && TMath::Abs(motherpdg21) != 111) {
934 fhTest_invmass_MCcutAll->Fill(8, invmass);
935 }
936 if (grandmotherId11 == grandmotherId21) { fhTest_invmass_MCcutAll->Fill(9, invmass); }
937 }
938 if ((motherId11 == motherId12 && motherId21 != motherId22)
939 || (motherId11 != motherId12 && motherId21 == motherId22)) {
940 fhTest_invmass_MCcutAll->Fill(10, invmass);
941 }
942 if (motherId11 != motherId12 && motherId21 != motherId22) {
943 fhTest_invmass_MCcutAll->Fill(11, invmass);
944 if (sameGrandmothersSum == 12) fhTest_invmass_MCcutAll->Fill(13, invmass);
945 if (sameGrandmothersSum == 6) fhTest_invmass_MCcutAll->Fill(14, invmass);
946 if (sameGrandmothersSum == 4) {
947 fhTest_invmass_MCcutAll->Fill(15, invmass);
948 if (grandmotherId11 < 0 || grandmotherId12 < 0 || grandmotherId21 < 0 || grandmotherId22 < 0) {
949 fhTest_invmass_MCcutAll->Fill(16, invmass);
950 }
951 if (grandmotherId11 == grandmotherId12) { fhTest_invmass_MCcutAll->Fill(17, invmass); }
952 if (grandmotherId11 != grandmotherId12) { fhTest_invmass_MCcutAll->Fill(18, invmass); }
953 }
954 if (sameGrandmothersSum == 2) fhTest_invmass_MCcutAll->Fill(19, invmass);
955 if (sameGrandmothersSum == 0) fhTest_invmass_MCcutAll->Fill(20, invmass);
956 fhTest_peakCheck1->Fill(sameGrandmothersSum);
957 fhTest_peakCheck2->Fill(sameMothersSum);
958 if ((grandmotherId11 < 0 || grandmotherId12 < 0 || grandmotherId21 < 0 || grandmotherId22 < 0)
959 && sameGrandmothersSum == 12) {
960 fhTest_peakCheck3->Fill(grandmotherId11);
961 fhTest_peakCheck3->Fill(grandmotherId12);
962 fhTest_peakCheck3->Fill(grandmotherId21);
963 fhTest_peakCheck3->Fill(grandmotherId22);
964 }
965 }
966 }
967 }
968 }
969 fhTest_ReconstructedPi0PerEvent->Fill(nofPi0);
970}
971
972
973LmvmKinePar CbmAnaConversionTest::CalculateKinematicParamsReco(const TVector3 electron1, const TVector3 electron2)
974{
975 LmvmKinePar params;
976
977 Double_t energyP = TMath::Sqrt(electron1.Mag2() + M2E);
978 TLorentzVector lorVecP(electron1, energyP);
979
980 Double_t energyM = TMath::Sqrt(electron2.Mag2() + M2E);
981 TLorentzVector lorVecM(electron2, energyM);
982
983 TVector3 momPair = electron1 + electron2;
984 Double_t energyPair = energyP + energyM;
985 Double_t ptPair = momPair.Perp();
986 Double_t pzPair = momPair.Pz();
987 Double_t yPair = 0.5 * TMath::Log((energyPair + pzPair) / (energyPair - pzPair));
988 Double_t anglePair = lorVecM.Angle(lorVecP.Vect());
989 Double_t theta = 180. * anglePair / TMath::Pi();
990 Double_t minv = 2. * TMath::Sin(anglePair / 2.) * TMath::Sqrt(electron1.Mag() * electron2.Mag());
991
992 params.fMomentumMag = momPair.Mag();
993 params.fPt = ptPair;
994 params.fRapidity = yPair;
995 params.fMinv = minv;
996 params.fAngle = theta;
997 return params;
998}
999
1000
1001Double_t CbmAnaConversionTest::Invmass_4particlesRECO(const TVector3 part1, const TVector3 part2, const TVector3 part3,
1002 const TVector3 part4)
1003// calculation of invariant mass from four electrons/positrons
1004{
1005 Double_t energy1 = TMath::Sqrt(part1.Mag2() + M2E);
1006 TLorentzVector lorVec1(part1, energy1);
1007
1008 Double_t energy2 = TMath::Sqrt(part2.Mag2() + M2E);
1009 TLorentzVector lorVec2(part2, energy2);
1010
1011 Double_t energy3 = TMath::Sqrt(part3.Mag2() + M2E);
1012 TLorentzVector lorVec3(part3, energy3);
1013
1014 Double_t energy4 = TMath::Sqrt(part4.Mag2() + M2E);
1015 TLorentzVector lorVec4(part4, energy4);
1016
1017 TLorentzVector sum;
1018 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
1019
1020 return sum.Mag();
1021}
1022
1023
1024void CbmAnaConversionTest::CombineElectrons_STSonly()
1025{
1026 Int_t nof = fVector_momenta.size();
1027 cout << "CbmAnaConversionTest: CombineElectrons_STSonly, nof - " << nof << endl;
1028 Int_t nofPhotons = 0;
1029 if (nof >= 2) {
1030 for (int a = 0; a < nof - 1; a++) {
1031 for (int b = a + 1; b < nof; b++) {
1032 Int_t check1 = (fVector_gt[a]->GetParamLast()->GetQp()
1033 > 0); // positive or negative charge (qp = charge over momentum ratio)
1034 Int_t check2 = (fVector_gt[b]->GetParamLast()->GetQp() > 0);
1035 Int_t test = check1 + check2;
1036 if (test != 1) continue; // need one electron and one positron
1037 if (fVector_gtIndex[a] == fVector_gtIndex[b]) continue;
1038
1039 //LmvmKinePar params1 = CalculateKinematicParamsReco(fVector_electronRICH_momenta[a], fVector_electronRICH_momenta[b]);
1040 CbmAnaConversionKinematicParams paramsTest =
1041 CbmAnaConversionKinematicParams::KinematicParams_2particles_Reco(fVector_momenta[a], fVector_momenta[b]);
1042
1043 // opening angle cut depending on pt of e+e- pair
1044 Double_t openingAngleCut = 1.8 - 0.6 * paramsTest.fPt;
1045 Double_t invMassCut = 0.03;
1046
1047 Int_t IsPhoton_openingAngle1 = (paramsTest.fAngle < openingAngleCut);
1048 Int_t IsPhoton_invMass1 = (paramsTest.fMinv < invMassCut);
1049
1050
1051 if (IsPhoton_openingAngle1 && IsPhoton_invMass1) {
1052 nofPhotons++;
1053 vector<int> pair; // = {a, b};
1054 pair.push_back(a);
1055 pair.push_back(b);
1056 fVector_reconstructedPhotons_STSonly.push_back(pair);
1057 //fhElectrons_invmass_cut->Fill(params1.fMinv);
1058
1059 // event mixing arrays
1060 vector<TVector3> pairmomenta;
1061 pairmomenta.push_back(fVector_momenta[a]);
1062 pairmomenta.push_back(fVector_momenta[b]);
1063 fMixedTest_STSonly_photons.push_back(pairmomenta);
1064 fMixedTest_STSonly_eventno.push_back(globalEventNo);
1065 vector<bool> pair_hasRichInd;
1066 //pair_hasRichInd.push_back(HasRichInd(fVector_gtIndex[a], a) );
1067 //pair_hasRichInd.push_back(HasRichInd(fVector_gtIndex[b], b) );
1068 pair_hasRichInd.push_back(fVector_withRichSignal[a]);
1069 pair_hasRichInd.push_back(fVector_withRichSignal[b]);
1070 fMixedTest_STSonly_hasRichInd.push_back(pair_hasRichInd);
1071 }
1072 }
1073 }
1074 }
1075}
1076
1077
1078void CbmAnaConversionTest::CombinePhotons_STSonly()
1079{
1080 Int_t nof = fVector_reconstructedPhotons_STSonly.size();
1081 cout << "CbmAnaConversionTest: CombinePhotons_STSonly, nof - " << nof << endl;
1082 Int_t nofPi0 = 0;
1083 Int_t nofPi0_0 = 0;
1084 Int_t nofPi0_1 = 0;
1085 Int_t nofPi0_2 = 0;
1086 Int_t nofPi0_3 = 0;
1087 Int_t nofPi0_4 = 0;
1088
1089
1090 if (nof >= 2) {
1091 for (int a = 0; a < nof - 1; a++) {
1092 for (int b = a + 1; b < nof; b++) {
1093 Int_t electron11 = fVector_reconstructedPhotons_STSonly[a][0]; // track with STS signal
1094 Int_t electron12 = fVector_reconstructedPhotons_STSonly[a][1]; // track with STS signal
1095 Int_t electron21 = fVector_reconstructedPhotons_STSonly[b][0]; // track with STS signal
1096 Int_t electron22 = fVector_reconstructedPhotons_STSonly[b][1]; // track with STS signal
1097
1098 Int_t gtIndex11 = fVector_gtIndex[electron11];
1099 Int_t gtIndex12 = fVector_gtIndex[electron12];
1100 Int_t gtIndex21 = fVector_gtIndex[electron21];
1101 Int_t gtIndex22 = fVector_gtIndex[electron22];
1102
1103 //Int_t richIndex11 = fVector_richIndex[electron11];
1104 //Int_t richIndex12 = fVector_richIndex[electron12];
1105 //Int_t richIndex21 = fVector_richIndex[electron21];
1106 //Int_t richIndex22 = fVector_richIndex[electron22];
1107
1108 if (gtIndex11 == gtIndex12 || gtIndex11 == gtIndex21 || gtIndex11 == gtIndex22 || gtIndex12 == gtIndex21
1109 || gtIndex12 == gtIndex22 || gtIndex21 == gtIndex22) {
1110 //if(electron12 == electron21 || electron12 == electron22) {
1111 //cout << "CbmAnaConversionTest - CombinePhotons_STSonly: Test_DoubleIndex." << endl;
1112 continue;
1113 }
1114
1115
1116 CbmAnaConversionKinematicParams paramsTest = CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(
1117 fVector_momenta[electron11], fVector_momenta[electron12], fVector_momenta[electron21],
1118 fVector_momenta[electron22]);
1119
1120 Double_t invmass = paramsTest.fMinv;
1121
1122
1123 //Int_t nofRICHindices = (HasRichInd(gtIndex11, electron11)) + (HasRichInd(gtIndex12, electron12)) + (HasRichInd(gtIndex21, electron21)) + (HasRichInd(gtIndex22, electron22));
1124 Int_t nofRICHindices = fVector_withRichSignal[electron11] + fVector_withRichSignal[electron12]
1125 + fVector_withRichSignal[electron21] + fVector_withRichSignal[electron22];
1126
1127 if (nofRICHindices == 0) {
1128 fhTest_invmass_RICHindex0->Fill(invmass);
1129 nofPi0_0++;
1130 }
1131 if (nofRICHindices == 1) {
1132 fhTest_invmass_RICHindex1->Fill(invmass);
1133 nofPi0_1++;
1134 }
1135 if (nofRICHindices == 2) {
1136 fhTest_invmass_RICHindex2->Fill(invmass);
1137 nofPi0_2++;
1138 }
1139 if (nofRICHindices == 3) {
1140 fhTest_invmass_RICHindex3->Fill(invmass);
1141 nofPi0_3++;
1142 //cout << "CbmAnaConversionTest: CombinePhotons_STSonly, photon combined! " << invmass << "/" << gtIndex11 << "/" << gtIndex12 << "/" << gtIndex21 << "/" << gtIndex22 << endl;
1143 }
1144 if (nofRICHindices == 4) {
1145 fhTest_invmass_RICHindex4->Fill(invmass);
1146 nofPi0_4++;
1147 }
1148
1149 nofPi0++;
1150 }
1151 }
1152 }
1153 cout << "CbmAnaConversionTest: CombinePhotons_STSonly, nofPi0 - " << nofPi0 << "/" << nofPi0_0 << "/" << nofPi0_1
1154 << "/" << nofPi0_2 << "/" << nofPi0_3 << "/" << nofPi0_4 << endl;
1155}
1156
1157
1158Bool_t CbmAnaConversionTest::HasRichInd(Int_t gtIndex, Int_t arrayIndex)
1159{
1160 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(gtIndex);
1161 if (NULL == gTrack) return false;
1162 int stsInd = gTrack->GetStsTrackIndex();
1163 int richInd = gTrack->GetRichRingIndex();
1164
1165 if (stsInd < 0) return false;
1166 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
1167 if (stsTrack == NULL) return false;
1168 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
1169 if (stsMatch == NULL) return false;
1170 if (stsMatch->GetNofLinks() <= 0) return false;
1171 int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
1172 if (stsMcTrackId < 0) return false;
1173 CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
1174 if (mcTrack1 == NULL) return false;
1175
1176 if (richInd < 0) return false;
1177 CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
1178 if (richMatch == NULL) return false;
1179 int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
1180 if (richMcTrackId < 0) return false;
1181 CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
1182 if (mcTrack2 == NULL) return false;
1183
1184 /*
1185 // calculate refitted momenta at primary vertex
1186 TVector3 refittedMomentum;
1187 CbmL1PFFitter fPFFitter;
1188 vector<CbmStsTrack> stsTracks;
1189 stsTracks.resize(1);
1190 stsTracks[0] = *stsTrack;
1191 vector<CbmL1PFFitter::PFFieldRegion> vField;
1192 vector<float> chiPrim;
1193 fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
1194 //cand.chi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
1195 //cand.chi2Prim = chiPrim[0];
1196 const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();
1197 vtxTrack->Momentum(refittedMomentum);
1198 //float result_chi = chiPrim[0];
1199
1200
1201
1202 // Doing refit of momenta with electron assumption
1203 CbmL1PFFitter fPFFitter_electron;
1204 vector<CbmStsTrack> stsTracks_electron;
1205 stsTracks_electron.resize(1);
1206 stsTracks_electron[0] = *stsTrack;
1207 vector<CbmL1PFFitter::PFFieldRegion> vField_electron;
1208 vector<float> chiPrim_electron;
1209 vector<int> pidHypo_electron;
1210 pidHypo_electron.push_back(11);
1211 fPFFitter_electron.Fit(stsTracks_electron, pidHypo_electron);
1212 fPFFitter_electron.GetChiToVertex(stsTracks_electron, vField_electron, chiPrim_electron, fKFVertex, 3e6);
1213
1214 TVector3 refittedMomentum_electron;
1215 const FairTrackParam* vtxTrack_electron = stsTracks_electron[0].GetParamFirst();
1216 vtxTrack_electron->Momentum(refittedMomentum_electron);
1217*/
1218
1219
1220 Bool_t electron_rich =
1221 CbmLitGlobalElectronId::GetInstance().IsRichElectron(gtIndex, fVector_momenta[arrayIndex].Mag());
1222 return electron_rich;
1223
1224 //return true;
1225}
1226
1227
1228void CbmAnaConversionTest::MixedEventTest_3p1()
1229// combines photons from two different events, taken from each time 200 events
1230{
1231 Int_t nof = fMixedTest_3p1_photons.size();
1232 //cout << "CbmAnaConversionRecoFull: MixedEventTest4 - nof entries " << nof << endl;
1233 for (Int_t a = 0; a < nof - 1; a++) {
1234 for (Int_t b = a + 1; b < nof; b++) {
1235 if (fMixedTest_3p1_eventno[a] == fMixedTest_3p1_eventno[b])
1236 continue; // to make sure that the photons are from two different events
1237 if (fMixedTest_3p1_combined[a] + fMixedTest_3p1_combined[b] != 1)
1238 continue; // to make sure, that the combination 3+1 is used
1239
1240 TVector3 e11 = fMixedTest_3p1_photons[a][0];
1241 TVector3 e12 = fMixedTest_3p1_photons[a][1];
1242 TVector3 e21 = fMixedTest_3p1_photons[b][0];
1243 TVector3 e22 = fMixedTest_3p1_photons[b][1];
1244
1245
1246 CbmAnaConversionKinematicParams params =
1247 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(e11, e12, e21, e22);
1248 fhTest_eventMixing_3p1->Fill(params.fMinv);
1249 //cout << "CbmAnaConversionRecoFull: MixedEventTest4(), event filled!, part" << endl;
1250
1251 if ((fMixedTest_3p1_combined[a] == 1 && e11.Mag() < 0.6)
1252 || (fMixedTest_3p1_combined[b] == 1 && e21.Mag() < 0.6)) {
1253 fhTest_eventMixing_3p1_pCut->Fill(params.fMinv);
1254 }
1255
1256
1257 Double_t OpeningAngleGG = CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma(e11, e12, e21, e22);
1258 if (OpeningAngleGG > 15 && OpeningAngleGG < 45) fhTest_eventMixing_3p1_GGcut->Fill(params.fMinv);
1259
1260
1261 //Double_t ANNe11 = fMixedTest_3p1_ann[a][0];
1262 Double_t ANNe12 = fMixedTest_3p1_ann[a][1];
1263 Double_t ANNe21 = fMixedTest_3p1_ann[b][0];
1264 Double_t ANNe22 = fMixedTest_3p1_ann[b][1];
1265
1266 if (ANNe12 > -1 && ANNe21 > -1 && ANNe22 > -1) fhTest_eventMixing_3p1_ANNcuts->Fill(1, params.fMinv);
1267 if (ANNe12 > -0.9 && ANNe21 > -0.9 && ANNe22 > -0.9) fhTest_eventMixing_3p1_ANNcuts->Fill(2, params.fMinv);
1268 if (ANNe12 > -0.8 && ANNe21 > -0.8 && ANNe22 > -0.8) fhTest_eventMixing_3p1_ANNcuts->Fill(3, params.fMinv);
1269 if (ANNe12 > -0.7 && ANNe21 > -0.7 && ANNe22 > -0.7) fhTest_eventMixing_3p1_ANNcuts->Fill(4, params.fMinv);
1270 if (ANNe12 > -0.6 && ANNe21 > -0.6 && ANNe22 > -0.6) fhTest_eventMixing_3p1_ANNcuts->Fill(5, params.fMinv);
1271 if (ANNe12 > -0.5 && ANNe21 > -0.5 && ANNe22 > -0.5) fhTest_eventMixing_3p1_ANNcuts->Fill(6, params.fMinv);
1272 if (ANNe12 > -0.0 && ANNe21 > -0.0 && ANNe22 > -0.0) fhTest_eventMixing_3p1_ANNcuts->Fill(7, params.fMinv);
1273 }
1274 }
1275}
1276
1277
1278void CbmAnaConversionTest::MixedEventTest_STSonly()
1279{
1280 Int_t nof = fMixedTest_STSonly_photons.size();
1281 for (Int_t a = 0; a < nof - 1; a++) {
1282 for (Int_t b = a + 1; b < nof; b++) {
1283 if (fMixedTest_STSonly_eventno[a] == fMixedTest_STSonly_eventno[b]) continue;
1284
1285 TVector3 e11 = fMixedTest_STSonly_photons[a][0];
1286 TVector3 e12 = fMixedTest_STSonly_photons[a][1];
1287 TVector3 e21 = fMixedTest_STSonly_photons[b][0];
1288 TVector3 e22 = fMixedTest_STSonly_photons[b][1];
1289
1290 Bool_t IsRich11 = fMixedTest_STSonly_hasRichInd[a][0];
1291 Bool_t IsRich12 = fMixedTest_STSonly_hasRichInd[a][1];
1292 Bool_t IsRich21 = fMixedTest_STSonly_hasRichInd[b][0];
1293 Bool_t IsRich22 = fMixedTest_STSonly_hasRichInd[b][1];
1294 Int_t IsRichSum = IsRich11 + IsRich12 + IsRich21 + IsRich22;
1295
1296
1297 CbmAnaConversionKinematicParams params =
1298 CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco(e11, e12, e21, e22);
1299
1300 if (IsRichSum == 2) { fhTest_eventMixing_STSonly_2p2->Fill(params.fMinv); }
1301 if (IsRichSum == 3) { fhTest_eventMixing_STSonly_3p1->Fill(params.fMinv); }
1302 if (IsRichSum == 4) { fhTest_eventMixing_STSonly_4p0->Fill(params.fMinv); }
1303 }
1304 }
1305}
1306
1307
1308Double_t CbmAnaConversionTest::ElectronANNvalue(Int_t globalTrackIndex, Double_t momentum)
1309{
1310 if (NULL == fGlobalTracks || NULL == fRichRings) return -2;
1311 //CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(globalTrackIndex);
1312 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(globalTrackIndex));
1313 Int_t richId = globalTrack->GetRichRingIndex();
1314 if (richId < 0) return -2;
1315 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(richId));
1316 if (NULL == ring) return -2;
1317
1318 Double_t ann = CbmRichElectronIdAnn::GetInstance().CalculateAnnValue(globalTrackIndex, momentum);
1319 return ann;
1320}
M2E
#define M2E
Definition CbmAnaJpsiKinematicParams.h:14
CbmStsTrack.h
Data class for STS tracks.
Int_t
int Int_t
Definition RootTypesDef.h:16
Bool_t
bool Bool_t
Definition RootTypesDef.h:15
CbmAnaConversionCutSettings::CalcChiCut
static Double_t CalcChiCut(Double_t pt)
Definition CbmAnaConversionCutSettings.h:80
CbmAnaConversionGlobalFunctions::OpeningAngleBetweenGamma
static Double_t OpeningAngleBetweenGamma(const TVector3 part11, const TVector3 part12, const TVector3 part21, const TVector3 part22)
Definition CbmAnaConversionGlobalFunctions.h:80
CbmAnaConversionKinematicParams::KinematicParams_2particles_Reco
static CbmAnaConversionKinematicParams KinematicParams_2particles_Reco(const TVector3 electron1, const TVector3 electron2)
Definition CbmAnaConversionKinematicParams.h:128
CbmAnaConversionKinematicParams::KinematicParams_4particles_Reco
static CbmAnaConversionKinematicParams KinematicParams_4particles_Reco(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
Definition CbmAnaConversionKinematicParams.h:87
CbmAnaConversionTest::fVector_electronRICH_reconstructedPhotons
std::vector< std::vector< int > > fVector_electronRICH_reconstructedPhotons
Definition CbmAnaConversionTest.h:137
CbmAnaConversionTest::fVector_electronRICH_mcIndex
std::vector< int > fVector_electronRICH_mcIndex
Definition CbmAnaConversionTest.h:135
CbmAnaConversionTest::fMixedTest_3p1_eventno
std::vector< int > fMixedTest_3p1_eventno
Definition CbmAnaConversionTest.h:160
CbmAnaConversionTest::CalculateKinematicParamsReco
LmvmKinePar CalculateKinematicParamsReco(const TVector3 electron1, const TVector3 electron2)
Definition CbmAnaConversionTest.cxx:973
CbmAnaConversionTest::fMixedTest_3p1_combined
std::vector< int > fMixedTest_3p1_combined
Definition CbmAnaConversionTest.h:161
CbmAnaConversionTest::fVector_reconstructedPhotons_FromSTSandRICH
std::vector< std::vector< int > > fVector_reconstructedPhotons_FromSTSandRICH
Definition CbmAnaConversionTest.h:130
CbmAnaConversionTest::fVector_chi
std::vector< double > fVector_chi
Definition CbmAnaConversionTest.h:125
CbmAnaConversionTest::fVector_electronRICH_gtIndex
std::vector< int > fVector_electronRICH_gtIndex
Definition CbmAnaConversionTest.h:134
CbmAnaConversionTest::fElectrons_gtid
std::vector< int > fElectrons_gtid
Definition CbmAnaConversionTest.h:80
CbmAnaConversionTest::fVector_richIndex
std::vector< int > fVector_richIndex
Definition CbmAnaConversionTest.h:127
CbmAnaConversionTest::fElectrons_pi0mcid
std::vector< int > fElectrons_pi0mcid
Definition CbmAnaConversionTest.h:83
CbmAnaConversionTest::fVector_electronRICH_gt
std::vector< CbmGlobalTrack * > fVector_electronRICH_gt
Definition CbmAnaConversionTest.h:133
CbmAnaConversionTest::fVector_AllMomenta
std::vector< TVector3 > fVector_AllMomenta
Definition CbmAnaConversionTest.h:119
CbmAnaConversionTest::fVector_reconstructedPhotons_STSonly
std::vector< std::vector< int > > fVector_reconstructedPhotons_STSonly
Definition CbmAnaConversionTest.h:140
CbmAnaConversionTest::fMixedTest_STSonly_hasRichInd
std::vector< std::vector< bool > > fMixedTest_STSonly_hasRichInd
Definition CbmAnaConversionTest.h:153
CbmAnaConversionTest::Invmass_4particlesRECO
Double_t Invmass_4particlesRECO(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
Definition CbmAnaConversionTest.cxx:1001
CbmAnaConversionTest::fElectrons_same
std::vector< int > fElectrons_same
Definition CbmAnaConversionTest.h:84
CbmAnaConversionTest::fMixedTest_STSonly_photons
std::vector< std::vector< TVector3 > > fMixedTest_STSonly_photons
Definition CbmAnaConversionTest.h:151
CbmAnaConversionTest::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition CbmAnaConversionTest.h:73
CbmAnaConversionTest::fVector_gtIndex
std::vector< int > fVector_gtIndex
Definition CbmAnaConversionTest.h:126
CbmAnaConversionTest::fMixedTest_STSonly_eventno
std::vector< int > fMixedTest_STSonly_eventno
Definition CbmAnaConversionTest.h:152
CbmAnaConversionTest::fMixedTest_3p1_photons
std::vector< std::vector< TVector3 > > fMixedTest_3p1_photons
Definition CbmAnaConversionTest.h:159
CbmAnaConversionTest::ElectronANNvalue
Double_t ElectronANNvalue(Int_t globalTrackIndex, Double_t momentum)
Definition CbmAnaConversionTest.cxx:1308
CbmAnaConversionTest::fVector_gt
std::vector< CbmGlobalTrack * > fVector_gt
Definition CbmAnaConversionTest.h:123
CbmAnaConversionTest::fVector_withRichSignal
std::vector< bool > fVector_withRichSignal
Definition CbmAnaConversionTest.h:129
CbmAnaConversionTest::fVector_mcIndex
std::vector< int > fVector_mcIndex
Definition CbmAnaConversionTest.h:128
CbmAnaConversionTest::HasRichInd
Bool_t HasRichInd(Int_t gtIndex, Int_t arrayIndex)
Definition CbmAnaConversionTest.cxx:1158
CbmAnaConversionTest::fHistoList_test
std::vector< TH1 * > fHistoList_test
Definition CbmAnaConversionTest.h:78
CbmAnaConversionTest::fMixedTest_3p1_ann
std::vector< std::vector< Double_t > > fMixedTest_3p1_ann
Definition CbmAnaConversionTest.h:162
CbmAnaConversionTest::fElectrons_richInd
std::vector< int > fElectrons_richInd
Definition CbmAnaConversionTest.h:82
CbmAnaConversionTest::fElectrons_mcid
std::vector< int > fElectrons_mcid
Definition CbmAnaConversionTest.h:81
CbmAnaConversionTest::CalcInvMass
Double_t CalcInvMass(Int_t e1, Int_t e2, Int_t e3, Int_t e4)
Definition CbmAnaConversionTest.cxx:421
CbmAnaConversionTest::fVector_momenta
std::vector< TVector3 > fVector_momenta
Definition CbmAnaConversionTest.h:124
CbmAnaConversionTest::fRichRingMatches
TClonesArray * fRichRingMatches
Definition CbmAnaConversionTest.h:70
CbmAnaConversionTest::fVector_electronRICH_momenta
std::vector< TVector3 > fVector_electronRICH_momenta
Definition CbmAnaConversionTest.h:136
CbmGlobalTrack::GetStsTrackIndex
int32_t GetStsTrackIndex() const
Definition CbmGlobalTrack.h:44
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition CbmL1PFFitter.cxx:501
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, const std::vector< CbmMvdHit > &vMvdHits, const std::vector< CbmStsHit > &vStsHits, const std::vector< int > &pidHypo)
Definition CbmL1PFFitter.cxx:172
CbmLitGlobalElectronId::IsRichElectron
Bool_t IsRichElectron(Int_t globalTrackIndex, Double_t momentum)
Identify electron in RICH detector.
Definition CbmLitGlobalElectronId.cxx:64
CbmLitGlobalElectronId::SetRichAnnCut
void SetRichAnnCut(Double_t par)
Set cut on RICH ANN output value.
Definition CbmLitGlobalElectronId.h:115
CbmLitGlobalElectronId::GetInstance
static CbmLitGlobalElectronId & GetInstance()
Definition CbmLitGlobalElectronId.h:30
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMCTrack::Get4Momentum
void Get4Momentum(TLorentzVector &momentum) const
Definition CbmMCTrack.h:173
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
CbmRichElectronIdAnn::CalculateAnnValue
double CalculateAnnValue(int globalTrackIndex, double momentum)
Calculate output value of the ANN.
Definition CbmRichElectronIdAnn.cxx:85
CbmRichElectronIdAnn::GetInstance
static CbmRichElectronIdAnn & GetInstance()
Definition CbmRichElectronIdAnn.h:47
LmvmKinePar::fAngle
Double_t fAngle
Definition LmvmKinePar.h:23
LmvmKinePar::fMinv
Double_t fMinv
Definition LmvmKinePar.h:22
LmvmKinePar::fPt
Double_t fPt
Definition LmvmKinePar.h:20
LmvmKinePar::fMomentumMag
Double_t fMomentumMag
Definition LmvmKinePar.h:19
LmvmKinePar::fRapidity
Double_t fRapidity
Definition LmvmKinePar.h:21
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21