20 #include "FairRootManager.h"
21 #include "KFParticleTopoReconstructor.h"
22 #include "KFTopoPerformance.h"
23 #include "TDatabasePDG.h"
26 #define M2E 2.6112004954086e-7
32 : fKFMcParticles(NULL)
35 , fStsTrackMatches(NULL)
37 , fKFparticleFinderQA(NULL)
42 , particlecounter_2daughters(0)
43 , particlecounter_3daughters(0)
44 , particlecounter_4daughters(0)
45 , particlecounter_all(0)
46 , fhPi0_NDaughters(NULL)
47 , fNofGeneratedPi0_allEvents(0)
48 , fNofPi0_kfparticle_allEvents(0)
50 , fNofPi0_kfparticle(0)
55 , fHistoList_kfparticle()
60 , fhInvMassPi0WithFullReco(NULL)
61 , fhInvMass2Gammas(NULL)
62 , fhInvMass2Gammas_cut(NULL)
63 , fhKF_particlevector(NULL)
64 , fhKF_trackvector(NULL)
66 , fhKF_NofPi0_signal(NULL)
67 , fhKF_NofPi0_trackvector(NULL)
69 , fhKF_invmass_fullReco(NULL)
77 FairRootManager* ioman = FairRootManager::Instance();
79 Fatal(
"CbmAnaConversionKF::Init",
"RootManager not instantised!");
82 fKFMcParticles = (TClonesArray*) ioman->GetObject(
"KFMCParticles");
84 Fatal(
"CbmAnaConversionKF::Init",
"No KFMCParticles array!");
87 fMcTracks = (TClonesArray*) ioman->GetObject(
"MCTrack");
89 Fatal(
"CbmAnaConversionKF::Init",
"No MCTrack array!");
92 fStsTracks = (TClonesArray*) ioman->GetObject(
"StsTrack");
94 Fatal(
"CbmAnaConversionKF::Init",
"No StsTrack array!");
99 Fatal(
"CbmAnaConversionKF::Init",
"No StsTrackMatch array!");
117 "fhPi0_NDaughters",
"fhPi0_NDaughters;number of daughers;#", 4, 0.5, 4.5);
118 fhPi0Ratio =
new TH1D(
"fhPi0Ratio",
"fhPi0Ratio; ratio;#", 1000, 0., 0.02);
119 fhPi0_mass =
new TH1D(
"fhPi0_mass",
"fhPi0_mass;mass;#", 500, 0., 0.5);
125 "fhInvMassPi0WithFullReco;invmass;#",
132 new TH1D(
"fhInvMass2Gammas",
"fhInvMass2Gammas;invmass;#", 400, 0., 2.);
135 "fhInvMass2Gammas_cut",
"fhInvMass2Gammas_cut;invmass;#", 400, 0., 2.);
139 new TH1D(
"fhKF_particlevector",
"fhKF_particlevector;;#", 10, 0., 10.);
145 new TH1D(
"fhKF_trackvector",
"fhKF_trackvector;;#", 10, 0., 10.);
151 fhKF_NofPi0 =
new TH1D(
"fhKF_NofPi0",
"fhKF_NofPi0;nof;#", 10, -0.5, 9.5);
154 new TH1D(
"fhKF_NofPi0_signal",
"fhKF_NofPi0_signal;nof;#", 10, -0.5, 9.5);
157 "fhKF_NofPi0_trackvector",
"fhKF_NofPi0_trackvector;nof;#", 10, -0.5, 9.5);
162 "fhKF_invmass_fullReco",
"fhKF_invmass_fullReco;invmass;#", 400, 0., 2.);
169 gDirectory->mkdir(
"KF");
170 gDirectory->cd(
"KF");
174 gDirectory->cd(
"..");
176 cout <<
"CbmAnaConversionKF: Realtime - " <<
fTime << endl;
206 cout <<
"CbmAnaConversionKF: kf works" << endl;
208 cout <<
"CbmAnaConversionKF: kf works not" << endl;
385 cout <<
"CbmAnaConversionKF: test nof " << nofparticles << endl;
386 for (
int i = 0;
i < nofparticles;
i++) {
389 if (pdg == 111) { cout <<
"CbmAnaConversionKF: test successful!" << endl; }
393 vector<vector<int>> ids;
394 const vector<KFParticle>& particles =
396 for (
unsigned int iPart = 0; iPart <
fSignalIds.size(); iPart++) {
397 if (particles[
fSignalIds[iPart]].GetPDG() != 111)
continue;
400 const KFParticle& pi0 = particles[
fSignalIds[iPart]];
401 vector<int> electrons;
402 for (
int iGamma = 0; iGamma < pi0.NDaughters(); iGamma++) {
403 const int GammaID = pi0.DaughterIds()[iGamma];
404 const KFParticle& Gamma = particles[GammaID];
405 for (
int iElectron = 0; iElectron < Gamma.NDaughters(); iElectron++) {
406 int ElectronID = Gamma.DaughterIds()[iElectron];
407 const KFParticle& Electron = particles[ElectronID];
408 int STStrackID = Electron.DaughterIds()[0];
409 electrons.push_back(STStrackID);
412 ids.push_back(electrons);
415 if (ids.size() > 0) {
416 cout <<
"NEW TEST: (sts ids) ";
417 for (
unsigned int i = 0;
i < ids.size();
i++) {
418 for (
int j = 0; j < 4; j++) {
419 cout <<
" " << ids[
i][j];
428 for (
unsigned int i = 0;
i < ids.size();
i++) {
429 for (
int j = 0; j < 4; j++) {
431 if (stsTrack == NULL)
return;
434 if (stsMatch == NULL)
return;
436 if (stsMcTrackId < 0)
return;
446 cout <<
" " << stsMcTrackId <<
"/" <<
mcTracks[j]->GetPdgCode()
447 <<
"(motherid" <<
mcTracks[j]->GetMotherId() <<
",motherpdg"
449 << grandmother->
GetPdgCode() <<
",grandmotherid"
456 cout <<
"mass: " << mass << endl;
467 TLorentzVector lorVec1;
470 TLorentzVector lorVec2;
473 TLorentzVector lorVec3;
476 TLorentzVector lorVec4;
481 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
482 cout <<
"mc: \t" << sum.Px() <<
" / " << sum.Py() <<
" / " << sum.Pz()
483 <<
" / " << sum.E() <<
"\t => mag = " << sum.Mag() << endl;
490 Int_t pi0counter_trackvector = 0;
491 const KFPTrackVector* kftrackvector;
492 kftrackvector =
fKFtopo->GetTracks();
493 Int_t kftv_size = kftrackvector->Size();
494 cout <<
"CbmAnaConversionKF: size of kftrackvector: " << kftv_size << endl;
496 kfvector_int vector_pdgs = kftrackvector->PDG();
498 cout <<
"CbmAnaConversionKF: pdgs in trackvector: ";
499 for (
int i = 0;
i < kftv_size;
i++) {
500 cout << vector_pdgs[
i] <<
" / ";
501 if (TMath::Abs(vector_pdgs[
i]) == 111) {
502 pi0counter_trackvector++;
524 cout <<
"CbmAnaConversionKF: size of particlevector: " << pv_size << endl;
525 cout <<
"CbmAnaConversionKF: size of particleMatch: " <<
particleMatch.size()
528 Int_t electroncounter = 0;
529 Int_t pi0counter = 0;
530 Int_t pi0counter_signal = 0;
531 Int_t gammacounter = 0;
532 for (
int i = 0;
i < pv_size;
i++) {
556 cout <<
"CbmAnaConversionKF: nof electrons in particlevector: "
557 << electroncounter << endl;
558 cout <<
"CbmAnaConversionKF: nof pi0 in particlevector: " << pi0counter
560 cout <<
"CbmAnaConversionKF: nof gamma in particlevector: " << gammacounter
568 for (
int a = 0; a < nof - 1; a++) {
569 for (
int b = a + 1; b < nof; b++) {
572 Double_t openingAngle = electron1.GetAngle(electron2);
574 TLorentzVector lorentzE1;
575 lorentzE1.SetPxPyPzE(electron1.GetPx(),
579 TLorentzVector lorentzE2;
580 lorentzE2.SetPxPyPzE(electron2.GetPx(),
584 TLorentzVector sum = lorentzE1 + lorentzE2;
585 Double_t invmass = sum.Mag();
587 if (openingAngle < 1 && invmass < 0.03) {
602 for (
int a = 0; a < nof - 1; a++) {
603 for (
int b = a + 1; b < nof; b++) {
623 for (
int a = 0; a < nof - 3; a++) {
624 for (
int b = a + 1; b < nof - 2; b++) {
625 for (
int c = b + 1; c < nof - 1; c++) {
626 for (
int d = c + 1;
d < nof;
d++) {
631 Int_t test1234 = check1 + check2 + check3 + check4;
641 particle1, particle2, particle3, particle4);
645 Double_t invmass_cut = 0.02;
647 if (particle1.GetZ() == particle2.GetZ()
648 && particle3.GetZ() == particle4.GetZ()) {
651 if (invmass_12 < invmass_cut && invmass_34 < invmass_cut) {
655 if (particle1.GetZ() == particle3.GetZ()
656 && particle2.GetZ() == particle4.GetZ()) {
659 if (invmass_13 < invmass_cut && invmass_24 < invmass_cut) {
663 if (particle1.GetZ() == particle4.GetZ()
664 && particle2.GetZ() == particle3.GetZ()) {
667 if (invmass_14 < invmass_cut && invmass_23 < invmass_cut) {
720 for (
int a = 0; a < nof - 1; a++) {
721 for (
int b = a + 1; b < nof; b++) {
725 Double_t openingAngle =
730 if ((TMath::Abs(particle1.GetZ() - particle2.GetZ()) < 0.005)
731 && (openingAngle < 5)) {
746 TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
747 Double_t energy1 = TMath::Sqrt(momentum1.Mag2() +
M2E);
748 TLorentzVector lorVec1(momentum1, energy1);
750 TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
751 Double_t energy2 = TMath::Sqrt(momentum2.Mag2() +
M2E);
752 TLorentzVector lorVec2(momentum2, energy2);
754 TVector3 momentum3(part3.GetPx(), part3.GetPy(), part3.GetPz());
755 Double_t energy3 = TMath::Sqrt(momentum3.Mag2() +
M2E);
756 TLorentzVector lorVec3(momentum3, energy3);
758 TVector3 momentum4(part4.GetPx(), part4.GetPy(), part4.GetPz());
759 Double_t energy4 = TMath::Sqrt(momentum4.Mag2() +
M2E);
760 TLorentzVector lorVec4(momentum4, energy4);
764 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
772 TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
773 Double_t energy1 = TMath::Sqrt(momentum1.Mag2());
774 TLorentzVector lorVec1(momentum1, energy1);
776 TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
777 Double_t energy2 = TMath::Sqrt(momentum2.Mag2());
778 TLorentzVector lorVec2(momentum2, energy2);
782 sum = lorVec1 + lorVec2;
790 TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
791 Double_t energy1 = TMath::Sqrt(momentum1.Mag2() +
M2E);
792 TLorentzVector lorVec1(momentum1, energy1);
794 TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
795 Double_t energy2 = TMath::Sqrt(momentum2.Mag2() +
M2E);
796 TLorentzVector lorVec2(momentum2, energy2);
800 sum = lorVec1 + lorVec2;
808 TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
809 Double_t energy1 = TMath::Sqrt(momentum1.Mag2());
810 TLorentzVector lorVec1(momentum1, energy1);
812 TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
813 Double_t energy2 = TMath::Sqrt(momentum2.Mag2());
814 TLorentzVector lorVec2(momentum2, energy2);
817 Double_t angleBetweenPhotons = lorVec1.Angle(lorVec2.Vect());
818 Double_t theta = 180. * angleBetweenPhotons /
TMath::Pi();