28#include "FairRootManager.h"
30#include "TDirectory.h"
43 , fGlobalTracks(nullptr)
45 , fStsTrackMatches(nullptr)
46 , El_Photon_Eta_refmomentum()
47 , El_Photon_Eta_MCtrack()
49 , El_Photon_Pion_Eta_refmomentum()
50 , El_Photon_Pion_Eta_MCtrack()
51 , El_Photon_Pion_Eta_Id()
52 , Pion_Eta_refmomentum()
55 , All_El_refmomentum()
58 , All_Pion_refmomentum()
67 , ECPGA_leptons_RefMom()
70 , ECPGA_pions_RefMom()
73 , InvMass_eta_gg_mc(nullptr)
74 , InvMass_eta_gg_reffited(nullptr)
75 , InvMassPhoton_eta_gg_mc(nullptr)
76 , InvMassPhoton_eta_gg_reffited(nullptr)
77 , OpeningAnglePhoton_eta_gg_mc(nullptr)
78 , OpeningAnglePhoton_eta_gg_reffited(nullptr)
79 , OpeningAngle_eta_gg_between_gg_mc(nullptr)
80 , OpeningAngle_eta_gg_between_gg_reffited(nullptr)
81 , InvMass_eta_gg_allcombinations_mc(nullptr)
82 , InvMass_eta_gg_allcombinations_reffited(nullptr)
84 , InvMass_eta_gg_reco_aftercuts(nullptr)
85 , rap_vs_pt_eta_gg_reco_aftercuts(nullptr)
86 , rap_vs_pt_NOTeta_gg_reco_aftercuts(nullptr)
87 , fHistoList_eta_ppg()
88 , InvMass_eta_ppg_mc(nullptr)
89 , InvMass_eta_ppg_reffited(nullptr)
90 , InvMassPhoton_eta_ppg_mc(nullptr)
91 , InvMassPhoton_eta_ppg_reffited(nullptr)
92 , OpeningAnglePhoton_eta_ppg_mc(nullptr)
93 , OpeningAnglePhoton_eta_ppg_reffited(nullptr)
94 , InvMass_eta_ppg_allcombinations_mc(nullptr)
95 , InvMass_eta_ppg_allcombinations_reffited(nullptr)
96 , Pion_P_fromEta_reco(nullptr)
97 , Pion_P_elsewhere_reco(nullptr)
98 , Pion_Pt_fromEta_reco(nullptr)
99 , Pion_Pt_elsewhere_reco(nullptr)
100 , OA_betweenPions_fromEta_mc(nullptr)
101 , OA_betweenPions_fromEta_reco(nullptr)
102 , OA_betweenPions_fromEta_reco_wrongcombinations(nullptr)
103 , EMT_eta_ppg(nullptr)
104 , EMT_eta_three_body(nullptr)
105 , InvMass_eta_ppg_reco_aftercuts(nullptr)
106 , rap_vs_pt_eta_ppg_reco_aftercuts(nullptr)
107 , rap_vs_pt_NOTeta_ppg_reco_aftercuts(nullptr)
108 , fHistoList_eta_ppp()
109 , InvMass_eta_ppp_mc(nullptr)
110 , InvMass_eta_ppp_reffited(nullptr)
111 , InvMass_eta_Npion_mc(nullptr)
112 , InvMass_eta_Npion_reffited(nullptr)
113 , InvMass_eta_ppp_allcombinations_mc(nullptr)
114 , InvMass_eta_ppp_allcombinations_reffited(nullptr)
116 , EMT_gg_pair_momenta()
118 , EMT_ppg_ee_pair_momenta()
120 , EMT_ppg_pp_pair_momenta()
121 , EMT_ppg_positive_pion_Event()
122 , EMT_ppg_positive_pion_momenta()
123 , EMT_ppg_negative_pion_Event()
124 , EMT_ppg_negative_pion_momenta()
134 FairRootManager* ioman = FairRootManager::Instance();
135 if (
nullptr == ioman) { Fatal(
"CbmKresEtaMCAnalysis::Init",
"RootManager not instantised!"); }
137 fMcTracks = (TClonesArray*) ioman->GetObject(
"MCTrack");
138 if (
nullptr ==
fMcTracks) { Fatal(
"CbmKresEtaMCAnalysis::Init",
"No MCTrack array!"); }
140 fGlobalTracks = (TClonesArray*) ioman->GetObject(
"GlobalTrack");
141 if (
nullptr ==
fGlobalTracks) { Fatal(
"CbmKresEtaMCAnalysis::Init",
"No GlobalTrack array!"); }
143 fStsTracks = (TClonesArray*) ioman->GetObject(
"StsTrack");
144 if (
nullptr ==
fStsTracks) { Fatal(
"CbmKresEtaMCAnalysis::Init",
"No StsTrack array!"); }
147 if (
nullptr ==
fStsTrackMatches) { Fatal(
"CbmKresEtaMCAnalysis::Init",
"No StsTrackMatch array!"); }
156 Int_t nofMcTracks =
fMcTracks->GetEntriesFast();
157 for (
int i = 0; i < nofMcTracks; i++) {
159 if (mctrack ==
nullptr)
continue;
200 for (Int_t iTr = 0; iTr < ngTracks; iTr++) {
202 if (
nullptr == gTrack)
continue;
207 if (stsInd < 0)
continue;
209 if (stsTrack ==
nullptr)
continue;
211 if (stsMatch ==
nullptr)
continue;
214 if (McTrackId < 0)
continue;
216 if (mcTrack ==
nullptr)
continue;
219 if (TMath::Abs(pdg) != 11 && TMath::Abs(pdg) != 211)
continue;
236 if (TMath::Abs(pdg) == 11) {
241 if (TMath::Abs(pdg) == 211) {
249 if (motherId == -1)
continue;
251 if (mcMotherTrack ==
nullptr)
continue;
253 if (TMath::Abs(pdg) == 11 && mcMotherTrack->
GetPdgCode() == 22) {
256 if (mcGrTrack ==
nullptr)
continue;
273 if (TMath::Abs(pdg) == 211 && mcMotherTrack->
GetPdgCode() == 221) {
298 if (Event % 1000 == 0) {
304 if (Event % 10 == 0) {
322 Double_t openingAngle;
323 TLorentzVector gamma1;
324 TLorentzVector gamma2;
345 Double_t angle = gamma1.Angle(gamma2.Vect());
346 openingAngle = 180. * angle / TMath::Pi();
352 TVector3 electron3, TVector3 electron4)
354 double M2El = 2.6112004954086e-7;
355 Double_t energy1 = TMath::Sqrt(electron1.Mag2() + M2El);
356 TLorentzVector lorVec1(electron1, energy1);
358 Double_t energy2 = TMath::Sqrt(electron2.Mag2() + M2El);
359 TLorentzVector lorVec2(electron2, energy2);
361 Double_t energy3 = TMath::Sqrt(electron3.Mag2() + M2El);
362 TLorentzVector lorVec3(electron3, energy3);
364 Double_t energy4 = TMath::Sqrt(electron4.Mag2() + M2El);
365 TLorentzVector lorVec4(electron4, energy4);
367 TLorentzVector lorPhoton1 = lorVec1 + lorVec2;
368 TLorentzVector lorPhoton2 = lorVec3 + lorVec4;
370 Double_t angleBetweenPhotons = lorPhoton1.Angle(lorPhoton2.Vect());
371 Double_t theta = 180. * angleBetweenPhotons / TMath::Pi();
382 if (MC.size() < 4)
return;
384 for (
size_t i = 0; i < MC.size(); i++) {
385 for (
size_t j = i + 1; j < MC.size(); j++) {
386 for (
size_t k = j + 1; k < MC.size(); k++) {
387 for (
size_t l = k + 1; l < MC.size(); l++) {
389 int pdg1 = MC.at(i)->GetPdgCode();
390 int pdg2 = MC.at(j)->GetPdgCode();
391 int pdg3 = MC.at(k)->GetPdgCode();
392 int pdg4 = MC.at(l)->GetPdgCode();
394 if (pdg1 + pdg2 + pdg3 + pdg4 != 0)
continue;
395 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11 || TMath::Abs(pdg3) != 11 || TMath::Abs(pdg4) != 11)
398 int motherId1 = MC.at(i)->GetMotherId();
399 int motherId2 = MC.at(j)->GetMotherId();
400 int motherId3 = MC.at(k)->GetMotherId();
401 int motherId4 = MC.at(l)->GetMotherId();
403 int mcId1 = Id.at(i);
404 int mcId2 = Id.at(j);
405 int mcId3 = Id.at(k);
406 int mcId4 = Id.at(l);
408 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
410 if (motherId1 == -1 || motherId2 == -1 || motherId3 == -1 || motherId4 == -1)
continue;
422 if (mcMotherPdg1 != 22 || mcMotherPdg2 != 22 || mcMotherPdg3 != 22 || mcMotherPdg4 != 22)
continue;
429 if (grandmotherId1 == -1)
continue;
432 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
433 && grandmotherId1 == grandmotherId4) {
436 int mcGrandmotherPdg1 = GrMother1->
GetPdgCode();
437 if (mcGrandmotherPdg1 != 221)
continue;
440 Double_t InvMass_reco =
442 cout <<
"Decay eta -> gamma gamma -> e+e- e+e- detected!\t\t mc mass: " << InvMass_true
443 <<
"\t, reco mass: " << InvMass_reco << endl;
444 cout <<
"motherids: " << motherId1 <<
"/" << motherId2 <<
"/" << motherId3 <<
"/" << motherId4 << endl;
445 cout <<
"grandmotherid: " << grandmotherId1 <<
"/" << grandmotherId2 <<
"/" << grandmotherId3 <<
"/"
446 << grandmotherId4 << endl;
447 cout <<
"pdgs: " << mcGrandmotherPdg1 <<
"-->" << mcMotherPdg1 <<
"/" << mcMotherPdg2 <<
"/" << mcMotherPdg3
448 <<
"/" << mcMotherPdg4 <<
"-->" << pdg1 <<
"/" << pdg2 <<
"/" << pdg3 <<
"/" << pdg4 << endl;
450 gg[0]->Fill(InvMass_true);
451 gg[1]->Fill(InvMass_reco);
453 cout <<
"\t \t mc true info: " << endl;
454 cout <<
"particle 1: \t" << MC.at(i)->GetPdgCode() <<
";\t pt = " << MC.at(i)->GetPt()
455 <<
";\t X = " << MC.at(i)->GetStartX() <<
";\t Y = " << MC.at(i)->GetStartY()
456 <<
";\t Z = " << MC.at(i)->GetStartZ() <<
";\t E = " << MC.at(i)->GetEnergy() << endl;
457 cout <<
"particle 2: \t" << MC.at(j)->GetPdgCode() <<
";\t pt = " << MC.at(j)->GetPt()
458 <<
";\t X = " << MC.at(j)->GetStartX() <<
";\t Y = " << MC.at(j)->GetStartY()
459 <<
";\t Z = " << MC.at(j)->GetStartZ() <<
";\t E = " << MC.at(j)->GetEnergy() << endl;
460 cout <<
"particle 3: \t" << MC.at(k)->GetPdgCode() <<
";\t pt = " << MC.at(k)->GetPt()
461 <<
";\t X = " << MC.at(k)->GetStartX() <<
";\t Y = " << MC.at(k)->GetStartY()
462 <<
";\t Z = " << MC.at(k)->GetStartZ() <<
";\t E = " << MC.at(k)->GetEnergy() << endl;
463 cout <<
"particle 4: \t" << MC.at(l)->GetPdgCode() <<
";\t pt = " << MC.at(l)->GetPt()
464 <<
";\t X = " << MC.at(l)->GetStartX() <<
";\t Y = " << MC.at(l)->GetStartY()
465 <<
";\t Z = " << MC.at(l)->GetStartZ() <<
";\t E = " << MC.at(l)->GetEnergy() << endl;
467 Double_t OpeningAngle1_mc = 0;
468 Double_t OpeningAngle2_mc = 0;
469 Double_t OpeningAngle1_refitted = 0;
470 Double_t OpeningAngle2_refitted = 0;
471 Double_t InvMass_realg1_mc = 0;
472 Double_t InvMass_realg2_mc = 0;
473 Double_t InvMass_realg1_refitted = 0;
474 Double_t InvMass_realg2_refitted = 0;
476 Double_t openingAngleBetweenGammasReco = 0;
478 if (motherId1 == motherId2) {
481 gg[4]->Fill(OpeningAngle1_mc);
482 gg[4]->Fill(OpeningAngle2_mc);
485 gg[5]->Fill(OpeningAngle1_refitted);
486 gg[5]->Fill(OpeningAngle2_refitted);
490 gg[2]->Fill(InvMass_realg1_mc);
491 gg[2]->Fill(InvMass_realg2_mc);
494 gg[3]->Fill(InvMass_realg1_refitted);
495 gg[3]->Fill(InvMass_realg2_refitted);
496 openingAngleBetweenGammasReco =
499 if (motherId1 == motherId3) {
502 gg[4]->Fill(OpeningAngle1_mc);
503 gg[4]->Fill(OpeningAngle2_mc);
506 gg[5]->Fill(OpeningAngle1_refitted);
507 gg[5]->Fill(OpeningAngle2_refitted);
511 gg[2]->Fill(InvMass_realg1_mc);
512 gg[2]->Fill(InvMass_realg2_mc);
515 gg[3]->Fill(InvMass_realg1_refitted);
516 gg[3]->Fill(InvMass_realg2_refitted);
517 openingAngleBetweenGammasReco =
520 if (motherId1 == motherId4) {
523 gg[4]->Fill(OpeningAngle1_mc);
524 gg[4]->Fill(OpeningAngle2_mc);
527 gg[5]->Fill(OpeningAngle1_refitted);
528 gg[5]->Fill(OpeningAngle2_refitted);
532 gg[2]->Fill(InvMass_realg1_mc);
533 gg[2]->Fill(InvMass_realg2_mc);
536 gg[3]->Fill(InvMass_realg1_refitted);
537 gg[3]->Fill(InvMass_realg2_refitted);
538 openingAngleBetweenGammasReco =
541 Double_t openingAngleBetweenGammas =
543 gg[6]->Fill(openingAngleBetweenGammas);
544 gg[7]->Fill(openingAngleBetweenGammasReco);
554 vector<Int_t> , vector<TVector3> RefMomEl,
555 vector<CbmMCTrack*> MCEl, vector<Int_t> ,
559 if (MCPion.size() < 2 || MCEl.size() < 2)
return;
563 for (
size_t i = 0; i < MCEl.size(); i++) {
564 for (
size_t j = i + 1; j < MCEl.size(); j++) {
565 int pdg1 = MCEl.at(i)->GetPdgCode();
566 int pdg2 = MCEl.at(j)->GetPdgCode();
568 if (pdg1 + pdg2 != 0)
continue;
569 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11)
continue;
571 int motherId1 = MCEl.at(i)->GetMotherId();
572 int motherId2 = MCEl.at(j)->GetMotherId();
574 if (motherId1 != motherId2)
continue;
575 if (motherId1 == -1 || motherId2 == -1)
continue;
583 if (grandmotherId1 == -1)
continue;
586 int mcGrandmotherPdg1 = GrMother1->
GetPdgCode();
587 if (mcGrandmotherPdg1 != 221)
continue;
590 for (
size_t k = 0; k < MCPion.size(); k++) {
591 for (
size_t l = k + 1; l < MCPion.size(); l++) {
592 int pdg3 = MCPion.at(k)->GetPdgCode();
593 int pdg4 = MCPion.at(l)->GetPdgCode();
595 if (pdg3 + pdg4 != 0)
continue;
596 if (TMath::Abs(pdg3) != 211 || TMath::Abs(pdg4) != 211)
continue;
598 int motherId3 = MCPion.at(k)->GetMotherId();
599 int motherId4 = MCPion.at(l)->GetMotherId();
601 if (motherId3 != motherId4)
continue;
602 if (motherId3 == -1 || motherId4 == -1)
continue;
603 if (motherId3 != grandmotherId1)
continue;
605 Double_t InvMass_true =
608 RefMomPion.at(k), RefMomPion.at(l));
609 cout <<
"Decay eta -> p+p- gamma -> p+p- e+e- detected!\t\t mc mass: " << InvMass_true
610 <<
"\t, reco mass: " << InvMass_reco << endl;
611 cout <<
"motherids: " << motherId1 <<
"/" << motherId2 <<
"/" << motherId3 <<
"/" << motherId4 << endl;
612 cout <<
"grandmotherid: " << grandmotherId1 << endl;
613 cout <<
"pdgs: " << mcGrandmotherPdg1 <<
"-->" << mother1->
GetPdgCode() <<
"/" << pdg3 <<
"/" << pdg4 <<
"-->"
614 << pdg1 <<
"/" << pdg2 <<
"/" << pdg3 <<
"/" << pdg4 << endl;
616 ppg[0]->Fill(InvMass_true);
617 ppg[1]->Fill(InvMass_reco);
619 cout <<
"\t \t mc true info: " << endl;
620 cout <<
"particle 1: \t" << MCEl.at(i)->GetPdgCode() <<
";\t pt = " << MCEl.at(i)->GetPt()
621 <<
";\t X = " << MCEl.at(i)->GetStartX() <<
";\t Y = " << MCEl.at(i)->GetStartY()
622 <<
";\t Z = " << MCEl.at(i)->GetStartZ() <<
";\t E = " << MCEl.at(i)->GetEnergy() << endl;
623 cout <<
"particle 2: \t" << MCEl.at(j)->GetPdgCode() <<
";\t pt = " << MCEl.at(j)->GetPt()
624 <<
";\t X = " << MCEl.at(j)->GetStartX() <<
";\t Y = " << MCEl.at(j)->GetStartY()
625 <<
";\t Z = " << MCEl.at(j)->GetStartZ() <<
";\t E = " << MCEl.at(j)->GetEnergy() << endl;
626 cout <<
"particle 3: \t" << MCPion.at(k)->GetPdgCode() <<
";\t pt = " << MCPion.at(k)->GetPt()
627 <<
";\t X = " << MCPion.at(k)->GetStartX() <<
";\t Y = " << MCPion.at(k)->GetStartY()
628 <<
";\t Z = " << MCPion.at(k)->GetStartZ() <<
";\t E = " << MCPion.at(k)->GetEnergy() << endl;
629 cout <<
"particle 4: \t" << MCPion.at(l)->GetPdgCode() <<
";\t pt = " << MCPion.at(l)->GetPt()
630 <<
";\t X = " << MCPion.at(l)->GetStartX() <<
";\t Y = " << MCPion.at(l)->GetStartY()
631 <<
";\t Z = " << MCPion.at(l)->GetStartZ() <<
";\t E = " << MCPion.at(l)->GetEnergy() << endl;
638 ppg[2]->Fill(InvMass_g_mc);
639 ppg[3]->Fill(InvMass_g_reffited);
640 ppg[4]->Fill(OpeningAngle_g_mc);
641 ppg[5]->Fill(OpeningAngle_g_refitted);
650 vector<Int_t> IdNeutral, vector<TVector3> RefMomPion,
651 vector<CbmMCTrack*> MCPion, vector<Int_t> ,
654 if (MCNeutral.size() < 4 || MCPion.size() < 2)
return;
656 int fDecayedParticlePdg = 221;
660 for (
size_t i = 0; i < MCNeutral.size(); i++) {
661 for (
size_t j = i + 1; j < MCNeutral.size(); j++) {
662 for (
size_t k = j + 1; k < MCNeutral.size(); k++) {
663 for (
size_t l = k + 1; l < MCNeutral.size(); l++) {
665 int pdg1 = MCNeutral.at(i)->GetPdgCode();
666 int pdg2 = MCNeutral.at(j)->GetPdgCode();
667 int pdg3 = MCNeutral.at(k)->GetPdgCode();
668 int pdg4 = MCNeutral.at(l)->GetPdgCode();
670 if (pdg1 + pdg2 + pdg3 + pdg4 != 0)
continue;
671 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11 || TMath::Abs(pdg3) != 11 || TMath::Abs(pdg4) != 11)
674 int motherId1 = MCNeutral.at(i)->GetMotherId();
675 int motherId2 = MCNeutral.at(j)->GetMotherId();
676 int motherId3 = MCNeutral.at(k)->GetMotherId();
677 int motherId4 = MCNeutral.at(l)->GetMotherId();
679 int mcId1 = IdNeutral.at(i);
680 int mcId2 = IdNeutral.at(j);
681 int mcId3 = IdNeutral.at(k);
682 int mcId4 = IdNeutral.at(l);
684 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
686 if (motherId1 == -1 || motherId2 == -1 || motherId3 == -1 || motherId4 == -1)
continue;
698 if (mcMotherPdg1 != 22 || mcMotherPdg2 != 22 || mcMotherPdg3 != 22 || mcMotherPdg4 != 22)
continue;
705 if (grandmotherId1 == -1)
continue;
707 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
708 && grandmotherId1 == grandmotherId4) {
710 int mcGrandmotherPdg1 = GrMother1->
GetPdgCode();
711 if (mcGrandmotherPdg1 != 111)
continue;
715 if (TopTrack->
GetPdgCode() != fDecayedParticlePdg)
continue;
718 for (
size_t m = 0; m < MCPion.size(); m++) {
719 for (
size_t s = m + 1; s < MCPion.size(); s++) {
720 int pdg5 = MCPion.at(m)->GetPdgCode();
721 int pdg6 = MCPion.at(s)->GetPdgCode();
723 if (pdg5 + pdg6 != 0)
continue;
724 if (TMath::Abs(pdg5) != 211 || TMath::Abs(pdg6) != 211)
continue;
726 int motherId5 = MCPion.at(m)->GetMotherId();
727 int motherId6 = MCPion.at(s)->GetMotherId();
729 if (motherId5 != motherId6)
continue;
730 if (motherId5 == -1 || motherId6 == -1)
continue;
731 if (motherId5 != GrMother1->
GetMotherId())
continue;
734 MCNeutral.at(i), MCNeutral.at(j), MCNeutral.at(k), MCNeutral.at(l), MCPion.at(m), MCPion.at(s));
736 RefMomNeutral.at(i), RefMomNeutral.at(j), RefMomNeutral.at(k), RefMomNeutral.at(l), RefMomPion.at(m),
738 cout <<
"Decay eta -> p+p- p0 -> p+p- e+e-e+e- detected!\t\t "
740 << InvMass_six_true <<
"\t, reco mass: " << InvMass_six_reco << endl;
741 cout <<
"motherids for neutral: " << motherId1 <<
"/" << motherId2 <<
"/" << motherId3 <<
"/"
742 << motherId4 << endl;
743 cout <<
"grandmotherid: " << GrMother1->
GetMotherId() <<
"/" << motherId5 <<
"/" << motherId6 << endl;
744 cout <<
"pdgs: " << TopTrack->
GetPdgCode() <<
"-->" << GrMother1->
GetPdgCode() <<
"/" << pdg5 <<
"/"
747 ppp[0]->Fill(InvMass_six_true);
748 ppp[1]->Fill(InvMass_six_reco);
750 cout <<
"\t \t mc true info: " << endl;
751 cout <<
"particle 1: \t" << MCNeutral.at(i)->GetPdgCode() <<
";\t pt = " << MCNeutral.at(i)->GetPt()
752 <<
";\t X = " << MCNeutral.at(i)->GetStartX() <<
";\t Y = " << MCNeutral.at(i)->GetStartY()
753 <<
";\t Z = " << MCNeutral.at(i)->GetStartZ() <<
";\t E = " << MCNeutral.at(i)->GetEnergy()
755 cout <<
"particle 2: \t" << MCNeutral.at(j)->GetPdgCode() <<
";\t pt = " << MCNeutral.at(j)->GetPt()
756 <<
";\t X = " << MCNeutral.at(j)->GetStartX() <<
";\t Y = " << MCNeutral.at(j)->GetStartY()
757 <<
";\t Z = " << MCNeutral.at(j)->GetStartZ() <<
";\t E = " << MCNeutral.at(j)->GetEnergy()
759 cout <<
"particle 3: \t" << MCNeutral.at(k)->GetPdgCode() <<
";\t pt = " << MCNeutral.at(k)->GetPt()
760 <<
";\t X = " << MCNeutral.at(k)->GetStartX() <<
";\t Y = " << MCNeutral.at(k)->GetStartY()
761 <<
";\t Z = " << MCNeutral.at(k)->GetStartZ() <<
";\t E = " << MCNeutral.at(k)->GetEnergy()
763 cout <<
"particle 4: \t" << MCNeutral.at(l)->GetPdgCode() <<
";\t pt = " << MCNeutral.at(l)->GetPt()
764 <<
";\t X = " << MCNeutral.at(l)->GetStartX() <<
";\t Y = " << MCNeutral.at(l)->GetStartY()
765 <<
";\t Z = " << MCNeutral.at(l)->GetStartZ() <<
";\t E = " << MCNeutral.at(l)->GetEnergy()
767 cout <<
"particle 5: \t" << MCPion.at(m)->GetPdgCode() <<
";\t pt = " << MCPion.at(m)->GetPt()
768 <<
";\t X = " << MCPion.at(m)->GetStartX() <<
";\t Y = " << MCPion.at(m)->GetStartY()
769 <<
";\t Z = " << MCPion.at(m)->GetStartZ() <<
";\t E = " << MCPion.at(m)->GetEnergy() << endl;
770 cout <<
"particle 6: \t" << MCPion.at(s)->GetPdgCode() <<
";\t pt = " << MCPion.at(s)->GetPt()
771 <<
";\t X = " << MCPion.at(s)->GetStartX() <<
";\t Y = " << MCPion.at(s)->GetStartY()
772 <<
";\t Z = " << MCPion.at(s)->GetStartZ() <<
";\t E = " << MCPion.at(s)->GetEnergy() << endl;
775 MCNeutral.at(k), MCNeutral.at(l));
777 RefMomNeutral.at(i), RefMomNeutral.at(j), RefMomNeutral.at(k), RefMomNeutral.at(l));
778 ppp[2]->Fill(InvMass_true);
779 ppp[3]->Fill(InvMass_reco);
791 vector<TVector3> RefMom, vector<CbmMCTrack*> MC,
792 vector<Int_t> Id, vector<TH1*> gg)
794 if (MC.size() < 4)
return;
800 for (
size_t i = 0; i < MC.size(); i++) {
801 for (
size_t j = i + 1; j < MC.size(); j++) {
802 int pdg1 = MC.at(i)->GetPdgCode();
803 int pdg2 = MC.at(j)->GetPdgCode();
804 if (pdg1 + pdg2 != 0)
continue;
805 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11)
continue;
810 if (OpeningAngle_g_refitted > OpeningAngleCut || InvMass_g_refitted > GammaInvMassCut)
continue;
819 fMCId.push_back(Id.at(i));
820 fMCId.push_back(Id.at(j));
832 if (
EDGA_MC.size() < 2)
return;
833 for (
size_t gamma1 = 0; gamma1 <
EDGA_RefMom.size(); gamma1++) {
834 for (
size_t gamma2 = gamma1 + 1; gamma2 <
EDGA_RefMom.size(); gamma2++) {
846 int mcId1 =
EDGA_Id[gamma1][0];
847 int mcId2 =
EDGA_Id[gamma1][1];
848 int mcId3 =
EDGA_Id[gamma2][0];
849 int mcId4 =
EDGA_Id[gamma2][1];
851 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
855 if (openingAngleBetweenGammasReco < 10 || openingAngleBetweenGammasReco > 40)
continue;
861 gg[8]->Fill(InvMass_true);
862 gg[9]->Fill(InvMass_reco);
879 if (GrandMothTrack1->
GetPdgCode() == 221) CorrectEta = 1;
883 if (CorrectEta == 1) {
895 int Event, vector<TVector3> RefMomPion,
896 vector<CbmMCTrack*> MCPion, vector<Int_t> ,
897 vector<TVector3> RefMomEl, vector<CbmMCTrack*> MCEl,
898 vector<Int_t> IdEl, vector<TH1*> ppg, vector<TH1*> ppp)
900 if (MCPion.size() < 2 || MCEl.size() < 2)
return;
906 for (
size_t i = 0; i < MCEl.size(); i++) {
907 for (
size_t j = i + 1; j < MCEl.size(); j++) {
908 int pdg1 = MCEl.at(i)->GetPdgCode();
909 int pdg2 = MCEl.at(j)->GetPdgCode();
910 if (pdg1 + pdg2 != 0)
continue;
911 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11)
continue;
916 if (OpeningAngle_g_refitted > OpeningAngleCut || InvMass_g_refitted > GammaInvMassCut)
continue;
925 fMCId.push_back(IdEl.at(i));
926 fMCId.push_back(IdEl.at(j));
940 for (
size_t k = 0; k < MCPion.size(); k++) {
941 if (MCPion.at(k)->GetPdgCode() == 211) {
949 for (
size_t l = k + 1; l < MCPion.size(); l++) {
950 int pdg3 = MCPion.at(k)->GetPdgCode();
951 int pdg4 = MCPion.at(l)->GetPdgCode();
952 if (pdg3 + pdg4 != 0)
continue;
953 if (TMath::Abs(pdg3) != 211 || TMath::Abs(pdg4) != 211)
continue;
956 double Pions_angle_reco =
959 if (Pions_angle_reco > 20)
continue;
992 int correctPionspair = 0;
995 if (EtaTrack->
GetPdgCode() == 221) correctPionspair = 1;
998 if (correctPionspair == 1) {
1030 ppg[6]->Fill(InvMass_true);
1031 ppg[7]->Fill(InvMass_reco);
1040 if (EtaTrack->
GetPdgCode() == 221) CorrectEta = 1;
1044 if (CorrectEta == 1) {
1075 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
1082 if (InvMass_pi_reco < 0.12 || InvMass_pi_reco > 0.15)
continue;
1084 for (
size_t k = 0; k < MCPion.size(); k++) {
1085 for (
size_t l = k + 1; l < MCPion.size(); l++) {
1086 int pdg5 = MCPion.at(k)->GetPdgCode();
1087 int pdg6 = MCPion.at(l)->GetPdgCode();
1089 if (pdg5 + pdg6 != 0)
continue;
1090 if (TMath::Abs(pdg5) != 211 || TMath::Abs(pdg6) != 211)
continue;
1092 double InvMass_six_true =
1094 double InvMass_six_reco =
1096 double Pions_angle_reco =
1099 if (Pions_angle_reco > 20)
continue;
1101 ppp[4]->Fill(InvMass_six_true);
1102 ppp[5]->Fill(InvMass_six_reco);
1113 cout <<
"Mixing for eta->gg channel - nof entries " << nof << endl;
1114 for (Int_t a = 0; a < nof - 1; a++) {
1115 for (Int_t b = a + 1; b < nof; b++) {
1123 if (openingAngleBetweenGammasReco < 10 || openingAngleBetweenGammasReco > 40)
continue;
1135 cout <<
"Mixing for eta->(pp)g channel - nof lepton pairs = " << nof_leptons <<
"; nof pion pairs = " << nof_pions
1137 for (Int_t a = 0; a < nof_leptons; a++) {
1138 for (Int_t b = 0; b < nof_pions; b++) {
1158 cout <<
"Mixing 3 bodies for eta-> p p g channel - nof photons = " << nof_photons
1159 <<
"; nof +pions = " << nof_positive <<
"; nof -pions = " << nof_negative << endl;
1161 for (Int_t p = 0; p < nof_positive; p++) {
1162 for (Int_t m = 0; m < nof_negative; m++) {
1169 if (Pions_angle_reco > 20)
continue;
1171 for (Int_t a = 0; a < nof_photons; a++) {
1188 gDirectory->mkdir(
"EtaMCAnalysis");
1189 gDirectory->cd(
"EtaMCAnalysis");
1191 gDirectory->mkdir(
"eta_gg");
1192 gDirectory->cd(
"eta_gg");
1196 gDirectory->cd(
"..");
1198 gDirectory->mkdir(
"eta_ppg");
1199 gDirectory->cd(
"eta_ppg");
1203 gDirectory->cd(
"..");
1205 gDirectory->mkdir(
"eta_ppp");
1206 gDirectory->cd(
"eta_ppp");
1210 gDirectory->cd(
"..");
1213 gDirectory->cd(
"..");
1219 InvMass_eta_gg_mc =
new TH1D(
"InvMass_eta_gg_mc",
"InvMass_eta_gg_mc; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1222 new TH1D(
"InvMass_eta_gg_reffited",
"InvMass_eta_gg_reffited; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1225 new TH1D(
"InvMassPhoton_eta_gg_mc",
"InvMassPhoton_eta_gg_mc; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1228 "InvMassPhoton_eta_gg_reffited",
"InvMassPhoton_eta_gg_reffited; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1231 new TH1D(
"OpeningAnglePhoton_eta_gg_mc",
"OpeningAnglePhoton_eta_gg_mc (between e+e- from #gamma); angle [deg];#",
1235 "OpeningAnglePhoton_eta_gg_reffited (between e+e- from #gamma); "
1240 "OpeningAngle_eta_gg_between_gg_mc (between #gamma#gamma from "
1241 "#eta); angle [deg];#",
1245 "OpeningAngle_eta_gg_between_gg_reffited (between #gamma#gamma "
1246 "from #eta); angle [deg];#",
1250 new TH1D(
"InvMass_eta_gg_allcombinations_mc",
"InvMass_eta_gg_allcombinations_mc; invariant mass in GeV/c^{2};#",
1254 new TH1D(
"InvMass_eta_gg_allcombinations_reffited",
1255 "InvMass_eta_gg_allcombinations_reffited; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1257 EMT_eta_gg =
new TH1D(
"EMT_eta_gg",
"EMT_eta_gg; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1260 "InvMass_eta_gg_reco_aftercuts",
"InvMass_eta_gg_reco_aftercuts; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1263 new TH2D(
"rap_vs_pt_eta_gg_reco_aftercuts",
"rap_vs_pt_eta_gg_reco_aftercuts; rapidity y; p_{t} in GeV/c ", 10, 0.,
1267 new TH2D(
"rap_vs_pt_NOTeta_gg_reco_aftercuts",
"rap_vs_pt_NOTeta_gg_reco_aftercuts; rapidity y; p_{t} in GeV/c ",
1268 10, 0., 4., 5, 0., 2.);
1274 new TH1D(
"InvMass_eta_ppg_mc",
"InvMass_eta_ppg_mc; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1277 new TH1D(
"InvMass_eta_ppg_reffited",
"InvMass_eta_ppg_reffited; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1280 new TH1D(
"InvMassPhoton_eta_ppg_mc",
"InvMassPhoton_eta_ppg_mc; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1283 "InvMassPhoton_eta_ppg_reffited",
"InvMassPhoton_eta_ppg_reffited; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1286 new TH1D(
"OpeningAnglePhoton_eta_ppg_mc",
"OpeningAnglePhoton_eta_ppg_mc (between e+e- from #gamma); angle [deg];#",
1290 "OpeningAnglePhoton_eta_ppg_reffited (between e+e- from #gamma); "
1295 new TH1D(
"InvMass_eta_ppg_allcombinations_mc",
"InvMass_eta_ppg_allcombinations_mc; invariant mass in GeV/c^{2};#",
1299 new TH1D(
"InvMass_eta_ppg_allcombinations_reffited",
1300 "InvMass_eta_ppg_allcombinations_reffited; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1302 Pion_P_fromEta_reco =
new TH1D(
"Pion_P_fromEta_reco",
"Pion_P_fromEta_reco; pions P in GeV/c;#", 500, 0.0, 5.0);
1304 Pion_P_elsewhere_reco =
new TH1D(
"Pion_P_elsewhere_reco",
"Pion_P_elsewhere_reco; pions P in GeV/c;#", 500, 0.0, 5.0);
1307 new TH1D(
"Pion_Pt_fromEta_reco",
"Pion_Pt_fromEta_reco; pions P_{t} in GeV/c;#", 500, 0.0, 5.0);
1310 new TH1D(
"Pion_Pt_elsewhere_reco",
"Pion_Pt_elsewhere_reco; pions P_{t} in GeV/c;#", 500, 0.0, 5.0);
1313 "OA_betweenPions_fromEta_mc (between #p^{+}#p^{-} from #eta); "
1318 "OA_betweenPions_fromEta_reco (between #p^{+}#p^{-} from #eta); "
1323 "OA_betweenPions_fromEta_reco_wrongcombinations (between "
1324 "#p^{+}#p^{-} from #eta); angle [deg];#",
1327 EMT_eta_ppg =
new TH1D(
"EMT_eta_ppg",
"EMT_eta_ppg; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1330 new TH1D(
"EMT_eta_three_body",
"EMT_eta_three_body; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1333 "InvMass_eta_ppg_reco_aftercuts",
"InvMass_eta_ppg_reco_aftercuts; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1336 new TH2D(
"rap_vs_pt_eta_ppg_reco_aftercuts",
"rap_vs_pt_eta_ppg_reco_aftercuts; rapidity y; p_{t} in GeV/c ", 10,
1340 new TH2D(
"rap_vs_pt_NOTeta_ppg_reco_aftercuts",
"rap_vs_pt_NOTeta_ppg_reco_aftercuts; rapidity y; p_{t} in GeV/c ",
1341 10, 0., 4., 5, 0., 2.);
1347 new TH1D(
"InvMass_eta_ppp_mc",
"InvMass_eta_ppp_mc; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1350 new TH1D(
"InvMass_eta_ppp_reffited",
"InvMass_eta_ppp_reffited; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1353 new TH1D(
"InvMass_eta_Npion_mc",
"InvMass_eta_Npion_mc; invariant mass in GeV/c^{2};#", 30, 0.0, 0.3);
1356 new TH1D(
"InvMass_eta_Npion_reffited",
"InvMass_eta_Npion_reffited; invariant mass in GeV/c^{2};#", 30, 0., 0.3);
1359 new TH1D(
"InvMass_eta_ppp_allcombinations_mc",
"InvMass_eta_ppp_allcombinations_mc; invariant mass in GeV/c^{2};#",
1363 new TH1D(
"InvMass_eta_ppp_allcombinations_reffited",
1364 "InvMass_eta_ppp_allcombinations_reffited; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
Data class for STS tracks.
int32_t GetStsTrackIndex() const
TH1D * InvMass_eta_gg_reffited
vector< TH1 * > fHistoList_eta_gg
vector< CbmMCTrack * > El_Photon_Pion_Eta_MCtrack
Double_t CalculateOpeningAngleBetweenGammas_Reco(TVector3 electron1, TVector3 electron2, TVector3 electron3, TVector3 electron4)
virtual ~CbmKresEtaMCAnalysis()
TH1D * InvMass_eta_gg_allcombinations_reffited
TClonesArray * fGlobalTracks
TH1D * Pion_Pt_elsewhere_reco
Double_t CalculateOpeningAngleBetweenGammas_MC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2, CbmMCTrack *mctrack3, CbmMCTrack *mctrack4)
vector< int > El_Photon_Eta_Id
std::vector< std::vector< TVector3 > > EMT_ppg_pp_pair_momenta
vector< TVector3 > El_Photon_Eta_refmomentum
TH1D * InvMass_eta_ppp_mc
vector< CbmMCTrack * > fMCtracks
vector< CbmMCTrack * > All_El_MCtrack
TH1D * InvMass_eta_gg_allcombinations_mc
TH1D * Pion_Pt_fromEta_reco
void EtaChargedPionsGammaAnalysis(vector< TVector3 > RefMomPion, vector< CbmMCTrack * > MCPion, vector< Int_t > IdPion, vector< TVector3 > RefMomEl, vector< CbmMCTrack * > MCEl, vector< Int_t > IdEl, vector< TH1 * > ppg)
TH1D * InvMassPhoton_eta_ppg_mc
TH1D * OpeningAnglePhoton_eta_ppg_mc
TClonesArray * fStsTrackMatches
TH1D * InvMass_eta_ppg_allcombinations_mc
TH1D * OA_betweenPions_fromEta_reco
TH1D * InvMass_eta_ppg_mc
std::vector< int > EMT_ppg_pp_Event
TH1D * InvMass_eta_Npion_reffited
TH1D * InvMass_eta_Npion_mc
vector< TVector3 > Pion_Eta_refmomentum
void EtaChargedPionsGammaAnalysis_plusBG(double OpeningAngleCut, double GammaInvMassCut, int Event, vector< TVector3 > RefMomPion, vector< CbmMCTrack * > MCPion, vector< Int_t > IdPion, vector< TVector3 > RefMomEl, vector< CbmMCTrack * > MCEl, vector< Int_t > IdEl, vector< TH1 * > ppg, vector< TH1 * > ppp)
void Exec(int Event, double OpeningAngleCut, double GammaInvMassCut)
std::vector< std::vector< TVector3 > > ECPGA_pions_RefMom
std::vector< std::vector< CbmMCTrack * > > ECPGA_leptons_MC
TH1D * OpeningAnglePhoton_eta_gg_reffited
TH1D * EMT_eta_three_body
std::vector< int > EMT_ppg_positive_pion_Event
TH1D * OpeningAngle_eta_gg_between_gg_reffited
vector< int > Pion_Eta_Id
std::vector< TVector3 > EMT_ppg_positive_pion_momenta
vector< CbmMCTrack * > Pion_Eta_MCtrack
std::vector< std::vector< CbmMCTrack * > > EDGA_MC
std::vector< int > EMT_gg_Event
TH1D * InvMass_eta_gg_reco_aftercuts
void EtaPosNegNeutralPionsAnalysis(vector< TVector3 > RefMomNeutral, vector< CbmMCTrack * > MCNeutral, vector< Int_t > IdNeutral, vector< TVector3 > RefMomPion, vector< CbmMCTrack * > MCPion, vector< Int_t > IdPion, vector< TH1 * > ppp)
std::vector< std::vector< int > > EDGA_Id
std::vector< std::vector< TVector3 > > EMT_ppg_ee_pair_momenta
vector< TVector3 > frefmomenta
TH2D * rap_vs_pt_NOTeta_ppg_reco_aftercuts
TH1D * InvMass_eta_ppp_allcombinations_reffited
TH2D * rap_vs_pt_eta_gg_reco_aftercuts
TH1D * InvMass_eta_ppg_allcombinations_reffited
void EtaDoubleGammaAnalysis(vector< TVector3 > RefMom, vector< CbmMCTrack * > MC, vector< Int_t > Id, vector< TH1 * > gg)
vector< TVector3 > El_Photon_Pion_Eta_refmomentum
TH1D * InvMassPhoton_eta_gg_mc
TH1D * InvMassPhoton_eta_gg_reffited
std::vector< std::vector< TVector3 > > ECPGA_leptons_RefMom
std::vector< std::vector< TVector3 > > EMT_gg_pair_momenta
TH1D * OA_betweenPions_fromEta_reco_wrongcombinations
TH1D * InvMassPhoton_eta_ppg_reffited
TH1D * OpeningAnglePhoton_eta_gg_mc
std::vector< std::vector< CbmMCTrack * > > ECPGA_pions_MC
vector< TVector3 > All_El_refmomentum
TClonesArray * fStsTracks
void EtaDoubleGammaAnalysis_plusBG(double OpeningAngleCut, double GammaInvMassCut, int Event, vector< TVector3 > RefMom, vector< CbmMCTrack * > MC, vector< Int_t > Id, vector< TH1 * > gg)
vector< int > All_Pion_Id
TH2D * rap_vs_pt_NOTeta_gg_reco_aftercuts
TH1D * InvMass_eta_ppg_reco_aftercuts
TH1D * InvMass_eta_ppp_reffited
std::vector< TVector3 > EMT_ppg_negative_pion_momenta
vector< CbmMCTrack * > All_Pion_MCtrack
std::vector< std::vector< TVector3 > > EDGA_RefMom
vector< TH1 * > fHistoList_eta_ppg
TH1D * OpeningAnglePhoton_eta_ppg_reffited
std::vector< int > EMT_ppg_ee_Event
vector< TH1 * > fHistoList_eta_ppp
vector< CbmMCTrack * > El_Photon_Eta_MCtrack
TH1D * OpeningAngle_eta_gg_between_gg_mc
TH1D * OA_betweenPions_fromEta_mc
TH1D * Pion_P_fromEta_reco
std::vector< int > EMT_ppg_negative_pion_Event
TH1D * InvMass_eta_ppg_reffited
vector< int > El_Photon_Pion_Eta_Id
vector< TVector3 > All_Pion_refmomentum
TH1D * Pion_P_elsewhere_reco
std::vector< std::vector< int > > ECPGA_leptons_Id
TH1D * InvMass_eta_ppp_allcombinations_mc
TH2D * rap_vs_pt_eta_ppg_reco_aftercuts
static TVector3 FitToVertexAndGetChi(CbmStsTrack *stsTrack, double x, double y, double z, double &chi)
static LmvmKinePar CalculateKinematicParams_2el_2pions(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
static Double_t CalculateOpeningAngle_Reco(TVector3 electron1, TVector3 electron2)
static double Invmass_2particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2)
static Double_t CalculateOpeningAngle_MC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2)
static Double_t CalculateOpeningAngleBetweenPions_MC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2)
calculate opening angle between two pions using MCtrue momenta
static Double_t CalculateOpeningAngleBetweenPions_Reco(TVector3 electron1, TVector3 electron2)
calculate opening angle between two pions using reconstructed momenta
static LmvmKinePar CalculateKinematicParams_4particles(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
static double Invmass_6particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4, const CbmMCTrack *mctrack5, const CbmMCTrack *mctrack6)
static double Invmass_2particles_RECO(const TVector3 part1, const TVector3 part2)
static double Invmass_4particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
static double Invmass_4particles_RECO(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
static double Invmass_4el_2pions_RECO(const TVector3 part1El, const TVector3 part2El, const TVector3 part3El, const TVector3 part4El, const TVector3 part5Pion, const TVector3 part6Pion)
static double Invmass_2el_2pions_RECO(const TVector3 part1El, const TVector3 part2El, const TVector3 part3Pion, const TVector3 part4Pion)
int32_t GetMotherId() const
int32_t GetPdgCode() const
void Get4Momentum(TLorentzVector &momentum) const
int32_t GetNofLinks() const
const CbmLink & GetMatchedLink() const