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