diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index f7c7551413e47cd51ee2955f29fe38101eca41a9..f7277bdd7df1fd7588b7913a89114206cfe01dec 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -233,19 +233,11 @@ void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, doubl
fXThetaXHist->GetRandom2(X,ThetaX);
fYPhiYHist->GetRandom2(Y,PhiY);
}
- // Direction
- double Xdir = sin(ThetaX); // cos(90-x) = sin(x)
- double Ydir = sin(PhiY);
- double Zdir = sqrt(1-Xdir*Xdir-Ydir*Ydir); // alpha^2 + beta^2 + gamma^2 = 1
- // Stretch factor to extend unitary vector from ZEmission to ZProfile
- double S = fZProfile-fZEmission ;
- TVector3 BeamDir(Xdir*S/Zdir,Ydir*S/Zdir,S);
-
- Xdir = BeamDir.X();
- Ydir = BeamDir.Y();
- // POSITION/DIRECTION AT Z EMISSION//
- X = X0 - Xdir;
- Y = Y0 - Ydir;
+ // Backward propagtion from ZProfile (in general on target)
+ // to ZEmission (position where beam is generated in simulation.)
+ double dZ = fZProfile-fZEmission ; //(Z1-Z0)
+ X= X0-tan(ThetaX)*dZ;
+ Y= Y0-tan(PhiY)*dZ;
Z = fZEmission;
}
diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index b19851ca6e57f04a99f0d9bd7eae04f764df1235..587919787ea6cbb1859a616b0989b37aabe9ad31 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -78,7 +78,11 @@ QFS::QFS(){
fshootB=false;
fshoot1=true;
fshoot2=true;
- fisotropic = true;
+ fIsotropic = true;
+
+ fCondEcmPos = false;
+ fCondOffshellMassPos = false;
+ fIsAllowed = false;
fUseExInGeant4=true;
@@ -220,6 +224,9 @@ Particle QFS::GetParticle(string name, NPL::InputParser parser){
////////////////////////////////////////////////////////////////////////////////////////////
void QFS::CalculateVariables(){
+ fCondEcmPos = false;
+ fCondOffshellMassPos = false;
+
if(fBeamEnergy < 0)
fBeamEnergy = 0 ;
@@ -250,7 +257,8 @@ void QFS::CalculateVariables(){
// Off-shell mass of the bound nucleon from E conservation
// in virtual dissociation of A -> B + a
double buffer = mA*mA + mB*mB - 2*mA*sqrt(mB*mB+Pa.Mag2()) ;
- if(buffer<=0) { cout<<"ERROR off shell mass ma_off=\t"<<buffer<<endl; return;}
+ if(buffer>0) {fCondOffshellMassPos = true;}
+ else{/*cout<<"ERROR off shell mass ma_off=\t"<<buffer<<endl; Dump();*/ return;}
ma_off = sqrt(buffer);
//deduced total energies of "a" and "B" in restframe of A
@@ -279,7 +287,8 @@ void QFS::CalculateVariables(){
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) {fCondEcmPos = true;}
+ if(fEcm<=0) { /*cout<<"ERROR Ecm negative =\t"<<fEcm<<endl;Dump();*/ return;}
ECM_a = (s + ma_off*ma_off - mT*mT)/(2*sqrt(s));
ECM_T = (s + mT*mT - ma_off*ma_off)/(2*sqrt(s));
@@ -298,8 +307,11 @@ void QFS::CalculateVariables(){
void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEnergyLab1, double& ThetaLab2, double& PhiLab2, double& KineticEnergyLab2){
+ fIsAllowed=false;
CalculateVariables();
-
+ if(fCondOffshellMassPos && fCondEcmPos) fIsAllowed= true;
+ else return;
+
double thetaCM_1 = fThetaCM;
double thetaCM_2 = M_PI - thetaCM_1;
double phiCM_2 = fPhiCM;
@@ -501,17 +513,6 @@ TGraph* QFS::GetPhi2VsPhi1(double AngleStep_CM){
return(fPhi2VsPhi1);
}
-///////////////////////////////////////////////////////////////////////////////
-// Check whenever the reaction is allowed at a given energy
-bool QFS::IsAllowed(){//double Energy){
- //double AvailableEnergy = Energy + fParticle1.Mass() + fParticle2.Mass();
- //double RequiredEnergy = fParticle3.Mass() + fParticle4.Mass();
-
- //if(AvailableEnergy>RequiredEnergy)
- return true;
- //else
- // return false;
-}
////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index 0afde4ccda55ae51395e2116f19b584d1c115b59..c5954ef76ee1cd8eb3dfde21efaf66651ef3c4a1 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -83,7 +83,11 @@ namespace NPL{
TVector3 fInternalMomentum; // Internal momentum of the removed cluster
TH1F* fPerpMomentumHist; // Perpendicular momentum distribution in beam at rest frame
TH1F* fParMomentumHist; // Parallel momentum distribution in beam at rest frame
- double fisotropic;
+ double fIsotropic;
+
+ bool fCondEcmPos;
+ bool fCondOffshellMassPos;
+ bool fIsAllowed;
TGraph* fTheta2VsTheta1;
TGraph* fPhi2VsPhi1;
@@ -103,7 +107,6 @@ namespace NPL{
void TestR3B();
void KineRelativistic(double &ThetaLab1, double &PhiLab1, double &KineticEnergyLab1, double &ThetaLab2, double &PhiLab2, double &KineticEnergyLab2);
TVector3 ShootInternalMomentum();
- bool IsAllowed();
void Dump();
private: // intern precompute variable
@@ -169,6 +172,7 @@ namespace NPL{
void SetPhiCM(const double& angle) {fPhiCM = angle;}
void SetInternalMomentum(const TVector3& mom) {fInternalMomentum = mom;}
void SetMomentumSigma(const double& sigma) {fMomentumSigma = sigma;}
+ void SetIsAllowed(const bool& val) {fIsAllowed = val;}
void SetPerpMomentumHist (TH1F* PerpMomentumHist)
{delete fPerpMomentumHist; fPerpMomentumHist = PerpMomentumHist;}
void SetParMomentumHist (TH1F* ParMomentumHist)
@@ -187,7 +191,8 @@ namespace NPL{
double GetThetaCM() const {return fThetaCM;}
double GetPhiCM() const {return fPhiCM;}
double GetMomentumSigma() const {return fMomentumSigma;}
- bool GetUseExInGeant4() const { return fUseExInGeant4; }
+ bool IsAllowed() const {return fIsAllowed;}
+ bool GetUseExInGeant4() const { return fUseExInGeant4; }
double GetExcitationA() const {return fExcitationA;}
double GetExcitationB() const {return fExcitationB;}
TVector3 GetInternalMomentum() const {return fInternalMomentum;}
diff --git a/NPSimulation/EventGenerator/EventGeneratorBeam.cc b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
index d169d083067a2a4f8f3c1217c4e8a81957db65c4..4e33ce3cf296bcfc2ae0d443f2409589355d3ae9 100644
--- a/NPSimulation/EventGenerator/EventGeneratorBeam.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
@@ -97,12 +97,13 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
double InitialBeamEnergy, x0, y0, z0, Beam_thetaX, Beam_phiY;
m_Beam->GenerateRandomEvent(InitialBeamEnergy, x0, y0, z0, Beam_thetaX, Beam_phiY);
- //Set the direction cosines: alpha=90-Beam_thetaX, beta=90-Beam_phiY
- double Xdir = sin(Beam_thetaX); // cos(90-x) = sin(x)
- double Ydir = sin(Beam_phiY);
- double Zdir = sqrt(1-Xdir*Xdir-Ydir*Ydir); // alpha^2 + beta^2 + gamma^2 = 1
+ double Xdir = tan(Beam_thetaX); //tan(thetax)= px/pz
+ double Ydir = tan(Beam_phiY); //tan(phiy)= py/pz
+ double Zdir = 1; // fix pz=1 arbitrarily
G4ThreeVector BeamDir(Xdir,Ydir,Zdir);
+ BeamDir = BeamDir.unit();
G4ThreeVector BeamPos(x0,y0,z0);
+
///////////////////////////////////////////////////////
///// Add the Beam particle to the particle Stack /////
///////////////////////////////////////////////////////
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 189d397dc3e91dbeb312c14bdf742b88424c6673..1a8ba7f3edb150c20c14cf7562837c173ced425d 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -183,13 +183,13 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
// If the condition is met, the event is generated
if (m_shoot && m_length < m_StepSize) {
if(m_ReactionType==QFS){
- if ( m_QFS.IsAllowed() ) {
+ //if ( m_QFS.IsAllowed() ) {
return true;
- }
- else{
- m_shoot=false;
- std::cout << "QFS not allowed" << std::endl;
- }
+ //}
+ //else{
+ // m_shoot=false;
+ // std::cout << "QFS not allowed" << std::endl;
+ //}
}
else if(m_ReactionType==TwoBody){
@@ -475,14 +475,6 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
// Set the Energy of the reaction
m_QFS.SetBeamEnergy(reac_energy);
- double Beam_theta = pdirection.theta();
- double Beam_phi = pdirection.phi();
-
-
- /////////////////////////////////////////////////////////////////
- ///// Angles for emitted particles following Cross Section //////
- ///// Angles are in the beam frame //////
- /////////////////////////////////////////////////////////////////
// Shoot and Set a Random ThetaCM
double costheta = G4RandFlat::shoot() * 2 - 1;
@@ -491,45 +483,55 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_QFS.SetThetaCM(theta);
m_QFS.SetPhiCM(phi);
- // Shoot and set internal momentum for the removed cluster
- // if a momentum Sigma is given then shoot in 3 indep. Gaussians
- // if input files are given for distributions use them instead
- m_QFS.ShootInternalMomentum();
-
- //////////////////////////////////////////////////
- ///// Momentum and angles from kinematics /////
- //////////////////////////////////////////////////
- // Variable where to store results
+ // Lab frame variables where to store results
double Theta1, Phi1, TKE1, Theta2, Phi2, TKE2, ThetaB, PhiB, TKEB;
- // Compute Kinematic using previously defined ThetaCM
- m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
+ int j=0;
+ m_QFS.SetIsAllowed(false);
+ while(!m_QFS.IsAllowed()){
+ // Shoot internal momentum for the removed cluster
+ // if a momentum Sigma is given then shoot in 3 indep. Gaussians
+ // if input files are given for distributions use them instead
+ m_QFS.ShootInternalMomentum();
- G4ThreeVector U(1, 0, 0);
- G4ThreeVector V(0, 1, 0);
- G4ThreeVector ZZ(0, 0, 1);
+ // Go from CM to Lab
+ m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
+
+ j++;
+ if(j>100) cout<< "ERROR: too many iteration and QFS kinematical conditions not allowed"<<endl;
+ }
+
+ //---------------------------------------------------------
+ //Rotations to switch from frame with beam along Z to world
+ //---------------------------------------------------------
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+ G4ThreeVector ux(1, 0, 0);
+ G4ThreeVector uy(0, 1, 0);
+ G4ThreeVector uz(0, 0, 1);
// Momentum in beam and world frame for light particle 1
G4ThreeVector momentum_kine1_beam(sin(Theta1) * cos(Phi1),
sin(Theta1) * sin(Phi1), cos(Theta1));
G4ThreeVector momentum_kine1_world = momentum_kine1_beam;
- momentum_kine1_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
- momentum_kine1_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+ momentum_kine1_world.rotate(Beam_theta, uy); // rotation of Beam_theta around Y axis
+ momentum_kine1_world.rotate(Beam_phi, uz); // rotation of Beam_phi around Z axis
// Momentum in beam and world frame for light particle 2
G4ThreeVector momentum_kine2_beam(sin(Theta2) * cos(Phi2),
sin(Theta2) * sin(Phi2), cos(Theta2));
G4ThreeVector momentum_kine2_world = momentum_kine2_beam;
- momentum_kine2_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
- momentum_kine2_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+ momentum_kine2_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ momentum_kine2_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
// Momentum in beam and world frame for heavy residual
//
//G4ThreeVector momentum_kineB_beam(sin(ThetaB) * cos(PhiB + pi),
// sin(ThetaB) * sin(PhiB + pi), cos(ThetaB));
//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
+ //momentum_kineB_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ //momentum_kineB_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
TLorentzVector* P_A = m_QFS.GetEnergyImpulsionLab_A();
TLorentzVector* P_B = m_QFS.GetEnergyImpulsionLab_B();
@@ -538,12 +540,13 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
momentum_kineB_beam = momentum_kineB_beam.unit();
TKEB = P_B->Energy() - m_QFS.GetParticleB()->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
+ momentum_kineB_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ momentum_kineB_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
ThetaB = P_B->Angle(P_A->Vect());
- if (ThetaB < 0) ThetaB += M_PI;
- PhiB = M_PI + P_B->Vect().Phi();
+ //if (ThetaB < 0) ThetaB += M_PI;
+ //PhiB = M_PI + P_B->Vect().Phi();
+ PhiB = P_B->Vect().Phi();
if (fabs(PhiB) < 1e-6) PhiB = 0;
@@ -602,44 +605,45 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_ReactionConditions->SetVertexPositionX(localPosition.x());
m_ReactionConditions->SetVertexPositionY(localPosition.y());
m_ReactionConditions->SetVertexPositionZ(localPosition.z());
+
m_ReactionConditions->SetBeamEmittanceTheta(
PrimaryTrack->GetMomentumDirection().theta() / deg);
m_ReactionConditions->SetBeamEmittancePhi(
PrimaryTrack->GetMomentumDirection().phi() / deg);
m_ReactionConditions->SetBeamEmittanceThetaX(
- PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+ PrimaryTrack->GetMomentumDirection().perpPart(uy).angle(uz) / deg);
m_ReactionConditions->SetBeamEmittancePhiY(
- PrimaryTrack->GetMomentumDirection().angle(V) / deg);
+ PrimaryTrack->GetMomentumDirection().perpPart(ux).angle(uz) / deg);
// Names 1,2 and B//
m_ReactionConditions->SetParticleName(Light1Name->GetParticleName());
m_ReactionConditions->SetParticleName(Light2Name->GetParticleName());
m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
- // Angle 1,2 and B //
+ // Angle 1,2 and B in fram with beam along Z axis//
m_ReactionConditions->SetTheta(Theta1 / deg);
m_ReactionConditions->SetTheta(Theta2 / deg);
m_ReactionConditions->SetTheta(ThetaB / deg);
m_ReactionConditions->SetPhi(Phi1 / deg);
m_ReactionConditions->SetPhi(Phi2 / deg);
m_ReactionConditions->SetPhi(PhiB / deg);
- // Energy 1,2 and B //
+ // Total Kinetic Energy 1,2 and B //
m_ReactionConditions->SetKineticEnergy(TKE1);
m_ReactionConditions->SetKineticEnergy(TKE2);
m_ReactionConditions->SetKineticEnergy(TKEB);
// 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());
+ //Excitation energy of reaction product B
m_ReactionConditions->SetExcitationEnergy4(m_QFS.GetExcitationB());
- // Momuntum X 1,2 and B //
+ // Momentum Dir. X 1,2 and B in world frame //
m_ReactionConditions->SetMomentumDirectionX(momentum_kine1_world.x());
m_ReactionConditions->SetMomentumDirectionX(momentum_kine2_world.x());
m_ReactionConditions->SetMomentumDirectionX(momentum_kineB_world.x());
- // Momuntum Y 1,2 and B //
+ // Momentum Dir. Y 1,2 and B in world frame//
m_ReactionConditions->SetMomentumDirectionY(momentum_kine1_world.y());
m_ReactionConditions->SetMomentumDirectionY(momentum_kine2_world.y());
m_ReactionConditions->SetMomentumDirectionY(momentum_kineB_world.y());
- // Momuntum Z 1,2 and B //
+ // Momentum Dir. Z 1,2 and B in world frame//
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine1_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine2_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kineB_world.z());