diff --git a/Inputs/DetectorConfiguration/hira_upgrade.detector b/Inputs/DetectorConfiguration/hira_upgrade.detector
index df05ab9fc281e93480b277aae94b4bfe72a9db9e..845a13982a9df8dc337c08af2960fc4f8ffeedde 100644
--- a/Inputs/DetectorConfiguration/hira_upgrade.detector
+++ b/Inputs/DetectorConfiguration/hira_upgrade.detector
@@ -15,122 +15,171 @@ HIRAArray
%%%%%%% Telescope0 %%%%%%%
HiraTelescope
- A= -177.657 64.6929 291.606
- B= -223.829 29.3441 262.173
- C= -211.837 -6.11417 273.759
- D= -165.941 29.2929 303.557
- ThinSi_DE= 1
+ A= -156.305 -173.834 263.552
+ B= -188.852 -203.91 215.761
+ C= -163.032 -232.465 208.619
+ D= -130.464 -202.409 256.407
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope1 %%%%%%%
HiraTelescope
- A= -165.967 134.197 274.275
- B= -216.232 99.9158 250.91
- C= -208.327 65.8351 268.555
- D= -158.327 100.269 292.267
- ThinSi_DE= 1
+ A= -173.728 -100.878 289.399
+ B= -210.345 -134.992 247.647
+ C= -188.391 -167.381 246.239
+ D= -151.755 -133.29 287.992
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope2 %%%%%%%
HiraTelescope
- A= -146.452 197.836 245.343
- B= -200.596 166.121 227.729
- C= -196.529 134.907 251.064
- D= -142.626 166.862 268.991
- ThinSi_DE= 1
+ A= -181.571 -22.7927 301.035
+ B= -222.662 -59.2104 265.919
+ C= -204.954 -93.7872 270.812
+ D= -163.848 -57.3939 305.933
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
-%%%%%%% Telescope4 %%%%%%%
+%%%%%%% Telescope3 %%%%%%%
HiraTelescope
- A= -247.076 64.6929 235.688
- B= -284.058 29.3441 195.308
- C= -275.473 -6.11417 209.604
- D= -238.853 29.2929 250.264
- ThinSi_DE= 1
+ A= -179.436 56.4514 297.868
+ B= -225.175 19.5814 269.647
+ C= -211.881 -15.4248 281.087
+ D= -166.129 21.4205 309.318
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
-%%%%%%% Telescope5 %%%%%%%
+%%%%%%% Telescope4 %%%%%%%
HiraTelescope
- A= -231.3 134.197 221.974
- B= -273.805 99.9158 186.395
- C= -270.736 65.8351 205.485
- D= -228.576 100.269 241.33
- ThinSi_DE= 1
+ A= -229.126 -173.834 203.451
+ B= -245.577 -203.91 148.019
+ C= -218.809 -232.465 149.078
+ D= -202.337 -202.409 204.513
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
-%%%%%%% Telescope6 %%%%%%%
+%%%%%%% Telescope5 %%%%%%%
HiraTelescope
- A= -204.962 197.836 199.079
- B= -252.702 166.121 168.051
- C= -254.812 134.907 191.643
- D= -207.386 166.862 222.911
- ThinSi_DE= 1
+ A= -253.591 -100.878 222.766
+ B= -275.758 -134.992 171.848
+ C= -254.417 -167.381 177.191
+ D= -232.233 -133.29 228.116
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
-%%%%%%% Telescope7 %%%%%%%
+%%%%%%% Telescope6 %%%%%%%
HiraTelescope
- A= -169.213 252.829 168.003
- B= -221.672 225.065 141.077
- C= -228.398 198.083 168.682
- D= -176.21 226.163 195.81
- ThinSi_DE= 1
+ A= -264.604 -22.7927 231.462
+ B= -293.053 -59.2104 185.502
+ C= -277.676 -93.7872 195.554
+ D= -249.214 -57.3939 241.523
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope7 %%%%%%%
HiraTelescope
- A= -125.617 296.772 130.106
- B= -182.071 274.173 106.652
- C= -192.649 252.602 137.605
- D= -136.411 275.579 161.212
- ThinSi_DE= 1
+ A= -261.607 56.4514 229.095
+ B= -296.582 19.5814 188.289
+ C= -287.402 -15.4248 203.233
+ D= -252.418 21.4205 244.053
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope8 %%%%%%%
HiraTelescope
- A= -305.189 64.6929 153.152
- B= -328.749 29.3441 103.724
- C= -324.718 -6.11417 119.904
- D= -301.587 29.2929 169.495
- ThinSi_DE= 1
+ A= -292.743 -139.329 137.835
+ B= -296.229 -171.569 81.2998
+ C= -273.828 -202.179 91.0104
+ D= -270.325 -169.961 147.556
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope9 %%%%%%%
HiraTelescope
- A= -286.092 134.197 144.649
- B= -316.337 99.9158 98.1979
- C= -318.984 65.8351 117.351
- D= -289.147 100.269 163.956
- ThinSi_DE= 1
+ A= -314.059 -63.4009 146.481
+ B= -324.663 -98.8582 92.8326
+ C= -309.021 -132.523 105.284
+ D= -298.405 -97.0896 158.945
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope10 %%%%%%%
HiraTelescope
- A= -254.211 197.836 130.455
- B= -290.793 166.121 86.8249
- C= -299.709 134.907 108.77
- D= -263.497 166.862 152.536
- ThinSi_DE= 1
+ A= -319.069 15.7508 148.513
+ B= -337.278 -21.1212 97.9493
+ C= -328.857 -56.1293 113.33
+ D= -310.641 -19.282 163.908
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
%%%%%%% Telescope11 %%%%%%%
HiraTelescope
- A= -210.939 252.829 111.189
- B= -253.233 225.065 70.1018
- C= -267.736 198.083 94.5344
- D= -225.76 226.163 135.735
- ThinSi_DE= 1
+ A= -307.519 94.1016 143.828
+ B= -333.433 57.6897 96.3898
+ C= -332.328 23.1183 114.738
+ D= -306.411 59.5059 162.192
+ ThinSi_DE= 0
ThickSi_E= 1
CsI= 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Microball
+ RING1= 0
+ RING2= 0
+ RING3= 1
+ RING4= 1
+ RING5= 1
+ RING6= 1
+ RING7= 1
+ RING8= 1
+ RING9= 1
+ %DISABLE_CRYSTAL 19
+ DISABLE_CRYSTAL 20
+ DISABLE_CRYSTAL 21
+ DISABLE_CRYSTAL 22
+ %DISABLE_CRYSTAL 31
+ DISABLE_CRYSTAL 32
+ DISABLE_CRYSTAL 33
+ DISABLE_CRYSTAL 34
+ %DISABLE_CRYSTAL 43
+ DISABLE_CRYSTAL 44
+ DISABLE_CRYSTAL 45
+ DISABLE_CRYSTAL 46
+ DETECTOR_FLIP= 0
+ INCLUDE_CHAMBER= 0
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+NeutronWall
+ THETA= 50
+ PHI= 0
+ R= 4000
+ BARS= 25
+ VETOWALL= 1
+ VWDISTANCE= 429
+ OVERLAP= 3
+ VWMATERIAL= BC400
+ NWMATERIAL= NE213
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+NeutronWall
+ THETA= 50
+ PHI= 0
+ R= 4573
+ BARS= 25
+ VETOWALL= 0
+ VWDISTANCE= 300
+ OVERLAP= 3
+ VWMATERIAL= BC400
+ NWMATERIAL= NE213
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/NPLib/Core/NPGlobalSystemOfUnits.h b/NPLib/Core/NPGlobalSystemOfUnits.h
index cbce24d95d96d6e22d283e5906017262fdfcbdf7..3e53094b3b924d2d22d61250101af91455b42b08 100644
--- a/NPLib/Core/NPGlobalSystemOfUnits.h
+++ b/NPLib/Core/NPGlobalSystemOfUnits.h
@@ -79,6 +79,7 @@ using NPUNITS::megahertz;
using NPUNITS::ns;
using NPUNITS::s;
using NPUNITS::ms;
+using NPUNITS::us;
using NPUNITS::eplus;
using NPUNITS::e_SI;
using NPUNITS::coulomb;
diff --git a/NPLib/Core/NPInputParser.cxx b/NPLib/Core/NPInputParser.cxx
index 1a8f1070fd1882b92f978b7a7f45033c70eab7a8..af1c2b5572b37e8f695702aaca7e58b5e813f8a3 100644
--- a/NPLib/Core/NPInputParser.cxx
+++ b/NPLib/Core/NPInputParser.cxx
@@ -580,8 +580,8 @@ double NPL::ApplyUnit(double value, std::string unit){
return value*NPUNITS::ns;
}
- else if(unit=="s"){
- return value*NPUNITS::s;
+ else if(unit=="us"){
+ return value*NPUNITS::us;
}
else if(unit=="ms"){
diff --git a/NPLib/Core/NPSystemOfUnits.h b/NPLib/Core/NPSystemOfUnits.h
index 2c5c615ec1f1328566b75898c4013be9dc15dbe1..d6513b20982f7c315e2de09104b38929edf6a0f0 100644
--- a/NPLib/Core/NPSystemOfUnits.h
+++ b/NPLib/Core/NPSystemOfUnits.h
@@ -110,6 +110,7 @@ namespace NPUNITS {
// symbols
static const double ns = nanosecond;
+ static const double us = microsecond;
static const double s = second;
static const double ms = millisecond;
diff --git a/NPLib/Detectors/Chio/TChio_anData.cxx b/NPLib/Detectors/Chio/TChio_anData.cxx
index 12d3fbb9d184e60b21fdb904d5af84f8cb4b29d5..7c38f0d888fc56a64b20d51c1caeadfe104a62d5 100644
--- a/NPLib/Detectors/Chio/TChio_anData.cxx
+++ b/NPLib/Detectors/Chio/TChio_anData.cxx
@@ -27,28 +27,19 @@ using namespace std;
ClassImp(TChio_anData)
-TChio_anData::TChio_anData()
-{
- // Default constructor
-
- // (E)
- fChio_an_Energy.clear();
- fChio_an_Energy_pileup.clear();
-}
+TChio_anData::TChio_anData(){
+ }
-TChio_anData::~TChio_anData()
-{
+TChio_anData::~TChio_anData(){
}
-
-void TChio_anData::Clear()
-{
- // (E)
+void TChio_anData::Clear(){
fChio_an_Energy.clear();
- fChio_an_Energy_pileup.clear();
+ fChio_an_Time.clear();
+ fChio_an_DetectorNbr.clear();
}
@@ -56,14 +47,12 @@ void TChio_anData::Clear()
void TChio_anData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
-
// Chio_an
// (E)
cout << "Chio_an_MultE = " << fChio_an_Energy.size() << endl;
- for (UShort_t i = 0; i < fChio_an_Energy.size(); i++)
- cout << " Energy: " << fChio_an_Energy[i] << endl;
-
- cout << "Chio_an_MultE pileup = " << fChio_an_Energy_pileup.size() << endl;
- for (UShort_t i = 0; i < fChio_an_Energy_pileup.size(); i++)
- cout << " Energy: " << fChio_an_Energy_pileup[i] << endl;
+ for (UShort_t i = 0; i < fChio_an_Energy.size(); i++){
+ cout << " Detector: " << fChio_an_DetectorNbr[i] ;
+ cout << " Energy: " << fChio_an_Energy[i];
+ cout << " Time: " << fChio_an_Time[i] << endl;
+ }
}
diff --git a/NPLib/Detectors/Chio/TChio_anData.h b/NPLib/Detectors/Chio/TChio_anData.h
index 808e6647765a96e6a7e55b7d0a89c44b1c61a94c..ce0160bba3e92f858322adfaa7a3858eaf2515b3 100644
--- a/NPLib/Detectors/Chio/TChio_anData.h
+++ b/NPLib/Detectors/Chio/TChio_anData.h
@@ -28,30 +28,34 @@ using namespace std;
class TChio_anData : public TObject {
private:
// ADC
- vector<UShort_t> fChio_an_Energy;
- vector<UShort_t> fChio_an_Energy_pileup;
+ vector<double> fChio_an_Energy;
+ vector<double> fChio_an_Time;
+ vector<double> fChio_an_DetectorNbr;
public:
TChio_anData();
virtual ~TChio_anData();
void Clear();
- void Clear(const Option_t*) {};
+ void Clear(const Option_t*) {};
void Dump() const;
///////////////////// GETTERS ////////////////////////
// (E)
- UShort_t GetMultE() {return fChio_an_Energy.size();}
- UShort_t GetEnergy(Int_t i) {return fChio_an_Energy.at(i);}
- UShort_t GetMultE_pileup() {return fChio_an_Energy_pileup.size();}
- UShort_t GetEnergy_pileup(Int_t i) {return fChio_an_Energy_pileup.at(i);}
+ UShort_t GetMult() {return fChio_an_Energy.size();}
+ UShort_t GetEnergy(unsigned int i) {return fChio_an_Energy[i];}
+ UShort_t GetTime(unsigned int i) {return fChio_an_Time[i];}
+ UShort_t GetDetectorNbr(unsigned i) {return fChio_an_DetectorNbr[i];}
///////////////////// SETTERS ////////////////////////
// (E)
- void SetEnergy(UShort_t E) {fChio_an_Energy.push_back(E);}
- void SetEnergy_pileup(UShort_t E) {fChio_an_Energy_pileup.push_back(E);}
+ inline void SetEnergyAndTime(double E, double T, unsigned int Det) {
+ fChio_an_Energy.push_back(E);
+ fChio_an_Time.push_back(T);
+ fChio_an_DetectorNbr.push_back(Det);
+ }
- ClassDef(TChio_anData,1) // Chio_anData structure
+ ClassDef(TChio_anData,2) // Chio_anData structure
};
#endif
diff --git a/NPLib/Detectors/Chio/TChio_anPhysics.cxx b/NPLib/Detectors/Chio/TChio_anPhysics.cxx
index 09edc31e779cfb6c72c19573d79ccea5661ee5f7..4d060a2a74e1efe340128275136f422aa850bb8b 100644
--- a/NPLib/Detectors/Chio/TChio_anPhysics.cxx
+++ b/NPLib/Detectors/Chio/TChio_anPhysics.cxx
@@ -101,9 +101,9 @@ void TChio_anPhysics::AddParameterToCalibrationManager()
void TChio_anPhysics::InitializeRootInputRaw()
{
TChain* inputChain = RootInput::getInstance()->GetChain() ;
- inputChain->SetBranchStatus("CHIO_AN" , true) ;
+ inputChain->SetBranchStatus("ChioAn" , true) ;
inputChain->SetBranchStatus("fChio_an_*" , true) ;
- inputChain->SetBranchAddress("CHIO_AN" , &EventData) ;
+ inputChain->SetBranchAddress("ChioAn" , &EventData) ;
}
// Activated associated Branches and link it to the private member DetectorPhysics address
diff --git a/NPLib/Detectors/Chio/TChio_digData.cxx b/NPLib/Detectors/Chio/TChio_digData.cxx
index 5bc21ed61fb5c3b2f9c35f2e60bac1daf2f5671c..6d6fa2d36d51085aa1f76b3f3ba025b2e840ddec 100644
--- a/NPLib/Detectors/Chio/TChio_digData.cxx
+++ b/NPLib/Detectors/Chio/TChio_digData.cxx
@@ -26,43 +26,39 @@ using namespace std;
ClassImp(TChio_digData)
-
+////////////////////////////////////////////////////////////////////////////////
TChio_digData::TChio_digData()
{
- // Default constructor
-
- // (E)
- fChio_dig_Energy.clear();
- //sum
- fChio_dig_Sum = 0;
}
-
+////////////////////////////////////////////////////////////////////////////////
TChio_digData::~TChio_digData()
{
}
-
+////////////////////////////////////////////////////////////////////////////////
void TChio_digData::Clear()
{
- // (E)
fChio_dig_Energy.clear();
- //sum
- fChio_dig_Sum = 0;
-}
+ fChio_dig_Time.clear();
+}
+////////////////////////////////////////////////////////////////////////////////
void TChio_digData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
// (E)
cout << "Chio_dig_Esize = " << fChio_dig_Energy.size() << endl;
- for (UShort_t i = 0; i < fChio_dig_Energy.size(); i++)
- cout << " Energy: " << fChio_dig_Energy[i] << endl;
- //sum
- cout<<"Chio_dig_sum" << fChio_dig_Sum << endl;
+}
+////////////////////////////////////////////////////////////////////////////////
+TGraph* TChio_digData::GetEnergyAsGraph(){
+ TGraph* res = new TGraph (fChio_dig_Energy.size(),&fChio_dig_Time[0],&fChio_dig_Energy[0]);
+ res->Sort();
+ return res;
+
}
diff --git a/NPLib/Detectors/Chio/TChio_digData.h b/NPLib/Detectors/Chio/TChio_digData.h
index 3c1069acc41100d8c7aeea1e9619caf17caceb08..fe1d56585550544b4a5f4dfa2df0187d9cdf4b7d 100644
--- a/NPLib/Detectors/Chio/TChio_digData.h
+++ b/NPLib/Detectors/Chio/TChio_digData.h
@@ -3,16 +3,16 @@
#include <vector>
#include "TObject.h"
+#include "TGraph.h"
using namespace std;
class TChio_digData : public TObject {
private:
// ADC
- vector<UShort_t> fChio_dig_Energy;
+ vector<double> fChio_dig_Energy;
+ vector<double> fChio_dig_Time;
- //sum
- Int_t fChio_dig_Sum;
public:
TChio_digData();
@@ -24,19 +24,18 @@ class TChio_digData : public TObject {
///////////////////// GETTERS ////////////////////////
// (E)
- UShort_t GetEsize() {return fChio_dig_Energy.size();}
- UShort_t GetEnergy(Int_t i) {return fChio_dig_Energy.at(i);}
-
- //sum
- Int_t GetSum() {return fChio_dig_Sum;}
+ inline unsigned int GetEsize() {return fChio_dig_Energy.size();}
+ inline vector<double> GetEnergy() {return fChio_dig_Energy;}
+ TGraph* GetEnergyAsGraph();
///////////////////// SETTERS ////////////////////////
- //sum
- void SetSum(UShort_t E) {fChio_dig_Sum += (Int_t)E;}
// (E)
- void SetEnergy(UShort_t E) {fChio_dig_Energy.push_back(E);}
+ void SetEnergy(vector<double>& E) {fChio_dig_Energy=E;}
+ void AddEnergyPoint(double& E,double& T){
+ fChio_dig_Energy.push_back(E);
+ fChio_dig_Time.push_back(T);}
- ClassDef(TChio_digData,1) // Chio_digData structure
+ ClassDef(TChio_digData,2) // Chio_digData structure
};
#endif
diff --git a/NPLib/Detectors/Chio/TChio_digPhysics.cxx b/NPLib/Detectors/Chio/TChio_digPhysics.cxx
index 1c9c2e64db0d800b85262b540bbbaa05bfe2c907..ab0923212f1633dcffff615dc7d4a2196b0ea731 100644
--- a/NPLib/Detectors/Chio/TChio_digPhysics.cxx
+++ b/NPLib/Detectors/Chio/TChio_digPhysics.cxx
@@ -105,9 +105,9 @@ void TChio_digPhysics::AddParameterToCalibrationManager()
void TChio_digPhysics::InitializeRootInputRaw()
{
TChain* inputChain = RootInput::getInstance()->GetChain() ;
- inputChain->SetBranchStatus("CHIO_DIG" , true) ;
+ inputChain->SetBranchStatus("ChioDig" , true) ;
inputChain->SetBranchStatus("fChio_dig_*" , true) ;
- inputChain->SetBranchAddress("CHIO_DIG" , &EventData) ;
+ inputChain->SetBranchAddress("ChioDig" , &EventData) ;
}
// Activated associated Branches and link it to the private member DetectorPhysics address
@@ -243,12 +243,7 @@ void TChio_digPhysics::BuildSimplePhysicalEvent()
// this value is defined by RC cuircuit in pre-amplifier.
// ask technitian for R and C values
- // read digitizer data into "Energy"
- for(int ch=0;ch<Nch;ch++){
- Energy.push_back((double)EventData->GetEnergy(ch));
- // Energy[ch] = (int)EventData->GetEnergy(ch);
- // cout << EventData->GetEnergy(ch) << endl;
- }
+ Energy = EventData->GetEnergy();
if(Nch!=0){
int ch = 0;
diff --git a/NPLib/Detectors/GASPARD/GaspardTrackerTrapezoid.cxx b/NPLib/Detectors/GASPARD/GaspardTrackerTrapezoid.cxx
index 9402c2fe543bcd78fafe31796200bd6ebc41aec5..82596e18970b93e32e7fba13744cc6e525c6c17c 100644
--- a/NPLib/Detectors/GASPARD/GaspardTrackerTrapezoid.cxx
+++ b/NPLib/Detectors/GASPARD/GaspardTrackerTrapezoid.cxx
@@ -30,10 +30,13 @@
#include <cmath>
#include <stdlib.h>
+// NPLib
+#include "NPSystemOfUnits.h"
// Gaspard
#include "TGaspardTrackerPhysics.h"
-
+//Root
+#include"TRotation.h"
GaspardTrackerTrapezoid::GaspardTrackerTrapezoid(map<int, GaspardTrackerModule*> &Module,
TGaspardTrackerPhysics* EventPhysics)
: m_ModuleTest(Module),
@@ -87,7 +90,7 @@ void GaspardTrackerTrapezoid::ReadConfiguration(NPL::InputParser parser){
double Phi = blocks[i]->GetDouble("PHI","deg");
vector<double> beta = blocks[i]->GetVectorDouble("BETA","deg");
- AddModule(R,Theta,Phi,beta[0],beta[1],beta[2]);
+ AddModule(Theta,Phi,R,beta[0],beta[1],beta[2]);
m_ModuleTest[m_index["Trapezoid"] + m_NumberOfModule] = this;
}
}
@@ -289,10 +292,6 @@ void GaspardTrackerTrapezoid::AddModule(double theta,
{
m_NumberOfModule++;
- // convert from degree to radian:
- theta *= M_PI/180.;
- phi *= M_PI/180.;
-
// Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
TVector3 U ;
// Vector V on Module Face (parallele to X Strip)
@@ -312,21 +311,21 @@ void GaspardTrackerTrapezoid::AddModule(double theta,
W = C.Unit();
V = W.Cross(YperpW);
- U = W.Cross(U);
+ U = V.Cross(W);
U = U.Unit();
V = V.Unit();
- U.Rotate( beta_u * M_PI/180. , U ) ;
- V.Rotate( beta_u * M_PI/180. , U ) ;
-
- U.Rotate( beta_v * M_PI/180. , V ) ;
- V.Rotate( beta_v * M_PI/180. , V ) ;
-
- U.Rotate( beta_w * M_PI/180. , W ) ;
- V.Rotate( beta_w * M_PI/180. , W ) ;
+ C=W*10*NPUNITS::mm+C;
+ TRotation R;
+ R.Rotate(beta_u,U);
+ R.Rotate(beta_v,V);
+ R.Rotate(beta_w,W);
+ //W=R*W;
+ V=R*V;
+ U=R*U;
vector<double> lineX;
vector<double> lineY;
vector<double> lineZ;
@@ -338,19 +337,20 @@ void GaspardTrackerTrapezoid::AddModule(double theta,
double X, Y, Z;
// Moving C to the 1.1 corner:
- C.SetX( C.X() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.X() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.X() ) ;
- C.SetY( C.Y() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Y() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Y() ) ;
- C.SetZ( C.Z() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Z() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Z() ) ;
-
+ TVector3 Strip_1_1;
+ Strip_1_1.SetX( C.X() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.X() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.X() ) ;
+ Strip_1_1.SetY( C.Y() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Y() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Y() ) ;
+ Strip_1_1.SetZ( C.Z() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Z() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Z() ) ;
+
for (int i = 0; i < m_NumberOfStripsX; i++) {
lineX.clear();
lineY.clear();
lineZ.clear();
for (int j = 0; j < m_NumberOfStripsY; j++) {
- X = C.X() + m_StripPitchX*U.X()*i + m_StripPitchY*V.X()*j ;
- Y = C.Y() + m_StripPitchX*U.Y()*i + m_StripPitchY*V.Y()*j ;
- Z = C.Z() + m_StripPitchX*U.Z()*i + m_StripPitchY*V.Z()*j ;
+ X = Strip_1_1.X() + m_StripPitchX*U.X()*i + m_StripPitchY*V.X()*j ;
+ Y = Strip_1_1.Y() + m_StripPitchX*U.Y()*i + m_StripPitchY*V.Y()*j ;
+ Z = Strip_1_1.Z() + m_StripPitchX*U.Z()*i + m_StripPitchY*V.Z()*j ;
lineX.push_back(X);
lineY.push_back(Y);
diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc
index 4ed7fa3ee315ade8898132a74658c6cdc1f21610..a02d0cdee388aa7bf2cb7180891552ee83eccdc0 100644
--- a/NPSimulation/Detectors/Chio/Chio.cc
+++ b/NPSimulation/Detectors/Chio/Chio.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009-2106 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2016 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 *
@@ -52,7 +52,7 @@
// NPTool header
#include "Chio.hh"
-#include "CalorimeterScorers.hh"
+#include "DriftElectronScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
#include "NPSDetectorFactory.hh"
@@ -61,6 +61,10 @@
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
+// ROOT
+#include "TH1D.h"
+#include "TF1.h"
+
using namespace std;
using namespace CLHEP;
@@ -68,11 +72,11 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Chio_NS{
// Energy and time Resolution
-// const double EnergyThreshold = 0.1*MeV;
- //const double ResoTime = 4.5*ns ;
- // const double ResoEnergy = 1.0*MeV ;
+ const double ChargeThreshold = 1;
+ const double ResoTime = 4.5*ns ;
+ // const double ResoEnergy = 1.0*MeV ;
//const double Radius = 50*mm ;
- // const double Width = 100*mm ;
+ // const double Width = 100*mm ;
const double Thickness = 300*mm ;
const string Material = "BC400";
}
@@ -81,7 +85,9 @@ namespace Chio_NS{
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Chio Specific Method
Chio::Chio(){
- m_Event = new TChio_anData() ;
+ m_Event_an = new TChio_anData() ;
+ m_Event_dig = new TChio_digData() ;
+
m_ChioScorer = 0;
m_SquareDetector = 0;
m_CylindricalDetector = 0;
@@ -125,27 +131,48 @@ G4LogicalVolume* Chio::BuildDetector(){
G4Box* sGas = new G4Box("Chio_Gas",6*cm*0.5,
6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5);
+ // Frish grid
+ G4Box* sGrid = new G4Box("Chio_Grid",1*um*0.5,
+ 6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5);
+
+ // Cathode
+ G4Box* sCathode = new G4Box("Chio_Cathode",1*um*0.5,
+ 6*cm*0.5,12*cm*0.5-2*12*micrometer*0.5);
+
+
// Entrance/Exit windows
G4Box* sWindows = new G4Box("Chio_Windows",7*cm*0.5,
7*cm*0.5,12*micrometer*0.5);
G4Material* Fe= MaterialManager::getInstance()->GetMaterialFromLibrary("Fe");
+ G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+
G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin);
G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
-
+
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
- MPT->AddConstProperty("DE_YIELD",1e-3);
-// MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
- MPT->AddConstProperty("DE_ABSLENGTH",1*km);
- MPT->AddConstProperty("DE_DRIFTSPEED",8e-3*mm/ns);
+ MPT->AddConstProperty("DE_YIELD",1e-2);
+ // MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
+ MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
+ MPT->AddConstProperty("DE_DRIFTSPEED",11*cm/microsecond);
MPT->AddConstProperty("DE_TRANSVERSALSPREAD",6e-5*mm2/ns);
MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-5*mm2/ns);
CF4->SetMaterialPropertiesTable(MPT);
+
+ G4MaterialPropertiesTable* MPT2 = new G4MaterialPropertiesTable();
+ MPT2->AddConstProperty("DE_YIELD",1);
+ MPT2->AddConstProperty("DE_AMPLIFICATION",2);
+ MPT2->AddConstProperty("DE_ABSLENGTH",1*pc);
+
+ Al->SetMaterialPropertiesTable(MPT2);
+
m_SquareDetector = new G4LogicalVolume(sChamber,Fe,"logic_Chio_Box",0,0,0);
G4LogicalVolume* logicGas = new G4LogicalVolume(sGas,CF4,"logic_Gas",0,0,0);
+ G4LogicalVolume* logicGrid = new G4LogicalVolume(sGrid,Al,"logic_Grid",0,0,0);
+ G4LogicalVolume* logicCathode = new G4LogicalVolume(sCathode,Fe,"logic_Cathode",0,0,0);
G4LogicalVolume* logicWindows = new G4LogicalVolume(sWindows,Mylar,"logic_Windows",0,0,0);
G4RotationMatrix* Rot = new G4RotationMatrix();
@@ -154,6 +181,14 @@ G4LogicalVolume* Chio::BuildDetector(){
logicGas,
"ChioGas",m_SquareDetector,false,0);
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(-2.5*cm,0,0)),
+ logicGrid,
+ "ChioGrid",logicGas,false,0);
+
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)),
+ logicCathode,
+ "ChioCathode",logicGas,false,0);
+
new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(0,0,6*cm-6*micrometer)),
logicWindows,
"ChioExitWindows",m_SquareDetector,false,0);
@@ -162,14 +197,25 @@ G4LogicalVolume* Chio::BuildDetector(){
logicWindows,
"ChioEntranceWindows",m_SquareDetector,false,0);
+ G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,-1000*volt/cm,0.0));
+ // Create an equation of motion for this field
+ G4EqMagElectricField* Equation = new G4EqMagElectricField(field);
+ G4MagIntegratorStepper* Stepper = new G4ClassicalRK4( Equation, 8 );
+
+ // Get the global field manager
+ G4FieldManager* FieldManager= new G4FieldManager();
+ // Set this field to the global field manager
+ FieldManager->SetDetectorField(field );
+ logicGas->SetFieldManager(FieldManager,true);
-/* G4Region* DriftRegion = new G4Region("DriftRegion");
- DriftRegion->AddRootLogicalVolume(logicGas);
- G4ProductionCuts* cuts = new G4ProductionCuts;
- cuts->SetProductionCut(1e-9*micrometer); // same cuts for gamma, e- and e+
- DriftRegion->SetProductionCuts(cuts);
-
- new DriftElectron("DriftElectron",DriftRegion);*/
+ G4MagInt_Driver* IntgrDriver = new G4MagInt_Driver(0.1*mm,
+ Stepper,
+ Stepper->GetNumberOfVariables() );
+
+ G4ChordFinder* ChordFinder = new G4ChordFinder(IntgrDriver);
+ FieldManager->SetChordFinder( ChordFinder );
+
+ logicCathode->SetSensitiveDetector(m_ChioScorer);
m_SquareDetector->SetVisAttributes(m_VisChamber);
logicGas->SetVisAttributes(m_VisGas);
logicWindows->SetVisAttributes(m_VisWindows);
@@ -258,40 +304,106 @@ void Chio::ConstructDetector(G4LogicalVolume* world){
void Chio::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("Chio")){
- pTree->Branch("Chio", "TChioData", &m_Event) ;
+ if(!pTree->FindBranch("ChioAn")){
+ pTree->Branch("ChioAn", "TChio_anData", &m_Event_an) ;
}
- pTree->SetBranchAddress("Chio", &m_Event) ;
+ pTree->SetBranchAddress("ChioAn", &m_Event_an) ;
+
+ ///////////////
+ if(!pTree->FindBranch("ChioDig")){
+ pTree->Branch("ChioDig", "TChio_digData", &m_Event_dig) ;
+ }
+ pTree->SetBranchAddress("ChioDig", &m_Event_dig) ;
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void Chio::ReadSensitive(const G4Event* event){
- m_Event->Clear();
+ m_Event_an->Clear();
///////////
- // Calorimeter scorer
- NPS::HitsMap<G4double*>* CaloHitMap;
- std::map<G4int, G4double**>::iterator Calo_itr;
-
- G4int CaloCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Calorimeter");
- CaloHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CaloCollectionID));
- /*
- // Loop on the Calo map
- for (Calo_itr = CaloHitMap->GetMap()->begin() ; Calo_itr != CaloHitMap->GetMap()->end() ; Calo_itr++){
-
- G4double* Info = *(Calo_itr->second);
- double Energy = RandGauss::shoot(Info[0],Chio_NS::ResoEnergy);
- if(Energy>Chio_NS::EnergyThreshold){
- double Time = RandGauss::shoot(Info[1],Chio_NS::ResoTime);
- int DetectorNbr = (int) Info[2];
- //m_Event->SetEnergy(DetectorNbr,Energy);
- //m_Event->SetTime(DetectorNbr,Time);
+ // Cathode analogic scorer
+ NPS::HitsMap<G4double*>* CathodeHitMap;
+ std::map<G4int, G4double**>::iterator Cathode_itr;
+
+ G4int CathodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Cathode_an");
+ CathodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CathodeCollectionID));
+
+ // Loop on the Cathode map
+ for (Cathode_itr = CathodeHitMap->GetMap()->begin() ; Cathode_itr != CathodeHitMap->GetMap()->end() ; Cathode_itr++){
+ G4double* Info = *(Cathode_itr->second);
+ double Count= Info[0];
+ if(Count>Chio_NS::ChargeThreshold-1){
+ double Time = RandGauss::shoot(Info[1],Chio_NS::ResoTime);
+ m_Event_an->SetEnergyAndTime(Count,Time,Info[2]);
+ }
+ }
+ // clear map for next event
+ CathodeHitMap->clear();
+
+ m_Event_dig->Clear();
+
+ ///////////
+ CathodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Cathode_dig");
+ CathodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CathodeCollectionID));
+
+ // Loop on the Cathode map
+ TH1D* h = new TH1D("h","h",25000,0,25000);
+ for (Cathode_itr = CathodeHitMap->GetMap()->begin() ; Cathode_itr != CathodeHitMap->GetMap()->end() ; Cathode_itr++){
+ G4double* Info = *(Cathode_itr->second);
+ if(Info[0]){
+ h->Fill(Info[1],Info[0]);
+ }
}
- }*/
+
+ vector<double> E,T;
+ for(int i = 0 ; i < h->GetNbinsX() ; i++){
+ double count = h->GetBinContent(i);
+ double time = h->GetBinCenter(i);
+ if(count)
+ // m_Event_dig->AddEnergyPoint(count,time);
+ E.push_back(count);
+ T.push_back(time+500);
+ }
+
+ SimulateDigitizer(E,T,1.40*microsecond,0,8750,25,5);
+
+ delete h;
// clear map for next event
- CaloHitMap->clear();
+ CathodeHitMap->clear();
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Chio::SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop, double step,double noise){
+
+ static string formula;
+ formula= "";
+ static string Es,Ts,var,cond;
+ static string fall;
+ fall=std::to_string(fallTime);
+
+ for(unsigned int i = 0 ; i < E.size() ; i++){
+ if(E[i]!=0 && T[i]!=0){
+ Es = std::to_string(E[i]);
+ Ts = std::to_string(T[i]);
+ cond = ")*(x>"+Ts+")+";
+ var = "(x-"+Ts+")";
+ formula += Es+"*-1*exp(-"+var+"/"+fall+cond;
+ }
+ }
+ formula+="0";
+ TF1* f = new TF1("f",formula.c_str(),start,stop);
+ unsigned int size = (stop-start)/step;
+ for(unsigned int i = 0 ; i < size ; i++){
+ double time = start+i*step;
+ double energy = f->Eval(time)+noise*(1-2*G4UniformRand());
+ m_Event_dig->AddEnergyPoint(energy,time);
+ }
+
+ delete f;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -305,10 +417,13 @@ void Chio::InitializeScorers() {
return ;
// Otherwise the scorer is initialised
- vector<int> level; level.push_back(0);
- G4VPrimitiveScorer* Calorimeter= new CALORIMETERSCORERS::PS_Calorimeter("Calorimeter",level, 0) ;
+ G4VPrimitiveScorer* Cathode_an= new DRIFTELECTRONSCORERS::PS_DECathode("Cathode_an",0) ;
+ G4VPrimitiveScorer* Cathode_dig= new DRIFTELECTRONSCORERS::PS_DEDigitizer("Cathode_dig",0) ;
+
//and register it to the multifunctionnal detector
- m_ChioScorer->RegisterPrimitive(Calorimeter);
+ m_ChioScorer->RegisterPrimitive(Cathode_an);
+ m_ChioScorer->RegisterPrimitive(Cathode_dig);
+
G4SDManager::GetSDMpointer()->AddNewDetector(m_ChioScorer) ;
}
diff --git a/NPSimulation/Detectors/Chio/Chio.hh b/NPSimulation/Detectors/Chio/Chio.hh
index 65ab3a9cf14e5c8f98b55dd65f0fe71f3ddb0563..6aa7be66b4ba30f7abc09b9864b4cdf0872f6d60 100644
--- a/NPSimulation/Detectors/Chio/Chio.hh
+++ b/NPSimulation/Detectors/Chio/Chio.hh
@@ -1,7 +1,7 @@
#ifndef Chio_h
#define Chio_h 1
/*****************************************************************************
- * Copyright (C) 2009-XYEARX this file is part of the NPTool Project *
+ * Copyright (C) 2009-2017 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 *
@@ -35,6 +35,8 @@ using namespace std;
// NPTool header
#include "NPSVDetector.hh"
#include "TChio_anData.h"
+#include "TChio_digData.h"
+
#include "NPInputParser.h"
class Chio : public NPS::VDetector{
@@ -84,6 +86,7 @@ class Chio : public NPS::VDetector{
public: // Scorer
// Initialize all Scorer used by the MUST2Array
void InitializeScorers() ;
+ void SimulateDigitizer(vector<double> E, vector<double> T, double fallTime,double start,double stop,double step,double noise);
// Associated Scorer
G4MultiFunctionalDetector* m_ChioScorer ;
@@ -91,7 +94,8 @@ class Chio : public NPS::VDetector{
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
private:
- TChio_anData* m_Event ;
+ TChio_anData* m_Event_an ;
+ TChio_digData* m_Event_dig ;
////////////////////////////////////////////////////
///////////////Private intern Data//////////////////
diff --git a/NPSimulation/Detectors/GASPARD/GaspardTracker.cc b/NPSimulation/Detectors/GASPARD/GaspardTracker.cc
index b8840de6fa156c82b6a6f85b40c068ce609c9cfa..78d58edbba0676fb3f872db23ede91545e94b59d 100644
--- a/NPSimulation/Detectors/GASPARD/GaspardTracker.cc
+++ b/NPSimulation/Detectors/GASPARD/GaspardTracker.cc
@@ -199,15 +199,9 @@ void GaspardTracker::InitializeScorers()
// Called at in the EventAction::EndOfEventAction
void GaspardTracker::ReadSensitive(const G4Event* event)
{
- // Before looping on each sub-detector, clear the static variable
- // ms_InterCoord
- // This is done on the first element of the m_Modules vector.
- // This should be done here since this variable (of type TIneractionCoordinates)
- // deals with multiplicity of events > 1.
- m_Modules[0]->GetInterCoordPointer()->Clear();
-
- // We do the same for the static variable ms_Event
- m_Modules[0]->GetEventPointer()->Clear();
+ // Before looping on each sub-detector, clear the static variable ms_Event
+ if(m_Modules.size()>0)
+ m_Modules[0]->GetEventPointer()->Clear();
// loop on sub-detectors belonging to GaspardTracker
int nbDetectors = m_Modules.size();
diff --git a/NPSimulation/Detectors/GASPARD/GaspardTrackerTrapezoid.cc b/NPSimulation/Detectors/GASPARD/GaspardTrackerTrapezoid.cc
index f3e93267cdecbaa844d1b5bf9646d6ce098da5c6..a95ce30d5f89be461fa9a04936fbe513c853b035 100644
--- a/NPSimulation/Detectors/GASPARD/GaspardTrackerTrapezoid.cc
+++ b/NPSimulation/Detectors/GASPARD/GaspardTrackerTrapezoid.cc
@@ -396,9 +396,9 @@ void GaspardTrackerTrapezoid::ConstructDetector(G4LogicalVolume* world)
// w perpendicular to (u,v) plan and pointing ThirdStage
// Phi is angle between X axis and projection in (X,Y) plan
// Theta is angle between position vector and z axis
- G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
- G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
- G4double wZ = m_R[i] * cos(Theta / rad);
+ G4double wX = m_R[i] * sin(Theta) * cos(Phi);
+ G4double wY = m_R[i] * sin(Theta) * sin(Phi);
+ G4double wZ = m_R[i] * cos(Theta);
MMw = G4ThreeVector(wX, wY, wZ);
// vector corresponding to the center of the module
@@ -406,9 +406,9 @@ void GaspardTrackerTrapezoid::ConstructDetector(G4LogicalVolume* world)
// vector parallel to one axis of silicon plane
// in fact, this is vector u
- G4double ii = cos(Theta / rad) * cos(Phi / rad);
- G4double jj = cos(Theta / rad) * sin(Phi / rad);
- G4double kk = -sin(Theta / rad);
+ G4double ii = cos(Theta) * cos(Phi);
+ G4double jj = cos(Theta) * sin(Phi);
+ G4double kk = -sin(Theta);
G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
MMw = MMw.unit();
diff --git a/NPSimulation/Detectors/Hira/Hira.cc b/NPSimulation/Detectors/Hira/Hira.cc
index eb340f716562db4b3bcf5dffb3e4ecc6ae1ffeca..dd9dd922e5575a23164be69c62db454743cfdd67 100644
--- a/NPSimulation/Detectors/Hira/Hira.cc
+++ b/NPSimulation/Detectors/Hira/Hira.cc
@@ -93,32 +93,32 @@ Hira::~Hira(){
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
void Hira::ReadConfiguration(NPL::InputParser parser ){
- vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("HiraTelescope");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " Telescope found " << endl;
-
- // Cartesian Case
- vector<string> cart = {"A","B","C","D","ThickSi_E","ThinSi_DE","CsI"};
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Hira Telescope " << i+1 << endl;
- G4ThreeVector A = NPS::ConvertVector(blocks[i]->GetTVector3("A","mm"));
- G4ThreeVector B = NPS::ConvertVector(blocks[i]->GetTVector3("B","mm"));
- G4ThreeVector C = NPS::ConvertVector(blocks[i]->GetTVector3("C","mm"));
- G4ThreeVector D = NPS::ConvertVector(blocks[i]->GetTVector3("D","mm"));
- m_build_ThinSi = blocks[i]->GetInt("ThinSi_DE");
- m_build_ThickSi = blocks[i]->GetInt("ThickSi_E");
- m_build_CsI = blocks[i]->GetInt("CsI");
-
- AddTelescope(A,B,C,D) ;
- }
-
- else{
- cout << "ERROR: Missing token for M2Telescope blocks, check your input file" << endl;
- exit(1);
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("HiraTelescope");
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " Telescope found " << endl;
+
+ // Cartesian Case
+ vector<string> cart = {"A","B","C","D","ThickSi_E","ThinSi_DE","CsI"};
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(cart)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Hira Telescope " << i+1 << endl;
+ G4ThreeVector A = NPS::ConvertVector(blocks[i]->GetTVector3("A","mm"));
+ G4ThreeVector B = NPS::ConvertVector(blocks[i]->GetTVector3("B","mm"));
+ G4ThreeVector C = NPS::ConvertVector(blocks[i]->GetTVector3("C","mm"));
+ G4ThreeVector D = NPS::ConvertVector(blocks[i]->GetTVector3("D","mm"));
+ m_build_ThinSi = blocks[i]->GetInt("ThinSi_DE");
+ m_build_ThickSi = blocks[i]->GetInt("ThickSi_E");
+ m_build_CsI = blocks[i]->GetInt("CsI");
+
+ AddTelescope(A,B,C,D) ;
+ }
+
+ else{
+ cout << "ERROR: Missing token for M2Telescope blocks, check your input file" << endl;
+ exit(1);
+ }
}
- }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -290,7 +290,7 @@ void Hira::InitializeScorers(){
void Hira::InitializeRootOutput(){
TTree *pTree = RootOutput::getInstance()->GetTree();
if(!pTree->FindBranch("Hira")){
- pTree->Branch("Hira", "THiraData", &m_EventHira) ;
+ pTree->Branch("Hira", "THiraData", &m_EventHira) ;
}
// This insure that the object are correctly bind in case of
// a redifinition of the geometry in the simulation
@@ -453,7 +453,7 @@ void Hira::VolumeMaker(G4int DetectorNumber,
"CsI_Cristal",
m_LogicCluster,
false,
- CrystalNbr,true);
+ CrystalNbr,false);
delete rotM;
}
}
@@ -462,67 +462,66 @@ void Hira::VolumeMaker(G4int DetectorNumber,
// 3x3 CsI crystal //
/////////////////////
/*double pTheta_1 = -5.6*deg;
- double pPhi_1 = 46.5*deg;
- G4Trap* solidCsI1 = new G4Trap(NameCsI, 0.5*CsIThickness, pTheta_1, pPhi_1, 0.5*23.73, 0.5*23.73, 0.5*23.73, 0, 0.5*32.53, 0.5*32.53, 0.5*32.53, 0);
-
- cout << "Theta1= " << pTheta_1*180/3.1415 << endl;
- cout << "Phi1= " << pPhi_1*180/3.1415 << endl;
-
- double pTheta_2 = -4.0*deg;
- double pPhi_2 = 0*deg;
- G4Trap* solidCsI2 = new G4Trap(NameCsI, 0.5*CsIThickness, pTheta_2, pPhi_2, 0.5*22.55, 0.5*23.73, 0.5*23.73, 0, 0.5*28.05, 0.5*32.53, 0.5*32.53, 0);
-
- G4Trd* solidCsI3 = new G4Trd("SolidCluster", 0.5*22.55,0.5*28.05,0.5*22.55,0.5*28.05, 0.5*CsIThickness);
-
- G4RotationMatrix* rotM1 = new G4RotationMatrix;
- m_LogicCsICrystal = new G4LogicalVolume(solidCsI1, m_MaterialCsI, "logicCsICrystal", 0, 0, 0);
- m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
- m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
- double Xpos = 27;
- double Ypos = 27;
- Pos = G4ThreeVector(-Xpos,-Ypos,0);
- rotM1->rotateZ(0*deg);
- new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal1",m_LogicCluster,false,1,true);
-
- rotM1->rotateZ(-90*deg);
- Pos = G4ThreeVector(-Xpos,Ypos,0);
- new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal2",m_LogicCluster,false,3,true);
-
- rotM1->rotateZ(-90*deg);
- Pos = G4ThreeVector(Xpos,Ypos,0);
- new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal7",m_LogicCluster,false,7,true);
-
- rotM1->rotateZ(-90*deg);
- Pos = G4ThreeVector(Xpos,-Ypos,0);
- new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal9",m_LogicCluster,false,9,true);
-
- G4RotationMatrix* rotM3 = new G4RotationMatrix;
- m_LogicCsICrystal = new G4LogicalVolume(solidCsI3, m_MaterialCsI, "logicCsICrystal", 0, 0, 0);
- m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
- m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
- Pos = G4ThreeVector(0.,0.,0);
- rotM3->rotateZ(0);
- new G4PVPlacement(G4Transform3D(*rotM3,Pos),m_LogicCsICrystal,"CsI_Cristal5",m_LogicCluster,false,5,true);
-
- G4RotationMatrix* rotM2 = new G4RotationMatrix;
- m_LogicCsICrystal = new G4LogicalVolume(solidCsI2, m_MaterialCsI, "logicCsICrystal", 0, 0, 0);
- m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
- m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
- Pos = G4ThreeVector(-Xpos,0.,0);
- rotM2->rotateZ(0*deg);
- new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal2",m_LogicCluster,false,2,true);
-
- Pos = G4ThreeVector(Xpos,0.,0);
- rotM2->rotateZ(180*deg);
- new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal8",m_LogicCluster,false,8,true);
-
- Pos = G4ThreeVector(0,-Ypos,0);
- rotM2->rotateZ(-90*deg);
- new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal4",m_LogicCluster,false,4,true);
-
- Pos = G4ThreeVector(0,Ypos,0);
- rotM2->rotateZ(180*deg);
- new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal6",m_LogicCluster,false,6,true);*/
+ double pPhi_1 = 46.5*deg;
+ G4Trap* solidCsI1 = new G4Trap(NameCsI, 0.5*CsIThickness, pTheta_1, pPhi_1, 0.5*23.73, 0.5*23.73, 0.5*23.73, 0, 0.5*32.53, 0.5*32.53, 0.5*32.53, 0);
+
+ cout << "Theta1= " << pTheta_1*180/3.1415 << endl;
+ cout << "Phi1= " << pPhi_1*180/3.1415 << endl;
+
+ double pTheta_2 = -4.0*deg;
+ double pPhi_2 = 0*deg;
+ G4Trap* solidCsI2 = new G4Trap(NameCsI, 0.5*CsIThickness, pTheta_2, pPhi_2, 0.5*22.55, 0.5*23.73, 0.5*23.73, 0, 0.5*28.05, 0.5*32.53, 0.5*32.53, 0);
+
+ G4Trd* solidCsI3 = new G4Trd("SolidCluster", 0.5*22.55,0.5*28.05,0.5*22.55,0.5*28.05, 0.5*CsIThickness);
+
+ G4RotationMatrix* rotM1 = new G4RotationMatrix;
+ m_LogicCsICrystal = new G4LogicalVolume(solidCsI1, m_MaterialCsI, "logicCsICrystal", 0, 0, 0);
+ m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
+ m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
+ double Xpos = 27;
+ double Ypos = 27;
+ Pos = G4ThreeVector(-Xpos,-Ypos,0);
+ rotM1->rotateZ(0*deg);
+ new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal1",m_LogicCluster,false,1,true);
+
+ rotM1->rotateZ(-90*deg);
+ Pos = G4ThreeVector(-Xpos,Ypos,0);
+ new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal2",m_LogicCluster,false,3,true);
+
+ rotM1->rotateZ(-90*deg);
+ Pos = G4ThreeVector(Xpos,Ypos,0);
+ new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal7",m_LogicCluster,false,7,true);
+
+ rotM1->rotateZ(-90*deg);
+ Pos = G4ThreeVector(Xpos,-Ypos,0);
+ new G4PVPlacement(G4Transform3D(*rotM1,Pos),m_LogicCsICrystal,"CsI_Cristal9",m_LogicCluster,false,9,true);
+
+ G4RotationMatrix* rotM3 = new G4RotationMatrix;
+ m_LogicCsICrystal = new G4LogicalVolume(solidCsI3, m_MaterialCsI, "logicCsICrystal", 0, 0, 0);
+ m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
+ m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
+ Pos = G4ThreeVector(0.,0.,0);
+ rotM3->rotateZ(0);
+ new G4PVPlacement(G4Transform3D(*rotM3,Pos),m_LogicCsICrystal,"CsI_Cristal5",m_LogicCluster,false,5,true);
+
+ G4RotationMatrix* rotM2 = new G4RotationMatrix;
+ m_LogicCsICrystal = new G4LogicalVolume(solidCsI2, m_MaterialCsI, "logicCsICrystal", 0, 0, 0);
+ m_LogicCsICrystal->SetSensitiveDetector(m_CsIScorer);
+ m_LogicCsICrystal->SetVisAttributes(m_CsIVisAtt);
+ Pos = G4ThreeVector(-Xpos,0.,0);
+ rotM2->rotateZ(0*deg);
+ new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal2",m_LogicCluster,false,2,true);
+
+ Pos = G4ThreeVector(Xpos,0.,0);
+ rotM2->rotateZ(180*deg);
+ new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal8",m_LogicCluster,false,8,true);
+ Pos = G4ThreeVector(0,-Ypos,0);
+ rotM2->rotateZ(-90*deg);
+ new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal4",m_LogicCluster,false,4,true);
+
+ Pos = G4ThreeVector(0,Ypos,0);
+ rotM2->rotateZ(180*deg);
+ new G4PVPlacement(G4Transform3D(*rotM2,Pos),m_LogicCsICrystal,"CsI_Cristal6",m_LogicCluster,false,6,true);*/
}
}
@@ -560,5 +559,3 @@ extern"C" {
proxy_nps_hira p_nps_hira;
}
-
-
diff --git a/NPSimulation/Detectors/Hira/Hira.hh b/NPSimulation/Detectors/Hira/Hira.hh
index e2215f7df2851f484b563172b833dadf3056ffdf..bd4ed7f0cb705ad8306d6d88b25e7224e82c1953 100644
--- a/NPSimulation/Detectors/Hira/Hira.hh
+++ b/NPSimulation/Detectors/Hira/Hira.hh
@@ -233,4 +233,4 @@ public:
static NPS::VDetector* Construct();
};
-#endif
+#endif
\ No newline at end of file
diff --git a/NPSimulation/Process/G4DEAmplification.cc b/NPSimulation/Process/G4DEAmplification.cc
index a171c318799cfb361f50751e682361000d1d9044..8e8b3e272e50f841873ba81d9c2944666840d9ff 100644
--- a/NPSimulation/Process/G4DEAmplification.cc
+++ b/NPSimulation/Process/G4DEAmplification.cc
@@ -97,19 +97,17 @@ G4DEAmplification::~G4DEAmplification(){}
G4VParticleChange*
G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
{
+
+// if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" )
+// return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
+
// Get the primary track
aParticleChange.Initialize(aTrack);
- // Check that the parent process is not amplification
- if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" )
- return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
-
const G4Material* aMaterial = aTrack.GetMaterial();
G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
-
-
G4ThreeVector x0 = pPreStepPoint->GetPosition();
G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
G4double t0 = pPreStepPoint->GetGlobalTime();
@@ -121,24 +119,14 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
- if(!aMaterialPropertiesTable->ConstPropertyExists("DE_PAIRENERGY") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_AMPLIFICATION") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_DRIFTSPEED") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_TRANSVERSALSPREAD") ||
- !aMaterialPropertiesTable->ConstPropertyExists("DE_LONGITUDINALSPREAD"))
+ if(!aMaterialPropertiesTable->ConstPropertyExists("DE_AMPLIFICATION") ||
+ !aMaterialPropertiesTable->ConstPropertyExists("DE_YIELD"))
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
- G4double pair_energy = 0;
- pair_energy =
- aMaterialPropertiesTable->GetConstProperty("DE_PAIRENERGY");
+ G4double pair_energy = pPreStepPoint->GetKineticEnergy();
- G4double amplification=0;
- amplification=
- aMaterialPropertiesTable->GetConstProperty("DE_AMPLIFICATION");
- G4double Yield = 0;
- Yield=
- aMaterialPropertiesTable->GetConstProperty("DE_YIELD");
+ G4double amplification = aMaterialPropertiesTable->GetConstProperty("DE_AMPLIFICATION");
+ G4double Yield = aMaterialPropertiesTable->GetConstProperty("DE_YIELD");
G4int number_electron = amplification*Yield;
//if no electron leave
@@ -148,7 +136,6 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
}
aParticleChange.SetNumberOfSecondaries(number_electron);
-
// Create the secondary tracks
for(G4int i = 0 ; i < number_electron ; i++){
// Random direction at creation
@@ -158,42 +145,48 @@ G4DEAmplification::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4ThreeVector p;
p.setRThetaPhi(1,theta,phi);
- // Random Position along the step with matching time
- G4double rand = G4UniformRand();
- G4ThreeVector pos = x0 + rand * aStep.GetDeltaPosition();
- G4double time = t0+ rand* aStep.GetDeltaTime();
-
G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),p, pair_energy);
- G4Track* aSecondaryTrack = new G4Track(particle,time,pos);
+ G4Track* aSecondaryTrack = new G4Track(particle,t0,x0);
aSecondaryTrack->SetTouchableHandle(
aStep.GetPreStepPoint()->GetTouchableHandle());
-
+ aSecondaryTrack->SetCreatorProcess(this);
aSecondaryTrack->SetParentID(aTrack.GetTrackID());
aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
aParticleChange.AddSecondary(aSecondaryTrack);
}
-
- return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
+ // to kill the primary track to avoid it being re-amplified
+ aParticleChange.ProposeEnergy(0);
+ // The original track is left intact
+ return &aParticleChange;
}
// GetMeanFreePath
// ---------------
//
-G4double G4DEAmplification::GetMeanFreePath(const G4Track& ,
+G4double G4DEAmplification::GetMeanFreePath(const G4Track& aTrack,
G4double ,
G4ForceCondition* condition)
{
+ if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" ){
+ *condition = NotForced;
+ return DBL_MAX;
+ }
*condition = StronglyForced;
return DBL_MAX;
}
// GetMeanLifeTime
// ---------------
//
-G4double G4DEAmplification::GetMeanLifeTime(const G4Track& ,
+G4double G4DEAmplification::GetMeanLifeTime(const G4Track& aTrack,
G4ForceCondition* condition)
{
+ if(aTrack.GetCreatorProcess()->GetProcessName()=="DEAmplification" ){
+ *condition = NotForced;
+ return DBL_MAX;
+ }
+
*condition = StronglyForced;
return DBL_MAX;
}
diff --git a/NPSimulation/Process/G4DETransport.cc b/NPSimulation/Process/G4DETransport.cc
index 87578f5d8ec9a1fccf4dd6fec8b4f606c1e1544c..aa7a75fbcc7b230761583c9672dc93b44866a29e 100644
--- a/NPSimulation/Process/G4DETransport.cc
+++ b/NPSimulation/Process/G4DETransport.cc
@@ -49,6 +49,8 @@
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4DETransport.hh"
+#include "G4FieldManager.hh"
+#include "G4ElectroMagneticField.hh"
/////////////////////////
// Class Implementation
@@ -74,6 +76,8 @@ G4DETransport::G4DETransport(const G4String& processName, G4ProcessType type)
}
SetProcessSubType(fDETransport);
+ G4TransportationManager* TransportationManager = G4TransportationManager::GetTransportationManager();
+ m_SafetyHelper = TransportationManager->GetSafetyHelper();
}
// G4DETransport::G4DETransport(const G4OpAbsorpton &right)
@@ -88,7 +92,6 @@ G4DETransport::~G4DETransport(){}
////////////
// Methods
-////////////
// PostStepDoIt
// -------------
@@ -97,6 +100,7 @@ G4VParticleChange*
G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
{
aParticleChange.Initialize(aTrack);
+
G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
G4ThreeVector x0 = pPreStepPoint->GetPosition();
@@ -105,6 +109,15 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4double Energy = pPreStepPoint->GetKineticEnergy();
+ // The time scale is imposed by the distance travelled
+ G4double step_length = aStep.GetDeltaPosition().mag();
+
+ // allow internal relocation of the track by the kernel taking a 0 length intermediate step
+ // suppress also parasite infinetismal step that slow down the tracking
+ if(step_length<100*micrometer){
+ return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
+ }
+
// Get the material table
const G4Material* aMaterial = aTrack.GetMaterial();
@@ -123,33 +136,55 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
// Electron follow the field direction
- G4ThreeVector driftDir(0,1,0);
+ // The field direction is taken from the field manager
+ G4double* fieldArr = new G4double[6];
+ G4double Point[4]={x0.x(),x0.y(),x0.z(),t0};
+ G4FieldManager* fMng = pPreStepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume()->
+ GetFieldManager();
+
+ G4ElectroMagneticField* field = (G4ElectroMagneticField*)fMng->GetDetectorField();
+ field->GetFieldValue(Point,fieldArr) ;
+
+ // Electron move opposite to the field direction, hance the minus sign
+ G4ThreeVector driftDir(-fieldArr[3],-fieldArr[4],-fieldArr[5]);
+ // Normalised the drift direction
+ driftDir = driftDir.unit();
G4double v_drift = aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED");
- // The time scale is imposed by the drift speed
- G4double step = aStep.GetDeltaTime();
- if(!step)
- return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
- // Should be equal to delta pos y
- G4double step_length = step*v_drift;
G4double sigmaTrans = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift);
G4double sigmaLong = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift);
G4double d_long = G4RandGauss::shoot(0,sigmaLong);
G4double d_trans = G4RandGauss::shoot(0,sigmaTrans);
-
+ G4double d_drift = step_length+d_long;
+
G4ThreeVector trans(G4RandGauss::shoot(0,d_trans),0,0);
trans.rotateY(twopi*G4UniformRand());
- G4ThreeVector d_Pos = (step_length+d_long)*driftDir+trans;
- // Random Position along the step with matching time
+ G4ThreeVector d_Pos = (d_drift)*driftDir+trans;
+
+ // Should be equal to delta length
+ G4double step = (d_drift)/v_drift;
+
+ // New particle Position with matching time
G4ThreeVector pos = x0 + d_Pos;
G4double time = t0 + step;
+
+ // Garanty that the electron does not jump outside the current volume
+ G4double safety = m_SafetyHelper->ComputeSafety(x0,d_Pos.mag());
+
+ // If the distance travelled if above safety, the step is not taken
+ if(d_Pos.mag()>safety){
+ return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
+ }
aParticleChange.ProposeMomentumDirection(d_Pos.unit());
aParticleChange.ProposeEnergy(Energy);
aParticleChange.ProposePosition(pos);
aParticleChange.ProposeGlobalTime(time);
+ aParticleChange.ProposeLocalTime(time);
+ aParticleChange.ProposeVelocity(v_drift/c_light);
+ m_SafetyHelper->ReLocateWithinVolume(pos);
return &aParticleChange;
}
@@ -161,5 +196,8 @@ G4double G4DETransport::GetMeanFreePath(const G4Track& ,
G4double ,
G4ForceCondition* )
{
- return 0.5*mm;
+ // Typical distance after which the electron position should be reevaluated
+ // to take into account the diffusivity of the charge cloud inside the
+ // medium
+ return 1*mm;
}
diff --git a/NPSimulation/Process/G4DETransport.hh b/NPSimulation/Process/G4DETransport.hh
index 6505959f1168e290aa501396888450cd5a2a6b16..d885efe8bb961e53ff986337814a062fd3a04599 100644
--- a/NPSimulation/Process/G4DETransport.hh
+++ b/NPSimulation/Process/G4DETransport.hh
@@ -59,7 +59,8 @@
#include "G4DynamicParticle.hh"
#include "G4Material.hh"
#include "G4DriftElectron.hh"
-
+#include "G4TransportationManager.hh"
+class G4SafetyHelper;
// Class Description:
// Discrete Process -- Bulk transport of Drift Electron.
// Class inherits publicly from G4VDiscreteProcess
@@ -104,14 +105,16 @@ public:
G4double GetMeanFreePath(const G4Track& aTrack,
G4double ,
G4ForceCondition* );
- // Returns the transport length for bulk transport of drift
- // electron in media with a specified attenuation length.
+ // Returns the typical travel length for transport of drift
+ // electron in media
G4VParticleChange* PostStepDoIt(const G4Track& aTrack,
const G4Step& aStep);
// This is the method implementing bulk transport of drift
// electron.
+private:
+ G4SafetyHelper* m_SafetyHelper;
};
////////////////////
diff --git a/NPSimulation/Process/G4DriftElectron.cc b/NPSimulation/Process/G4DriftElectron.cc
index 04a09993aa1671b95de427f945f2272758cd75f0..daadc4d2d54ae3bee4248b2881951dfcc6085f81 100644
--- a/NPSimulation/Process/G4DriftElectron.cc
+++ b/NPSimulation/Process/G4DriftElectron.cc
@@ -53,6 +53,9 @@ G4DriftElectron* G4DriftElectron::Definition()
// search in particle table]
G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* anInstance = pTable->FindParticle(name);
+ // Note: The charge of drift electron is null so they are not afected by the
+ // electro magnetic field directly. This is taken care of "manually" by the
+ // G4DETransport process that get the local field direction
if (anInstance ==0)
{
// create particle
@@ -65,14 +68,15 @@ G4DriftElectron* G4DriftElectron::Definition()
// stable lifetime decay table
// shortlived subType anti_encoding
anInstance = new G4ParticleDefinition(
- name, electron_mass_c2, 0.0*MeV, -1.*eplus,
+ name, electron_mass_c2, 0.0*MeV, 0,
1, 0, 0,
0, 0, 0,
"lepton", 1, 0, 11,
true, -1.0, NULL,
- false, "e"
+ false, "driftion"
);
}
+
theInstance = reinterpret_cast<G4DriftElectron*>(anInstance);
return theInstance;
}
diff --git a/NPSimulation/Process/G4DriftElectronProcessIndex.hh b/NPSimulation/Process/G4DriftElectronProcessIndex.hh
index 9593a740b382df3dbf6c7cac5d2a45907b816b80..b82ce051fdcb0e430f322cd889884371d4f9c499 100644
--- a/NPSimulation/Process/G4DriftElectronProcessIndex.hh
+++ b/NPSimulation/Process/G4DriftElectronProcessIndex.hh
@@ -66,10 +66,10 @@ inline
G4String G4DriftElectronProcessName(G4int processNumber)
{
switch ( processNumber ) {
- case kDEIonization: return "Ionization";
- case kDEAmplification: return "Amplification";
- case kDEAbsorption: return "Absorption";
- default: return "NoProcess";
+ case kDEIonization: return "DEIonization";
+ case kDEAmplification: return "DEAmplification";
+ case kDEAbsorption: return "DEAbsorption";
+ default: return "DENoProcess";
}
}
diff --git a/NPSimulation/Process/G4IonizationWithDE.cc b/NPSimulation/Process/G4IonizationWithDE.cc
index 2856283b4ad1de21a22170a4a407f7344feef3f3..b4a2abe6cd2b0c4edbfb4ad66606c9835ae2275b 100644
--- a/NPSimulation/Process/G4IonizationWithDE.cc
+++ b/NPSimulation/Process/G4IonizationWithDE.cc
@@ -49,7 +49,8 @@
#include "G4SystemOfUnits.hh"
#include "G4ParticleTypes.hh"
#include "G4EmProcessSubType.hh"
-
+#include "G4ElectroMagneticField.hh"
+#include "G4FieldManager.hh"
#include "G4IonizationWithDE.hh"
#include <iostream>
using namespace std;
@@ -163,14 +164,16 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
- G4double pair_energy=0;
- pair_energy=
+ G4double v_drift=
+ aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED");
+ G4double pair_energy=
aMaterialPropertiesTable->GetConstProperty("DE_PAIRENERGY");
G4double IonizationWithDEYield = 0;
IonizationWithDEYield=
aMaterialPropertiesTable->GetConstProperty("DE_YIELD");
G4int number_electron = IonizationWithDEYield*TotalEnergyDeposit/pair_energy;
+ number_electron = G4Poisson(number_electron);
//if no electron leave
if(number_electron<1){
aParticleChange.SetNumberOfSecondaries(0);
@@ -181,12 +184,20 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
// Create the secondary tracks
for(G4int i = 0 ; i < number_electron ; i++){
- // Random direction at creation
- G4double cost = 1-2*G4UniformRand();
- G4double theta = acos(cost);
- G4double phi = twopi*G4UniformRand();
- G4ThreeVector p;
- p.setRThetaPhi(1,theta,phi);
+ // Electron follow the field direction
+ // The field direction is taken from the field manager
+ G4double* fieldArr = new G4double[6];
+ G4double Point[4]={x0.x(),x0.y(),x0.z(),t0};
+ G4FieldManager* fMng = pPreStepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume()->
+ GetFieldManager();
+
+ G4ElectroMagneticField* field = (G4ElectroMagneticField*)fMng->GetDetectorField();
+ field->GetFieldValue(Point,fieldArr) ;
+
+ // Electron move opposite to the field direction, hance the minus sign
+ G4ThreeVector p(-fieldArr[3],-fieldArr[4],-fieldArr[5]);
+ // Normalised the drift direction
+ p = p.unit();
// Random Position along the step with matching time
G4double rand = G4UniformRand();
@@ -195,9 +206,7 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
G4DynamicParticle* particle = new G4DynamicParticle(G4DriftElectron::DriftElectron(),p, pair_energy);
G4Track* aSecondaryTrack = new G4Track(particle,time,pos);
-
- aSecondaryTrack->SetTouchableHandle(
- aStep.GetPreStepPoint()->GetTouchableHandle());
+ aSecondaryTrack->SetVelocity(v_drift/c_light);
aSecondaryTrack->SetParentID(aTrack.GetTrackID());
aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
diff --git a/NPSimulation/Process/PhysicsList.cc b/NPSimulation/Process/PhysicsList.cc
index 6cc59d9a53f4eecd52b04bb0ab50851e4bda9ee2..cf08ee4881d969ba626ae8c5bbe1c2e33d2e191d 100644
--- a/NPSimulation/Process/PhysicsList.cc
+++ b/NPSimulation/Process/PhysicsList.cc
@@ -42,7 +42,7 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList(){
m_EmList = "Option4";
defaultCutValue = 1*mm;//0.2*mm;
opticalPhysicsList = NULL;
-
+ driftElectronPhysicsList = NULL;
ReadConfiguration("PhysicsListOption.txt");
G4LossTableManager::Instance();
SetVerboseLevel(0);
@@ -116,8 +116,8 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList(){
opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true);
}
+ // Drift electron for gazeous detector simulation
if(m_DriftElectronPhysics){
-
driftElectronPhysicsList = new G4DriftElectronPhysics(0);
driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6);
}
@@ -213,7 +213,7 @@ void PhysicsList::ConstructParticle(){
}
- if(driftElectronPhysicsList){
+ if(m_DriftElectronPhysics){
((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructParticle();
}
@@ -291,7 +291,7 @@ void PhysicsList::AddStepMax(){
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddParametrisation(){
-
+/*
G4FastSimulationManagerProcess* drift =
new G4FastSimulationManagerProcess("DriftElectron");
@@ -307,7 +307,7 @@ void PhysicsList::AddParametrisation(){
if(particle->GetParticleName()=="e-")
pmanager->AddDiscreteProcess(drift);
- }
+ }*/
}
diff --git a/NPSimulation/Scorers/CMakeLists.txt b/NPSimulation/Scorers/CMakeLists.txt
index ed220c926f9a921baa5df8d650afd3fb5ad7ab5a..f54c2bdbf5d8e506ce9ab80ed3998e087a52f42a 100644
--- a/NPSimulation/Scorers/CMakeLists.txt
+++ b/NPSimulation/Scorers/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSScorers SHARED NPSHitsMap.hh CalorimeterScorers.cc SiliconScorers.cc PhotoDiodeScorers.cc ObsoleteGeneralScorers.cc)
+add_library(NPSScorers SHARED NPSHitsMap.hh CalorimeterScorers.cc SiliconScorers.cc PhotoDiodeScorers.cc ObsoleteGeneralScorers.cc DriftElectronScorers.cc )
target_link_libraries(NPSScorers ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates)
diff --git a/NPSimulation/Scorers/DriftElectronScorers.cc b/NPSimulation/Scorers/DriftElectronScorers.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4c0548f425673d72e9cd3bbfe81496dfd4905aa3
--- /dev/null
+++ b/NPSimulation/Scorers/DriftElectronScorers.cc
@@ -0,0 +1,198 @@
+ /*****************************************************************************
+ * Copyright (C) 2009-2016 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: morfouac@nscl.msu.edu *
+ * *
+ * Creation Date : January 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Scorer for Drift Electron collector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * This new type of scorer to count drift electron in the cathode *
+ * *
+ *****************************************************************************/
+#include "DriftElectronScorers.hh"
+#include "G4UnitsTable.hh"
+using namespace DRIFTELECTRONSCORERS ;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_DECathode::PS_DECathode(G4String name,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+ m_Index = -1 ;
+ m_Level = depth;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_DECathode::~PS_DECathode(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_DECathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+
+ // contain Energy Time, DetNbr, StripFront and StripBack
+ G4double* Infos = new G4double[9];
+ Infos[0] = 0;
+ Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
+ m_Position = aStep->GetPreStepPoint()->GetPosition();
+
+ // Interaction coordinates (used to fill the InteractionCoordinates branch)
+ Infos[2] = m_Position.x();
+ Infos[3] = m_Position.y();
+ Infos[4] = m_Position.z();
+ Infos[5] = m_Position.theta();
+ Infos[6] = m_Position.phi();
+ Infos[7] = m_DetectorNumber;
+
+ m_Index = m_DetectorNumber * 1e3 ;
+ G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName();
+
+ if(PID=="driftelectron"){
+ Infos[0] = 1;
+ }
+
+ // Check if the particle has interact before, if yes, add up the number of electron.
+ map<G4int, G4double**>::iterator it;
+ it= EvtMap->GetMap()->find(m_Index);
+ if(it!=EvtMap->GetMap()->end()){
+ G4double* dummy = *(it->second);
+ Infos[0]+=dummy[0];
+ Infos[1]=dummy[1];
+ }
+
+ EvtMap->set(m_Index, Infos);
+ return TRUE;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DECathode::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DECathode::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DECathode::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DECathode::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DECathode::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_DEDigitizer::PS_DEDigitizer(G4String name,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+ m_Index = -1 ;
+ m_Level = depth;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_DEDigitizer::~PS_DEDigitizer(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_DEDigitizer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+
+ // contain Energy Time, DetNbr, StripFront and StripBack
+ G4double* Infos = new G4double[9];
+ Infos[0] = 0;
+ Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
+ m_Position = aStep->GetPreStepPoint()->GetPosition();
+
+ // Interaction coordinates (used to fill the InteractionCoordinates branch)
+ Infos[2] = m_Position.x();
+ Infos[3] = m_Position.y();
+ Infos[4] = m_Position.z();
+ Infos[5] = m_Position.theta();
+ Infos[6] = m_Position.phi();
+ Infos[7] = m_DetectorNumber;
+
+ m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber*1e6 ;
+ G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName();
+
+ if(PID=="driftelectron"){
+ Infos[0] = 1;
+ }
+
+ // Check if the particle has interact before, if yes, add up the number of electron.
+ map<G4int, G4double**>::iterator it;
+ it= EvtMap->GetMap()->find(m_Index);
+ if(it!=EvtMap->GetMap()->end()){
+ G4double* dummy = *(it->second);
+ Infos[0]+=dummy[0];
+ Infos[1]=dummy[1];
+ }
+
+ EvtMap->set(m_Index, Infos);
+ return TRUE;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
diff --git a/NPSimulation/Scorers/DriftElectronScorers.hh b/NPSimulation/Scorers/DriftElectronScorers.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6b49c8a67eb960252bd1ab8dd05d4914bf6025e0
--- /dev/null
+++ b/NPSimulation/Scorers/DriftElectronScorers.hh
@@ -0,0 +1,95 @@
+#ifndef DECathodeScorers_h
+#define SiliconScorers_h 1
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : February 2013 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * File old the scorer specific to the Silicon Detector *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * This new style of scorer is aim to become the standard way of doing scorer*
+ * in NPTool. *
+ *The index is build using the TrackID, Detector Number and Strip Number. *
+ *The scorer Hold Energy and time together *
+ *Only one scorer is needed for a detector *
+ *****************************************************************************/
+#include "G4VPrimitiveScorer.hh"
+#include "NPSHitsMap.hh"
+
+#include <map>
+using namespace std;
+using namespace CLHEP;
+
+namespace DRIFTELECTRONSCORERS {
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ class PS_DECathode : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_DECathode(G4String name,G4int depth=0);
+ ~PS_DECathode();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ NPS::HitsMap<G4double*>* EvtMap;
+
+ private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
+ G4int m_Index;
+ G4int m_Level;
+ G4int m_DetectorNumber;
+ G4ThreeVector m_Position;
+
+ };
+////////////////////////////////////////////////////////////////////////////////
+ class PS_DEDigitizer : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_DEDigitizer(G4String name,G4int depth=0);
+ ~PS_DEDigitizer();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ NPS::HitsMap<G4double*>* EvtMap;
+
+ private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
+ G4int m_Index;
+ G4int m_Level;
+ G4int m_DetectorNumber;
+ G4ThreeVector m_Position;
+
+ };
+
+
+}
+
+#endif
diff --git a/Projects/BeLise/9Li.reaction b/Projects/BeLise/9Li.reaction
index 02db44e328bc5c5c8ab41d514f4d8ec1ae645b49..8e29608a76dc505d47dbe5903c3eec6d4e211967 100644
--- a/Projects/BeLise/9Li.reaction
+++ b/Projects/BeLise/9Li.reaction
@@ -7,12 +7,12 @@ Beam
SigmaEnergy= 20 MeV
SigmaThetaX= 0.6138 deg
SigmaPhiY= 0.3812 deg
- SigmaX= 5.1 mm
- SigmaY= 8.5 mm
+ SigmaX= 6.216 mm
+ SigmaY= 6.109 mm
MeanThetaX= 0 deg
MeanPhiY= 0 deg
- MeanX= 2.5 mm
- MeanY= -0.5 mm
+ MeanX= 0 mm
+ MeanY= 0 mm
%EnergyProfilePath=
%XThetaXProfilePath=
%YPhiYProfilePath=
@@ -24,8 +24,8 @@ TwoBodyReaction
Light= 3He
Heavy= 9Li
ExcitationEnergyLight= 0.0 MeV
- ExcitationEnergyHeavy= 0 MeV
- CrossSectionPath= 11Li(d,3He)10He.txt CS10He
+ ExcitationEnergyHeavy= 0.0 MeV
+ CrossSectionPath= 11Li(d,3He)10He.txt CS2
ShootLight= 1
ShootHeavy= 1
diff --git a/Projects/BeLise/Analysis.cxx b/Projects/BeLise/Analysis.cxx
index c8a8a07a7219eeac42150b405990af3e39e26106..5d0a7676115caa8266da6fbd60a5e092c34d0ac2 100644
--- a/Projects/BeLise/Analysis.cxx
+++ b/Projects/BeLise/Analysis.cxx
@@ -78,7 +78,8 @@ void Analysis::TreatEvent(){
// Beam energy is measured using F3 and F2 plastic TOF
double BeamEnergy = Rand.Gaus(Initial->GetIncidentInitialKineticEnergy(),4.5);
//BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
-
+ OriginalThetaLab = Initial->GetThetaCM(0);
+ OriginalELab = Initial->GetKineticEnergy(0);
myReaction->SetBeamEnergy(BeamEnergy);
//////////////////////////// LOOP on MUST2 Hit //////////////////
for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
@@ -174,6 +175,10 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalELab",&OriginalELab,"OriginalELab/D");
+
+
}
////////////////////////////////////////////////////////////////////////////////
@@ -189,6 +194,8 @@ void Analysis::ReInitValue(){
ELab = -1000;
ThetaLab = -1000;
ThetaCM = -1000;
+ OriginalThetaLab = -1000;
+ OriginalELab = -1000;
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/BeLise/Analysis.h b/Projects/BeLise/Analysis.h
index 1e546814cbe2a02b632230569772c525fdfe321d..8e3cc283c4cbf8bb1ced7d115b1db8bef10af010 100644
--- a/Projects/BeLise/Analysis.h
+++ b/Projects/BeLise/Analysis.h
@@ -64,7 +64,8 @@ class Analysis: public NPL::VAnalysis{
double Si_Y_M2;
double ZTarget;
double TargetThickness;
-
+ double OriginalThetaLab;
+ double OriginalELab;
NPL::EnergyLoss He3CD2 ;
NPL::EnergyLoss He3Al ;
NPL::EnergyLoss He3Si ;
diff --git a/Projects/BeLise/PhysicsListOption.txt b/Projects/BeLise/PhysicsListOption.txt
index d85f8d5dbca0295f25f81da9b192609ea6bf0523..8c522e9f8b5c39436811edde1a5f9eb64fc2da75 100644
--- a/Projects/BeLise/PhysicsListOption.txt
+++ b/Projects/BeLise/PhysicsListOption.txt
@@ -1,5 +1,5 @@
EmPhysicsList Option4
-DefaultCutOff 100
+DefaultCutOff 1000
IonBinaryCascadePhysics 0
NPIonInelasticPhysics 0
EmExtraPhysics 0
diff --git a/Projects/BeLise/chio.detector b/Projects/BeLise/chio.detector
index 0166a75518a63e8321c785c8c156a7c604d2755c..7c249a657b9cb7893bea3471564c537c5067d0bc 100644
--- a/Projects/BeLise/chio.detector
+++ b/Projects/BeLise/chio.detector
@@ -9,6 +9,7 @@ Target
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chio
- POS= 0 0 350 mm
- Shape= Square
+ POS= 0 0 400 mm
+ Shape= Square
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/Hira10/Analysis.cxx b/Projects/Hira10/Analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d0ab7b29fc57267db57bcf5bd5b1878c3438210b
--- /dev/null
+++ b/Projects/Hira10/Analysis.cxx
@@ -0,0 +1,335 @@
+ /*****************************************************************************
+ * Copyright (C) 2009-2014 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: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include<iostream>
+using namespace std;
+#include"Analysis.h"
+#include"NPAnalysisFactory.h"
+#include"NPDetectorManager.h"
+#include"NPOptionManager.h"
+#include"RootOutput.h"
+#include"RootInput.h"
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis(){
+}
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init(){
+ Hira = (THiraPhysics*) m_DetectorManager->GetDetector("HIRAArray");
+ InitialConditions=new TInitialConditions();
+ InitOutputBranch();
+ InitInputBranch();
+ Rand = TRandom3();
+ //TransferReaction= new NPL::Reaction("8B(p,4He)5Be@400");
+ //TransferReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+
+ TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
+ // Energy loss table: the G4Table are generated by the simulation
+ Triton_CD2 = EnergyLoss("triton_CD2.G4table","G4Table",100 );
+ Triton_CH2 = EnergyLoss("triton_CD2.G4table","G4Table",100 );
+ Deuteron_CH2 = EnergyLoss("deuteron_CH2.G4table","G4Table",100 );
+ He3_CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
+ He4_CH2 = EnergyLoss("alpha_CH2.G4table","G4Table",100 );
+
+ f_proton = new TF1("f_proton","1 - TMath::Exp(-(x-0.0918309)/27.7746)",0,10);
+ f_deuton = new TF1("f_deuton","1 - TMath::Exp(-(x-0.0434552)/21.134)",0,10);
+ f_triton = new TF1("f_triton","1 - TMath::Exp(-(x-0.0353106)/17.9059)",0,10);
+ f_3He = new TF1("f_3He","1 - TMath::Exp(-(x-0.0463954)/20.8906)",0,10);
+
+ // Energy loss table: the G4Table are generated by the simulation
+ Proton_CsI = EnergyLoss("proton_CsI.G4table","G4Table",100 );
+ Deuton_CsI = EnergyLoss("deuteron_CsI.G4table","G4Table",100 );
+ Triton_CsI = EnergyLoss("triton_CsI.G4table","G4Table",100 );
+ He3_CsI = EnergyLoss("He3_CsI.G4table","G4Table",100 );
+
+ proton = new NPL::Nucleus("1H");
+ deuton = new NPL::Nucleus("2H");
+ triton = new NPL::Nucleus("3H");
+ helium3 = new NPL::Nucleus("3He");
+ alpha = new NPL::Nucleus("4He");
+ beam = new NPL::Nucleus("48Ca");
+ target = new NPL::Nucleus("124Sn");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent(){
+ // Reinitiate calculated variable
+ ReInitValue();
+
+ double BeamEnergy = 120*beam->GetA();
+
+ TVector3 fImpulsionLab_beam = TVector3(0,0,sqrt(BeamEnergy*BeamEnergy + 2*BeamEnergy*beam->Mass()));
+ TLorentzVector fEnergyImpulsionLab_beam = TLorentzVector(fImpulsionLab_beam,beam->Mass()+BeamEnergy);
+
+ TVector3 fImpulsionLab_target = TVector3(0,0,0);
+ TLorentzVector fEnergyImpulsionLab_target = TLorentzVector(fImpulsionLab_target,target->Mass());
+
+ TLorentzVector fEnergyImpulsionLab_total;
+ fEnergyImpulsionLab_total = fEnergyImpulsionLab_beam + fEnergyImpulsionLab_target;
+
+ double BetaCM = fEnergyImpulsionLab_total.Beta();
+
+ InitialEnergy = InitialConditions->GetKineticEnergy(0);
+ double EDelta = 2.0;
+
+ if(Hira->EventMultiplicity==1) InitialEnergy_Hira = InitialEnergy;
+ else InitialEnergy_Hira = -100;
+ double phi_in = acos(InitialConditions->GetMomentumDirectionX(0)/sin(InitialConditions->GetThetaCM(0)*deg));
+
+ ECM_initial = proton->GetEnergyCM(InitialEnergy, InitialConditions->GetThetaCM(0)*deg, phi_in, BetaCM);
+ if(Hira->ThickSi_E.size()==1){
+ ECM_initial_Hira = proton->GetEnergyCM(InitialEnergy_Hira, InitialConditions->GetThetaCM(0)*deg, phi_in, BetaCM);
+ }
+ else ECM_initial_Hira = -100;
+ //TransferReaction->SetBeamEnergy(BeamEnergy);
+ //////////////////////////// LOOP on Hira Hit //////////////////
+ if(Hira -> EventMultiplicity == 1){
+ for(unsigned int countHira = 0 ; countHira < Hira->EventMultiplicity ; countHira++){
+ event += 1;
+ TelescopeNumber = Hira->TelescopeNumber[countHira];
+
+ TargetThickness = m_DetectorManager->GetTargetThickness();
+
+ X = Hira->GetPositionOfInteraction(countHira).X();
+ Y = Hira->GetPositionOfInteraction(countHira).Y();
+ Z = Hira->GetPositionOfInteraction(countHira).Z();
+
+ TVector3 PositionOnHira = TVector3(X,Y,Z);
+ TVector3 ZUnit = TVector3(0,0,1);
+
+ X_target = InitialConditions->GetIncidentPositionX();// + Rand.Gaus(0,2.);
+ Y_target = InitialConditions->GetIncidentPositionY();// + Rand.Gaus(0,2.);
+ double Z_target = InitialConditions->GetIncidentPositionZ();
+
+ //TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
+ TVector3 PositionOnTarget = TVector3(0,0,0);
+ TVector3 HitDirection = PositionOnHira - PositionOnTarget;
+ TVector3 HitDirectionUnit = HitDirection.Unit();
+ TVector3 BeamDirection = TVector3(0,0,1);
+
+ ThetaLab = BeamDirection.Angle(HitDirection);
+ PhiLab = HitDirection.Phi();
+
+
+ E_ThickSi = Hira->ThickSi_E[countHira];
+ //E_ThinSi = Hira->ThinSi_E[countHira];
+
+ ELab = E_ThickSi;
+ if(Hira->CsI_E[countHira]>0){
+ E_CsI = Hira->CsI_E[countHira];
+ ELab += E_CsI;
+
+ PID = pow(E_ThickSi+E_CsI,1.78)-pow(E_CsI,1.78);
+
+ //Try to simulate the nuclear reaction loss
+ //ThicknessCsI = Proton_CsI.EvaluateMaterialThickness(0*MeV, Hira->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+ //ThicknessCsI = Deuton_CsI.EvaluateMaterialThickness(0*MeV, Hira->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+ //ThicknessCsI = Triton_CsI.EvaluateMaterialThickness(0*MeV, Hira->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+
+
+ //double eval = f_proton->Eval(ThicknessCsI/10);
+ //double eval = f_deuton->Eval(ThicknessCsI/10);
+ //double eval = f_triton->Eval(ThicknessCsI/10);
+ //double Random_value = Rand.Uniform(0,1);
+
+ }
+
+ if(fabs(InitialEnergy-ELab)>EDelta){
+ ELab = -100;
+ }
+
+
+ ////////////////////////////////////////////////////////////////////////
+ // Calculation of the center of mass energy depending on the particle //
+ ////////////////////////////////////////////////////////////////////////
+ ECM = proton->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ThetaCM = proton->GetThetaCM(ELab, ThetaLab, PhiLab, BetaCM)/deg;
+ /*if(PID>100 && PID<175){
+ ECM = proton->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ThetaCM = proton->GetThetaCM(ELab, ThetaLab, PhiLab, BetaCM)/deg;
+ Particle = 1;
+ }
+ else if(PID>175 && PID<265){
+ ECM = deuton->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ThetaCM = deuton->GetThetaCM(ELab, ThetaLab, PhiLab, BetaCM)/deg;
+ Particle = 2;
+ }
+ else if(PID>265 && PID<350){
+ ECM = triton->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ThetaCM = triton->GetThetaCM(ELab, ThetaLab, PhiLab, BetaCM)/deg;
+ Particle = 3;
+ }
+ else if(PID>1100 && PID<1370){
+ ECM = helium3->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ThetaCM = helium3->GetThetaCM(ELab, ThetaLab, PhiLab, BetaCM)/deg;
+ Particle = 4;
+ }
+ else if(PID>1370 && PID<1700){
+ ECM = alpha->GetEnergyCM(ELab, ThetaLab, PhiLab, BetaCM);
+ ThetaCM = alpha->GetThetaCM(ELab, ThetaLab, PhiLab, BetaCM)/deg;
+ Particle = 5;
+ }
+ else{
+ ECM = -100;
+ ThetaCM = -100;
+ Particle = -1;
+ }*/
+ ////////////////////////////////////////////////////////////////////////
+
+ double ImpulsionLab = sqrt(ELab*ELab + 2*ELab*proton->Mass());
+ double ImpulsionLabX = ImpulsionLab*sin(ThetaLab)*cos(PhiLab);
+ double ImpulsionLabY = ImpulsionLab*sin(ThetaLab)*sin(PhiLab);
+ double ImpulsionLabZ = ImpulsionLab*cos(ThetaLab);
+
+ TVector3 VImpulsionLAB = TVector3(ImpulsionLabX, ImpulsionLabY, ImpulsionLabZ);
+ TLorentzVector LVEnergyImpulsionLAB = TLorentzVector(VImpulsionLAB,proton->Mass()+ELab);
+
+ TLorentzVector LVEnergyImpulsionCM = LVEnergyImpulsionLAB;
+ LVEnergyImpulsionCM.Boost(0,0,-BetaCM);
+
+ Rapidity_Lab = LVEnergyImpulsionLAB.Rapidity();
+ Rapidity_CM = LVEnergyImpulsionCM.Rapidity();
+ Pper = LVEnergyImpulsionCM.Pt();
+
+ /*double Px_Lab = LVEnergyImpulsionLAB.Px();
+ double Py_Lab = LVEnergyImpulsionLAB.Py();
+ double Pz_Lab = LVEnergyImpulsionLAB.Pz();
+
+ double Px_CM = LVEnergyImpulsionCM.Px();
+ double Py_CM = LVEnergyImpulsionCM.Py();
+ double Pz_CM = LVEnergyImpulsionCM.Pz();
+
+ double Pper_Lab = sqrt(Px_Lab*Px_Lab+Py_Lab*Py_Lab);
+ double Pper_CM = sqrt(Px_CM*Px_CM+Py_CM*Py_CM);
+
+ double test_Lab = 0.5*log((proton->Mass()+ELab + Pz_Lab)/(proton->Mass()+ELab - Pz_Lab));
+ double test_CM = 0.5*log((proton->Mass()+ECM + Pz_CM)/(proton->Mass()+ECM - Pz_CM));
+
+ cout << "********************************" << endl;
+ cout << Pper_Lab << " " << Pper_CM << endl;
+ cout << Rapidity_Lab << " " << test_Lab << endl;
+ cout << Rapidity_CM << " " << test_CM << endl;
+ cout << endl;*/
+
+ ThetaLab = ThetaLab/deg;
+ PhiLab = PhiLab/deg;
+
+ }//end loop Hira
+ }//end if Hira
+
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End(){
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("ThicknessCsI",&ThicknessCsI,"ThicknessCsI/D");
+ RootOutput::getInstance()->GetTree()->Branch( "ThetaLab" , &ThetaLab , "ThetaLab/D" );
+ RootOutput::getInstance()->GetTree()->Branch( "PhiLab" , &PhiLab , "PhiLab/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM,"ThetaCM/D") ;
+ RootOutput::getInstance()->GetTree()->Branch( "E_ThinSi" , &E_ThinSi , "E_ThinSi/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "E_ThickSi" , &E_ThickSi , "E_ThickSi/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "E_CsI" , &E_CsI , "E_CsI/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "ELab" , &ELab , "ELab/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "ECM" , &ECM , "ECM/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "ECM_initial" , &ECM_initial , "ECM_initial/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "ECM_initial_Hira" , &ECM_initial_Hira , "ECM_initial_Hira/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy,"ExcitationEnergy/D") ;
+ RootOutput::getInstance()->GetTree()->Branch( "X" , &X , "X/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "Y" , &Y , "Y/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "Z" , &Z , "Z/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "TelescopeNumber" , &TelescopeNumber , "TelescopeNumber/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "X_target" , &X_target , "X_target/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "Y_target" , &Y_target , "Y_target/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch("InitialEnergy",&InitialEnergy,"InitialEnergy/D");
+ RootOutput::getInstance()->GetTree()->Branch("InitialEnergy_Hira",&InitialEnergy_Hira,"InitialEnergy_Hira/D");
+ RootOutput::getInstance()->GetTree()->Branch("Rapidity_CM",&Rapidity_CM,"Rapidity_CM/D");
+ RootOutput::getInstance()->GetTree()->Branch("Rapidity_Lab",&Rapidity_Lab,"Rapidity_Lab/D");
+ RootOutput::getInstance()->GetTree()->Branch("Pper",&Pper,"Pper/D");
+ RootOutput::getInstance()->GetTree()->Branch("PID",&PID,"PID/D");
+ RootOutput::getInstance()->GetTree()->Branch("Particle",&Particle,"Particle/I");
+ //RootOutput::getInstance()->GetTree()-> Branch("InteractionCoordinates","TInteractionCoordinates",&InteractionCoordinates);
+ //RootOutput::getInstance()->GetTree()->Branch("InitialConditions","TInitialConditions",&InitialConditions);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+ ThicknessCsI= -100;
+ E_ThinSi = -100;
+ E_ThickSi = -100;
+ E_CsI = -100;
+ ELab = -100;
+ ECM_initial = -100;
+ ECM_initial_Hira = -100;
+ ECM = -100;
+ ThetaLab = -100;
+ PhiLab = -100;
+ ThetaCM = -100;
+ ExcitationEnergy = -100;
+ X = -100;
+ Y = -100;
+ Z = -100;
+ TelescopeNumber = -1;
+ X_target = -100;
+ Y_target = -100;
+ InitialEnergy = -100;
+ InitialEnergy_Hira = -100;
+ Rapidity_Lab = -100;
+ Rapidity_CM = -100;
+ Pper = -100;
+ PID = -100;
+ Particle = -1;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct(){
+ return (NPL::VAnalysis*) new Analysis();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
+
+ proxy p;
+}
\ No newline at end of file
diff --git a/Projects/Hira10/Analysis.h b/Projects/Hira10/Analysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..595ed65c99ea916db95624b07aeefbfb4ec871a8
--- /dev/null
+++ b/Projects/Hira10/Analysis.h
@@ -0,0 +1,115 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 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: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include"NPVAnalysis.h"
+#include"THiraPhysics.h"
+#include "TInitialConditions.h"
+#include "TInteractionCoordinates.h"
+#include "NPEnergyLoss.h"
+#include "NPReaction.h"
+#include "TRandom3.h"
+#include "TF1.h"
+#include "TLorentzVector.h"
+
+class Analysis: public NPL::VAnalysis{
+public:
+ Analysis();
+ ~Analysis();
+
+public:
+ void Init();
+ void TreatEvent();
+ void End();
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+ double event;
+ double good_event;
+ double ProtonEnergy;
+
+private:
+ double ThicknessCsI;
+ double TargetThickness;
+ double ExcitationEnergy;
+ double ELab;
+ double ECM;
+ double ECM_initial;
+ double ECM_initial_Hira;
+ double E_ThinSi;
+ double E_ThickSi;
+ double E_CsI;
+ double PhiLab;
+ double ThetaLab;
+ double ThetaLab_simu;
+ double ThetaCM;
+ double X,Y,Z;
+ double TelescopeNumber;
+ double EnergyThreshold;
+ double X_target;
+ double Y_target;
+ double InitialEnergy;
+ double InitialEnergy_Hira;
+ double Rapidity_Lab;
+ double Rapidity_CM;
+ double Pper;
+ double PID;
+ int Particle;
+
+
+
+ NPL::Reaction* TransferReaction;
+
+ // intermediate variable
+ TRandom3 Rand;
+
+
+ TF1* f_proton;
+ TF1* f_deuton;
+ TF1* f_triton;
+ TF1* f_3He;
+
+ NPL::EnergyLoss Proton_CsI ;
+ NPL::EnergyLoss Deuton_CsI ;
+ NPL::EnergyLoss Triton_CsI ;
+ NPL::EnergyLoss Triton_CD2 ;
+ NPL::EnergyLoss Triton_CH2 ;
+ NPL::EnergyLoss Deuteron_CH2 ;
+ NPL::EnergyLoss He3_CD2;
+ NPL::EnergyLoss He3_CsI;
+ NPL::EnergyLoss He4_CH2;
+
+ NPL::Nucleus* proton;
+ NPL::Nucleus* deuton;
+ NPL::Nucleus* triton;
+ NPL::Nucleus* helium3;
+ NPL::Nucleus* alpha;
+ NPL::Nucleus* beam;
+ NPL::Nucleus* target;
+
+ THiraPhysics* Hira;
+ TInitialConditions* InitialConditions;
+ TInteractionCoordinates* InteractionCoordinates;
+};
+#endif
\ No newline at end of file
diff --git a/Projects/Hira10/CMakeLists.txt b/Projects/Hira10/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6806d778cdd5240538c93278f2ac5dcce3429083
--- /dev/null
+++ b/Projects/Hira10/CMakeLists.txt
@@ -0,0 +1,2 @@
+cmake_minimum_required (VERSION 2.8)
+include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
diff --git a/Projects/Hira10/RunToTreat.txt b/Projects/Hira10/RunToTreat.txt
new file mode 100755
index 0000000000000000000000000000000000000000..98c0bbfb0b9a849f23bb10d4cdb2f148224b725f
--- /dev/null
+++ b/Projects/Hira10/RunToTreat.txt
@@ -0,0 +1,13 @@
+TTreeName
+ SimulatedTree
+RootFileName
+ ../../Outputs/Simulation/hira_acceptance.root
+ %../../Outputs/Simulation/hiraU_deutonpbuu_nucl.root
+ %../../Outputs/Simulation/hiraU_tritonpbuu_nucl.root
+ %../../Outputs/Simulation/hiraU_3Hepbuu_nucl.root
+ %../../Outputs/Simulation/hiraU_alphapbuu_nucl.root
+
+
+ %../../Outputs/Simulation/hiraU_protonpbuu_nucl.root
+ %../../Outputs/Simulation/hiraU_flat_p_nucl.root
+ %../../Outputs/Simulation/hiraU_protonsource_nucl.root
diff --git a/Projects/Hira10/cmake_install.cmake b/Projects/Hira10/cmake_install.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..3d5e9eb942c7e9093d7839d764ffe76d68b61ffb
--- /dev/null
+++ b/Projects/Hira10/cmake_install.cmake
@@ -0,0 +1,39 @@
+# Install script for directory: /Users/pierremorfouace/Physics/NPTool/nptool/Projects/Hira10
+
+# Set the install prefix
+if(NOT DEFINED CMAKE_INSTALL_PREFIX)
+ set(CMAKE_INSTALL_PREFIX "/usr/local")
+endif()
+string(REGEX REPLACE "/$" "" CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}")
+
+# Set the install configuration name.
+if(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME)
+ if(BUILD_TYPE)
+ string(REGEX REPLACE "^[^A-Za-z0-9_]+" ""
+ CMAKE_INSTALL_CONFIG_NAME "${BUILD_TYPE}")
+ else()
+ set(CMAKE_INSTALL_CONFIG_NAME "Release")
+ endif()
+ message(STATUS "Install configuration: \"${CMAKE_INSTALL_CONFIG_NAME}\"")
+endif()
+
+# Set the component getting installed.
+if(NOT CMAKE_INSTALL_COMPONENT)
+ if(COMPONENT)
+ message(STATUS "Install component: \"${COMPONENT}\"")
+ set(CMAKE_INSTALL_COMPONENT "${COMPONENT}")
+ else()
+ set(CMAKE_INSTALL_COMPONENT)
+ endif()
+endif()
+
+if(CMAKE_INSTALL_COMPONENT)
+ set(CMAKE_INSTALL_MANIFEST "install_manifest_${CMAKE_INSTALL_COMPONENT}.txt")
+else()
+ set(CMAKE_INSTALL_MANIFEST "install_manifest.txt")
+endif()
+
+string(REPLACE ";" "\n" CMAKE_INSTALL_MANIFEST_CONTENT
+ "${CMAKE_INSTALL_MANIFEST_FILES}")
+file(WRITE "/Users/pierremorfouace/Physics/NPTool/nptool/Projects/Hira10/${CMAKE_INSTALL_MANIFEST}"
+ "${CMAKE_INSTALL_MANIFEST_CONTENT}")
diff --git a/Projects/Hira10/macro/.DS_Store b/Projects/Hira10/macro/.DS_Store
new file mode 100644
index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6
Binary files /dev/null and b/Projects/Hira10/macro/.DS_Store differ
diff --git a/Projects/Hira10/macro/ShowAcceptance.C b/Projects/Hira10/macro/ShowAcceptance.C
new file mode 100644
index 0000000000000000000000000000000000000000..9acabeb15d060c4c77be3ce5721ee1ae26937f40
--- /dev/null
+++ b/Projects/Hira10/macro/ShowAcceptance.C
@@ -0,0 +1,37 @@
+TChain* chain;
+
+/////////////////////////////////////////////////
+void LoadChain(){
+ chain = new TChain("PhysicsTree");
+ chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hira_acceptance.root");
+ return;
+}
+
+/////////////////////////////////////////////////
+void ShowAcceptance(){
+ gROOT->SetStyle("pierre_style");
+
+ LoadChain();
+
+ TCanvas *c1 = new TCanvas("c1","c1",800,800);
+ c1->cd();
+ chain->Draw("ThetaLab:ECM>>h1(220,0,220,180,0,180)","","colz");
+ TH2F *h1 = (TH2F*) gDirectory->FindObjectAny("h1");
+ h1->GetXaxis()->SetTitle("E_{c.m.} (MeV)");
+ h1->GetYaxis()->SetTitle("#theta_{lab} (deg)");
+ h1->GetXaxis()->CenterTitle();
+ h1->GetYaxis()->CenterTitle();
+ h1->Draw("colz");
+
+ TCanvas *c2 = new TCanvas("c2","c2",800,800);
+ c2->cd();
+ chain->Draw("ThetaCM:ECM>>h2(220,0,220,180,0,180)","","colz");
+ TH2F *h2 = (TH2F*) gDirectory->FindObjectAny("h2");
+ h2->GetXaxis()->SetTitle("E_{c.m.} (MeV)");
+ h2->GetYaxis()->SetTitle("#theta_{c.m.} (deg)");
+ h2->GetXaxis()->CenterTitle();
+ h2->GetYaxis()->CenterTitle();
+ h2->Draw("colz");
+
+
+}
diff --git a/Projects/Hira10/macro/ShowEcm.C b/Projects/Hira10/macro/ShowEcm.C
new file mode 100644
index 0000000000000000000000000000000000000000..46f87afb1cf3492f873f0351072c686fe2eabf16
--- /dev/null
+++ b/Projects/Hira10/macro/ShowEcm.C
@@ -0,0 +1,83 @@
+TChain* chain1=NULL ;
+TChain* chain2=NULL ;
+
+TCanvas* c1 = NULL;
+TH1F* heff = NULL;
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadChain(){
+ chain1 = new TChain("PhysicsTree");
+ chain1->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hiraU_protonpbuu_multi.root");
+
+ chain2 = new TChain("PhysicsTree");
+ //chain2->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hiraU_pid_nucl.root");
+ //chain2->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hiraU_flat_p_nucl.root");
+ //chain2->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hiraU_protonsource_nucl.root");
+ chain2->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/test.root");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadEfficiency(){
+ //TFile* f1 = new TFile("/Users/pierremorfouace/Physics/NPTool/nptool/Projects/Hira_upgrade/Efficiency/Efficiency_HiraU_flat.root");
+
+ TFile* f1 = new TFile("/Users/pierremorfouace/Physics/NPTool/nptool/Projects/Hira_upgrade/Efficiency_HiraU.root");
+ heff = (TH1F*) f1->FindObjectAny("heff_CM");
+ //heff->Sumw2();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void ShowEcm()
+{
+ gROOT->SetStyle("nptool");
+ LoadChain();
+ LoadEfficiency();
+
+ TFile *f1 = new TFile("/Users/pierremorfouace/Physics/NPTool/nptool/Inputs/EventGenerator/sn112e120_p_energyCM.root");
+ TH1F* hecm_pbuu_sn112 = (TH1F*) f1->FindObjectAny("henergy");
+ hecm_pbuu_sn112->SetMarkerColor(2);
+ hecm_pbuu_sn112->SetLineColor(2);
+ hecm_pbuu_sn112->Rebin(1);
+ hecm_pbuu_sn112->Scale(62);
+
+
+ /*chain1->Draw("ECM>>hECM1(50,0,300)","InitialEnergy>0","E1");
+ TH1F* hECM1 = (TH1F*) gDirectory->FindObjectAny("hECM1");
+ hECM1->SetMarkerColor(4);
+ hECM1->SetLineColor(4);
+ hECM1->SetMarkerStyle(8);*/
+ //hECM1->Sumw2();
+
+ chain2->Draw("ECM>>hECM2(50,0,300)","Particle==1 && ELab<200","E1");
+ TH1F* hECM2 = (TH1F*) gDirectory->FindObjectAny("hECM2");
+ hECM2->SetMarkerColor(2);
+ hECM2->SetLineColor(2);
+ hECM2->SetMarkerStyle(8);
+ //hECM2->Sumw2();
+
+ chain2->Draw("ECM_initial>>hECMi(50,0,300)","","E1");
+ TH1F* hECMi = (TH1F*) gDirectory->FindObjectAny("hECMi");
+ hECMi->SetMarkerColor(1);
+ hECMi->SetLineColor(1);
+ hECMi->SetMarkerStyle(8);
+ hECMi->Scale(1);
+
+ TH1F *hECM_correct = new TH1F("hECM_correct","hECM_correct",50,0,300);
+ hECM_correct->Divide(hECM2,heff,1,1);
+ hECM_correct->SetMarkerColor(8);
+ hECM_correct->SetLineColor(8);
+ hECM_correct->SetMarkerStyle(8);
+ hECM_correct->Scale(1.3);
+
+
+
+ //c1 = new TCanvas("c1","c1",800,600);
+ //c1->SetLogy();
+
+ //hECM1->Draw("E1same");
+
+ hECMi->Draw("E1");
+ hECM_correct->Draw("E1same");
+ hECM2->Draw("E1same");
+ hecm_pbuu_sn112->Draw("E1same");
+
+}
\ No newline at end of file
diff --git a/Projects/Hira10/macro/ShowPlot.C b/Projects/Hira10/macro/ShowPlot.C
new file mode 100644
index 0000000000000000000000000000000000000000..c67e429120c78fa7a535e4241a35b26437edb411
--- /dev/null
+++ b/Projects/Hira10/macro/ShowPlot.C
@@ -0,0 +1,234 @@
+TChain* chain1=NULL ;
+TChain* chain2=NULL ;
+
+TCanvas* c1 = NULL;
+TH1F* heff = NULL;
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadChain(){
+ chain1 = new TChain("PhysicsTree");
+ chain1->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hiraU_protonpbuu_multi.root");
+
+ chain2 = new TChain("PhysicsTree");
+ chain2->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/hiraU_protonpbuu_nucl.root");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadEfficiency(){
+ TFile* f1 = new TFile("/Users/pierremorfouace/Physics/NPTool/nptool/Projects/Hira_upgrade/Efficiency/Efficiency_HiraU_flat.root");
+ heff = (TH1F*) f1->FindObjectAny("heff_CM");
+ //heff->Sumw2();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void ShowPlot()
+{
+ gROOT->SetStyle("nptool");
+ LoadChain();
+ LoadEfficiency();
+
+ TFile *f1 = new TFile("/Users/pierremorfouace/Physics/NPTool/nptool/Inputs/EventGenerator/sn112e120_p_energyCM.root");
+ TH1F* hecm_pbuu_sn112 = (TH1F*) f1->FindObjectAny("henergy");
+ hecm_pbuu_sn112->SetMarkerColor(2);
+ hecm_pbuu_sn112->SetLineColor(2);
+ hecm_pbuu_sn112->Rebin(1);
+ hecm_pbuu_sn112->Scale(31);
+
+ TCanvas *c1 = new TCanvas("c1","c1",1200,600);
+ c1->Divide(2,1);
+
+ // ELab //
+ c1->cd(1);
+ TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
+ pad1->SetBottomMargin(0); // Upper and lower plot are joined
+ pad1->SetGridx(); // Vertical grid
+ pad1->Draw(); // Draw the upper pad: pad1
+ pad1->cd()->SetLogy(); // pad1 becomes the current pad
+
+ chain1->Draw("ELab>>hELab1(20,0,200)","InitialEnergy>0","E1");
+ TH1F* hELab1 = (TH1F*) gDirectory->FindObjectAny("hELab1");
+ hELab1->SetMarkerColor(4);
+ hELab1->SetLineColor(4);
+ hELab1->SetMarkerStyle(8);
+ hELab1->SetName("multi");
+ hELab1->GetXaxis()->SetTitle("E_{lab} (MeV)");
+ hELab1->GetYaxis()->SetTitle("counts (a.u.)");
+ hELab1->GetXaxis()->SetRangeUser(0,220);
+
+ chain2->Draw("InitialEnergy_Hira>>hEi(20,0,200)","InitialEnergy>0","E1");
+ TH1F* hEi = (TH1F*) gDirectory->FindObjectAny("hEi");
+ hEi->SetMarkerColor(1);
+ hEi->SetLineColor(1);
+ hEi->SetMarkerStyle(8);
+ hEi->SetName("normal");
+ hEi->GetXaxis()->SetTitle("E_{lab} (MeV)");
+ hEi->GetYaxis()->SetTitle("counts (a.u.)");
+ hEi->GetXaxis()->SetRangeUser(0,220);
+
+ chain2->Draw("ELab>>hELab2(20,0,200)","InitialEnergy_Hira>0","E1");
+ TH1F* hELab2 = (TH1F*) gDirectory->FindObjectAny("hELab2");
+ hELab2->SetMarkerColor(2);
+ hELab2->SetLineColor(2);
+ hELab2->SetMarkerStyle(8);
+ hELab2->SetName("nucl");
+
+ hEi->Draw("E1");
+ hELab1->Draw("E1same");
+ hELab2->Draw("E1same");
+
+ hELab1->GetYaxis()->SetTitleSize(25);
+ hELab1->GetYaxis()->SetTitleFont(43);
+ hELab1->GetYaxis()->SetTitleOffset(1.55);
+
+
+ TLegend *leg = new TLegend(0.6,0.7,0.9,0.95);
+ leg->SetBorderSize(1);
+ leg->SetLineColor(1);
+ leg->AddEntry("normal","Original spectrum","lp");
+ leg->AddEntry("multi","Multiple scattering","lp");
+ leg->AddEntry("nucl","Nuclear reaction loss","lp");
+ leg->Draw();
+
+ c1->cd(1);
+ TPad *pad2 = new TPad("pad2", "pad2", 0.0, 0.05, 1, 0.3);
+ pad2->SetTopMargin(0);
+ pad2->SetBottomMargin(0.36);
+ pad2->SetGridx();
+ pad2->Draw();
+ pad2->cd(); // pad2 becomes the current pad
+ // Define the ratio plot
+ TH1F *hratio_lab = (TH1F*)hELab2->Clone("hratio_lab");
+ hratio_lab->Sumw2();
+ hratio_lab->Divide(hEi);
+ hratio_lab->SetMarkerColor(2);
+ hratio_lab->SetLineColor(2);
+ hratio_lab->SetMarkerStyle(8);
+ hratio_lab->SetMarkerSize(1);
+ hratio_lab->GetXaxis()->SetRangeUser(0,220);
+ hratio_lab->GetYaxis()->SetRangeUser(0.5,1.05);
+ hratio_lab->Draw("E1");
+
+ hratio_lab->GetXaxis()->SetTitle("E_{lab} (MeV)");
+ hratio_lab->GetXaxis()->CenterTitle();
+ hratio_lab->GetYaxis()->SetTitle("ratio");
+ hratio_lab->GetYaxis()->SetNdivisions(505);
+ hratio_lab->GetYaxis()->SetTitleSize(25);
+ hratio_lab->GetYaxis()->SetTitleFont(43);
+ hratio_lab->GetYaxis()->SetTitleOffset(1.55);
+ hratio_lab->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
+ hratio_lab->GetYaxis()->SetLabelSize(20);
+
+ // X axis ratio plot settings
+ hratio_lab->GetXaxis()->SetTitleSize(25);
+ hratio_lab->GetXaxis()->SetTitleFont(43);
+ hratio_lab->GetXaxis()->SetTitleOffset(4.);
+ hratio_lab->GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
+ hratio_lab->GetXaxis()->SetLabelSize(20);
+
+
+ TH1F *hratio_lab2 = (TH1F*)hELab1->Clone("hratio_lab");
+ hratio_lab2->Sumw2();
+ hratio_lab2->Divide(hEi);
+ hratio_lab2->SetMarkerColor(4);
+ hratio_lab2->SetLineColor(4);
+ hratio_lab2->SetMarkerStyle(8);
+ hratio_lab2->SetMarkerSize(1);
+ hratio_lab2->Draw("E1same");
+
+
+
+ // ECM //
+ c1->cd(2);
+ TPad *pad3 = new TPad("pad3", "pad3", 0, 0.3, 1, 1.0);
+ pad3->SetBottomMargin(0); // Upper and lower plot are joined
+ pad3->SetGridx(); // Vertical grid
+ pad3->Draw(); // Draw the upper pad: pad1
+ pad3->cd()->SetLogy(); // pad1 becomes the current pad
+
+ chain1->Draw("ECM>>hECM1(20,0,200)","","E1");
+ TH1F* hECM1 = (TH1F*) gDirectory->FindObjectAny("hECM1");
+ hECM1->SetMarkerColor(4);
+ hECM1->SetLineColor(4);
+ hECM1->SetMarkerStyle(8);
+ hECM1->GetXaxis()->SetTitle("E_{c.m.} (MeV)");
+ hECM1->GetYaxis()->SetTitle("counts (a.u.)");
+ hECM1->GetXaxis()->SetRangeUser(0,200);
+
+ chain2->Draw("ECM_initial_Hira>>hECMi(20,0,200)","InitialEnergy>0 && InitialEnergy<200","E1");
+ TH1F* hECMi = (TH1F*) gDirectory->FindObjectAny("hECMi");
+ hECMi->SetMarkerColor(1);
+ hECMi->SetLineColor(1);
+ hECMi->SetMarkerStyle(8);
+ hECMi->SetName("normal");
+ hECMi->GetXaxis()->SetTitle("E_{c.m.} (MeV)");
+ hECMi->GetYaxis()->SetTitle("counts (a.u.)");
+ hECMi->GetXaxis()->SetRangeUser(0,200);
+ hECMi->GetYaxis()->SetRangeUser(10,2e5);
+
+ chain2->Draw("ECM>>hECM2(20,0,200)","","E1");
+ TH1F* hECM2 = (TH1F*) gDirectory->FindObjectAny("hECM2");
+ hECM2->SetMarkerColor(2);
+ hECM2->SetLineColor(2);
+ hECM2->SetMarkerStyle(8);
+ //hECM2->Sumw2();
+
+ hECM1->GetYaxis()->SetTitleSize(25);
+ hECM1->GetYaxis()->SetTitleFont(43);
+ hECM1->GetYaxis()->SetTitleOffset(1.55);
+
+ hECMi->Draw("E1");
+ hECM1->Draw("E1same");
+ hECM2->Draw("E1same");
+
+ cout << "Mean value of the energy in the c.m" << endl;
+ cout << "Original: " << hECMi->GetMean() << endl;
+ cout << "Multiple: " << hECM1->GetMean() << endl;
+ cout << "Nuclear: " << hECM2->GetMean() << endl;
+
+ c1->cd(2);
+ TPad *pad4 = new TPad("pad4", "pad4", 0.0, 0.05, 1, 0.3);
+ pad4->SetTopMargin(0);
+ pad4->SetBottomMargin(0.36);
+ pad4->SetGridx();
+ pad4->Draw();
+ pad4->cd(); // pad4 becomes the current pad
+ // Define the ratio plot
+ TH1F *hratio_cm = (TH1F*)hECM2->Clone("hratio_cm");
+ hratio_cm->Sumw2();
+ hratio_cm->Divide(hECMi);
+ hratio_cm->SetMarkerColor(2);
+ hratio_cm->SetLineColor(2);
+ hratio_cm->SetMarkerStyle(8);
+ hratio_cm->SetMarkerSize(1);
+ hratio_cm->GetXaxis()->SetRangeUser(0,200);
+ hratio_cm->GetYaxis()->SetRangeUser(0.5,1.05);
+ hratio_cm->Draw("E1");
+
+ hratio_cm->GetXaxis()->SetTitle("E_{c.m.} (MeV)");
+ hratio_cm->GetXaxis()->CenterTitle();
+ hratio_cm->GetYaxis()->SetTitle("ratio");
+ hratio_cm->GetYaxis()->SetNdivisions(505);
+ hratio_cm->GetYaxis()->SetTitleSize(25);
+ hratio_cm->GetYaxis()->SetTitleFont(43);
+ hratio_cm->GetYaxis()->SetTitleOffset(1.55);
+ hratio_cm->GetYaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
+ hratio_cm->GetYaxis()->SetLabelSize(20);
+
+ // X axis ratio plot settings
+ hratio_cm->GetXaxis()->SetTitleSize(25);
+ hratio_cm->GetXaxis()->SetTitleFont(43);
+ hratio_cm->GetXaxis()->SetTitleOffset(4.);
+ hratio_cm->GetXaxis()->SetLabelFont(43); // Absolute font size in pixel (precision 3)
+ hratio_cm->GetXaxis()->SetLabelSize(20);
+
+ TH1F *hratio_cm2 = (TH1F*)hECM1->Clone("hratio_cm");
+ hratio_cm2->Sumw2();
+ hratio_cm2->Divide(hECMi);
+ hratio_cm2->SetMarkerColor(4);
+ hratio_cm2->SetLineColor(4);
+ hratio_cm2->SetMarkerStyle(8);
+ hratio_cm2->SetMarkerSize(1);
+ hratio_cm2->Draw("E1same");
+
+
+}
\ No newline at end of file
diff --git a/Projects/Hira10/macro/ShowResults.C b/Projects/Hira10/macro/ShowResults.C
new file mode 100644
index 0000000000000000000000000000000000000000..b44a8bff902163548c2e3985db0a9ea55e963598
--- /dev/null
+++ b/Projects/Hira10/macro/ShowResults.C
@@ -0,0 +1,129 @@
+#include"NPReaction.h"
+
+
+TChain* chain=NULL ;
+TCanvas* c1 = NULL;
+TCanvas* c2 = NULL;
+TCanvas* c3 = NULL;
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadChain(){
+ chain = new TChain("PhysicsTree");
+ chain->Add("../../Outputs/Analysis/hiraU_protonpbuu_nucl.root");
+ //chain->Add("../../Outputs/Analysis/hiraU_flat_p_nucl.root");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void ShowResults(){
+ gROOT->SetStyle("nptool");
+
+ LoadChain();
+
+ c1 = new TCanvas("Example1","Example1",800,800);
+ c1->Divide(2,2);
+
+ // Light Particle ID //
+ // E-DE
+ c1->cd(1);
+ chain->Draw("E_ThickSi:E_CsI>>hIDE(1000,0,220,1000,0,20)","","colz");
+ TH1F* hIDE = (TH1F*) gDirectory->FindObjectAny("hIDE");
+ hIDE->GetXaxis()->SetTitle("E_{CsI} (MeV)");
+ hIDE->GetYaxis()->SetTitle("E_{Si} (MeV)");
+
+ // TotalEnergy
+ int bin=50;
+ double Emin = 0;
+ double Emax = 300;
+ c1->cd(2);
+
+ chain->Draw(Form("InitialEnergy_Hira>>hEi(%d,%f,%f)",bin,Emin,Emax),"InitialEnergy>0","E1");
+ TH1F* hEi = (TH1F*) gDirectory->FindObjectAny("hEi");
+ hEi->GetXaxis()->SetTitle("E_{lab} (MeV)");
+ hEi->GetYaxis()->SetTitle("counts");
+ hEi->SetMarkerColor(2);
+ hEi->SetLineColor(2);
+ hEi->Sumw2();
+
+ chain->Draw(Form("ELab>>hEd(%d,%f,%f)",bin,Emin,Emax),"InitialEnergy>0","E1same");
+ TH1F* hEd = (TH1F*) gDirectory->FindObjectAny("hEd");
+ hEd->SetMarkerColor(4);
+ hEd->SetLineColor(4);
+ hEd->Sumw2();
+
+
+ hEi->Draw("E1");
+ hEd->Draw("E1same");
+
+
+ c1->cd(3);
+ TH1F* heff = new TH1F("heff","heff",bin,Emin,Emax);
+ heff->Sumw2();
+ heff->GetXaxis()->SetTitle("E_{lab} (MeV)");
+ heff->GetYaxis()->SetTitle("Eff (%)");
+ heff->GetYaxis()->SetRangeUser(0,100);
+ heff->GetXaxis()->SetRangeUser(0,200);
+ heff->SetMarkerColor(4);
+ heff->SetLineColor(4);
+
+ heff->Divide(hEd,hEi,1,1);
+ heff->Draw("E1");
+
+ TH1F* heff_nucl = new TH1F("heff_nucl","heff_nucl",bin,Emin,Emax);
+ heff_nucl->Sumw2();
+ heff_nucl->SetMarkerColor(8);
+ heff_nucl->SetLineColor(8);
+
+
+ /*heff->Fit("pol1");
+
+ TF1* f = heff->GetFunction("pol1");
+ f->SetLineWidth(2);
+ f->SetLineColor(kOrange-3);
+ f->SetNpx(1000);*/
+
+ c1->cd(4);
+ chain->Draw("ELab:ThetaLab>>htheta(160,0,120,500,0,200)","InitialEnergy>0","colz");
+ TH1F* htheta = (TH1F*) gDirectory->FindObjectAny("htheta");
+ htheta->GetXaxis()->SetTitle("#theta_{lab} (deg)");
+ htheta->GetYaxis()->SetTitle("E_{lab} (MeV)");
+
+
+ c2 = new TCanvas("CM","CM",800,800);
+ c2->Divide(2,1);
+ c2->cd(1)->SetLogy();
+ chain->Draw("ECM_initial_Hira>>hECMi(50,0,300)","InitialEnergy>0 && InitialEnergy<200","E1");
+ TH1F* hECMi = (TH1F*) gDirectory->FindObjectAny("hECMi");
+ hECMi->GetXaxis()->SetTitle("E_{CM} (MeV)");
+ hECMi->SetMarkerColor(2);
+ hECMi->SetLineColor(2);
+ hECMi->Sumw2();
+
+ chain->Draw("ECM>>hECM(50,0,300)","ELab<200 && ELab>0","E1same");
+ TH1F* hECM = (TH1F*) gDirectory->FindObjectAny("hECM");
+ hECM->SetMarkerColor(4);
+ hECM->SetLineColor(4);
+ hECM->Sumw2();
+
+ hECMi->Draw("E1");
+ hECM->Draw("E1same");
+
+ c2->cd(2);
+ chain->Draw("ECM:ThetaLab>>hECMtheta(500,0,80,500,0,120)","","colz");
+ TH1F* hECMtheta = (TH1F*) gDirectory->FindObjectAny("hECMtheta");
+ hECMtheta->GetXaxis()->SetTitle("#theta_{lab} (deg)");
+ hECMtheta->GetYaxis()->SetTitle("E_{CM} (MeV)");
+
+
+ c3 = new TCanvas("Eff","Eff",0,0,800,800);
+ c3->cd();
+ TH1F* heff_CM = new TH1F("heff_CM","heff_CM",50,0,300);
+ heff_CM->Sumw2();
+ heff_CM->GetXaxis()->SetTitle("E_{CM} (MeV)");
+ heff_CM->GetYaxis()->SetTitle("Eff");
+ heff_CM->Divide(hECM,hECMi,1,1);
+ heff_CM->Draw("E1");
+
+ heff_CM->SaveAs("Efficiency_HiraU.root");
+
+
+}
diff --git a/Projects/MUGAST_LISE/68Nidp.reaction b/Projects/MUGAST_LISE/68Nidp.reaction
new file mode 100644
index 0000000000000000000000000000000000000000..9d1be70f582d8abf38409e5cf0a82a52294a07c2
--- /dev/null
+++ b/Projects/MUGAST_LISE/68Nidp.reaction
@@ -0,0 +1,33 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 68Ni
+ Energy= 1700
+ SigmaEnergy= 10
+ ExcitationEnergy= 0 MeV
+ SigmaThetaX= 3 mrad
+ SigmaPhiY= 3 mrad
+ SigmaX= 25 mm
+ SigmaY= 25 mm
+ MeanThetaX= 1 mrad
+ MeanPhiY= 1 mrad
+ MeanX= 0 mm
+ MeanY= 0 mm
+ %EnergyProfilePath= eLise.txt EL
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 68Ni
+ Target= 2H
+ Light= 1H
+ Heavy= 69Ni
+ ExcitationEnergyLight= 0.0 MeV
+ ExcitationEnergyHeavy= 1.0 MeV
+ CrossSectionPath= fe61_f5_2_250.txt h
+ ShootLight= 1
+ ShootHeavy= 0
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
diff --git a/Projects/MUGAST_LISE/Analysis.cxx b/Projects/MUGAST_LISE/Analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d18aca94d9e02a98a8c3d964c1dae2fd60693155
--- /dev/null
+++ b/Projects/MUGAST_LISE/Analysis.cxx
@@ -0,0 +1,285 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2014 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: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include<iostream>
+using namespace std;
+
+#include "Analysis.h"
+#include "NPFunction.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPOptionManager.h"
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis(){
+}
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init() {
+ // initialize input and output branches
+ InitOutputBranch();
+ InitInputBranch();
+ myInit = new TInitialConditions();
+ // get MUST2 and Gaspard objects
+ M2 = (TMust2Physics*) m_DetectorManager -> GetDetector("M2Telescope");
+ GD = (GaspardTracker*) m_DetectorManager -> GetDetector("GaspardTracker");
+
+ // get reaction information
+ myReaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ OriginalBeamEnergy = myReaction.GetBeamEnergy();
+
+ // target thickness
+ TargetThickness = m_DetectorManager->GetTargetThickness();
+ string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+ // Cryo target case
+ WindowsThickness = m_DetectorManager->GetWindowsThickness();
+ string WindowsMaterial = m_DetectorManager->GetWindowsMaterial();
+
+ // energy losses
+ string light=NPL::ChangeNameToG4Standard(myReaction.GetNucleus3().GetName());
+ string beam=NPL::ChangeNameToG4Standard(myReaction.GetNucleus1().GetName());
+
+ LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
+ LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100);
+ LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
+ BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
+
+ if(WindowsThickness){
+ BeamWindow= new NPL::EnergyLoss(beam+"_"+WindowsMaterial+".G4table","G4Table",100);
+ LightWindow= new NPL::EnergyLoss(light+"_"+WindowsMaterial+".G4table","G4Table",100);
+ }
+
+ else{
+ BeamWindow= NULL;
+ LightWindow=NULL;
+ }
+
+ // initialize various parameters
+ Rand = TRandom3();
+ DetectorNumber = 0;
+ ThetaNormalTarget = 0;
+ ThetaM2Surface = 0;
+ Si_E_M2 = 0;
+ CsI_E_M2 = 0;
+ Energy = 0;
+ ThetaGDSurface = 0;
+ X = 0;
+ Y = 0;
+ Z = 0;
+ dE = 0;
+ dTheta=0;
+ BeamDirection = TVector3(0,0,1);
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent() {
+ // Reinitiate calculated variable
+ ReInitValue();
+ //double zImpact = Rand.Uniform(-TargetThickness*0.5,TargetThickness*0.5);
+ double zImpact = 0 ;
+ BeamImpact = TVector3(0,0,zImpact);
+ // determine beam energy for a randomized interaction point in target
+ // 1% FWHM randominastion (E/100)/2.35
+ myReaction.SetBeamEnergy(Rand.Gaus(myInit->GetIncidentFinalKineticEnergy(),myInit->GetIncidentFinalKineticEnergy()/235));
+ OriginalThetaLab = myInit->GetThetaLab_WorldFrame(0);
+ OriginalELab = myInit->GetKineticEnergy(0);
+
+ //////////////////////////// LOOP on MUST2 //////////////////
+ for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
+ /************************************************/
+ //Part 0 : Get the usefull Data
+ // MUST2
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ X = M2 -> GetPositionOfInteraction(countMust2).X();
+ Y = M2 -> GetPositionOfInteraction(countMust2).Y();
+ Z = M2 -> GetPositionOfInteraction(countMust2).Z();
+
+ ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2= M2->CsI_E[countMust2];
+
+ // if CsI
+ if(CsI_E_M2>0 ){
+ // The energy in CsI is calculate form dE/dx Table because
+ Energy = CsI_E_M2;
+ Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
+ Energy+=Si_E_M2;
+ }
+
+ else
+ Energy = Si_E_M2;
+
+ // Evaluate energy using the thickness
+ ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
+ // Target Correction
+ ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2.-zImpact, ThetaNormalTarget);
+
+ if(LightWindow)
+ ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = myReaction.ReconstructRelativistic( ELab , ThetaLab );
+ ThetaLab=ThetaLab/deg;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = myReaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ /************************************************/
+ }//end loop MUST2
+
+ ////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////// LOOP on GASPARD //////////////////
+ if(GD->GetEnergyDeposit()>0){
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaGDSurface = 0;
+ ThetaNormalTarget = 0;
+ TVector3 HitDirection = GD -> GetPositionOfInteraction() - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ X = GD -> GetPositionOfInteraction().X();
+ Y = GD -> GetPositionOfInteraction().Y();
+ Z = GD -> GetPositionOfInteraction().Z();
+
+ ThetaGDSurface = HitDirection.Angle( TVector3(0,0,1) ) ;
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Energy = GD->GetEnergyDeposit();
+ // Target Correction
+ ELab = LightTarget.EvaluateInitialEnergy( Energy ,TargetThickness*0.5-zImpact, ThetaNormalTarget);
+
+ if(LightWindow)
+ ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
+
+ dE= myInit->GetKineticEnergy(0) - ELab ;
+ dTheta = (myInit->GetThetaLab_WorldFrame(0)*deg-ThetaLab)/deg ;
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = myReaction.ReconstructRelativistic( ELab , ThetaLab );
+
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = myReaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ ThetaLab=ThetaLab/deg;
+
+ /************************************************/
+ }//end loop GASPARD
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End(){
+}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("Run",&Run,"Run/I");
+ RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D");
+ RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D");
+ RootOutput::getInstance()->GetTree()->Branch("Z",&Z,"Z/D");
+ RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalELab",&OriginalELab,"OriginalELab/D");
+
+
+ RootOutput::getInstance()->GetTree()->Branch("dTheta",&dTheta,"dTheta/D");
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ RootInput::getInstance()->GetChain()->SetBranchStatus("Run",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Run",&Run);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&myInit);
+}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+ Ex = -1000 ;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
+ X = -1000;
+ Y = -1000;
+ Z = -1000;
+ dE= -1000;
+ dTheta= -1000;
+ OriginalThetaLab = -1000;
+ OriginalELab = -1000;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the AnalysisFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct(){
+ return (NPL::VAnalysis*) new Analysis();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+class proxy_analysis{
+ public:
+ proxy_analysis(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+};
+
+proxy_analysis p_analysis;
+}
+
diff --git a/Projects/MUGAST_LISE/Analysis.h b/Projects/MUGAST_LISE/Analysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..5b401c8b9ae850de6edb928dfa00a0e451666dc4
--- /dev/null
+++ b/Projects/MUGAST_LISE/Analysis.h
@@ -0,0 +1,97 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 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: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include"NPVAnalysis.h"
+#include"NPEnergyLoss.h"
+#include"NPReaction.h"
+#include"RootOutput.h"
+#include"RootInput.h"
+#include "TMust2Physics.h"
+#include "GaspardTracker.h"
+#include "TInitialConditions.h"
+#include <TRandom3.h>
+#include <TVector3.h>
+#include <TMath.h>
+
+class Analysis: public NPL::VAnalysis{
+ public:
+ Analysis();
+ ~Analysis();
+
+ public:
+ void Init();
+ void TreatEvent();
+ void End();
+
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+ private:
+ double Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+ NPL::Reaction myReaction;
+ // Energy loss table: the G4Table are generated by the simulation
+ NPL::EnergyLoss LightTarget;
+ NPL::EnergyLoss LightAl;
+ NPL::EnergyLoss LightSi;
+ NPL::EnergyLoss BeamCD2;
+ NPL::EnergyLoss* BeamWindow;
+ NPL::EnergyLoss* LightWindow;
+
+ double TargetThickness ;
+ double WindowsThickness;
+ // Beam Energy
+ double OriginalBeamEnergy ; // AMEV
+ double FinalBeamEnergy;
+
+ // intermediate variable
+ TVector3 BeamDirection;
+ TVector3 BeamImpact;
+ TRandom3 Rand ;
+ int Run;
+ int DetectorNumber ;
+ double ThetaNormalTarget;
+ double ThetaM2Surface ;
+ double Si_E_M2 ;
+ double CsI_E_M2 ;
+ double Energy ;
+ double OriginalThetaLab;
+ double OriginalELab;
+ double ThetaGDSurface ;
+ double X ;
+ double Y ;
+ double Z ;
+
+ double dE;
+ double dTheta;
+ // Branches and detectors
+ TMust2Physics* M2;
+ GaspardTracker* GD;
+ TInitialConditions* myInit ;
+
+};
+#endif
diff --git a/Projects/MUGAST_LISE/CMakeLists.txt b/Projects/MUGAST_LISE/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..22c74affdfc45019bdda2594f8439c52d4ab97ec
--- /dev/null
+++ b/Projects/MUGAST_LISE/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required (VERSION 2.8)
+# Setting the policy to match Cmake version
+cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
+# include the default NPAnalysis cmake file
+include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
diff --git a/Projects/MUGAST_LISE/DetectorConfiguration/MUGAST_LISE.detector b/Projects/MUGAST_LISE/DetectorConfiguration/MUGAST_LISE.detector
new file mode 100644
index 0000000000000000000000000000000000000000..60f8969ccd92c5a9b71b0abc05199b1d956b81c4
--- /dev/null
+++ b/Projects/MUGAST_LISE/DetectorConfiguration/MUGAST_LISE.detector
@@ -0,0 +1,116 @@
+%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+Target
+ THICKNESS= 9.5
+ RADIUS= 7.5 mm
+ MATERIAL= CD2
+ Angle= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+AnnularS1
+ Z= -100
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 0 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 45 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= -1
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 90 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 135 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 180 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 225 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 270 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GaspardTracker Trapezoid
+ R= -90 mm
+ THETA= 60 deg
+ PHI= 315 deg
+ BETA= -60 0 -90 deg
+ FIRSTSTAGE= 1
+ SECONDSTAGE= 0
+ THIRDSTAGE= 0
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Tiara Barrel
+ InnerBarrel= 0
+ OuterBarrel= 0
+ Chamber= 0
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X1_Y1= 115.12 -11.35 1540.94 mm
+ X1_Y128= 24.33 -11.60 1900.81 mm
+ X128_Y1= 103.64 -103.42 1206.51 mm
+ X128_Y128= 12.89 -103.90 1602.29 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
diff --git a/Projects/MUGAST_LISE/alpha.source b/Projects/MUGAST_LISE/alpha.source
new file mode 100644
index 0000000000000000000000000000000000000000..c2fb3781ce066b164f6e5882c4490356031ecc00
--- /dev/null
+++ b/Projects/MUGAST_LISE/alpha.source
@@ -0,0 +1,17 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy are given in MeV , Position in mm %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Isotropic
+ EnergyLow= 5 MeV
+ EnergyHigh= 5 MeV
+ HalfOpenAngleMin= 90 deg
+ HalfOpenAngleMax= 180 deg
+ x0= 0 mm
+ y0= 0 mm
+ z0= 0 mm
+ Particle= alpha
+ ExcitationEnergy= 0 MeV
+
+% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
diff --git a/Projects/MUGAST_LISE/show.cxx b/Projects/MUGAST_LISE/show.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6522307216c39cc8babc9596a59365754a91925e
--- /dev/null
+++ b/Projects/MUGAST_LISE/show.cxx
@@ -0,0 +1,19 @@
+void show(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/PhysicsTree.root","READ");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ tree->Draw("ELab:ThetaLab>>h(360,0,180,360,0,20)","ELab>0","colz");
+ NPL::Reaction r("68Ni(d,p)69Ni@680");
+ r.GetKinematicLine3()->Draw("c");
+ TCanvas*c =new TCanvas();
+ c->Divide(2,1);
+ c->cd(1);
+ tree->Draw("ELab:OriginalELab>>hE(500,0,16,500,0,16)","ELab>0","colz");
+ TLine* lineE = new TLine(0,0,16,16);
+ lineE->Draw();
+ c->cd(2);
+ tree->Draw("ThetaLab:OriginalThetaLab>>hT(500,100,180,500,100,180)","ELab>0","colz");
+ TLine* lineT = new TLine(100,100,180,180);
+ lineT->Draw();
+}
diff --git a/Projects/macros/GeometricalEfficiency.C b/Projects/macros/GeometricalEfficiency.C
index b5c7fef1e7e8f264d53dcfee8cd646c6a7f2e9ef..26b0d5daed665f4985aa44f8d2679adeaa98aac6 100644
--- a/Projects/macros/GeometricalEfficiency.C
+++ b/Projects/macros/GeometricalEfficiency.C
@@ -88,19 +88,20 @@ void GeometricalEfficiency(const char * fname = "myResult.root"){
}
hEmittTheta->Sumw2();
- hEmittThetaCM->Sumw2();
- hDetecTheta->Sumw2();
- hDetecThetaCM->Sumw2();
+ hEmittThetaCM->Sumw2();
+ hDetecTheta->Sumw2();
+ hDetecThetaCM->Sumw2();
- TCanvas *c0 = new TCanvas("c0", "Distrib",800,800);
- hEmittTheta->Draw("");
+ TCanvas *c0 = new TCanvas("c0", "Distrib",800,800);
+ hEmittTheta->Draw("");
+ hDetecTheta->SetMarkerColor(kAzure+7);
hDetecTheta->Draw("same");
// efficiency in lab frame in %
TCanvas *c = new TCanvas("c", "efficiency",800,800);
c->SetTopMargin(0.01);
c->SetRightMargin(0.03);
- TH1F *hEfficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 90);
+ TH1F *hEfficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
hEfficiency->Divide(hDetecTheta, hEmittTheta, 1, 1);
hEfficiency->GetXaxis()->SetTitle("#Theta (deg)");
hEfficiency->GetYaxis()->SetTitle("#epsilon");