From b4cc6749f124e5a1d99f1897c7c6d6e8d3d55267 Mon Sep 17 00:00:00 2001
From: fysquare <freddy.flavigny@laposte.net>
Date: Thu, 14 May 2020 15:01:25 +0200
Subject: [PATCH] For QFS: Add Internal Momentum branch in SimulatedTree
---
Examples/Example5/CheckSimu.C | 41 +++++++++++++++++++++++++--
NPLib/Physics/NPQFS.cxx | 32 ++++++++++-----------
NPLib/Physics/NPQFS.h | 6 +++-
NPLib/Physics/TReactionConditions.cxx | 6 ++++
NPLib/Physics/TReactionConditions.h | 3 ++
NPSimulation/Process/BeamReaction.cc | 8 +++++-
6 files changed, 75 insertions(+), 21 deletions(-)
diff --git a/Examples/Example5/CheckSimu.C b/Examples/Example5/CheckSimu.C
index 237fe38a1..188924f87 100644
--- a/Examples/Example5/CheckSimu.C
+++ b/Examples/Example5/CheckSimu.C
@@ -58,6 +58,9 @@ void CheckSimu(const char * fname = "Example5"){
// Declare histograms
// for emitted particle
TH1F *hEmThetaCM = new TH1F("hEmThetaCM", "Light Ejectile Theta CM",180, 0, 180);
+ TH1F *hEmInternalMomX = new TH1F("hEmInternalMomX", "Internal Momentum (X) of removed cluster",500, -500, 500);
+ TH1F *hEmInternalMomY = new TH1F("hEmInternalMomY", "Internal Momentum (Y) of removed cluster",500, -500, 500);
+ TH1F *hEmInternalMomZ = new TH1F("hEmInternalMomZ", "Internal Momentum (Z) of removed cluster",500, -500, 500);
TH1F *hEmTheta1IF = new TH1F("hEmTheta1IF", "Ejectile1 Theta (reaction frame)", 180, 0, 180);
TH1F *hEmPhi1IF = new TH1F("hEmPhi1IF", "Ejectile1 Phi (reaction frame)", 360, 0, 360);
TH1F *hEmTheta1WF = new TH1F("hEmTheta1WF", "Ejectile1 Theta (world frame)", 180, 0, 180);
@@ -94,6 +97,9 @@ void CheckSimu(const char * fname = "Example5"){
// Fill histos
// ejected particles
hEmThetaCM -> Fill(reacCond->GetThetaCM());
+ hEmInternalMomX -> Fill(reacCond->GetInternalMomentum().X());
+ hEmInternalMomY -> Fill(reacCond->GetInternalMomentum().Y());
+ hEmInternalMomZ -> Fill(reacCond->GetInternalMomentum().Z());
hEmTheta1IF -> Fill(reacCond->GetThetaLab_BeamFrame(0));
hEmTheta1WF -> Fill(reacCond->GetThetaLab_WorldFrame(0));
@@ -122,8 +128,8 @@ void CheckSimu(const char * fname = "Example5"){
// Display emmitted paricles histograms
- canvas2 = new TCanvas("canvas2", "Emmited particles properties in reaction frame",1000,600);
- canvas2->Divide(3,2);
+ canvas2 = new TCanvas("canvas2", "Emmited particles properties in reaction frame",1000,1000);
+ canvas2->Divide(3,3);
canvas2->cd(1);
hEmThetaCM->SetXTitle("#theta_{c.m.}");
@@ -154,17 +160,48 @@ void CheckSimu(const char * fname = "Example5"){
Kine->SetLineColor(kRed);
Kine->SetLineStyle(2);
Kine->Draw("csame");
+ TGraph* Kine2 = qfs.GetTheta2VsTheta1(10);
+ Kine2->SetLineColor(kRed);
+ Kine2->SetMarkerColor(kRed);
+ Kine2->SetMarkerStyle(20);
+ Kine2->SetMarkerSize(1.3);
+ Kine2->Draw("Psame");
+
+
+
canvas2->cd(5);
hEmPhi1VsPhi2->Draw("colz");
hEmPhi1VsPhi2->SetXTitle("#phi_{1} (deg)");
hEmPhi1VsPhi2->SetYTitle("#phi_{2} (deg)");
+ TGraph* KinePhi = qfs.GetPhi2VsPhi1(1);
+ KinePhi->SetMarkerColor(kRed);
+ KinePhi->SetMarkerSize(0.4);
+ KinePhi->Draw("Psame");
+
canvas2->cd(6);
hEmOpAngle->Draw();
hEmOpAngle->SetXTitle("Opening angle (1-2) (deg)");
hEmOpAngle->SetYTitle("Counts");
+ canvas2->cd(7);
+ hEmInternalMomX->Draw();
+ hEmInternalMomX->SetXTitle("P_{x} (MeV/c)");
+ hEmInternalMomX->SetYTitle("Counts");
+
+ canvas2->cd(8);
+ hEmInternalMomY->Draw();
+ hEmInternalMomY->SetXTitle("P_{y} (MeV/c)");
+ hEmInternalMomY->SetYTitle("Counts");
+
+ canvas2->cd(9);
+ hEmInternalMomZ->Draw();
+ hEmInternalMomZ->SetXTitle("P_{z} (MeV/c)");
+ hEmInternalMomZ->SetYTitle("Counts");
+
+
+
}
diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index defdd5675..ed46fd18a 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -74,6 +74,7 @@ QFS::QFS(){
fExcitationA = 0;
fExcitationB = 0;
fMomentumSigma = 0;
+ fInternalMomentum = {0., 0.,0. };
fshootB=false;
fshoot1=true;
fshoot2=true;
@@ -215,11 +216,12 @@ void QFS::CalculateVariables(){
//Internal momentum of removed cluster/nucleon
//gRandom->SetSeed(0);
- //double mom_sigma = 0; // MeV/c
- Pa.SetX(gRandom->Gaus(0.,fMomentumSigma));
- Pa.SetY(gRandom->Gaus(0.,fMomentumSigma));
- Pa.SetZ(gRandom->Gaus(0.,fMomentumSigma));
-
+ //Pa.SetX(gRandom->Gaus(0.,fMomentumSigma));
+ //Pa.SetY(gRandom->Gaus(0.,fMomentumSigma));
+ //Pa.SetZ(gRandom->Gaus(0.,fMomentumSigma));
+ Pa.SetX(fInternalMomentum.X());
+ Pa.SetY(fInternalMomentum.Y());
+ Pa.SetZ(fInternalMomentum.Z());
//Internal momentum of heavy recoil after removal
PB.SetXYZ( (-Pa.X()) , (-Pa.Y()) , (-Pa.Z()) );
@@ -341,7 +343,6 @@ void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEn
////////////////////////////////////////////////////////////////////////////////////////////
void QFS::Dump(){
-
cout<<endl;
cout<<"------------------------------------"<<endl;
cout<<"------------ DUMP QFS --------------"<<endl;
@@ -396,8 +397,6 @@ void QFS::Dump(){
cout<<"Phi1:\t"<<fEnergyImpulsionLab_1.Vect().Phi()*180./TMath::Pi()<<endl;
cout<<"Phi2:\t"<<fEnergyImpulsionLab_2.Vect().Phi()*180./TMath::Pi()<<endl;
-
-
}
////////////////////////////////////////////////////////////////////////////////////////////
@@ -414,6 +413,7 @@ TGraph* QFS::GetTheta2VsTheta1(double AngleStep_CM){
vector<double> vx;
vector<double> vy;
double theta1,phi1,E1,theta2,phi2,E2;
+ SetPhiCM(20.*TMath::Pi()/180.);
for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
SetThetaCM(angle*TMath::Pi()/180.);
@@ -433,15 +433,13 @@ TGraph* QFS::GetPhi2VsPhi1(double AngleStep_CM){
vector<double> vx;
vector<double> vy;
double theta1,phi1,E1,theta2,phi2,E2;
+ SetThetaCM(20.*TMath::Pi()/180.);
- for (double theta=0 ; theta <= 180 ; theta+=AngleStep_CM){
- for (double angle=-180 ; angle <= 180 ; angle+=AngleStep_CM){
- SetThetaCM(theta*TMath::Pi()/180.);
- SetPhiCM(angle*TMath::Pi()/180.);
- KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
- vx.push_back(phi1*180./M_PI);
- vy.push_back(phi2*180./M_PI);
- }
+ for (double angle=-180 ; angle <= 180 ; angle+=AngleStep_CM){
+ SetPhiCM(angle*TMath::Pi()/180.);
+ KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
+ vx.push_back(phi1*180./M_PI);
+ vy.push_back(phi2*180./M_PI);
}
fPhi2VsPhi1 = new TGraph(vx.size(),&vx[0],&vy[0]);
@@ -538,7 +536,7 @@ void QFS::CalculateVariablesOld(){
s = ma_off*ma_off + mT*mT + 2*mT*Ea_lab ;
fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
fEcm = sqrt(s) - m1 -m2;
- if(fEcm<=0) { cout<<"ERROR Ecm negative =\t"<<fEcm<<endl;Dump(); return;}
+ if(fEcm<=0) { cout<<"ERROR Ecm negative =\t"<<fEcm<<endl;DumpFormatted(); return;}
vector<double> theta1;
vector<double> theta2;
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index 60af1e0ef..fde81d853 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -79,7 +79,8 @@ namespace NPL{
double fBeamEnergy; // Beam energy in MeV
double fExcitationA; // Excitation energy in MeV of the beam, useful for isomers
double fExcitationB; // Excitation energy in MeV of beam-like ejectile
- double fMomentumSigma; // Width of the momentum distribution of the removed cluster (sigma)
+ double fMomentumSigma; // Width of the momentum distribution (sigma)
+ TVector3 fInternalMomentum; // Internal momentum of the removed cluster
double fisotropic;
TGraph* fTheta2VsTheta1;
@@ -162,6 +163,7 @@ namespace NPL{
void SetBeamEnergy(const double& eBeam) {fBeamEnergy = eBeam;}
void SetThetaCM(const double& angle) {fThetaCM = angle;}
void SetPhiCM(const double& angle) {fPhiCM = angle;}
+ void SetInternalMomentum(const TVector3& mom) {fInternalMomentum = mom;}
void SetMomentumSigma(const double& sigma) {fMomentumSigma = sigma;}
//GETTERS
@@ -174,6 +176,8 @@ namespace NPL{
bool GetShoot2() const {return fshoot2;}
bool GetShootB() const {return fshootB;}
double GetThetaCM() const {return fThetaCM;}
+ double GetMomentumSigma() const {return fMomentumSigma;}
+ TVector3 GetInternalMomentum() const {return fInternalMomentum;}
TLorentzVector* GetEnergyImpulsionLab_A() {return &fEnergyImpulsionLab_A;}
TLorentzVector* GetEnergyImpulsionLab_T() {return &fEnergyImpulsionLab_T;}
diff --git a/NPLib/Physics/TReactionConditions.cxx b/NPLib/Physics/TReactionConditions.cxx
index 31e8a1d1d..bd312dc2d 100644
--- a/NPLib/Physics/TReactionConditions.cxx
+++ b/NPLib/Physics/TReactionConditions.cxx
@@ -49,6 +49,7 @@ void TReactionConditions::Clear(){
fRC_Vertex_Position_Z = -1;
fRC_ThetaCM = -1;
+ fRC_Internal_Momentum = {-1, -1; -1};
// emmitted particles
fRC_Particle_Name.clear();
@@ -81,6 +82,11 @@ void TReactionConditions::Dump() const{
<< fRC_Vertex_Position_X << " ; "
<< fRC_Vertex_Position_Y << " ; "
<< fRC_Vertex_Position_Z << ")" << endl;
+ cout << "\t If QFS, internal momentum, otherwise unused (-1) : ( "
+ << fRC_Internal_Momentum.X() << " ; "
+ << fRC_Internal_Momentum.Y() << " ; "
+ << fRC_Internal_Momentum.Z() << ")" << endl;
+
// emmitted particle
unsigned int size = fRC_Particle_Name.size();
diff --git a/NPLib/Physics/TReactionConditions.h b/NPLib/Physics/TReactionConditions.h
index 551c875ac..dc4a40436 100644
--- a/NPLib/Physics/TReactionConditions.h
+++ b/NPLib/Physics/TReactionConditions.h
@@ -59,6 +59,7 @@ private:
double fRC_ExcitationEnergy3;
double fRC_ExcitationEnergy4;
double fRC_ThetaCM;
+ TVector3 fRC_Internal_Momentum;
// emmitted particles
vector<string> fRC_Particle_Name;
vector<double> fRC_Theta;
@@ -92,6 +93,7 @@ public:
void SetExcitationEnergy3 (const double& Ex) {fRC_ExcitationEnergy3=Ex;};//!
void SetExcitationEnergy4 (const double& Ex) {fRC_ExcitationEnergy4=Ex;};//!
void SetThetaCM (const double & ThetaCM) {fRC_ThetaCM = ThetaCM;}//!
+ void SetInternalMomentum (const TVector3& IntMom) {fRC_Internal_Momentum = IntMom;}//!
// emmitted particles
void SetParticleName (const string & Particle_Name) {fRC_Particle_Name.push_back(Particle_Name);}//!
@@ -118,6 +120,7 @@ public:
double GetExcitation3 () const {return fRC_ExcitationEnergy3 ;}//!
double GetExcitation4 () const {return fRC_ExcitationEnergy4 ;}//!
double GetThetaCM () const {return fRC_ThetaCM;}//!
+ TVector3 GetInternalMomentum () const {return fRC_InternalMomentum;}//!
// emmitted particles
int GetParticleMultiplicity() const {return fRC_Kinetic_Energy.size();}//!
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 446ffc9d6..dcfe2ce82 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -455,9 +455,14 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
//m_QFS.ShootRandomPhiCM();
double theta = RandFlat::shoot() * pi;
double phi = RandFlat::shoot() * 2. * pi - pi; //rand in [-pi,pi]
+ double momX = gRandom->Gaus(0.,m_QFS.GetMomentumSigma());
+ double momY = gRandom->Gaus(0.,m_QFS.GetMomentumSigma());
+ double momZ = gRandom->Gaus(0.,m_QFS.GetMomentumSigma());
+ TVector3 mom(momX, momY, momZ);
m_QFS.SetThetaCM(theta);
m_QFS.SetPhiCM(phi);
+ m_QFS.SetInternalMomentum(mom);
//////////////////////////////////////////////////
///// Momentum and angles from kinematics /////
@@ -563,8 +568,9 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_ReactionConditions->SetKineticEnergy(TKE1);
m_ReactionConditions->SetKineticEnergy(TKE2);
m_ReactionConditions->SetKineticEnergy(TKEB);
- // ThetaCM and Ex//
+ // ThetaCM, Ex and Internal Momentum of removed cluster//
m_ReactionConditions->SetThetaCM(m_QFS.GetThetaCM() / deg);
+ m_ReactionConditions->SetInternalMomentum(m_QFS.GetInternalMomentum());
//m_ReactionConditions->SetExcitationEnergy3(m_QFS.GetExcitation3());
//m_ReactionConditions->SetExcitationEnergy4(m_QFS.GetExcitation4());
// Momuntum X 3 and 4 //
--
GitLab