From 154c10b14df6967313dce9a44ff0b970392be619 Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Mon, 24 Mar 2014 11:56:26 +0000
Subject: [PATCH] * Adding support for Ex distribution in Simulation - Similar
mecanism to the CS system - Allow user to provide a random Ex distribution
in txt or root file format * Changing Display of Gaspard, Paris and Must2
detector for the GMAP_LOI preparation
---
Inputs/EventGenerator/10He.reaction | 4 +-
NPLib/Physics/NPReaction.cxx | 30 +-
NPLib/Physics/NPReaction.h | 297 +-
NPSimulation/GASPARD/GaspardTrackerSquare.cc | 12 +-
.../GASPARD/GaspardTrackerTrapezoid.cc | 11 +-
NPSimulation/MUST2/MUST2Array.cc | 2419 ++++++++---------
NPSimulation/Paris/ParisCluster.cc | 15 +-
.../src/EventGeneratorTwoBodyReaction.cc | 3 +
8 files changed, 1404 insertions(+), 1387 deletions(-)
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 32a211b93..4e92aed97 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -36,8 +36,8 @@ TwoBodyReaction
Light= 3He
Heavy= 10He
ExcitationEnergyLight= 0.0
- ExcitationEnergyHeavy= 3.0
- CrossSectionPath= flat.txt CS10He
+ ExcitationEnergyHeavy= 1.3
+ CrossSectionPath= 11Li(d,3He)10He.txt CS10He
ShootLight= 1
ShootHeavy= 1
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index b26c96698..b25a373b6 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -88,6 +88,8 @@ Reaction::Reaction(){
++offset;
fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
+ fExcitationEnergyHist = NULL;
+
fshoot3=true;
fshoot4=true;
@@ -151,11 +153,11 @@ Reaction::Reaction(string reaction){
++offset;
fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
+ fCrossSectionHist = NULL;
+
fshoot3=true;
fshoot4=true;
-
-
initializePrecomputeVariable();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -176,6 +178,9 @@ Reaction::~Reaction(){
if(fCrossSectionHist)
delete fCrossSectionHist;
+
+ if(fExcitationEnergyHist)
+ delete fExcitationEnergyHist;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -208,6 +213,13 @@ double Reaction::ShootRandomThetaCM(){
SetThetaCM( theta=fCrossSectionHist->GetRandom()*deg );
return theta;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Reaction::ShootRandomExcitationEnergy(){
+ if(fExcitationEnergyHist){
+ SetExcitation4(fExcitationEnergyHist->GetRandom());
+ }
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
double &ThetaLab4, double &KineticEnergyLab4){
@@ -308,6 +320,7 @@ void Reaction::ReadConfigurationFile(string Path){
bool check_Heavy = false ;
bool check_ExcitationEnergy3 = false ;
bool check_ExcitationEnergy4 = false ;
+ bool check_ExcitationEnergyDistribution = false;
bool check_CrossSectionPath = false ;
bool check_shoot3 = false ;
bool check_shoot4 = false;
@@ -385,7 +398,16 @@ void Reaction::ReadConfigurationFile(string Path){
fExcitation4 = atof(DataBuffer.c_str()) * MeV;
if(fVerboseLevel==1) cout << "Excitation Energy Nuclei 4: " << fExcitation4 / MeV << " MeV" << endl;
}
-
+
+ else if (DataBuffer=="ExcitationEnergyDistribution="){
+ check_ExcitationEnergyDistribution = true;
+ string FileName,HistName;
+ ReactionFile >> FileName >> HistName;
+ if(fVerboseLevel==1) cout << "Reading Excitation Energy Distribution file: " << FileName << endl;
+ fExcitationEnergyHist = Read1DProfile(FileName, HistName );
+ fExcitation4 = 0 ;
+ }
+
else if (DataBuffer== "CrossSectionPath=") {
check_CrossSectionPath = true ;
string FileName,HistName;
@@ -442,7 +464,7 @@ void Reaction::ReadConfigurationFile(string Path){
///////////////////////////////////////////////////
// If all Token found toggle out
if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergy3
- && check_ExcitationEnergy4 && check_CrossSectionPath && check_shoot4 && check_shoot3)
+ && (check_ExcitationEnergy4 || check_ExcitationEnergyDistribution ) && check_CrossSectionPath && check_shoot4 && check_shoot3)
ReadingStatus = false;
}
}
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 72a2b9654..2c4130c5d 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -53,169 +53,172 @@ using namespace NPL;
using namespace std;
namespace NPL{
-
+
class Reaction{
-
- public: // Constructors and Destructors
- Reaction();
- // This constructor aim to provide a fast way to instantiate a reaction without input file
- // The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV
- Reaction(string);
- virtual ~Reaction();
-
- public: // Various Method
- void ReadConfigurationFile(string Path);
-
- private:
- int fVerboseLevel;
-
- private: // use for Monte Carlo simulation
- bool fshoot3;
- bool fshoot4;
-
- private: // use to display the kinematical line
- TGraph* fKineLine3 ;
- TGraph* fKineLine4 ;
- TGraph* fTheta3VsTheta4;
- TGraph* fLineBrho3;
- TGraph* fAngleLine;
- private:
- Beam* fNuclei1; // Beam
- Nucleus* fNuclei2; // Target
- Nucleus* fNuclei3; // Light ejectile
- Nucleus* fNuclei4; // Heavy ejectile
- double fQValue; // Q-value in MeV
- double fBeamEnergy; // Beam energy in MeV
- double fThetaCM; // Center-of-mass angle in radian
- double fExcitation3; // Excitation energy in MeV
- double fExcitation4; // Excitation energy in MeV
- TH1F* fCrossSectionHist;
-
- public:
+
+ public: // Constructors and Destructors
+ Reaction();
+ // This constructor aim to provide a fast way to instantiate a reaction without input file
+ // The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV
+ Reaction(string);
+ virtual ~Reaction();
+
+ public: // Various Method
+ void ReadConfigurationFile(string Path);
+
+ private:
+ int fVerboseLevel;
+
+ private: // use for Monte Carlo simulation
+ bool fshoot3;
+ bool fshoot4;
+
+ private: // use to display the kinematical line
+ TGraph* fKineLine3 ;
+ TGraph* fKineLine4 ;
+ TGraph* fTheta3VsTheta4;
+ TGraph* fLineBrho3;
+ TGraph* fAngleLine;
+ private:
+ Beam* fNuclei1; // Beam
+ Nucleus* fNuclei2; // Target
+ Nucleus* fNuclei3; // Light ejectile
+ Nucleus* fNuclei4; // Heavy ejectile
+ double fQValue; // Q-value in MeV
+ double fBeamEnergy; // Beam energy in MeV
+ double fThetaCM; // Center-of-mass angle in radian
+ double fExcitation3; // Excitation energy in MeV
+ double fExcitation4; // Excitation energy in MeV
+ TH1F* fCrossSectionHist; // Differential cross section in CM frame
+ TH1F* fExcitationEnergyHist; // Distribution of Excitation energy
+ public:
// Getters and Setters
- void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
- void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();}
- void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
- void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
+ void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
+ void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();}
+ void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
+ void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
// For retro compatibility
- void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
- void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
- void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
- void SetCrossSectionHist (TH1F* CrossSectionHist) {delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;}
-
- double GetBeamEnergy() const {return fBeamEnergy;}
- double GetThetaCM() const {return fThetaCM;}
- double GetExcitation3() const {return fExcitation3;}
- double GetExcitation4() const {return fExcitation4;}
- double GetQValue() const {return fQValue;}
- Nucleus* GetNucleus1() const {return fNuclei1;}
- Nucleus* GetNucleus2() const {return fNuclei2;}
- Nucleus* GetNucleus3() const {return fNuclei3;}
- Nucleus* GetNucleus4() const {return fNuclei4;}
- TH1F* GetCrossSectionHist() const {return fCrossSectionHist;}
- int GetVerboseLevel() const {return fVerboseLevel;}
- bool GetShoot3() const {return fshoot3;}
- bool GetShoot4() const {return fshoot4;}
-
-
- public:
+ void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
+ void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
+ void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
+ void SetCrossSectionHist (TH1F* CrossSectionHist) {delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;}
+
+ double GetBeamEnergy() const {return fBeamEnergy;}
+ double GetThetaCM() const {return fThetaCM;}
+ double GetExcitation3() const {return fExcitation3;}
+ double GetExcitation4() const {return fExcitation4;}
+ double GetQValue() const {return fQValue;}
+ Nucleus* GetNucleus1() const {return fNuclei1;}
+ Nucleus* GetNucleus2() const {return fNuclei2;}
+ Nucleus* GetNucleus3() const {return fNuclei3;}
+ Nucleus* GetNucleus4() const {return fNuclei4;}
+ TH1F* GetCrossSectionHist() const {return fCrossSectionHist;}
+ int GetVerboseLevel() const {return fVerboseLevel;}
+ bool GetShoot3() const {return fshoot3;}
+ bool GetShoot4() const {return fshoot4;}
+
+
+ public:
// Modify the CS histo to so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
- void SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax);
-
- private: // intern precompute variable
- void initializePrecomputeVariable();
- double m1;
- double m2;
- double m3;
- double m4;
-
+ void SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax);
+
+ private: // intern precompute variable
+ void initializePrecomputeVariable();
+ double m1;
+ double m2;
+ double m3;
+ double m4;
+
// Lorents Vector
- TLorentzVector fEnergyImpulsionLab_1;
- TLorentzVector fEnergyImpulsionLab_2;
- TLorentzVector fEnergyImpulsionLab_3;
- TLorentzVector fEnergyImpulsionLab_4;
- TLorentzVector fTotalEnergyImpulsionLab;
-
- TLorentzVector fEnergyImpulsionCM_1;
- TLorentzVector fEnergyImpulsionCM_2;
- TLorentzVector fEnergyImpulsionCM_3;
- TLorentzVector fEnergyImpulsionCM_4;
- TLorentzVector fTotalEnergyImpulsionCM;
-
+ TLorentzVector fEnergyImpulsionLab_1;
+ TLorentzVector fEnergyImpulsionLab_2;
+ TLorentzVector fEnergyImpulsionLab_3;
+ TLorentzVector fEnergyImpulsionLab_4;
+ TLorentzVector fTotalEnergyImpulsionLab;
+
+ TLorentzVector fEnergyImpulsionCM_1;
+ TLorentzVector fEnergyImpulsionCM_2;
+ TLorentzVector fEnergyImpulsionCM_3;
+ TLorentzVector fEnergyImpulsionCM_4;
+ TLorentzVector fTotalEnergyImpulsionCM;
+
// Impulsion Vector3
- TVector3 fImpulsionLab_1;
- TVector3 fImpulsionLab_2;
- TVector3 fImpulsionLab_3;
- TVector3 fImpulsionLab_4;
-
- TVector3 fImpulsionCM_1;
- TVector3 fImpulsionCM_2;
- TVector3 fImpulsionCM_3;
- TVector3 fImpulsionCM_4;
-
+ TVector3 fImpulsionLab_1;
+ TVector3 fImpulsionLab_2;
+ TVector3 fImpulsionLab_3;
+ TVector3 fImpulsionLab_4;
+
+ TVector3 fImpulsionCM_1;
+ TVector3 fImpulsionCM_2;
+ TVector3 fImpulsionCM_3;
+ TVector3 fImpulsionCM_4;
+
// CM Energy composante & CM impulsion norme
- Double_t ECM_1;
- Double_t ECM_2;
- Double_t ECM_3;
- Double_t ECM_4;
- Double_t pCM_1;
- Double_t pCM_2;
- Double_t pCM_3;
- Double_t pCM_4;
-
+ Double_t ECM_1;
+ Double_t ECM_2;
+ Double_t ECM_3;
+ Double_t ECM_4;
+ Double_t pCM_1;
+ Double_t pCM_2;
+ Double_t pCM_3;
+ Double_t pCM_4;
+
// Mandelstam variable
- Double_t s;
-
+ Double_t s;
+
// Center of Mass Kinematic
- Double_t BetaCM;
-
- public:
- TLorentzVector GetEnergyImpulsionLab_1() const {return fEnergyImpulsionLab_1;}
- TLorentzVector GetEnergyImpulsionLab_2() const {return fEnergyImpulsionLab_2;}
- TLorentzVector GetEnergyImpulsionLab_3() const {return fEnergyImpulsionLab_3;}
- TLorentzVector GetEnergyImpulsionLab_4() const {return fEnergyImpulsionLab_4;}
-
-
- public: // Kinematics
- // Check that the reaction is alowed
- bool CheckKinematic();
-
+ Double_t BetaCM;
+
+ public:
+ TLorentzVector GetEnergyImpulsionLab_1() const {return fEnergyImpulsionLab_1;}
+ TLorentzVector GetEnergyImpulsionLab_2() const {return fEnergyImpulsionLab_2;}
+ TLorentzVector GetEnergyImpulsionLab_3() const {return fEnergyImpulsionLab_3;}
+ TLorentzVector GetEnergyImpulsionLab_4() const {return fEnergyImpulsionLab_4;}
+
+
+ public: // Kinematics
+ // Check that the reaction is alowed
+ bool CheckKinematic();
+
// Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
- double ShootRandomThetaCM();
-
+ double ShootRandomThetaCM();
+
+ // Use fExcitationEnergyHist to shoot a Random Excitation energy and set fExcitation4 to this value
+ void ShootRandomExcitationEnergy();
+
// Compute ThetaLab and EnergyLab for product of reaction
- void KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
- double &ThetaLab4, double &KineticEnergyLab4);
-
+ void KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
+ double &ThetaLab4, double &KineticEnergyLab4);
+
// Return Excitation Energy
- double ReconstructRelativistic(double EnergyLab, double ThetaLab);
-
+ double ReconstructRelativistic(double EnergyLab, double ThetaLab);
+
// Return ThetaCM
// EnergyLab: energy measured in the laboratory frame
// ExcitationEnergy: excitation energy previously calculated.
- double EnergyLabToThetaCM(double EnergyLab, double ThetaLab);
-
- void SetNuclei3(double EnergyLab, double ThetaLab);
-
- TGraph* GetKinematicLine3(double AngleStep_CM=1);
- TGraph* GetKinematicLine4(double AngleStep_CM=1);
- TGraph* GetBrhoLine3(double AngleStep_CM=1);
- TGraph* GetThetaLabVersusThetaCM(double AngleStep_CM=1);
- TGraph* GetTheta3VsTheta4(double AngleStep_CM=1);
- void PrintKinematic();
- double GetP_CM_1() {return pCM_1;}
- double GetP_CM_2() {return pCM_2;}
- double GetP_CM_3() {return pCM_3;}
- double GetP_CM_4() {return pCM_4;}
- double GetE_CM_1() {return ECM_1;}
- double GetE_CM_2() {return ECM_2;}
- double GetE_CM_3() {return ECM_3;}
- double GetE_CM_4() {return ECM_4;}
-
+ double EnergyLabToThetaCM(double EnergyLab, double ThetaLab);
+
+ void SetNuclei3(double EnergyLab, double ThetaLab);
+
+ TGraph* GetKinematicLine3(double AngleStep_CM=1);
+ TGraph* GetKinematicLine4(double AngleStep_CM=1);
+ TGraph* GetBrhoLine3(double AngleStep_CM=1);
+ TGraph* GetThetaLabVersusThetaCM(double AngleStep_CM=1);
+ TGraph* GetTheta3VsTheta4(double AngleStep_CM=1);
+ void PrintKinematic();
+ double GetP_CM_1() {return pCM_1;}
+ double GetP_CM_2() {return pCM_2;}
+ double GetP_CM_3() {return pCM_3;}
+ double GetP_CM_4() {return pCM_4;}
+ double GetE_CM_1() {return ECM_1;}
+ double GetE_CM_2() {return ECM_2;}
+ double GetE_CM_3() {return ECM_3;}
+ double GetE_CM_4() {return ECM_4;}
+
// Print private paremeter
- void Print() const;
-
-ClassDef(Reaction,0)
+ void Print() const;
+
+ ClassDef(Reaction,0)
};
}
diff --git a/NPSimulation/GASPARD/GaspardTrackerSquare.cc b/NPSimulation/GASPARD/GaspardTrackerSquare.cc
index dac055ae8..9e07060a1 100644
--- a/NPSimulation/GASPARD/GaspardTrackerSquare.cc
+++ b/NPSimulation/GASPARD/GaspardTrackerSquare.cc
@@ -234,8 +234,10 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber,
new G4PVPlacement(G4Transform3D(*MMrot, MMpos), logicGPDSquare, Name, world, false, 0);
- logicGPDSquare->SetVisAttributes(G4VisAttributes::Invisible);
- if (m_non_sensitive_part_visiualisation) logicGPDSquare->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
+
+ G4VisAttributes* SquareVisAtt = new G4VisAttributes(G4Colour(0.90, 0.90, 0.90));
+ SquareVisAtt->SetForceWireframe(true);
+ logicGPDSquare->SetVisAttributes(SquareVisAtt);
//Place two marker to identify the u and v axis on silicon face:
//marker are placed a bit before the silicon itself so they don't perturbate simulation
@@ -276,7 +278,7 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber,
logicFirstStage->SetSensitiveDetector(m_FirstStageScorer);
///Visualisation of FirstStage Strip
- G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)); // blue
+ G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); // blue
logicFirstStage->SetVisAttributes(FirstStageVisAtt);
}
@@ -296,7 +298,7 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber,
logicSecondStage->SetSensitiveDetector(m_SecondStageScorer);
///Visualisation of SecondStage Strip
- G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
+ G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)) ;
logicSecondStage->SetVisAttributes(SecondStageVisAtt) ;
}
@@ -316,7 +318,7 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber,
logicThirdStage->SetSensitiveDetector(m_ThirdStageScorer);
///Visualisation of Third Stage
- G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.9, 0.0)); // green
+ G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); // green
logicThirdStage->SetVisAttributes(ThirdStageVisAtt);
}
}
diff --git a/NPSimulation/GASPARD/GaspardTrackerTrapezoid.cc b/NPSimulation/GASPARD/GaspardTrackerTrapezoid.cc
index 5094093f4..790830ce3 100644
--- a/NPSimulation/GASPARD/GaspardTrackerTrapezoid.cc
+++ b/NPSimulation/GASPARD/GaspardTrackerTrapezoid.cc
@@ -235,8 +235,9 @@ void GaspardTrackerTrapezoid::VolumeMaker(G4int TelescopeNumber ,
new G4PVPlacement(G4Transform3D(*MMrot, MMpos), logicGPDTrapezoid, Name, world, false, 0);
- logicGPDTrapezoid->SetVisAttributes(G4VisAttributes::Invisible);
- if (m_non_sensitive_part_visiualisation) logicGPDTrapezoid->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
+ G4VisAttributes* TrapezoideVisAtt = new G4VisAttributes(G4Colour(0.90, 0.90, 0.90));
+ TrapezoideVisAtt->SetForceWireframe(true);
+ logicGPDTrapezoid->SetVisAttributes(TrapezoideVisAtt);
//Place two marker to identify the u and v axis on silicon face:
//marker are placed a bit before the silicon itself so they don't perturbate simulation
@@ -286,7 +287,7 @@ void GaspardTrackerTrapezoid::VolumeMaker(G4int TelescopeNumber ,
logicFirstStage->SetSensitiveDetector(m_FirstStageScorer);
///Visualisation of FirstStage Strip
- G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)); // blue
+ G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); // blue
logicFirstStage->SetVisAttributes(FirstStageVisAtt);
}
@@ -315,7 +316,7 @@ void GaspardTrackerTrapezoid::VolumeMaker(G4int TelescopeNumber ,
logicSecondStage->SetSensitiveDetector(m_SecondStageScorer);
///Visualisation of SecondStage Strip
- G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5));
+ G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
logicSecondStage->SetVisAttributes(SecondStageVisAtt);
}
@@ -344,7 +345,7 @@ void GaspardTrackerTrapezoid::VolumeMaker(G4int TelescopeNumber ,
logicThirdStage->SetSensitiveDetector(m_ThirdStageScorer);
///Visualisation of Third Stage
- G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.9, 0.0)); // red
+ G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); // red
logicThirdStage->SetVisAttributes(ThirdStageVisAtt);
}
}
diff --git a/NPSimulation/MUST2/MUST2Array.cc b/NPSimulation/MUST2/MUST2Array.cc
index d611e88b4..834c272d9 100644
--- a/NPSimulation/MUST2/MUST2Array.cc
+++ b/NPSimulation/MUST2/MUST2Array.cc
@@ -59,492 +59,490 @@ using namespace MUST2 ;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// MUST2Array Specific Method
-MUST2Array::MUST2Array()
-{
- m_Event = new TMust2Data() ;
- InitializeMaterial();
+MUST2Array::MUST2Array(){
+ m_Event = new TMust2Data() ;
+ InitializeMaterial();
}
-MUST2Array::~MUST2Array()
-{
- delete m_MaterialAluminium;
- delete m_MaterialIron;
- delete m_MaterialCsI;
- delete m_MaterialVacuum ;
- delete m_MaterialMyl;
- delete m_MaterialHarvar;
+MUST2Array::~MUST2Array(){
+ delete m_MaterialAluminium;
+ delete m_MaterialIron;
+ delete m_MaterialCsI;
+ delete m_MaterialVacuum ;
+ delete m_MaterialMyl;
+ delete m_MaterialHarvar;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MUST2Array::AddTelescope( G4ThreeVector X1_Y1,
- G4ThreeVector X128_Y1,
- G4ThreeVector X1_Y128,
- G4ThreeVector X128_Y128,
- bool wSi,
- bool wSiLi,
- bool wCsI)
+ G4ThreeVector X128_Y1,
+ G4ThreeVector X1_Y128,
+ G4ThreeVector X128_Y128,
+ bool wSi,
+ bool wSiLi,
+ bool wCsI)
{
- m_DefinitionType.push_back(true);
-
- m_X1_Y1 .push_back(X1_Y1);
- m_X128_Y1 .push_back(X128_Y1);
- m_X1_Y128 .push_back(X1_Y128);
- m_X128_Y128 .push_back(X128_Y128);
- m_wSi.push_back(wSi);
- m_wSiLi.push_back(wSiLi);
- m_wCsI.push_back(wCsI);
-
- m_R.push_back(0);
- m_Theta.push_back(0);
- m_Phi.push_back(0);
- m_beta_u.push_back(0);
- m_beta_v.push_back(0);
- m_beta_w.push_back(0);
+ m_DefinitionType.push_back(true);
+
+ m_X1_Y1 .push_back(X1_Y1);
+ m_X128_Y1 .push_back(X128_Y1);
+ m_X1_Y128 .push_back(X1_Y128);
+ m_X128_Y128 .push_back(X128_Y128);
+ m_wSi.push_back(wSi);
+ m_wSiLi.push_back(wSiLi);
+ m_wCsI.push_back(wCsI);
+
+ m_R.push_back(0);
+ m_Theta.push_back(0);
+ m_Phi.push_back(0);
+ m_beta_u.push_back(0);
+ m_beta_v.push_back(0);
+ m_beta_w.push_back(0);
}
void MUST2Array::AddTelescope( G4double R,
- G4double Theta,
- G4double Phi,
- G4double beta_u,
- G4double beta_v,
- G4double beta_w,
- bool wSi,
- bool wSiLi,
- bool wCsI)
+ G4double Theta,
+ G4double Phi,
+ G4double beta_u,
+ G4double beta_v,
+ G4double beta_w,
+ bool wSi,
+ bool wSiLi,
+ bool wCsI)
{
- G4ThreeVector empty = G4ThreeVector(0, 0, 0);
+ G4ThreeVector empty = G4ThreeVector(0, 0, 0);
- m_DefinitionType.push_back(false);
+ m_DefinitionType.push_back(false);
- m_R.push_back(R);
- m_Theta.push_back(Theta);
- m_Phi.push_back(Phi);
- m_beta_u.push_back(beta_u);
- m_beta_v.push_back(beta_v);
- m_beta_w.push_back(beta_w);
- m_wSi.push_back(wSi);
- m_wSiLi.push_back(wSiLi);
- m_wCsI.push_back(wCsI);
+ m_R.push_back(R);
+ m_Theta.push_back(Theta);
+ m_Phi.push_back(Phi);
+ m_beta_u.push_back(beta_u);
+ m_beta_v.push_back(beta_v);
+ m_beta_w.push_back(beta_w);
+ m_wSi.push_back(wSi);
+ m_wSiLi.push_back(wSiLi);
+ m_wCsI.push_back(wCsI);
- m_X1_Y1 .push_back(empty) ;
- m_X128_Y1 .push_back(empty) ;
- m_X1_Y128 .push_back(empty) ;
- m_X128_Y128 .push_back(empty) ;
+ m_X1_Y1 .push_back(empty) ;
+ m_X128_Y1 .push_back(empty) ;
+ m_X1_Y128 .push_back(empty) ;
+ m_X128_Y128 .push_back(empty) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MUST2Array::VolumeMaker( G4int TelescopeNumber,
- G4ThreeVector MMpos,
- G4RotationMatrix* MMrot,
- bool wSi,
- bool wSiLi,
- bool wCsI,
- G4LogicalVolume* world)
+ G4ThreeVector MMpos,
+ G4RotationMatrix* MMrot,
+ bool wSi,
+ bool wSiLi,
+ bool wCsI,
+ G4LogicalVolume* world)
{
- G4double NbrTelescopes = TelescopeNumber;
- G4String DetectorNumber;
- std::ostringstream Number;
- Number << NbrTelescopes ;
- DetectorNumber = Number.str();
+ G4double NbrTelescopes = TelescopeNumber;
+ G4String DetectorNumber;
+ std::ostringstream Number;
+ Number << NbrTelescopes ;
+ DetectorNumber = Number.str();
+
+ ////////////////////////////////////////////////////////////////
+ ////////////// Starting Volume Definition //////////////////////
+ ////////////////////////////////////////////////////////////////
+ G4Trd* solidMM = new G4Trd("MUST2Telescope" + DetectorNumber, 0.5*FaceFront, 0.5*FaceBack, 0.5*FaceFront, 0.5*FaceBack, 0.5*Length);
+ G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, m_MaterialIron, "MUST2Telescope" + DetectorNumber, 0, 0, 0);
+ G4String Name = "MUST2Telescope" + DetectorNumber ;
-////////////////////////////////////////////////////////////////
-////////////// Starting Volume Definition //////////////////////
-////////////////////////////////////////////////////////////////
+ new G4PVPlacement( G4Transform3D(*MMrot, MMpos),
+ logicMM ,
+ Name ,
+ world ,
+ false ,
+ 0 );
+ if (m_non_sensitive_part_visiualisation){
+ G4VisAttributes* FrameVisAtt = new G4VisAttributes(G4Colour(0.80, 0.80, 0.80));
+ FrameVisAtt->SetForceWireframe(true);
+ logicMM->SetVisAttributes(FrameVisAtt) ;
+ }
+ else logicMM->SetVisAttributes(G4VisAttributes::Invisible) ;
- G4Trd* solidMM = new G4Trd("MUST2Telescope" + DetectorNumber, 0.5*FaceFront, 0.5*FaceBack, 0.5*FaceFront, 0.5*FaceBack, 0.5*Length);
- G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, m_MaterialIron, "MUST2Telescope" + DetectorNumber, 0, 0, 0);
+ G4ThreeVector positionVacBox = G4ThreeVector(0, 0, VacBox_PosZ);
- G4String Name = "MUST2Telescope" + DetectorNumber ;
+ G4Trd* solidVacBox = new G4Trd("solidVacBox", 0.5*SiliconFace, 0.5*CsIFaceFront, 0.5*SiliconFace, 0.5*CsIFaceFront, 0.5*VacBoxThickness);
+ G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0);
- new G4PVPlacement( G4Transform3D(*MMrot, MMpos) ,
- logicMM ,
- Name ,
- world ,
- false ,
- 0 );
+ new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, 0);
+ logicVacBox->SetVisAttributes(G4VisAttributes::Invisible);
+
+ G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)) ;
+ ////////////////////////////////////////////////////////////////
+ /////////////////Si Strip Construction//////////////////////////
+ ////////////////////////////////////////////////////////////////
- if (m_non_sensitive_part_visiualisation) logicMM->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90))) ;
- else logicMM->SetVisAttributes(G4VisAttributes::Invisible) ;
+ if (wSi) {
+ G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, AluStripFront_PosZ);
+ G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripBack_PosZ);
- G4ThreeVector positionVacBox = G4ThreeVector(0, 0, VacBox_PosZ);
+ G4Box* solidAluStrip = new G4Box("AluBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*AluStripThickness);
+ G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0);
- G4Trd* solidVacBox = new G4Trd("solidVacBox", 0.5*SiliconFace, 0.5*CsIFaceFront, 0.5*SiliconFace, 0.5*CsIFaceFront, 0.5*VacBoxThickness);
- G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0);
+ new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicMM, false, 0);
+ new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicMM, false, 0);
- new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, 0);
+ logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible);
- logicVacBox->SetVisAttributes(G4VisAttributes::Invisible);
+ G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ);
-////////////////////////////////////////////////////////////////
-/////////////////Si Strip Construction//////////////////////////
-////////////////////////////////////////////////////////////////
+ G4Box* solidSilicon = new G4Box("solidSilicon", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconThickness);
+ G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0);
- if (wSi) {
- G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, AluStripFront_PosZ);
- G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripBack_PosZ);
+ new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicMM, false, 0);
- G4Box* solidAluStrip = new G4Box("AluBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*AluStripThickness);
- G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0);
- new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicMM, false, 0);
- new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicMM, false, 0);
+ ///Set Silicon strip sensible
+ logicSilicon->SetSensitiveDetector(m_StripScorer);
+
+ ///Visualisation of Silicon Strip
+ logicSilicon->SetVisAttributes(SiliconVisAtt) ;
+ }
+
+ ////////////////////////////////////////////////////////////////
+ //////////////////// SiLi Construction ////////////////////////
+ ////////////////////////////////////////////////////////////////
- logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible);
+ if (wSiLi) {
+ G4double SiLiSpace = 8 * mm;
+ G4RotationMatrix* rotSiLi = new G4RotationMatrix(0,0,0);
+ G4Box* solidSiLi = new G4Box("SiLi", 0.5*SiliconFace+0.5*SiLiSpace, 0.5*SiliconFace, 0.5*SiLiThickness);
+ G4LogicalVolume* logicSiLi = new G4LogicalVolume(solidSiLi, m_MaterialAluminium, Name + "_SiLi" , 0, 0, 0);
- G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ);
+ logicSiLi->SetVisAttributes(G4VisAttributes::Invisible);
- G4Box* solidSilicon = new G4Box("solidSilicon", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconThickness);
- G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0);
+ new G4PVPlacement( G4Transform3D(*rotSiLi, G4ThreeVector(0,0,0) ) ,
+ logicSiLi ,
+ Name + "_SiLi",
+ logicVacBox ,
+ false ,
+ 0);
- new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicMM, false, 0);
+ // SiLi are placed inside of the VacBox...
+ // Left/Right define when looking to detector from Si to CsI
+ G4double SiLi_HighY_Upper = 19.86 * mm;
+ G4double SiLi_HighY_Center = 25.39 * mm;
+ G4double SiLi_WidthX_Left = 22.85 * mm;
+ G4double SiLi_WidthX_Right = 24.9 * mm;
+ G4double SiLi_ShiftX = 0.775 * mm;
+ // SiLi are organized by two group of 8 Up(9 to 15) and Down(1 to 8).
+ G4ThreeVector ShiftSiLiUp = G4ThreeVector(-0.25 * SiliconFace - 0.5 * SiLiSpace, 0, 0) ;
+ G4ThreeVector ShiftSiLiDown = G4ThreeVector(0.25 * SiliconFace + 0.5 * SiLiSpace, 0, 0) ;
- ///Set Silicon strip sensible
- logicSilicon->SetSensitiveDetector(m_StripScorer);
+ // SiLi : left side of SiLi detector
+ G4Box* solidSiLi_LT = new G4Box("SiLi_LT" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
+ G4Box* solidSiLi_RT = new G4Box("SiLi_RT" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
+ G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
+ G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
+ G4Box* solidSiLi_LB = new G4Box("SiLi_LB" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
+ G4Box* solidSiLi_RB = new G4Box("SiLi_RB" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
+ G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
+ G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
- ///Visualisation of Silicon Strip
- G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ;
- logicSilicon->SetVisAttributes(SiliconVisAtt) ;
- }
+ G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume(solidSiLi_LT , m_MaterialSilicon , "SiLi_LT" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume(solidSiLi_RT , m_MaterialSilicon , "SiLi_RT" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume(solidSiLi_LC1 , m_MaterialSilicon , "SiLi_LC1" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume(solidSiLi_RC1 , m_MaterialSilicon , "SiLi_RC1" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume(solidSiLi_LB , m_MaterialSilicon , "SiLi_LB" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume(solidSiLi_RB , m_MaterialSilicon , "SiLi_RB" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume(solidSiLi_LC2 , m_MaterialSilicon , "SiLi_LC2" , 0 , 0 , 0);
+ G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume(solidSiLi_RC2 , m_MaterialSilicon , "SiLi_RC2" , 0 , 0 , 0);
-////////////////////////////////////////////////////////////////
-//////////////////// SiLi Construction ////////////////////////
-////////////////////////////////////////////////////////////////
+ G4double interSiLi = 0.5 * mm;
- if (wSiLi) {
- G4double SiLiSpace = 8 * mm;
+ // Top
+ G4ThreeVector positionSiLi_LT_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
+ 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
+ 0);
- G4RotationMatrix* rotSiLi = new G4RotationMatrix(0,0,0);
+ positionSiLi_LT_up += ShiftSiLiUp ;
- // G4Box* solidSiLi = new G4Box("SiLi", 0.5*SiLiFaceX, 0.5*SiLiFaceY, 0.5*SiLiThickness);
-
- G4Box* solidSiLi = new G4Box("SiLi", 0.5*SiliconFace+0.5*SiLiSpace, 0.5*SiliconFace, 0.5*SiLiThickness);
-
- G4LogicalVolume* logicSiLi = new G4LogicalVolume(solidSiLi, m_MaterialAluminium, Name + "_SiLi" , 0, 0, 0);
-
- new G4PVPlacement( G4Transform3D(*rotSiLi, G4ThreeVector(0,0,0) ) ,
- logicSiLi ,
- Name + "_SiLi",
- logicVacBox ,
- false ,
- 0);
-
- // SiLi are placed inside of the VacBox...
- // Left/Right define when looking to detector from Si to CsI
- G4double SiLi_HighY_Upper = 19.86 * mm;
- G4double SiLi_HighY_Center = 25.39 * mm;
- G4double SiLi_WidthX_Left = 22.85 * mm;
- G4double SiLi_WidthX_Right = 24.9 * mm;
- G4double SiLi_ShiftX = 0.775 * mm;
-
- // SiLi are organized by two group of 8 Up(9 to 15) and Down(1 to 8).
- G4ThreeVector ShiftSiLiUp = G4ThreeVector(-0.25 * SiliconFace - 0.5 * SiLiSpace, 0, 0) ;
- G4ThreeVector ShiftSiLiDown = G4ThreeVector(0.25 * SiliconFace + 0.5 * SiLiSpace, 0, 0) ;
-
- // SiLi : left side of SiLi detector
- G4Box* solidSiLi_LT = new G4Box("SiLi_LT" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
- G4Box* solidSiLi_RT = new G4Box("SiLi_RT" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
- G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
- G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
- G4Box* solidSiLi_LB = new G4Box("SiLi_LB" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
- G4Box* solidSiLi_RB = new G4Box("SiLi_RB" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness);
- G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
- G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness);
-
- G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume(solidSiLi_LT , m_MaterialSilicon , "SiLi_LT" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume(solidSiLi_RT , m_MaterialSilicon , "SiLi_RT" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume(solidSiLi_LC1 , m_MaterialSilicon , "SiLi_LC1" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume(solidSiLi_RC1 , m_MaterialSilicon , "SiLi_RC1" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume(solidSiLi_LB , m_MaterialSilicon , "SiLi_LB" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume(solidSiLi_RB , m_MaterialSilicon , "SiLi_RB" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume(solidSiLi_LC2 , m_MaterialSilicon , "SiLi_LC2" , 0 , 0 , 0);
- G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume(solidSiLi_RC2 , m_MaterialSilicon , "SiLi_RC2" , 0 , 0 , 0);
-
- G4double interSiLi = 0.5 * mm;
-
- // Top
- G4ThreeVector positionSiLi_LT_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
- 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
- 0);
-
- positionSiLi_LT_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_RT_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
- 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
- 0);
-
- positionSiLi_RT_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_LC1_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
- 0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
- 0);
-
- positionSiLi_LC1_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_RC1_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
- 0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
- 0);
-
- positionSiLi_RC1_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_LB_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
- -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi ,
- 0);
-
- positionSiLi_LB_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_RB_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
- -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi ,
- 0);
-
- positionSiLi_RB_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_LC2_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
- -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
- 0);
-
- positionSiLi_LC2_up += ShiftSiLiUp ;
-
- G4ThreeVector positionSiLi_RC2_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
- -0.5 * SiLi_HighY_Center - 0.5 * interSiLi ,
- 0);
-
- positionSiLi_RC2_up += ShiftSiLiUp ;
-
-
- // Down
- G4ThreeVector positionSiLi_LT_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
- 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
- 0);
-
- positionSiLi_LT_down += ShiftSiLiDown ;
-
- G4ThreeVector positionSiLi_RT_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
- 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
- 0);
-
- positionSiLi_RT_down += ShiftSiLiDown ;
-
- G4ThreeVector positionSiLi_LC1_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
- 0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
- 0);
-
- positionSiLi_LC1_down += ShiftSiLiDown ;
-
- G4ThreeVector positionSiLi_RC1_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
- 0.5 * SiLi_HighY_Center + 0.5 * interSiLi ,
- 0);
+ G4ThreeVector positionSiLi_RT_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
+ 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
+ 0);
+
+ positionSiLi_RT_up += ShiftSiLiUp ;
+
+ G4ThreeVector positionSiLi_LC1_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
+ 0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
+ 0);
+
+ positionSiLi_LC1_up += ShiftSiLiUp ;
+
+ G4ThreeVector positionSiLi_RC1_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
+ 0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
+ 0);
+
+ positionSiLi_RC1_up += ShiftSiLiUp ;
+
+ G4ThreeVector positionSiLi_LB_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
+ -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi ,
+ 0);
+
+ positionSiLi_LB_up += ShiftSiLiUp ;
+
+ G4ThreeVector positionSiLi_RB_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
+ -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi ,
+ 0);
+
+ positionSiLi_RB_up += ShiftSiLiUp ;
+
+ G4ThreeVector positionSiLi_LC2_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
+ -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
+ 0);
+
+ positionSiLi_LC2_up += ShiftSiLiUp ;
+
+ G4ThreeVector positionSiLi_RC2_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
+ -0.5 * SiLi_HighY_Center - 0.5 * interSiLi ,
+ 0);
+
+ positionSiLi_RC2_up += ShiftSiLiUp ;
+
+
+ // Down
+ G4ThreeVector positionSiLi_LT_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
+ 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
+ 0);
+
+ positionSiLi_LT_down += ShiftSiLiDown ;
- positionSiLi_RC1_down += ShiftSiLiDown ;
+ G4ThreeVector positionSiLi_RT_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
+ 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi,
+ 0);
- G4ThreeVector positionSiLi_LB_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
- -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi,
- 0);
+ positionSiLi_RT_down += ShiftSiLiDown ;
- positionSiLi_LB_down += ShiftSiLiDown ;
+ G4ThreeVector positionSiLi_LC1_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
+ 0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
+ 0);
- G4ThreeVector positionSiLi_RB_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
- -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi,
- 0);
+ positionSiLi_LC1_down += ShiftSiLiDown ;
- positionSiLi_RB_down += ShiftSiLiDown ;
+ G4ThreeVector positionSiLi_RC1_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
+ 0.5 * SiLi_HighY_Center + 0.5 * interSiLi ,
+ 0);
- G4ThreeVector positionSiLi_LC2_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
- -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
- 0);
+ positionSiLi_RC1_down += ShiftSiLiDown ;
- positionSiLi_LC2_down += ShiftSiLiDown ;
+ G4ThreeVector positionSiLi_LB_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
+ -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi,
+ 0);
- G4ThreeVector positionSiLi_RC2_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
- -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
- 0);
-
- positionSiLi_RC2_down += ShiftSiLiDown ;
-
-
- // up
- new G4PVPlacement(0 , positionSiLi_LT_up , logicSiLi_LT , Name + "_SiLi_Pad9" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_RT_up , logicSiLi_RT , Name + "_SiLi_Pad10" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" , logicSiLi , false , 0);
- new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" , logicSiLi , false , 0);
+ positionSiLi_LB_down += ShiftSiLiDown ;
- new G4PVPlacement(0 , positionSiLi_LB_up , logicSiLi_LB , Name + "_SiLi_Pad16" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_RB_up , logicSiLi_RB , Name + "_SiLi_Pad15" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" , logicSiLi , false , 0);
- new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" , logicSiLi , false , 0);
+ G4ThreeVector positionSiLi_RB_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
+ -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi,
+ 0);
-
- // down
- new G4PVPlacement(0 , positionSiLi_LT_down , logicSiLi_LT , Name + "_SiLi_Pad2" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_RT_down , logicSiLi_RT , Name + "_SiLi_Pad1" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" , logicSiLi , false , 0) ;
+ positionSiLi_RB_down += ShiftSiLiDown ;
- new G4PVPlacement(0 , positionSiLi_LB_down , logicSiLi_LB , Name + "_SiLi_Pad7" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_RB_down , logicSiLi_RB , Name + "_SiLi_Pad8" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" , logicSiLi , false , 0) ;
- new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" , logicSiLi , false , 0) ;
+ G4ThreeVector positionSiLi_LC2_down = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
+ -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
+ 0);
+ positionSiLi_LC2_down += ShiftSiLiDown ;
+ G4ThreeVector positionSiLi_RC2_down = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
+ -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
+ 0);
- logicSiLi->SetVisAttributes(G4VisAttributes(G4Colour(1, 1., 1.)));
+ positionSiLi_RC2_down += ShiftSiLiDown ;
- // Set SiLi sensible
- logicSiLi_LT->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_RT->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_LC1->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_RC1->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_LB->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_RB->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_LC2->SetSensitiveDetector(m_SiLiScorer);
- logicSiLi_RC2->SetSensitiveDetector(m_SiLiScorer);
+ // up
+ new G4PVPlacement(0 , positionSiLi_LT_up , logicSiLi_LT , Name + "_SiLi_Pad9" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_RT_up , logicSiLi_RT , Name + "_SiLi_Pad10" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" , logicSiLi , false , 0);
+ new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" , logicSiLi , false , 0);
- // Mark blue a SiLi to see telescope orientation
- logicSiLi_LT->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
- logicSiLi_RT->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
- logicSiLi_LC1->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
- logicSiLi_RC1->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
+ new G4PVPlacement(0 , positionSiLi_LB_up , logicSiLi_LB , Name + "_SiLi_Pad16" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_RB_up , logicSiLi_RB , Name + "_SiLi_Pad15" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" , logicSiLi , false , 0);
+ new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" , logicSiLi , false , 0);
- logicSiLi_LB->SetVisAttributes(G4VisAttributes(G4Colour(0, 0, 1.)));
- logicSiLi_RB->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
- logicSiLi_LC2->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
- logicSiLi_RC2->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0)));
-
-
- delete rotSiLi;
- }
-////////////////////////////////////////////////////////////////
-//////////////////// CsI Construction//////////////////////////
-////////////////////////////////////////////////////////////////
+ // down
+ new G4PVPlacement(0 , positionSiLi_LT_down , logicSiLi_LT , Name + "_SiLi_Pad2" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_RT_down , logicSiLi_RT , Name + "_SiLi_Pad1" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" , logicSiLi , false , 0) ;
- if (wCsI) {
-
- G4ThreeVector positionCsI = G4ThreeVector(0, 0, CsI_PosZ);
- G4Trd* solidCsI = new G4Trd("csI", 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIThickness);
+ new G4PVPlacement(0 , positionSiLi_LB_down , logicSiLi_LB , Name + "_SiLi_Pad7" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_RB_down , logicSiLi_RB , Name + "_SiLi_Pad8" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" , logicSiLi , false , 0) ;
+ new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" , logicSiLi , false , 0) ;
- G4LogicalVolume* logicCsI = new G4LogicalVolume(solidCsI, m_MaterialAluminium, Name + "_CsI_Mylar", 0, 0, 0);
- new G4PVPlacement(0, positionCsI, logicCsI, Name + "_CsI_Mylar", logicMM, false, 0);
+ // Set SiLi sensible
+ logicSiLi_LT->SetSensitiveDetector(m_SiLiScorer);
+ logicSiLi_RT->SetSensitiveDetector(m_SiLiScorer);
+ logicSiLi_LC1->SetSensitiveDetector(m_SiLiScorer);
+ logicSiLi_RC1->SetSensitiveDetector(m_SiLiScorer);
- G4ThreeVector positionMylarCsI = G4ThreeVector(0, 0, MylarCsIThickness * 0.5 - CsIThickness * 0.5);
+ logicSiLi_LB->SetSensitiveDetector(m_SiLiScorer);
+ logicSiLi_RB->SetSensitiveDetector(m_SiLiScorer);
+ logicSiLi_LC2->SetSensitiveDetector(m_SiLiScorer);
+ logicSiLi_RC2->SetSensitiveDetector(m_SiLiScorer);
- G4Box* solidMylarCsI = new G4Box("MylarCsIBox", 0.5*CsIFaceFront, 0.5*CsIFaceFront, 0.5*MylarCsIThickness);
- G4LogicalVolume* logicMylarCsI = new G4LogicalVolume(solidMylarCsI, m_MaterialMyl, Name + "_CsI_Mylar", 0, 0, 0);
+ // Mark blue a SiLi to see telescope orientation
+ G4VisAttributes* SiLiVisAtt = new G4VisAttributes(G4Colour(0.3, 1, 0.3));
- new G4PVPlacement(0, positionMylarCsI, logicMylarCsI, Name + "_CsI_Mylar", logicCsI, false, 0);
+ logicSiLi_LT->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_RT->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_LC1->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_RC1->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_LB->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_RB->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_LC2->SetVisAttributes(SiLiVisAtt);
+ logicSiLi_RC2->SetVisAttributes(SiLiVisAtt);
- logicCsI->SetVisAttributes(G4VisAttributes::Invisible);
- logicMylarCsI->SetVisAttributes(G4VisAttributes::Invisible);
- // Cristal1
- G4Trap* solidCristal1 = new G4Trap("Cristal1", 40.*mm / 2., 6.693896*deg, 41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal1 = new G4LogicalVolume(solidCristal1, m_MaterialCsI, Name + "_CsI_Cristal1", 0, 0, 0);
+ delete rotSiLi;
+ }
- // Cristal2
- G4Trap* solidCristal2 = new G4Trap("Cristal2", 40.*mm / 2., 17.8836*deg, (74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal2 = new G4LogicalVolume(solidCristal2, m_MaterialCsI, Name + "_CsI_Cristal2", 0, 0, 0);
+ ////////////////////////////////////////////////////////////////
+ //////////////////// CsI Construction//////////////////////////
+ ////////////////////////////////////////////////////////////////
- // Cristal3
- G4Trap* solidCristal3 = new G4Trap("Cristal3", 40.*mm / 2., 18.243*deg, 13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal3 = new G4LogicalVolume(solidCristal3, m_MaterialCsI, Name + "_CsI_Cristal3", 0, 0, 0);
+ if (wCsI) {
- // Cristal4
+ G4ThreeVector positionCsI = G4ThreeVector(0, 0, CsI_PosZ);
+ G4Trd* solidCsI = new G4Trd("csI", 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIThickness);
- G4Trap* solidCristal4 = new G4Trap("Cristal4", 40.*mm / 2., 24.0482*deg, 44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal4 = new G4LogicalVolume(solidCristal4, m_MaterialCsI, Name + "_CsI_Cristal4", 0, 0, 0);
+ G4LogicalVolume* logicCsI = new G4LogicalVolume(solidCsI, m_MaterialAluminium, Name + "_CsI_Mylar", 0, 0, 0);
+ new G4PVPlacement(0, positionCsI, logicCsI, Name + "_CsI_Mylar", logicMM, false, 0);
+ G4ThreeVector positionMylarCsI = G4ThreeVector(0, 0, MylarCsIThickness * 0.5 - CsIThickness * 0.5);
- // Cristal1s
+ G4Box* solidMylarCsI = new G4Box("MylarCsIBox", 0.5*CsIFaceFront, 0.5*CsIFaceFront, 0.5*MylarCsIThickness);
+ G4LogicalVolume* logicMylarCsI = new G4LogicalVolume(solidMylarCsI, m_MaterialMyl, Name + "_CsI_Mylar", 0, 0, 0);
- G4Trap* solidCristal1s = new G4Trap("Cristal1s", 40.*mm / 2., 6.693896*deg, -41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal1s = new G4LogicalVolume(solidCristal1s, m_MaterialCsI, Name + "_CsI_Cristal1s", 0, 0, 0);
+ new G4PVPlacement(0, positionMylarCsI, logicMylarCsI, Name + "_CsI_Mylar", logicCsI, false, 0);
- // Cristal2s
- G4Trap* solidCristal2s = new G4Trap("Cristal2s", 40.*mm / 2., 17.8836*deg, -(74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal2s = new G4LogicalVolume(solidCristal2s, m_MaterialCsI, Name + "_CsI_Cristal2s", 0, 0, 0);
+ logicCsI->SetVisAttributes(G4VisAttributes::Invisible);
+ logicMylarCsI->SetVisAttributes(G4VisAttributes::Invisible);
- // Cristal3s
+ // Cristal1
+ G4Trap* solidCristal1 = new G4Trap("Cristal1", 40.*mm / 2., 6.693896*deg, 41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal1 = new G4LogicalVolume(solidCristal1, m_MaterialCsI, Name + "_CsI_Cristal1", 0, 0, 0);
- G4Trap* solidCristal3s = new G4Trap("Cristal3s", 40.*mm / 2., 18.243*deg, -13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal3s = new G4LogicalVolume(solidCristal3s, m_MaterialCsI, Name + "_CsI_Cristal3s", 0, 0, 0);
+ // Cristal2
+ G4Trap* solidCristal2 = new G4Trap("Cristal2", 40.*mm / 2., 17.8836*deg, (74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal2 = new G4LogicalVolume(solidCristal2, m_MaterialCsI, Name + "_CsI_Cristal2", 0, 0, 0);
- // Cristal4s
+ // Cristal3
+ G4Trap* solidCristal3 = new G4Trap("Cristal3", 40.*mm / 2., 18.243*deg, 13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal3 = new G4LogicalVolume(solidCristal3, m_MaterialCsI, Name + "_CsI_Cristal3", 0, 0, 0);
- G4Trap* solidCristal4s = new G4Trap("Cristal4s", 40.*mm / 2., 24.0482*deg, -44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg);
- G4LogicalVolume* logicCristal4s = new G4LogicalVolume(solidCristal4s, m_MaterialCsI, Name + "_CsI_Cristal4s", 0, 0, 0);
- //to see telescope orientation
- G4LogicalVolume* logicCristal4sbis = new G4LogicalVolume(solidCristal4s, m_MaterialCsI, Name + "_CsI_Cristal4s", 0, 0, 0);
+ // Cristal4
+ G4Trap* solidCristal4 = new G4Trap("Cristal4", 40.*mm / 2., 24.0482*deg, 44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal4 = new G4LogicalVolume(solidCristal4, m_MaterialCsI, Name + "_CsI_Cristal4", 0, 0, 0);
- G4double XEdge1 = 16.96 * mm + DistInterCsI * 0.5;
- G4double YEdge1 = 15.01 * mm + DistInterCsI * 0.5;
- G4double XEdge2 = 50.63 * mm + DistInterCsI * 1.5;
- G4double YEdge2 = 48.67 * mm + DistInterCsI * 1.5;
- G4ThreeVector positionCristal1 = G4ThreeVector(-XEdge1, YEdge1, 0 * mm);
- G4ThreeVector positionCristal2 = G4ThreeVector(-XEdge1, YEdge2, 0);
- G4ThreeVector positionCristal3 = G4ThreeVector(-XEdge2, YEdge1, 0);
- G4ThreeVector positionCristal4 = G4ThreeVector(-XEdge2, YEdge2, 0);
+ // Cristal1s
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1, logicCristal1, Name + "_CsI_Cristal1", logicCsI, false, 1);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2, logicCristal2, Name + "_CsI_Cristal2", logicCsI, false, 2);
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3, logicCristal3, Name + "_CsI_Cristal3", logicCsI, false, 3);
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4, logicCristal4, Name + "_CsI_Cristal4", logicCsI, false, 4);
+ G4Trap* solidCristal1s = new G4Trap("Cristal1s", 40.*mm / 2., 6.693896*deg, -41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal1s = new G4LogicalVolume(solidCristal1s, m_MaterialCsI, Name + "_CsI_Cristal1s", 0, 0, 0);
+ // Cristal2s
- G4ThreeVector positionCristal1b = G4ThreeVector(XEdge1, -YEdge1, 0 * mm);
- G4ThreeVector positionCristal2b = G4ThreeVector(XEdge1, -YEdge2, 0);
- G4ThreeVector positionCristal3b = G4ThreeVector(XEdge2, -YEdge1, 0);
- G4ThreeVector positionCristal4b = G4ThreeVector(XEdge2, -YEdge2, 0);
+ G4Trap* solidCristal2s = new G4Trap("Cristal2s", 40.*mm / 2., 17.8836*deg, -(74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal2s = new G4LogicalVolume(solidCristal2s, m_MaterialCsI, Name + "_CsI_Cristal2s", 0, 0, 0);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1b, logicCristal1, Name + "_CsI_Cristal5", logicCsI, false, 5);
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal2b, logicCristal2, Name + "_CsI_Cristal6", logicCsI, false, 6);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3b, logicCristal3, Name + "_CsI_Cristal7", logicCsI, false, 7);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4b, logicCristal4, Name + "_CsI_Cristal8", logicCsI, false, 8);
+ // Cristal3s
- G4ThreeVector positionCristal1s = G4ThreeVector(-XEdge1, -YEdge1, 0 * mm);
- G4ThreeVector positionCristal2s = G4ThreeVector(-XEdge1, -YEdge2, 0);
- G4ThreeVector positionCristal3s = G4ThreeVector(-XEdge2, -YEdge1, 0);
- G4ThreeVector positionCristal4s = G4ThreeVector(-XEdge2, -YEdge2, 0);
+ G4Trap* solidCristal3s = new G4Trap("Cristal3s", 40.*mm / 2., 18.243*deg, -13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal3s = new G4LogicalVolume(solidCristal3s, m_MaterialCsI, Name + "_CsI_Cristal3s", 0, 0, 0);
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1s, logicCristal1s, Name + "_CsI_Cristal9", logicCsI, false, 9);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2s, logicCristal2s, Name + "_CsI_Cristal10", logicCsI, false, 10);
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3s, logicCristal3s, Name + "_CsI_Cristal11", logicCsI, false, 11);
- new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4s, logicCristal4sbis, Name + "_CsI_Cristal12", logicCsI, false, 12);
+ // Cristal4s
- G4ThreeVector positionCristal1sb = G4ThreeVector(XEdge1, YEdge1, 0 * mm);
- G4ThreeVector positionCristal2sb = G4ThreeVector(XEdge1, YEdge2, 0);
- G4ThreeVector positionCristal3sb = G4ThreeVector(XEdge2, YEdge1, 0);
- G4ThreeVector positionCristal4sb = G4ThreeVector(XEdge2, YEdge2, 0);
+ G4Trap* solidCristal4s = new G4Trap("Cristal4s", 40.*mm / 2., 24.0482*deg, -44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg);
+ G4LogicalVolume* logicCristal4s = new G4LogicalVolume(solidCristal4s, m_MaterialCsI, Name + "_CsI_Cristal4s", 0, 0, 0);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1sb, logicCristal1s, Name + "_CsI_Cristal13", logicCsI, false, 13);
- new G4PVPlacement(Rotation(180, 0, 0), positionCristal2sb, logicCristal2s, Name + "_CsI_Cristal14", logicCsI, false, 14);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3sb, logicCristal3s, Name + "_CsI_Cristal15", logicCsI, false, 15);
- new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4sb, logicCristal4s, Name + "_CsI_Cristal16", logicCsI, false, 16);
+ G4double XEdge1 = 16.96 * mm + DistInterCsI * 0.5;
+ G4double YEdge1 = 15.01 * mm + DistInterCsI * 0.5;
+ G4double XEdge2 = 50.63 * mm + DistInterCsI * 1.5;
+ G4double YEdge2 = 48.67 * mm + DistInterCsI * 1.5;
- ///Set CsI sensible
- logicCristal1->SetSensitiveDetector(m_CsIScorer);
- logicCristal2->SetSensitiveDetector(m_CsIScorer);
- logicCristal3->SetSensitiveDetector(m_CsIScorer);
- logicCristal4->SetSensitiveDetector(m_CsIScorer);
+ G4ThreeVector positionCristal1 = G4ThreeVector(-XEdge1, YEdge1, 0 * mm);
+ G4ThreeVector positionCristal2 = G4ThreeVector(-XEdge1, YEdge2, 0);
+ G4ThreeVector positionCristal3 = G4ThreeVector(-XEdge2, YEdge1, 0);
+ G4ThreeVector positionCristal4 = G4ThreeVector(-XEdge2, YEdge2, 0);
- logicCristal1s->SetSensitiveDetector(m_CsIScorer);
- logicCristal2s->SetSensitiveDetector(m_CsIScorer);
- logicCristal3s->SetSensitiveDetector(m_CsIScorer);
- logicCristal4s->SetSensitiveDetector(m_CsIScorer);
- logicCristal4sbis->SetSensitiveDetector(m_CsIScorer);
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1, logicCristal1, Name + "_CsI_Cristal1", logicCsI, false, 1);
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2, logicCristal2, Name + "_CsI_Cristal2", logicCsI, false, 2);
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3, logicCristal3, Name + "_CsI_Cristal3", logicCsI, false, 3);
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4, logicCristal4, Name + "_CsI_Cristal4", logicCsI, false, 4);
- //Mark blue a CsI corner crystal to see telescope orientation
- logicCristal4sbis->SetVisAttributes(G4VisAttributes(G4Colour(0, 0, 1.)));
- }
+
+ G4ThreeVector positionCristal1b = G4ThreeVector(XEdge1, -YEdge1, 0 * mm);
+ G4ThreeVector positionCristal2b = G4ThreeVector(XEdge1, -YEdge2, 0);
+ G4ThreeVector positionCristal3b = G4ThreeVector(XEdge2, -YEdge1, 0);
+ G4ThreeVector positionCristal4b = G4ThreeVector(XEdge2, -YEdge2, 0);
+
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1b, logicCristal1, Name + "_CsI_Cristal5", logicCsI, false, 5);
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal2b, logicCristal2, Name + "_CsI_Cristal6", logicCsI, false, 6);
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3b, logicCristal3, Name + "_CsI_Cristal7", logicCsI, false, 7);
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4b, logicCristal4, Name + "_CsI_Cristal8", logicCsI, false, 8);
+
+ G4ThreeVector positionCristal1s = G4ThreeVector(-XEdge1, -YEdge1, 0 * mm);
+ G4ThreeVector positionCristal2s = G4ThreeVector(-XEdge1, -YEdge2, 0);
+ G4ThreeVector positionCristal3s = G4ThreeVector(-XEdge2, -YEdge1, 0);
+ G4ThreeVector positionCristal4s = G4ThreeVector(-XEdge2, -YEdge2, 0);
+
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1s, logicCristal1s, Name + "_CsI_Cristal9", logicCsI, false, 9);
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2s, logicCristal2s, Name + "_CsI_Cristal10", logicCsI, false, 10);
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3s, logicCristal3s, Name + "_CsI_Cristal11", logicCsI, false, 11);
+ new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4s, logicCristal4s, Name + "_CsI_Cristal12", logicCsI, false, 12);
+
+ G4ThreeVector positionCristal1sb = G4ThreeVector(XEdge1, YEdge1, 0 * mm);
+ G4ThreeVector positionCristal2sb = G4ThreeVector(XEdge1, YEdge2, 0);
+ G4ThreeVector positionCristal3sb = G4ThreeVector(XEdge2, YEdge1, 0);
+ G4ThreeVector positionCristal4sb = G4ThreeVector(XEdge2, YEdge2, 0);
+
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1sb, logicCristal1s, Name + "_CsI_Cristal13", logicCsI, false, 13);
+ new G4PVPlacement(Rotation(180, 0, 0), positionCristal2sb, logicCristal2s, Name + "_CsI_Cristal14", logicCsI, false, 14);
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3sb, logicCristal3s, Name + "_CsI_Cristal15", logicCsI, false, 15);
+ new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4sb, logicCristal4s, Name + "_CsI_Cristal16", logicCsI, false, 16);
+
+ ///Set CsI sensible
+ logicCristal1->SetSensitiveDetector(m_CsIScorer);
+ logicCristal2->SetSensitiveDetector(m_CsIScorer);
+ logicCristal3->SetSensitiveDetector(m_CsIScorer);
+ logicCristal4->SetSensitiveDetector(m_CsIScorer);
+
+ logicCristal1s->SetSensitiveDetector(m_CsIScorer);
+ logicCristal2s->SetSensitiveDetector(m_CsIScorer);
+ logicCristal3s->SetSensitiveDetector(m_CsIScorer);
+ logicCristal4s->SetSensitiveDetector(m_CsIScorer);
+
+ //Mark blue a CsI corner crystal to see telescope orientation
+ G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1, 0.5, 0));
+ logicCristal1 ->SetVisAttributes(CsIVisAtt);
+ logicCristal2 ->SetVisAttributes(CsIVisAtt);
+ logicCristal3 ->SetVisAttributes(CsIVisAtt);
+ logicCristal4 ->SetVisAttributes(CsIVisAtt);
+ logicCristal1s ->SetVisAttributes(CsIVisAtt);
+ logicCristal2s ->SetVisAttributes(CsIVisAtt);
+ logicCristal3s ->SetVisAttributes(CsIVisAtt);
+ logicCristal4s ->SetVisAttributes(CsIVisAtt);
+
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -554,872 +552,859 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber,
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
-void MUST2Array::ReadConfiguration(string Path)
-{
- ifstream ConfigFile ;
- ConfigFile.open(Path.c_str()) ;
- string LineBuffer ;
- string DataBuffer ;
-
- // A:X1_Y1 --> X:1 Y:1
- // B:X128_Y1 --> X:128 Y:1
- // C:X1_Y128 --> X:1 Y:128
- // D:X128_Y128 --> X:128 Y:128
-
- G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
- G4ThreeVector A , B , C , D ;
- G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
- int SI = 0 , SILI = 0 , CSI = 0 ;
-
- bool check_A = false ;
- bool check_C = false ;
- bool check_B = false ;
- bool check_D = false ;
-
- bool check_Theta = false ;
- bool check_Phi = false ;
- bool check_R = false ;
-// bool check_beta = false ;
-
- bool check_SI = false ;
- bool check_SILI = false ;
- bool check_CSI = false ;
- bool check_VIS = false ;
-
- bool ReadingStatus = false ;
-
-
- while (!ConfigFile.eof()) {
-
- getline(ConfigFile, LineBuffer);
-
- // If line is a Start Up MUST2 bloc, Reading toggle to true
- if (LineBuffer.compare(0, 11, "M2Telescope") == 0)
- {
- G4cout << "///" << G4endl ;
- G4cout << "Telescope found: " << G4endl ;
- ReadingStatus = true ;
-
- }
-
- // Else don't toggle to Reading Block Status
- else ReadingStatus = false ;
-
- // Reading Block
- while(ReadingStatus)
- {
-
- ConfigFile >> DataBuffer ;
- // Comment Line
- if(DataBuffer.compare(0, 1, "%") == 0) {
- ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );
-
- }
-
- // Finding another telescope (safety), toggle out
- else if (DataBuffer.compare(0, 11, "M2Telescope") == 0) {
- cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
- ReadingStatus = false ;
- }
-
- // Position method
-
- else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
- check_A = true;
- ConfigFile >> DataBuffer ;
- Ax = atof(DataBuffer.c_str()) ;
- Ax = Ax * mm ;
- ConfigFile >> DataBuffer ;
- Ay = atof(DataBuffer.c_str()) ;
- Ay = Ay * mm ;
- ConfigFile >> DataBuffer ;
- Az = atof(DataBuffer.c_str()) ;
- Az = Az * mm ;
-
- A = G4ThreeVector(Ax, Ay, Az);
- cout << "X1 Y1 corner position : " << A << endl;
-
- }
-
-
- else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
- check_B = true;
- ConfigFile >> DataBuffer ;
- Bx = atof(DataBuffer.c_str()) ;
- Bx = Bx * mm ;
- ConfigFile >> DataBuffer ;
- By = atof(DataBuffer.c_str()) ;
- By = By * mm ;
- ConfigFile >> DataBuffer ;
- Bz = atof(DataBuffer.c_str()) ;
- Bz = Bz * mm ;
-
- B = G4ThreeVector(Bx, By, Bz);
- cout << "X128 Y1 corner position : " << B << endl;
-
- }
-
-
- else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
- check_C = true;
- ConfigFile >> DataBuffer ;
- Cx = atof(DataBuffer.c_str()) ;
- Cx = Cx * mm ;
- ConfigFile >> DataBuffer ;
- Cy = atof(DataBuffer.c_str()) ;
- Cy = Cy * mm ;
- ConfigFile >> DataBuffer ;
- Cz = atof(DataBuffer.c_str()) ;
- Cz = Cz * mm ;
-
- C = G4ThreeVector(Cx, Cy, Cz);
- cout << "X1 Y128 corner position : " << C << endl;
-
- }
-
- else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
- check_D = true;
- ConfigFile >> DataBuffer ;
- Dx = atof(DataBuffer.c_str()) ;
- Dx = Dx * mm ;
- ConfigFile >> DataBuffer ;
- Dy = atof(DataBuffer.c_str()) ;
- Dy = Dy * mm ;
- ConfigFile >> DataBuffer ;
- Dz = atof(DataBuffer.c_str()) ;
- Dz = Dz * mm ;
-
- D = G4ThreeVector(Dx, Dy, Dz);
- cout << "X128 Y128 corner position : " << D << endl;
-
- }
-
- // End Position Method
-
- // Angle method
- else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
- check_Theta = true;
- ConfigFile >> DataBuffer ;
- Theta = atof(DataBuffer.c_str()) ;
- Theta = Theta * deg;
- cout << "Theta: " << Theta / deg << endl;
-
- }
-
- //Angle method
- else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
- check_Phi = true;
- ConfigFile >> DataBuffer ;
- Phi = atof(DataBuffer.c_str()) ;
- Phi = Phi * deg;
- cout << "Phi: " << Phi / deg << endl;
-
- }
-
- //Angle method
- else if (DataBuffer.compare(0, 2, "R=") == 0) {
- check_R = true;
- ConfigFile >> DataBuffer ;
- R = atof(DataBuffer.c_str()) ;
- R = R * mm;
- cout << "R: " << R / mm << endl;
-
- }
-
- //Angle method
- else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
-// check_beta = true;
- ConfigFile >> DataBuffer ;
- beta_u = atof(DataBuffer.c_str()) ;
- beta_u = beta_u * deg ;
- ConfigFile >> DataBuffer ;
- beta_v = atof(DataBuffer.c_str()) ;
- beta_v = beta_v * deg ;
- ConfigFile >> DataBuffer ;
- beta_w = atof(DataBuffer.c_str()) ;
- beta_w = beta_w * deg ;
- G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
-
- }
-
- // End Angle Method
-
-
- // Common part
- else if (DataBuffer.compare(0, 3, "SI=") == 0) {
- check_SI = true ;
- ConfigFile >> DataBuffer;
- SI = atof(DataBuffer.c_str()) ;
- G4cout << " SI= " << SI << G4endl ;
-
- }
-
-
- else if (DataBuffer.compare(0, 5, "SILI=") == 0) {
- check_SILI = true ;
- ConfigFile >> DataBuffer;
- SILI = atof(DataBuffer.c_str()) ;
- G4cout << " SILI= " << SILI << G4endl ;
-
- }
-
-
- else if (DataBuffer.compare(0, 4, "CSI=") == 0) {
- check_CSI = true ;
- ConfigFile >> DataBuffer;
- CSI = atof(DataBuffer.c_str()) ;
- G4cout << " CSI= " << CSI << G4endl ;
-
- }
-
-
- else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
- check_VIS = true ;
- ConfigFile >> DataBuffer;
- if (DataBuffer.compare(0, 3, "all") == 0)
- {
- m_non_sensitive_part_visiualisation = true;
- G4cout << " VIS= all" << G4endl ;
- }
-
- else
- G4cout << " VIS= Sensible Only" << G4endl ;
-
-
- }
-
- // If no MUST2 Token and no comment, toggle out
- else
- {ReadingStatus = false; G4cout << "other token " << DataBuffer << G4endl ;}
-
-
-
- /////////////////////////////////////////////////
- // If All necessary information there, toggle out
- if ( ( check_SI && check_SILI && check_CSI && check_VIS ) && ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R) ) )
- {
- ReadingStatus = false;
-
- ///Add The previously define telescope
- //With position method
- if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R))
- {
- AddTelescope( A,
- B,
- C,
- D,
- SI == 1,
- SILI == 1,
- CSI == 1);
- }
-
- //with angle method
- else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D))
- {
- AddTelescope( R,
- Theta,
- Phi,
- beta_u,
- beta_v,
- beta_w,
- SI == 1,
- SILI == 1,
- CSI == 1);
- }
-
-
- check_A = false ;
- check_C = false ;
- check_B = false ;
- check_D = false ;
-
- check_Theta = false ;
- check_Phi = false ;
- check_R = false ;
-// check_beta = false ;
-
- check_SI = false ;
- check_SILI = false ;
- check_CSI = false ;
- check_VIS = false ;
-
- }
- }
- }
+void MUST2Array::ReadConfiguration(string Path){
+ ifstream ConfigFile ;
+ ConfigFile.open(Path.c_str()) ;
+ string LineBuffer ;
+ string DataBuffer ;
+
+ // A:X1_Y1 --> X:1 Y:1
+ // B:X128_Y1 --> X:128 Y:1
+ // C:X1_Y128 --> X:1 Y:128
+ // D:X128_Y128 --> X:128 Y:128
+
+ G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz ;
+ G4ThreeVector A , B , C , D ;
+ G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0 ;
+ int SI = 0 , SILI = 0 , CSI = 0 ;
+
+ bool check_A = false ;
+ bool check_C = false ;
+ bool check_B = false ;
+ bool check_D = false ;
+
+ bool check_Theta = false ;
+ bool check_Phi = false ;
+ bool check_R = false ;
+ // bool check_beta = false ;
+
+ bool check_SI = false ;
+ bool check_SILI = false ;
+ bool check_CSI = false ;
+ bool check_VIS = false ;
+
+ bool ReadingStatus = false ;
+
+
+ while (!ConfigFile.eof()) {
+
+ getline(ConfigFile, LineBuffer);
+
+ // If line is a Start Up MUST2 bloc, Reading toggle to true
+ if (LineBuffer.compare(0, 11, "M2Telescope") == 0)
+ {
+ G4cout << "///" << G4endl ;
+ G4cout << "Telescope found: " << G4endl ;
+ ReadingStatus = true ;
+
+ }
+
+ // Else don't toggle to Reading Block Status
+ else ReadingStatus = false ;
+
+ // Reading Block
+ while(ReadingStatus)
+ {
+
+ ConfigFile >> DataBuffer ;
+ // Comment Line
+ if(DataBuffer.compare(0, 1, "%") == 0) {
+ ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );
+
+ }
+
+ // Finding another telescope (safety), toggle out
+ else if (DataBuffer.compare(0, 11, "M2Telescope") == 0) {
+ cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
+ ReadingStatus = false ;
+ }
+
+ // Position method
+
+ else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
+ check_A = true;
+ ConfigFile >> DataBuffer ;
+ Ax = atof(DataBuffer.c_str()) ;
+ Ax = Ax * mm ;
+ ConfigFile >> DataBuffer ;
+ Ay = atof(DataBuffer.c_str()) ;
+ Ay = Ay * mm ;
+ ConfigFile >> DataBuffer ;
+ Az = atof(DataBuffer.c_str()) ;
+ Az = Az * mm ;
+
+ A = G4ThreeVector(Ax, Ay, Az);
+ cout << "X1 Y1 corner position : " << A << endl;
+
+ }
+
+
+ else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
+ check_B = true;
+ ConfigFile >> DataBuffer ;
+ Bx = atof(DataBuffer.c_str()) ;
+ Bx = Bx * mm ;
+ ConfigFile >> DataBuffer ;
+ By = atof(DataBuffer.c_str()) ;
+ By = By * mm ;
+ ConfigFile >> DataBuffer ;
+ Bz = atof(DataBuffer.c_str()) ;
+ Bz = Bz * mm ;
+
+ B = G4ThreeVector(Bx, By, Bz);
+ cout << "X128 Y1 corner position : " << B << endl;
+
+ }
+
+
+ else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
+ check_C = true;
+ ConfigFile >> DataBuffer ;
+ Cx = atof(DataBuffer.c_str()) ;
+ Cx = Cx * mm ;
+ ConfigFile >> DataBuffer ;
+ Cy = atof(DataBuffer.c_str()) ;
+ Cy = Cy * mm ;
+ ConfigFile >> DataBuffer ;
+ Cz = atof(DataBuffer.c_str()) ;
+ Cz = Cz * mm ;
+
+ C = G4ThreeVector(Cx, Cy, Cz);
+ cout << "X1 Y128 corner position : " << C << endl;
+
+ }
+
+ else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
+ check_D = true;
+ ConfigFile >> DataBuffer ;
+ Dx = atof(DataBuffer.c_str()) ;
+ Dx = Dx * mm ;
+ ConfigFile >> DataBuffer ;
+ Dy = atof(DataBuffer.c_str()) ;
+ Dy = Dy * mm ;
+ ConfigFile >> DataBuffer ;
+ Dz = atof(DataBuffer.c_str()) ;
+ Dz = Dz * mm ;
+
+ D = G4ThreeVector(Dx, Dy, Dz);
+ cout << "X128 Y128 corner position : " << D << endl;
+
+ }
+
+ // End Position Method
+
+ // Angle method
+ else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
+ check_Theta = true;
+ ConfigFile >> DataBuffer ;
+ Theta = atof(DataBuffer.c_str()) ;
+ Theta = Theta * deg;
+ cout << "Theta: " << Theta / deg << endl;
+
+ }
+
+ //Angle method
+ else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
+ check_Phi = true;
+ ConfigFile >> DataBuffer ;
+ Phi = atof(DataBuffer.c_str()) ;
+ Phi = Phi * deg;
+ cout << "Phi: " << Phi / deg << endl;
+
+ }
+
+ //Angle method
+ else if (DataBuffer.compare(0, 2, "R=") == 0) {
+ check_R = true;
+ ConfigFile >> DataBuffer ;
+ R = atof(DataBuffer.c_str()) ;
+ R = R * mm;
+ cout << "R: " << R / mm << endl;
+
+ }
+
+ //Angle method
+ else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
+ // check_beta = true;
+ ConfigFile >> DataBuffer ;
+ beta_u = atof(DataBuffer.c_str()) ;
+ beta_u = beta_u * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_v = atof(DataBuffer.c_str()) ;
+ beta_v = beta_v * deg ;
+ ConfigFile >> DataBuffer ;
+ beta_w = atof(DataBuffer.c_str()) ;
+ beta_w = beta_w * deg ;
+ G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ;
+
+ }
+
+ // End Angle Method
+
+
+ // Common part
+ else if (DataBuffer.compare(0, 3, "SI=") == 0) {
+ check_SI = true ;
+ ConfigFile >> DataBuffer;
+ SI = atof(DataBuffer.c_str()) ;
+ G4cout << " SI= " << SI << G4endl ;
+
+ }
+
+
+ else if (DataBuffer.compare(0, 5, "SILI=") == 0) {
+ check_SILI = true ;
+ ConfigFile >> DataBuffer;
+ SILI = atof(DataBuffer.c_str()) ;
+ G4cout << " SILI= " << SILI << G4endl ;
+
+ }
+
+
+ else if (DataBuffer.compare(0, 4, "CSI=") == 0) {
+ check_CSI = true ;
+ ConfigFile >> DataBuffer;
+ CSI = atof(DataBuffer.c_str()) ;
+ G4cout << " CSI= " << CSI << G4endl ;
+
+ }
+
+
+ else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
+ check_VIS = true ;
+ ConfigFile >> DataBuffer;
+ if (DataBuffer.compare(0, 3, "all") == 0)
+ {
+ m_non_sensitive_part_visiualisation = true;
+ G4cout << " VIS= all" << G4endl ;
+ }
+
+ else
+ G4cout << " VIS= Sensible Only" << G4endl ;
+
+
+ }
+
+ // If no MUST2 Token and no comment, toggle out
+ else
+ {ReadingStatus = false; G4cout << "other token " << DataBuffer << G4endl ;}
+
+
+
+ /////////////////////////////////////////////////
+ // If All necessary information there, toggle out
+ if ( ( check_SI && check_SILI && check_CSI && check_VIS ) && ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R) ) )
+ {
+ ReadingStatus = false;
+
+ ///Add The previously define telescope
+ //With position method
+ if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R))
+ {
+ AddTelescope( A,
+ B,
+ C,
+ D,
+ SI == 1,
+ SILI == 1,
+ CSI == 1);
+ }
+
+ //with angle method
+ else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D))
+ {
+ AddTelescope( R,
+ Theta,
+ Phi,
+ beta_u,
+ beta_v,
+ beta_w,
+ SI == 1,
+ SILI == 1,
+ CSI == 1);
+ }
+
+
+ check_A = false ;
+ check_C = false ;
+ check_B = false ;
+ check_D = false ;
+
+ check_Theta = false ;
+ check_Phi = false ;
+ check_R = false ;
+ // check_beta = false ;
+
+ check_SI = false ;
+ check_SILI = false ;
+ check_CSI = false ;
+ check_VIS = false ;
+
+ }
+ }
+ }
}
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
-void MUST2Array::ConstructDetector(G4LogicalVolume* world)
-{
- G4RotationMatrix* MMrot = NULL ;
- G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;
- G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ;
- bool Si = true ;
- bool SiLi = true ;
- bool CsI = true ;
-
- G4int NumberOfTelescope = m_DefinitionType.size() ;
-
- for (G4int i = 0 ; i < NumberOfTelescope ; i++) {
- // By Point
- if (m_DefinitionType[i]) {
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // w perpendicular to (u,v) plan and pointing CsI
- MMu = m_X128_Y1[i] - m_X1_Y1[i] ;
- MMu = MMu.unit() ;
-
- MMv = m_X1_Y128[i] - m_X1_Y1[i] ;
- MMv = MMv.unit() ;
-
- MMw = MMv.cross(MMu) ;
- // if (MMw.z() > 0)MMw = MMv.cross(MMu) ;
- MMw = MMw.unit() ;
-
- MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4 ;
-
- // Passage Matrix from Lab Referential to Telescope Referential
- MMrot = new G4RotationMatrix(MMv, MMu, MMw);
- MMpos = MMw * Length * 0.5 + MMCenter ;
- }
+void MUST2Array::ConstructDetector(G4LogicalVolume* world){
+ G4RotationMatrix* MMrot = NULL ;
+ G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;
+ G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ;
+ bool Si = true ;
+ bool SiLi = true ;
+ bool CsI = true ;
+
+ G4int NumberOfTelescope = m_DefinitionType.size() ;
+
+ for (G4int i = 0 ; i < NumberOfTelescope ; i++) {
+ // By Point
+ if (m_DefinitionType[i]) {
+ // (u,v,w) unitary vector associated to telescope referencial
+ // (u,v) // to silicon plan
+ // w perpendicular to (u,v) plan and pointing CsI
+ MMu = m_X128_Y1[i] - m_X1_Y1[i] ;
+ MMu = MMu.unit() ;
+
+ MMv = m_X1_Y128[i] - m_X1_Y1[i] ;
+ MMv = MMv.unit() ;
+
+ MMw = MMv.cross(MMu) ;
+ // if (MMw.z() > 0)MMw = MMv.cross(MMu) ;
+ MMw = MMw.unit() ;
+
+ MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4 ;
+
+ // Passage Matrix from Lab Referential to Telescope Referential
+ MMrot = new G4RotationMatrix(MMv, MMu, MMw);
+ MMpos = MMw * Length * 0.5 + MMCenter ;
+ }
- // By Angle
- else {
- G4double Theta = m_Theta[i] ;
- G4double Phi = m_Phi[i] ;
-
- // (u,v,w) unitary vector associated to telescope referencial
- // (u,v) // to silicon plan
- // 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);
- MMw = G4ThreeVector(wX, wY, wZ) ;
-
- // vector corresponding to the center of the module
- G4ThreeVector CT = MMw;
-
- // vector parallel to one axis of silicon plane
- G4double ii = cos(Theta / rad) * cos(Phi / rad);
- G4double jj = cos(Theta / rad) * sin(Phi / rad);
- G4double kk = -sin(Theta / rad);
- G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
-
- MMw = MMw.unit();
- MMu = MMw.cross(Y);
- MMv = MMw.cross(MMu);
- MMv = MMv.unit();
- MMu = MMu.unit();
-
- // Passage Matrix from Lab Referential to Telescope Referential
- // MUST2
- MMrot = new G4RotationMatrix(MMu, MMv, MMw);
- // Telescope is rotate of Beta angle around MMv axis.
- MMrot->rotate(m_beta_u[i], MMu);
- MMrot->rotate(m_beta_v[i], MMv);
- MMrot->rotate(m_beta_w[i], MMw);
- // translation to place Telescope
- MMpos = MMw * Length * 0.5 + CT ;
- }
+ // By Angle
+ else {
+ G4double Theta = m_Theta[i] ;
+ G4double Phi = m_Phi[i] ;
+
+ // (u,v,w) unitary vector associated to telescope referencial
+ // (u,v) // to silicon plan
+ // 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);
+ MMw = G4ThreeVector(wX, wY, wZ) ;
+
+ // vector corresponding to the center of the module
+ G4ThreeVector CT = MMw;
+
+ // vector parallel to one axis of silicon plane
+ G4double ii = cos(Theta / rad) * cos(Phi / rad);
+ G4double jj = cos(Theta / rad) * sin(Phi / rad);
+ G4double kk = -sin(Theta / rad);
+ G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
+
+ MMw = MMw.unit();
+ MMu = MMw.cross(Y);
+ MMv = MMw.cross(MMu);
+ MMv = MMv.unit();
+ MMu = MMu.unit();
+
+ // Passage Matrix from Lab Referential to Telescope Referential
+ // MUST2
+ MMrot = new G4RotationMatrix(MMu, MMv, MMw);
+ // Telescope is rotate of Beta angle around MMv axis.
+ MMrot->rotate(m_beta_u[i], MMu);
+ MMrot->rotate(m_beta_v[i], MMv);
+ MMrot->rotate(m_beta_w[i], MMw);
+ // translation to place Telescope
+ MMpos = MMw * Length * 0.5 + CT ;
+ }
- Si = m_wSi[i] ;
- SiLi = m_wSiLi[i] ;
- CsI = m_wCsI[i] ;
+ Si = m_wSi[i] ;
+ SiLi = m_wSiLi[i] ;
+ CsI = m_wCsI[i] ;
- VolumeMaker(i + 1 , MMpos , MMrot , Si , SiLi , CsI , world);
- }
+ VolumeMaker(i + 1 , MMpos , MMrot , Si , SiLi , CsI , world);
+ }
- delete MMrot ;
+ delete MMrot ;
}
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
-void MUST2Array::InitializeRootOutput()
-{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("MUST2", "TMust2Data", &m_Event) ;
+void MUST2Array::InitializeRootOutput(){
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ pTree->Branch("MUST2", "TMust2Data", &m_Event) ;
}
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
-void MUST2Array::ReadSensitive(const G4Event* event)
-{
- G4String DetectorNumber ;
- m_Event->Clear() ;
-
-//////////////////////////////////////////////////////////////////////////////////////
-//////////////////////// Used to Read Event Map of detector //////////////////////////
-//////////////////////////////////////////////////////////////////////////////////////
-
-// Si
- std::map<G4int, G4int*>::iterator DetectorNumber_itr;
- std::map<G4int, G4double*>::iterator Energy_itr;
- std::map<G4int, G4double*>::iterator Time_itr;
- std::map<G4int, G4int*>::iterator X_itr;
- std::map<G4int, G4int*>::iterator Y_itr;
- std::map<G4int, G4double*>::iterator Pos_X_itr;
- std::map<G4int, G4double*>::iterator Pos_Y_itr;
- std::map<G4int, G4double*>::iterator Pos_Z_itr;
- std::map<G4int, G4double*>::iterator Ang_Theta_itr;
- std::map<G4int, G4double*>::iterator Ang_Phi_itr;
-
- G4THitsMap<G4int>* DetectorNumberHitMap;
- G4THitsMap<G4double>* EnergyHitMap;
- G4THitsMap<G4double>* TimeHitMap;
- G4THitsMap<G4int>* XHitMap;
- G4THitsMap<G4int>* YHitMap;
- G4THitsMap<G4double>* PosXHitMap;
- G4THitsMap<G4double>* PosYHitMap;
- G4THitsMap<G4double>* PosZHitMap;
- G4THitsMap<G4double>* AngThetaHitMap;
- G4THitsMap<G4double>* AngPhiHitMap;
-
-// Si(Li)
- std::map<G4int, G4double*>::iterator SiLiEnergy_itr;
- std::map<G4int, G4int*>::iterator SiLiPadNbr_itr;
- G4THitsMap<G4double>* SiLiEnergyHitMap;
- G4THitsMap<G4int>* SiLiPadNbrHitMap;
-
-// CsI
- std::map<G4int, G4double*>::iterator CsIEnergy_itr;
- std::map<G4int, G4int*>::iterator CsICristalNbr_itr;
- G4THitsMap<G4double>* CsIEnergyHitMap;
- G4THitsMap<G4int>* CsICristalNbrHitMap;
-//////////////////////////////////////////////////////////////////////////////////////
-//////////////////////////////////////////////////////////////////////////////////////
-
- // Read the Scorer associate to the Silicon Strip
-
- //Detector Number
- G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/DetectorNumber");
- DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID));
- DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin();
-
- //Energy
- G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripEnergy") ;
- EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
- Energy_itr = EnergyHitMap->GetMap()->begin();
-
- //Time of Flight
- G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripTime");
- TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID));
- Time_itr = TimeHitMap->GetMap()->begin();
-
- //Strip Number X
- G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberX") ;
- XHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID));
- X_itr = XHitMap->GetMap()->begin();
-
- //Strip Number Y
- G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberY");
- YHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID));
- Y_itr = YHitMap->GetMap()->begin();
-
- //Interaction Coordinate X
- G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordX");
- PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID));
- Pos_X_itr = PosXHitMap->GetMap()->begin();
-
- //Interaction Coordinate Y
- G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordY");
- PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID));
- Pos_Y_itr = PosYHitMap->GetMap()->begin();
-
- //Interaction Coordinate Z
- G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordZ");
- PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID));
- Pos_Z_itr = PosXHitMap->GetMap()->begin();
-
- //Interaction Coordinate Angle Theta
- G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngTheta");
- AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID));
- Ang_Theta_itr = AngThetaHitMap->GetMap()->begin();
-
- //Interaction Coordinate Angle Phi
- G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngPhi");
- AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID));
- Ang_Phi_itr = AngPhiHitMap->GetMap()->begin();
-
- // Read the Scorer associate to the SiLi
- //Energy
- G4int SiLiEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiEnergy");
- SiLiEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiEnergyCollectionID));
- SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin();
-
- G4int SiLiPadNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiPadNbr");
- SiLiPadNbrHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(SiLiPadNbrCollectionID));
- SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin();
-
- // Read the Scorer associate to the CsI crystal
- //Energy
- G4int CsIEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIEnergy");
- CsIEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIEnergyCollectionID));
- CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin();
-
- G4int CsICristalNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsICristalNbr");
- CsICristalNbrHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICristalNbrCollectionID));
- CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin();
-
-
-/////////////////////
-
- G4int sizeN = DetectorNumberHitMap->entries();
- G4int sizeE = EnergyHitMap->entries();
- G4int sizeT = TimeHitMap->entries();
- G4int sizeX = XHitMap->entries();
- G4int sizeY = YHitMap->entries();
-
- // Loop on Telescope Number entry
- for (G4int l = 0 ; l < sizeN ; l++) {
- G4int N = *(DetectorNumber_itr->second);
- G4int NTrackID = DetectorNumber_itr->first - N;
-
-
- if (N > 0) {
-
- m_Event->SetMMStripXEDetectorNbr(N) ;
- m_Event->SetMMStripYEDetectorNbr(N) ;
- m_Event->SetMMStripXTDetectorNbr(N) ;
- m_Event->SetMMStripYTDetectorNbr(N) ;
-
- // Energy
- Energy_itr = EnergyHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < sizeE ; h++) {
- G4int ETrackID = Energy_itr->first - N;
- G4double E = *(Energy_itr->second);
-
- if (ETrackID == NTrackID) {
- m_Event->SetMMStripXEEnergy(RandGauss::shoot(E, ResoStrip));
- m_Event->SetMMStripYEEnergy(RandGauss::shoot(E, ResoStrip));
- }
-
- Energy_itr++;
- }
-
-
- // Time
- Time_itr = TimeHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < sizeT ; h++) {
- G4int TTrackID = Time_itr->first - N;
- G4double T = *(Time_itr->second);
-
- if (TTrackID == NTrackID) {
- m_Event->SetMMStripXTTime(RandGauss::shoot(T, ResoTimeMust)) ;
- m_Event->SetMMStripYTTime(RandGauss::shoot(T, ResoTimeMust)) ;
- }
-
- Time_itr++;
- }
-
-
- // X
- X_itr = XHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < sizeX ; h++) {
- G4int XTrackID = X_itr->first - N;
- G4int X = *(X_itr->second);
- if (XTrackID == NTrackID) {
- m_Event->SetMMStripXEStripNbr(X);
- m_Event->SetMMStripXTStripNbr(X);
- }
-
- X_itr++;
- }
-
- // Y
- Y_itr = YHitMap->GetMap()->begin() ;
- for (G4int h = 0 ; h < sizeY ; h++) {
- G4int YTrackID = Y_itr->first - N ;
- G4int Y = *(Y_itr->second);
- if (YTrackID == NTrackID) {
- m_Event->SetMMStripYEStripNbr(Y);
- m_Event->SetMMStripYTStripNbr(Y);
- }
-
- Y_itr++;
- }
-
- // Pos X
- Pos_X_itr = PosXHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < PosXHitMap->entries() ; h++) {
- G4int PosXTrackID = Pos_X_itr->first - N ;
- G4double PosX = *(Pos_X_itr->second) ;
- if (PosXTrackID == NTrackID) {
- ms_InterCoord->SetDetectedPositionX(PosX) ;
- }
- Pos_X_itr++;
- }
-
- // Pos Y
- Pos_Y_itr = PosYHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < PosYHitMap->entries() ; h++) {
- G4int PosYTrackID = Pos_Y_itr->first - N ;
- G4double PosY = *(Pos_Y_itr->second) ;
- if (PosYTrackID == NTrackID) {
- ms_InterCoord->SetDetectedPositionY(PosY) ;
- }
- Pos_Y_itr++;
- }
-
- // Pos Z
- Pos_Z_itr = PosZHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < PosZHitMap->entries() ; h++) {
- G4int PosZTrackID = Pos_Z_itr->first - N ;
- G4double PosZ = *(Pos_Z_itr->second) ;
- if (PosZTrackID == NTrackID) {
- ms_InterCoord->SetDetectedPositionZ(PosZ) ;
- }
- Pos_Z_itr++;
- }
-
- // Angle Theta
- Ang_Theta_itr = AngThetaHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < AngThetaHitMap->entries(); h++) {
- G4int AngThetaTrackID = Ang_Theta_itr->first - N ;
- G4double AngTheta = *(Ang_Theta_itr->second) ;
- if (AngThetaTrackID == NTrackID) {
- ms_InterCoord->SetDetectedAngleTheta(AngTheta) ;
- }
- Ang_Theta_itr++;
- }
-
- // Angle Phi
- Ang_Phi_itr = AngPhiHitMap->GetMap()->begin();
- for (G4int h = 0 ; h < AngPhiHitMap->entries() ; h++) {
- G4int AngPhiTrackID = Ang_Phi_itr->first - N ;
- G4double AngPhi = *(Ang_Phi_itr->second) ;
- if (AngPhiTrackID == NTrackID) {
- ms_InterCoord->SetDetectedAnglePhi(AngPhi) ;
- }
- Ang_Phi_itr++;
- }
-
-
- // Si(Li)
- SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ;
- SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin() ;
-
- for (G4int h = 0 ; h < SiLiEnergyHitMap->entries() ; h++) {
- G4int SiLiEnergyTrackID = SiLiEnergy_itr->first -N ;
- G4double SiLiEnergy = *(SiLiEnergy_itr->second) ;
-
- G4int PadNbr = *(SiLiPadNbr_itr->second) ;
-
- if (SiLiEnergyTrackID == NTrackID) {
- m_Event->SetMMSiLiEEnergy(RandGauss::shoot(SiLiEnergy, ResoSiLi)) ;
- m_Event->SetMMSiLiEPadNbr(PadNbr);
- m_Event->SetMMSiLiTPadNbr(PadNbr);
- m_Event->SetMMSiLiTTime(1);
- m_Event->SetMMSiLiTDetectorNbr(N);
- m_Event->SetMMSiLiEDetectorNbr(N);
- }
-
- SiLiEnergy_itr++ ;
- SiLiPadNbr_itr++ ;
- }
-
-
- // CsI
- CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ;
- CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin() ;
-
- for (G4int h = 0 ; h < CsIEnergyHitMap->entries() ; h++) {
- G4int CsIEnergyTrackID = CsIEnergy_itr->first-N ;
- G4double CsIEnergy = *(CsIEnergy_itr->second) ;
-
- G4int CristalNumber = *(CsICristalNbr_itr->second) ;
- if (CsIEnergyTrackID == NTrackID) {
- m_Event->SetMMCsIEEnergy(RandGauss::shoot(CsIEnergy, ResoCsI*sqrt(CsIEnergy)));
- m_Event->SetMMCsIECristalNbr(CristalNumber);
- m_Event->SetMMCsITCristalNbr(CristalNumber);
- m_Event->SetMMCsITTime(1);
- m_Event->SetMMCsITDetectorNbr(N);
- m_Event->SetMMCsIEDetectorNbr(N);
- }
-
- CsIEnergy_itr++ ;
- CsICristalNbr_itr++ ;
- }
+void MUST2Array::ReadSensitive(const G4Event* event){
+ G4String DetectorNumber ;
+ m_Event->Clear() ;
+
+ //////////////////////////////////////////////////////////////////////////////////////
+ //////////////////////// Used to Read Event Map of detector //////////////////////////
+ //////////////////////////////////////////////////////////////////////////////////////
+
+ // Si
+ std::map<G4int, G4int*>::iterator DetectorNumber_itr;
+ std::map<G4int, G4double*>::iterator Energy_itr;
+ std::map<G4int, G4double*>::iterator Time_itr;
+ std::map<G4int, G4int*>::iterator X_itr;
+ std::map<G4int, G4int*>::iterator Y_itr;
+ std::map<G4int, G4double*>::iterator Pos_X_itr;
+ std::map<G4int, G4double*>::iterator Pos_Y_itr;
+ std::map<G4int, G4double*>::iterator Pos_Z_itr;
+ std::map<G4int, G4double*>::iterator Ang_Theta_itr;
+ std::map<G4int, G4double*>::iterator Ang_Phi_itr;
+
+ G4THitsMap<G4int>* DetectorNumberHitMap;
+ G4THitsMap<G4double>* EnergyHitMap;
+ G4THitsMap<G4double>* TimeHitMap;
+ G4THitsMap<G4int>* XHitMap;
+ G4THitsMap<G4int>* YHitMap;
+ G4THitsMap<G4double>* PosXHitMap;
+ G4THitsMap<G4double>* PosYHitMap;
+ G4THitsMap<G4double>* PosZHitMap;
+ G4THitsMap<G4double>* AngThetaHitMap;
+ G4THitsMap<G4double>* AngPhiHitMap;
+
+ // Si(Li)
+ std::map<G4int, G4double*>::iterator SiLiEnergy_itr;
+ std::map<G4int, G4int*>::iterator SiLiPadNbr_itr;
+ G4THitsMap<G4double>* SiLiEnergyHitMap;
+ G4THitsMap<G4int>* SiLiPadNbrHitMap;
+
+ // CsI
+ std::map<G4int, G4double*>::iterator CsIEnergy_itr;
+ std::map<G4int, G4int*>::iterator CsICristalNbr_itr;
+ G4THitsMap<G4double>* CsIEnergyHitMap;
+ G4THitsMap<G4int>* CsICristalNbrHitMap;
+ //////////////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////////////////
+
+ // Read the Scorer associate to the Silicon Strip
+
+ //Detector Number
+ G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/DetectorNumber");
+ DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID));
+ DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin();
+
+ //Energy
+ G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripEnergy") ;
+ EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
+ Energy_itr = EnergyHitMap->GetMap()->begin();
+
+ //Time of Flight
+ G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripTime");
+ TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID));
+ Time_itr = TimeHitMap->GetMap()->begin();
+
+ //Strip Number X
+ G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberX") ;
+ XHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID));
+ X_itr = XHitMap->GetMap()->begin();
+
+ //Strip Number Y
+ G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberY");
+ YHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID));
+ Y_itr = YHitMap->GetMap()->begin();
+
+ //Interaction Coordinate X
+ G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordX");
+ PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID));
+ Pos_X_itr = PosXHitMap->GetMap()->begin();
+
+ //Interaction Coordinate Y
+ G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordY");
+ PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID));
+ Pos_Y_itr = PosYHitMap->GetMap()->begin();
+
+ //Interaction Coordinate Z
+ G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordZ");
+ PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID));
+ Pos_Z_itr = PosXHitMap->GetMap()->begin();
+
+ //Interaction Coordinate Angle Theta
+ G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngTheta");
+ AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID));
+ Ang_Theta_itr = AngThetaHitMap->GetMap()->begin();
+
+ //Interaction Coordinate Angle Phi
+ G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngPhi");
+ AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID));
+ Ang_Phi_itr = AngPhiHitMap->GetMap()->begin();
+
+ // Read the Scorer associate to the SiLi
+ //Energy
+ G4int SiLiEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiEnergy");
+ SiLiEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiEnergyCollectionID));
+ SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin();
+
+ G4int SiLiPadNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiPadNbr");
+ SiLiPadNbrHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(SiLiPadNbrCollectionID));
+ SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin();
+
+ // Read the Scorer associate to the CsI crystal
+ //Energy
+ G4int CsIEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIEnergy");
+ CsIEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIEnergyCollectionID));
+ CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin();
+
+ G4int CsICristalNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsICristalNbr");
+ CsICristalNbrHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICristalNbrCollectionID));
+ CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin();
+
+
+ /////////////////////
+
+ G4int sizeN = DetectorNumberHitMap->entries();
+ G4int sizeE = EnergyHitMap->entries();
+ G4int sizeT = TimeHitMap->entries();
+ G4int sizeX = XHitMap->entries();
+ G4int sizeY = YHitMap->entries();
+
+ // Loop on Telescope Number entry
+ for (G4int l = 0 ; l < sizeN ; l++) {
+ G4int N = *(DetectorNumber_itr->second);
+ G4int NTrackID = DetectorNumber_itr->first - N;
+
+
+ if (N > 0) {
+
+ m_Event->SetMMStripXEDetectorNbr(N) ;
+ m_Event->SetMMStripYEDetectorNbr(N) ;
+ m_Event->SetMMStripXTDetectorNbr(N) ;
+ m_Event->SetMMStripYTDetectorNbr(N) ;
+
+ // Energy
+ Energy_itr = EnergyHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < sizeE ; h++) {
+ G4int ETrackID = Energy_itr->first - N;
+ G4double E = *(Energy_itr->second);
+
+ if (ETrackID == NTrackID) {
+ m_Event->SetMMStripXEEnergy(RandGauss::shoot(E, ResoStrip));
+ m_Event->SetMMStripYEEnergy(RandGauss::shoot(E, ResoStrip));
+ }
+
+ Energy_itr++;
+ }
+ // Time
+ Time_itr = TimeHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < sizeT ; h++) {
+ G4int TTrackID = Time_itr->first - N;
+ G4double T = *(Time_itr->second);
+ if (TTrackID == NTrackID) {
+ m_Event->SetMMStripXTTime(RandGauss::shoot(T, ResoTimeMust)) ;
+ m_Event->SetMMStripYTTime(RandGauss::shoot(T, ResoTimeMust)) ;
}
- DetectorNumber_itr++;
+ Time_itr++;
+ }
+
+
+ // X
+ X_itr = XHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < sizeX ; h++) {
+ G4int XTrackID = X_itr->first - N;
+ G4int X = *(X_itr->second);
+ if (XTrackID == NTrackID) {
+ m_Event->SetMMStripXEStripNbr(X);
+ m_Event->SetMMStripXTStripNbr(X);
+ }
+
+ X_itr++;
+ }
+
+ // Y
+ Y_itr = YHitMap->GetMap()->begin() ;
+ for (G4int h = 0 ; h < sizeY ; h++) {
+ G4int YTrackID = Y_itr->first - N ;
+ G4int Y = *(Y_itr->second);
+ if (YTrackID == NTrackID) {
+ m_Event->SetMMStripYEStripNbr(Y);
+ m_Event->SetMMStripYTStripNbr(Y);
+ }
+
+ Y_itr++;
+ }
+
+ // Pos X
+ Pos_X_itr = PosXHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < PosXHitMap->entries() ; h++) {
+ G4int PosXTrackID = Pos_X_itr->first - N ;
+ G4double PosX = *(Pos_X_itr->second) ;
+ if (PosXTrackID == NTrackID) {
+ ms_InterCoord->SetDetectedPositionX(PosX) ;
+ }
+ Pos_X_itr++;
+ }
+
+ // Pos Y
+ Pos_Y_itr = PosYHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < PosYHitMap->entries() ; h++) {
+ G4int PosYTrackID = Pos_Y_itr->first - N ;
+ G4double PosY = *(Pos_Y_itr->second) ;
+ if (PosYTrackID == NTrackID) {
+ ms_InterCoord->SetDetectedPositionY(PosY) ;
+ }
+ Pos_Y_itr++;
+ }
+
+ // Pos Z
+ Pos_Z_itr = PosZHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < PosZHitMap->entries() ; h++) {
+ G4int PosZTrackID = Pos_Z_itr->first - N ;
+ G4double PosZ = *(Pos_Z_itr->second) ;
+ if (PosZTrackID == NTrackID) {
+ ms_InterCoord->SetDetectedPositionZ(PosZ) ;
+ }
+ Pos_Z_itr++;
+ }
+
+ // Angle Theta
+ Ang_Theta_itr = AngThetaHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < AngThetaHitMap->entries(); h++) {
+ G4int AngThetaTrackID = Ang_Theta_itr->first - N ;
+ G4double AngTheta = *(Ang_Theta_itr->second) ;
+ if (AngThetaTrackID == NTrackID) {
+ ms_InterCoord->SetDetectedAngleTheta(AngTheta) ;
+ }
+ Ang_Theta_itr++;
+ }
+
+ // Angle Phi
+ Ang_Phi_itr = AngPhiHitMap->GetMap()->begin();
+ for (G4int h = 0 ; h < AngPhiHitMap->entries() ; h++) {
+ G4int AngPhiTrackID = Ang_Phi_itr->first - N ;
+ G4double AngPhi = *(Ang_Phi_itr->second) ;
+ if (AngPhiTrackID == NTrackID) {
+ ms_InterCoord->SetDetectedAnglePhi(AngPhi) ;
+ }
+ Ang_Phi_itr++;
+ }
+
+ // Si(Li)
+ SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ;
+ SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin() ;
+
+ for (G4int h = 0 ; h < SiLiEnergyHitMap->entries() ; h++) {
+ G4int SiLiEnergyTrackID = SiLiEnergy_itr->first -N ;
+ G4double SiLiEnergy = *(SiLiEnergy_itr->second) ;
+
+ G4int PadNbr = *(SiLiPadNbr_itr->second) ;
+
+ if (SiLiEnergyTrackID == NTrackID) {
+ m_Event->SetMMSiLiEEnergy(RandGauss::shoot(SiLiEnergy, ResoSiLi)) ;
+ m_Event->SetMMSiLiEPadNbr(PadNbr);
+ m_Event->SetMMSiLiTPadNbr(PadNbr);
+ m_Event->SetMMSiLiTTime(1);
+ m_Event->SetMMSiLiTDetectorNbr(N);
+ m_Event->SetMMSiLiEDetectorNbr(N);
+ }
+
+ SiLiEnergy_itr++ ;
+ SiLiPadNbr_itr++ ;
+ }
+
+ // CsI
+ CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ;
+ CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin() ;
+
+ for (G4int h = 0 ; h < CsIEnergyHitMap->entries() ; h++) {
+ G4int CsIEnergyTrackID = CsIEnergy_itr->first-N ;
+ G4double CsIEnergy = *(CsIEnergy_itr->second) ;
+
+ G4int CristalNumber = *(CsICristalNbr_itr->second) ;
+ if (CsIEnergyTrackID == NTrackID) {
+ m_Event->SetMMCsIEEnergy(RandGauss::shoot(CsIEnergy, ResoCsI*sqrt(CsIEnergy)));
+ m_Event->SetMMCsIECristalNbr(CristalNumber);
+ m_Event->SetMMCsITCristalNbr(CristalNumber);
+ m_Event->SetMMCsITTime(1);
+ m_Event->SetMMCsITDetectorNbr(N);
+ m_Event->SetMMCsIEDetectorNbr(N);
+ }
+
+ CsIEnergy_itr++ ;
+ CsICristalNbr_itr++ ;
+ }
+
}
-
- // clear map for next event
- DetectorNumberHitMap ->clear();
- EnergyHitMap ->clear() ;
- TimeHitMap ->clear() ;
- XHitMap ->clear() ;
- YHitMap ->clear() ;
- SiLiEnergyHitMap ->clear() ;
- SiLiPadNbrHitMap ->clear() ;
- CsIEnergyHitMap ->clear() ;
- CsICristalNbrHitMap ->clear() ;
- PosXHitMap ->clear() ;
- PosYHitMap ->clear() ;
- PosZHitMap ->clear() ;
- AngThetaHitMap ->clear() ;
- AngPhiHitMap ->clear() ;
-
+
+ DetectorNumber_itr++;
+ }
+
+ // clear map for next event
+ DetectorNumberHitMap ->clear();
+ EnergyHitMap ->clear() ;
+ TimeHitMap ->clear() ;
+ XHitMap ->clear() ;
+ YHitMap ->clear() ;
+ SiLiEnergyHitMap ->clear() ;
+ SiLiPadNbrHitMap ->clear() ;
+ CsIEnergyHitMap ->clear() ;
+ CsICristalNbrHitMap ->clear() ;
+ PosXHitMap ->clear() ;
+ PosYHitMap ->clear() ;
+ PosZHitMap ->clear() ;
+ AngThetaHitMap ->clear() ;
+ AngPhiHitMap ->clear() ;
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void MUST2Array::InitializeScorers()
- {
- // Silicon Associate Scorer
- m_StripScorer = new G4MultiFunctionalDetector("MUST2_StripScorer");
-
- G4VPrimitiveScorer* DetNbr = new GENERALSCORERS::PSDetectorNumber("DetectorNumber","MUST2Telescope", 0);
- G4VPrimitiveScorer* Energy = new GENERALSCORERS::PSEnergy("StripEnergy","MUST2Telescope", 0);
- G4VPrimitiveScorer* TOF = new GENERALSCORERS::PSTOF("StripTime","MUST2Telescope", 0);
-
- G4VPrimitiveScorer* StripPositionX = new PSStripNumberX("StripNumberX", 0, SiliconFace, 128);
- G4VPrimitiveScorer* StripPositionY = new PSStripNumberY("StripNumberY", 0, SiliconFace, 128);
-
- G4VPrimitiveScorer* InteractionCoordinatesX = new GENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","MUST2Telescope", 0);
- G4VPrimitiveScorer* InteractionCoordinatesY = new GENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","MUST2Telescope", 0);
- G4VPrimitiveScorer* InteractionCoordinatesZ = new GENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","MUST2Telescope", 0);
- G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new GENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","MUST2Telescope", 0);
- G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new GENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","MUST2Telescope", 0) ;
-
-
- //and register it to the multifunctionnal detector
- m_StripScorer->RegisterPrimitive(DetNbr);
- m_StripScorer->RegisterPrimitive(Energy);
- m_StripScorer->RegisterPrimitive(TOF);
- m_StripScorer->RegisterPrimitive(StripPositionX);
- m_StripScorer->RegisterPrimitive(StripPositionY);
- m_StripScorer->RegisterPrimitive(InteractionCoordinatesX);
- m_StripScorer->RegisterPrimitive(InteractionCoordinatesY);
- m_StripScorer->RegisterPrimitive(InteractionCoordinatesZ);
- m_StripScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
- m_StripScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);
-
- // SiLi Associate Scorer
- m_SiLiScorer = new G4MultiFunctionalDetector("MUST2_SiLiScorer");
- G4VPrimitiveScorer* SiLiEnergy = new GENERALSCORERS::PSEnergy("SiLiEnergy","MUST2Telescope", 0) ;
- G4VPrimitiveScorer* SiLiPadNbr = new PSPadOrCristalNumber("SiLiPadNbr",0) ;
- m_SiLiScorer->RegisterPrimitive(SiLiEnergy);
- m_SiLiScorer->RegisterPrimitive(SiLiPadNbr);
-
- // CsI Associate Scorer
- m_CsIScorer = new G4MultiFunctionalDetector("MUST2_CsIScorer");
- G4VPrimitiveScorer* CsIEnergy = new GENERALSCORERS::PSEnergy("CsIEnergy","MUST2Telescope", 0) ;
- G4VPrimitiveScorer* CsICristalNbr = new PSPadOrCristalNumber("CsICristalNbr",0) ;
- m_CsIScorer->RegisterPrimitive(CsIEnergy) ;
- m_CsIScorer->RegisterPrimitive(CsICristalNbr) ;
-
- // Add All Scorer to the Global Scorer Manager
- G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer) ;
- G4SDManager::GetSDMpointer()->AddNewDetector(m_SiLiScorer) ;
- G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ;
- }
+void MUST2Array::InitializeScorers() {
+ // Silicon Associate Scorer
+ m_StripScorer = new G4MultiFunctionalDetector("MUST2_StripScorer");
+
+ G4VPrimitiveScorer* DetNbr = new GENERALSCORERS::PSDetectorNumber("DetectorNumber","MUST2Telescope", 0);
+ G4VPrimitiveScorer* Energy = new GENERALSCORERS::PSEnergy("StripEnergy","MUST2Telescope", 0);
+ G4VPrimitiveScorer* TOF = new GENERALSCORERS::PSTOF("StripTime","MUST2Telescope", 0);
+
+ G4VPrimitiveScorer* StripPositionX = new PSStripNumberX("StripNumberX", 0, SiliconFace, 128);
+ G4VPrimitiveScorer* StripPositionY = new PSStripNumberY("StripNumberY", 0, SiliconFace, 128);
+
+ G4VPrimitiveScorer* InteractionCoordinatesX = new GENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","MUST2Telescope", 0);
+ G4VPrimitiveScorer* InteractionCoordinatesY = new GENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","MUST2Telescope", 0);
+ G4VPrimitiveScorer* InteractionCoordinatesZ = new GENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","MUST2Telescope", 0);
+ G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new GENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","MUST2Telescope", 0);
+ G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new GENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","MUST2Telescope", 0) ;
+
+
+ //and register it to the multifunctionnal detector
+ m_StripScorer->RegisterPrimitive(DetNbr);
+ m_StripScorer->RegisterPrimitive(Energy);
+ m_StripScorer->RegisterPrimitive(TOF);
+ m_StripScorer->RegisterPrimitive(StripPositionX);
+ m_StripScorer->RegisterPrimitive(StripPositionY);
+ m_StripScorer->RegisterPrimitive(InteractionCoordinatesX);
+ m_StripScorer->RegisterPrimitive(InteractionCoordinatesY);
+ m_StripScorer->RegisterPrimitive(InteractionCoordinatesZ);
+ m_StripScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
+ m_StripScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);
+
+ // SiLi Associate Scorer
+ m_SiLiScorer = new G4MultiFunctionalDetector("MUST2_SiLiScorer");
+ G4VPrimitiveScorer* SiLiEnergy = new GENERALSCORERS::PSEnergy("SiLiEnergy","MUST2Telescope", 0) ;
+ G4VPrimitiveScorer* SiLiPadNbr = new PSPadOrCristalNumber("SiLiPadNbr",0) ;
+ m_SiLiScorer->RegisterPrimitive(SiLiEnergy);
+ m_SiLiScorer->RegisterPrimitive(SiLiPadNbr);
+
+ // CsI Associate Scorer
+ m_CsIScorer = new G4MultiFunctionalDetector("MUST2_CsIScorer");
+ G4VPrimitiveScorer* CsIEnergy = new GENERALSCORERS::PSEnergy("CsIEnergy","MUST2Telescope", 0) ;
+ G4VPrimitiveScorer* CsICristalNbr = new PSPadOrCristalNumber("CsICristalNbr",0) ;
+ m_CsIScorer->RegisterPrimitive(CsIEnergy) ;
+ m_CsIScorer->RegisterPrimitive(CsICristalNbr) ;
+
+ // Add All Scorer to the Global Scorer Manager
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_SiLiScorer) ;
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ;
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void MUST2Array::InitializeMaterial()
- {
-
- ////////////////////////////////////////////////////////////////
-/////////////////Element Definition ///////////////////////////
-////////////////////////////////////////////////////////////////
- G4String symbol ;
- G4double density = 0. , a = 0, z = 0 ;
- G4int ncomponents = 0, natoms = 0 ;
-
- G4Element* H = new G4Element("Hydrogen" , symbol = "H" , z = 1 , a = 1.01 * g / mole);
- G4Element* C = new G4Element("Carbon" , symbol = "C" , z = 6 , a = 12.011 * g / mole);
- G4Element* N = new G4Element("Nitrogen" , symbol = "N" , z = 7 , a = 14.01 * g / mole);
- G4Element* O = new G4Element("Oxigen" , symbol = "O" , z = 8 , a = 16.00 * g / mole);
- G4Element* I = new G4Element("Iode" , symbol = "I" , z = 53 , a = 126.9 * g / mole);
- G4Element* Cs = new G4Element("Cesium" , symbol = "Cs" , z = 55 , a = 132.9 * g / mole);
-
- G4Element* Co = new G4Element("Cobalt" , symbol = "Co" , z = 27 , a = 58.933 * g / mole);
- G4Element* Cr = new G4Element("Cromium" , symbol = "Cr" , z = 24 , a = 51.996 * g / mole);
- G4Element* Ni = new G4Element("Nickel" , symbol = "Ni" , z = 28 , a = 58.69 * g / mole);
- G4Element* Fe = new G4Element("Iron" , symbol = "Fe" , z = 26 , a = 55.847 * g / mole);
- G4Element* W = new G4Element("Tungsten" , symbol = "W" , z = 74 , a = 183.5 * g / mole);
-
-////////////////////////////////////////////////////////////////
-/////////////////Material Definition ///////////////////////////
-////////////////////////////////////////////////////////////////
-
- // Si
- a = 28.0855 * g / mole;
- density = 2.321 * g / cm3;
- m_MaterialSilicon = new G4Material("Si", z = 14., a, density);
-
- // Al
- density = 2.702 * g / cm3;
- a = 26.98 * g / mole;
- m_MaterialAluminium = new G4Material("Aluminium", z = 13., a, density);
-
- // Iron
- density = 7.874 * g / cm3;
- a = 55.847 * g / mole;
- m_MaterialIron = new G4Material("Iron", z = 26., a, density);
-
- // CsI
- density = 4.51 * g / cm3;
- m_MaterialCsI = new G4Material("CsI", density, ncomponents = 2);
- m_MaterialCsI->AddElement(Cs , natoms = 1);
- m_MaterialCsI->AddElement(I , natoms = 1);
-
- // Vacuum
- density = 0.000000001 * mg / cm3;
- m_MaterialVacuum = new G4Material("Vacuum", density, ncomponents = 2);
- m_MaterialVacuum->AddElement(N, .7);
- m_MaterialVacuum->AddElement(O, .3);
-
- // Mylar
- density = 1.397 * g / cm3;
- m_MaterialMyl = new G4Material("Mylar", density, ncomponents = 3);
- m_MaterialMyl->AddElement(C, natoms = 10);
- m_MaterialMyl->AddElement(H, natoms = 8);
- m_MaterialMyl->AddElement(O, natoms = 4);
-
- // Havar
- m_MaterialHarvar = new G4Material("Havar", 8.3*g / cm3, 5);
- m_MaterialHarvar->AddElement(Co , 42);
- m_MaterialHarvar->AddElement(Cr , 20);
- m_MaterialHarvar->AddElement(Ni , 13);
- m_MaterialHarvar->AddElement(Fe , 19);
- m_MaterialHarvar->AddElement(W , 1);
-
-
- }
+void MUST2Array::InitializeMaterial(){
+
+ ////////////////////////////////////////////////////////////////
+ /////////////////Element Definition ///////////////////////////
+ ////////////////////////////////////////////////////////////////
+ G4String symbol ;
+ G4double density = 0. , a = 0, z = 0 ;
+ G4int ncomponents = 0, natoms = 0 ;
+
+ G4Element* H = new G4Element("Hydrogen" , symbol = "H" , z = 1 , a = 1.01 * g / mole);
+ G4Element* C = new G4Element("Carbon" , symbol = "C" , z = 6 , a = 12.011 * g / mole);
+ G4Element* N = new G4Element("Nitrogen" , symbol = "N" , z = 7 , a = 14.01 * g / mole);
+ G4Element* O = new G4Element("Oxigen" , symbol = "O" , z = 8 , a = 16.00 * g / mole);
+ G4Element* I = new G4Element("Iode" , symbol = "I" , z = 53 , a = 126.9 * g / mole);
+ G4Element* Cs = new G4Element("Cesium" , symbol = "Cs" , z = 55 , a = 132.9 * g / mole);
+
+ G4Element* Co = new G4Element("Cobalt" , symbol = "Co" , z = 27 , a = 58.933 * g / mole);
+ G4Element* Cr = new G4Element("Cromium" , symbol = "Cr" , z = 24 , a = 51.996 * g / mole);
+ G4Element* Ni = new G4Element("Nickel" , symbol = "Ni" , z = 28 , a = 58.69 * g / mole);
+ G4Element* Fe = new G4Element("Iron" , symbol = "Fe" , z = 26 , a = 55.847 * g / mole);
+ G4Element* W = new G4Element("Tungsten" , symbol = "W" , z = 74 , a = 183.5 * g / mole);
+
+ ////////////////////////////////////////////////////////////////
+ /////////////////Material Definition ///////////////////////////
+ ////////////////////////////////////////////////////////////////
+
+ // Si
+ a = 28.0855 * g / mole;
+ density = 2.321 * g / cm3;
+ m_MaterialSilicon = new G4Material("Si", z = 14., a, density);
+
+ // Al
+ density = 2.702 * g / cm3;
+ a = 26.98 * g / mole;
+ m_MaterialAluminium = new G4Material("Aluminium", z = 13., a, density);
+
+ // Iron
+ density = 7.874 * g / cm3;
+ a = 55.847 * g / mole;
+ m_MaterialIron = new G4Material("Iron", z = 26., a, density);
+
+ // CsI
+ density = 4.51 * g / cm3;
+ m_MaterialCsI = new G4Material("CsI", density, ncomponents = 2);
+ m_MaterialCsI->AddElement(Cs , natoms = 1);
+ m_MaterialCsI->AddElement(I , natoms = 1);
+
+ // Vacuum
+ density = 0.000000001 * mg / cm3;
+ m_MaterialVacuum = new G4Material("Vacuum", density, ncomponents = 2);
+ m_MaterialVacuum->AddElement(N, .7);
+ m_MaterialVacuum->AddElement(O, .3);
+
+ // Mylar
+ density = 1.397 * g / cm3;
+ m_MaterialMyl = new G4Material("Mylar", density, ncomponents = 3);
+ m_MaterialMyl->AddElement(C, natoms = 10);
+ m_MaterialMyl->AddElement(H, natoms = 8);
+ m_MaterialMyl->AddElement(O, natoms = 4);
+
+ // Havar
+ m_MaterialHarvar = new G4Material("Havar", 8.3*g / cm3, 5);
+ m_MaterialHarvar->AddElement(Co , 42);
+ m_MaterialHarvar->AddElement(Cr , 20);
+ m_MaterialHarvar->AddElement(Ni , 13);
+ m_MaterialHarvar->AddElement(Fe , 19);
+ m_MaterialHarvar->AddElement(W , 1);
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4RotationMatrix* Rotation(double tetaX, double tetaY, double tetaZ)
-{
- double PI = 3.141592653589793238;
- double radX = tetaX * PI / 180.;
- double radY = tetaY * PI / 180.;
- double radZ = tetaZ * PI / 180.;
-
- G4ThreeVector col1 = G4ThreeVector(cos(radZ) * cos(radY),
- -sin(radZ) * cos(radY),
- -sin(radY));
- G4ThreeVector col2 = G4ThreeVector(sin(radZ) * cos(radX) - cos(radZ) * sin(radY) * sin(radX),
- cos(radZ) * cos(radX) + sin(radZ) * sin(radY) * sin(radX),
- -cos(radY) * sin(radX));
- G4ThreeVector col3 = G4ThreeVector(sin(radZ) * sin(radX) + cos(radZ) * sin(radY) * sin(radX),
- cos(radZ) * sin(radX) - sin(radZ) * sin(radY) * cos(radX),
- cos(radY) * cos(radX));
-
- return (new G4RotationMatrix(col1, col2, col3));
+G4RotationMatrix* Rotation(double tetaX, double tetaY, double tetaZ){
+ double PI = 3.141592653589793238;
+ double radX = tetaX * PI / 180.;
+ double radY = tetaY * PI / 180.;
+ double radZ = tetaZ * PI / 180.;
+
+ G4ThreeVector col1 = G4ThreeVector(cos(radZ) * cos(radY),
+ -sin(radZ) * cos(radY),
+ -sin(radY));
+ G4ThreeVector col2 = G4ThreeVector(sin(radZ) * cos(radX) - cos(radZ) * sin(radY) * sin(radX),
+ cos(radZ) * cos(radX) + sin(radZ) * sin(radY) * sin(radX),
+ -cos(radY) * sin(radX));
+ G4ThreeVector col3 = G4ThreeVector(sin(radZ) * sin(radX) + cos(radZ) * sin(radY) * sin(radX),
+ cos(radZ) * sin(radX) - sin(radZ) * sin(radY) * cos(radX),
+ cos(radY) * cos(radX));
+
+ return (new G4RotationMatrix(col1, col2, col3));
}
diff --git a/NPSimulation/Paris/ParisCluster.cc b/NPSimulation/Paris/ParisCluster.cc
index 2ee93847d..43f265315 100644
--- a/NPSimulation/Paris/ParisCluster.cc
+++ b/NPSimulation/Paris/ParisCluster.cc
@@ -212,15 +212,16 @@ void ParisCluster::VolumeMaker(G4int DetecNumber,
false ,
0);
- logicParisCluster->SetVisAttributes(G4VisAttributes::Invisible);
- if (m_non_sensitive_part_visiualisation) logicParisCluster->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90)));
-
+ G4VisAttributes* ClusterVisAtt = new G4VisAttributes(G4Colour(0.90, 0.90, 0.90));
+ ClusterVisAtt->SetForceWireframe(true);
+ logicParisCluster->SetVisAttributes(ClusterVisAtt);
// Phoswich volume construction (empty volume)
G4ThreeVector positionPhoSwStage = G4ThreeVector(LaBr3Face+InterSpace, LaBr3Face+InterSpace, 0.);
G4Box* solidPhoSwStage = new G4Box(Name+"_PhoSwStage", 0.5*LaBr3Face, 0.5*LaBr3Face, 0.5*Length);
G4LogicalVolume* logicPhoSwStage = new G4LogicalVolume(solidPhoSwStage, Vacuum, "logicPhoSwStage", 0, 0, 0);
-
+ logicPhoSwStage->SetVisAttributes(G4VisAttributes::Invisible);
+
new G4PVPlacement(0,
positionPhoSwStage,
logicPhoSwStage,
@@ -230,7 +231,7 @@ void ParisCluster::VolumeMaker(G4int DetecNumber,
// 0);
positionPhoSwStage = G4ThreeVector(0., LaBr3Face+InterSpace, 0.);
- G4PVPlacement(0,
+ new G4PVPlacement(0,
positionPhoSwStage,
logicPhoSwStage,
Name+"_PhoSwStage",
@@ -323,7 +324,7 @@ void ParisCluster::VolumeMaker(G4int DetecNumber,
logicLaBr3Stage->SetSensitiveDetector(m_LaBr3StageScorer);
// Visualisation of LaBr3Stage Strip
- G4VisAttributes* LaBr3VisAtt = new G4VisAttributes(G4Colour(0., 0., 1.));
+ G4VisAttributes* LaBr3VisAtt = new G4VisAttributes(G4Colour(0., 0.5, 1.));
logicLaBr3Stage->SetVisAttributes(LaBr3VisAtt);
// CsI or NaI
@@ -345,7 +346,7 @@ void ParisCluster::VolumeMaker(G4int DetecNumber,
logicCsIStage->SetSensitiveDetector(m_CsIStageScorer);
// Visualisation of CsIStage Strip
- G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1., 0., 0.));
+ G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1., 0.25, 0.));
logicCsIStage->SetVisAttributes(CsIVisAtt);
}
diff --git a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
index 0a421521c..7cc5671b0 100644
--- a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
+++ b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
@@ -126,6 +126,9 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
+ // Generate the excitation energy if a distribution is given
+ m_Reaction->ShootRandomExcitationEnergy();
+
G4ParticleDefinition* HeavyName
= G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation4()*MeV);
--
GitLab