From 57a57f343612fce664367cc105ef368d24f0084b Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Tue, 3 Nov 2009 12:56:50 +0000
Subject: [PATCH] * Slow down the beam to the interaction point in the target
+ Target.{hh,cc} class now has a method returning the target material +
VEventGenerator.{hh,cc} class now returns the effective thickness of
target seen by the beam before interaction (method CalculateBeamInteraction)
+ Propagate these changes in the different EvenGenerators
* Update Gaspard documentation
---
INSTALL | 2 +-
NPAnalysis/Gaspard/RunToTreat.txt | 4 +-
NPDocumentation/Gaspard.tex | 4 +-
NPSimulation/include/Target.hh | 13 ++---
NPSimulation/include/VEventGenerator.hh | 7 ++-
NPSimulation/src/EventGeneratorBeam.cc | 4 +-
NPSimulation/src/EventGeneratorTransfert.cc | 24 +++++++++-
.../src/EventGeneratorTransfertToResonance.cc | 4 +-
NPSimulation/src/VEventGenerator.cc | 47 ++++++++++++-------
TODO | 7 ++-
10 files changed, 79 insertions(+), 37 deletions(-)
diff --git a/INSTALL b/INSTALL
index 35f9cbaf7..2235f0b7d 100644
--- a/INSTALL
+++ b/INSTALL
@@ -29,7 +29,7 @@ I) REQUIREMENTS
II) WORKING CONFIGURATIONS
The NPTool package has been mainly tested (and developped) with
the two following configurations:
- 1) Mac OS X (10.5.7) + G4 (4.9.2) + CLHEP (2.0.4.2) + ROOT (5.22/00 and 5.24/00) + GSL ()
+ 1) Mac OS X (10.5.7) + G4 (4.9.2 and 4.9.3.b1) + CLHEP (2.0.4.2) + ROOT (5.22/00 and 5.24/00) + GSL ()
2) Linux Fedora (kernel 2.6.29) + G4 (4.9.1p3) + CLHEP (2.0.4.2) + ROOT (5.22/00) + GSL (1.12)
Please, report any working or non-working configuration.
diff --git a/NPAnalysis/Gaspard/RunToTreat.txt b/NPAnalysis/Gaspard/RunToTreat.txt
index e1e642cec..5e6a33893 100644
--- a/NPAnalysis/Gaspard/RunToTreat.txt
+++ b/NPAnalysis/Gaspard/RunToTreat.txt
@@ -1,5 +1,5 @@
TTreeName
SimulatedTree
RootFileName
-% ../../Outputs/Simulation/mySimul.root
- ../../Outputs/Simulation/sn132dp_10MeVA_T1_B0_E0.root
+ ../../Outputs/Simulation/mySimul.root
+% ../../Outputs/Simulation/sn132dp_10MeVA_T1_B0_E0.root
diff --git a/NPDocumentation/Gaspard.tex b/NPDocumentation/Gaspard.tex
index 5f0ccce02..e486f4bc7 100755
--- a/NPDocumentation/Gaspard.tex
+++ b/NPDocumentation/Gaspard.tex
@@ -169,10 +169,10 @@ contains three classes:
\subsection{Adding a new detector shape to Gaspard}
A special class (GaspardTrackerDummyShape) has been created to show how
to define a new module for the Gaspard Tracker. This class describes a
-simple 5x5 cm2 square telescope made of three layers of silicon which
+simple 5x5 cm$^2$ square telescope made of three layers of silicon which
has been used for some preliminary studies of the tracker. So, when
considering adding a new module to the Gaspard Tracker, please do not use
-this class but create you own instead. However, for the explanations the
+this class but create your own instead. However, for the explanations the
GaspardTrackerDummyShape case will be considered.
When adding a new detector you need to follow several steps:
diff --git a/NPSimulation/include/Target.hh b/NPSimulation/include/Target.hh
index 5fcee069b..92da2be70 100644
--- a/NPSimulation/include/Target.hh
+++ b/NPSimulation/include/Target.hh
@@ -79,12 +79,13 @@ public:
public:
- G4double GetTargetThickness() {return m_TargetThickness;}
- 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 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;}
private:
diff --git a/NPSimulation/include/VEventGenerator.hh b/NPSimulation/include/VEventGenerator.hh
index 136487077..f5d040a10 100644
--- a/NPSimulation/include/VEventGenerator.hh
+++ b/NPSimulation/include/VEventGenerator.hh
@@ -11,7 +11,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : January 2009 *
- * Last update : *
+ * Last update : 03/11/2009 *
*---------------------------------------------------------------------------*
* Decription: *
* All event generator added in the project should derive from this virtual *
@@ -23,6 +23,8 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
+ * + 03/11/09: Adding EffectiveThiknessBeforeInteraction in the *
+ * CalculateBeamInteraction() method (N. de Sereville) *
* *
*****************************************************************************/
// C++ header
@@ -64,7 +66,8 @@ public:
double MeanPosY, double SigmaPosY, double MeanPosPhi, double SigmaPosPhi,
Target* target,
G4ThreeVector &InterCoord, double &AngleEmittanceTheta, double &AngleEmittancePhi,
- double &AngleIncidentTheta, double &AngleIncidentPhi);
+ 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 664ccf1d5..c0418d722 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -180,11 +180,13 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent, G4ParticleGun* particle
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);
+ Beam_theta, Beam_phi,
+ EffectiveThickness);
// write vertex position to ROOT file
G4double x0 = InterCoord.x();
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 87bba8e11..2f4d42251 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -30,6 +30,7 @@
// G4 headers
#include "G4ParticleTable.hh"
+#include "G4EmCalculator.hh"
#include "G4ParticleGun.hh"
#include "G4RotationMatrix.hh"
@@ -320,18 +321,25 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
//////////////////////////////////////////////////
//////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);
+
///////////////////////////////////////////////////////////////////////
///// Calculate the incident beam direction as well as the vertex /////
///// of interaction in target /////
@@ -339,11 +347,13 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
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);
+ Beam_theta, Beam_phi,
+ EffectiveThickness);
// write vertex position to ROOT file
G4double x0 = InterCoord.x();
@@ -380,9 +390,19 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
///// Angles for emitted particles following Cross Section //////
///// Angles are in the beam frame //////
/////////////////////////////////////////////////////////////////
- // Beam incident energy
+ // 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
+ const G4int NbLayers = 50;
+ 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
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index eeb2e8afc..790eb9322 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -397,6 +397,7 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
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);
@@ -405,7 +406,8 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
0, m_SigmaY, 0, m_SigmaPhiY,
m_Target,
InterCoord, Beam_thetaX, Beam_phiY,
- Beam_theta, Beam_phi);
+ Beam_theta, Beam_phi,
+ EffectiveThickness);
else
InterCoord = G4ThreeVector(0,0,0);
diff --git a/NPSimulation/src/VEventGenerator.cc b/NPSimulation/src/VEventGenerator.cc
index 6ab94ebb3..4dc2aac3d 100644
--- a/NPSimulation/src/VEventGenerator.cc
+++ b/NPSimulation/src/VEventGenerator.cc
@@ -9,7 +9,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : January 2009 *
- * Last update : *
+ * Last update : 03/11/2009 *
*---------------------------------------------------------------------------*
* Decription: *
* All event generator added in the project should derive from this virtual *
@@ -21,6 +21,8 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
+ * + 03/11/09: Adding EffectiveThiknessBeforeInteraction in the *
+ * CalculateBeamInteraction() method (N. de Sereville) *
* *
*****************************************************************************/
#include "VEventGenerator.hh"
@@ -28,8 +30,8 @@
// C++ headers
#include "cmath"
-// ROOT headers
-//#include"TMath.h"
+// G4 headers
+#include "G4UnitsTable.hh"
// CLHEP headers
#include "Randomize.hh"
@@ -51,7 +53,8 @@ void VEventGenerator::CalculateBeamInteraction(double MeanPosX, double SigmaPosX
double MeanPosY, double SigmaPosY, double MeanPosPhi, double SigmaPosPhi,
Target* target,
G4ThreeVector &InterCoord, double &AngleEmittanceTheta, double &AngleEmittancePhi,
- double &AngleIncidentTheta, double &AngleIncidentPhi)
+ double &AngleIncidentTheta, double &AngleIncidentPhi,
+ double &EffectiveTargetThicknessBeforeInteraction)
{
// target parameters
double TargetThickness = target->GetTargetThickness();
@@ -71,6 +74,7 @@ void VEventGenerator::CalculateBeamInteraction(double MeanPosX, double SigmaPosX
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
@@ -78,18 +82,6 @@ void VEventGenerator::CalculateBeamInteraction(double MeanPosX, double SigmaPosX
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);
@@ -99,6 +91,29 @@ void VEventGenerator::CalculateBeamInteraction(double MeanPosX, double SigmaPosX
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);
}
diff --git a/TODO b/TODO
index 5bb0b745c..a7a00e5e1 100644
--- a/TODO
+++ b/TODO
@@ -22,22 +22,21 @@
TODO for NPTool:
----------------
- + Build a dummy detector as a simple example of how NPTool is working (Adrien)
+ + Retrieve energyloss tables from G4 to be used for slowing the beam
+ and in the analysis.
+ Target work (Nicolas)
- + Make homogeneous treatment of materials definition wrt other detectors
- (Adrien)
+ Take correctly into account the emittance when the target is tilted.
+ IsInsideTarget method for debugging purposes
+ Add a dedicated class to deal with materials (see example from G4 tutorial)
(still under debate)
+ Split physic list and give the possibility to choose which package to use
+ Add support for messengers
+ + Scorers: CopyNumber v.s. DetectorNumber? (to be checked)
+ Build NPTool with the Autotool/Automake
+ Add documentation (Adrien + Nicolas)
TODO for NPAnalysis:
--------------------
- + Add calibration manager class from Julien (Adrien)
+ Use only one class for the analysis per detector instead of two currently
(ROOT feature in private member definition ( //!))
--
GitLab