diff --git a/NPLib/Detectors/PISTA/TPISTAPhysics.cxx b/NPLib/Detectors/PISTA/TPISTAPhysics.cxx
index 7266d3d536f5b3d8325076a69385e997e9362e56..c0893cc90e0c4512e5ce865c93031b880fa948ee 100644
--- a/NPLib/Detectors/PISTA/TPISTAPhysics.cxx
+++ b/NPLib/Detectors/PISTA/TPISTAPhysics.cxx
@@ -67,6 +67,15 @@ void TPISTAPhysics::AddDetector(TVector3){
void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
m_NumberOfDetectors++;
+ //double Height = 118; // mm
+ double Height = 61.8; // mm
+ //double Base = 95; // mm
+ double Base = 78.1; // mm
+ double NumberOfStrips = 128;
+ double StripPitchHeight = Height / NumberOfStrips; // mm
+ double StripPitchBase = Base / NumberOfStrips; // mm
+
+
// Vector U on detector face (parallel to Y strips) Y strips are along X axis
TVector3 U;
// Vector V on detector face (parallel to X strips)
@@ -78,7 +87,7 @@ void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
C = TVector3(R*sin(Theta)*cos(Phi),
R*sin(Theta)*sin(Phi),
- R*cos(Theta));
+ Height*0.5+R*cos(Theta));
TVector3 P = TVector3(cos(Theta)*cos(Phi),
cos(Theta)*sin(Phi),
@@ -91,12 +100,6 @@ void TPISTAPhysics::AddDetector(double R, double Theta, double Phi){
U = U.Unit();
V = V.Unit();
- double Height = 118; // mm
- double Base = 95; // mm
- double NumberOfStrips = 128;
- double StripPitchHeight = Height / NumberOfStrips; // mm
- double StripPitchBase = Base / NumberOfStrips; // mm
-
vector<double> lineX;
vector<double> lineY;
vector<double> lineZ;
diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc
index e66bc6168a3e9aa0269523fcf3bf3732af3b4bd2..cb9e2933807bc7907dd4bd378f419ab27cbc6be3 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -278,6 +278,35 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
return material;
}
+ else if(Name == "235U"){
+ if(!density)
+ density = 19.1 * g / cm3;
+ G4Element* U235 = new G4Element("U235","U235",1);
+
+ G4Isotope* isotope = new G4Isotope("235U",92,235);
+ U235->AddIsotope(isotope,1);
+
+ G4Material* material = new G4Material("NPS_" + Name, density, 1, kStateSolid, 293.15*kelvin);
+ material->AddElement(U235,1);
+ m_Material[Name] = material;
+ return material;
+ }
+
+ else if(Name == "238U"){
+ if(!density)
+ density = 19.1 * g / cm3;
+ G4Element* U238 = new G4Element("U238","U238",1);
+
+ G4Isotope* isotope = new G4Isotope("238U",92,238);
+ U238->AddIsotope(isotope,1);
+
+ G4Material* material = new G4Material("NPS_" + Name, density, 1, kStateSolid, 293.15*kelvin);
+ material->AddElement(U238,1);
+ m_Material[Name] = material;
+ return material;
+ }
+
+
else if (Name == "Gd") {
if (!density)
density = 7.90 * g / cm3;
diff --git a/NPSimulation/Detectors/PISTA/PISTA.cc b/NPSimulation/Detectors/PISTA/PISTA.cc
index 0a02fb280cd0201f8d7fa0f46c2ee74e747111bf..420612ff28788053c7732d3c59f08ae9dadc0567 100644
--- a/NPSimulation/Detectors/PISTA/PISTA.cc
+++ b/NPSimulation/Detectors/PISTA/PISTA.cc
@@ -63,9 +63,12 @@ namespace PISTA_NS{
const double ResoEnergy = 0.015*MeV ;
// Trapezoid dimension
- const double TrapezoidBaseLarge = 95*mm;
- const double TrapezoidBaseSmall = 45*mm;
- const double TrapezoidHeight = 118*mm;
+ //const double TrapezoidBaseLarge = 95*mm;
+ const double TrapezoidBaseLarge = 78.1*mm;
+ //const double TrapezoidBaseSmall = 45*mm;
+ const double TrapezoidBaseSmall = 43.3*mm;
+ //const double TrapezoidHeight = 118*mm;
+ const double TrapezoidHeight = 61.8*mm;
const double TrapezoidLength = 1*cm;
const double FirstStageThickness = 100*um;
const double SecondStageThickness = 1*mm;
@@ -226,7 +229,7 @@ void PISTA::ConstructDetector(G4LogicalVolume* world){
G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ;
G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ;
- G4double wZ = m_R[i] * cos(m_Theta[i] ) ;
+ G4double wZ = TrapezoidHeight*0.5 + m_R[i] * cos(m_Theta[i] ) ;
G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ;
// So the face of the detector is at R instead of the middle
Det_pos+=Det_pos.unit()*PISTA_NS::TrapezoidLength*0.5;
diff --git a/Projects/PISTA/Analysis.cxx b/Projects/PISTA/Analysis.cxx
index ed8955b67a5196df8a1ff48574375837dddae3a2..e01382896d6ca5ed8be1a5efe048d68754402ff9 100644
--- a/Projects/PISTA/Analysis.cxx
+++ b/Projects/PISTA/Analysis.cxx
@@ -39,6 +39,7 @@ void Analysis::Init(){
PISTA= (TPISTAPhysics*) m_DetectorManager->GetDetector("PISTA");
InitialConditions = new TInitialConditions();
ReactionConditions = new TReactionConditions();
+ InteractionCoordinates = new TInteractionCoordinates();
InitOutputBranch();
InitInputBranch();
Rand = TRandom3();
@@ -60,14 +61,25 @@ void Analysis::TreatEvent(){
OriginalElab = ReactionConditions->GetKineticEnergy(0);
OriginalBeamEnergy = ReactionConditions->GetBeamEnergy();
+ int mult = InteractionCoordinates->GetDetectedMultiplicity();
+ if(mult>0){
+ for(int i=0; i<mult; i++){
+ double Xpista = InteractionCoordinates->GetDetectedPositionX(i);
+ double Ypista = InteractionCoordinates->GetDetectedPositionY(i);
+ double Zpista = InteractionCoordinates->GetDetectedPositionZ(i);
+ R = sqrt(Xpista*Xpista + Ypista*Ypista + Zpista*Zpista);
+ }
+ }
XTarget = InitialConditions->GetIncidentPositionX();
YTarget = InitialConditions->GetIncidentPositionY();
ZTarget = InitialConditions->GetIncidentPositionZ();
TVector3 BeamDirection = InitialConditions->GetBeamDirection();
TVector3 BeamPosition(XTarget,YTarget,ZTarget);
+ //TVector3 PositionOnTarget(0,0,0);
+ TVector3 PositionOnTarget(Rand.Gaus(XTarget, 0.6/2.35), Rand.Gaus(YTarget, 0.6/2.35), 0);
- BeamEnergy = InitialConditions->GetIncidentInitialKineticEnergy();
+ BeamEnergy = 1428.;//InitialConditions->GetIncidentInitialKineticEnergy();
BeamEnergy = U238C.Slow(BeamEnergy,TargetThickness*0.5,0);
Transfer->SetBeamEnergy(BeamEnergy);
@@ -75,8 +87,7 @@ void Analysis::TreatEvent(){
for(unsigned int i = 0; i<PISTA->EventMultiplicity; i++){
if(PISTA->E.size()>0){
double Energy = PISTA->DE[i] + PISTA->E[i];
-
- TVector3 HitDirection = PISTA->GetPositionOfInteraction(i);
+ TVector3 HitDirection = PISTA->GetPositionOfInteraction(i)-PositionOnTarget;
//ThetaLab = HitDirection.Angle(BeamDirection);
ThetaLab = HitDirection.Angle(TVector3(0,0,1));
@@ -107,6 +118,7 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/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");
}
////////////////////////////////////////////////////////////////////////////////
@@ -117,6 +129,9 @@ void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("fRC_*",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("InteractionCoordinates",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fDetected_*",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&InteractionCoordinates);
}
////////////////////////////////////////////////////////////////////////////////
@@ -134,6 +149,7 @@ void Analysis::ReInitValue(){
YTarget = -1000;
ZTarget = -1000;
OriginalThetaLab = -1000;
+ R = -1000;
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/PISTA/Analysis.h b/Projects/PISTA/Analysis.h
index 3ec3cf61893310cc6e5a272963e05a2bee31b792..f5b182ca138522491ca781ec07f754d1b87583af 100644
--- a/Projects/PISTA/Analysis.h
+++ b/Projects/PISTA/Analysis.h
@@ -25,6 +25,7 @@
#include "TPISTAPhysics.h"
#include "TInitialConditions.h"
#include "TReactionConditions.h"
+#include "TInteractionCoordinates.h"
#include "NPEnergyLoss.h"
#include "NPReaction.h"
#include "TRandom3.h"
@@ -46,6 +47,7 @@ class Analysis: public NPL::VAnalysis{
private:
double OriginalBeamEnergy;
double BeamEnergy;
+ double R;
double XTarget;
double YTarget;
double ZTarget;
@@ -68,6 +70,7 @@ class Analysis: public NPL::VAnalysis{
private:
TPISTAPhysics* PISTA;
+ TInteractionCoordinates* InteractionCoordinates;
TInitialConditions* InitialConditions;
TReactionConditions* ReactionConditions;
diff --git a/Projects/PISTA/PISTA.detector b/Projects/PISTA/PISTA.detector
index bc33042c2a7a66e9f543302264bf1f56b9298a04..d815893669198c4db742f88ce2ee87578f35b6b1 100644
--- a/Projects/PISTA/PISTA.detector
+++ b/Projects/PISTA/PISTA.detector
@@ -9,42 +9,42 @@ Target
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 315 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 270 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 225 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 180 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 135 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 90 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 45 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PISTA
- R= 100 mm
+ R= 85 mm
THETA= 60 deg
PHI= 0 deg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%