From fefe788042d36d1b76e324eb9d24ecffdb7e698d Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Mon, 27 Jul 2020 16:38:54 +0200
Subject: [PATCH] Advancing on strasse missing mass
---
NPLib/Physics/NPQFS.cxx | 5 +--
NPLib/Physics/TReactionConditions.h | 9 ++++-
NPSimulation/Process/BeamReaction.cc | 10 ++---
Projects/Strasse/Analysis.cxx | 43 +++++++++++++---------
Projects/Strasse/Analysis.h | 2 +-
Projects/Strasse/reaction/C12_p2p.reaction | 10 ++---
6 files changed, 45 insertions(+), 34 deletions(-)
diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 663cc53e0..0eca56b43 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -215,7 +215,7 @@ void QFS::CalculateVariables(){
//cout<<"---- COMPUTE ------"<<endl;
// cout<<"--CM--"<<endl;
- mA = fNucleiA.Mass(); // Beam mass in MeV
+ mA = fNucleiA.Mass(); // Beam mass in MeV
mT = fNucleiT.Mass(); // Target mass in MeV
mB = fNucleiB.Mass(); // Heavy residual mass in MeV
m1 = mT; // scattered target nucleon (same mass);
@@ -307,7 +307,7 @@ void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEn
pCM_2*sin(thetaCM_2)*sin(phiCM_2),
pCM_2*cos(thetaCM_2),
ECM_2);
-
+
fEnergyImpulsionCM_1 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_2;
//-- Boost in the direction of the moving cluster "a" --//
@@ -346,7 +346,6 @@ void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEn
// Kinetic Energy in the lab frame
KineticEnergyLab1 = fEnergyImpulsionLab_1.E() - m1;
KineticEnergyLab2 = fEnergyImpulsionLab_2.E() - m2;
-
// test for total energy conversion
//if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_1.E()+fEnergyImpulsionLab_2.E())) > 1e-6)
// cout << "Problem for energy conservation" << endl;
diff --git a/NPLib/Physics/TReactionConditions.h b/NPLib/Physics/TReactionConditions.h
index f4e5c613f..93c70537e 100644
--- a/NPLib/Physics/TReactionConditions.h
+++ b/NPLib/Physics/TReactionConditions.h
@@ -104,7 +104,12 @@ public:
void SetMomentumDirectionX (const double & Momentum_Direction_X) {fRC_Momentum_Direction_X.push_back(Momentum_Direction_X);}//!
void SetMomentumDirectionY (const double & Momentum_Direction_Y) {fRC_Momentum_Direction_Y.push_back(Momentum_Direction_Y);}//!
void SetMomentumDirectionZ (const double & Momentum_Direction_Z) {fRC_Momentum_Direction_Z.push_back(Momentum_Direction_Z);}//!
- void SetMomentum (const TVector3 & Momentum) {fRC_Momentum.push_back(Momentum);}//!
+ void SetMomentum (const TVector3 & Momentum) {
+ Momentum.Unit();
+ SetMomentumDirectionX(Momentum.X());
+ SetMomentumDirectionY(Momentum.Y());
+ SetMomentumDirectionZ(Momentum.Z());
+ }//!
///////////////////// GETTERS ////////////////////////
// Beam parameter
@@ -133,7 +138,7 @@ public:
double GetMomentumDirectionX (const int &i) const {return fRC_Momentum_Direction_X[i];}//!
double GetMomentumDirectionY (const int &i) const {return fRC_Momentum_Direction_Y[i];}//!
double GetMomentumDirectionZ (const int &i) const {return fRC_Momentum_Direction_Z[i];}//!
- TVector3 GetParticleMomentum (const int &i) const {return fRC_Momentum[i];}//!
+ TVector3 GetParticleMomentum (const int &i) const {return TVector3(fRC_Momentum_Direction_X[i],fRC_Momentum_Direction_Y[i],fRC_Momentum_Direction_Z[i]).Unit();}//!
TVector3 GetBeamDirection () const ;
TVector3 GetParticleDirection (const int i) const ;
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 909451f77..2a5903217 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -482,10 +482,10 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
TLorentzVector* P_A = m_QFS.GetEnergyImpulsionLab_A();
TLorentzVector* P_B = m_QFS.GetEnergyImpulsionLab_B();
-
+
G4ThreeVector momentum_kineB_beam( P_B->Px(), P_B->Py(), P_B->Pz() );
momentum_kineB_beam = momentum_kineB_beam.unit();
- TKEB = m_QFS.GetEnergyImpulsionLab_B()->Energy() - m_QFS.GetNucleusB()->Mass();
+ TKEB = P_B->Energy() - m_QFS.GetNucleusB()->Mass();
G4ThreeVector momentum_kineB_world = momentum_kineB_beam;
momentum_kineB_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
momentum_kineB_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
@@ -549,15 +549,15 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_ReactionConditions->SetInternalMomentum(m_QFS.GetInternalMomentum());
//m_ReactionConditions->SetExcitationEnergy3(m_QFS.GetExcitation3());
//m_ReactionConditions->SetExcitationEnergy4(m_QFS.GetExcitation4());
- // Momuntum X 3 and 4 //
+ // Momuntum X 1,2 and B //
m_ReactionConditions->SetMomentumDirectionX(momentum_kine1_world.x());
m_ReactionConditions->SetMomentumDirectionX(momentum_kine2_world.x());
m_ReactionConditions->SetMomentumDirectionX(momentum_kineB_world.x());
- // Momuntum Y 3 and 4 //
+ // Momuntum Y 1,2 and B //
m_ReactionConditions->SetMomentumDirectionY(momentum_kine1_world.y());
m_ReactionConditions->SetMomentumDirectionY(momentum_kine2_world.y());
m_ReactionConditions->SetMomentumDirectionY(momentum_kineB_world.y());
- // Momuntum Z 3 and 4 //
+ // Momuntum Z 1,2 and B //
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine1_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine2_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kineB_world.z());
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index d22b19e08..92dc2e844 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -27,6 +27,7 @@ using namespace std;
#include"NPOptionManager.h"
#include"NPFunction.h"
#include"NPTrackingUtility.h"
+#include"NPPhysicalConstants.h"
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -40,19 +41,18 @@ void Analysis::Init(){
DC= new TInteractionCoordinates;
RC= new TReactionConditions;
-
InitOutputBranch();
InitInputBranch();
Strasse = (TStrassePhysics*) m_DetectorManager -> GetDetector("Strasse");
Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana");
// reaction properties
- myQFS = new NPL::QFS();
- myQFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ m_QFS = new NPL::QFS();
+ m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
// reaction properties
myBeam = new NPL::Beam();
myBeam->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
- InitialBeamEnergy = myBeam->GetEnergy() * myBeam->GetA();
+ InitialBeamEnergy = myBeam->GetEnergy()/* * myBeam->GetA()*/;
// target thickness
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -62,8 +62,7 @@ void Analysis::Init(){
protonTarget = NPL::EnergyLoss("proton_"+TargetMaterial+".G4table","G4Table",100);
protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100);
protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100);
-
-
+ LV_T.SetVectM(TVector3(0,0,0),NPUNITS::proton_mass_c2);
}
////////////////////////////////////////////////////////////////////////////////
@@ -109,20 +108,28 @@ void Analysis::TreatEvent(){
if(i1!=i2){
E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]);
E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]);
+ double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),0);
+ //double TA = RC->GetBeamEnergy();
+ // setting up Lorentz Vector from measured trajectories and energies
+ TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetNucleusA()->Mass()))); // for like there is no BDC
+ Proton1=E1*Proton1.Unit();
+ Proton2=E2*Proton2.Unit();
+
+ LV_A.SetVectM(PA,m_QFS->GetNucleusA()->Mass());
+ double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2));
+ double P2= sqrt(E2*(E2+2*NPUNITS::proton_mass_c2));
+
+ LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2);
+ LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2);
+
+ // computing Ex from Missing Mass
+ LV_B = LV_A + LV_T - LV_p1 - LV_p2;
+ //LV_B = RC->GetParticleMomentum(2);
+ Ex = LV_B.M() - m_QFS->GetNucleusB()->Mass();
+ cout << Ex << " " << m_QFS->GetNucleusB()->Mass() << endl;
}
+
}
- //double thickness_before = 0;
- //double EA_vertex = BeamTarget.Slow(InitialBeamEnergy,thickness_before,0);
-
- // setting up Lorentz Vector from measured trajectories and energies
- //LV_A.SetVect(PA); LV_p1.SetE(EA_vertex);
- //LV_p1.SetVect(P1); LV_p1.SetE(E1);
- //LV_p2.SetVect(P2); LV_p1.SetE(E2);
-
- // computing Ex from Missing Mass
- //double EB = LV_A.E() + LV_T.E() - LV_p1.E() - LV_p2.E();
- //TVector3 PB = LV_A.Vect() + LV_p1.Vect() - LV_p2.Vect();
- //Ex = TMath::Sqrt( EB*EB - PB.Mag2() ) - myQFS->GetNucleusB()->Mass();
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h
index 08494fc00..52cf9bea4 100644
--- a/Projects/Strasse/Analysis.h
+++ b/Projects/Strasse/Analysis.h
@@ -73,7 +73,7 @@ class Analysis: public NPL::VAnalysis{
TLorentzVector LV_p2;
NPL::Beam* myBeam;
- NPL::QFS* myQFS;
+ NPL::QFS* m_QFS;
// Energy loss table: the G4Table are generated by the simulation
EnergyLoss BeamTarget;
EnergyLoss protonTarget;
diff --git a/Projects/Strasse/reaction/C12_p2p.reaction b/Projects/Strasse/reaction/C12_p2p.reaction
index 47d81151a..97a681ccc 100755
--- a/Projects/Strasse/reaction/C12_p2p.reaction
+++ b/Projects/Strasse/reaction/C12_p2p.reaction
@@ -4,10 +4,10 @@ Beam
Particle= 12C
Energy= 4800 MeV
SigmaEnergy= 0 MeV
- SigmaThetaX= 0.1 deg
- SigmaPhiY= 0.1 deg
- SigmaX= 5 mm
- SigmaY= 5 mm
+ SigmaThetaX= 0 deg
+ SigmaPhiY= 0 deg
+ SigmaX= 0 mm
+ SigmaY= 0 mm
MeanThetaX= 0 deg
MeanPhiY= 0 deg
MeanX= 0 mm
@@ -22,7 +22,7 @@ QFSReaction
Heavy= 11B
ExcitationEnergyBeam= 0.0 MeV
ExcitationEnergyHeavy= 0.0 MeV
- MomentumSigma= 50.0
+ MomentumSigma= 0.0
ShootHeavy= 1
ShootLight= 1
--
GitLab