Skip to content
Snippets Groups Projects
Commit 6d7e4f28 authored by Morfouace's avatar Morfouace
Browse files

Merge branch 'NPTool.2.dev' of https://gitlab.in2p3.fr/np/nptool into NPTool.2.dev

parents 5ee5b8a3 e26a4f52
No related branches found
No related tags found
No related merge requests found
Pipeline #70545 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Examples/Example5/CheckSimu.C
+ 39
− 2
View file @ 6d7e4f28
......@@ -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");
}
NPLib/Physics/NPQFS.cxx
+ 14
− 16
View file @ 6d7e4f28
......@@ -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(0.*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(0.*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]);
......
NPLib/Physics/NPQFS.h
+ 5
− 1
View file @ 6d7e4f28
......@@ -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;}
......
NPLib/Physics/TReactionConditions.cxx
+ 6
− 0
View file @ 6d7e4f28
......@@ -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();
......
NPLib/Physics/TReactionConditions.h
+ 3
− 0
View file @ 6d7e4f28
......@@ -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_Internal_Momentum;}//!
// emmitted particles
int GetParticleMultiplicity() const {return fRC_Kinetic_Energy.size();}//!
......
NPSimulation/Detectors/PISTA/PISTA.hh
+ 1
− 1
View file @ 6d7e4f28
......@@ -103,7 +103,7 @@ class PISTA : public NPS::VDetector{
vector<double> m_Phi;
// Visualisation Attribute
G4VisAttributes* m_VisTrap;
//G4VisAttributes* m_VisTrap;
// Needed for dynamic loading of the library
public:
......
NPSimulation/Process/BeamReaction.cc
+ 7
− 1
View file @ 6d7e4f28
......@@ -455,9 +455,14 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
//m_QFS.ShootRandomPhiCM();
double theta = G4RandFlat::shoot() * pi;
double phi = G4RandFlat::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 //
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment