diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector
index efc9692767cd8e6c6c691e153d9615d0757508b1..0c25c3e988018a4baad82e44162fc84aac8650e2 100644
--- a/Inputs/DetectorConfiguration/Riken_65mm.detector
+++ b/Inputs/DetectorConfiguration/Riken_65mm.detector
@@ -84,16 +84,25 @@ SILI= 0
CSI= 1
VIS= all
%%%%%%% Telescope 6 %%%%%%%
-M2Telescope
-X1_Y1= 155.77 50.23 8.18
-X1_Y128= 155.77 -50.23 8.18
-X128_Y1= 133.17 50.23 -89.7
-X128_Y128= 133.17 -50.23 -89.7
-SI= 1
-SILI= 1
-CSI= 0
-VIS= all
-
+%M2Telescope
+%X1_Y1= 155.77 50.23 8.18
+%X1_Y128= 155.77 -50.23 8.18
+%X128_Y1= 133.17 50.23 -89.7
+%X128_Y128= 133.17 -50.23 -89.7
+%SI= 1
+%SILI= 1
+%CSI= 0
+%VIS= all
+%% Alternate Telescope 6 %%
+M2Telescope
+ THETA= -130
+ PHI= 0
+ R= 150
+ BETA= 0 0 0
+ SI= 1
+ SILI= 1
+ CSI= 0
+ VIS= all
%%%%%%% Telescope 7 %%%%%%%
M2Telescope
X1_Y1= 27.07 50.23 -154.49
@@ -104,7 +113,7 @@ SI= 1
SILI= 1
CSI= 0
VIS= all
-
+
%%%%%%%%%%%%%%%%%%%%
AddThinSi
%%%%%%%%% Det 1 %%%%%%%%
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 8e2f05efa48fca6883c88dcdea97f4f0a89a992f..b840502107be8c8016561b15b8255f256248f175 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -2,7 +2,7 @@
%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Beam energy given in MeV ; Excitation in MeV
-Transfert
+TransfertToResonance
Beam= 11Li
Target= 2H
Light= 3He
@@ -10,13 +10,16 @@ Transfert
ExcitationEnergy= 1.0
BeamEnergy= 550
BeamEnergySpread= 0
- BeamFWHMX= 6.232
- BeamFWHMY= 9.069
- BeamEmmitanceTheta= 0.01208
- BeamEmmitancePhi= 0.01681
+ SigmaThetaX= 0.6921330164
+ SigmaPhiY= 0.963142053
+ SigmaX= 6.232
+ SigmaY= 9.069
+ ResonanceDecayZ= 2
+ ResonanceDecayA= 8
CrossSectionPath= 11Li(d,3He)10He.txt
- ShootLight= 1
+ ShootLight= 0
ShootHeavy= 1
+ ShootDecayProduct= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%0.69u
diff --git a/NPAnalysis/10He_Riken/Analysis b/NPAnalysis/10He_Riken/Analysis
index 6c0904b1217ac65c52f5f7efcf964871d826d228..7a15c0ec1ee2c7c99032279403bdbfc030232674 100755
Binary files a/NPAnalysis/10He_Riken/Analysis and b/NPAnalysis/10He_Riken/Analysis differ
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index a09db085a49b944b53913ca04d20a004871cd329..b51e713f1c3a6e1c1131ddf4364b57db83f84a74 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -68,16 +68,17 @@ int main(int argc,char** argv)
myDetector -> BuildPhysicalEvent() ;
E = M2 -> GetEnergyDeposit();
-
- XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1);
- YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
+ XTarget = Init->GetICPositionX(0);
+ YTarget = Init->GetICPositionY(0);
+// XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1);
+// YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
TVector3 HitDirection = M2 -> GetPositionOfInteraction() - TVector3(XTarget,YTarget,0);
- BeamTheta = RandomEngine.Gaus( Init->GetICIncidentAngleTheta(0)*deg , 2*deg ) ;
- BeamPhi = RandomEngine.Gaus( Init->GetICIncidentAnglePhi(0)*deg , 2*deg ) ;
+// BeamTheta = RandomEngine.Gaus( Init->GetICIncidentAngleTheta(0)*deg , 2*deg ) ;
+// BeamPhi = RandomEngine.Gaus( Init->GetICIncidentAnglePhi(0)*deg , 2*deg ) ;
-// BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
+ BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ;
// Angle between beam and particle
@@ -98,7 +99,7 @@ int main(int argc,char** argv)
20*micrometer , // Target Thickness at 0 degree
ThetaMM2Surface );
- // E = E + ThinSi ;
+// E = E + ThinSi ;
E= He3StripAl.EvaluateInitialEnergy( E , // Energy of the detected particle
0.4*micrometer , // Target Thickness at 0 degree
diff --git a/NPLib/MUST2/Must2Array.cxx b/NPLib/MUST2/Must2Array.cxx
index b627d4524f000f8bdea56abc9f951702fbae34bd..0c740c710ab3cef6cedfac91508b7f4aa3c8f533 100644
--- a/NPLib/MUST2/Must2Array.cxx
+++ b/NPLib/MUST2/Must2Array.cxx
@@ -239,9 +239,9 @@ void MUST2Array::ReadConfiguration(string Path)
//with angle method
else if ( check_Theta && check_Phi && check_R && check_beta )
{
- AddTelescope( R ,
- Theta ,
+ AddTelescope( Theta ,
Phi ,
+ R ,
beta_u ,
beta_v ,
beta_w );
@@ -504,14 +504,14 @@ void MUST2Array::AddTelescope( double theta ,
for( int j = 0 ; j < 128 ; j++ )
{
- X = C.X() + StripPitch * ( U.X()*i + V.X()*j ) ;
- Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j ) ;
- Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j ) ;
-
- lineX.push_back(X) ;
- lineY.push_back(Y) ;
- lineZ.push_back(Z) ;
-
+ X = C.X() + StripPitch * ( U.X()*i + V.X()*j ) ;
+ Y = C.Y() + StripPitch * ( U.Y()*i + V.Y()*j ) ;
+ Z = C.Z() + StripPitch * ( U.Z()*i + V.Z()*j ) ;
+
+ lineX.push_back(X) ;
+ lineY.push_back(Y) ;
+ lineZ.push_back(Z) ;
+
}
OneTelescopeStripPositionX.push_back(lineX) ;
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index 0d7bb48efeac3841eb9a95993746a5eb8f322b61..1c4d652cb23cfb52e65b9a18a897a16699146db0 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -23,7 +23,6 @@
// STL
#include <vector>
-
int main(int argc, char** argv)
{
diff --git a/NPSimulation/include/EventGeneratorTransfert.hh b/NPSimulation/include/EventGeneratorTransfert.hh
index f66562a6f5b39af8a53e8aab96e2a233597da910..09c3229bcd19f77956f9380ee79a085559c3d73d 100644
--- a/NPSimulation/include/EventGeneratorTransfert.hh
+++ b/NPSimulation/include/EventGeneratorTransfert.hh
@@ -26,8 +26,10 @@
// C++ headers
#include <string>
-// NPTool headers
+// NPSimulation
#include "VEventGenerator.hh"
+
+// NPLib
#include "TInitialConditions.h"
// NPLib header
@@ -52,8 +54,8 @@ class EventGeneratorTransfert : public VEventGenerator
double BeamEnergy , // Beam Energy
double ExcitationEnergy , // Excitation of Heavy Nuclei
double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
+ double SigmaX ,
+ double SigmaY ,
double SigmaThetaX ,
double SigmaPhiY ,
bool ShootLight ,
@@ -118,8 +120,8 @@ class EventGeneratorTransfert : public VEventGenerator
double BeamEnergy , // Beam Energy
double ExcitationEnergy , // Excitation of Heavy Nuclei
double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
+ double SigmaX ,
+ double SigmaY ,
double SigmaThetaX ,
double SigmaPhiY ,
bool ShootLight ,
diff --git a/NPSimulation/include/EventGeneratorTransfertToResonance.hh b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
index 312a3239c469882783f75322ac4b7fd6a7835a99..348ee5ed17289ce43132b4d0bc4e838d1c4ce997 100644
--- a/NPSimulation/include/EventGeneratorTransfertToResonance.hh
+++ b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
@@ -16,17 +16,22 @@
* Decription: *
* This event Generator is used to simulated two body TransfertReaction. *
* A Phase Space calculation is then performed to decay the Heavy product. *
+ * The TGenPhaseSpace from ROOT is used to calculate a phase space decay *
+ * with flat distribution *
*---------------------------------------------------------------------------*
* Comment: *
- * This class is not yet operationnal. *
+ * *
* *
*****************************************************************************/
// C++ header
#include <string>
-// NPTool header
+// NPSimulation
#include "VEventGenerator.hh"
+// NPLib
+#include "TInitialConditions.h"
+
// NPLib header
#include "NPReaction.h"
@@ -37,118 +42,114 @@ using namespace NPL ;
class EventGeneratorTransfertToResonance : public VEventGenerator
{
-public: // Constructors and Destructors
-
- // Default constructor used to allocate memory
- EventGeneratorTransfertToResonance();
-
- // This is the constructor to be used
- EventGeneratorTransfertToResonance(string name1 , //Beam nuclei
- string name2 , //Target nuclei
- string name3 , //Product of reaction
- string name4 , //Product of reaction
- double BeamEnergy , //Beam Energy
- double ExcitationEnergy , //Excitation of Heavy Nuclei
- double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path); //Path of the differentiel Cross Section
-
- // Default Destructor
- virtual ~EventGeneratorTransfertToResonance();
-
-public: // Inherit from VEventGenerator class
- void ReadConfiguration(string) ;
- void GenerateEvent(G4Event*, G4ParticleGun*) ;
-
- void SetTargetCoordinate(G4double X, G4double Y, G4double Z) {
- m_TargetX = X;
- m_TargetY = Y;
- m_TargetZ = Z;
- }
-
- void SetTargetThickness(double TargetThickness) {
- m_TargetThickness = TargetThickness ;
- }
-
- void SetTargetRadius(double TargetRadius) {
- m_TargetRadius = TargetRadius ;
- }
-
-private: // Particle Shoot Option
- bool m_ShootLight ;
- bool m_ShootHeavy ;
- bool m_ShootDecayProduct ;
-
-private: // TTree to store initial value of beam and reaction
- double m_InitialLightEnergy ;
- double m_InitialLightTheta ;
- double m_InitialBeamEnergy ;
- double m_InitialBeamTheta ;
- double m_InitialBeamX ;
- double m_InitialBeamY ;
- double m_InitialThetaCM ;
-
-private: // Beam Parameter
- double m_BeamEnergySpread ;
- double m_BeamFWHMX ;
- double m_BeamFWHMY ;
- double m_BeamEmmitanceTheta ;
- double m_BeamEmmitancePhi ;
-
-private: // Target Parameter
- double m_TargetThickness ;
- double m_TargetRadius ;
- double m_TargetX ;
- double m_TargetY ;
- double m_TargetZ ;
-
-private: // Reaction
- Reaction* m_Reaction ;
-
-private: // Resonance decay channel
- int m_ResonanceDecayZ ;
- int m_ResonanceDecayA ;
-
- //Other
- void Print() const ;
- void InitializeRootOutput() ;
- void ResonanceDecay(G4int parentZ ,
- G4int parentA ,
- G4double EnergyHeavy ,
- G4double ThetaHeavy ,
- G4double PhiHeavy ,
- G4double x0 ,
- G4double y0 ,
- G4double z0 ,
- G4Event* anEvent ,
- G4ParticleGun* particleGun);
-
-
- void SetEverything(string name1 , //Beam nuclei
- string name2 , //Target nuclei
- string name3 , //Product of reaction
- string name4 , //Product of reaction
- double BeamEnergy , //Beam Energy
- double ExcitationEnergy , //Excitation of Heavy Nuclei
- double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi ,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path); //Path of the differentiel Cross Section
+ public: // Constructors and Destructors
+
+ // Default constructor used to allocate memory
+ EventGeneratorTransfertToResonance();
+
+ // This is the constructor to be used
+ EventGeneratorTransfertToResonance(string name1 , //Beam nuclei
+ string name2 , //Target nuclei
+ string name3 , //Product of reaction
+ string name4 , //Product of reaction
+ double BeamEnergy , //Beam Energy
+ double ExcitationEnergy , //Excitation of Heavy Nuclei
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY,
+ int ResonanceDecayZ ,
+ int ResonanceDecayA ,
+ bool ShootLight ,
+ bool ShootHeavy ,
+ bool ShootDecayProduct ,
+ string Path); //Path of the differentiel Cross Section
+
+ // Default Destructor
+ virtual ~EventGeneratorTransfertToResonance();
+
+ public: // Inherit from VEventGenerator class
+ void ReadConfiguration(string) ;
+ void GenerateEvent(G4Event*, G4ParticleGun*) ;
+
+ void SetTargetCoordinate(G4double X, G4double Y, G4double Z) {
+ m_TargetX = X;
+ m_TargetY = Y;
+ m_TargetZ = Z;
+ }
+
+ void SetTargetThickness(double TargetThickness) {
+ m_TargetThickness = TargetThickness ;
+ }
+
+ void SetTargetRadius(double TargetRadius) {
+ m_TargetRadius = TargetRadius ;
+ }
+
+ private: // Particle Shoot Option
+ bool m_ShootLight ;
+ bool m_ShootHeavy ;
+ bool m_ShootDecayProduct ;
+
+ private: // TTree to store initial value of beam and reaction
+ TInitialConditions* m_InitConditions;
+
+ private: // Beam Parameter
+ double m_BeamEnergy ;
+ double m_BeamEnergySpread ;
+ double m_SigmaX ;
+ double m_SigmaY ;
+ double m_SigmaThetaX ;
+ double m_SigmaPhiY ;
+
+ private: // Target Parameter
+ double m_TargetThickness ;
+ double m_TargetRadius ;
+ double m_TargetX ;
+ double m_TargetY ;
+ double m_TargetZ ;
+
+ private: // Reaction
+ Reaction* m_Reaction ;
+
+ private: // Resonance decay channel
+ int m_ResonanceDecayZ ;
+ int m_ResonanceDecayA ;
+
+ //Other
+ void Print() const ;
+ void InitializeRootOutput() ;
+ void ResonanceDecay( G4double EnergyHeavy ,
+ G4double ThetaHeavy ,
+ G4double PhiHeavy ,
+ G4double x0 ,
+ G4double y0 ,
+ G4double z0 ,
+ G4Event* anEvent ,
+ G4ParticleGun* particleGun);
+
+ // This method return a random Vector of dimension N and magnitude R
+ // The return distribution populate uniformely the surface of the N-Sphere of radius R
+ vector<double> PhaseSpaceUniformGenerator( int N , double R);
+
+ void SetEverything(string name1 , //Beam nuclei
+ string name2 , //Target nuclei
+ string name3 , //Product of reaction
+ string name4 , //Product of reaction
+ double BeamEnergy , //Beam Energy
+ double ExcitationEnergy , //Excitation of Heavy Nuclei
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY ,
+ int ResonanceDecayZ ,
+ int ResonanceDecayA ,
+ bool ShootLight ,
+ bool ShootHeavy ,
+ bool ShootDecayProduct ,
+ string Path); //Path of the differentiel Cross Section
};
#endif
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 12a9205e554c38723eab62ef7e6bc12d8c834da3..569c489fc1633b003b7567abd2375d6f645c4a9a 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -464,9 +464,6 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
particleGun->GeneratePrimaryVertex(anEvent);
}
if (m_ShootHeavy) { // Case of recoil particle
- //
- // The case of recoil particle has not been tested with the new version of the event generator
- //
// Particle type
particleGun->SetParticleDefinition(HeavyName);
// Particle energy
@@ -484,10 +481,6 @@ 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
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
index 0a1fb9e2616f8f72f5b8940e8ea411fa4bf1d2ae..df57cbdde4252d73ca295d129cea40b15601682f 100644
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
@@ -14,9 +14,11 @@
* Decription: *
* This event Generator is used to simulated two body TransfertReaction. *
* A Phase Space calculation is then performed to decay the Heavy product. *
+ * The TGenPhaseSpace from ROOT is used to calculate a phase space decay *
+ * with flat distribution *
*---------------------------------------------------------------------------*
* Comment: *
- * This class is not yet operationnal. *
+ * *
* *
*****************************************************************************/
@@ -30,6 +32,7 @@
// G4 headers
#include "G4ParticleTable.hh"
#include "G4ParticleGun.hh"
+#include "G4RotationMatrix.hh"
// G4 headers including CLHEP headers
// for generating random numbers
@@ -39,6 +42,9 @@
#include "EventGeneratorTransfertToResonance.hh"
#include "RootOutput.h"
+//Root Headers
+#include "TGenPhaseSpace.h"
+
using namespace std;
using namespace CLHEP;
@@ -48,21 +54,21 @@ using namespace CLHEP;
EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance()
{
//------------- Default Constructor -------------
-
+ m_InitConditions = new TInitialConditions() ;
m_Reaction = new Reaction() ;
- m_BeamFWHMX = 0 ;
- m_BeamFWHMY = 0 ;
- m_BeamEmmitanceTheta = 0 ;
- m_BeamEmmitancePhi = 0 ;
- m_ResonanceDecayZ = 0 ;
- m_ResonanceDecayA = 0 ;
+ m_SigmaX = 0 ;
+ m_SigmaY = 0 ;
+ m_SigmaThetaX = 0 ;
+ m_SigmaPhiY = 0 ;
+ m_ResonanceDecayZ = 0 ;
+ m_ResonanceDecayA = 0 ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance()
{
//------------- Default Destructor ------------
-
+ delete m_InitConditions;
delete m_Reaction ;
}
@@ -74,10 +80,10 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string
double BeamEnergy ,
double ExcitationEnergy ,
double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY ,
int ResonanceDecayZ ,
int ResonanceDecayA ,
bool ShootLight ,
@@ -94,10 +100,10 @@ EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string
BeamEnergy , //Beam Energy
ExcitationEnergy , //Excitation of Heavy Nuclei
BeamEnergySpread ,
- BeamFWHMX ,
- BeamFWHMY ,
- BeamEmmitanceTheta ,
- BeamEmmitancePhi ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY ,
ResonanceDecayZ ,
ResonanceDecayA ,
ShootLight ,
@@ -113,13 +119,7 @@ void EventGeneratorTransfertToResonance::InitializeRootOutput()
{
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("ThetaCM" , &m_InitialThetaCM , "ThetaCM/D");
- pTree->Branch("InitialLightEnergy", &m_InitialLightEnergy , "InitialLightEnergy/D");
- pTree->Branch("InitialLightTheta" , &m_InitialLightTheta , "InitialLightTheta/D");
- pTree->Branch("InitialBeamEnergy" , &m_InitialBeamEnergy , "InitialBeamEnergy/D");
- pTree->Branch("InitialBeamTheta" , &m_InitialBeamTheta , "InitialBeamTheta/D");
- pTree->Branch("InitialBeamX" , &m_InitialBeamX , "InitialBeamX/D");
- pTree->Branch("InitialBeamY" , &m_InitialBeamY , "InitialBeamY/D");
+ pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
}
@@ -140,7 +140,7 @@ void EventGeneratorTransfertToResonance::ReadConfiguration(string Path)
////////Reaction Setting needs///////
string Beam, Target, Heavy, Light, CrossSectionPath ;
- G4double BeamEnergy = 0 , ExcitationEnergy = 0 , BeamEnergySpread = 0 , BeamFWHMX = 0 , BeamFWHMY = 0 , BeamEmmitanceTheta = 0 , BeamEmmitancePhi=0, ResonanceDecayZ = 0 , ResonanceDecayA = 0 ;
+ G4double BeamEnergy = 0 , ExcitationEnergy = 0 , BeamEnergySpread = 0 , SigmaX = 0 , SigmaY = 0 , SigmaThetaX = 0 , SigmaPhiY=0, ResonanceDecayZ = 0 , ResonanceDecayA = 0 ;
bool ShootLight = false ;
bool ShootHeavy = false ;
bool ShootDecayProduct = false ;
@@ -237,32 +237,32 @@ while(ReadingStatus){
G4cout << "Beam Energy Spread " << BeamEnergySpread / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "BeamFWHMX=") == 0) {
+ else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
check_FWHMX = true ;
ReactionFile >> DataBuffer;
- BeamFWHMX = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM X " << BeamFWHMX << " mm" << G4endl;
+ SigmaX = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Beam FWHM X " << SigmaX << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "BeamFWHMY=") == 0) {
+ else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
check_FWHMY = true ;
ReactionFile >> DataBuffer;
- BeamFWHMY = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM Y " << BeamFWHMX << " mm" << G4endl;
+ SigmaY = atof(DataBuffer.c_str()) * mm;
+ G4cout << "Beam FWHM Y " << SigmaX << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 19, "BeamEmmitanceTheta=") == 0) {
+ else if (DataBuffer.compare(0, 19, "SigmaThetaX=") == 0) {
check_EmmitanceTheta = true ;
ReactionFile >> DataBuffer;
- BeamEmmitanceTheta = atof(DataBuffer.c_str()) * rad;
- G4cout << "Beam Emmitance Theta " << BeamEmmitanceTheta / deg << " deg" << G4endl;
+ SigmaThetaX = atof(DataBuffer.c_str()) * rad;
+ G4cout << "Beam Emmitance Theta " << SigmaThetaX / deg << " deg" << G4endl;
}
- else if (DataBuffer.compare(0, 17, "BeamEmmitancePhi=") == 0) {
+ else if (DataBuffer.compare(0, 17, "SigmaPhiY=") == 0) {
check_EmmitancePhi = true ;
ReactionFile >> DataBuffer;
- BeamEmmitancePhi = atof(DataBuffer.c_str()) * rad;
- G4cout << "Beam Emmitance Phi " << BeamEmmitancePhi / deg << " deg" << G4endl;
+ SigmaPhiY = atof(DataBuffer.c_str()) * rad;
+ G4cout << "Beam Emmitance Phi " << SigmaPhiY / deg << " deg" << G4endl;
}
else if (DataBuffer.compare(0, 17, "ResonanceDecayZ=") == 0) {
@@ -335,10 +335,10 @@ while(ReadingStatus){
BeamEnergy ,
ExcitationEnergy ,
BeamEnergySpread ,
- BeamFWHMX ,
- BeamFWHMY ,
- BeamEmmitanceTheta ,
- BeamEmmitancePhi ,
+ SigmaX ,
+ SigmaY ,
+ SigmaThetaX ,
+ SigmaPhiY ,
ResonanceDecayZ ,
ResonanceDecayA ,
ShootLight ,
@@ -353,10 +353,12 @@ while(ReadingStatus){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
{
+ // Clear contents of Precedent event (now stored in ROOTOutput)
+ m_InitConditions->Clear();
+
//////////////////////////////////////////////////
//////Define the kind of particle to shoot////////
//////////////////////////////////////////////////
-
G4int LightZ = m_Reaction->GetNucleus3()->GetZ() ;
G4int LightA = m_Reaction->GetNucleus3()->GetA() ;
@@ -366,380 +368,346 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part
G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
- //////////////////////////////////////////////////
- /////Choose ThetaCM following Cross Section //////
- //////////////////////////////////////////////////
-
- //Calling RandGeneral fonction, which taking a double array as argument
- CLHEP::RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize() ) ;
- G4double ThetaCM = CrossSectionShoot.shoot() * (180 * deg) ;
- //G4double ThetaCM = 0;
- //Set the Theta angle of reaction
- m_InitialThetaCM = ThetaCM;
-
- //////////////////////////////////////////////////
- ////Momentum and direction set with kinematics////
- //////////////////////////////////////////////////
-
- //Variable where to store results.
- G4double ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy;
- //Compute Kinematic using previously defined ThetaCM
- m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
- //to write thetaCM in the root file
- m_InitialLightTheta = ThetaLight / deg ;
- m_InitialLightEnergy = EnergyLight / MeV ;
-
- // Vertex position and beam angle
- // 11Li Beam@Riken
- G4double x0 = 1000 * cm ;
- G4double y0 = 1000 * cm ;
- //shoot inside the target.
+ 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 Beam_thetaX = 0 ;
+ G4double Beam_phiY = 0 ;
+
+ //shoot inside the target with correlated angle
if (m_TargetRadius != 0) {
while (sqrt(x0*x0 + y0*y0) > m_TargetRadius) {
- x0 = RandGauss::shoot(0, m_BeamFWHMX / 2.35) ;
- y0 = RandGauss::shoot(0, m_BeamFWHMY / 2.35) ;
+ RandomGaussian2D(0, 0, m_SigmaX, m_SigmaThetaX, x0, Beam_thetaX);
+ RandomGaussian2D(0, 0, m_SigmaY, m_SigmaPhiY , y0, Beam_phiY );
}
}
-
else {
- x0 = 0 ;
- y0 = 0 ;
+ RandomGaussian2D(0,0,0,m_SigmaThetaX,x0,Beam_thetaX);
+ RandomGaussian2D(0,0,0,m_SigmaPhiY ,y0,Beam_phiY );
}
-
- // FIXME a changer pour prendre en compte correctement l'emmitance
- G4double Beam_theta = RandGauss::shoot(0, m_BeamEmmitanceTheta) ;
- // must shoot inside the target.
- G4double z0 = (-m_TargetThickness / 2 + CLHEP::RandFlat::shoot() * m_TargetThickness) ;
-
- // Move to the target
- x0 += m_TargetX ;
- y0 += m_TargetY ;
- z0 += m_TargetZ ;
+ // write emittance angles to ROOT file
+ m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
+ m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
- // Store initial value
- m_InitialBeamX = x0 ;
- m_InitialBeamY = y0 ;
- m_InitialBeamTheta = Beam_theta ;
+ // 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);
+
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+ G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ cos(Beam_theta) * sin(Beam_phi),
+ -sin(Beam_theta));
+ G4ThreeVector col2(-sin(Beam_phi),
+ cos(Beam_phi),
+ 0);
+ G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ sin(Beam_theta) * sin(Beam_phi),
+ cos(Beam_theta));
+ G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+ // Beam incident energy
+ G4double NominalBeamEnergy = m_BeamEnergy;
+ G4double IncidentBeamEnergy = RandGauss::shoot(NominalBeamEnergy, m_BeamEnergySpread / 2.35);
+ m_Reaction->SetBeamEnergy(IncidentBeamEnergy);
+ m_InitConditions->SetICIncidentEnergy(IncidentBeamEnergy / MeV);
+ // Angles
+ RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
+ G4double ThetaCM = CrossSectionShoot.shoot() * (180*deg);
+ G4double phi = RandFlat::shoot() * 2*pi;
+ // write angles to ROOT file
+ m_InitConditions->SetICEmittedAngleThetaCM(ThetaCM / deg);
+ m_InitConditions->SetICEmittedAnglePhiIncidentFrame(phi / deg);
//////////////////////////////////////////////////
- /////Now define everything for light particle/////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
//////////////////////////////////////////////////
- particleGun->SetParticleDefinition(LightName);
-
- G4double theta = ThetaLight + Beam_theta ;
- G4double phi = CLHEP::RandFlat::shoot() * 2 * pi ;
- G4double particle_energy = EnergyLight ;
-
- // Direction of particle, energy and laboratory angle
- G4double momentum_x = sin(theta) * cos(phi) ;
- G4double momentum_y = sin(theta) * sin(phi) ;
- G4double momentum_z = cos(theta) ;
+ // Variable where to store results
+ G4double ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy;
+ // Set the Theta angle of reaction
+ m_Reaction->SetThetaCM(ThetaCM);
+ // Compute Kinematic using previously defined ThetaCM
+ m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kineLight_beam(sin(ThetaLight) * cos(phi),
+ sin(ThetaLight) * sin(phi),
+ cos(ThetaLight));
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kineHeavy_beam(sin(ThetaHeavy) * cos(phi),
+ sin(ThetaHeavy) * sin(phi),
+ cos(ThetaHeavy));
+
+ // write angles/energy to ROOT file
+ m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(ThetaLight / deg);
+ m_InitConditions->SetICEmittedEnergy(EnergyLight/MeV);
+ //must shoot inside the target.
+ G4double z0 = (-m_TargetThickness / 2 + RandFlat::shoot() * m_TargetThickness);
- //Set the gun to shoot
- particleGun->SetParticleMomentumDirection(G4ThreeVector(momentum_x, momentum_y, momentum_z)) ;
- particleGun->SetParticleEnergy(particle_energy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Move to the target
+ x0 += m_TargetX ;
+ y0 += m_TargetY ;
+ z0 += m_TargetZ ;
- //Shoot the light particle
- if (m_ShootLight) particleGun->GeneratePrimaryVertex(anEvent) ;
+ // write vertex position to ROOT file
+ m_InitConditions->SetICPositionX(x0);
+ m_InitConditions->SetICPositionY(y0);
+ m_InitConditions->SetICPositionZ(z0);
//////////////////////////////////////////////////
- /////Now define everything for heavy particle/////
+ ///////// Set up everything for shooting /////////
//////////////////////////////////////////////////
-
- theta = ThetaHeavy + Beam_theta ;
- particle_energy = EnergyHeavy ;
-
- if (m_ShootHeavy) ResonanceDecay( HeavyZ ,
- HeavyA ,
- EnergyHeavy ,
- theta ,
- phi ,
- x0 ,
- y0 ,
- z0 ,
- anEvent ,
- particleGun );
+ if (m_ShootLight) { // Case of light particle
+ // Particle type
+ particleGun->SetParticleDefinition(LightName);
+ // Particle energy
+ particleGun->SetParticleEnergy(EnergyLight);
+ // Particle vertex position
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
+ // Particle direction
+ // Kinematical angles in the beam frame are transformed
+ // to the "world" frame
+ G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineLight_beam;
+ // get theta and phi in the world frame
+ G4double theta_world = momentum_kine_world.theta();
+ G4double phi_world = momentum_kine_world.phi();
+ if (phi_world < 1e-6) phi_world += 2*pi;
+ // write angles in ROOT file
+ m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
+ m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
+ //Set the gun to shoot
+ particleGun->SetParticleMomentumDirection(momentum_kine_world);
+ //Shoot the light particle
+ particleGun->GeneratePrimaryVertex(anEvent);
+ }
+
+ if (m_ShootHeavy) { // Case of recoil particle
+ // Particle type
+ particleGun->SetParticleDefinition(HeavyName);
+ // Particle energy
+ particleGun->SetParticleEnergy(EnergyHeavy);
+ // Particle vertex position
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
+ // Particle direction
+ // Kinematical angles in the beam frame are transformed
+ // to the "world" frame
+ G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineHeavy_beam;
+ // get theta and phi in the world frame
+ G4double theta_world = momentum_kine_world.theta();
+ G4double phi_world = momentum_kine_world.phi();
+ if (phi_world < 1e-6) phi_world += 2*pi;
+
+ EventGeneratorTransfertToResonance::ResonanceDecay( EnergyHeavy ,
+ theta_world ,
+ phi_world ,
+ x0 ,
+ y0 ,
+ z0 ,
+ anEvent ,
+ particleGun );
+
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::ResonanceDecay(G4int parentZ ,
- G4int parentA ,
- G4double EnergyHeavy ,
- G4double ThetaHeavy ,
- G4double PhiHeavy ,
- G4double x0 ,
- G4double y0 ,
- G4double z0 ,
- G4Event* anEvent ,
- G4ParticleGun* particleGun)
+void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy ,
+ G4double ThetaHeavy ,
+ G4double PhiHeavy ,
+ G4double x0 ,
+ G4double y0 ,
+ G4double z0 ,
+ G4Event* anEvent ,
+ G4ParticleGun* particleGun )
{
- //FIXME
+ G4double parentZ = m_Reaction->GetNucleus4()->GetZ() ;
+ G4double parentA = m_Reaction->GetNucleus4()->GetA() ;
+
G4int NumberOfNeutrons = (parentA - parentZ) - (m_ResonanceDecayA - m_ResonanceDecayZ) ;
G4int NumberOfProtons = parentZ - m_ResonanceDecayZ ;
G4int NumberOfDecayProducts = NumberOfNeutrons + NumberOfProtons ;
- /*G4int NumberOfNeutrons = 1 ;
- G4int NumberOfProtons = 0 ;
- G4int NumberOfDecayProducts = NumberOfNeutrons + NumberOfProtons ;*/
-
if (NumberOfNeutrons < 0 || NumberOfProtons < 0) {
G4cout << "Error input for Resonance decay" << G4endl;
return;
}
-
+
else {
- //Obtain Mass of daughter Nuclei
- G4ParticleDefinition* parent
- = G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, m_Reaction->GetExcitation()) ;
- G4ParticleDefinition* daughter
- = G4ParticleTable::GetParticleTable()->GetIon(m_ResonanceDecayZ, m_ResonanceDecayA, 0.) ;
- G4ParticleDefinition* neutron
- = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
- G4ParticleDefinition* proton
- = G4ParticleTable::GetParticleTable()->FindParticle("proton");
-
- G4double M = parent -> GetPDGMass() ;
- G4double md = daughter -> GetPDGMass() ;
-
- G4double mn = neutron -> GetPDGMass() ;
- G4double mp = proton -> GetPDGMass() ;
-
- // FIXME Hexaneutron
- /* mn = 6*(neutron-> GetPDGMass())-100*keV ;
- G4double E = (M*M - mn*mn + md*md) / (2*M) ;
- G4double p = sqrt(( M*M - (md + mn)*(md + mn) ) * ( M*M - (md - mn)*(md - mn) )) / (2*M) ;*/
-
- // G4double cosThetaCM = CLHEP::RandFlat::shoot()*2-1 ;
- // G4double ThetaCM = acos(cosThetaCM) ;
- // G4double PhiCM = CLHEP::RandFlat::shoot() * 2 * pi ;
-
- // FIXME Hexaneutron
- /*
- G4double DaughterMomentumX = sin(ThetaCM)*cos(PhiCM) ;
- G4double DaughterMomentumY = sin(ThetaCM)*sin(PhiCM) ;
- G4double DaughterMomentumZ = cos(ThetaCM) ;*/
-
- G4double Q = m_Reaction->GetExcitation() + M - md - mn * NumberOfNeutrons - mp * NumberOfProtons ;
-
- vector<G4double> DecayProductsMomentumCM ;
- vector<G4double> DecayProductsMomentumXCM ;
- vector<G4double> DecayProductsMomentumYCM ;
- vector<G4double> DecayProductsMomentumZCM ;
- vector<G4double> DecayProductsThetaCM ;
- vector<G4double> DecayProductsPhiCM ;
-
- G4double AvaibleEnergy = Q ;
- //Calculate Energy of Heavy
- G4double MomentumCMHeavy = CLHEP::RandFlat::shoot() * AvaibleEnergy ;
- AvaibleEnergy = AvaibleEnergy - MomentumCMHeavy ;
-
- G4double MomentumCM = 0 ;
- // Calculate DecayProducts Properties
- for (G4int i = 0 ; i < NumberOfDecayProducts ; i++) {
- if (i != NumberOfDecayProducts - 1) MomentumCM = CLHEP::RandFlat::shoot() * AvaibleEnergy ;
- else MomentumCM = AvaibleEnergy ;
-
- AvaibleEnergy = AvaibleEnergy - MomentumCM ;
- DecayProductsMomentumCM.push_back(MomentumCM) ;
-
- G4double cosThetaCM = CLHEP::RandFlat::shoot() * 2 - 1 ;
- G4double ThetaCM = acos(cosThetaCM) ;
- DecayProductsThetaCM.push_back(ThetaCM) ;
-
- G4double PhiCM = CLHEP::RandFlat::shoot() * 2 * pi ;
- DecayProductsPhiCM.push_back(PhiCM) ;
-
- DecayProductsMomentumXCM.push_back(sin(ThetaCM)*cos(PhiCM)) ;
- DecayProductsMomentumYCM.push_back(sin(ThetaCM)*sin(PhiCM)) ;
- DecayProductsMomentumZCM.push_back(cos(ThetaCM)) ;
- }
-
- // Applied conservation of Momentum to daughter nuclei
- G4double DaughterMomentumX = 0 ;
- G4double DaughterMomentumY = 0 ;
- G4double DaughterMomentumZ = 0 ;
-
- for (G4int i = 0 ; i < NumberOfDecayProducts ; i++) {
- DaughterMomentumX = (DaughterMomentumX - DecayProductsMomentumCM[i] * DecayProductsMomentumXCM[i]) ;
- DaughterMomentumY = (DaughterMomentumY - DecayProductsMomentumCM[i] * DecayProductsMomentumYCM[i]) ;
- DaughterMomentumZ = (DaughterMomentumZ - DecayProductsMomentumCM[i] * DecayProductsMomentumZCM[i]) ;
- }
-
- // Daughter MomentumCM
- G4double p =
- sqrt(DaughterMomentumX * DaughterMomentumX + DaughterMomentumY * DaughterMomentumY + DaughterMomentumZ * DaughterMomentumZ) ;
-
- // Daughter EnergyCm
- G4double E = sqrt(p * p + md * md) ;
- DaughterMomentumX = DaughterMomentumX / p ;
- DaughterMomentumY = DaughterMomentumY / p ;
- DaughterMomentumZ = DaughterMomentumZ / p ;
-
- // CM to Lab //
-
- // Initial Lab to CM Momentum
- G4double InitialE = sqrt(EnergyHeavy * EnergyHeavy + M * M) ;
- G4double InitialMomentum = EnergyHeavy ;
- G4double InitialMomentumX = sin(ThetaHeavy) * cos(PhiHeavy) ;
- G4double InitialMomentumY = sin(ThetaHeavy) * sin(PhiHeavy) ;
- G4double InitialMomentumZ = cos(ThetaHeavy) ;
-
- // Beta and Gamma CM to Lab
- /* G4double betaX = (DaughterMomentumX*p - InitialMomentumX*InitialMomentum)/E ;
- G4double betaY = (DaughterMomentumY*p - InitialMomentumY*InitialMomentum)/E ;
- G4double betaZ = (DaughterMomentumZ*p - InitialMomentumZ*InitialMomentum)/E ;
- G4double beta = sqrt (betaX*betaX + betaY*betaY + betaZ*betaZ ) ;
- G4double beta2 = beta*beta ;
- G4double gamma = 1 / ( sqrt(1 - beta2 ) ) ;*/
-
- G4double betaX = -(InitialMomentumX * InitialMomentum) / InitialE ;
- G4double betaY = -(InitialMomentumY * InitialMomentum) / InitialE ;
- G4double betaZ = -(InitialMomentumZ * InitialMomentum) / InitialE ;
- G4double beta = sqrt(betaX * betaX + betaY * betaY + betaZ * betaZ) ;
- G4double beta2 = beta * beta ;
- G4double gamma = 1 / (sqrt(1 - beta2)) ;
- // Calculate everything for heavy particule
-
- /* G4double NewEnergy =
- gamma*E
- - betaX*gamma*DaughterMomentumX*p
- - betaY*gamma*DaughterMomentumY*p
- - betaZ*gamma*DaughterMomentumZ*p;*/
-
- G4double NewMomentumX =
- - betaX * gamma * E
- + DaughterMomentumX * p + (gamma - 1) * (betaX * betaX) / (beta2) * DaughterMomentumX * p
- + (gamma - 1) * (betaX * betaY) / (beta2) * DaughterMomentumY * p
- + (gamma - 1) * (betaX * betaZ) / (beta2) * DaughterMomentumZ * p ;
-
- G4double NewMomentumY =
- - betaY * gamma * E
- + (gamma - 1) * (betaY * betaX) / (beta2) * DaughterMomentumX * p
- + DaughterMomentumY * p + (gamma - 1) * (betaY * betaY) / (beta2) * DaughterMomentumY * p
- + (gamma - 1) * (betaY * betaZ) / (beta2) * DaughterMomentumZ * p ;
-
- G4double NewMomentumZ =
- - betaZ * gamma * E
- + (gamma - 1) * (betaZ * betaX) / (beta2) * DaughterMomentumX * p
- + (gamma - 1) * (betaZ * betaY) / (beta2) * DaughterMomentumY * p
- + DaughterMomentumZ * p + (gamma - 1) * (betaZ * betaZ) / (beta2) * DaughterMomentumZ * p ;
-
- G4double
- NewEnergy = sqrt(NewMomentumX * NewMomentumX + NewMomentumY * NewMomentumY + NewMomentumZ * NewMomentumZ) ;
- NewMomentumX = NewMomentumX / NewEnergy ;
- NewMomentumY = NewMomentumY / NewEnergy ;
- NewMomentumZ = NewMomentumZ / NewEnergy ;
-
- //Set the gun to shoot
- particleGun->SetParticleDefinition(daughter) ;
- particleGun->SetParticleMomentumDirection(G4ThreeVector(NewMomentumX, NewMomentumY, NewMomentumZ)) ;
- particleGun->SetParticleEnergy(NewEnergy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
-
- // Shoot the Daugter
- particleGun->GeneratePrimaryVertex(anEvent) ;
-
- if (m_ShootDecayProduct) {
- G4int jj = 0 ;
- for (; jj < NumberOfNeutrons ; jj++) {
- p = sqrt(DecayProductsMomentumXCM[jj] * DecayProductsMomentumXCM[jj]
- + DecayProductsMomentumYCM[jj] * DecayProductsMomentumYCM[jj]
- + DecayProductsMomentumZCM[jj] * DecayProductsMomentumZCM[jj]) ;
-
- E = sqrt(p * p + mn * mn) ;
-
- NewMomentumX =
- - betaX * gamma * E
- + DecayProductsMomentumXCM[jj] * p + (gamma - 1) * (betaX * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaX * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaX * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewMomentumY =
- - betaY * gamma * E
- + (gamma - 1) * (betaY * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + DecayProductsMomentumYCM[jj] * p + (gamma - 1) * (betaY * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaY * betaZ) / (beta2) * DecayProductsMomentumYCM[jj] * p ;
-
- NewMomentumZ =
- - betaZ * gamma * E
- + (gamma - 1) * (betaZ * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaZ * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + DecayProductsMomentumZCM[jj] * p + (gamma - 1) * (betaZ * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewEnergy = sqrt(NewMomentumX * NewMomentumX + NewMomentumY * NewMomentumY + NewMomentumZ * NewMomentumZ) ;
- NewMomentumX = NewMomentumX / NewEnergy ;
- NewMomentumY = NewMomentumY / NewEnergy ;
- NewMomentumZ = NewMomentumZ / NewEnergy ;
- //Set the gun to shoot
- particleGun->SetParticleDefinition(neutron) ;
- particleGun->SetParticleMomentumDirection(G4ThreeVector(NewMomentumX, NewMomentumY, NewMomentumZ)) ;
- particleGun->SetParticleEnergy(NewEnergy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
- particleGun->GeneratePrimaryVertex(anEvent) ;
- }
-
- for (; jj < NumberOfProtons ; jj++) {
- p = sqrt(DecayProductsMomentumXCM[jj] * DecayProductsMomentumXCM[jj]
- + DecayProductsMomentumYCM[jj] * DecayProductsMomentumYCM[jj]
- + DecayProductsMomentumZCM[jj] * DecayProductsMomentumZCM[jj]) ;
-
- E = sqrt(p * p + mp * mp) ;
-
- NewMomentumX =
- - betaX * gamma * E
- + DecayProductsMomentumXCM[jj] * p + (gamma - 1) * (betaX * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaX * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaX * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewMomentumY =
- - betaY * gamma * E
- + (gamma - 1) * (betaY * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + DecayProductsMomentumYCM[jj] * p + (gamma - 1) * (betaY * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + (gamma - 1) * (betaY * betaZ) / (beta2) * DecayProductsMomentumYCM[jj] * p ;
-
- NewMomentumZ =
- - betaZ * gamma * E
- + (gamma - 1) * (betaZ * betaX) / (beta2) * DecayProductsMomentumXCM[jj] * p
- + (gamma - 1) * (betaZ * betaY) / (beta2) * DecayProductsMomentumYCM[jj] * p
- + DecayProductsMomentumZCM[jj] * p + (gamma - 1) * (betaZ * betaZ) / (beta2) * DecayProductsMomentumZCM[jj] * p ;
-
- NewEnergy = sqrt(NewMomentumX * NewMomentumX + NewMomentumY * NewMomentumY + NewMomentumZ * NewMomentumZ) ;
- NewMomentumX = NewMomentumX / NewEnergy ;
- NewMomentumY = NewMomentumY / NewEnergy ;
- NewMomentumZ = NewMomentumZ / NewEnergy ;
- //Set the gun to shoot
- particleGun->SetParticleDefinition(neutron) ;
- particleGun->SetParticleMomentumDirection(G4ThreeVector(NewMomentumX, NewMomentumY, NewMomentumZ)) ;
- particleGun->SetParticleEnergy(NewEnergy) ;
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
- particleGun->GeneratePrimaryVertex(anEvent) ;
- }
- }
- }
+ //Obtain Mass of daughter Nuclei
+ G4ParticleDefinition* parent
+ = G4ParticleTable::GetParticleTable()->GetIon(parentZ, parentA, m_Reaction->GetExcitation()) ;
+
+ G4ParticleDefinition* daughter
+ = G4ParticleTable::GetParticleTable()->GetIon(m_ResonanceDecayZ, m_ResonanceDecayA, 0.) ;
+
+ G4ParticleDefinition* neutron
+ = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+
+ G4ParticleDefinition* proton
+ = G4ParticleTable::GetParticleTable()->FindParticle("proton");
+
+ G4double M = parent -> GetPDGMass() ;
+ G4double md = daughter -> GetPDGMass() ;
+ G4double mn = neutron -> GetPDGMass() ;
+ G4double mp = proton -> GetPDGMass() ;
+
+ G4double Q = M - md - mn * NumberOfNeutrons - mp * NumberOfProtons ;
+
+ vector<G4double> DecayProductsMomentumCM ;
+ vector<G4double> DecayProductsMomentumXCM ;
+ vector<G4double> DecayProductsMomentumYCM ;
+ vector<G4double> DecayProductsMomentumZCM ;
+ vector<G4double> DecayProductsThetaCM ;
+ vector<G4double> DecayProductsPhiCM ;
+
+ // Initial Lab Momentum
+ G4double InitialE = sqrt(EnergyHeavy * EnergyHeavy + M * M) ;
+ G4double InitialMomentumX = EnergyHeavy * sin(ThetaHeavy) * cos(PhiHeavy) ;
+ G4double InitialMomentumY = EnergyHeavy * sin(ThetaHeavy) * sin(PhiHeavy) ;
+ G4double InitialMomentumZ = EnergyHeavy * cos(ThetaHeavy) ;
+
+ TLorentzVector Initial = TLorentzVector(InitialMomentumX/GeV, InitialMomentumY/GeV, InitialMomentumZ/GeV,InitialE/GeV);
+
+ // Array of masses express in GeV/c2
+ double* masses = new double[NumberOfDecayProducts+1];
+
+ // Filling Array
+ masses[0] = md/GeV ;
+
+ int ll = 1 ;
+ for(int i = 0 ; i < NumberOfNeutrons ; i++)
+ {masses[ll] = mn/GeV ; ll++;}
+
+ for(int i = 0 ; i < NumberOfProtons ; i++)
+ {masses[ll] = mp/GeV ; ll++;}
+
+ // Instentiate a Phase Space Generator, with flat distrution
+
+
+ TGenPhaseSpace TPhaseSpace ;
+
+ if( !TPhaseSpace.SetDecay(Initial, NumberOfDecayProducts+1, masses,"Fermi") ) cout << "Warning: Phase Space Decay forbiden by kinematic";
+ TPhaseSpace.Generate() ;
+
+ TLorentzVector* daugterLV = TPhaseSpace.GetDecay(0);
+ G4ThreeVector Momentum = G4ThreeVector( daugterLV->X()*GeV ,
+ daugterLV->Y()*GeV ,
+ daugterLV->Z()*GeV );
+ double Energy = Momentum.mag() ;
+ Momentum.unit() ;
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(daughter) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+
+ if (m_ShootDecayProduct)
+ {
+ G4int jj = 1 ;
+ for ( int u = 0; u < NumberOfNeutrons ; u++)
+ {
+ TLorentzVector* neutronLV = TPhaseSpace.GetDecay(jj);
+
+ Momentum = G4ThreeVector( daugterLV->X()*GeV ,
+ daugterLV->Y()*GeV ,
+ daugterLV->Z()*GeV );
+ Energy = Momentum.mag() ;
+ Momentum.unit() ;
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(neutron) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+ jj++;
+ }
+
+ for ( int u = 0; u < NumberOfProtons ; u++)
+ {
+ TLorentzVector* protonLV = TPhaseSpace.GetDecay(jj);
+
+ Momentum = G4ThreeVector( daugterLV->X()*GeV ,
+ daugterLV->Y()*GeV ,
+ daugterLV->Z()*GeV );
+ Energy = Momentum.mag() ;
+ Momentum.unit() ;
+ //Set the gun to shoot
+ particleGun->SetParticleDefinition(proton) ;
+ particleGun->SetParticleMomentumDirection(Momentum) ;
+ particleGun->SetParticleEnergy(Energy) ;
+ particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)) ;
+ // Shoot the Daugter
+ particleGun->GeneratePrimaryVertex(anEvent) ;
+ jj++;
+ }
+
+ delete masses ;
+ }
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::SetEverything(string name1 ,
- string name2 ,
- string name3 ,
- string name4 ,
- double BeamEnergy ,
- double Excitation ,
- double BeamEnergySpread ,
- double BeamFWHMX ,
- double BeamFWHMY ,
- double BeamEmmitanceTheta ,
- double BeamEmmitancePhi ,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path)
+vector<double> EventGeneratorTransfertToResonance::PhaseSpaceUniformGenerator( int N , double R)
+ {
+ vector<double> V ;
+ V.reserve(N) ;
+ double Norme ;
+ double Buffer ;
+
+ for(int i = 0 ; i< N ; i++)
+ {
+ V.push_back( Buffer = CLHEP::RandFlat::shoot() );
+ Norme += Buffer*Buffer ;
+ }
+
+ Norme = sqrt(Norme) ;
+
+ for(int i = 0 ; i< N ; i++)
+ V[i] = V[i] / Norme ;
+
+ return V;
+ }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorTransfertToResonance::SetEverything( string name1 ,
+ string name2 ,
+ string name3 ,
+ string name4 ,
+ double BeamEnergy ,
+ double Excitation ,
+ double BeamEnergySpread ,
+ double SigmaX ,
+ double SigmaY ,
+ double SigmaThetaX ,
+ double SigmaPhiY ,
+ int ResonanceDecayZ ,
+ int ResonanceDecayA ,
+ bool ShootLight ,
+ bool ShootHeavy ,
+ bool ShootDecayProduct ,
+ string Path )
{
//------------- Constructor with nuclei names and beam energy ------------
@@ -750,12 +718,12 @@ void EventGeneratorTransfertToResonance::SetEverything(string name1
BeamEnergy,
Excitation,
Path) ;
-
+ m_BeamEnergy = BeamEnergy;
m_BeamEnergySpread = BeamEnergySpread ;
- m_BeamFWHMX = BeamFWHMX ;
- m_BeamFWHMY = BeamFWHMY ;
- m_BeamEmmitanceTheta = BeamEmmitanceTheta ;
- m_BeamEmmitancePhi = BeamEmmitancePhi ;
+ m_SigmaX = SigmaX ;
+ m_SigmaY = SigmaY ;
+ m_SigmaThetaX = SigmaThetaX ;
+ m_SigmaPhiY = SigmaPhiY ;
m_ResonanceDecayZ = ResonanceDecayZ ;
m_ResonanceDecayA = ResonanceDecayA ;
m_ShootLight = ShootLight ;
@@ -764,3 +732,5 @@ void EventGeneratorTransfertToResonance::SetEverything(string name1
}
+
+
diff --git a/NPSimulation/vis.mac b/NPSimulation/vis.mac
index 7130261ec20bc96c8fddd2d03740a8b7316d5944..e8345fb35fdaef1effc82ccf8f853be6f7e0cc3b 100644
--- a/NPSimulation/vis.mac
+++ b/NPSimulation/vis.mac
@@ -9,7 +9,7 @@
# choose a graphic system
#/vis/open OGLIX
#/vis/open OGLSX
-#/vis/open VRML2FILE
+/vis/open VRML2FILE
/vis/scene/create
/vis/drawVolume
/vis/viewer/set/viewpointThetaPhi 0 0 deg