From ce1e792c93754f385d96ae940f25dea6280fb1a7 Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Wed, 18 Jul 2012 13:23:13 +0000
Subject: [PATCH] Change for new NPReaction Class by Pierre
---
.../TComptonTelescopeDataDict.cxx | 2 +-
.../TComptonTelescopeProcessDataDict.cxx | 2 +-
NPLib/Eurogam/TEurogamDataDict.cxx | 2 +-
NPLib/Physics/NPNucleus.cxx | 108 ++-
NPLib/Physics/NPNucleus.h | 97 ++-
NPLib/Physics/NPReaction.cxx | 415 ++++++---
NPLib/Physics/NPReaction.h | 163 ++--
NPLib/Vamos/TVamosCHIODataDict.cxx | 2 +-
NPLib/Vamos/TVamosDCDataDict.cxx | 2 +-
NPLib/Vamos/TVamosFingerDataDict.cxx | 2 +-
NPLib/Vamos/TVamosPlasticDataDict.cxx | 2 +-
.../EventGeneratorTransfertToResonance.hh | 151 ----
NPSimulation/src/EventGeneratorTransfert.cc | 42 +-
.../src/EventGeneratorTransfertToResonance.cc | 823 ------------------
14 files changed, 590 insertions(+), 1223 deletions(-)
delete mode 100644 NPSimulation/include/EventGeneratorTransfertToResonance.hh
delete mode 100644 NPSimulation/src/EventGeneratorTransfertToResonance.cc
diff --git a/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx b/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx
index 395a14e4c..0b23325ed 100644
--- a/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx
+++ b/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:43:12 2012
+// File generated by rootcint at Wed Jul 18 14:41:32 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx b/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx
index 1ed66ad98..7e2d1a123 100644
--- a/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx
+++ b/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:43:14 2012
+// File generated by rootcint at Wed Jul 18 14:41:34 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPLib/Eurogam/TEurogamDataDict.cxx b/NPLib/Eurogam/TEurogamDataDict.cxx
index 688f67b98..670ac9d75 100644
--- a/NPLib/Eurogam/TEurogamDataDict.cxx
+++ b/NPLib/Eurogam/TEurogamDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:43:17 2012
+// File generated by rootcint at Wed Jul 18 14:41:36 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPLib/Physics/NPNucleus.cxx b/NPLib/Physics/NPNucleus.cxx
index 5dec2f6ca..bc4d608ab 100644
--- a/NPLib/Physics/NPNucleus.cxx
+++ b/NPLib/Physics/NPNucleus.cxx
@@ -9,7 +9,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : febuary 2009 *
- * Last update : *
+ * Last update : may 2012 morfouac@ipno.in2p3.fr *
*---------------------------------------------------------------------------*
* Decription: *
* This class manage a nucleus, data are read in the nubtab03.asc file *
@@ -24,9 +24,20 @@
// C++ headers
#include <iostream>
#include <fstream>
-#include <stdlib.h>
+#include <string>
+#include <cstdlib>
+#include <sstream>
+//#include <stdlib.h>
+#include <cmath>
+#include <vector>
+
#include "NPNucleus.h"
+#include "Constant.h"
+
+// Use CLHEP System of unit and Physical Constant
+//#include "CLHEP/Units/GlobalSystemOfUnits.h"
+//#include "CLHEP/Units/PhysicalConstants.h"
using namespace std;
using namespace NPL;
@@ -118,7 +129,7 @@ void Nucleus::Extract(const char* line)
string ligne = line;
// name of the isotope
- string s_name = ligne.substr(11,6);
+ string s_name = ligne.substr(11,7);
fName = s_name.data();
// charge and mass
@@ -165,6 +176,9 @@ void Nucleus::Extract(const char* line)
if (s_spinparity.find("17/2") != string::npos) fSpin = 8.5 ;
if (s_spinparity.find("19/2") != string::npos) fSpin = 9.5 ;
if (s_spinparity.find("21/2") != string::npos) fSpin = 10.5 ;
+ //cout << fName << endl;
+ GetNucleusName();
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -178,11 +192,93 @@ void Nucleus::Print() const
cout << "Jpi = " << fSpinParity << " (" << fSpin << fParity << ")" << endl;
}
-string Nucleus::GetName() const
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::GetNucleusName()
+{
+fNucleusName=string(fName,6);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::EnergyToBrho()
+{
+ fBrho = sqrt(pow(fKineticEnergy,2) + 2*fKineticEnergy*Mass()) * 1e6 * echarge / C / (GetZ() * echarge);
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::EnergyToTof() // second/meter
+{
+ fTimeOfFlight = 1/sqrt(1-(Mass()*Mass())/(fKineticEnergy+Mass())/(fKineticEnergy+Mass()))/C;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::TofToEnergy()
+{
+ double Energy = sqrt( Mass()*Mass()/(1-pow((1/(C*fTimeOfFlight)),2)) );
+ fKineticEnergy = Energy - Mass();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::BrhoToEnergy()
+{
+ fKineticEnergy = sqrt( pow((fBrho*C*GetZ()*echarge)/(1e6*echarge),2) + pow(Mass(),2) ) - Mass();
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::EnergyToBeta()
+{
+ fBeta = sqrt(pow(fKineticEnergy,2) + 2*fKineticEnergy*Mass())/(fKineticEnergy + Mass());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::BetaToEnergy()
+{
+ fKineticEnergy = Mass()/sqrt(1-pow(fBeta,2)) - Mass();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::BetaToGamma()
+{
+ fGamma = 1/sqrt(1-pow(fBeta,2));
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Nucleus::BetaToVelocity() // in cm/ns
{
-string Name=string(fName,1);
-return Name;
+ fVelocity = C*fBeta*1e-7;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+double Nucleus::DopplerCorrection(double EnergyLabGamma, double ThetaLabGamma)
+{
+ double EnergyGammaCorrected = EnergyLabGamma*(1-fBeta*cos(ThetaLabGamma))/( sqrt(1-pow(fBeta,2)) );
+
+ return EnergyGammaCorrected;
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/NPLib/Physics/NPNucleus.h b/NPLib/Physics/NPNucleus.h
index a8e5f0001..d60dedf7d 100644
--- a/NPLib/Physics/NPNucleus.h
+++ b/NPLib/Physics/NPNucleus.h
@@ -24,6 +24,7 @@
*****************************************************************************/
#define uma 931.49432
#include <string>
+#include "TLorentzVector.h"
using namespace std;
namespace NPL
@@ -37,34 +38,80 @@ namespace NPL
~Nucleus();
private :
- const char* fName; // Nucleus name
- int fCharge; // Nucleus charge
- int fAtomicWeight; // Nucleus atomic weight
- double fMassExcess; // Nucleus mass excess in keV
- const char* fSpinParity; // Nucleus spin and parity
- double fSpin; // Nucleus spin
- const char* fParity; // Nucleus parity
-
+ //intrinsic properties
+ const char* fName; // Nucleus name
+ string fNucleusName;
+ int fCharge; // Nucleus charge
+ int fAtomicWeight; // Nucleus atomic weight
+ double fMassExcess; // Nucleus mass excess in keV
+ const char* fSpinParity; // Nucleus spin and parity
+ double fSpin; // Nucleus spin
+ const char* fParity; // Nucleus parity
+ //kinematic properties
+ double fKineticEnergy;
+ double fBeta;
+ double fGamma;
+ double fBrho;
+ double fTimeOfFlight;
+ double fVelocity;
+ TLorentzVector fEnergyImpulsion;
+
+ public:
+ void EnergyToBrho();
+ void EnergyToTof();
+ void BetaToVelocity();
+ void BrhoToEnergy();
+ void BrhoToTof() {BrhoToEnergy(); EnergyToTof();}
+ void TofToEnergy();
+ void TofToBrho() {TofToEnergy(); EnergyToBrho();}
+ void EnergyToBeta();
+ void BetaToEnergy();
+ void BetaToGamma();
+ double DopplerCorrection(double EnergyLabGamma, double ThetaLabGamma);
+
+
protected :
- void Extract(const char* line);
+ void Extract(const char* line);
public :
- string GetName() const;
- int GetZ() const {return fCharge;}
- int GetA() const {return fAtomicWeight;}
- double GetMassExcess() const {return fMassExcess;}
- const char* GetSpinParity() const {return fSpinParity;}
- double GetSpin() const {return fSpin;}
- const char* GetParity() const {return fParity;}
- void SetName(const char* name) {fName = name;}
- void SetZ(int charge) {fCharge = charge;}
- void SetA(int mass) {fAtomicWeight = mass;}
- void SetMassExcess(double massexcess) {fMassExcess = massexcess;}
- void SetSpinParity(const char* spinparity) {fSpinParity = spinparity;}
- void SetSpin(double spin) {fSpin = spin;}
- void SetParity(const char* parity) {fParity = parity;}
-
- // Nuclear mass in MeV
+ void GetNucleusName();
+ string GetName() const {return fNucleusName;}
+ int GetZ() const {return fCharge;}
+ int GetA() const {return fAtomicWeight;}
+ double GetMassExcess() const {return fMassExcess;}
+ const char* GetSpinParity() const {return fSpinParity;}
+ double GetSpin() const {return fSpin;}
+ const char* GetParity() const {return fParity;}
+ double GetEnergy() const {return fKineticEnergy;}
+ double GetBrho() const {return fBrho;}
+ double GetTimeOfFlight() const {return fTimeOfFlight;}
+ double GetBeta() const {return fBeta;}
+ double GetGamma() const {return fGamma;}
+ double GetVelocity() const {return fVelocity;}
+ TLorentzVector GetEnergyImpulsion() const {return fEnergyImpulsion;}
+ void SetName(const char* name) {fName = name;}
+ void SetZ(int charge) {fCharge = charge;}
+ void SetA(int mass) {fAtomicWeight = mass;}
+ void SetMassExcess(double massexcess) {fMassExcess = massexcess;}
+ void SetSpinParity(const char* spinparity) {fSpinParity = spinparity;}
+ void SetSpin(double spin) {fSpin = spin;}
+ void SetParity(const char* parity) {fParity = parity;}
+ void SetKineticEnergy(double energy) {fKineticEnergy = energy; EnergyToBrho(); EnergyToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
+ void SetBrho(double brho) {fBrho = brho; BrhoToEnergy(); BrhoToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
+ void SetTimeOfFlight(double tof) {fTimeOfFlight = tof; TofToEnergy(); TofToBrho(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
+ void SetEnergyImpulsion(TLorentzVector P) {fEnergyImpulsion = P;
+ fKineticEnergy = P.E() - Mass();
+ EnergyToBrho();
+ EnergyToTof();
+ EnergyToBeta();
+ BetaToGamma();
+ BetaToVelocity();}
+ void SetBeta(double beta) {fBeta = beta; BetaToGamma(); BetaToEnergy(); EnergyToTof(); EnergyToBrho();BetaToVelocity();}
+
+
+
+
+ // Nuclear mass in MeV
double Mass() const {return (fAtomicWeight*uma + fMassExcess/1000.);}
void Print() const ;
};
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index 65a1c54a8..5f0f99a20 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -48,6 +48,7 @@
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
+
using namespace NPL;
@@ -62,8 +63,8 @@ Reaction::Reaction()
fNuclei4 = new Nucleus();
fBeamEnergy = 0;
fThetaCM = 0;
- fExcitationLight = 0;
- fExcitationHeavy = 0;
+ fExcitation3 = 0;
+ fExcitation4 = 0;
fQValue = 0;
fCrossSectionAngleMin = 0;
fCrossSectionAngleMax = 180;
@@ -91,8 +92,8 @@ void Reaction::SetEveryThing(string name1, string name2, string name3, string na
fNuclei4 = new Nucleus(name4);
fBeamEnergy = BeamEnergy;
fThetaCM = 0;
- fExcitationLight = ExcitationEnergyLight;
- fExcitationHeavy = ExcitationEnergyHeavy;
+ fExcitation3 = ExcitationEnergyLight;
+ fExcitation4 = ExcitationEnergyHeavy;
fQValue = (fNuclei1->GetMassExcess() + fNuclei2->GetMassExcess()
- fNuclei3->GetMassExcess() - fNuclei4->GetMassExcess()) / 1000;
@@ -144,9 +145,9 @@ void Reaction::SetEveryThing(string name1, string name2, string name3, string na
fCrossSectionAngleMax = thetamax;
CSFile.close();
- CrossSectionSize = CrossSectionBuffer.size();
- CrossSection = new double[CrossSectionSize] ;
- for(int i = 0 ; i <CrossSectionSize ; i++ ) CrossSection[i] = CrossSectionBuffer[i];
+ fCrossSectionSize = CrossSectionBuffer.size();
+ fCrossSection = new double[fCrossSectionSize] ;
+ for(int i = 0 ; i < fCrossSectionSize ; i++ ) fCrossSection[i] = CrossSectionBuffer[i];
initializePrecomputeVariable();
}
@@ -161,106 +162,103 @@ Reaction::~Reaction()
delete fNuclei4;
}
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
bool Reaction::CheckKinematic()
{
- // Check if kinematics is allowed
-
- // case of inverse kinematics
- double theta = fThetaCM;
- if (m1 > m2) theta = M_PI - fThetaCM;
-
- // total and kinetic energies in the lab
- double W3lab = W3cm * G * (1 + B*beta3cm*cos(theta));
- double W4lab = W4cm * G * (1 + B*beta4cm*cos(theta + M_PI));
- // test for total energy conversion
- if (fabs(WtotLab - (W3lab+W4lab)) > 1e-6)
- return false;
-
- else return true;
+ double theta = fThetaCM;
+ if (m1 > m2) theta = M_PI - fThetaCM;
+
+ fEnergyImpulsionCM_3 = TLorentzVector(pCM_3*sin(theta),0,pCM_3*cos(theta),ECM_3);
+ fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
+
+ fEnergyImpulsionLab_3 = fEnergyImpulsionCM_3;
+ fEnergyImpulsionLab_3.Boost(0,0,BetaCM);
+ fEnergyImpulsionLab_4 = fEnergyImpulsionCM_4;
+ fEnergyImpulsionLab_4.Boost(0,0,BetaCM);
+
+ if ( fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_3.E()+fEnergyImpulsionLab_4.E()))> 1e-6){
+ cout << "Problem with energy conservation" << endl;
+ return false;
+ }
+
+ else{
+ //cout << "Kinematic OK" << endl;
+ return true;
+ }
+
}
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Reaction::KineRelativistic(double &ThetaLab3, double &EnergieLab3,
- double &ThetaLab4, double &EnergieLab4) const
+void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
+ double &ThetaLab4, double &KineticEnergyLab4)
{
-// 2-body relativistic kinematics: direct + inverse
-// EnergieLab3,4 : lab energy in MeV of the 2 ejectiles
-// ThetaLab3,4 : angles in rad
-
- // case of inverse kinematics
- double theta = fThetaCM;
- if (m1 > m2) theta = M_PI - fThetaCM;
-
- // lab angles
- ThetaLab3 = atan(sin(theta) / (cos(theta) + K3) / G);
- if (ThetaLab3 < 0) ThetaLab3 += M_PI;
- ThetaLab4 = atan(sin(M_PI + theta) / (cos(M_PI + theta) + K4) / G);
- if (fabs(ThetaLab3) < 1e-6) ThetaLab3 = 0;
- ThetaLab4 = fabs(ThetaLab4);
- if (fabs(ThetaLab4) < 1e-6) ThetaLab4 = 0;
-
- // total and kinetic energies in the lab
- double W3lab = W3cm * G * (1 + B*beta3cm*cos(theta));
- double W4lab = W4cm * G * (1 + B*beta4cm*cos(theta + M_PI));
- // test for total energy conversion
- if (fabs(WtotLab - (W3lab+W4lab)) > 1e-6)
- cout << "Problem for energy conservation" << endl;
- EnergieLab3 = W3lab - m3;
- EnergieLab4 = W4lab - m4;
-
+ // 2-body relativistic kinematics: direct + inverse
+ // EnergieLab3,4 : lab energy in MeV of the 2 ejectiles
+ // ThetaLab3,4 : angles in rad
+
+ // case of inverse kinematics
+ double theta = fThetaCM;
+ if (m1 > m2) theta = M_PI - fThetaCM;
+
+ fEnergyImpulsionCM_3 = TLorentzVector(pCM_3*sin(theta),0,pCM_3*cos(theta),ECM_3);
+ fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
+
+ fEnergyImpulsionLab_3 = fEnergyImpulsionCM_3;
+ fEnergyImpulsionLab_3.Boost(0,0,BetaCM);
+ fEnergyImpulsionLab_4 = fEnergyImpulsionCM_4;
+ fEnergyImpulsionLab_4.Boost(0,0,BetaCM);
+
+ // Angle in the lab frame
+ ThetaLab3 = fEnergyImpulsionLab_3.Angle(fEnergyImpulsionLab_1.Vect());
+ if (ThetaLab3 < 0) ThetaLab3 += M_PI;
+
+ ThetaLab4 = fEnergyImpulsionLab_4.Angle(fEnergyImpulsionLab_1.Vect());
+ if (fabs(ThetaLab3) < 1e-6) ThetaLab3 = 0;
+ ThetaLab4 = fabs(ThetaLab4);
+ if (fabs(ThetaLab4) < 1e-6) ThetaLab4 = 0;
+
+ // Kinetic Energy in the lab frame
+ KineticEnergyLab3 = fEnergyImpulsionLab_3.E() - m3;
+ KineticEnergyLab4 = fEnergyImpulsionLab_4.E() - m4;
+
+ // test for total energy conversion
+ if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_3.E()+fEnergyImpulsionLab_4.E())) > 1e-6)
+ cout << "Problem for energy conservation" << endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab) const
+double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab)
{
// EnergyLab in MeV
// ThetaLab in rad
- double P1 = sqrt(2*m1*fBeamEnergy+(fBeamEnergy*fBeamEnergy));
- double P3 = sqrt(2*m3*EnergyLab+(EnergyLab*EnergyLab));
- double P4 = sqrt(P1*P1+P3*P3-(2*P1*P3*cos(ThetaLab)));
- double E4 = fBeamEnergy+m1+m2-(EnergyLab+m3);
- double m4e = sqrt((E4*E4)-(P4*P4));
- double Eex= m4e-fNuclei4->Mass();
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3*E3 - m3*m3);
+
+ fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3*sin(ThetaLab),0,p_Lab_3*cos(ThetaLab),E3);
+ fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
+
+ double Eex = fEnergyImpulsionLab_4.Mag() - fNuclei4->Mass();
- return Eex;
+ return Eex;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//Return ThetaCM
-double Reaction::EnergyLabToThetaCM( double EnergyLab , double ExcitationEnergy ) const
+double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab)
{
- if(ExcitationEnergy == -500) ExcitationEnergy = fExcitationHeavy;
-
- double E1 = (fBeamEnergy+m1) ;
- double E3 = (EnergyLab+m3) ;
-
- // Compute Mandelstan variable
- double s = 2*m2*E1 + m1*m1 + m2*m2 ;
- double u = -2*m2*E3 + m2*m2 + m3*m3 ;
- // Compute CM impulsion:
- //before reaction
- double P2CM = ( sqrt( ( s-(m1-m2)*(m1-m2) )*( s-(m1+m2)*(m1+m2) ) ) ) / (2*sqrt(s)) ;
- // after reaction
- double P3CM = ( sqrt( ( s-(m3-m4)*(m3-m4) )*( s-(m3+m4)*(m3+m4) ) ) ) / (2*sqrt(s)) ;
-
- // Compute CM Energy
- double E2CM = (s + m2*m2 -m1*m1)/(2*sqrt(s)) ;
- double E3CM = (s + m3*m3 -m4*m4)/(2*sqrt(s)) ;
-
- double u0 = m2*m2 + m3*m3 - 2*(E2CM*E3CM + P2CM*P3CM) ;
-
- double Pi = 3.141592654 ;
- double ThetaCM = Pi - acos ( 1-(u-u0)/(2*P2CM*P3CM) ) ;
-
- return(ThetaCM);
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3*E3 - m3*m3);
+
+ fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3*sin(ThetaLab),0,p_Lab_3*cos(ThetaLab),E3);
+ fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3;
+ fEnergyImpulsionCM_3.Boost(0,0,-BetaCM);
+
+ double ThetaCM = CLHEP::pi - fEnergyImpulsionCM_1.Angle(fEnergyImpulsionCM_3.Vect());
+
+ return(ThetaCM);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -273,12 +271,14 @@ void Reaction::Print() const
<< fBeamEnergy << " MeV"
<< endl ;
- cout << "Exc Light = " << fExcitationLight << " MeV" << endl;
- cout << "Exc Heavy = " << fExcitationHeavy << " MeV" << endl;
+ cout << "Exc Nuclei 3 = " << fExcitation3 << " MeV" << endl;
+ cout << "Exc Nuclei 4 = " << fExcitation4 << " MeV" << endl;
cout << "Qgg = " << fQValue << " MeV" << endl;
}
-
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::ReadConfigurationFile(string Path)
{
////////General Reading needs////////
@@ -287,15 +287,15 @@ void Reaction::ReadConfigurationFile(string Path)
////////Reaction Setting needs///////
string Beam, Target, Heavy, Light, CrossSectionPath ;
- double BeamEnergy = 0 , ExcitationEnergyLight = 0, ExcitationEnergyHeavy = 0;
+ double BeamEnergy = 0 , ExcitationEnergy3 = 0, ExcitationEnergy4 = 0;
double CSHalfOpenAngleMin = 0, CSHalfOpenAngleMax = 0;
bool ReadingStatus = false ;
bool check_Beam = false ;
bool check_Target = false ;
bool check_Light = false ;
bool check_Heavy = false ;
- bool check_ExcitationEnergyLight = false ;
- bool check_ExcitationEnergyHeavy = false ;
+ bool check_ExcitationEnergy3 = false ;
+ bool check_ExcitationEnergy4 = false ;
bool check_BeamEnergy = false ;
bool check_CrossSectionPath = false ;
@@ -361,18 +361,18 @@ void Reaction::ReadConfigurationFile(string Path)
cout << "Heavy " << Heavy << endl;
}
- else if (DataBuffer=="ExcitationEnergyLight=") {
- check_ExcitationEnergyLight = true ;
+ else if (DataBuffer=="ExcitationEnergy3=" || DataBuffer=="ExcitationEnergyLight=") {
+ check_ExcitationEnergy3 = true ;
ReactionFile >> DataBuffer;
- ExcitationEnergyLight = atof(DataBuffer.c_str()) * MeV;
- cout << "Excitation Energy Light" << ExcitationEnergyLight / MeV << " MeV" << endl;
+ ExcitationEnergy3 = atof(DataBuffer.c_str()) * MeV;
+ cout << "Excitation Energy Nuclei 3: " << ExcitationEnergy3 / MeV << " MeV" << endl;
}
- else if (DataBuffer=="ExcitationEnergyHeavy=") {
- check_ExcitationEnergyHeavy = true ;
+ else if (DataBuffer=="ExcitationEnergy4=" || DataBuffer=="ExcitationEnergyHeavy=") {
+ check_ExcitationEnergy4 = true ;
ReactionFile >> DataBuffer;
- ExcitationEnergyHeavy = atof(DataBuffer.c_str()) * MeV;
- cout << "Excitation Energy Heavy" << ExcitationEnergyHeavy / MeV << " MeV" << endl;
+ ExcitationEnergy4 = atof(DataBuffer.c_str()) * MeV;
+ cout << "Excitation Energy Nuclei 4: " << ExcitationEnergy4 / MeV << " MeV" << endl;
}
else if (DataBuffer=="BeamEnergy=") {
@@ -383,7 +383,7 @@ void Reaction::ReadConfigurationFile(string Path)
}
else if (DataBuffer== "CrossSectionPath=") {
- check_CrossSectionPath = true ;
+ check_CrossSectionPath = true ;
ReactionFile >> CrossSectionPath;
cout << "Cross Section File: " << CrossSectionPath << endl ;
}
@@ -408,42 +408,203 @@ void Reaction::ReadConfigurationFile(string Path)
///////////////////////////////////////////////////
// If all Token found toggle out
- if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergyLight
- && check_ExcitationEnergyHeavy && check_BeamEnergy && check_CrossSectionPath)
+ if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergy3
+ && check_ExcitationEnergy4 && check_BeamEnergy && check_CrossSectionPath)
ReadingStatus = false;
}
}
- SetEveryThing(Beam, Target, Light, Heavy,BeamEnergy,ExcitationEnergyLight, ExcitationEnergyHeavy,CrossSectionPath, CSHalfOpenAngleMin, CSHalfOpenAngleMax);
+ SetEveryThing(Beam, Target, Light, Heavy,BeamEnergy,ExcitationEnergy3, ExcitationEnergy4,CrossSectionPath, CSHalfOpenAngleMin, CSHalfOpenAngleMax);
ReactionFile.close();
}
-
+
+////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::initializePrecomputeVariable()
{
m1 = fNuclei1->Mass();
m2 = fNuclei2->Mass();
- m3 = fNuclei3->Mass() + fExcitationLight;
- m4 = fNuclei4->Mass() + fExcitationHeavy;
-
- // center-of-mass velocity
- WtotLab = (fBeamEnergy + m1) + m2;
- P1 = sqrt(pow(fBeamEnergy,2) + 2*m1*fBeamEnergy);
- B = P1 / WtotLab;
- G = 1 / sqrt(1 - pow(B,2));
-
- // total energy of the ejectiles in the center-of-mass
- W3cm = (pow(WtotLab,2) + pow(G,2)*(pow(m3,2) - pow(m4,2)))
- / (2 * G * WtotLab);
- W4cm = (pow(WtotLab,2) + pow(G,2)*(pow(m4,2) - pow(m3,2)))
- / (2 * G * WtotLab);
-
- // velocity of the ejectiles in the center-of-mass
- beta3cm = sqrt(1 - pow(m3,2)/pow(W3cm,2));
- beta4cm = sqrt(1 - pow(m4,2)/pow(W4cm,2));
-
- // Constants of the kinematics
- K3 = B / beta3cm;
- K4 = B / beta4cm;
- }
+ m3 = fNuclei3->Mass() + fExcitation3;
+ m4 = fNuclei4->Mass() + fExcitation4;
+
+ s = m1*m1 + m2*m2 + 2*m2*(fBeamEnergy + m1);
+
+ fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
+
+ ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
+ ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
+ ECM_3 = (s + m3*m3 - m4*m4)/(2*sqrt(s));
+ ECM_4 = (s + m4*m4 - m3*m3)/(2*sqrt(s));
+
+ pCM_3 = sqrt(ECM_3*ECM_3 - m3*m3);
+ pCM_4 = sqrt(ECM_4*ECM_4 - m4*m4);
+
+ fImpulsionLab_1 = TVector3(0,0,sqrt(fBeamEnergy*fBeamEnergy + 2*fBeamEnergy*m1));
+ fImpulsionLab_2 = TVector3(0,0,0);
+
+ fEnergyImpulsionLab_1 = TLorentzVector(fImpulsionLab_1,m1+fBeamEnergy);
+ fEnergyImpulsionLab_2 = TLorentzVector(fImpulsionLab_2,m2);
+
+ fTotalEnergyImpulsionLab = fEnergyImpulsionLab_1 + fEnergyImpulsionLab_2;
+
+ BetaCM = fTotalEnergyImpulsionLab.Beta();
+
+ fEnergyImpulsionCM_1 = fEnergyImpulsionLab_1;
+ fEnergyImpulsionCM_1.Boost(0,0,-BetaCM);
+
+ fEnergyImpulsionCM_2 = fEnergyImpulsionLab_2;
+ fEnergyImpulsionCM_2.Boost(0,0,-BetaCM);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::SetNuclei3(double EnergyLab, double ThetaLab)
+{
+ double p3 = sqrt(pow(EnergyLab,2) + 2*m3*EnergyLab);
+
+ fEnergyImpulsionLab_3 = TLorentzVector(p3*sin(ThetaLab),0,p3*cos(ThetaLab),EnergyLab+m3);
+ fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
+
+ fNuclei3->SetEnergyImpulsion(fEnergyImpulsionLab_3);
+ fNuclei4->SetEnergyImpulsion(fEnergyImpulsionLab_4);
+
+ fThetaCM = EnergyLabToThetaCM(EnergyLab, ThetaLab);
+ fExcitation4 = ReconstructRelativistic(EnergyLab, ThetaLab);
+}
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetKinematicLine3()
+{
+ int size = 360;
+ double x[size];
+ double y[size];
+ double theta3,E3,theta4,E4,Brho;
+
+ for (int i = 0; i < size; ++i)
+ {
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei3->SetKineticEnergy(E3);
+ Brho = fNuclei3->GetBrho();
+
+ x[i] = theta3/deg;
+ y[i] = E3;
+ }
+ TGraph* KineLine3 = new TGraph(size,x,y);
+ KineLine3->SetTitle("Kinematic Line of particule 3");
+ return(KineLine3);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetKinematicLine4()
+{
+ int size = 360;
+ double x[size];
+ double y[size];
+ double theta3,E3,theta4,E4,Brho;
+
+ for (int i = 0; i < size; ++i)
+ {
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei4->SetKineticEnergy(E4);
+ Brho = fNuclei4->GetBrho();
+
+ x[i] = theta4/deg;
+ y[i] = E4;
+ }
+ TGraph* KineLine4= new TGraph(size,x,y);
+ KineLine4->SetTitle("Kinematic Line of particule 4");
+ return(KineLine4);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetBrhoLine3()
+{
+ int size = 360;
+ double x[size];
+ double y[size];
+ double theta3,E3,theta4,E4;
+ double Brho;
+
+ for (int i = 0; i < size; ++i)
+ {
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+ fNuclei3->SetKineticEnergy(E3);
+ Brho = fNuclei3->GetBrho();
+
+ x[i] = theta3/deg;
+ y[i] = Brho;
+ }
+ TGraph* LineBrho3= new TGraph(size,x,y);
+ return(LineBrho3);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetThetaLabVersusThetaCM()
+{
+ int size = 360;
+ double x[size];
+ double y[size];
+ double theta3,E3,theta4,E4;
+
+ for (int i = 0; i < size; ++i)
+ {
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+
+ x[i] = fThetaCM/deg;
+ y[i] = theta3/deg;
+ }
+ TGraph* AngleLine= new TGraph(size,x,y);
+ return(AngleLine);
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::PrintKinematic()
+{
+ int size = 360;
+ double theta3,E3,theta4,E4,Brho3,Brho4;
+
+ cout << endl;
+ cout << "*********************** Print Kinematic ***********************" << endl;
+ cout << "ThetaCM" << " " << "ThetaLab" << " " << "EnergyLab3" << " " << "Brho3" << " " << "EnergyLab4" << " " << "Brho4" << endl;
+ for (int i = 0; i < size; ++i)
+ {
+ SetThetaCM(((double)i)/2*deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+
+ fNuclei3->SetKineticEnergy(E3);
+ Brho3 = fNuclei3->GetBrho();
+
+ fNuclei4->SetKineticEnergy(E4);
+ Brho4 = fNuclei4->GetBrho();
+
+ cout << (double)i/2 << " " << theta3/deg << " " << E3 << " " << Brho3 << " " << E4 << " " << Brho4 << endl;
+ }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 4937e14f5..1e80b060a 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -39,8 +39,15 @@
#include "NPNucleus.h"
+// ROOT header
+#include "TLorentzVector.h"
+#include "TLorentzRotation.h"
+#include "TVector3.h"
+#include "TGraph.h"
+#include "TCanvas.h"
+#include "TH2F.h"
+
using namespace std;
-//class Nucleus;
namespace NPL
{
@@ -57,85 +64,115 @@ namespace NPL
void ReadConfigurationFile(string Path);
private:
- Nucleus *fNuclei1; // Beam
- Nucleus *fNuclei2; // Target
- Nucleus *fNuclei3; // Light ejectile
- Nucleus *fNuclei4; // Heavy ejectile
- double fQValue; // Q-value in MeV
- double fBeamEnergy; // Beam energy in MeV
- double fThetaCM; // Center-of-mass angle in radian
- double fExcitation; // Excitation energy in MeV
- double fExcitationLight; // Excitation energy in MeV
- double fExcitationHeavy; // Excitation energy in MeV
- double* CrossSection; // Differential CrossSection
- int CrossSectionSize; // Size of array containing Differention CrossSection
- double fCrossSectionAngleMin; // Minimum angle of the differential cross-section given by the user
- double fCrossSectionAngleMax; // Maximum angle of the differential cross-section given by the user
+ Nucleus *fNuclei1; // Beam
+ Nucleus *fNuclei2; // Target
+ Nucleus *fNuclei3; // Light ejectile
+ Nucleus *fNuclei4; // Heavy ejectile
+ double fQValue; // Q-value in MeV
+ double fBeamEnergy; // Beam energy in MeV
+ double fThetaCM; // Center-of-mass angle in radian
+ double fExcitation3; // Excitation energy in MeV
+ double fExcitation4; // Excitation energy in MeV
+ double* fCrossSection; // Differential CrossSection
+ int fCrossSectionSize; // Size of array containing Differention CrossSection
+ double fCrossSectionAngleMin; // Minimum angle of the differential cross-section given by the user
+ double fCrossSectionAngleMax; // Maximum angle of the differential cross-section given by the user
public:
// Getters and Setters
void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();}
- void SetExcitationLight(double exci) {fExcitationLight = exci; initializePrecomputeVariable();}
- void SetExcitationHeavy(double exci) {fExcitationHeavy = exci; initializePrecomputeVariable();}
+ void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
+ void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;}
- double GetExcitationHeavy() const {return fExcitationHeavy;}
- double GetExcitation() const {return fExcitationHeavy;}
- double GetExcitationLight() const {return fExcitationLight;}
+ double GetExcitation3() const {return fExcitation3;}
+ double GetExcitation4() const {return fExcitation4;}
double GetQValue() const {return fQValue;}
Nucleus* GetNucleus1() const {return fNuclei1;}
Nucleus* GetNucleus2() const {return fNuclei2;}
Nucleus* GetNucleus3() const {return fNuclei3;}
Nucleus* GetNucleus4() const {return fNuclei4;}
- double* GetCrossSection() const {return CrossSection;}
- int GetCrossSectionSize() const {return CrossSectionSize;}
+ double* GetCrossSection() const {return fCrossSection;}
+ int GetCrossSectionSize() const {return fCrossSectionSize;}
double GetCrossSectionAngleMin() const {return fCrossSectionAngleMin;}
double GetCrossSectionAngleMax() const {return fCrossSectionAngleMax;}
- private: // intern precompute variable
- void initializePrecomputeVariable();
- double m1;
- double m2;
- double m3;
- double m4;
-
- // center-of-mass velocity
- double WtotLab;
- double P1;
- double B;
- double G;
-
- // total energy of the ejectiles in the center-of-mass
- double W3cm;
- double W4cm;
-
- // velocity of the ejectiles in the center-of-mass
- double beta3cm;
- double beta4cm;
-
- // Constants of the kinematics
- double K3;
- double K4;
+
+ private: // intern precompute variable
+ void initializePrecomputeVariable();
+ double m1;
+ double m2;
+ double m3;
+ double m4;
+
+ // Lorents Vector
+ TLorentzVector fEnergyImpulsionLab_1;
+ TLorentzVector fEnergyImpulsionLab_2;
+ TLorentzVector fEnergyImpulsionLab_3;
+ TLorentzVector fEnergyImpulsionLab_4;
+ TLorentzVector fTotalEnergyImpulsionLab;
+
+ TLorentzVector fEnergyImpulsionCM_1;
+ TLorentzVector fEnergyImpulsionCM_2;
+ TLorentzVector fEnergyImpulsionCM_3;
+ TLorentzVector fEnergyImpulsionCM_4;
+ TLorentzVector fTotalEnergyImpulsionCM;
+
+ // Impulsion Vector3
+ TVector3 fImpulsionLab_1;
+ TVector3 fImpulsionLab_2;
+ TVector3 fImpulsionLab_3;
+ TVector3 fImpulsionLab_4;
+
+ TVector3 fImpulsionCM_1;
+ TVector3 fImpulsionCM_2;
+ TVector3 fImpulsionCM_3;
+ TVector3 fImpulsionCM_4;
+
+ // CM Energy composante & CM impulsion norme
+ Double_t ECM_1;
+ Double_t ECM_2;
+ Double_t ECM_3;
+ Double_t ECM_4;
+ Double_t pCM_3;
+ Double_t pCM_4;
+
+ // Mandelstam variable
+ Double_t s;
+
+ // Center of Mass Kinematic
+ Double_t BetaCM;
+
public: // Kinematics
- // Check that the reaction is alowed
- bool CheckKinematic();
-
- // Compute ThetaLab and EnergyLab for product of reaction
- void KineRelativistic(double &ThetaLab3, double &EnergieLab3,
- double &ThetaLab4, double &EnergieLab4) const;
-
- // Return Excitation Energy
- double ReconstructRelativistic(double EnergyLab, double ThetaLab) const;
-
- // Return ThetaCM
- // EnergyLab: energy measured in the laboratory frame
- // ExcitationEnergy: excitation energy previously calculated. If no argument given, fExcitation is used
- double EnergyLabToThetaCM(double EnergyLab, double ExcitationEnergy = -500) const;
-
- // Print private paremeter
- void Print() const;
+ // Check that the reaction is alowed
+ bool CheckKinematic();
+
+ // Compute ThetaLab and EnergyLab for product of reaction
+ void KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
+ double &ThetaLab4, double &KineticEnergyLab4);
+
+ // Return Excitation Energy
+ double ReconstructRelativistic(double EnergyLab, double ThetaLab);
+
+
+ // Return ThetaCM
+ // EnergyLab: energy measured in the laboratory frame
+ // ExcitationEnergy: excitation energy previously calculated.
+ double EnergyLabToThetaCM(double EnergyLab, double ThetaLab);
+
+ void SetNuclei3(double EnergyLab, double ThetaLab);
+
+
+ TGraph* GetKinematicLine3();
+ TGraph* GetKinematicLine4();
+ TGraph* GetBrhoLine3();
+ TGraph* GetThetaLabVersusThetaCM();
+ void PrintKinematic();
+
+ // Print private paremeter
+ void Print() const;
};
}
#endif
diff --git a/NPLib/Vamos/TVamosCHIODataDict.cxx b/NPLib/Vamos/TVamosCHIODataDict.cxx
index df0790581..f99e9be06 100644
--- a/NPLib/Vamos/TVamosCHIODataDict.cxx
+++ b/NPLib/Vamos/TVamosCHIODataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:44:03 2012
+// File generated by rootcint at Wed Jul 18 14:42:27 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPLib/Vamos/TVamosDCDataDict.cxx b/NPLib/Vamos/TVamosDCDataDict.cxx
index b3377069a..350ba7911 100644
--- a/NPLib/Vamos/TVamosDCDataDict.cxx
+++ b/NPLib/Vamos/TVamosDCDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:44:04 2012
+// File generated by rootcint at Wed Jul 18 14:42:28 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPLib/Vamos/TVamosFingerDataDict.cxx b/NPLib/Vamos/TVamosFingerDataDict.cxx
index 1990684d7..8edb1bab5 100644
--- a/NPLib/Vamos/TVamosFingerDataDict.cxx
+++ b/NPLib/Vamos/TVamosFingerDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:44:02 2012
+// File generated by rootcint at Wed Jul 18 14:42:27 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPLib/Vamos/TVamosPlasticDataDict.cxx b/NPLib/Vamos/TVamosPlasticDataDict.cxx
index 76df2db6f..1de81741e 100644
--- a/NPLib/Vamos/TVamosPlasticDataDict.cxx
+++ b/NPLib/Vamos/TVamosPlasticDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 11:44:01 2012
+// File generated by rootcint at Wed Jul 18 14:42:26 2012
// Do NOT change. Changes will be lost next time file is generated
//
diff --git a/NPSimulation/include/EventGeneratorTransfertToResonance.hh b/NPSimulation/include/EventGeneratorTransfertToResonance.hh
deleted file mode 100644
index 8b2b5e941..000000000
--- a/NPSimulation/include/EventGeneratorTransfertToResonance.hh
+++ /dev/null
@@ -1,151 +0,0 @@
-#ifndef EventGeneratorTransfertToResonance_H
-#define EventGeneratorTransfertToResonance_H
-/*****************************************************************************
- * Copyright (C) 2009-2010 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
- * *
- * Creation Date : January 2009 *
- * Last update : January 2011 *
- *---------------------------------------------------------------------------*
- * 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: *
- * + 20/01/2011: Add support for excitation energy for light ejectile *
- * (N. de Sereville) *
- * *
- * *
- *****************************************************************************/
-// C++ header
-#include <string>
-
-// NPSimulation
-#include "VEventGenerator.hh"
-#include "Target.hh"
-
-// NPLib
-#include "TInitialConditions.h"
-
-// NPLib header
-#include "NPReaction.h"
-
-using namespace std;
-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 ExcitationEnergyLight , //Excitation of Light Nuclei
- double ExcitationEnergyHeavy , //Excitation of Heavy Nuclei
- double BeamEnergySpread ,
- double SigmaX ,
- double SigmaY ,
- double SigmaThetaX ,
- double SigmaPhiY ,
- double ResonanceWidth ,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path ,
- double CSThetaMin ,
- double CSThetaMax); //Path of the differentiel Cross Section
-
- // Default Destructor
- virtual ~EventGeneratorTransfertToResonance();
-
- public: // Inherit from VEventGenerator class
- void ReadConfiguration(string);
- void GenerateEvent(G4Event*, G4ParticleGun*);
- void SetTarget(Target* Target);
-
- 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
- Target* m_Target;
-
- private: // Reaction
- Reaction* m_Reaction ;
-
- private: // Resonance decay channel
- double m_ResonanceWidth ;
- double m_ResonanceMean ;
- int m_ResonanceDecayZ ;
- int m_ResonanceDecayA ;
-
- // When the Phase Space Generator is called, the weight of the current configuration is return and stored in this variable
- // Spectrum then need to be weighted by this paramater to be realistic
- // NB: This procedure avoid long calculation time of the rejection methods previously instantiate and therefore allow simulation of manybody phase space decay
- double m_EventWeight ;
-
- //Other
- void InitializeRootOutput() ;
- void ResonanceDecay( 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 ExcitationEnergyLight , //Excitation of Light Nuclei
- double ExcitationEnergyHeavy , //Excitation of Heavy Nuclei
- double BeamEnergySpread ,
- double SigmaX ,
- double SigmaY ,
- double SigmaThetaX ,
- double SigmaPhiY ,
- double ResonanceWidth ,
- int ResonanceDecayZ ,
- int ResonanceDecayA ,
- bool ShootLight ,
- bool ShootHeavy ,
- bool ShootDecayProduct ,
- string Path ,
- double CSThetaMin ,
- double CSThetaMax); //Path of the differentiel Cross Section
-};
-
-#endif
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 435b5f687..90c42c4b9 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -199,109 +199,109 @@ void EventGeneratorTransfert::ReadConfiguration(string Path)
//Search for comment Symbol %
if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
- else if (DataBuffer.compare(0, 5, "Beam=") == 0) {
+ else if (DataBuffer=="Beam=") {
check_Beam = true ;
ReactionFile >> DataBuffer;
Beam = DataBuffer;
G4cout << "Beam " << Beam << G4endl;
}
- else if (DataBuffer.compare(0, 7, "Target=") == 0) {
+ else if (DataBuffer=="Target=") {
check_Target = true ;
ReactionFile >> DataBuffer;
Target = DataBuffer;
G4cout << "Target " << Target << G4endl;
}
- else if (DataBuffer.compare(0, 6, "Light=") == 0) {
+ else if (DataBuffer=="Light=") {
check_Light = true ;
ReactionFile >> DataBuffer;
Light = DataBuffer;
G4cout << "Light " << Light << G4endl;
}
- else if (DataBuffer.compare(0, 6, "Heavy=") == 0) {
+ else if (DataBuffer=="Heavy=") {
check_Heavy = true ;
ReactionFile >> DataBuffer;
Heavy = DataBuffer;
G4cout << "Heavy " << Heavy << G4endl;
}
- else if (DataBuffer.compare(0, 22, "ExcitationEnergyLight=") == 0) {
+ else if (DataBuffer=="ExcitationEnergy3=" || DataBuffer=="ExcitationEnergyLight=") {
check_ExcitationEnergyLight = true ;
ReactionFile >> DataBuffer;
ExcitationEnergyLight = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Excitation Energy Light" << ExcitationEnergyLight / MeV << " MeV" << G4endl;
+ G4cout << "Excitation Energy Nuclei 3: " << ExcitationEnergyLight / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 22, "ExcitationEnergyHeavy=") == 0) {
+ else if (DataBuffer=="ExcitationEnergy4=" || DataBuffer=="ExcitationEnergyHeavy=") {
check_ExcitationEnergyHeavy = true ;
ReactionFile >> DataBuffer;
ExcitationEnergyHeavy = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Excitation Energy Heavy" << ExcitationEnergyHeavy / MeV << " MeV" << G4endl;
+ G4cout << "Excitation Energy Nuclei 4: " << ExcitationEnergyHeavy / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "BeamEnergy=") == 0) {
+ else if (DataBuffer=="BeamEnergy=") {
check_BeamEnergy = true ;
ReactionFile >> DataBuffer;
BeamEnergy = atof(DataBuffer.c_str()) * MeV;
G4cout << "Beam Energy " << BeamEnergy / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 17, "BeamEnergySpread=") == 0) {
+ else if (DataBuffer=="BeamEnergySpread=") {
check_BeamEnergySpread = true ;
ReactionFile >> DataBuffer;
BeamEnergySpread = atof(DataBuffer.c_str()) * MeV;
G4cout << "Beam Energy Spread " << BeamEnergySpread / MeV << " MeV" << G4endl;
}
- else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
+ else if (DataBuffer=="SigmaX=") {
check_FWHMX = true ;
ReactionFile >> DataBuffer;
SigmaX = atof(DataBuffer.c_str()) * mm;
G4cout << "Beam FWHM X " << SigmaX << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
+ else if (DataBuffer=="SigmaY=") {
check_FWHMY = true ;
ReactionFile >> DataBuffer;
SigmaY = atof(DataBuffer.c_str()) * mm;
G4cout << "Beam FWHM Y " << SigmaX << " mm" << G4endl;
}
- else if (DataBuffer.compare(0, 12, "SigmaThetaX=") == 0) {
+ else if (DataBuffer=="SigmaThetaX=") {
check_EmmitanceTheta = true ;
ReactionFile >> DataBuffer;
SigmaThetaX = atof(DataBuffer.c_str()) * deg;
G4cout << "Beam Emmitance Theta " << SigmaThetaX / deg << " deg" << G4endl;
}
- else if (DataBuffer.compare(0, 10, "SigmaPhiY=") == 0) {
+ else if (DataBuffer=="SigmaPhiY=") {
check_EmmitancePhi = true ;
ReactionFile >> DataBuffer;
SigmaPhiY = atof(DataBuffer.c_str()) * deg;
G4cout << "Beam Emmitance Phi " << SigmaPhiY / deg << " deg" << G4endl;
}
- else if (DataBuffer.compare(0, 17, "CrossSectionPath=") == 0) {
+ else if (DataBuffer=="CrossSectionPath=") {
check_CrossSectionPath = true ;
ReactionFile >> CrossSectionPath;
G4cout << "Cross Section File: " << CrossSectionPath << G4endl ;
}
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
+ else if (DataBuffer=="HalfOpenAngleMin=") {
ReactionFile >> DataBuffer;
CSHalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
G4cout << "HalfOpenAngleMin " << CSHalfOpenAngleMin / deg << " degree" << G4endl;
}
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMax=") == 0) {
+ else if (DataBuffer=="HalfOpenAngleMax=") {
ReactionFile >> DataBuffer;
CSHalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
G4cout << "HalfOpenAngleMax " << CSHalfOpenAngleMax / deg << " degree" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "ShootLight=") == 0) {
+ else if (DataBuffer=="ShootLight=") {
check_ShootLight = true ;
ReactionFile >> DataBuffer;
if (atof(DataBuffer.c_str()) == 1) ShootLight = true ;
@@ -309,7 +309,7 @@ void EventGeneratorTransfert::ReadConfiguration(string Path)
else G4cout << "Shoot Light particle : no" << G4endl;
}
- else if (DataBuffer.compare(0, 11, "ShootHeavy=") == 0) {
+ else if (DataBuffer=="ShootHeavy=") {
check_ShootHeavy = true ;
ReactionFile >> DataBuffer;
if (atof(DataBuffer.c_str()) == 1) ShootHeavy = true ;
@@ -383,14 +383,14 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent)
G4int LightA = m_Reaction->GetNucleus3()->GetA() ;
G4ParticleDefinition* LightName
- = G4ParticleTable::GetParticleTable()->GetIon(LightZ, LightA, m_Reaction->GetExcitationLight()*MeV);
+ = G4ParticleTable::GetParticleTable()->GetIon(LightZ, LightA, m_Reaction->GetExcitation3()*MeV);
// Recoil
G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
G4ParticleDefinition* HeavyName
- = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitationHeavy()*MeV);
+ = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation4()*MeV);
// Beam
G4int BeamZ = m_Reaction->GetNucleus1()->GetZ();
diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc
deleted file mode 100644
index dd4c1a3d5..000000000
--- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc
+++ /dev/null
@@ -1,823 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-2010 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
- * *
- * Creation Date : January 2009 *
- * Last update : January 2011 *
- *---------------------------------------------------------------------------*
- * 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: *
- * + 20/01/2011: Add support for excitation energy for light ejectile *
- * (N. de Sereville) *
- * *
- * *
- *****************************************************************************/
-
-// C++ headers
-#include <iostream>
-#include <fstream>
-#include <limits>
-
-// G4 header defining G4 types
-#include "globals.hh"
-
-// G4 headers
-#include "G4ParticleTable.hh"
-#include "G4ParticleGun.hh"
-#include "G4RotationMatrix.hh"
-
-// G4 headers including CLHEP headers
-// for generating random numbers
-#include "Randomize.hh"
-
-// NPTool headers
-#include "EventGeneratorTransfertToResonance.hh"
-#include "RootOutput.h"
-
-//Root Headers
-#include "TGenPhaseSpace.h"
-
-//CLHEP
-#include "CLHEP/Random/RandBreitWigner.h"
-
-using namespace std;
-using namespace CLHEP;
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance()
-{
- //------------- Default Constructor -------------
- m_InitConditions = new TInitialConditions();
- m_Reaction = new Reaction();
- m_Target = new Target();
- m_SigmaX = 0;
- m_SigmaY = 0;
- m_SigmaThetaX = 0;
- m_SigmaPhiY = 0;
- m_ResonanceDecayZ = 0;
- m_ResonanceDecayA = 0;
- m_EventWeight = 0;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorTransfertToResonance::~EventGeneratorTransfertToResonance()
-{
- //------------- Default Destructor ------------
- delete m_InitConditions;
- delete m_Reaction;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::SetTarget(Target* Target)
- {
- if(Target!=0)
- m_Target = Target;
- }
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorTransfertToResonance::EventGeneratorTransfertToResonance( string name1,
- string name2,
- string name3,
- string name4,
- double BeamEnergy,
- double ExcitationEnergyLight,
- double ExcitationEnergyHeavy,
- double BeamEnergySpread,
- double SigmaX,
- double SigmaY,
- double SigmaThetaX,
- double SigmaPhiY,
- double ResonanceWidth,
- int ResonanceDecayZ,
- int ResonanceDecayA,
- bool ShootLight,
- bool ShootHeavy,
- bool ShootDecayProduct,
- string Path,
- double CSThetaMin,
- double CSThetaMax)
-{
- //------------- Constructor with nuclei names and beam energy ------------
-
- SetEverything( name1,
- name2,
- name3,
- name4,
- BeamEnergy,
- ExcitationEnergyLight,
- ExcitationEnergyHeavy,
- BeamEnergySpread,
- SigmaX,
- SigmaY,
- SigmaThetaX,
- SigmaPhiY,
- ResonanceWidth,
- ResonanceDecayZ,
- ResonanceDecayA,
- ShootLight,
- ShootHeavy,
- ShootDecayProduct,
- Path,
- CSThetaMin,
- CSThetaMax);
-
- m_EventWeight = 0;
-
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::InitializeRootOutput()
-{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
- pTree->Branch("EventWeight",&m_EventWeight,"EventWeigt/D");
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Inherit from VEventGenerator
-
-void EventGeneratorTransfertToResonance::ReadConfiguration(string Path)
-{
-////////General Reading needs////////
- string LineBuffer;
- string DataBuffer;
-
-////////Reaction Setting needs///////
- string Beam, Target, Heavy, Light, CrossSectionPath ;
- G4double BeamEnergy = 0, ExcitationEnergyLight = 0, ExcitationEnergyHeavy = 0;
- G4double BeamEnergySpread = 0 , SigmaX = 0 , SigmaY = 0 , SigmaThetaX = 0 , SigmaPhiY=0, ResonanceWidth = 0 ,ResonanceDecayZ = 0 , ResonanceDecayA = 0 ;
- G4double CSHalfOpenAngleMin = 0; G4double CSHalfOpenAngleMax = 180;
- bool ShootLight = false ;
- bool ShootHeavy = false ;
- bool ShootDecayProduct = false ;
-
- bool ReadingStatus = false ;
- bool check_Beam = false ;
- bool check_Target = false ;
- bool check_Light = false ;
- bool check_Heavy = false ;
- bool check_ExcitationEnergyLight = false ;
- bool check_ExcitationEnergyHeavy = false ;
- bool check_BeamEnergy = false ;
- bool check_BeamEnergySpread = false ;
- bool check_FWHMX = false ;
- bool check_FWHMY = false ;
- bool check_EmmitanceTheta = false ;
- bool check_EmmitancePhi = false ;
- bool check_CrossSectionPath = false ;
- bool check_ShootLight = false ;
- bool check_ShootHeavy = false ;
- bool check_ResonanceWidth = false ;
- bool check_ResonanceDecayZ = false ;
- bool check_ResonanceDecayA = false ;
- bool check_ShootDecayProduct = false ;
-
-//////////////////////////////////////////////////////////////////////////////////////////
- ifstream ReactionFile;
- ReactionFile.open(Path.c_str());
-
- if (ReactionFile.is_open()) {} else {
- return;
- }
-
- while (!ReactionFile.eof()) {
- //Pick-up next line
- getline(ReactionFile, LineBuffer);
-
-
-
- if (LineBuffer.compare(0, 20, "TransfertToResonance") == 0) { ReadingStatus = true ;}
-
-
- while(ReadingStatus){
-
- ReactionFile >> DataBuffer;
-
- //Search for comment Symbol %
- if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
-
- else if (DataBuffer.compare(0, 5, "Beam=") == 0) {
- check_Beam = true ;
- ReactionFile >> DataBuffer;
- Beam = DataBuffer;
- G4cout << "Beam " << Beam << G4endl;
- }
-
- else if (DataBuffer.compare(0, 7, "Target=") == 0) {
- check_Target = true ;
- ReactionFile >> DataBuffer;
- Target = DataBuffer;
- G4cout << "Target " << Target << G4endl;
- }
-
- else if (DataBuffer.compare(0, 6, "Light=") == 0) {
- check_Light = true ;
- ReactionFile >> DataBuffer;
- Light = DataBuffer;
- G4cout << "Light " << Light << G4endl;
- }
-
- else if (DataBuffer.compare(0, 6, "Heavy=") == 0) {
- check_Heavy = true ;
- ReactionFile >> DataBuffer;
- Heavy = DataBuffer;
- G4cout << "Heavy " << Heavy << G4endl;
- }
-
- else if (DataBuffer.compare(0, 22, "ExcitationEnergyLight=") == 0) {
- check_ExcitationEnergyLight = true ;
- ReactionFile >> DataBuffer;
- ExcitationEnergyLight = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Excitation Energy Light " << ExcitationEnergyLight / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 22, "ExcitationEnergyHeavy=") == 0) {
- check_ExcitationEnergyHeavy = true ;
- ReactionFile >> DataBuffer;
- ExcitationEnergyHeavy = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Excitation Energy Heavy " << ExcitationEnergyHeavy / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 11, "BeamEnergy=") == 0) {
- check_BeamEnergy = true ;
- ReactionFile >> DataBuffer;
- BeamEnergy = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Beam Energy " << BeamEnergy / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "BeamEnergySpread=") == 0) {
- check_BeamEnergySpread = true ;
- ReactionFile >> DataBuffer;
- BeamEnergySpread = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Beam Energy Spread " << BeamEnergySpread / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
- check_FWHMX = true ;
- ReactionFile >> DataBuffer;
- SigmaX = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM X " << SigmaX << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
- check_FWHMY = true ;
- ReactionFile >> DataBuffer;
- SigmaY = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM Y " << SigmaX << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 19, "SigmaThetaX=") == 0) {
- check_EmmitanceTheta = true ;
- ReactionFile >> DataBuffer;
- SigmaThetaX = atof(DataBuffer.c_str()) * deg;
- G4cout << "Beam Emmitance Theta " << SigmaThetaX / deg << " deg" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "SigmaPhiY=") == 0) {
- check_EmmitancePhi = true ;
- ReactionFile >> DataBuffer;
- SigmaPhiY = atof(DataBuffer.c_str()) * deg;
- G4cout << "Beam Emmitance Phi " << SigmaPhiY / deg << " deg" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 15, "ResonanceWidth=") == 0) {
- check_ResonanceWidth = true ;
- ReactionFile >> DataBuffer;
- ResonanceWidth = atof(DataBuffer.c_str())*MeV;
- G4cout << "Resonance Width " << ResonanceWidth/MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "ResonanceDecayZ=") == 0) {
- check_ResonanceDecayZ = true ;
- ReactionFile >> DataBuffer;
- ResonanceDecayZ = atof(DataBuffer.c_str());
- G4cout << "ResonanceDecayZ " << ResonanceDecayZ << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "ResonanceDecayA=") == 0) {
- check_ResonanceDecayA = true ;
- ReactionFile >> DataBuffer;
- ResonanceDecayA = atof(DataBuffer.c_str());
- G4cout << "ResonanceDecayA " << ResonanceDecayA << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "CrossSectionPath=") == 0) {
- check_CrossSectionPath = true ;
- ReactionFile >> CrossSectionPath;
- G4cout << "Cross Section File: " << CrossSectionPath << G4endl ;
- }
-
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
- ReactionFile >> DataBuffer;
- CSHalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
- G4cout << "HalfOpenAngleMin " << CSHalfOpenAngleMin / deg << " degree" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMax=") == 0) {
- ReactionFile >> DataBuffer;
- CSHalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
- G4cout << "HalfOpenAngleMax " << CSHalfOpenAngleMax / deg << " degree" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 11, "ShootLight=") == 0) {
- check_ShootLight = true ;
- ReactionFile >> DataBuffer;
- if (atof(DataBuffer.c_str()) == 1) ShootLight = true ;
- if (ShootLight) G4cout << "Shoot Light particle : yes" << G4endl;
- else G4cout << "Shoot Light particle : no" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 11, "ShootHeavy=") == 0) {
- check_ShootHeavy = true ;
- ReactionFile >> DataBuffer;
- if (atof(DataBuffer.c_str()) == 1) ShootHeavy = true ;
- if (ShootHeavy) G4cout << "Shoot Heavy particle : yes" << G4endl;
- else G4cout << "Shoot Heavy particle : no" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 18, "ShootDecayProduct=") == 0) {
- check_ShootDecayProduct = true ;
- ReactionFile >> DataBuffer;
- if (atof(DataBuffer.c_str()) == 1) ShootDecayProduct = true ;
- if (ShootDecayProduct) G4cout << "Shoot DecayProducts from decay : yes" << G4endl;
- else G4cout << "Shoot DecayProducts from decay : no" << G4endl;
- }
-
-
- ///////////////////////////////////////////////////
- // If no Transfert Token and no comment, toggle out
- else
- {ReadingStatus = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
-
- ///////////////////////////////////////////////////
- // If all Token found toggle out
- if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergyLight && check_ExcitationEnergyHeavy
- && check_BeamEnergy && check_BeamEnergySpread && check_FWHMX && check_FWHMY && check_EmmitanceTheta
- && check_EmmitancePhi && check_CrossSectionPath && check_ShootLight && check_ShootHeavy
- && check_ResonanceWidth && check_ResonanceDecayZ && check_ResonanceDecayA && check_ShootDecayProduct)
- ReadingStatus = false ;
-
- }
-
-
- }
-
- SetEverything( Beam,
- Target,
- Light,
- Heavy,
- BeamEnergy,
- ExcitationEnergyLight,
- ExcitationEnergyHeavy,
- BeamEnergySpread,
- SigmaX,
- SigmaY,
- SigmaThetaX,
- SigmaPhiY,
- ResonanceWidth,
- ResonanceDecayZ,
- ResonanceDecayA,
- ShootLight,
- ShootHeavy,
- ShootDecayProduct,
- CrossSectionPath,
- CSHalfOpenAngleMin,
- CSHalfOpenAngleMax);
-
- ReactionFile.close();
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
-{
- // Initialize event weight to one.
- m_EventWeight = 1 ;
-
- // If first time, write the DeDx table
- if(anEvent->GetEventID()==0)
- {
- //-------------- Before living, wrtie the DeDx Table -------------------
-
- G4int LightZx = m_Reaction->GetNucleus3()->GetZ() ;
- G4int LightAx = m_Reaction->GetNucleus3()->GetA() ;
-
- G4int BeamZx = m_Reaction->GetNucleus1()->GetZ() ;
- G4int BeamAx = m_Reaction->GetNucleus1()->GetA() ;
-
- if(m_Target!=0)
- {
- m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZx,LightAx, 0.) ,0, m_BeamEnergy+4*m_BeamEnergySpread);
- m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(BeamZx,BeamAx, 0.) ,0, m_BeamEnergy+4*m_BeamEnergySpread);
- }
-
- }
-
- // Clear contents of Precedent event (now stored in ROOTOutput)
- m_InitConditions->Clear();
-
- //////////////////////////////////////////////////
- //////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.);
-
- // Shoot the Resonance energy following the mean and width value
- // EXX should always be more than specific heat of the reaction
- // double EXX = RandBreitWigner::shoot(m_ResonanceMean,m_ResonanceWidth) ;
- double EXX = RandGauss::shoot(m_ResonanceMean,m_ResonanceWidth) ;
- m_Reaction->SetExcitationHeavy( EXX );
-
- while ( m_Reaction->CheckKinematic()==false )
- {
-// EXX = RandBreitWigner::shoot(m_ResonanceMean,m_ResonanceWidth) ;
- EXX = RandGauss::shoot(m_ResonanceMean,m_ResonanceWidth) ;
- m_Reaction->SetExcitationHeavy( EXX );
- }
- // 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 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 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();
- 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);
-
- // 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 //////
- /////////////////////////////////////////////////////////////////
- // Angles
- RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
- G4double ThetaCM = (m_Reaction->GetCrossSectionAngleMin() + CrossSectionShoot.shoot() * (m_Reaction->GetCrossSectionAngleMax() - m_Reaction->GetCrossSectionAngleMin())) * deg;
- G4double phi = RandFlat::shoot() * 2*pi;
- // write angles to ROOT file
- m_InitConditions->SetICEmittedAngleThetaCM(ThetaCM / deg);
- m_InitConditions->SetICEmittedAnglePhiIncidentFrame(phi / deg);
-
- //////////////////////////////////////////////////
- ///// Momentum and angles from kinematics /////
- ///// Angles are in the beam frame /////
- //////////////////////////////////////////////////
- // 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+pi),
- sin(ThetaHeavy) * sin(phi+pi),
- cos(ThetaHeavy));
-
- //////////////////////////////////////////////////
- ///////// Set up everything for shooting /////////
- //////////////////////////////////////////////////
- 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);
- // write angles/energy to ROOT file
- m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(ThetaLight / deg);
- m_InitConditions->SetICEmittedEnergy(EnergyLight/MeV);
- //Set the gun to shoot
- particleGun->SetParticleMomentumDirection(momentum_kine_world);
- //Shoot the light particle
- particleGun->GeneratePrimaryVertex(anEvent);
- }
-
- // Case of recoil particle
- // 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;
-
- if(m_ShootHeavy || m_ShootDecayProduct)
- EventGeneratorTransfertToResonance::ResonanceDecay( EnergyHeavy,
- theta_world,
- phi_world,
- x0,
- y0,
- z0,
- anEvent,
- particleGun);
-
-
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy,
- G4double ThetaHeavy,
- G4double PhiHeavy,
- G4double x0,
- G4double y0,
- G4double z0,
- G4Event* anEvent,
- G4ParticleGun* particleGun)
-{
- 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 ;
-
- 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, 0 );
-
- 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() + m_Reaction->GetExcitation()*MeV;
- G4double md = daughter -> GetPDGMass() ;
- G4double mn = neutron -> GetPDGMass() ;
- G4double mp = proton -> GetPDGMass() ;
-
- // Check that we are above threshold:
- // If the Resonnance go below the threshold, decay is forced at thereshold
- if (M < md + NumberOfNeutrons*mn + NumberOfProtons*mp)
- {
- double NewExx = parent -> GetPDGMass() - md - NumberOfNeutrons*mn - NumberOfProtons*mp ;
- M = parent -> GetPDGMass() + NewExx ;
- }
-
-
- G4double InitialE = EnergyHeavy + M ;
- G4double InitialMomentumX = sqrt( InitialE*InitialE - M*M) * sin(ThetaHeavy) * cos(PhiHeavy) ;
- G4double InitialMomentumY = sqrt( InitialE*InitialE - M*M) * sin(ThetaHeavy) * sin(PhiHeavy) ;
- G4double InitialMomentumZ = sqrt( InitialE*InitialE - M*M) * 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) ) cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles "<<endl;
-
- // Generate an event and store is weight in the Output Tree
- m_EventWeight = TPhaseSpace.Generate() ;
-
-
- TLorentzVector* daugterLV ;
- double Energy;
- G4ThreeVector Momentum;
- if (m_ShootHeavy)
- {
- daugterLV = TPhaseSpace.GetDecay(0);
-
- Momentum = G4ThreeVector( daugterLV->X()*GeV ,
- daugterLV->Y()*GeV ,
- daugterLV->Z()*GeV );
-
- Energy = daugterLV->E()*GeV-md ;
- 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) ;
-
- // get theta and phi in the world frame
- G4double theta_world = Momentum.theta();
- G4double phi_world = Momentum.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);
- m_InitConditions->SetICEmittedEnergy(Energy);
- }
-
- if (m_ShootDecayProduct)
- {
- G4int jj = 1 ;
- for ( int u = 0; u < NumberOfNeutrons ; u++)
- {
- daugterLV = TPhaseSpace.GetDecay(jj);
-
- Momentum = G4ThreeVector( daugterLV->X()*GeV,
- daugterLV->Y()*GeV,
- daugterLV->Z()*GeV);
- Energy = daugterLV->E()*GeV-mn ;
- 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) ;
-
- // get theta and phi in the world frame
- G4double theta_world = Momentum.theta();
- G4double phi_world = Momentum.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);
- m_InitConditions->SetICEmittedEnergy(Energy);
-
- jj++;
- }
-
- for ( int u = 0; u < NumberOfProtons ; u++)
- {
- daugterLV = TPhaseSpace.GetDecay(jj);
-
- Momentum = G4ThreeVector( daugterLV->X()*GeV,
- daugterLV->Y()*GeV,
- daugterLV->Z()*GeV);
- Energy = daugterLV->E()*GeV-mp ;
- 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);
-
- // get theta and phi in the world frame
- G4double theta_world = Momentum.theta();
- G4double phi_world = Momentum.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);
- m_InitConditions->SetICEmittedEnergy(Energy);
- jj++;
- }
-
-
- }
-
- delete masses ;
- }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfertToResonance::SetEverything( string name1,
- string name2,
- string name3,
- string name4,
- double BeamEnergy,
- double ExcitationLight,
- double ExcitationHeavy,
- double BeamEnergySpread,
- double SigmaX,
- double SigmaY,
- double SigmaThetaX,
- double SigmaPhiY,
- double ResonanceWidth,
- int ResonanceDecayZ,
- int ResonanceDecayA,
- bool ShootLight,
- bool ShootHeavy,
- bool ShootDecayProduct,
- string Path,
- double CSThetaMin,
- double CSThetaMax)
-{
- //------------- Constructor with nuclei names and beam energy ------------
-
- m_Reaction = new Reaction ( name1 ,
- name2 ,
- name3 ,
- name4 ,
- BeamEnergy,
- ExcitationLight,
- ExcitationHeavy,
- Path,
- CSThetaMin,
- CSThetaMax) ;
-
- m_BeamEnergy = BeamEnergy;
- m_BeamEnergySpread = BeamEnergySpread;
- m_SigmaX = SigmaX;
- m_SigmaY = SigmaY;
- m_SigmaThetaX = SigmaThetaX;
- m_SigmaPhiY = SigmaPhiY;
- m_ResonanceWidth = ResonanceWidth;
- m_ResonanceMean = ExcitationHeavy;
- m_ResonanceDecayZ = ResonanceDecayZ;
- m_ResonanceDecayA = ResonanceDecayA;
- m_ShootLight = ShootLight;
- m_ShootHeavy = ShootHeavy;
- m_ShootDecayProduct = ShootDecayProduct;
-
-}
-
--
GitLab