From 52fa66214a8d980ea4b432f61f5975628a5a62d8 Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Thu, 17 Sep 2009 17:43:42 +0000
Subject: [PATCH] * Add the CalculateBeamInteraction() method in
VEventGenerator.{h,cxx}
---
.../gaspardTestSpheric.detector | 2 +-
Inputs/EventGenerator/132Sndp.reaction | 4 +-
NPSimulation/include/VEventGenerator.hh | 20 +--
NPSimulation/src/EventGeneratorBeam.cc | 72 ++++-------
NPSimulation/src/EventGeneratorTransfert.cc | 55 +++-----
.../src/EventGeneratorTransfertToResonance.cc | 53 +++-----
NPSimulation/src/VEventGenerator.cc | 118 +++++++++++++-----
7 files changed, 163 insertions(+), 161 deletions(-)
diff --git a/Inputs/DetectorConfiguration/gaspardTestSpheric.detector b/Inputs/DetectorConfiguration/gaspardTestSpheric.detector
index 6896d3de6..0b46aa52e 100644
--- a/Inputs/DetectorConfiguration/gaspardTestSpheric.detector
+++ b/Inputs/DetectorConfiguration/gaspardTestSpheric.detector
@@ -10,7 +10,7 @@ GeneralTarget
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Target
THICKNESS= 9.7
- ANGLE= 40
+ ANGLE= 0
RADIUS= 12
MATERIAL= CD2
X= 0
diff --git a/Inputs/EventGenerator/132Sndp.reaction b/Inputs/EventGenerator/132Sndp.reaction
index 0ff5e1462..0618fe3f6 100644
--- a/Inputs/EventGenerator/132Sndp.reaction
+++ b/Inputs/EventGenerator/132Sndp.reaction
@@ -12,8 +12,8 @@ Transfert
BeamEnergySpread= 0
SigmaX= 0.851
SigmaY= 0.851
- SigmaThetaX= 0
- SigmaPhiY= 0
+ SigmaThetaX= 1
+ SigmaPhiY= 1
CrossSectionPath= sn132dp_gs_10AMeV.txt
ShootLight= 1
ShootHeavy= 0
diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index 52b4d1497..136487077 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -28,12 +28,13 @@
// C++ header
#include <string>
-// G4 header defining G$ types
+// G4 header defining G4 types
#include "globals.hh"
// G4 headers
#include "G4ParticleGun.hh"
#include "G4Event.hh"
+#include "G4ThreeVector.hh"
// NPTool headers
#include "Target.hh"
@@ -56,14 +57,17 @@ public:
// Used in some case to generate event inside the target
virtual void SetTarget(Target*) {};
- virtual void SetTargetThickness(G4double) {};
- virtual void SetTargetAngle(G4double) {};
- virtual void SetTargetRadius(G4double) {};
- virtual void SetTargetCoordinate(G4double, G4double, G4double) {};
-
- // Used to simulate beam emmitance effect
- void RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y, double NumberOfSigma=10000);
+ // Used to calculate the incident beam direction (taking into account
+ // the emittance) and the vertex of interaction in target
+ void CalculateBeamInteraction(double MeanPosX, double SigmaPosX, double MeanPosTheta, double SigmaPosTheta,
+ double MeanPosY, double SigmaPosY, double MeanPosPhi, double SigmaPosPhi,
+ Target* target,
+ G4ThreeVector &InterCoord, double &AngleEmittanceTheta, double &AngleEmittancePhi,
+ double &AngleIncidentTheta, double &AngleIncidentPhi);
+
+ // Used to simulate beam emmitance effect
+ void RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y, double NumberOfSigma = 10000);
};
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index 1b6410dc6..b3071b08a 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -173,58 +173,36 @@ void EventGeneratorBeam::ReadConfiguration(string Path)
void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun)
{
m_InitConditions->Clear();
- // Vertex position and beam angle
- G4double x0 = 1000 * cm ;
- G4double y0 = 1000 * cm ;
- G4double z0 = 0;
- G4double Beam_thetaX = 0 ;
- G4double Beam_phiY = 0 ;
- G4double TargetThick = 0;
- //shoot inside the target with correlated angle
- if (m_Target->GetTargetRadius() != 0)
- {
- while (sqrt(x0*x0 + y0*y0) > m_Target->GetTargetRadius())
- {
- RandomGaussian2D(0 , 0 , m_SigmaX , m_SigmaThetaX , x0 , Beam_thetaX );
- RandomGaussian2D(0 , 0 , m_SigmaY , m_SigmaPhiY , y0 , Beam_phiY );
- G4double dz = x0 * tan(m_Target->GetTargetAngle());
- TargetThick = m_Target->GetTargetThickness() / cos(m_Target->GetTargetAngle());
- z0 = dz + (-TargetThick / 2 + RandFlat::shoot() * TargetThick);
- }
- }
-
- else
- {
- RandomGaussian2D( 0 , 0 , 0 , m_SigmaThetaX , x0 , Beam_thetaX );
- RandomGaussian2D( 0 , 0 , 0 , m_SigmaPhiY , y0 , Beam_phiY );
- TargetThick = m_Target->GetTargetThickness() / cos(m_Target->GetTargetAngle());
- z0 = (-TargetThick / 2 + RandFlat::shoot() * TargetThick);
- }
-
- m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
- m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
-
-
- // Calculate Angle in spherical coordinate, passing by the direction vector dir
- G4double Xdir = cos( pi/2. - Beam_thetaX ) ;
- G4double Ydir = cos( pi/2. - Beam_phiY ) ;
- G4double Zdir = sin( pi/2. - Beam_thetaX ) + sin( pi/2. - Beam_phiY) ;
-
- G4double Beam_theta = acos ( Zdir / sqrt( Xdir*Xdir + Ydir*Ydir + Zdir*Zdir ) );
- G4double Beam_phi = atan2( Ydir , Xdir );
+ ///////////////////////////////////////////////////////////////////////
+ ///// Calculate the incident beam direction as well as the vertex /////
+ ///// of interaction in target /////
+ ///////////////////////////////////////////////////////////////////////
+ G4ThreeVector InterCoord;
+ G4double Beam_thetaX = 0, Beam_phiY = 0;
+ G4double Beam_theta = 0, Beam_phi = 0;
+ CalculateBeamInteraction(0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ m_Target,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi);
+
+ // write vertex position to ROOT file
+ G4double x0 = InterCoord.x();
+ G4double y0 = InterCoord.y();
+ G4double z0 = InterCoord.z();
+ m_InitConditions->SetICPositionX(x0);
+ m_InitConditions->SetICPositionY(y0);
+ m_InitConditions->SetICPositionZ(z0);
+
+ // write emittance angles to ROOT file
+ m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
+ m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
- // Move to the target
- x0 += m_Target->GetTargetX() ;
- y0 += m_Target->GetTargetY() ;
- z0 += m_Target->GetTargetZ() ;
-
// Store initial value
m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
- m_InitConditions->SetICPositionX(x0);
- m_InitConditions->SetICPositionY(y0);
- m_InitConditions->SetICPositionZ(z0);
+
//////////////////////////////////////////////////
/////Now define everything for light particle/////
//////////////////////////////////////////////////
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index a38a7129b..87bba8e11 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -68,7 +68,6 @@ EventGeneratorTransfert::~EventGeneratorTransfert()
//------------- Default Destructor ------------
delete m_InitConditions;
delete m_Reaction;
- delete m_Target;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -333,37 +332,23 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
G4ParticleDefinition* HeavyName
= G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation()*MeV);
- // Vertex position and beam angle inte world frame
- G4double x0 = 1000 * cm ;
- G4double y0 = 1000 * cm ;
- G4double z0 = 0;
- G4double Beam_thetaX = 0 ;
- G4double Beam_phiY = 0 ;
- G4double TargetThick = 0;
-
- // Shoot inside the target with correlated angle
- if (m_Target->GetTargetRadius() != 0) {
- while (sqrt(x0*x0 + y0*y0) > m_Target->GetTargetRadius()) {
- RandomGaussian2D(0, 0, m_SigmaX, m_SigmaThetaX, x0, Beam_thetaX);
- RandomGaussian2D(0, 0, m_SigmaY, m_SigmaPhiY , y0, Beam_phiY );
- }
- G4double dz = x0 * tan(m_Target->GetTargetAngle());
- TargetThick = m_Target->GetTargetThickness() / cos(m_Target->GetTargetAngle());
- z0 = dz + (-TargetThick / 2 + RandFlat::shoot() * TargetThick);
- }
- else {
- RandomGaussian2D(0,0,0,m_SigmaThetaX,x0,Beam_thetaX);
- RandomGaussian2D(0,0,0,m_SigmaPhiY ,y0,Beam_phiY );
- TargetThick = m_Target->GetTargetThickness() / cos(m_Target->GetTargetAngle());
- z0 = (-TargetThick / 2 + RandFlat::shoot() * TargetThick);
- }
-
- // Move to the target
- x0 += m_Target->GetTargetX() ;
- y0 += m_Target->GetTargetY() ;
- z0 += m_Target->GetTargetZ() ;
+ ///////////////////////////////////////////////////////////////////////
+ ///// Calculate the incident beam direction as well as the vertex /////
+ ///// of interaction in target /////
+ ///////////////////////////////////////////////////////////////////////
+ G4ThreeVector InterCoord;
+ G4double Beam_thetaX = 0, Beam_phiY = 0;
+ G4double Beam_theta = 0, Beam_phi = 0;
+ CalculateBeamInteraction(0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ m_Target,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi);
// write vertex position to ROOT file
+ G4double x0 = InterCoord.x();
+ G4double y0 = InterCoord.y();
+ G4double z0 = InterCoord.z();
m_InitConditions->SetICPositionX(x0);
m_InitConditions->SetICPositionY(y0);
m_InitConditions->SetICPositionZ(z0);
@@ -372,16 +357,6 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
- // Calculate Angle in spherical coordinate, passing by the direction vector dir
- G4double Xdir = cos(pi/2. - Beam_thetaX);
- G4double Ydir = cos(pi/2. - Beam_phiY );
- G4double Zdir = sin(pi/2. - Beam_thetaX) + sin(pi/2. - Beam_phiY);
-
- G4double Beam_theta = acos(Zdir / sqrt(Xdir*Xdir + Ydir*Ydir + Zdir*Zdir)) * rad;
- G4double Beam_phi = atan2(Ydir, Xdir) * rad;
- if (Beam_phi < 0) Beam_phi += 2*pi;
- if (Beam_theta < 1e-6) Beam_phi = 0;
-
// write angles to ROOT file
m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index 46261eb7b..d3836acf4 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -391,44 +391,31 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
m_Reaction->SetExcitation( EXX );
- // Vertex position and beam angle inte world frame
- G4double x0 = 1000 * cm ;
- G4double y0 = 1000 * cm ;
- G4double z0 =0;
- G4double Beam_thetaX = 0 ;
- G4double Beam_phiY = 0 ;
- G4double TargetThick = 0;
-
- //shoot inside the target with correlated angle
- if (m_Target->GetTargetRadius() != 0) {
- while (sqrt(x0*x0 + y0*y0) > m_Target->GetTargetRadius()) {
- RandomGaussian2D(0, 0, m_SigmaX, m_SigmaThetaX, x0, Beam_thetaX);
- RandomGaussian2D(0, 0, m_SigmaY, m_SigmaPhiY , y0, Beam_phiY );
- }
- G4double dz = x0 * tan(m_Target->GetTargetAngle());
- TargetThick = m_Target->GetTargetThickness() / cos(m_Target->GetTargetAngle());
- z0 = dz + (-TargetThick / 2 + RandFlat::shoot() * TargetThick);
- }
- else {
- RandomGaussian2D(0,0,0,m_SigmaThetaX,x0,Beam_thetaX);
- RandomGaussian2D(0,0,0,m_SigmaPhiY ,y0,Beam_phiY );
- TargetThick = m_Target->GetTargetThickness() / cos(m_Target->GetTargetAngle());
- z0 = (-TargetThick / 2 + RandFlat::shoot() * TargetThick);
- }
+ ///////////////////////////////////////////////////////////////////////
+ ///// Calculate the incident beam direction as well as the vertex /////
+ ///// of interaction in target /////
+ ///////////////////////////////////////////////////////////////////////
+ G4ThreeVector InterCoord;
+ G4double Beam_thetaX = 0, Beam_phiY = 0;
+ G4double Beam_theta = 0, Beam_phi = 0;
+ CalculateBeamInteraction(0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ m_Target,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi);
+
+ // write vertex position to ROOT file
+ G4double x0 = InterCoord.x();
+ G4double y0 = InterCoord.y();
+ G4double z0 = InterCoord.z();
+ m_InitConditions->SetICPositionX(x0);
+ m_InitConditions->SetICPositionY(y0);
+ m_InitConditions->SetICPositionZ(z0);
// write emittance angles to ROOT file
m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
- // Calculate Angle in spherical coordinate, passing by the direction vector dir
- G4double Xdir = cos(pi/2. - Beam_thetaX);
- G4double Ydir = cos(pi/2. - Beam_phiY );
- G4double Zdir = sin(pi/2. - Beam_thetaX) + sin(pi/2. - Beam_phiY);
-
- G4double Beam_theta = acos(Zdir / sqrt(Xdir*Xdir + Ydir*Ydir + Zdir*Zdir)) * rad;
- G4double Beam_phi = atan2(Ydir, Xdir) * rad;
- if (Beam_phi < 0) Beam_phi += 2*pi;
-
// write angles to ROOT file
m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index 4c91b8696..6ab94ebb3 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -25,41 +25,99 @@
*****************************************************************************/
#include "VEventGenerator.hh"
+// C++ headers
+#include "cmath"
+
// ROOT headers
-#include"TMath.h"
+//#include"TMath.h"
// CLHEP headers
#include "Randomize.hh"
+
VEventGenerator::VEventGenerator()
- {
- }
+{
+}
+
+
VEventGenerator::~VEventGenerator()
- {
- }
-
-void VEventGenerator::RandomGaussian2D(double MeanX,double MeanY,double SigmaX,double SigmaY,double &X,double &Y,double NumberOfSigma)
- {
-
- if(SigmaX!=0)
- {
- X=2* NumberOfSigma *SigmaX;
- while(X>NumberOfSigma*SigmaX)
- X = RandGauss::shoot( MeanX , SigmaX);
- double a = NumberOfSigma * SigmaX/2 ;
- double b = NumberOfSigma * SigmaY/2 ;
-
- double SigmaYPrim = b * sqrt( 1 - X*X / (a*a) ) ;
- SigmaYPrim = 2*SigmaYPrim / NumberOfSigma ;
-
- Y = RandGauss::shoot( MeanY , SigmaYPrim) ;
- }
-
- else
- {
- X = MeanX ;
- Y = RandGauss::shoot( MeanY , SigmaY) ;
- }
-
- }
+{
+}
+
+
+
+void VEventGenerator::CalculateBeamInteraction(double MeanPosX, double SigmaPosX, double MeanPosTheta, double SigmaPosTheta,
+ double MeanPosY, double SigmaPosY, double MeanPosPhi, double SigmaPosPhi,
+ Target* target,
+ G4ThreeVector &InterCoord, double &AngleEmittanceTheta, double &AngleEmittancePhi,
+ double &AngleIncidentTheta, double &AngleIncidentPhi)
+{
+ // target parameters
+ double TargetThickness = target->GetTargetThickness();
+ double TargetRadius = target->GetTargetRadius();
+ double TargetAngle = target->GetTargetAngle();
+
+ // beam interaction parameters
+ double x0 = 1000 * cm;
+ double y0 = 1000 * cm;
+ double z0 = 0 * cm;
+ double dz = 0 * cm;
+
+ // calculate emittance parameters (x,theta) and (y,phi)
+ if (TargetRadius != 0) { // case of finite target dimensions
+ while (sqrt(x0*x0 + y0*y0) > TargetRadius) {
+ RandomGaussian2D(MeanPosX, MeanPosTheta, SigmaPosX, SigmaPosTheta, x0, AngleEmittanceTheta);
+ RandomGaussian2D(MeanPosY, MeanPosPhi, SigmaPosY, SigmaPosPhi, y0, AngleEmittancePhi);
+ }
+ // in case target is tilted, correct the z-position of interaction
+ dz = x0 * tan(TargetAngle);
+ }
+ else { // if no target radius is given consider a point-like target
+ RandomGaussian2D(0, 0, 0, SigmaPosTheta, x0, AngleEmittanceTheta);
+ RandomGaussian2D(0, 0, 0, SigmaPosPhi, y0, AngleEmittancePhi);
+ }
+
+ // correct for the target angle wrt the beam axis
+ // this simple correction is only valid if the beam is parallel to the beam axis
+ // should be improved in a next version
+ TargetThickness /= cos(TargetAngle);
+ z0 = dz + (-TargetThickness / 2 + RandFlat::shoot() * TargetThickness);
+
+ // Move to the target position
+ x0 += target->GetTargetX();
+ y0 += target->GetTargetY();
+ z0 += target->GetTargetZ();
+ InterCoord = G4ThreeVector(x0, y0, z0);
+
+ // Calculate incident angle in spherical coordinate, passing by the direction vector dir
+ double Xdir = sin(AngleEmittanceTheta);
+ double Ydir = sin(AngleEmittancePhi);
+ double Zdir = cos(AngleEmittanceTheta) + cos(AngleEmittancePhi);
+
+ AngleIncidentTheta = acos(Zdir / sqrt(Xdir*Xdir + Ydir*Ydir + Zdir*Zdir)) * rad;
+ AngleIncidentPhi = atan2(Ydir, Xdir) * rad;
+ if (AngleIncidentPhi < 0) AngleIncidentPhi += 2*pi;
+ if (AngleIncidentTheta < 1e-6) AngleIncidentPhi = 0;
+}
+
+
+
+void VEventGenerator::RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y, double NumberOfSigma)
+{
+ if (SigmaX != 0) {
+ X = 2 * NumberOfSigma*SigmaX;
+ while (X > NumberOfSigma*SigmaX) X = RandGauss::shoot(MeanX, SigmaX);
+
+ double a = NumberOfSigma * SigmaX/2;
+ double b = NumberOfSigma * SigmaY/2;
+ double SigmaYPrim = b * sqrt(1 - X*X / (a*a));
+
+ SigmaYPrim = 2*SigmaYPrim / NumberOfSigma;
+ Y = RandGauss::shoot(MeanY, SigmaYPrim);
+ }
+ else {
+ X = MeanX;
+ Y = RandGauss::shoot(MeanY, SigmaY);
+ }
+}
--
GitLab