diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index 6ad7d6b9ab381ecddb16dc446eb828093a30eeb4..f71dd7e45f9a83f1255d41943076d7e6dae2e507 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 80
- EnergyHigh= 80
+ EnergyLow= 20
+ EnergyHigh= 20
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 90
+ HalfOpenAngleMax= 5
x0= 0
y0= 0
z0= 0
diff --git a/NPAnalysis/Example1/Analysis.cxx b/NPAnalysis/Example1/Analysis.cxx
index febc1f59403b6a9c82320d04ffa931c7a5d09738..3581d19f045e15c2ca2252dcf050a4ae3d5f268d 100644
--- a/NPAnalysis/Example1/Analysis.cxx
+++ b/NPAnalysis/Example1/Analysis.cxx
@@ -63,7 +63,7 @@ void Analysis::Init(){
He3CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
He3Al = EnergyLoss("He3_Al.G4table","G4Table",10);
He3Si = EnergyLoss("He3_Si.G4table","G4Table",10);
- Li11CD2 = EnergyLoss("Li11[0.0]_CD2.G4table","G4Table",100);
+ //Li11CD2 = EnergyLoss("Li11[0.0]_CD2.G4table","G4Table",100);
}
////////////////////////////////////////////////////////////////////////////////
@@ -79,7 +79,7 @@ void Analysis::TreatEvent(){
TVector3 BeamDirection = Initial->GetBeamDirection();
// Beam energy is measured using F3 and F2 plastic TOF
double BeamEnergy = Rand.Gaus(Initial->GetIncidentInitialKineticEnergy(),4.5);
- BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
+ //BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
He10Reaction->SetBeamEnergy(BeamEnergy);
//////////////////////////// LOOP on MUST2 + SSSD Hit //////////////////
diff --git a/NPAnalysis/Example1/Analysis.h b/NPAnalysis/Example1/Analysis.h
index d7f368ed30c33f0fd6b14a0d33b0aba1db929bc9..10215aa6d91645e0ed52fde8a5c213cc9c1e2b43 100644
--- a/NPAnalysis/Example1/Analysis.h
+++ b/NPAnalysis/Example1/Analysis.h
@@ -22,7 +22,7 @@
* *
*****************************************************************************/
#include"NPVAnalysis.h"
-#include"TAnnularS1Physics.h"
+#include"THiraPhysics.h"
#include "TMust2Physics.h"
#include "TSSSDPhysics.h"
#include "TInitialConditions.h"
@@ -76,6 +76,5 @@ class Analysis: public NPL::VAnalysis{
TMust2Physics* M2;
TSSSDPhysics* SSSD;
TInitialConditions* Initial;
-
};
#endif
diff --git a/NPAnalysis/Hira/Analysis.cxx b/NPAnalysis/Hira/Analysis.cxx
index 9f03961962f9d140a92d65de5f3fb5a0ca457bf9..c95877e97133df7a28e433a82ab171a09f116875 100644
--- a/NPAnalysis/Hira/Analysis.cxx
+++ b/NPAnalysis/Hira/Analysis.cxx
@@ -33,9 +33,9 @@ int main(int argc, char** argv)
// Instantiate some Reaction
NPL::Reaction* TransfertReaction = new Reaction ;
- //TransfertReaction -> ReadConfigurationFile("34Ar_pd.reaction") ;
+ TransfertReaction -> ReadConfigurationFile("34Ar_pd.reaction") ;
//TransfertReaction -> ReadConfigurationFile("46Ar_pd.reaction") ;
- TransfertReaction -> ReadConfigurationFile("11Be_d3He.reaction") ;
+ //TransfertReaction -> ReadConfigurationFile("11Be_d3He.reaction") ;
//Get Detector pointer :
THiraPhysics* Hira = (THiraPhysics*) myDetector -> GetDetector("HIRAArray") ;
diff --git a/NPAnalysis/Hira/RunToTreat.txt b/NPAnalysis/Hira/RunToTreat.txt
index df496a3640190bc6ca8cd0ee12a35713d1587872..3a48e00e5e11c3504a43a2c7a5a18b4ae360e136 100644
--- a/NPAnalysis/Hira/RunToTreat.txt
+++ b/NPAnalysis/Hira/RunToTreat.txt
@@ -4,7 +4,7 @@ RootFileName
%/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/46Ar_pd_gs.root
%/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/46Ar_pd_1st.root
- /Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/11Be_d3He.root
+ /Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/test.root
diff --git a/NPAnalysis/newHira/Analysis.cxx b/NPAnalysis/newHira/Analysis.cxx
index a4384b8fa057454681917a5ba0bd8e7de8074ca8..3581d19f045e15c2ca2252dcf050a4ae3d5f268d 100644
--- a/NPAnalysis/newHira/Analysis.cxx
+++ b/NPAnalysis/newHira/Analysis.cxx
@@ -6,12 +6,13 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: *
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
* *
* Creation Date : march 2025 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
+ * Class describing the property of an Analysis object *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -24,6 +25,8 @@ using namespace std;
#include"NPAnalysisFactory.h"
#include"NPDetectorManager.h"
#include"NPOptionManager.h"
+#include"RootOutput.h"
+#include"RootInput.h"
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -33,148 +36,172 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
+ M2= (TMust2Physics*) m_DetectorManager->GetDetector("MUST2Array");
+ SSSD= (TSSSDPhysics*) m_DetectorManager->GetDetector("SSSD");
+ Initial=new TInitialConditions();
InitOutputBranch();
InitInputBranch();
-
- Hira= (THiraPhysics*) m_DetectorManager -> GetDetector("HIRAArray");
- LightCD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
- LightAl = EnergyLoss("proton_Al.G4table","G4Table",100);
- LightSi = EnergyLoss("proton_Si.G4table","G4Table",100);
- //BeamCD2 = EnergyLoss("Rb88_CD2.G4table","G4Table",100);
- myReaction = new NPL::Reaction();
- myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
- TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
- OriginalBeamEnergy = myReaction->GetBeamEnergy();
Rand = TRandom3();
-
- E_ThinSi = 0;
- E_ThickSi = 0;
- E_CsI = 0;
- ELab = 0;
- ThetaLab = 0;
- PhiLab = 0;
- ThetaLab_simu = 0;
- ThetaCM = 0;
- ExcitationEnergy = 0;
- X = 0;
- Y = 0;
- Z = 0;
- TelescopeNumber = 0;
- EnergyThreshold = 1;
-
+ He10Reaction= new NPL::Reaction();
+ He10Reaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ DetectorNumber = 0 ;
+ ThetaNormalTarget = 0 ;
+ ThetaM2Surface = 0;
+ X_M2 = 0 ;
+ Y_M2 = 0 ;
+ Z_M2 = 0 ;
+ Si_E_M2 = 0 ;
+ CsI_E_M2 = 0 ;
+ E_SSSD = 0 ;
+ Energy = 0;
+ E_M2 = 0;
+ Si_X_M2 = 0;
+ Si_Y_M2 = 0;
+ ZTarget = 0;
+ TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
+ // Energy loss table: the G4Table are generated by the simulation
+ He3CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
+ He3Al = EnergyLoss("He3_Al.G4table","G4Table",10);
+ He3Si = EnergyLoss("He3_Si.G4table","G4Table",10);
+ //Li11CD2 = EnergyLoss("Li11[0.0]_CD2.G4table","G4Table",100);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
// Reinitiate calculated variable
ReInitValue();
- if(Hira -> ThickSi_E.size() == 1){
- for(int countHira = 0 ; countHira < Hira->ThickSi_E.size() ; countHira++){
-
- TelescopeNumber = Hira->TelescopeNumber[countHira];
-
- TargetThickness = m_DetectorManager->GetTargetThickness();
+ // Get the Init information on beam position and energy
+ // and apply by hand the experimental resolution
+ // This is because the beam diagnosis are not simulated
+ // PPAC position resolution on target is assumed to be 1mm
+ double XTarget = Rand.Gaus(Initial->GetIncidentPositionX(),1);
+ double YTarget = Rand.Gaus(Initial->GetIncidentPositionY(),1);
+ TVector3 BeamDirection = Initial->GetBeamDirection();
+ // Beam energy is measured using F3 and F2 plastic TOF
+ double BeamEnergy = Rand.Gaus(Initial->GetIncidentInitialKineticEnergy(),4.5);
+ //BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
+
+ He10Reaction->SetBeamEnergy(BeamEnergy);
+ //////////////////////////// LOOP on MUST2 + SSSD Hit //////////////////
+ for(unsigned int countSSSD = 0 ; countSSSD < SSSD->Energy.size() ; countSSSD++){
+ for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
+ /************************************************/
+ //Part 0 : Get the usefull Data
+ // MUST2
+ int X = M2->Si_X[countMust2];
+ int Y = M2->Si_Y[countMust2];
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+ Si_X_M2 = X ;
+ Si_Y_M2 = Y ;
+ //SSSD
+ int SiNumber = SSSD->DetectorNumber[countSSSD];
+
+ /************************************************/
+ // Matching between Thin Si and MUST2, and Forward Telescope Only
+ if(TelescopeNumber==SiNumber && TelescopeNumber<5){
+ DetectorNumber = TelescopeNumber ;
+ /************************************************/
+ // Part 1 : Impact Angle
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+ if(XTarget>-1000 && YTarget>-1000){
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
+ Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
+ Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
+ }
- X = Hira->GetPositionOfInteraction(0).X();
- Y = Hira->GetPositionOfInteraction(0).Y();
- Z = Hira->GetPositionOfInteraction(0).Z();
+ else{
+ BeamDirection = TVector3(-1000,-1000,-1000);
+ ThetaM2Surface = -1000 ;
+ ThetaNormalTarget = -1000 ;
+ }
- TVector3 PositionOnHira = TVector3(X,Y,Z);
- TVector3 ZUnit = TVector3(0,0,1);
+ /************************************************/
+
+ /************************************************/
+
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2= M2->CsI_E[countMust2];
+ E_SSSD = SSSD->Energy[countSSSD];
+
+ // if CsI
+ if(CsI_E_M2>0 ){
+ // The energy in CsI is calculate form dE/dx Table because
+ // 20um resolution is poor
+ Energy =
+ He3Si.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer,
+ ThetaM2Surface, 0.01*MeV,
+ 450.*MeV,0.001*MeV ,1000);
+ E_M2=CsI_E_M2;
+ }
- double X_target = InitialConditions->GetIncidentPositionX();
- double Y_target = InitialConditions->GetIncidentPositionY();
- double Z_target = InitialConditions->GetIncidentPositionZ();
+ else
+ Energy = Si_E_M2;
- TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
- TVector3 HitDirection = PositionOnHira - PositionOnTarget;
- TVector3 HitDirectionUnit = HitDirection.Unit();
+ E_M2 += Si_E_M2;
- TVector3 BeamDirection = InitialConditions->GetBeamDirection();
- double XBeam = BeamDirection.X();
- double YBeam = BeamDirection.Y();
- double ZBeam = BeamDirection.Z();
+ // Evaluate energy using the thickness
+ ELab = He3Al.EvaluateInitialEnergy( Energy ,2*0.4*micrometer , ThetaM2Surface);
+ ELab = He3Si.EvaluateInitialEnergy( ELab ,20*micrometer , ThetaM2Surface);
+ ELab = He3Al.EvaluateInitialEnergy( ELab ,0.4*micrometer , ThetaM2Surface);
+ // Target Correction
+ ELab = He3CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ /************************************************/
- ThetaLab = BeamDirection.Angle(HitDirection);
- ThetaLab = ThetaLab/deg;
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = He10Reaction -> ReconstructRelativistic( ELab , ThetaLab );
+ /************************************************/
- PhiLab = HitDirection.Phi()/deg;
- E_ThickSi = Hira->ThickSi_E[countHira];
- E_ThinSi = Hira->ThinSi_E[countHira];
- //ELab = E_ThinSi + E_ThickSi;
- if(Hira->CsI_E.size() == 1){
- for(int countCsI =0; countCsI<Hira->CsI_E.size(); countCsI++){
- E_CsI = Hira->CsI_E[countCsI];
- //ELab += E_CsI;
- if(E_CsI>EnergyThreshold) ELab = E_ThinSi + E_ThickSi + E_CsI;
- }
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = He10Reaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ ThetaLab=ThetaLab/deg;
+ /************************************************/
+
}
- else ELab = -1000;
-
- // ********************** Angle in the CM frame *****************************
- myReaction -> SetNuclei3(ELab, ThetaLab*deg);
- ThetaCM = myReaction->GetThetaCM()/deg;
- ExcitationEnergy = myReaction->GetExcitation4();
- }
- }
+ } //end loop SSSD
+ }//end loop MUST2
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
+
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
-
- RootOutput::getInstance()->GetTree()->Branch( "ThetaLab" , &ThetaLab , "ThetaLab/D" );
- RootOutput::getInstance()->GetTree()->Branch( "PhiLab" , &PhiLab , "PhiLab/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM,"ThetaCM/D") ;
- RootOutput::getInstance()->GetTree()->Branch( "E_ThinSi" , &E_ThinSi , "E_ThinSi/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "E_ThickSi" , &E_ThickSi , "E_ThickSi/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "E_CsI" , &E_CsI , "E_CsI/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "ELab" , &ELab , "ELab/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy,"ExcitationEnergy/D") ;
- RootOutput::getInstance()->GetTree()->Branch( "X" , &X , "X/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "Y" , &Y , "Y/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "Z" , &Z , "Z/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "TelescopeNumber" , &TelescopeNumber , "TelescopeNumber/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&InteractionCoordinates);
- RootOutput::getInstance()->GetTree()->Branch("InitialConditions","TInitialConditions",&InitialConditions);
-
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
}
+
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
- //TInteractionCoordinate
- InteractionCoordinates = new TInteractionCoordinates();
- RootInput::getInstance()->GetChain()->SetBranchStatus("InteractionCoordinates",true);
- RootInput::getInstance()->GetChain()->SetBranchStatus("fDetected_*",true);
- RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&InteractionCoordinates);
-
- //TInitialConditions
- InitialConditions = new TInitialConditions();
- RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
- RootInput::getInstance()->GetChain()->SetBranchStatus("fIC_*",true);
- RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions);
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
}
-////////////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
- E_ThinSi = -1000;
- E_ThickSi = -1000;
- E_CsI = -1000;
- ELab = -1000;
- ThetaLab = -1000;
- PhiLab = -1000;
- ThetaCM = -1000;
- ExcitationEnergy = -1000;
- X = -1000;
- Y = -1000;
- Z = -1000;
- TelescopeNumber = -1;
+ Ex = -1000 ;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
}
-
////////////////////////////////////////////////////////////////////////////////
-// Construct Method to be pass to the AnalysisFactory //
+// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
return (NPL::VAnalysis*) new Analysis();
diff --git a/NPAnalysis/newHira/Analysis.h b/NPAnalysis/newHira/Analysis.h
index 463519dcafd043b5f18f20a144c975c8c260cede..10215aa6d91645e0ed52fde8a5c213cc9c1e2b43 100644
--- a/NPAnalysis/newHira/Analysis.h
+++ b/NPAnalysis/newHira/Analysis.h
@@ -22,17 +22,13 @@
* *
*****************************************************************************/
#include"NPVAnalysis.h"
-#include"NPEnergyLoss.h"
-#include"NPReaction.h"
-#include"RootOutput.h"
-#include"RootInput.h"
-#include "THiraPhysics.h"
+#include"THiraPhysics.h"
+#include "TMust2Physics.h"
+#include "TSSSDPhysics.h"
#include "TInitialConditions.h"
-#include "TInteractionCoordinates.h"
-#include <TRandom3.h>
-#include <TVector3.h>
-#include <TMath.h>
-
+#include "NPEnergyLoss.h"
+#include "NPReaction.h"
+#include "TRandom3.h"
class Analysis: public NPL::VAnalysis{
public:
Analysis();
@@ -42,40 +38,43 @@ class Analysis: public NPL::VAnalysis{
void Init();
void TreatEvent();
void End();
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
- void InitOutputBranch();
- void InitInputBranch();
- void ReInitValue();
- static NPL::VAnalysis* Construct();
-
private:
- double ExcitationEnergy;
- double ELab;
- double E_ThinSi;
- double E_ThickSi;
- double E_CsI;
- double PhiLab;
- double ThetaLab;
- double ThetaLab_simu;
- double ThetaCM;
- double X,Y,Z;
- double TelescopeNumber;
- double EnergyThreshold;
+ double Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+ NPL::Reaction* He10Reaction;
+
+ // intermediate variable
+ TRandom3 Rand;
+ int DetectorNumber;
+ double ThetaNormalTarget;
+ double ThetaM2Surface;
+ double X_M2;
+ double Y_M2;
+ double Z_M2;
+ double Si_E_M2;
+ double CsI_E_M2;
+ double E_SSSD;
+ double Energy = 0;
+ double E_M2 = 0;
+ double Si_X_M2 = 0;
+ double Si_Y_M2 = 0;
+ double ZTarget = 0;
+ double TargetThickness;
- NPL::Reaction* myReaction;
- TInitialConditions* InitialConditions ;
- TInteractionCoordinates* InteractionCoordinates;
+ NPL::EnergyLoss He3CD2 ;
+ NPL::EnergyLoss He3Al ;
+ NPL::EnergyLoss He3Si ;
+ NPL::EnergyLoss Li11CD2 ;
- // Energy loss table: the G4Table are generated by the simulation
- EnergyLoss LightCD2;
- EnergyLoss LightAl;
- EnergyLoss LightSi;
- EnergyLoss BeamCD2;
-
- double TargetThickness ;
- double OriginalBeamEnergy ;
- TRandom3 Rand ;
- // Hira object
- THiraPhysics* Hira;
+ TMust2Physics* M2;
+ TSSSDPhysics* SSSD;
+ TInitialConditions* Initial;
};
#endif
diff --git a/NPSimulation/Core/EventGeneratorIsotropic.cc b/NPSimulation/Core/EventGeneratorIsotropic.cc
index ce6ae242b39cf959110ab17e30c189ec843d1aaf..7c018a95d3d386ac39b8c6deb5aadf6173f6bdbd 100644
--- a/NPSimulation/Core/EventGeneratorIsotropic.cc
+++ b/NPSimulation/Core/EventGeneratorIsotropic.cc
@@ -165,6 +165,7 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path,int){
else if(m_particleName=="alpha") { m_particleName="4He" ; check_ExcitationEnergy = true ;}
else if(m_particleName=="gamma") { check_ExcitationEnergy = true ;}
else if(m_particleName=="neutron") { check_ExcitationEnergy = true ;}
+ else { check_ExcitationEnergy = true ;}
}
diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc
index 4c0e7cd6924f9fdd2fedb237fe4cb5002a66ab2e..3687fb48041db77a28c1a9563c61cc4bcd2b81e1 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -123,6 +123,7 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
m_Material[Name]=material;
return material;
}
+
else if(Name == "Kapton"){
G4Material* material = new G4Material(Name, 1.39*g/cm3,3);
@@ -316,7 +317,38 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
m_Material[Name]=material;
return material;
}
+
+ else if(Name == "P10_1atm"){
+ G4Material* material = new G4Material(Name, 1.74*mg/cm3,3); //@ 0K, 1 atm
+ material->AddElement(GetElementFromLibrary("Ar"),0.9222);
+ material->AddElement(GetElementFromLibrary("C"),0.0623);
+ material->AddElement(GetElementFromLibrary("H"),0.0155);
+ m_Material[Name]=material;
+ return material;
+ }
+
+ else if(Name == "P10"){
+ G4Material* material = new G4Material(Name, 0.57*mg/cm3,3); //@ 0K, 1/3 atm
+ material->AddElement(GetElementFromLibrary("Ar"),0.9222);
+ material->AddElement(GetElementFromLibrary("C"),0.0623);
+ material->AddElement(GetElementFromLibrary("H"),0.0155);
+ m_Material[Name]=material;
+ return material;
+ }
+ else if(Name == "Air_1atm"){ // 1 atm
+ G4Material* material = new G4Material("Air", 1.290*mg/cm3, 2);
+ material->AddElement(GetElementFromLibrary("N"), 0.7);
+ material->AddElement(GetElementFromLibrary("O"), 0.3);
+ }
+
+ else if(Name == "Air"){ // 1/3 atm
+ G4Material* material = new G4Material("Air", 1.290/3*mg/cm3, 2);
+ material->AddElement(GetElementFromLibrary("N"), 0.7);
+ material->AddElement(GetElementFromLibrary("O"), 0.3);
+ }
+
+
else{
G4cout << "ERROR: Material requested \""<< Name <<"\" is not available in the Material Librairy" << G4endl;
exit(1);
diff --git a/NPSimulation/Hira/Hira.cc b/NPSimulation/Hira/Hira.cc
index b01a6c7294bede52d28d51532082848ad7299b9b..00d7cd80a2f6c7ded0c4a4ff93554a43ca00c13e 100644
--- a/NPSimulation/Hira/Hira.cc
+++ b/NPSimulation/Hira/Hira.cc
@@ -530,9 +530,13 @@ void Hira::VolumeMaker(G4int DetectorNumber,
///////////////////////////////////////////////////
G4String NameCsI = "CsI"+DetNumber;
+ double X1 = (CsIXFront-CsIXBack)/2.;
+ double Y1 = (CsIYFront-CsIYBack)/2.;
+ double l = sqrt(pow(X1,2) + pow(Y1,2));
+
double pDz = 0.5*CsIThickness;
- double pTheta = -atan( (CsIXBack-CsIXFront)/(2*CsIThickness) );
- double pPhi = 0;
+ double pTheta = -atan( (l)/(CsIThickness) );
+ double pPhi = atan( X1/Y1 );
double pDy1 = 0.5*CsIYFront;
double pDx1 = 0.5*CsIXFront;
double pDx2 = 0.5*CsIXFront;
@@ -567,28 +571,33 @@ void Hira::VolumeMaker(G4int DetectorNumber,
m_LogicCluster,"Cluster",
m_logicMotherVolume,false,0);
- const G4double CsIXMiddle = CsIXFront + CsIThickness*tan(-pTheta);
+ const G4double CsIXMiddle = CsIXFront + (CsIThickness/2)*tan(-pTheta)*sin(pPhi);
+ const G4double CsIYMiddle = CsIYFront + (CsIThickness/2)*tan(-pTheta)*cos(pPhi);
const G4double DistInterCsIX = CsIXMiddle+DistInterCsI;
- const G4double DistInterCsIY = CsIYFront+DistInterCsI;
+ const G4double DistInterCsIY = CsIYMiddle+DistInterCsI;
G4ThreeVector Origin(-0.5*DistInterCsIX,-0.5*DistInterCsIY,0);
- G4RotationMatrix* rotM = new G4RotationMatrix;
- const G4double dangle = 180.*deg;
+ G4ThreeVector Pos;
+ const G4double dangle = 90.*deg;
// A cluster is a 2 by 2 aggregat of CsI crystal
unsigned int CsINbr = 1;
for(unsigned int i = 0 ; i < 2 ; i++){
- for(unsigned int j = 0 ; j <2 ; j++){
+ for(unsigned int j = 0 ; j < 2 ; j++){
+ G4RotationMatrix* rotM = new G4RotationMatrix;
unsigned int CrystalNbr = CsINbr++;
- if(i==0)rotM->rotateZ((i)*dangle);
- if(i==1)rotM->rotateZ((i+j)*dangle);
- G4ThreeVector Pos = Origin + G4ThreeVector(i*DistInterCsIX,j*DistInterCsIY,0);
-
+ if(i==0 && j==0)rotM->rotateZ(0);
+ if(i==1 && j==0)rotM->rotateZ(dangle);
+ if(i==0 && j==1)rotM->rotateZ(-dangle);
+ if(i==1 && j==1)rotM->rotateZ(2*dangle);
+ Pos = Origin + G4ThreeVector(i*DistInterCsIX,j*DistInterCsIY,0);
+
new G4PVPlacement(G4Transform3D(*rotM,Pos),
m_LogicCsICrystal,
"CsI_Cristal",
m_LogicCluster,
false,
CrystalNbr);
+ delete rotM;
}
}
}
diff --git a/NPSimulation/Hira/Hira.hh b/NPSimulation/Hira/Hira.hh
index f049bbe9f6f00d0cac8a8b4a55803001e6e1196f..112ae19a7be25acb566da699f46d5d6ca342cd95 100644
--- a/NPSimulation/Hira/Hira.hh
+++ b/NPSimulation/Hira/Hira.hh
@@ -71,8 +71,8 @@ namespace HIRA
const G4double CsIXFront = 33.*mm;
const G4double CsIXBack = 37.*mm;
const G4double CsIYFront = 33.*mm;
- const G4double CsIYBack = 33.*mm;
- const G4double DistInterCsI = 0.01*mm;
+ const G4double CsIYBack = 37.*mm;
+ const G4double DistInterCsI = 1.5*mm;
const G4double ClusterFaceFront = 7*cm;
const G4double ClusterFaceBack = 9*cm;