CbmRoot
CbmAnaConversionTest2.cxx
Go to the documentation of this file.
1 
9 
12 #include "CbmGlobalTrack.h"
13 #include "CbmL1PFFitter.h"
14 #include "CbmLitGlobalElectronId.h"
15 #include "CbmMCTrack.h"
16 #include "CbmRichPoint.h"
17 #include "CbmRichRing.h"
18 #include "CbmStsTrack.h"
19 #include "CbmTrackMatchNew.h"
20 #include "FairLogger.h"
21 #include "FairRootManager.h"
22 #include "L1Field.h"
23 
24 #include <algorithm>
25 #include <map>
26 
27 
28 using namespace std;
29 
30 
32  : fRichPoints(NULL)
33  , fRichRings(NULL)
34  , fRichRingMatches(NULL)
35  , fMcTracks(NULL)
36  , fStsTracks(NULL)
37  , fStsTrackMatches(NULL)
38  , fGlobalTracks(NULL)
39  , fPrimVertex(NULL)
40  , fKFVertex()
41  , fHistoList_test2()
42  , fhTest2_invmass_RICHindex0(NULL)
43  , fhTest2_invmass_RICHindex1(NULL)
44  , fhTest2_invmass_RICHindex2(NULL)
45  , fhTest2_invmass_RICHindex3(NULL)
46  , fhTest2_invmass_RICHindex4(NULL)
47  , fVector_gt()
48  , fVector_momenta()
49  , fVector_mctrack()
50  , fVector_gtIndex()
51  , fVector_richIndex()
52  , fhTest2_invmass_gee_mc(NULL)
53  , fhTest2_invmass_gee_refitted(NULL)
54  , fhTest2_invmass_gg_mc(NULL)
55  , fhTest2_invmass_gg_refitted(NULL)
56  , fhTest2_invmass_all_mc(NULL)
57  , fhTest2_invmass_all_refitted(NULL)
58  , fhTest2_pt_vs_rap_gee(NULL)
59  , fhTest2_pt_vs_rap_gg(NULL)
60  , fhTest2_pt_vs_rap_all(NULL)
61  , fhTest2_startvertexElectrons_gee(NULL)
62  , fhTest2_startvertexElectrons_gg(NULL)
63  , fhTest2_startvertexElectrons_all(NULL)
64  , fhTest2_2rich_invmass_gee_mc(NULL)
65  , fhTest2_2rich_invmass_gee_refitted(NULL)
66  , fhTest2_2rich_invmass_gg_mc(NULL)
67  , fhTest2_2rich_invmass_gg_refitted(NULL)
68  , fhTest2_2rich_invmass_all_mc(NULL)
69  , fhTest2_2rich_invmass_all_refitted(NULL)
70  , fhTest2_2rich_pt_vs_rap_gee(NULL)
71  , fhTest2_2rich_pt_vs_rap_gg(NULL)
72  , fhTest2_2rich_pt_vs_rap_all(NULL)
73  , fhTest2_electrons_pt_vs_p(NULL)
74  , fhTest2_electrons_pt_vs_p_withRICH(NULL)
75  , fhTest2_electrons_pt_vs_p_noRICH(NULL)
76  , fhTest2_3rich_electrons_theta_included(NULL)
77  , fhTest2_3rich_electrons_theta_missing(NULL)
78  , fhTest2_3rich_electrons_thetaVSp_included(NULL)
79  , fhTest2_3rich_electrons_thetaVSp_missing(NULL) {}
80 
82 
83 
85  FairRootManager* ioman = FairRootManager::Instance();
86  if (NULL == ioman) {
87  Fatal("CbmAnaConversion::Init", "RootManager not instantised!");
88  }
89 
90  fRichPoints = (TClonesArray*) ioman->GetObject("RichPoint");
91  if (NULL == fRichPoints) {
92  Fatal("CbmAnaConversion::Init", "No RichPoint array!");
93  }
94 
95  fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
96  if (NULL == fMcTracks) {
97  Fatal("CbmAnaConversion::Init", "No MCTrack array!");
98  }
99 
100  fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
101  if (NULL == fStsTracks) {
102  Fatal("CbmAnaConversion::Init", "No StsTrack array!");
103  }
104 
105  fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
106  if (NULL == fStsTrackMatches) {
107  Fatal("CbmAnaConversion::Init", "No StsTrackMatch array!");
108  }
109 
110  fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
111  if (NULL == fGlobalTracks) {
112  Fatal("CbmAnaConversion::Init", "No GlobalTrack array!");
113  }
114 
115  // Get pointer to PrimaryVertex object from IOManager if it exists
116  // The old name for the object is "PrimaryVertex" the new one
117  // "PrimaryVertex." Check first for the new name
118  fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex."));
119  if (nullptr == fPrimVertex) {
120  fPrimVertex = dynamic_cast<CbmVertex*>(ioman->GetObject("PrimaryVertex"));
121  }
122  if (nullptr == fPrimVertex) { LOG(fatal) << "No PrimaryVertex array!"; }
123 
124  fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
125  if (NULL == fRichRings) {
126  Fatal("CbmAnaConversion::Init", "No RichRing array!");
127  }
128 
129  fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
130  if (NULL == fRichRingMatches) {
131  Fatal("CbmAnaConversion::Init", "No RichRingMatch array!");
132  }
133 
134  InitHistos();
135 }
136 
137 
139  fHistoList_test2.clear();
140 
141  Double_t invmassSpectra_nof = 800;
142  Double_t invmassSpectra_start = -0.00125;
143  Double_t invmassSpectra_end = 1.99875;
144 
145 
146  fhTest2_invmass_RICHindex0 = new TH1D(
147  "fhTest2_invmass_RICHindex0",
148  "fhTest2_invmass_RICHindex0; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
149  600,
150  -0.0025,
151  2.9975);
153  fhTest2_invmass_RICHindex1 = new TH1D(
154  "fhTest2_invmass_RICHindex1",
155  "fhTest2_invmass_RICHindex1; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
156  600,
157  -0.0025,
158  2.9975);
160  fhTest2_invmass_RICHindex2 = new TH1D(
161  "fhTest2_invmass_RICHindex2",
162  "fhTest2_invmass_RICHindex2; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
163  600,
164  -0.0025,
165  2.9975);
167  fhTest2_invmass_RICHindex3 = new TH1D(
168  "fhTest2_invmass_RICHindex3",
169  "fhTest2_invmass_RICHindex3; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
170  600,
171  -0.0025,
172  2.9975);
174  fhTest2_invmass_RICHindex4 = new TH1D(
175  "fhTest2_invmass_RICHindex4",
176  "fhTest2_invmass_RICHindex4; invariant mass of 4 e^{#pm} in GeV/c^{2}; #",
177  600,
178  -0.0025,
179  2.9975);
181 
182 
183  fhTest2_invmass_gee_mc = new TH1D(
184  "fhTest2_invmass_gee_mc",
185  "fhTest2_invmass_gee_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
186  invmassSpectra_nof,
187  invmassSpectra_start,
188  invmassSpectra_end);
189  fhTest2_invmass_gee_refitted = new TH1D(
190  "fhTest2_invmass_gee_refitted",
191  "fhTest2_invmass_gee_refitted;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
192  invmassSpectra_nof,
193  invmassSpectra_start,
194  invmassSpectra_end);
196  new TH1D("fhTest2_invmass_gg_mc",
197  "fhTest2_invmass_gg_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
198  invmassSpectra_nof,
199  invmassSpectra_start,
200  invmassSpectra_end);
201  fhTest2_invmass_gg_refitted = new TH1D(
202  "fhTest2_invmass_gg_refitted",
203  "fhTest2_invmass_gg_refitted;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
204  invmassSpectra_nof,
205  invmassSpectra_start,
206  invmassSpectra_end);
207  fhTest2_invmass_all_mc = new TH1D(
208  "fhTest2_invmass_all_mc",
209  "fhTest2_invmass_all_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
210  invmassSpectra_nof,
211  invmassSpectra_start,
212  invmassSpectra_end);
213  fhTest2_invmass_all_refitted = new TH1D(
214  "fhTest2_invmass_all_refitted",
215  "fhTest2_invmass_all_refitted;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
216  invmassSpectra_nof,
217  invmassSpectra_start,
218  invmassSpectra_end);
225 
227  new TH2D("fhTest2_pt_vs_rap_gee",
228  "fhTest2_pt_vs_rap_gee;p_{t} in GeV/c; rapidity y",
229  240,
230  -2.,
231  10.,
232  270,
233  -2.,
234  7.);
236  new TH2D("fhTest2_pt_vs_rap_gg",
237  "fhTest2_pt_vs_rap_gg;p_{t} in GeV/c; rapidity y",
238  240,
239  -2.,
240  10.,
241  270,
242  -2.,
243  7.);
245  new TH2D("fhTest2_pt_vs_rap_all",
246  "fhTest2_pt_vs_rap_all;p_{t} in GeV/c; rapidity y",
247  240,
248  -2.,
249  10.,
250  270,
251  -2.,
252  7.);
256 
258  new TH1D("fhTest2_startvertexElectrons_gee",
259  "fhTest2_startvertexElectrons_gee;z in cm;#",
260  411,
261  -5.25,
262  200.25);
264  new TH1D("fhTest2_startvertexElectrons_gg",
265  "fhTest2_startvertexElectrons_gg;z in cm;#",
266  411,
267  -5.25,
268  200.25);
270  new TH1D("fhTest2_startvertexElectrons_all",
271  "fhTest2_startvertexElectrons_all;z in cm;#",
272  411,
273  -5.25,
274  200.25);
278 
279 
280  // 2 leptons in RICH
281  fhTest2_2rich_invmass_gee_mc = new TH1D(
282  "fhTest2_2rich_invmass_gee_mc",
283  "fhTest2_2rich_invmass_gee_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
284  invmassSpectra_nof,
285  invmassSpectra_start,
286  invmassSpectra_end);
288  new TH1D("fhTest2_2rich_invmass_gee_refitted",
289  "fhTest2_2rich_invmass_gee_refitted;invariant mass of 4 e^{#pm} "
290  "in GeV/c^{2};#",
291  invmassSpectra_nof,
292  invmassSpectra_start,
293  invmassSpectra_end);
294  fhTest2_2rich_invmass_gg_mc = new TH1D(
295  "fhTest2_2rich_invmass_gg_mc",
296  "fhTest2_2rich_invmass_gg_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
297  invmassSpectra_nof,
298  invmassSpectra_start,
299  invmassSpectra_end);
301  new TH1D("fhTest2_2rich_invmass_gg_refitted",
302  "fhTest2_2rich_invmass_gg_refitted;invariant mass of 4 e^{#pm} in "
303  "GeV/c^{2};#",
304  invmassSpectra_nof,
305  invmassSpectra_start,
306  invmassSpectra_end);
307  fhTest2_2rich_invmass_all_mc = new TH1D(
308  "fhTest2_2rich_invmass_all_mc",
309  "fhTest2_2rich_invmass_all_mc;invariant mass of 4 e^{#pm} in GeV/c^{2};#",
310  invmassSpectra_nof,
311  invmassSpectra_start,
312  invmassSpectra_end);
314  new TH1D("fhTest2_2rich_invmass_all_refitted",
315  "fhTest2_2rich_invmass_all_refitted;invariant mass of 4 e^{#pm} "
316  "in GeV/c^{2};#",
317  invmassSpectra_nof,
318  invmassSpectra_start,
319  invmassSpectra_end);
326 
327 
329  new TH2D("fhTest2_2rich_pt_vs_rap_gee",
330  "fhTest2_2rich_pt_vs_rap_gee;p_{t} in GeV/c; rapidity y",
331  240,
332  -2.,
333  10.,
334  270,
335  -2.,
336  7.);
338  new TH2D("fhTest2_2rich_pt_vs_rap_gg",
339  "fhTest2_2rich_pt_vs_rap_gg;p_{t} in GeV/c; rapidity y",
340  240,
341  -2.,
342  10.,
343  270,
344  -2.,
345  7.);
347  new TH2D("fhTest2_2rich_pt_vs_rap_all",
348  "fhTest2_2rich_pt_vs_rap_all;p_{t} in GeV/c; rapidity y",
349  240,
350  -2.,
351  10.,
352  270,
353  -2.,
354  7.);
358 
359 
360  // further tests
362  new TH2D("fhTest2_electrons_pt_vs_p",
363  "fhTest2_electrons_pt_vs_p;p_{t} in GeV/c; p in GeV/c",
364  240,
365  -2.,
366  10.,
367  360,
368  -2.,
369  16.);
372  new TH2D("fhTest2_electrons_pt_vs_p_withRICH",
373  "fhTest2_electrons_pt_vs_p_withRICH;p_{t} in GeV/c; p in GeV/c",
374  240,
375  -2.,
376  10.,
377  360,
378  -2.,
379  16.);
382  new TH2D("fhTest2_electrons_pt_vs_p_noRICH",
383  "fhTest2_electrons_pt_vs_p_noRICH;p_{t} in GeV/c; p in GeV/c",
384  240,
385  -2.,
386  10.,
387  360,
388  -2.,
389  16.);
391 
393  new TH1D("fhTest2_3rich_electrons_theta_included",
394  "fhTest2_3rich_electrons_theta_included;theta angle [deg];#",
395  90,
396  0,
397  90);
400  new TH1D("fhTest2_3rich_electrons_theta_missing",
401  "fhTest2_3rich_electrons_theta_missing;theta angle [deg];#",
402  90,
403  0,
404  90);
407  "fhTest2_3rich_electrons_thetaVSp_included",
408  "fhTest2_3rich_electrons_thetaVSp_included;theta angle [deg];p in GeV/c",
409  90,
410  0,
411  90,
412  540,
413  -2.,
414  16.);
417  "fhTest2_3rich_electrons_thetaVSp_missing",
418  "fhTest2_3rich_electrons_thetaVSp_missing;theta angle [deg];p in GeV/c",
419  90,
420  0,
421  90,
422  540,
423  -2.,
424  16.);
426 }
427 
428 
430  //gDirectory->cd("analysis-conversion");
431  gDirectory->mkdir("Test2");
432  gDirectory->cd("Test2");
433  for (UInt_t i = 0; i < fHistoList_test2.size(); i++) {
434  fHistoList_test2[i]->Write();
435  }
436  gDirectory->cd("..");
437 }
438 
439 
441  fVector_gt.clear();
442  fVector_momenta.clear();
443  fVector_mctrack.clear();
444  fVector_gtIndex.clear();
445  fVector_richIndex.clear();
446 
447 
448  Int_t ngTracks = fGlobalTracks->GetEntriesFast();
449  for (Int_t i = 0; i < ngTracks; i++) {
450  CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
451  if (NULL == gTrack) continue;
452  int stsInd = gTrack->GetStsTrackIndex();
453  int richInd = gTrack->GetRichRingIndex();
454 
455  if (stsInd < 0) continue;
456  CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
457  if (stsTrack == NULL) continue;
458  CbmTrackMatchNew* stsMatch =
459  (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
460  if (stsMatch == NULL) continue;
461  if (stsMatch->GetNofLinks() <= 0) continue;
462  int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
463  if (stsMcTrackId < 0) continue;
464  CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
465  if (mcTrack1 == NULL) continue;
466 
467 
468  // calculate refitted momenta at primary vertex
469  TVector3 refittedMomentum;
470  CbmL1PFFitter fPFFitter;
471  vector<CbmStsTrack> stsTracks;
472  stsTracks.resize(1);
473  stsTracks[0] = *stsTrack;
474  vector<L1FieldRegion> vField;
475  vector<float> chiPrim;
476  fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, fKFVertex, 3e6);
477  //cand.chi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
478  //cand.chi2Prim = chiPrim[0];
479  const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();
480  vtxTrack->Momentum(refittedMomentum);
481  //float result_chi = chiPrim[0];
482 
483 
484  // Doing refit of momenta with electron assumption
485  CbmL1PFFitter fPFFitter_electron;
486  vector<CbmStsTrack> stsTracks_electron;
487  stsTracks_electron.resize(1);
488  stsTracks_electron[0] = *stsTrack;
489  vector<L1FieldRegion> vField_electron;
490  vector<float> chiPrim_electron;
491  vector<int> pidHypo_electron;
492  pidHypo_electron.push_back(11);
493  fPFFitter_electron.Fit(stsTracks_electron, pidHypo_electron);
494  fPFFitter_electron.GetChiToVertex(
495  stsTracks_electron, vField_electron, chiPrim_electron, fKFVertex, 3e6);
496 
497  TVector3 refittedMomentum_electron;
498  const FairTrackParam* vtxTrack_electron =
499  stsTracks_electron[0].GetParamFirst();
500  vtxTrack_electron->Momentum(refittedMomentum_electron);
501  //float result_chi_electron = chiPrim_electron[0];
502  //float result_ndf_electron = stsTracks_electron[0].GetNDF();
503  //Double_t stsHits = stsTrack->GetNofHits();
504 
505 
506  Int_t pdg = mcTrack1->GetPdgCode();
507 
508  if (TMath::Abs(pdg) == 11) {
509  fhTest2_electrons_pt_vs_p->Fill(refittedMomentum_electron.Perp(),
510  refittedMomentum_electron.Mag());
511  if (richInd > 0) {
513  refittedMomentum_electron.Perp(), refittedMomentum_electron.Mag());
514  }
515  if (richInd < 0) {
516  fhTest2_electrons_pt_vs_p_noRICH->Fill(refittedMomentum_electron.Perp(),
517  refittedMomentum_electron.Mag());
518  }
519 
520 
521  fVector_momenta.push_back(refittedMomentum_electron);
522  fVector_mctrack.push_back(mcTrack1);
523  fVector_gtIndex.push_back(i);
524  fVector_gt.push_back(gTrack);
525 
526  if (richInd < 0) {
527  fVector_richIndex.push_back(0);
528  continue;
529  }
530  CbmTrackMatchNew* richMatch =
531  (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
532  if (richMatch == NULL) {
533  fVector_richIndex.push_back(0);
534  continue;
535  }
536  int richMcTrackId = richMatch->GetMatchedLink().GetIndex();
537  if (richMcTrackId < 0) {
538  fVector_richIndex.push_back(0);
539  continue;
540  }
541  CbmMCTrack* mcTrack2 = (CbmMCTrack*) fMcTracks->At(richMcTrackId);
542  if (mcTrack2 == NULL) {
543  fVector_richIndex.push_back(0);
544  continue;
545  }
546 
547  fVector_richIndex.push_back(richInd);
548  }
549  }
550 
553 }
554 
555 
557 // Calculating invariant mass of 4 ep/em, using MC data AND reconstructed momentum
558 {
559  cout << "CbmAnaConversionTest2: InvariantMassTest_3RICH - Start..." << endl;
560  cout << "CbmAnaConversionTest2: InvariantMassTest_3RICH - "
561  << fVector_mctrack.size() << endl;
562  int fill = 0;
563  if (fVector_mctrack.size() < 4) return;
564  for (unsigned int i = 0; i < fVector_mctrack.size() - 3; i++) {
565  if (i % 10 == 0)
566  cout << "CbmAnaConversionTest2: InvariantMassTest_3RICH - iteration i = "
567  << i << endl;
568  for (unsigned int j = i + 1; j < fVector_mctrack.size() - 2; j++) {
569  for (unsigned int k = j + 1; k < fVector_mctrack.size() - 1; k++) {
570  for (unsigned int l = k + 1; l < fVector_mctrack.size(); l++) {
571  if (fVector_mctrack[i]->GetPdgCode()
572  + fVector_mctrack[j]->GetPdgCode()
573  + fVector_mctrack[k]->GetPdgCode()
574  + fVector_mctrack[l]->GetPdgCode()
575  != 0)
576  continue;
577 
578  if (fVector_mctrack.size() != fVector_richIndex.size()) {
579  cout << "CbmAnaConversionTest2: InvariantMassTest_4epem - not "
580  "matching number of entries!"
581  << endl;
582  continue;
583  }
584 
585  Bool_t IsPi0 = false;
586 
587  // starting points of each electron (-> i.e. conversion points of gamma OR decay points of pi0, depending on decay channel)
588  TVector3 pi0start_i;
589  fVector_mctrack[i]->GetStartVertex(pi0start_i);
590  TVector3 pi0start_j;
591  fVector_mctrack[j]->GetStartVertex(pi0start_j);
592  TVector3 pi0start_k;
593  fVector_mctrack[k]->GetStartVertex(pi0start_k);
594  TVector3 pi0start_l;
595  fVector_mctrack[l]->GetStartVertex(pi0start_l);
596 
597 
598  Int_t richIndex1 = (fVector_richIndex[i] > 0);
599  Int_t richIndex2 = (fVector_richIndex[j] > 0);
600  Int_t richIndex3 = (fVector_richIndex[k] > 0);
601  Int_t richIndex4 = (fVector_richIndex[l] > 0);
602 
603  if (richIndex1 + richIndex2 + richIndex3 + richIndex4 != 3) continue;
604 
605 
606  int motherId1 = fVector_mctrack[i]->GetMotherId();
607  int motherId2 = fVector_mctrack[j]->GetMotherId();
608  int motherId3 = fVector_mctrack[k]->GetMotherId();
609  int motherId4 = fVector_mctrack[l]->GetMotherId();
610 
611 
612  if ((motherId1 == motherId2 && motherId3 == motherId4)
613  || (motherId1 == motherId3 && motherId2 == motherId4)
614  || (motherId1 == motherId4
615  && motherId2 == motherId3)) { // start 1
616 
617 
618  int grandmotherId1 = -1;
619  int grandmotherId2 = -1;
620  int grandmotherId3 = -1;
621  int grandmotherId4 = -1;
622 
623  int mcMotherPdg1 = -1;
624  // int mcMotherPdg2 = -1; (FU) not used
625  // int mcMotherPdg3 = -1; (FU) not used
626  // int mcMotherPdg4 = -1; (FU) not used
627  int mcGrandmotherPdg1 = -1;
628  // int mcGrandmotherPdg2 = -1;
629  // int mcGrandmotherPdg3 = -1;
630  // int mcGrandmotherPdg4 = -1;
631 
632  CbmMCTrack* grandmother1 = nullptr;
633 
634  if (motherId1 != -1) {
635  CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
636  if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
637  grandmotherId1 = mother1->GetMotherId();
638  if (grandmotherId1 != -1) {
639  grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
640  if (NULL != grandmother1)
641  mcGrandmotherPdg1 = grandmother1->GetPdgCode();
642  }
643  }
644  if (motherId2 != -1) {
645  CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
646  // if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode(); (FU) not used
647  grandmotherId2 = mother2->GetMotherId();
648  if (grandmotherId2 != -1) {
649  // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
650  // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
651  }
652  }
653  if (motherId3 != -1) {
654  CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
655  // if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode(); (FU) not used
656  grandmotherId3 = mother3->GetMotherId();
657  if (grandmotherId3 != -1) {
658  // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
659  // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
660  }
661  }
662  if (motherId4 != -1) {
663  CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
664  // if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode(); (FU) not used
665  grandmotherId4 = mother4->GetMotherId();
666  if (grandmotherId4 != -1) {
667  // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
668  // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
669  }
670  }
671 
672  if (grandmotherId1 == -1) continue;
673 
674  if (motherId1 == motherId2 && motherId3 == motherId4) {
675  if (NofDaughters(motherId1) != 2 || NofDaughters(motherId3) != 2)
676  continue;
677  if ((grandmotherId1 == motherId3 && mcGrandmotherPdg1 == 111)
678  || (motherId1 == grandmotherId3 && mcMotherPdg1 == 111)) {
679 
680  fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
681  fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
682  fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
683  fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
684  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
685  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
686  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
687  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
688  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
689  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1
690  || pi0start_k.Z() > 1 || pi0start_l.Z() > 1)
691  continue;
692 
693 
694  Double_t invmass_mc = 0;
695  Double_t invmass_reco = 0;
699  fVector_mctrack[j],
700  fVector_mctrack[k],
701  fVector_mctrack[l]);
702  CbmAnaConversionKinematicParams params_reco =
705  fVector_momenta[j],
706  fVector_momenta[k],
707  fVector_momenta[l]);
708  invmass_mc = params_mc.fMinv;
709  invmass_reco = params_reco.fMinv;
710 
711  fhTest2_invmass_gee_mc->Fill(invmass_mc);
712  fhTest2_invmass_gee_refitted->Fill(invmass_reco);
713  fhTest2_invmass_all_mc->Fill(invmass_mc);
714  fhTest2_invmass_all_refitted->Fill(invmass_reco);
715 
716  fhTest2_pt_vs_rap_gee->Fill(params_reco.fPt,
717  params_reco.fRapidity);
718  fhTest2_pt_vs_rap_all->Fill(params_reco.fPt,
719  params_reco.fRapidity);
720 
721  IsPi0 = true;
722  }
723  }
724  if (motherId1 == motherId3 && motherId2 == motherId4) {
725  if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2)
726  continue;
727  if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
728  || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
729 
730  fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
731  fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
732  fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
733  fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
734  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
735  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
736  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
737  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
738  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
739  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1
740  || pi0start_k.Z() > 1 || pi0start_l.Z() > 1)
741  continue;
742 
743 
744  Double_t invmass_mc = 0;
745  Double_t invmass_reco = 0;
749  fVector_mctrack[j],
750  fVector_mctrack[k],
751  fVector_mctrack[l]);
752  CbmAnaConversionKinematicParams params_reco =
755  fVector_momenta[j],
756  fVector_momenta[k],
757  fVector_momenta[l]);
758  invmass_mc = params_mc.fMinv;
759  invmass_reco = params_reco.fMinv;
760 
761  fhTest2_invmass_gee_mc->Fill(invmass_mc);
762  fhTest2_invmass_gee_refitted->Fill(invmass_reco);
763  fhTest2_invmass_all_mc->Fill(invmass_mc);
764  fhTest2_invmass_all_refitted->Fill(invmass_reco);
765 
766  fhTest2_pt_vs_rap_gee->Fill(params_reco.fPt,
767  params_reco.fRapidity);
768  fhTest2_pt_vs_rap_all->Fill(params_reco.fPt,
769  params_reco.fRapidity);
770 
771  IsPi0 = true;
772  }
773  }
774  if (motherId1 == motherId4 && motherId2 == motherId3) {
775  if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2)
776  continue;
777  if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
778  || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
779 
780  fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
781  fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
782  fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
783  fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
784  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
785  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
786  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
787  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
788  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
789  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1
790  || pi0start_k.Z() > 1 || pi0start_l.Z() > 1)
791  continue;
792 
793 
794  Double_t invmass_mc = 0;
795  Double_t invmass_reco = 0;
799  fVector_mctrack[j],
800  fVector_mctrack[k],
801  fVector_mctrack[l]);
802  CbmAnaConversionKinematicParams params_reco =
805  fVector_momenta[j],
806  fVector_momenta[k],
807  fVector_momenta[l]);
808  invmass_mc = params_mc.fMinv;
809  invmass_reco = params_reco.fMinv;
810 
811  fhTest2_invmass_gee_mc->Fill(invmass_mc);
812  fhTest2_invmass_gee_refitted->Fill(invmass_reco);
813  fhTest2_invmass_all_mc->Fill(invmass_mc);
814  fhTest2_invmass_all_refitted->Fill(invmass_reco);
815 
816  fhTest2_pt_vs_rap_gee->Fill(params_reco.fPt,
817  params_reco.fRapidity);
818  fhTest2_pt_vs_rap_all->Fill(params_reco.fPt,
819  params_reco.fRapidity);
820 
821  IsPi0 = true;
822  }
823  }
824 
825 
826  // ===================================================================================================
827  // HERE DECAY pi0 -> gamma gamma -> e+e- e+e-
828  if (grandmotherId1 == grandmotherId2
829  && grandmotherId1 == grandmotherId3
830  && grandmotherId1 == grandmotherId4) {
831  if (mcGrandmotherPdg1 != 111) continue; // 111 = pi0, 221 = eta
832 
833  TVector3 pi0start;
834  grandmother1->GetStartVertex(pi0start);
835 
836  fhTest2_startvertexElectrons_gg->Fill(pi0start_i.Z());
837  fhTest2_startvertexElectrons_gg->Fill(pi0start_j.Z());
838  fhTest2_startvertexElectrons_gg->Fill(pi0start_k.Z());
839  fhTest2_startvertexElectrons_gg->Fill(pi0start_l.Z());
840  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
841  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
842  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
843  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
844 
845 
846  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
847  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1
848  || pi0start_l.Z() > 1)
849  continue;
850 
851  Double_t invmass_mc = 0;
852  Double_t invmass_reco = 0;
856  fVector_mctrack[j],
857  fVector_mctrack[k],
858  fVector_mctrack[l]);
859  CbmAnaConversionKinematicParams params_reco =
862  fVector_momenta[j],
863  fVector_momenta[k],
864  fVector_momenta[l]);
865  invmass_mc = params_mc.fMinv;
866  invmass_reco = params_reco.fMinv;
867 
868  fhTest2_invmass_gg_mc->Fill(invmass_mc);
869  fhTest2_invmass_gg_refitted->Fill(invmass_reco);
870  fhTest2_invmass_all_mc->Fill(invmass_mc);
871  fhTest2_invmass_all_refitted->Fill(invmass_reco);
872 
873  fhTest2_pt_vs_rap_gg->Fill(params_reco.fPt,
874  params_reco.fRapidity);
875  fhTest2_pt_vs_rap_all->Fill(params_reco.fPt,
876  params_reco.fRapidity);
877 
878  IsPi0 = true;
879 
880  cout << "########################################################"
881  "##############"
882  << endl;
883 
884  fill++;
885  }
886 
887 
888  if (richIndex1 <= 0 && IsPi0) {
890  fVector_momenta[i].Theta() * 180 / TMath::Pi());
892  fVector_momenta[j].Theta() * 180 / TMath::Pi());
894  fVector_momenta[k].Theta() * 180 / TMath::Pi());
896  fVector_momenta[l].Theta() * 180 / TMath::Pi());
898  fVector_momenta[i].Theta() * 180 / TMath::Pi(),
899  fVector_momenta[i].Mag());
901  fVector_momenta[j].Theta() * 180 / TMath::Pi(),
902  fVector_momenta[j].Mag());
904  fVector_momenta[k].Theta() * 180 / TMath::Pi(),
905  fVector_momenta[k].Mag());
907  fVector_momenta[l].Theta() * 180 / TMath::Pi(),
908  fVector_momenta[l].Mag());
909  }
910  if (richIndex2 <= 0 && IsPi0) {
912  fVector_momenta[j].Theta() * 180 / TMath::Pi());
914  fVector_momenta[i].Theta() * 180 / TMath::Pi());
916  fVector_momenta[k].Theta() * 180 / TMath::Pi());
918  fVector_momenta[l].Theta() * 180 / TMath::Pi());
920  fVector_momenta[j].Theta() * 180 / TMath::Pi(),
921  fVector_momenta[j].Mag());
923  fVector_momenta[i].Theta() * 180 / TMath::Pi(),
924  fVector_momenta[i].Mag());
926  fVector_momenta[k].Theta() * 180 / TMath::Pi(),
927  fVector_momenta[k].Mag());
929  fVector_momenta[l].Theta() * 180 / TMath::Pi(),
930  fVector_momenta[l].Mag());
931  }
932  if (richIndex3 <= 0 && IsPi0) {
934  fVector_momenta[k].Theta() * 180 / TMath::Pi());
936  fVector_momenta[i].Theta() * 180 / TMath::Pi());
938  fVector_momenta[j].Theta() * 180 / TMath::Pi());
940  fVector_momenta[l].Theta() * 180 / TMath::Pi());
942  fVector_momenta[k].Theta() * 180 / TMath::Pi(),
943  fVector_momenta[k].Mag());
945  fVector_momenta[i].Theta() * 180 / TMath::Pi(),
946  fVector_momenta[i].Mag());
948  fVector_momenta[j].Theta() * 180 / TMath::Pi(),
949  fVector_momenta[j].Mag());
951  fVector_momenta[l].Theta() * 180 / TMath::Pi(),
952  fVector_momenta[l].Mag());
953  }
954  if (richIndex4 <= 0 && IsPi0) {
956  fVector_momenta[l].Theta() * 180 / TMath::Pi());
958  fVector_momenta[i].Theta() * 180 / TMath::Pi());
960  fVector_momenta[j].Theta() * 180 / TMath::Pi());
962  fVector_momenta[k].Theta() * 180 / TMath::Pi());
964  fVector_momenta[l].Theta() * 180 / TMath::Pi(),
965  fVector_momenta[l].Mag());
967  fVector_momenta[i].Theta() * 180 / TMath::Pi(),
968  fVector_momenta[i].Mag());
970  fVector_momenta[j].Theta() * 180 / TMath::Pi(),
971  fVector_momenta[j].Mag());
973  fVector_momenta[k].Theta() * 180 / TMath::Pi(),
974  fVector_momenta[k].Mag());
975  }
976 
977  } // end 1
978  }
979  }
980  }
981  }
982 }
983 
984 
986 // Calculating invariant mass of 4 ep/em, using MC data AND reconstructed momentum
987 {
988  cout << "CbmAnaConversionTest2: InvariantMassTest_2RICH - Start..." << endl;
989  cout << "CbmAnaConversionTest2: InvariantMassTest_2RICH - "
990  << fVector_mctrack.size() << endl;
991  int fill = 0;
992  if (fVector_mctrack.size() < 4) return;
993  for (unsigned int i = 0; i < fVector_mctrack.size() - 3; i++) {
994  if (i % 10 == 0)
995  cout << "CbmAnaConversionTest2: InvariantMassTest_2RICH - iteration i = "
996  << i << endl;
997  for (unsigned int j = i + 1; j < fVector_mctrack.size() - 2; j++) {
998  for (unsigned int k = j + 1; k < fVector_mctrack.size() - 1; k++) {
999  for (unsigned int l = k + 1; l < fVector_mctrack.size(); l++) {
1000  if (fVector_mctrack[i]->GetPdgCode()
1001  + fVector_mctrack[j]->GetPdgCode()
1002  + fVector_mctrack[k]->GetPdgCode()
1003  + fVector_mctrack[l]->GetPdgCode()
1004  != 0)
1005  continue;
1006 
1007  if (fVector_mctrack.size() != fVector_richIndex.size()) {
1008  cout << "CbmAnaConversionTest2: InvariantMassTest_4epem - not "
1009  "matching number of entries!"
1010  << endl;
1011  continue;
1012  }
1013 
1014 
1015  // starting points of each electron (-> i.e. conversion points of gamma OR decay points of pi0, depending on decay channel)
1016  TVector3 pi0start_i;
1017  fVector_mctrack[i]->GetStartVertex(pi0start_i);
1018  TVector3 pi0start_j;
1019  fVector_mctrack[j]->GetStartVertex(pi0start_j);
1020  TVector3 pi0start_k;
1021  fVector_mctrack[k]->GetStartVertex(pi0start_k);
1022  TVector3 pi0start_l;
1023  fVector_mctrack[l]->GetStartVertex(pi0start_l);
1024 
1025 
1026  Int_t richIndex1 = (fVector_richIndex[i] > 0);
1027  Int_t richIndex2 = (fVector_richIndex[j] > 0);
1028  Int_t richIndex3 = (fVector_richIndex[k] > 0);
1029  Int_t richIndex4 = (fVector_richIndex[l] > 0);
1030 
1031  if (richIndex1 + richIndex2 + richIndex3 + richIndex4 != 2) continue;
1032 
1033  int motherId1 = fVector_mctrack[i]->GetMotherId();
1034  int motherId2 = fVector_mctrack[j]->GetMotherId();
1035  int motherId3 = fVector_mctrack[k]->GetMotherId();
1036  int motherId4 = fVector_mctrack[l]->GetMotherId();
1037 
1038 
1039  if ((motherId1 == motherId2 && motherId3 == motherId4)
1040  || (motherId1 == motherId3 && motherId2 == motherId4)
1041  || (motherId1 == motherId4 && motherId2 == motherId3)) {
1042 
1043 
1044  int grandmotherId1 = -1;
1045  int grandmotherId2 = -1;
1046  int grandmotherId3 = -1;
1047  int grandmotherId4 = -1;
1048 
1049  int mcMotherPdg1 = -1;
1050  // int mcMotherPdg2 = -1; (FU) not used
1051  // int mcMotherPdg3 = -1; (FU) not used
1052  // int mcMotherPdg4 = -1; (FU) not used
1053  int mcGrandmotherPdg1 = -1;
1054  // int mcGrandmotherPdg2 = -1;
1055  // int mcGrandmotherPdg3 = -1;
1056  // int mcGrandmotherPdg4 = -1;
1057 
1058  CbmMCTrack* grandmother1 = nullptr;
1059 
1060  if (motherId1 != -1) {
1061  CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
1062  if (NULL != mother1) mcMotherPdg1 = mother1->GetPdgCode();
1063  grandmotherId1 = mother1->GetMotherId();
1064  if (grandmotherId1 != -1) {
1065  grandmother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
1066  if (NULL != grandmother1)
1067  mcGrandmotherPdg1 = grandmother1->GetPdgCode();
1068  }
1069  }
1070  if (motherId2 != -1) {
1071  CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
1072  // if (NULL != mother2) mcMotherPdg2 = mother2->GetPdgCode(); (FU) not used
1073  grandmotherId2 = mother2->GetMotherId();
1074  if (grandmotherId2 != -1) {
1075  // CbmMCTrack* grandmother2 = (CbmMCTrack*) fMcTracks->At(grandmotherId2);
1076  // if (NULL != grandmother2) mcGrandmotherPdg2 = grandmother2->GetPdgCode();
1077  }
1078  }
1079  if (motherId3 != -1) {
1080  CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
1081  // if (NULL != mother3) mcMotherPdg3 = mother3->GetPdgCode(); (FU) not used
1082  grandmotherId3 = mother3->GetMotherId();
1083  if (grandmotherId3 != -1) {
1084  // CbmMCTrack* grandmother3 = (CbmMCTrack*) fMcTracks->At(grandmotherId3);
1085  // if (NULL != grandmother3) mcGrandmotherPdg3 = grandmother3->GetPdgCode();
1086  }
1087  }
1088  if (motherId4 != -1) {
1089  CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
1090  // if (NULL != mother4) mcMotherPdg4 = mother4->GetPdgCode(); (FU) not used
1091  grandmotherId4 = mother4->GetMotherId();
1092  if (grandmotherId4 != -1) {
1093  // CbmMCTrack* grandmother4 = (CbmMCTrack*) fMcTracks->At(grandmotherId4);
1094  // if (NULL != grandmother4) mcGrandmotherPdg4 = grandmother4->GetPdgCode();
1095  }
1096  }
1097 
1098  if (grandmotherId1 == -1) continue;
1099 
1100  if (motherId1 == motherId2 && motherId3 == motherId4) {
1101  if (NofDaughters(motherId1) != 2 || NofDaughters(motherId3) != 2)
1102  continue;
1103  if ((grandmotherId1 == motherId3 && mcGrandmotherPdg1 == 111)
1104  || (motherId1 == grandmotherId3 && mcMotherPdg1 == 111)) {
1105  /*
1106  fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
1107  fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
1108  fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
1109  fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
1110  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
1111  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
1112  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
1113  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
1114  */
1115  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1116  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1
1117  || pi0start_k.Z() > 1 || pi0start_l.Z() > 1)
1118  continue;
1119 
1120 
1121  Double_t invmass_mc = 0;
1122  Double_t invmass_reco = 0;
1126  fVector_mctrack[j],
1127  fVector_mctrack[k],
1128  fVector_mctrack[l]);
1129  CbmAnaConversionKinematicParams params_reco =
1132  fVector_momenta[j],
1133  fVector_momenta[k],
1134  fVector_momenta[l]);
1135  invmass_mc = params_mc.fMinv;
1136  invmass_reco = params_reco.fMinv;
1137 
1138  fhTest2_2rich_invmass_gee_mc->Fill(invmass_mc);
1139  fhTest2_2rich_invmass_gee_refitted->Fill(invmass_reco);
1140  fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
1141  fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
1142 
1143  fhTest2_2rich_pt_vs_rap_gee->Fill(params_reco.fPt,
1144  params_reco.fRapidity);
1145  fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt,
1146  params_reco.fRapidity);
1147  }
1148  }
1149  if (motherId1 == motherId3 && motherId2 == motherId4) {
1150  if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2)
1151  continue;
1152  if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
1153  || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
1154  /*
1155  fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
1156  fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
1157  fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
1158  fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
1159  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
1160  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
1161  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
1162  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
1163  */
1164  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1165  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1
1166  || pi0start_k.Z() > 1 || pi0start_l.Z() > 1)
1167  continue;
1168 
1169 
1170  Double_t invmass_mc = 0;
1171  Double_t invmass_reco = 0;
1175  fVector_mctrack[j],
1176  fVector_mctrack[k],
1177  fVector_mctrack[l]);
1178  CbmAnaConversionKinematicParams params_reco =
1181  fVector_momenta[j],
1182  fVector_momenta[k],
1183  fVector_momenta[l]);
1184  invmass_mc = params_mc.fMinv;
1185  invmass_reco = params_reco.fMinv;
1186 
1187  fhTest2_2rich_invmass_gee_mc->Fill(invmass_mc);
1188  fhTest2_2rich_invmass_gee_refitted->Fill(invmass_reco);
1189  fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
1190  fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
1191 
1192  fhTest2_2rich_pt_vs_rap_gee->Fill(params_reco.fPt,
1193  params_reco.fRapidity);
1194  fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt,
1195  params_reco.fRapidity);
1196  }
1197  }
1198  if (motherId1 == motherId4 && motherId2 == motherId3) {
1199  if (NofDaughters(motherId1) != 2 || NofDaughters(motherId2) != 2)
1200  continue;
1201  if ((grandmotherId1 == motherId2 && mcGrandmotherPdg1 == 111)
1202  || (motherId1 == grandmotherId2 && mcMotherPdg1 == 111)) {
1203  /*
1204  fhTest2_startvertexElectrons_gee->Fill(pi0start_i.Z());
1205  fhTest2_startvertexElectrons_gee->Fill(pi0start_j.Z());
1206  fhTest2_startvertexElectrons_gee->Fill(pi0start_k.Z());
1207  fhTest2_startvertexElectrons_gee->Fill(pi0start_l.Z());
1208  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
1209  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
1210  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
1211  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
1212  */
1213  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1214  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1
1215  || pi0start_k.Z() > 1 || pi0start_l.Z() > 1)
1216  continue;
1217 
1218 
1219  Double_t invmass_mc = 0;
1220  Double_t invmass_reco = 0;
1224  fVector_mctrack[j],
1225  fVector_mctrack[k],
1226  fVector_mctrack[l]);
1227  CbmAnaConversionKinematicParams params_reco =
1230  fVector_momenta[j],
1231  fVector_momenta[k],
1232  fVector_momenta[l]);
1233  invmass_mc = params_mc.fMinv;
1234  invmass_reco = params_reco.fMinv;
1235 
1236  fhTest2_2rich_invmass_gee_mc->Fill(invmass_mc);
1237  fhTest2_2rich_invmass_gee_refitted->Fill(invmass_reco);
1238  fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
1239  fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
1240 
1241  fhTest2_2rich_pt_vs_rap_gee->Fill(params_reco.fPt,
1242  params_reco.fRapidity);
1243  fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt,
1244  params_reco.fRapidity);
1245  }
1246  }
1247 
1248 
1249  // ===================================================================================================
1250  // HERE DECAY pi0 -> gamma gamma -> e+e- e+e-
1251  if (grandmotherId1 == grandmotherId2
1252  && grandmotherId1 == grandmotherId3
1253  && grandmotherId1 == grandmotherId4) {
1254  if (mcGrandmotherPdg1 != 111) continue; // 111 = pi0, 221 = eta
1255 
1256  TVector3 pi0start;
1257  grandmother1->GetStartVertex(pi0start);
1258  /*
1259  fhTest2_startvertexElectrons_gg->Fill(pi0start_i.Z());
1260  fhTest2_startvertexElectrons_gg->Fill(pi0start_j.Z());
1261  fhTest2_startvertexElectrons_gg->Fill(pi0start_k.Z());
1262  fhTest2_startvertexElectrons_gg->Fill(pi0start_l.Z());
1263  fhTest2_startvertexElectrons_all->Fill(pi0start_i.Z());
1264  fhTest2_startvertexElectrons_all->Fill(pi0start_j.Z());
1265  fhTest2_startvertexElectrons_all->Fill(pi0start_k.Z());
1266  fhTest2_startvertexElectrons_all->Fill(pi0start_l.Z());
1267  */
1268 
1269  // consider only electrons from the target (only then the momenta are correctly refitted/reconstructed)
1270  if (pi0start_i.Z() > 1 || pi0start_j.Z() > 1 || pi0start_k.Z() > 1
1271  || pi0start_l.Z() > 1)
1272  continue;
1273 
1274  Double_t invmass_mc = 0;
1275  Double_t invmass_reco = 0;
1278  fVector_mctrack[i],
1279  fVector_mctrack[j],
1280  fVector_mctrack[k],
1281  fVector_mctrack[l]);
1282  CbmAnaConversionKinematicParams params_reco =
1285  fVector_momenta[j],
1286  fVector_momenta[k],
1287  fVector_momenta[l]);
1288  invmass_mc = params_mc.fMinv;
1289  invmass_reco = params_reco.fMinv;
1290 
1291  fhTest2_2rich_invmass_gg_mc->Fill(invmass_mc);
1292  fhTest2_2rich_invmass_gg_refitted->Fill(invmass_reco);
1293  fhTest2_2rich_invmass_all_mc->Fill(invmass_mc);
1294  fhTest2_2rich_invmass_all_refitted->Fill(invmass_reco);
1295 
1296  fhTest2_2rich_pt_vs_rap_gg->Fill(params_reco.fPt,
1297  params_reco.fRapidity);
1298  fhTest2_2rich_pt_vs_rap_all->Fill(params_reco.fPt,
1299  params_reco.fRapidity);
1300 
1301  cout << "########################################################"
1302  "##############"
1303  << endl;
1304 
1305  fill++;
1306  }
1307  }
1308  }
1309  }
1310  }
1311  }
1312 }
1313 
1314 
1316  Int_t nofDaughters = 0;
1317  for (unsigned int i = 0; i < fVector_mctrack.size(); i++) {
1318  Int_t motherId_temp = fVector_mctrack[i]->GetMotherId();
1319  if (motherId == motherId_temp) nofDaughters++;
1320  }
1321 
1322 
1323  return nofDaughters;
1324 }
CbmRichPoint.h
CbmAnaConversionTest2::fRichPoints
TClonesArray * fRichPoints
Definition: CbmAnaConversionTest2.h:47
CbmAnaConversionTest2::fhTest2_invmass_gee_mc
TH1D * fhTest2_invmass_gee_mc
Definition: CbmAnaConversionTest2.h:77
CbmAnaConversionTest2::fVector_gt
std::vector< CbmGlobalTrack * > fVector_gt
Definition: CbmAnaConversionTest2.h:70
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmAnaConversionKinematicParams::fPt
Double_t fPt
Definition: CbmAnaConversionKinematicParams.h:19
CbmAnaConversionTest2::fhTest2_2rich_pt_vs_rap_gee
TH2D * fhTest2_2rich_pt_vs_rap_gee
Definition: CbmAnaConversionTest2.h:100
CbmAnaConversionTest2::fhTest2_3rich_electrons_thetaVSp_missing
TH2D * fhTest2_3rich_electrons_thetaVSp_missing
Definition: CbmAnaConversionTest2.h:112
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmAnaConversionTest2::fGlobalTracks
TClonesArray * fGlobalTracks
Definition: CbmAnaConversionTest2.h:53
CbmAnaConversionTest2::fhTest2_2rich_invmass_gg_mc
TH1D * fhTest2_2rich_invmass_gg_mc
Definition: CbmAnaConversionTest2.h:95
CbmAnaConversionTest2::fhTest2_invmass_RICHindex0
TH1D * fhTest2_invmass_RICHindex0
Definition: CbmAnaConversionTest2.h:62
CbmAnaConversionTest2::fhTest2_2rich_invmass_gee_refitted
TH1D * fhTest2_2rich_invmass_gee_refitted
Definition: CbmAnaConversionTest2.h:94
L1Field.h
CbmL1PFFitter
Definition: CbmL1PFFitter.h:31
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, std::vector< int > &pidHypo)
Definition: CbmL1PFFitter.cxx:81
CbmAnaConversionTest2::fVector_mctrack
std::vector< CbmMCTrack * > fVector_mctrack
Definition: CbmAnaConversionTest2.h:72
CbmAnaConversionTest2::fhTest2_electrons_pt_vs_p_noRICH
TH2D * fhTest2_electrons_pt_vs_p_noRICH
Definition: CbmAnaConversionTest2.h:108
CbmAnaConversionTest2::fhTest2_electrons_pt_vs_p
TH2D * fhTest2_electrons_pt_vs_p
Definition: CbmAnaConversionTest2.h:106
CbmGlobalTrack::GetRichRingIndex
Int_t GetRichRingIndex() const
Definition: CbmGlobalTrack.h:41
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmAnaConversionKinematicParams.h
CbmAnaConversionTest2.h
CbmGlobalTrack.h
CbmAnaConversionTest2::fhTest2_3rich_electrons_theta_included
TH1D * fhTest2_3rich_electrons_theta_included
Definition: CbmAnaConversionTest2.h:109
CbmAnaConversionTest2::fVector_richIndex
std::vector< int > fVector_richIndex
Definition: CbmAnaConversionTest2.h:74
CbmAnaConversionTest2::fPrimVertex
CbmVertex * fPrimVertex
Definition: CbmAnaConversionTest2.h:54
CbmRichRing.h
CbmAnaConversionTest2::fhTest2_invmass_all_mc
TH1D * fhTest2_invmass_all_mc
Definition: CbmAnaConversionTest2.h:81
CbmAnaConversionTest2::fhTest2_startvertexElectrons_gg
TH1D * fhTest2_startvertexElectrons_gg
Definition: CbmAnaConversionTest2.h:89
CbmAnaConversionTest2::fKFVertex
CbmKFVertex fKFVertex
Definition: CbmAnaConversionTest2.h:55
CbmAnaConversionTest2::InvariantMassTest_3RICH
void InvariantMassTest_3RICH()
Definition: CbmAnaConversionTest2.cxx:556
CbmAnaConversionTest2::fhTest2_pt_vs_rap_gee
TH2D * fhTest2_pt_vs_rap_gee
Definition: CbmAnaConversionTest2.h:84
CbmAnaConversionTest2::fhTest2_invmass_RICHindex3
TH1D * fhTest2_invmass_RICHindex3
Definition: CbmAnaConversionTest2.h:65
CbmAnaConversionCutSettings.h
CbmAnaConversionTest2::fhTest2_3rich_electrons_theta_missing
TH1D * fhTest2_3rich_electrons_theta_missing
Definition: CbmAnaConversionTest2.h:110
CbmAnaConversionTest2::fhTest2_pt_vs_rap_all
TH2D * fhTest2_pt_vs_rap_all
Definition: CbmAnaConversionTest2.h:86
CbmAnaConversionTest2::Finish
void Finish()
Definition: CbmAnaConversionTest2.cxx:429
CbmStsTrack.h
Data class for STS tracks.
CbmGlobalTrack::GetStsTrackIndex
Int_t GetStsTrackIndex() const
Definition: CbmGlobalTrack.h:38
CbmAnaConversionTest2::fhTest2_startvertexElectrons_gee
TH1D * fhTest2_startvertexElectrons_gee
Definition: CbmAnaConversionTest2.h:88
CbmAnaConversionTest2::Exec
void Exec()
Definition: CbmAnaConversionTest2.cxx:440
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
CbmAnaConversionTest2::InvariantMassTest_2RICH
void InvariantMassTest_2RICH()
Definition: CbmAnaConversionTest2.cxx:985
CbmAnaConversionTest2::fhTest2_invmass_gg_mc
TH1D * fhTest2_invmass_gg_mc
Definition: CbmAnaConversionTest2.h:79
CbmAnaConversionTest2::fhTest2_2rich_invmass_all_refitted
TH1D * fhTest2_2rich_invmass_all_refitted
Definition: CbmAnaConversionTest2.h:98
CbmAnaConversionTest2::InitHistos
void InitHistos()
Definition: CbmAnaConversionTest2.cxx:138
CbmAnaConversionTest2::fVector_gtIndex
std::vector< int > fVector_gtIndex
Definition: CbmAnaConversionTest2.h:73
CbmAnaConversionKinematicParams::fMinv
Double_t fMinv
Definition: CbmAnaConversionKinematicParams.h:21
CbmTrackMatchNew.h
CbmVertex
Definition: CbmVertex.h:26
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
CbmAnaConversionTest2::fHistoList_test2
std::vector< TH1 * > fHistoList_test2
Definition: CbmAnaConversionTest2.h:58
CbmAnaConversionTest2::CbmAnaConversionTest2
CbmAnaConversionTest2()
Definition: CbmAnaConversionTest2.cxx:31
CbmAnaConversionTest2::fMcTracks
TClonesArray * fMcTracks
Definition: CbmAnaConversionTest2.h:50
CbmL1PFFitter.h
CbmAnaConversionTest2::fhTest2_3rich_electrons_thetaVSp_included
TH2D * fhTest2_3rich_electrons_thetaVSp_included
Definition: CbmAnaConversionTest2.h:111
CbmAnaConversionTest2::fhTest2_invmass_gg_refitted
TH1D * fhTest2_invmass_gg_refitted
Definition: CbmAnaConversionTest2.h:80
CbmAnaConversionTest2::fhTest2_invmass_RICHindex2
TH1D * fhTest2_invmass_RICHindex2
Definition: CbmAnaConversionTest2.h:64
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition: CbmMCTrack.h:182
CbmAnaConversionKinematicParams::KinematicParams_4particles_MC
static CbmAnaConversionKinematicParams KinematicParams_4particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
Definition: CbmAnaConversionKinematicParams.h:30
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmAnaConversionTest2::fhTest2_electrons_pt_vs_p_withRICH
TH2D * fhTest2_electrons_pt_vs_p_withRICH
Definition: CbmAnaConversionTest2.h:107
CbmAnaConversionTest2::fhTest2_2rich_invmass_gg_refitted
TH1D * fhTest2_2rich_invmass_gg_refitted
Definition: CbmAnaConversionTest2.h:96
CbmAnaConversionTest2::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmAnaConversionTest2.h:49
CbmAnaConversionTest2::~CbmAnaConversionTest2
virtual ~CbmAnaConversionTest2()
Definition: CbmAnaConversionTest2.cxx:81
CbmAnaConversionTest2::fStsTracks
TClonesArray * fStsTracks
Definition: CbmAnaConversionTest2.h:51
CbmAnaConversionTest2::fhTest2_invmass_all_refitted
TH1D * fhTest2_invmass_all_refitted
Definition: CbmAnaConversionTest2.h:82
CbmAnaConversionTest2::fhTest2_2rich_invmass_all_mc
TH1D * fhTest2_2rich_invmass_all_mc
Definition: CbmAnaConversionTest2.h:97
CbmAnaConversionKinematicParams
Definition: CbmAnaConversionKinematicParams.h:16
CbmAnaConversionTest2::fhTest2_2rich_pt_vs_rap_gg
TH2D * fhTest2_2rich_pt_vs_rap_gg
Definition: CbmAnaConversionTest2.h:101
CbmGlobalTrack
Definition: CbmGlobalTrack.h:26
CbmAnaConversionTest2::fhTest2_startvertexElectrons_all
TH1D * fhTest2_startvertexElectrons_all
Definition: CbmAnaConversionTest2.h:90
CbmAnaConversionTest2::fhTest2_invmass_gee_refitted
TH1D * fhTest2_invmass_gee_refitted
Definition: CbmAnaConversionTest2.h:78
CbmMCTrack.h
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmAnaConversionTest2::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition: CbmAnaConversionTest2.h:52
CbmAnaConversionTest2::fhTest2_invmass_RICHindex4
TH1D * fhTest2_invmass_RICHindex4
Definition: CbmAnaConversionTest2.h:66
CbmAnaConversionTest2::fhTest2_pt_vs_rap_gg
TH2D * fhTest2_pt_vs_rap_gg
Definition: CbmAnaConversionTest2.h:85
CbmAnaConversionTest2::fRichRings
TClonesArray * fRichRings
Definition: CbmAnaConversionTest2.h:48
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmLitGlobalElectronId.h
CbmAnaConversionTest2::fhTest2_invmass_RICHindex1
TH1D * fhTest2_invmass_RICHindex1
Definition: CbmAnaConversionTest2.h:63
CbmAnaConversionTest2::Init
void Init()
Definition: CbmAnaConversionTest2.cxx:84
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmAnaConversionTest2::fhTest2_2rich_invmass_gee_mc
TH1D * fhTest2_2rich_invmass_gee_mc
Definition: CbmAnaConversionTest2.h:93
CbmAnaConversionTest2::NofDaughters
Int_t NofDaughters(Int_t motherId)
Definition: CbmAnaConversionTest2.cxx:1315
CbmAnaConversionTest2::fVector_momenta
std::vector< TVector3 > fVector_momenta
Definition: CbmAnaConversionTest2.h:71
CbmAnaConversionTest2::fhTest2_2rich_pt_vs_rap_all
TH2D * fhTest2_2rich_pt_vs_rap_all
Definition: CbmAnaConversionTest2.h:102
CbmAnaConversionKinematicParams::fRapidity
Double_t fRapidity
Definition: CbmAnaConversionKinematicParams.h:20