Merge branch 'NPTool.2.dev' of gitlab.in2p3.fr:np/nptool into NPTool.2.dev

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;

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;

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;

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;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

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;
 

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
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%