Skip to content
Snippets Groups Projects
Commit fba72618 authored by Pierre Morfouace's avatar Pierre Morfouace
Browse files

Updaintg pista simulation

parent 1e727c33
No related branches found
No related tags found
1 merge request!27Draft: [Epic] Preparation of the environement for the new GaseousDetectorScorers...
Hide whitespace changes
Inline Side-by-side
Showing
with 170 additions and 80 deletions
NPLib/Physics/TFissionConditions.h
+ 1
− 0
View file @ fba72618
......@@ -117,6 +117,7 @@ class TFissionConditions : public TObject{
int GetNeutronMultiplicity() const {return fFC_Neutron_Multiplicity;}//!
// emmitted particles
int GetFragmentMult() const {return fFC_Fragment_Z.size();};//!
string GetFragmentName (const int &i) const {return fFC_Fragment_Name[i];}//!
int GetFragmentZ (const int &i) const {return fFC_Fragment_Z[i];}//!
int GetFragmentA (const int &i) const {return fFC_Fragment_A[i];}//!
......
NPSimulation/Detectors/PISTA/PISTA.cc
+ 8
− 7
View file @ fba72618
......@@ -68,8 +68,8 @@ namespace PISTA_NS{
double TrapezoidBaseSmall = 41.0*mm;
double TrapezoidHeight = 57.7*mm;
double TrapezoidLength = 1*cm;
double FirstStageThickness = 300*um;
double SecondStageThickness = 1.5*mm;
double FirstStageThickness = 100*um;
double SecondStageThickness = 1.*mm;
double DistanceBetweenSi = 4*mm;
double FirstStageNbrOfStrips = 91;
double SecondStageNbrOfStrips = 57;
......@@ -118,8 +118,8 @@ void PISTA::AddDetector(G4ThreeVector A, G4ThreeVector B, G4ThreeVector C, G4Thr
TrapezoidBaseSmall = 41.0*mm;
TrapezoidHeight = 57.7*mm;
TrapezoidLength = 1*cm;
FirstStageThickness = 300*um;
SecondStageThickness = 1.5*mm;
FirstStageThickness = 100*um;
SecondStageThickness = 1.0*mm;
DistanceBetweenSi = 4*mm;
FirstStageNbrOfStrips = 91;
SecondStageNbrOfStrips = 57;
......@@ -157,8 +157,8 @@ void PISTA::AddDetector(double R, double Theta, double Phi){
TrapezoidBaseSmall = 41.0*mm;
TrapezoidHeight = 57.7*mm;
TrapezoidLength = 1*cm;
FirstStageThickness = 300*um;
SecondStageThickness = 1.5*mm;
FirstStageThickness = 100*um;
SecondStageThickness = 1.*mm;
DistanceBetweenSi = 4*mm;
FirstStageNbrOfStrips = 91;
SecondStageNbrOfStrips = 57;
......@@ -389,7 +389,8 @@ void PISTA::ReadSensitive(const G4Event* ){
int DetNbr = FirstStageScorer->GetDetectorWidth(i);
int StripFront = FirstStageScorer->GetStripWidth(i);
m_Event->SetPISTA_DE_DetectorNbr(DetNbr);
m_Event->SetPISTA_DE_StripNbr(StripFront);
if (m_version==1) m_Event->SetPISTA_DE_StripNbr(92-StripFront);
else m_Event->SetPISTA_DE_StripNbr(StripFront);
m_Event->SetPISTA_DE_StripEnergy(EnergyFront);
m_Event->SetPISTA_DE_StripTime(Time);
}
......
NPSimulation/Detectors/Vamos/Vamos.cc
+ 21
− 11
View file @ fba72618
......@@ -508,7 +508,7 @@ void Vamos::ReadConfiguration(NPL::InputParser parser) {
void Vamos::ConstructDetector(G4LogicalVolume* world) {
// Mother Volume of Vamos
G4double R = m_R + VamosVolumeThickness * 0.5;
/*G4double R = m_R + VamosVolumeThickness * 0.5;
G4double X = R * sin(m_Theta) * cos(Phi);
G4double Y = R * sin(m_Theta) * sin(Phi);
......@@ -549,10 +549,10 @@ void Vamos::ConstructDetector(G4LogicalVolume* world) {
G4LogicalVolume* MotherDetector = new G4LogicalVolume(
MotherDetectorSolid, VolumeMaterial, "MotherDetector", 0, 0, 0);
new G4PVPlacement(G4Transform3D(*Rot2, Det_pos2), MotherDetector,
*/
/*new G4PVPlacement(G4Transform3D(*Rot2, Det_pos2), MotherDetector,
"MotherDetector", MotherVolume, false, 0);
MotherDetector->SetVisAttributes(m_VisVolumeVamos);
MotherDetector->SetVisAttributes(m_VisVolumeVamos);*/
// Position the entry DCs and the magnets in the MotherVolume
......@@ -562,6 +562,7 @@ void Vamos::ConstructDetector(G4LogicalVolume* world) {
new G4PVPlacement(0, G4ThreeVector(0, 0, -VamosVolumeThickness * 0.5 + m_Z[1]),
BuildDC2(), "Entrance_DC2", MotherVolume, false, 2);*/
/*
new G4PVPlacement(
0,
G4ThreeVector(0, 0, (-VamosVolumeThickness + MagnetThickness) * 0.5 + 400),
......@@ -577,21 +578,30 @@ void Vamos::ConstructDetector(G4LogicalVolume* world) {
0,
G4ThreeVector(0, 0,
(-VamosVolumeThickness + MagnetThickness) * 0.5 + 1500 * mm),
BuildDipol(), "Vamos", MotherVolume, false, 0);
BuildDipol(), "Vamos", MotherVolume, false, 0);*/
// Position the system of detection at the end of Vamos in the sub Volume
/*new G4PVPlacement(0, m_PosCatcher, BuildBeamCatcher(), "BeamCatcher",
MotherDetector, false, 3);*/
G4RotationMatrix* RotVamos = new G4RotationMatrix();
RotVamos->rotateY(m_Theta);
G4double X1 = m_Z_TMW1 * sin(m_Theta) * cos(Phi);
G4double Y1 = m_Z_TMW1 * sin(m_Theta) * sin(Phi);
G4double Z1 = m_Z_TMW1 * cos(m_Theta);
G4ThreeVector Det_pos1 = G4ThreeVector(X1, Y1, Z1);
G4double X2 = m_Z_TMW2 * sin(m_Theta) * cos(Phi);
G4double Y2 = m_Z_TMW2 * sin(m_Theta) * sin(Phi);
G4double Z2 = m_Z_TMW2 * cos(m_Theta);
G4ThreeVector Det_pos2 = G4ThreeVector(X2, Y2, Z2);
new G4PVPlacement(
0,
G4ThreeVector(0, 0, -VamosVolumeThickness * 0.5 + TMW1_Thickness * 0.5),
BuildTMW1(), "TMW1", MotherVolume, false, 0);
G4Transform3D(*RotVamos, Det_pos1),
BuildTMW1(), "TMW1", world, false, 0);
new G4PVPlacement(
0,
G4ThreeVector(0, 0, -VamosVolumeThickness * 0.5 + TMW2_Thickness * 0.5 + m_Z_TMW2),
BuildTMW2(), "TMW2", MotherVolume, false, 0);
G4Transform3D(*RotVamos, Det_pos2),
BuildTMW2(), "TMW2", world, false, 0);
/* new G4PVPlacement(
......
Projects/PISTA/Analysis.cxx
+ 37
− 35
View file @ fba72618
......@@ -37,6 +37,7 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
PISTA= (TPISTAPhysics*) m_DetectorManager->GetDetector("PISTA");
FissionConditions = new TFissionConditions();
InitialConditions = new TInitialConditions();
ReactionConditions = new TReactionConditions();
InteractionCoordinates = new TInteractionCoordinates();
......@@ -46,7 +47,7 @@ void Analysis::Init(){
TargetThickness = m_DetectorManager->GetTargetThickness();
Transfer = new NPL::Reaction("238U(12C,12C)238U@1428");
Transfer = new NPL::Reaction("238U(12C,10Be)240Pu@1428");
// Energy loss table
//Be10C = EnergyLoss("EnergyLossTable/Be10_C.G4table","G4Table",100);
......@@ -57,20 +58,14 @@ void Analysis::Init(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
ReInitValue();
OriginalThetaLab = ReactionConditions->GetTheta(0);
OriginalElab = ReactionConditions->GetKineticEnergy(0);
OriginalBeamEnergy = ReactionConditions->GetBeamEnergy();
OriginalEx = ReactionConditions->GetExcitation4();
int mult = InteractionCoordinates->GetDetectedMultiplicity();
if(mult==1){
for(int i=0; i<mult; i++){
Xpista = InteractionCoordinates->GetDetectedPositionX(i);
Ypista = InteractionCoordinates->GetDetectedPositionY(i);
Zpista = InteractionCoordinates->GetDetectedPositionZ(i);
R = sqrt(Xpista*Xpista + Ypista*Ypista + Zpista*Zpista);
}
int FragmentMult = FissionConditions->GetFragmentMult();
for(int i=0; i<FragmentMult; i++){
double init_A = FissionConditions->GetFragmentA(i);
double init_Z = FissionConditions->GetFragmentZ(i);
init_FF_A.push_back(init_A);
init_FF_Z.push_back(init_Z);
}
XTarget = InitialConditions->GetIncidentPositionX();
YTarget = InitialConditions->GetIncidentPositionY();
ZTarget = InitialConditions->GetIncidentPositionZ();
......@@ -79,13 +74,11 @@ void Analysis::TreatEvent(){
TVector3 BeamPosition(XTarget,YTarget,ZTarget);
TVector3 PositionOnTarget(0,0,0);
//TVector3 PositionOnTarget(-1,0.5,0);
//TVector3 PositionOnTarget(Rand.Gaus(XTarget, 0.6/2.35), Rand.Gaus(YTarget, 0.6/2.35), 0);
//TVector3 PositionOnTarget(XTarget, YTarget, 0);
BeamEnergy = 1428.;//InitialConditions->GetIncidentInitialKineticEnergy();
BeamEnergy = U238C.Slow(BeamEnergy,TargetThickness*0.5,0);
//BeamEnergy = U238C.Slow(BeamEnergy,TargetThickness*0.5,0);
Transfer->SetBeamEnergy(BeamEnergy);
//Transfer->SetBeamEnergy(OriginalBeamEnergy);
//cout << PISTA->EventMultiplicity << endl;
if(PISTA->EventMultiplicity==1){
......@@ -121,45 +114,59 @@ void Analysis::TreatEvent(){
//Elab = Be10C.EvaluateInitialEnergy(Energy,TargetThickness*0.5,ThetaNormalTarget);
Elab = Energy;
Elab = C12C.EvaluateInitialEnergy(Energy,TargetThickness*0.5,ThetaNormalTarget);
OptimumEx = Transfer->ReconstructRelativistic(OriginalElab, OriginalThetaLab*deg);
Ex = Transfer->ReconstructRelativistic(Elab, ThetaLab);
ThetaCM = Transfer->EnergyLabToThetaCM(Elab, ThetaLab)/deg;
ThetaLab = ThetaLab/deg;
int mult = InteractionCoordinates->GetDetectedMultiplicity();
if(mult==2){
for(int j=0; j<mult; j++){
double Zint = InteractionCoordinates->GetDetectedPositionZ(j);
double A = InteractionCoordinates->GetA(j);
double Z = InteractionCoordinates->GetZ(j);
if(A!=10){
FF_A = A;
FF_Z = Z;
}
}
}
}
}
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("OriginalBeamEnergy",&OriginalBeamEnergy,"OriginalBeamEnergy/D");
RootOutput::getInstance()->GetTree()->Branch("OriginalEx",&OriginalEx,"OriginalEx/D");
RootOutput::getInstance()->GetTree()->Branch("BeamEnergy",&BeamEnergy,"BeamEnergy/D");
RootOutput::getInstance()->GetTree()->Branch("XTarget",&XTarget,"XTarget/D");
RootOutput::getInstance()->GetTree()->Branch("YTarget",&YTarget,"YTarget/D");
RootOutput::getInstance()->GetTree()->Branch("ZTarget",&ZTarget,"ZTarget/D");
RootOutput::getInstance()->GetTree()->Branch("OptimumEx",&OptimumEx,"OptimumEx/D");
RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
RootOutput::getInstance()->GetTree()->Branch("DeltaE",&DeltaE,"DeltaE/D");
RootOutput::getInstance()->GetTree()->Branch("Eres",&Eres,"Eres/D");
RootOutput::getInstance()->GetTree()->Branch("PID",&PID,"PID/D");
RootOutput::getInstance()->GetTree()->Branch("Elab",&Elab,"Elab/D");
RootOutput::getInstance()->GetTree()->Branch("OriginalElab",&OriginalElab,"OriginalElab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("R",&R,"R/D");
RootOutput::getInstance()->GetTree()->Branch("Xpista",&Xpista,"Xpista/D");
RootOutput::getInstance()->GetTree()->Branch("Ypista",&Ypista,"Ypista/D");
RootOutput::getInstance()->GetTree()->Branch("Zpista",&Zpista,"Zpista/D");
RootOutput::getInstance()->GetTree()->Branch("Xcalc",&Xcalc,"Xcalc/D");
RootOutput::getInstance()->GetTree()->Branch("Ycalc",&Ycalc,"Ycalc/D");
RootOutput::getInstance()->GetTree()->Branch("Zcalc",&Zcalc,"Zcalc/D");
RootOutput::getInstance()->GetTree()->Branch("Telescope",&Telescope,"Telescope/I");
RootOutput::getInstance()->GetTree()->Branch("FF_A",&FF_A,"FF_A/D");
RootOutput::getInstance()->GetTree()->Branch("FF_Z",&FF_Z,"FF_Z/D");
RootOutput::getInstance()->GetTree()->Branch("init_FF_Z",&init_FF_Z);
RootOutput::getInstance()->GetTree()->Branch("init_FF_A",&init_FF_A);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchStatus("FisisonConditions",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("fFC_*",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("FissionConditions",&FissionConditions);
RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("fIC_*",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions);
......@@ -173,32 +180,27 @@ void Analysis::InitInputBranch(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
OriginalBeamEnergy = -1000;
OriginalEx = -1000;
BeamEnergy = -1000;
OptimumEx = -1000;
Ex = -1000;
DeltaE = -1000;
Eres = -1000;
Elab = -1000;
OriginalElab = -1000;
OriginalThetaLab = -1000;
ThetaLab = -1000;
PhiLab = -1000;
ThetaCM = -1000;
XTarget = -1000;
YTarget = -1000;
ZTarget = -1000;
OriginalThetaLab = -1000;
R = -1000;
Xpista = -1000;
Ypista = -1000;
Zpista = -1000;
Xcalc = -1000;
Ycalc = -1000;
Zcalc = -1000;
PID = -1000;
Telescope = -1;
FF_A = -1;
FF_Z = -1;
init_FF_Z.clear();
init_FF_A.clear();
}
////////////////////////////////////////////////////////////////////////////////
......
Projects/PISTA/Analysis.h
+ 7
− 5
View file @ fba72618
......@@ -23,6 +23,7 @@
#include "NPVAnalysis.h"
#include "TPISTAPhysics.h"
#include "TFissionConditions.h"
#include "TInitialConditions.h"
#include "TReactionConditions.h"
#include "TInteractionCoordinates.h"
......@@ -45,8 +46,6 @@ class Analysis: public NPL::VAnalysis{
static NPL::VAnalysis* Construct();
private:
double OriginalBeamEnergy;
double OriginalEx;
double BeamEnergy;
double R;
double Xpista;
......@@ -58,20 +57,22 @@ class Analysis: public NPL::VAnalysis{
double XTarget;
double YTarget;
double ZTarget;
double OriginalElab;
double Elab;
double DeltaE;
double Eres;
double ThetaLab;
double PhiLab;
double ThetaCM;
double OptimumEx;
double Ex;
double PID;
double OriginalThetaLab;
int Telescope;
NPL::Reaction* Transfer;
double FF_A;
double FF_Z;
vector<double> init_FF_Z;
vector<double> init_FF_A;
TRandom3 Rand;
double ThetaNormalTarget;
double ThetaDetectorSurface;
......@@ -86,6 +87,7 @@ class Analysis: public NPL::VAnalysis{
TInteractionCoordinates* InteractionCoordinates;
TInitialConditions* InitialConditions;
TReactionConditions* ReactionConditions;
TFissionConditions* FissionConditions;
};
#endif
Projects/PISTA/GEF.source
+ 4
− 4
View file @ fba72618
......@@ -4,15 +4,15 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GEFReader
GEFversion= 2023.33
InputDataFile= /home/sofia/Physics/Codes/GEF_2024.3.3/out/240Pu_8MeV.lmd
GEFversion= 2023.13
InputDataFile= GEF/240Pu_14MeV.lmd
x0= 0
y0= 0
z0= 0 mm
TwoBodyReaction= 238U(12C,10Be)240Pu@1428
FissioningSystem= 240Pu
SigmaX= 0 mm
SigmaY= 0 mm
SigmaX= 1 mm
SigmaY= 1 mm
Particle= fragments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supported particle type: proton, neutron, gamma, fragments.
......
Projects/PISTA/PISTA.detector
+ 16
− 16
View file @ fba72618
......@@ -9,43 +9,43 @@ Target
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 315 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 270 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 225 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 180 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 135 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 90 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 45 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
R= 80 mm
THETA= 60 deg
R= 100 mm
THETA= 39 deg
PHI= 0 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......
Projects/PISTA/pista.detector
+ 2
− 2
View file @ fba72618
......@@ -14,12 +14,12 @@ Vamos
Theta= 20 deg
Vamos TMW1
Z = 0 mm
Z = 162.7 mm
Gas= iC4H10
Pressure= 0.01 bar
Vamos TMW2
Z = 109.3 mm
Z = 272 mm
Gas= iC4H10
Pressure= 0.01 bar
......
Projects/PISTA/pista_ana.detector 0 → 100644
+ 74
− 0
View file @ fba72618
%%%%%%%%%%%%%%%%%%%%
Target
THICKNESS= 0.44 micrometer
RADIUS= 5 mm
MATERIAL= C
ANGLE= 0 deg
X= 0 mm
Y= 0 mm
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= 36.71 96.1 50.7
POS_B= -35.92 95.8 51.3
POS_C= -19.35 56.5 93.4
POS_D= 20.8 56.6 93.1
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= -42.15 93.1 51.1
POS_B= -93.62 41.8 51.6
POS_C= -54.12 25.3 93.6
POS_D= -25.66 53.6 93.4
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= -95.98 35.6 51.3
POS_B= -95.63 -37.1 51.5
POS_C= -56.22 -20.7 93.7
POS_D= -56.39 19.5 93.6
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= -93.34 -43.5 50.9
POS_B= -41.88 -94.8 51.3
POS_C= -25.41 -55.2 93.2
POS_D= -53.84 -26.8 93.0
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= -35.66 -97.1 51.3
POS_B= 36.98 -96.9 51.1
POS_C= 20.9 -57.4 93.4
POS_D= -19.24 -57.6 93.5
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= 43.66 -94.4 50.9
POS_B= 95.19 -43.3 50.8
POS_C= 56.03 -26.6 93.1
POS_D= 27.54 -54.9 93.2
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= 96.87 -36.8 50.3
POS_B= 97.41 35.8 50.4
POS_C= 57.83 20.0 92.6
POS_D= 57.52 -20.1 92.6
%%%%%%%%%%%%%%%%%%%%
PISTA
VERSION= 1
POS_A= 94.39 42.0 50.8
POS_B= 42.94 93.3 50.7
POS_C= 26.47 54.1 93.0
POS_D= 54.9 25.7 93.1
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment