From 3337cfe73ec814fc03cb0f18951ce26910a4f221 Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Fri, 13 Nov 2009 00:42:11 +0000
Subject: [PATCH] * Now interaction coordinate and Energy loss of the beam are
compute using a target method, called in the different event generator
- Some test have been done, results are OK - Simulation is slightly slow
done by this process of slowing done the beam, mainly because an Energy loss
table is called every events...
---
.../DetectorConfiguration/Riken_65mm.detector | 7 +-
Inputs/DetectorConfiguration/e530.detector | 1 +
NPSimulation/include/Target.hh | 58 +++++---
NPSimulation/include/VEventGenerator.hh | 11 +-
NPSimulation/src/EventGeneratorBeam.cc | 23 +--
NPSimulation/src/EventGeneratorTransfert.cc | 50 +++----
.../src/EventGeneratorTransfertToResonance.cc | 39 ++++--
NPSimulation/src/Target.cc | 132 +++++++++++++++---
NPSimulation/src/VEventGenerator.cc | 99 -------------
9 files changed, 212 insertions(+), 208 deletions(-)
diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector
index 31d9cbda5..a0e5ee510 100644
--- a/Inputs/DetectorConfiguration/Riken_65mm.detector
+++ b/Inputs/DetectorConfiguration/Riken_65mm.detector
@@ -7,15 +7,16 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GeneralTarget
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%Target
- THICKNESS= 1
+Target
+ THICKNESS= 10000
RADIUS= 45
MATERIAL= CD2
+ ANGLE= 45
X= 0
Y= 0
Z= 0
-CryoTarget
+%CryoTarget
THICKNESS= 3000
RADIUS= 45
TEMPERATURE= 26
diff --git a/Inputs/DetectorConfiguration/e530.detector b/Inputs/DetectorConfiguration/e530.detector
index 444075b0c..fbe124019 100644
--- a/Inputs/DetectorConfiguration/e530.detector
+++ b/Inputs/DetectorConfiguration/e530.detector
@@ -22,6 +22,7 @@ Target
THICKNESS= 309.278350515
RADIUS= 7.5
MATERIAL= CD2
+ ANGLE= 45
X= 0
Y= 0
Z= 40
diff --git a/NPSimulation/include/Target.hh b/NPSimulation/include/Target.hh
index f34c7677e..3b24daf07 100644
--- a/NPSimulation/include/Target.hh
+++ b/NPSimulation/include/Target.hh
@@ -52,29 +52,41 @@ public:
public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(string Path);
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(string Path);
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world);
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world);
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput();
-
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event);
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput();
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event);
public:
- // method for debug purpose (still to be implemented)
- // This method should check if the results of the beam interaction within the target
- // (interaction coordinates) are well located inside the target volume
- bool IsInsideTarget() {return false;};
-
+ // method for debug purpose (still to be implemented)
+ // This method should check if the results of the beam interaction within the target
+ // (interaction coordinates) are well located inside the target volume
+ bool IsInsideTarget() {return false;};
+
+ // Used to calculate the incident beam direction (taking into account
+ // the emittance) and the vertex of interaction in target
+ // Also compute the energy lost by the beam in the target before interacting
+ void CalculateBeamInteraction(double MeanPosX, double SigmaPosX, double MeanPosTheta, double SigmaPosTheta,
+ double MeanPosY, double SigmaPosY, double MeanPosPhi, double SigmaPosPhi,
+ double IncidentBeamEnergy,
+ G4ParticleDefinition* BeamName,
+ G4ThreeVector &InterCoord, double &AngleEmittanceTheta, double &AngleEmittancePhi,
+ double &AngleIncidentTheta, double &AngleIncidentPhi,
+ double &FinalBeamEnergy);
+
+ // Used to simulate beam emmitance effect
+ void RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y, double NumberOfSigma = 10000);
public:
G4Material* GetMaterialFromLibrary(G4String MaterialName, G4double Temperature = 0, G4double Pressure = 0);
@@ -83,11 +95,11 @@ public:
public:
G4double GetTargetThickness() {return m_TargetThickness;}
G4Material* GetTargetMaterial() {return m_TargetMaterial;}
- G4double GetTargetRadius() {return m_TargetRadius;}
- G4double GetTargetAngle() {return m_TargetAngle;}
- G4double GetTargetX() {return m_TargetX;}
- G4double GetTargetY() {return m_TargetY;}
- G4double GetTargetZ() {return m_TargetZ;}
+ G4double GetTargetRadius() {return m_TargetRadius;}
+ G4double GetTargetAngle() {return m_TargetAngle;}
+ G4double GetTargetX() {return m_TargetX;}
+ G4double GetTargetY() {return m_TargetY;}
+ G4double GetTargetZ() {return m_TargetZ;}
G4int GetTargetNbLayers() {return m_TargetNbLayers;}
diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index f5d040a10..59bbd3bca 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -59,16 +59,7 @@ public:
// Used in some case to generate event inside the target
virtual void SetTarget(Target*) {};
-
- // 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,
- double &EffectiveTargetThicknessBeforeInteraction);
-
+
// 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 7388505fa..604ce59f8 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -179,19 +179,24 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
///////////////////////////////////////////////////////////////////////
///// Calculate the incident beam direction as well as the vertex /////
- ///// of interaction in target /////
+ ///// of interaction in target and Energy Loss of the beam within /////
+ ///// the target. /////
///////////////////////////////////////////////////////////////////////
G4ThreeVector InterCoord;
+
G4double Beam_thetaX = 0, Beam_phiY = 0;
G4double Beam_theta = 0, Beam_phi = 0;
- G4double EffectiveThickness = 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,
- EffectiveThickness);
-
+ G4double FinalBeamEnergy = 0 ;
+ G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
+
+ m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ InitialBeamEnergy,
+ m_particle,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi,
+ FinalBeamEnergy);
+
// write vertex position to ROOT file
G4double x0 = InterCoord.x();
G4double y0 = InterCoord.y();
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 0f7e15c94..af94442fb 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -31,7 +31,6 @@
// G4 headers
#include "G4ParticleTable.hh"
-#include "G4EmCalculator.hh"
#include "G4ParticleGun.hh"
#include "G4RotationMatrix.hh"
@@ -343,19 +342,27 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
///////////////////////////////////////////////////////////////////////
///// Calculate the incident beam direction as well as the vertex /////
- ///// of interaction in target /////
+ ///// of interaction in target and Energy Loss of the beam within /////
+ ///// the target. /////
///////////////////////////////////////////////////////////////////////
G4ThreeVector InterCoord;
+
G4double Beam_thetaX = 0, Beam_phiY = 0;
G4double Beam_theta = 0, Beam_phi = 0;
- G4double EffectiveThickness = 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,
- EffectiveThickness);
-
+ G4double FinalBeamEnergy = 0 ;
+ G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
+
+ m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ InitialBeamEnergy,
+ BeamName,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi,
+ FinalBeamEnergy);
+
+ m_Reaction->SetBeamEnergy(FinalBeamEnergy);
+ m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
+
// write vertex position to ROOT file
G4double x0 = InterCoord.x();
G4double y0 = InterCoord.y();
@@ -391,21 +398,7 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
///// Angles for emitted particles following Cross Section //////
///// Angles are in the beam frame //////
/////////////////////////////////////////////////////////////////
- // Beam nominal incident energy before interaction with the target
- G4double NominalBeamEnergy = m_BeamEnergy;
- G4double IncidentBeamEnergy = RandGauss::shoot(NominalBeamEnergy, m_BeamEnergySpread / 2.35);
- // Slowing down the beam to the interaction layer in the target
- // Number of Layers
- G4int NbLayers = m_Target->GetTargetNbLayers();
- G4EmCalculator emCalculator;
- for (G4int i = 0; i < NbLayers; i++) {
-// G4double dedx = emCalculator.GetDEDX(IncidentBeamEnergy, BeamName, m_Target->GetTargetMaterial());
- G4double dedx = emCalculator.ComputeTotalDEDX(IncidentBeamEnergy, BeamName, m_Target->GetTargetMaterial());
- G4double de = dedx * EffectiveThickness / NbLayers;
- IncidentBeamEnergy -= de;
- }
- m_Reaction->SetBeamEnergy(IncidentBeamEnergy);
- m_InitConditions->SetICIncidentEnergy(IncidentBeamEnergy / MeV);
+
// Angles
RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
G4double ThetaCM = CrossSectionShoot.shoot() * (180*deg);
@@ -458,16 +451,13 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
// write angles in ROOT file
m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
- // tests
-// G4cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXX" << G4endl;
-// G4cout << "cinematique dans ref world : " << G4endl;
-// G4cout << "\t" << momentum_kine_world << G4endl;
+
//Set the gun to shoot
particleGun->SetParticleMomentumDirection(momentum_kine_world);
//Shoot the light particle
particleGun->GeneratePrimaryVertex(anEvent);
}
- if (m_ShootHeavy) { // Case of recoil particle
+ if (m_ShootHeavy) { // Case of heavy particle
// Particle type
particleGun->SetParticleDefinition(HeavyName);
// Particle energy
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index f513b2cb0..c12b8b13f 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -371,18 +371,25 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
//////////////////////////////////////////////////
//////Define the kind of particle to shoot////////
//////////////////////////////////////////////////
+ // Light
G4int LightZ = m_Reaction->GetNucleus3()->GetZ() ;
G4int LightA = m_Reaction->GetNucleus3()->GetA() ;
G4ParticleDefinition* LightName
= G4ParticleTable::GetParticleTable()->GetIon(LightZ, LightA, 0.);
+ // Recoil
G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
G4ParticleDefinition* HeavyName
= G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation()*MeV);
+ // Beam
+ G4int BeamZ = m_Reaction->GetNucleus1()->GetZ();
+ G4int BeamA = m_Reaction->GetNucleus1()->GetA();
+ G4ParticleDefinition* BeamName = G4ParticleTable::GetParticleTable()->GetIon(BeamZ, BeamA, 0);
+
// Shoot the Resonance energy following the mean and width value
double EXX = -10 ;
@@ -391,27 +398,29 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
m_Reaction->SetExcitation( EXX );
- ///////////////////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////
///// Calculate the incident beam direction as well as the vertex /////
- ///// of interaction in target /////
+ ///// of interaction in target and Energy Loss of the beam within /////
+ ///// the target. /////
///////////////////////////////////////////////////////////////////////
G4ThreeVector InterCoord;
+
G4double Beam_thetaX = 0, Beam_phiY = 0;
G4double Beam_theta = 0, Beam_phi = 0;
- G4double EffectiveThickness = 0;
-
- if(m_SigmaX==0 && m_SigmaY==0 && m_SigmaThetaX == 0 && m_SigmaPhiY==0) InterCoord = G4ThreeVector(0,0,0);
+ G4double FinalBeamEnergy = 0 ;
+ G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
- else if( m_Target !=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,
- EffectiveThickness);
-
- else
- InterCoord = G4ThreeVector(0,0,0);
+ m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
+ 0, m_SigmaY, 0, m_SigmaPhiY,
+ InitialBeamEnergy,
+ BeamName,
+ InterCoord, Beam_thetaX, Beam_phiY,
+ Beam_theta, Beam_phi,
+ FinalBeamEnergy);
+
+ m_Reaction->SetBeamEnergy(FinalBeamEnergy);
+ m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
+
// write vertex position to ROOT file
G4double x0 = InterCoord.x();
diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc
index 798c028ee..4dc4efa24 100644
--- a/NPSimulation/src/Target.cc
+++ b/NPSimulation/src/Target.cc
@@ -21,7 +21,7 @@
* + 16/09/2009: Add support for positioning the target with an angle with *
* respect to the beam (N. de Sereville) *
* + 16/09/2009: Add CH2 material for targets (N. de Sereville) *
- * + 06/11/2009: Add new Token NBLAYERS defining the number of steps used *
+ * + 06/11/2009: Add new Token m_TargetNbLayers defining the number of steps used *
* to slow down the beam in the target (N. de Sereville) *
* *
*****************************************************************************/
@@ -41,7 +41,9 @@
#include "G4VPhysicalVolume.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
-
+#include "G4EmCalculator.hh"
+#include "Randomize.hh"
+using namespace CLHEP ;
// NPTool header
#include"Target.hh"
@@ -212,7 +214,7 @@ void Target::ReadConfiguration(string Path)
bool check_X = false ;
bool check_Y = false ;
bool check_Z = false ;
- bool check_NbLayers = false;
+ bool check_m_TargetNbLayers = false;
bool check_Temperature = false ;
bool check_Pressure = false ;
@@ -287,8 +289,8 @@ void Target::ReadConfiguration(string Path)
cout << m_TargetZ / mm << " )" << endl ;
}
- else if (DataBuffer.compare(0, 9, "NBLAYERS=") == 0) {
- check_NbLayers = true ;
+ else if (DataBuffer.compare(0, 9, "NbLayers=") == 0) {
+ check_m_TargetNbLayers = true ;
ConfigFile >> DataBuffer;
m_TargetNbLayers = atoi(DataBuffer.c_str());
cout << "Number of steps for slowing down the beam in target: " << m_TargetNbLayers << endl;
@@ -381,8 +383,8 @@ void Target::ReadConfiguration(string Path)
cout << m_TargetZ / mm << " )" << endl ;
}
- else if (DataBuffer.compare(0, 9, "NBLAYERS=") == 0) {
- check_NbLayers = true ;
+ else if (DataBuffer.compare(0, 9, "m_TargetNbLayers=") == 0) {
+ check_m_TargetNbLayers = true ;
ConfigFile >> DataBuffer;
m_TargetNbLayers = atoi(DataBuffer.c_str());
cout << "Number of steps for slowing down the beam in target: " << m_TargetNbLayers << endl;
@@ -457,20 +459,20 @@ void Target::ConstructDetector(G4LogicalVolume* world)
G4LogicalVolume* logicWindowsB = new G4LogicalVolume(solidWindowsB, m_WindowsMaterial, "logicTargetWindowsB");
PVPBuffer =
- new G4PVPlacement(0 ,
- TargetPos + G4ThreeVector(0., 0., 0.5*(m_TargetThickness + m_WindowsThickness)) ,
- logicWindowsF ,
- "Target Window Front" ,
- world ,
- false, 0);
+ new G4PVPlacement( 0 ,
+ TargetPos + G4ThreeVector(0., 0., 0.5*(m_TargetThickness + m_WindowsThickness)) ,
+ logicWindowsF ,
+ "Target Window Front" ,
+ world ,
+ false, 0 );
PVPBuffer =
- new G4PVPlacement(0 ,
- TargetPos + G4ThreeVector(0., 0., -0.5*(m_TargetThickness + m_WindowsThickness)) ,
- logicWindowsB ,
- "Target Window Back" ,
- world ,
- false, 0);
+ new G4PVPlacement( 0 ,
+ TargetPos + G4ThreeVector(0., 0., -0.5*(m_TargetThickness + m_WindowsThickness)) ,
+ logicWindowsB ,
+ "Target Window Back" ,
+ world ,
+ false, 0 );
G4VisAttributes* WindowsVisAtt = new G4VisAttributes(G4Colour(0.5, 1., 0.5));
logicWindowsF->SetVisAttributes(WindowsVisAtt);
@@ -489,4 +491,96 @@ void Target::InitializeRootOutput()
// Called at in the EventAction::EndOfEventAvtion
void Target::ReadSensitive(const G4Event*)
{}
+
+void Target::CalculateBeamInteraction( double MeanPosX, double SigmaPosX, double MeanPosTheta, double SigmaPosTheta,
+ double MeanPosY, double SigmaPosY, double MeanPosPhi, double SigmaPosPhi,
+ double IncidentBeamEnergy,
+ G4ParticleDefinition* BeamName,
+ G4ThreeVector &InterCoord, double &AngleEmittanceTheta, double &AngleEmittancePhi,
+ double &AngleIncidentTheta, double &AngleIncidentPhi,
+ double &FinalBeamEnergy)
+{
+
+ // target parameters
+ G4ThreeVector TargetNormal = G4ThreeVector( sin(m_TargetAngle) ,
+ 0 ,
+ cos(m_TargetAngle) );
+
+ // 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 (m_TargetRadius != 0) { // case of finite target dimensions
+ while (sqrt(x0*x0 + y0*y0) > m_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
+ // x is the vertical axis
+ dz = x0 * tan(m_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);
+ }
+
+ // 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);
+ G4ThreeVector BeamDir = G4ThreeVector(Xdir,Ydir,Zdir) ;
+
+ AngleIncidentTheta = BeamDir.theta() ;
+ AngleIncidentPhi = BeamDir.phi() ;
+ if (AngleIncidentPhi < 0) AngleIncidentPhi += 2*pi ;
+ if (AngleIncidentTheta < 1e-6) AngleIncidentPhi = 0 ;
+
+ // Calculation of effective target thickness and z-position of interaction
+ // when the target is tilted wrt the beam axis
+ double EffectiveThickness = m_TargetThickness / (BeamDir.unit()).dot(TargetNormal.unit());
+ double uniform = RandFlat::shoot();
+ z0 = dz + (-m_TargetThickness / 2 + uniform * m_TargetThickness);
+
+ // Calculate the effective thickness before interaction in target
+ // This is useful to slow down the beam
+ double EffectiveTargetThicknessBeforeInteraction = m_TargetThickness * uniform / cos(AngleIncidentTheta);
+
+ // Move to the target position
+ x0 += m_TargetX;
+ y0 += m_TargetY;
+ z0 += m_TargetZ;
+ InterCoord = G4ThreeVector(x0, y0, z0);
+
+ G4EmCalculator emCalculator;
+ for (G4int i = 0; i < m_TargetNbLayers; i++)
+ {
+ G4double dedx = emCalculator.ComputeTotalDEDX(IncidentBeamEnergy, BeamName, m_TargetMaterial);
+ G4double de = dedx * EffectiveTargetThicknessBeforeInteraction / m_TargetNbLayers;
+ IncidentBeamEnergy -= de;
+ }
+FinalBeamEnergy=IncidentBeamEnergy;
+}
+
+
+void Target::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);
+ }
+}
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index a1f0acb26..bf2b8e82b 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -47,102 +47,3 @@ VEventGenerator::~VEventGenerator()
{
}
-
-
-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,
- double &EffectiveTargetThicknessBeforeInteraction)
-{
- if (target != 0) {
- // 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
- // x is the vertical axis
- 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);
- }
-
- // 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;
-
- // Calculation of effective target thickness and z-position of interaction
- // when the target is tilted wrt the beam axis
- // * exact if target is perpendicular to the beam axis (target not tilted)
- // in any case of beam emittance
- // * exact if target is tilted wrt the beam axis and the beam is parallel
- // (no beam divergence)
- // * wrong if target is tilted and for general case of beam emittance
- // (should be fixed in next release). For small beam divergence this
- // should not be such a big problem
- TargetThickness /= cos(TargetAngle);
- double uniform = RandFlat::shoot();
- z0 = dz + (-TargetThickness / 2 + uniform * TargetThickness);
-
- // Calculate the effective thickness before interaction in target
- // This is useful to slow down the beam
- EffectiveTargetThicknessBeforeInteraction = TargetThickness * uniform / cos(AngleIncidentTheta);
-
- // Move to the target position
- x0 += target->GetTargetX();
- y0 += target->GetTargetY();
- z0 += target->GetTargetZ();
- InterCoord = G4ThreeVector(x0, y0, z0);
-
- } // end if target
- else {
- InterCoord = G4ThreeVector(0, 0, 0);
- AngleEmittanceTheta = 0;
- AngleEmittancePhi = 0;
- AngleIncidentTheta = 0;
- 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