From c45feafb149e65cf8f22e2c67633e7b0302c2d8a Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@cea.fr>
Date: Fri, 13 Jan 2023 16:40:07 +0100
Subject: [PATCH] Updating on Sofia for simulation with GLAD field
---
NPLib/Detectors/Sofia/CMakeLists.txt | 6 +-
NPLib/Detectors/Sofia/TGladMagnetPhysics.cxx | 84 +++++
NPLib/Detectors/Sofia/TGladMagnetPhysics.h | 113 ++++++
NPSimulation/Detectors/Sofia/CMakeLists.txt | 4 +
NPSimulation/Detectors/Sofia/Glad.cc | 353 ++++++++++++++++++
NPSimulation/Detectors/Sofia/Glad.hh | 132 +++++++
.../Detectors/Sofia/GladFieldPropagation.cc | 175 +++++++++
.../Detectors/Sofia/GladFieldPropagation.hh | 78 ++++
NPSimulation/Detectors/Sofia/SofTofW.cc | 100 +----
Projects/Sofia/Analysis.cxx | 131 ++++++-
Projects/Sofia/Analysis.h | 28 +-
Projects/Sofia/PhysicsListOption.txt | 2 +-
Projects/Sofia/sofia.detector | 38 +-
13 files changed, 1111 insertions(+), 133 deletions(-)
create mode 100644 NPLib/Detectors/Sofia/TGladMagnetPhysics.cxx
create mode 100644 NPLib/Detectors/Sofia/TGladMagnetPhysics.h
create mode 100644 NPSimulation/Detectors/Sofia/Glad.cc
create mode 100644 NPSimulation/Detectors/Sofia/Glad.hh
create mode 100644 NPSimulation/Detectors/Sofia/GladFieldPropagation.cc
create mode 100644 NPSimulation/Detectors/Sofia/GladFieldPropagation.hh
diff --git a/NPLib/Detectors/Sofia/CMakeLists.txt b/NPLib/Detectors/Sofia/CMakeLists.txt
index b30076cb7..1ead29e3c 100644
--- a/NPLib/Detectors/Sofia/CMakeLists.txt
+++ b/NPLib/Detectors/Sofia/CMakeLists.txt
@@ -18,13 +18,15 @@ add_custom_command(OUTPUT TSofTwimPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/sc
add_custom_command(OUTPUT TSofBeamIDDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofBeamID.h TSofBeamIDDict.cxx TSofBeamID.rootmap libNPSofia.dylib DEPENDS TSofBeamID.h)
+add_custom_command(OUTPUT TGladMagnetPhysicsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TGladMagnetPhysics.h TGladMagnetPhysicsDict.cxx TGladMagnetPhysics.rootmap libNPSofia.dylib DEPENDS TGladMagnetPhysics.h)
+
add_custom_command(OUTPUT TSofFissionFragmentDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TSofFissionFragment.h TSofFissionFragmentDict.cxx TSofFissionFragment.rootmap libNPSofia.dylib DEPENDS TSofFissionFragment.h)
add_custom_command(OUTPUT GladFieldMapDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh GladFieldMap.h GladFieldMapDict.cxx GladFieldMap.rootmap libNPSofia.dylib DEPENDS GladFieldMap.h)
-add_library(NPSofia SHARED TSofSciData.cxx TSofSciDataDict.cxx TSofSciPhysics.cxx TSofSciPhysicsDict.cxx TSofMwpcData.cxx TSofMwpcDataDict.cxx TSofMwpcPhysics.cxx TSofMwpcPhysicsDict.cxx TSofAtData.cxx TSofAtDataDict.cxx TSofAtPhysics.cxx TSofAtPhysicsDict.cxx TSofTrimData.cxx TSofTrimDataDict.cxx TSofTrimPhysics.cxx TSofTrimPhysicsDict.cxx TSofTwimData.cxx TSofTwimDataDict.cxx TSofTwimPhysics.cxx TSofTwimPhysicsDict.cxx TSofTofWData.cxx TSofTofWDataDict.cxx TSofTofWPhysics.cxx TSofTofWPhysicsDict.cxx TSofBeamID.cxx TSofBeamIDDict.cxx TSofFissionFragment.cxx TSofFissionFragmentDict.cxx GladFieldMap.cxx GladFieldMapDict.cxx)
+add_library(NPSofia SHARED TSofSciData.cxx TSofSciDataDict.cxx TSofSciPhysics.cxx TSofSciPhysicsDict.cxx TSofMwpcData.cxx TSofMwpcDataDict.cxx TSofMwpcPhysics.cxx TSofMwpcPhysicsDict.cxx TSofAtData.cxx TSofAtDataDict.cxx TSofAtPhysics.cxx TSofAtPhysicsDict.cxx TSofTrimData.cxx TSofTrimDataDict.cxx TSofTrimPhysics.cxx TSofTrimPhysicsDict.cxx TSofTwimData.cxx TSofTwimDataDict.cxx TSofTwimPhysics.cxx TSofTwimPhysicsDict.cxx TSofTofWData.cxx TSofTofWDataDict.cxx TSofTofWPhysics.cxx TSofTofWPhysicsDict.cxx TSofBeamID.cxx TSofBeamIDDict.cxx TSofFissionFragment.cxx TSofFissionFragmentDict.cxx GladFieldMap.cxx GladFieldMapDict.cxx TGladMagnetPhysics.cxx TGladMagnetPhysicsDict.cxx)
target_link_libraries(NPSofia ${ROOT_LIBRARIES} NPCore NPPhysics)
-install(FILES TSofSciData.h TSofSciPhysics.h TSofMwpcData.h TSofMwpcPhysics.h TSofAtData.h TSofAtPhysics.h TSofTrimData.h TSofTrimPhysics.h TSofTwimData.h TSofTwimPhysics.h TSofTofWData.h TSofTofWPhysics.h TSofBeamID.h TSofFissionFragment.h GladFieldMap.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES TSofSciData.h TSofSciPhysics.h TSofMwpcData.h TSofMwpcPhysics.h TSofAtData.h TSofAtPhysics.h TSofTrimData.h TSofTrimPhysics.h TSofTwimData.h TSofTwimPhysics.h TSofTofWData.h TSofTofWPhysics.h TSofBeamID.h TSofFissionFragment.h GladFieldMap.h TGladMagnetPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/Sofia/TGladMagnetPhysics.cxx b/NPLib/Detectors/Sofia/TGladMagnetPhysics.cxx
new file mode 100644
index 000000000..c8baa2b9f
--- /dev/null
+++ b/NPLib/Detectors/Sofia/TGladMagnetPhysics.cxx
@@ -0,0 +1,84 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : May 2021 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold GladMagnet treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+#include "TGladMagnetPhysics.h"
+
+// STL
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <stdlib.h>
+#include <limits>
+using namespace std;
+
+// NPL
+#include "NPOptionManager.h"
+#include "NPDetectorFactory.h"
+#include "NPSystemOfUnits.h"
+// ROOT
+using namespace NPUNITS;
+///////////////////////////////////////////////////////////////////////////
+
+ClassImp(TGladMagnetPhysics)
+ ///////////////////////////////////////////////////////////////////////////
+ TGladMagnetPhysics::TGladMagnetPhysics(){
+ }
+
+///////////////////////////////////////////////////////////////////////////
+void TGladMagnetPhysics::Clear(){
+}
+///////////////////////////////////////////////////////////////////////////
+
+//// Innherited from VDetector Class ////
+
+///////////////////////////////////////////////////////////////////////////
+void TGladMagnetPhysics::ReadConfiguration(NPL::InputParser parser){
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("SAMURAIMagnet");
+ // nothing to do
+}
+
+///////////////////////////////////////////////////////////////////////////
+map< string , TH1*> TGladMagnetPhysics::GetSpectra() {
+ map< string , TH1*> empty;
+ return empty;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VDetector* TGladMagnetPhysics::Construct(){
+ return (NPL::VDetector*) new TGladMagnetPhysics();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy_gladMagnet{
+ public:
+ proxy_gladMagnet(){
+ NPL::DetectorFactory::getInstance()->AddToken("Glad","Glad");
+ NPL::DetectorFactory::getInstance()->AddDetector("Glad",TGladMagnetPhysics::Construct);
+ }
+ };
+
+ proxy_gladMagnet p_gladMagnet;
+}
+
diff --git a/NPLib/Detectors/Sofia/TGladMagnetPhysics.h b/NPLib/Detectors/Sofia/TGladMagnetPhysics.h
new file mode 100644
index 000000000..cc604c890
--- /dev/null
+++ b/NPLib/Detectors/Sofia/TGladMagnetPhysics.h
@@ -0,0 +1,113 @@
+#ifndef TGladMagnetPhysics_H
+#define TGladMagnetPhysics_H
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : October 2021 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class is a front to load the lib when the magnet is called, nothing *
+ * is done *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+// STL
+#include <vector>
+#include <map>
+#include <string>
+
+// NPL
+#include "NPVDetector.h"
+#include "NPInputParser.h"
+#include "NPXmlParser.h"
+
+#if __cplusplus > 199711L && NPMULTITHREADING
+#include "NPDCReconstructionMT.h"
+#else
+#include "NPDCReconstruction.h"
+#endif
+
+// ROOT
+#include "TVector3.h"
+
+// Forward declaration
+//class TGladMagnetSpectra;
+
+class TGladMagnetPhysics : public TObject, public NPL::VDetector{
+ public:
+ TGladMagnetPhysics();
+ ~TGladMagnetPhysics() {};
+
+ public:
+ void Clear();
+ void Clear(const Option_t*) {};
+
+ // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
+ void ReadConfiguration(NPL::InputParser) ;
+
+
+ // Add Parameter to the CalibrationManger
+ void AddParameterToCalibrationManager() {;};
+
+ // Activated associated Branches and link it to the private member DetectorData address
+ // In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
+ void InitializeRootInputRaw(){;} ;
+
+ // Activated associated Branches and link it to the private member DetectorPhysics address
+ // In this method mother Branches (Detector) AND daughter leaf (parameter) have to be activated
+ void InitializeRootInputPhysics(){;} ;
+
+ // Create associated branches and associated private member DetectorPhysics address
+ void InitializeRootOutput(){;} ;
+
+ // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter.
+ void BuildPhysicalEvent(){;} ;
+
+ // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...).
+ // This method aimed to be used for analysis performed during experiment, when speed is requiered.
+ // NB: This method can eventually be the same as BuildPhysicalEvent.
+ void BuildSimplePhysicalEvent() {;};
+
+ // Same as above but for online analysis
+ void BuildOnlinePhysicalEvent() {BuildPhysicalEvent();};
+
+ // Those two method all to clear the Event Physics or Data
+ void ClearEventPhysics() {Clear();}
+ void ClearEventData() {Clear();}
+
+ // Method related to the TSpectra classes, aimed at providing a framework for online applications
+ // Instantiate the Spectra class and the histogramm throught it
+ void InitSpectra(){;};
+ // Fill the spectra hold by the spectra class
+ void FillSpectra(){;};
+ // Used for Online mainly, perform check on the histo and for example change their color if issues are found
+ void CheckSpectra(){;};
+ // Used for Online only, clear all the spectra hold by the Spectra class
+ void ClearSpectra(){;};
+ // Write Spectra to file
+ void WriteSpectra(){;};
+
+ private: // Spectra Class
+ // TGladMagnetSpectra* m_Spectra; // !
+
+ public: // Spectra Getter
+ std::map< std::string , TH1*> GetSpectra();
+
+ public: // Static constructor to be passed to the Detector Factory
+ static NPL::VDetector* Construct();
+ ClassDef(TGladMagnetPhysics,1) // GladMagnetPhysics structure
+};
+
+#endif
diff --git a/NPSimulation/Detectors/Sofia/CMakeLists.txt b/NPSimulation/Detectors/Sofia/CMakeLists.txt
index b9667c66f..0aa49ff74 100644
--- a/NPSimulation/Detectors/Sofia/CMakeLists.txt
+++ b/NPSimulation/Detectors/Sofia/CMakeLists.txt
@@ -1,4 +1,8 @@
add_library(NPSSofTofW SHARED SofTofW.cc)
add_library(NPSSofTwim SHARED SofTwim.cc)
+add_library(NPSGladFieldPropagation SHARED GladFieldPropagation.cc)
+add_library(NPSGlad SHARED Glad.cc)
target_link_libraries(NPSSofTofW NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPSofia)
target_link_libraries(NPSSofTwim NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPSofia)
+target_link_libraries(NPSGladFieldPropagation NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPSofia)
+target_link_libraries(NPSGlad NPSGladFieldPropagation NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPSofia)
diff --git a/NPSimulation/Detectors/Sofia/Glad.cc b/NPSimulation/Detectors/Sofia/Glad.cc
new file mode 100644
index 000000000..eb194a33f
--- /dev/null
+++ b/NPSimulation/Detectors/Sofia/Glad.cc
@@ -0,0 +1,353 @@
+/******************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ ******************************************************************************/
+
+/******************************************************************************
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr*
+ * *
+ * Creation Date : November 2020 *
+ * Last update : *
+ *----------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe a simple Glad setup for simulation *
+ * *
+ *----------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ ******************************************************************************/
+
+// C++ headers
+#include <sstream>
+#include <cmath>
+#include <limits>
+//G4 Geometry object
+#include "G4Tubs.hh"
+#include "G4Box.hh"
+#include "G4SubtractionSolid.hh"
+#include "G4UserLimits.hh"
+
+//G4 sensitive
+#include "G4SDManager.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+//G4 various object
+#include "G4Material.hh"
+#include "G4Transform3D.hh"
+#include "G4PVPlacement.hh"
+#include "G4VisAttributes.hh"
+#include "G4Colour.hh"
+#include "G4FieldManager.hh"
+#include "G4UniformMagField.hh"
+#include "G4TransportationManager.hh"
+#include "G4FastSimulationManager.hh"
+#include "G4Region.hh"
+
+// NPTool header
+#include "Glad.hh"
+#include "CalorimeterScorers.hh"
+#include "InteractionScorers.hh"
+#include "RootOutput.h"
+#include "MaterialManager.hh"
+#include "NPSDetectorFactory.hh"
+#include "NPOptionManager.h"
+#include "NPSHitsMap.hh"
+// CLHEP header
+#include "CLHEP/Random/RandGauss.h"
+
+//CAD Mesh
+#include "CADMesh.hh"
+
+using namespace std;
+using namespace CLHEP;
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+namespace Glad_NS{
+ const double GLAD_height = 2000*mm; // y
+ const double GLAD_width = 4000*mm; // x
+ const double GLAD_Depth = 6000*mm; //z
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Glad Specific Method
+Glad::Glad(){
+ m_Magnet= NULL;
+ m_GLAD_STL= NULL;
+ m_PropagationRegion= NULL;
+ m_GladScorer= NULL;
+ m_GLAD_TiltAngle= -14*deg;
+ m_Current= 2185;
+ m_StepSize= 50*mm;
+
+ m_GLAD_X = 0;
+ m_GLAD_Y = 0;
+ m_GLAD_Z = -1113.5*mm;
+
+ // RGB Color + Transparency
+ m_VisSquare = new G4VisAttributes(G4Colour(0.53, 0.81, 0.98, 1));
+ m_VisGLAD = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.9));
+ m_VisField = new G4VisAttributes(G4Colour(0.1, 0.6, 0.9, 0.5));
+ m_VisKapton = new G4VisAttributes(G4Colour(1, 0.4, 0., 0.6));
+
+}
+
+Glad::~Glad(){
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Glad::AddMagnet(G4ThreeVector POS, double Tilt_Angle, string fieldmap){
+ // Convert the POS value to R theta Phi as Spherical coordinate is easier in G4
+ m_R = POS.mag();
+ m_Theta = POS.theta();
+ m_Phi = POS.phi();
+
+ m_GLAD_TiltAngle = Tilt_Angle;
+ m_FieldMapFile = fieldmap;
+
+ return;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Glad::AddMagnet(double R, double Theta, double Phi, double Tilt_Angle, string fieldmap){
+ m_R = R;
+ m_Theta = Theta;
+ m_Phi = Phi;
+
+ m_GLAD_TiltAngle= Tilt_Angle;
+ m_FieldMapFile = fieldmap;
+
+ return;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Glad::BuildGLADFromSTL()
+{
+ if(!m_GLAD_STL){
+ string basepath = getenv("NPTOOL");
+ string path = basepath + "/NPSimulation/Detectors/Sofia/stl/GLAD_only.stl";
+
+ auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) path.c_str());
+ mesh->SetScale(mm);
+
+ G4Material* GLAD_Material = MaterialManager::getInstance()->GetMaterialFromLibrary("Inox");
+
+ auto cad_solid = mesh->GetSolid();
+ m_GLAD_STL = new G4LogicalVolume(cad_solid,GLAD_Material,"GLAD_Magnet",0,0,0);
+
+ m_GLAD_STL->SetVisAttributes(m_VisGLAD);
+ }
+
+ return m_GLAD_STL;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Glad::BuildMagnet()
+{
+ if(!m_Magnet){
+ // Shape - Box
+ G4Box* box = new G4Box("glad_Box",Glad_NS::GLAD_width*0.5,Glad_NS::GLAD_height*0.5,Glad_NS::GLAD_Depth*0.5);
+
+ // Material
+ G4Material* vac = MaterialManager::getInstance()->GetMaterialFromLibrary("Vaccuum");
+
+ // Logical Volume
+ m_Magnet = new G4LogicalVolume(box,vac,"logic_GLAD_box",0,0,0);
+ m_Magnet->SetVisAttributes(m_VisField);
+ }
+ return m_Magnet;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Virtual Method of NPS::VDetector class
+
+// Read stream at Configfile to pick-up parameters of detector (Position,...)
+// Called in DetecorConstruction::ReadDetextorConfiguration Method
+void Glad::ReadConfiguration(NPL::InputParser parser){
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Glad");
+
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// GLAD Magnet found ////" << endl;
+
+ vector<string> cart = {"POS","TILT_ANGLE","FIELDMAP","CURRENT","STEPSIZE"};
+ vector<string> sphe = {"R","Theta","Phi","TILT_ANGLE","FIELDMAP","CURRENT","STEPSIZE"};
+
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(cart)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Glad information: " << endl;
+
+ G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
+ m_GLAD_TiltAngle = blocks[i]->GetDouble("TILT_ANGLE","deg");
+ string fieldmap = blocks[i]->GetString("FIELDMAP");
+ m_Current = blocks[i]->GetDouble("CURRENT","A");
+ m_StepSize = blocks[i]->GetDouble("STEPSIZE","mm");
+
+ AddMagnet(Pos,m_GLAD_TiltAngle,fieldmap);
+ }
+ else if(blocks[i]->HasTokenList(sphe)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Glad information: " << endl;
+
+ double R = blocks[i]->GetDouble("R","mm");
+ double Theta = blocks[i]->GetDouble("Theta","deg");
+ double Phi = blocks[i]->GetDouble("Phi","deg");
+ m_GLAD_TiltAngle = blocks[i]->GetDouble("TILT_ANGLE","deg");
+ string fieldmap = blocks[i]->GetString("FIELDMAP");
+ m_Current = blocks[i]->GetDouble("CURRENT","A");
+ m_StepSize = blocks[i]->GetDouble("STEPSIZE","mm");
+
+ AddMagnet(R,Theta,Phi,m_GLAD_TiltAngle,fieldmap);
+ }
+ else{
+ cout << "ERROR: check your GLAD input file formatting " << endl;
+ exit(1);
+ }
+ }
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Construct detector and inialise sensitive part.
+// Called After DetecorConstruction::AddDetector Method
+void Glad::ConstructDetector(G4LogicalVolume* world){
+
+ G4double wX = m_R * sin(m_Theta) * cos(m_Phi);
+ G4double wY = m_R * sin(m_Theta) * sin(m_Phi);
+ G4double wZ = m_R * cos(m_Theta);
+ G4ThreeVector Mag_pos = G4ThreeVector(wX,wY,wZ);
+
+ G4ThreeVector u (cos(m_GLAD_TiltAngle), 0, -sin(m_GLAD_TiltAngle));
+ G4ThreeVector v (0,1,0);
+ G4ThreeVector w (sin(m_GLAD_TiltAngle), 0, cos(m_GLAD_TiltAngle));
+ G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w);
+
+ new G4PVPlacement(Rot, Mag_pos,
+ BuildMagnet(), "Glad", world, false, 0);
+
+ SetPropagationRegion();
+
+ /*G4RotationMatrix* RotGLAD = new G4RotationMatrix();
+ RotGLAD->rotateX(90*deg);
+ //RotGLAD->rotateZ(-(90-m_GLAD_TiltAngle)*deg);
+ RotGLAD->rotateZ(-90*deg);
+ RotGLAD->rotateZ(m_GLAD_TiltAngle);
+ new G4PVPlacement(RotGLAD, GLAD_pos,
+ BuildGLADFromSTL(),
+ "GLAD",
+ world, false, 0);
+ */
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Glad::SetPropagationRegion()
+{
+ if(!m_PropagationRegion){
+ m_PropagationRegion= new G4Region("NPGladFieldPropagation");
+ //m_PropagationRegion -> AddRootLogicalVolume(BuildPropvol());
+ m_PropagationRegion -> AddRootLogicalVolume(BuildMagnet());
+ m_PropagationRegion->SetUserLimits(new G4UserLimits(m_StepSize));
+ }
+
+ G4FastSimulationManager* mng = m_PropagationRegion->GetFastSimulationManager();
+ //To make sure no other models are present in the region
+ unsigned int size = m_PropagationModel.size();
+ for(unsigned int i = 0 ; i < size ; i++)
+ mng->RemoveFastSimulationModel(m_PropagationModel[i]);
+ m_PropagationModel.clear();
+
+ G4VFastSimulationModel* fsm;
+ fsm = new NPS::GladFieldPropagation("GladFieldPropagation", m_PropagationRegion);
+ ((NPS::GladFieldPropagation*) fsm)->SetGladEntrance(m_GLAD_X, m_GLAD_Y, m_GLAD_Z);
+ ((NPS::GladFieldPropagation*) fsm)->SetCurrent(m_Current);
+ ((NPS::GladFieldPropagation*) fsm)->SetStepSize(m_StepSize);
+ ((NPS::GladFieldPropagation*) fsm)->SetFieldMap(m_FieldMapFile);
+ /*if(m_Method == NPS::RungeKutta){
+ double r_max = sqrt(
+ Glad_NS::GLAD_Width * Glad_NS::GLAD_Width /4. +
+ Glad_NS::GLAD_Depth * Glad_NS::GLAD_Depth /4. );//sqrt(x^2 + z^2)
+ ((NPS::GladFieldPropagation*) fsm)->SetRmax(r_max);
+ }*/
+ m_PropagationModel.push_back(fsm);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Add Detector branch to the EventTree.
+// Called After DetecorConstruction::AddDetector Method
+void Glad::InitializeRootOutput(){
+ // RootOutput *pAnalysis = RootOutput::getInstance();
+ // TTree *pTree = pAnalysis->GetTree();
+ // if(!pTree->FindBranch("IdealGladData")){
+ // pTree->Branch("IdealGladData", "TGladIdealData", &m_Event) ;
+ // }
+ // pTree->SetBranchAddress("IdealGladData", &m_Event) ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Read sensitive part and fill the Root tree.
+// Called at in the EventAction::EndOfEventAvtion
+void Glad::ReadSensitive(const G4Event* ){
+ //m_Event->Clear();
+
+ CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_GladScorer->GetPrimitive(0);
+
+ unsigned int size = Scorer->GetMult();
+ for(unsigned int i = 0 ; i < size ; i++){
+ double energy = Scorer->GetEnergy(i);
+ }
+ Scorer->clear();
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////
+void Glad::InitializeScorers() {
+ // This check is necessary in case the geometry is reloaded
+ bool already_exist = false;
+ m_GladScorer = CheckScorer("GladScorer",already_exist) ;
+
+ if(already_exist)
+ return ;
+
+ // Otherwise the scorer is initialised
+ vector<int> level;
+ level.push_back(0);
+ //G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ;
+ G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",level, 0) ;
+ //and register it to the multifunctionnal detector
+ //m_GladScorer->RegisterPrimitive(Interaction);
+ m_GladScorer->RegisterPrimitive(Calorimeter);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_GladScorer) ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPS::VDetector* Glad::Construct(){
+ return (NPS::VDetector*) new Glad();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern"C" {
+ class proxy_nps_glad{
+ public:
+ proxy_nps_glad(){
+ NPS::DetectorFactory::getInstance()->AddToken("Glad","Glad");
+ NPS::DetectorFactory::getInstance()->AddDetector("Glad",Glad::Construct);
+ }
+ };
+
+ proxy_nps_glad p_nps_glad;
+}
diff --git a/NPSimulation/Detectors/Sofia/Glad.hh b/NPSimulation/Detectors/Sofia/Glad.hh
new file mode 100644
index 000000000..1c3f215f9
--- /dev/null
+++ b/NPSimulation/Detectors/Sofia/Glad.hh
@@ -0,0 +1,132 @@
+#ifndef Glad_h
+#define Glad_h 1
+/*****************************************************************************
+ * Copyright (C) 2009-2020 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr *
+ * *
+ * Creation Date : November 2020 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe Glad simulation *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// C++ header
+#include <string>
+#include <vector>
+using namespace std;
+
+// G4 headers
+#include "G4ThreeVector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4LogicalVolume.hh"
+#include "G4MultiFunctionalDetector.hh"
+#include "G4AssemblyVolume.hh"
+#include "G4VFastSimulationModel.hh"
+#include "G4VSolid.hh"
+
+// NPTool header
+#include "NPSVDetector.hh"
+#include "NPInputParser.h"
+
+#include "GladFieldPropagation.hh"
+
+class Glad : public NPS::VDetector{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+ public:
+ Glad() ;
+ virtual ~Glad() ;
+
+ ////////////////////////////////////////////////////
+ /////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+ public:
+ // Cartesian
+ void AddMagnet(G4ThreeVector POS, double Tilt_Angle, string fieldmap);
+ // Spherical
+ void AddMagnet(double R,double Theta,double Phi, double Tilt_Angle, string fieldmap);
+
+ G4LogicalVolume* BuildGLADFromSTL();
+ G4LogicalVolume* BuildMagnet();
+
+ private:
+ G4LogicalVolume* m_Magnet;
+ G4LogicalVolume* m_GLAD_STL;
+
+ ////////////////////////////////////////////////////
+ ////// Inherite from NPS::VDetector class /////////
+ ////////////////////////////////////////////////////
+ public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(NPL::InputParser) ;
+
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world) ;
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput() ;
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event) ;
+
+ public: // Scorer
+ void InitializeScorers() ;
+
+ // Set region were magnetic field is active
+ void SetPropagationRegion();
+
+ // Associated Scorer
+ G4MultiFunctionalDetector* m_GladScorer ;
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+ private:
+ //TGladData* m_Event ;
+ G4Region* m_PropagationRegion;
+ vector<G4VFastSimulationModel*> m_PropagationModel;
+
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+ private: // Geometry
+ // Detector Coordinate
+ double m_R;
+ double m_Theta;
+ double m_Phi;
+
+ // GLAD //
+ double m_GLAD_X;
+ double m_GLAD_Y;
+ double m_GLAD_Z;
+ double m_GLAD_TiltAngle;
+ double m_Current;
+ int m_StepSize;
+ string m_FieldMapFile;
+
+ // Visualisation Attribute
+ G4VisAttributes* m_VisSquare;
+ G4VisAttributes* m_VisGLAD;
+ G4VisAttributes* m_VisField;
+ G4VisAttributes* m_VisKapton;
+
+ // Needed for dynamic loading of the library
+ public:
+ static NPS::VDetector* Construct();
+};
+#endif
diff --git a/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc b/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc
new file mode 100644
index 000000000..55d366fd6
--- /dev/null
+++ b/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc
@@ -0,0 +1,175 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Elia Pilotto, Omar Nasr *
+ * contact address: pilottoelia@gmail.com, omar.nasr@etu.unicaen.fr *
+ * *
+ * Creation Date : September 2021 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+//C++ libraries
+#include <iostream>
+#include <string>
+#include <iomanip>
+#include <cmath>
+#include <Randomize.hh>
+#include <fstream>
+//G4 libraries
+#include "G4Electron.hh"
+#include "G4Gamma.hh"
+#include "G4EmCalculator.hh"
+#include "G4VPhysicalVolume.hh"
+#include "G4IonTable.hh"
+#include "G4UserLimits.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4PhysicalConstants.hh"
+//nptool libraries
+#include "NPFunction.h"
+#include "NPInputParser.h"
+#include "NPOptionManager.h"
+#include "NPSFunction.hh"
+//other
+#include "GladFieldPropagation.hh"
+#include "RootOutput.h"
+#include "TLorentzVector.h"
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::GladFieldPropagation::GladFieldPropagation(G4String modelName, G4Region* envelope)
+ : G4VFastSimulationModel(modelName, envelope) {
+
+ m_Map = NULL;
+ m_Initialized = false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::GladFieldPropagation::GladFieldPropagation(G4String modelName)
+ : G4VFastSimulationModel(modelName) {}
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::GladFieldPropagation::~GladFieldPropagation() {}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::GladFieldPropagation::AttachReactionConditions() {
+ // Reasssigned the branch address
+ /*
+ if (RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+ RootOutput::getInstance()->GetTree()->SetBranchAddress(
+ "ReactionConditions", &m_ReactionConditions);*/
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool NPS::GladFieldPropagation::IsApplicable(const G4ParticleDefinition& particleType) {
+ if (particleType.GetPDGCharge() == 0) return false;
+ return true;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool NPS::GladFieldPropagation::ModelTrigger(const G4FastTrack& fastTrack) {
+ return true;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::GladFieldPropagation::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
+ if(!m_Initialized){
+ m_Map = new GladFieldMap();
+
+ //Needed for the RungeKutta method
+ m_Map->SetCurrent(m_Current);
+ m_Map->SetLimit( 10.*m / m_StepSize );
+ m_Map->SetPropagationTimeInterval(0.8);
+ m_Map->SetGladEntrance(m_GladEntrance.X(), m_GladEntrance.Y(), m_GladEntrance.Z());
+
+ m_Map->LoadMap(m_FieldMap);
+ m_Initialized = true;
+ }
+
+ RungeKuttaPropagation(fastTrack, fastStep);
+
+ return;
+}
+
+/////////////////////////////////////////////////////////////////////////
+void NPS::GladFieldPropagation::RungeKuttaPropagation (const G4FastTrack& fastTrack, G4FastStep& fastStep){
+
+ static bool inside = false;//true if previous step is inside the volume
+ static vector<TVector3> trajectory;
+ static int count;//keeps track of the step reached inside trajectory
+
+ // Get the track info
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+
+ static G4VSolid* solid;//Stores the current solid
+ solid = PrimaryTrack->GetTouchable()->GetSolid();
+
+ G4ThreeVector localDir = fastTrack.GetPrimaryTrackLocalDirection();
+ G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
+ G4ThreeVector localMomentum = fastTrack.GetPrimaryTrackLocalMomentum();
+
+ double speed = PrimaryTrack->GetVelocity();
+ double charge = PrimaryTrack->GetParticleDefinition()->GetPDGCharge() / coulomb;
+ double ConF_p = 1 / (joule * c_light / (m/s)) ; // MeV/c to kg*m/s (SI units)
+
+ //Initially inside is false -> calculate trajectory
+ if (!inside){
+ count = 2;//skip first two positions as they are the same as the current position
+ double Brho = localMomentum.mag() * ConF_p / charge;
+ TVector3 pos (localPosition.x(), localPosition.y(), localPosition.z());
+ TVector3 dir (localDir.x(), localDir.y(), localDir.z());
+
+ trajectory.clear();
+ trajectory = m_Map->Propagate(Brho, pos, dir,true);
+
+ inside = true;
+ }
+
+ G4ThreeVector newPosition (trajectory[count].x(), trajectory[count].y(), trajectory[count].z());
+ //benchmark
+ //G4ThreeVector newDir = (newPosition - localPosition).unit();
+ //G4ThreeVector newMomentum = newDir * localMomentum.mag();
+
+ //Check if newPosition is not inside
+ if (solid->Inside(newPosition) != kInside){
+ inside = false;
+ G4ThreeVector toOut = solid->DistanceToOut(localPosition, localDir) * localDir;
+ newPosition = localPosition + toOut;
+ //benchmark
+ //newDir = (newPosition - localPosition).unit();
+ //newMomentum = newDir * localMomentum.mag();
+
+ //benchmark
+ //if(single_particle) PrintData(m_StepSize, newPosition, newMomentum, outRKsp);
+ //else PrintData(m_StepSize, newPosition, newMomentum, outRK);
+ //counter++;
+ }
+
+ //benchmark
+ G4ThreeVector newDir = (newPosition - localPosition).unit();
+ G4ThreeVector newMomentum = newDir * localMomentum.mag();
+
+
+ double time = PrimaryTrack->GetGlobalTime()+(newPosition - localPosition).mag()/speed;
+
+ fastStep.ProposePrimaryTrackFinalPosition( newPosition );
+ fastStep.SetPrimaryTrackFinalMomentum ( newMomentum );
+ fastStep.ProposePrimaryTrackFinalTime( time );
+
+ count++;
+ return;
+}
+
+
diff --git a/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh b/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh
new file mode 100644
index 000000000..7e833a2d8
--- /dev/null
+++ b/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh
@@ -0,0 +1,78 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre Morfouace *
+ * contact address: pierre.morfouace@cea.fr *
+ * *
+ * Creation Date : January 2023 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: Based On SamuraiFieldPropagation *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#ifndef GladFieldPropagation_h
+#define GladFieldPropagation_h
+
+#include "G4VFastSimulationModel.hh"
+#include "G4Abla.hh"
+#include "G4AblaInterface.hh"
+#include "G4Fragment.hh"
+#include "G4VPhysicalVolume.hh"
+#include "G4Region.hh"
+
+#include "GladFieldMap.h"
+
+#include <string>
+#include <vector>
+#include <fstream>
+#include <iomanip>
+
+namespace NPS{
+
+ class GladFieldPropagation : public G4VFastSimulationModel{
+ public:
+ GladFieldPropagation (G4String, G4Region*);
+ GladFieldPropagation (G4String);
+ ~GladFieldPropagation ();
+
+ public:
+ G4bool IsApplicable(const G4ParticleDefinition&);
+ G4bool ModelTrigger(const G4FastTrack &);
+ void DoIt(const G4FastTrack&, G4FastStep&);
+
+ private:
+ bool m_Initialized;
+ double m_StepSize;
+ double m_Rmax;
+ double m_Current;
+ TVector3 m_GladEntrance;
+ string m_FieldMap;
+ GladFieldMap* m_Map;
+
+ void RungeKuttaPropagation (const G4FastTrack& fastTrack, G4FastStep& fastStep);
+
+ public:
+ void AttachReactionConditions();
+
+ void SetStepSize(double step){m_StepSize=step;};
+ void SetFieldMap(string fieldmap){m_FieldMap=fieldmap;};
+ void SetRmax(double r_max){m_Rmax=r_max;};
+ void SetCurrent(double val){m_Current=val;}
+ void SetGladEntrance(double x, double y, double z){m_GladEntrance=TVector3(x,y,z);}
+ };
+}
+
+
+#endif
diff --git a/NPSimulation/Detectors/Sofia/SofTofW.cc b/NPSimulation/Detectors/Sofia/SofTofW.cc
index 093a4eaa2..078e9c02f 100644
--- a/NPSimulation/Detectors/Sofia/SofTofW.cc
+++ b/NPSimulation/Detectors/Sofia/SofTofW.cc
@@ -72,11 +72,6 @@ namespace SofTofW_NS{
const double tof_plastic_width = 32*mm;
const double tof_plastic_thickness = 0.5*mm;
const string Material = "BC400";
-
- const double GLAD_height = 1.5*m;
- const double GLAD_width = 10*m;
- const double GLAD_Leff = 2*m;
-
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -86,26 +81,15 @@ SofTofW::SofTofW(){
m_Event = new TSofTofWData() ;
m_TofScorer = 0;
m_PlasticTof = 0;
- m_GLAD= 0;
- m_GLAD_STL= 0;
m_VacuumPipe= 0;
m_TofWall = 0;
- m_Build_GLAD= 0;
- m_Build_MagneticField= 0;
- m_Build_VacuumPipe= 0;
m_VacuumPipeX= 0;
m_VacuumPipeY= 0;
m_VacuumPipeZ= 0;
- m_GLAD_MagField = 0;
- m_GLAD_DistanceFromTarget= 0;
- m_GLAD_TiltAngle= 7.;
- m_GLAD_Yoffset= 0.;
// RGB Color + Transparency
m_VisSquare = new G4VisAttributes(G4Colour(0.53, 0.81, 0.98, 1));
- m_VisGLAD = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.9));
- m_VisField = new G4VisAttributes(G4Colour(0.1, 0.6, 0.9, 0.9));
m_VisKapton = new G4VisAttributes(G4Colour(1, 0.4, 0., 0.6));
}
@@ -190,57 +174,6 @@ G4AssemblyVolume* SofTofW::BuildVacuumPipe(){
return m_VacuumPipe;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* SofTofW::BuildGLADFromSTL()
-{
- if(!m_GLAD_STL){
- string basepath = getenv("NPTOOL");
- string path = basepath + "/NPSimulation/Detectors/Sofia/stl/GLAD_only.stl";
-
- auto mesh = CADMesh::TessellatedMesh::FromSTL((char*) path.c_str());
- mesh->SetScale(mm);
-
- G4Material* GLAD_Material = MaterialManager::getInstance()->GetMaterialFromLibrary("Inox");
-
- auto cad_solid = mesh->GetSolid();
- m_GLAD_STL = new G4LogicalVolume(cad_solid,GLAD_Material,"GLAD_Magnet",0,0,0);
-
- m_GLAD_STL->SetVisAttributes(m_VisGLAD);
- }
-
- return m_GLAD_STL;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4AssemblyVolume* SofTofW::BuildMagneticField()
-{
- if(!m_GLAD){
- m_GLAD = new G4AssemblyVolume;
-
- // *** GLAD field *** //
- G4Box* box = new G4Box("glad_Box",SofTofW_NS::GLAD_width*0.5,SofTofW_NS::GLAD_height*0.5,SofTofW_NS::GLAD_Leff*0.5);
- G4Material* vac = MaterialManager::getInstance()->GetMaterialFromLibrary("Vaccuum");
- G4LogicalVolume* vol_field = new G4LogicalVolume(box,vac,"logic_GLAD_field",0,0,0);
- vol_field->SetVisAttributes(m_VisField);
- //vol_field->SetVisAttributes(G4VisAttributes::Invisible);
-
- G4UniformMagField* magField = new G4UniformMagField(G4ThreeVector(0,m_GLAD_MagField,0));
- //G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
- G4FieldManager* fieldMgr = new G4FieldManager(magField);
-
- //fieldMgr->SetDetectorField(magField);
-
- fieldMgr->CreateChordFinder(magField);
-
- vol_field->SetFieldManager(fieldMgr,true);
-
- G4ThreeVector Pos_field = G4ThreeVector(0,0,1*m);
- G4RotationMatrix* Rot_field = new G4RotationMatrix();
- Rot_field->rotateY(-m_GLAD_TiltAngle);
- m_GLAD->AddPlacedVolume(vol_field,Pos_field,Rot_field);
-
- }
- return m_GLAD;
-}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -253,8 +186,8 @@ void SofTofW::ReadConfiguration(NPL::InputParser parser){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
- vector<string> cart = {"POS","Build_GLAD","Build_VacuumPipe"};
- vector<string> sphe = {"R","Theta","Phi","Build_GLAD","Build_VacuumPipe"};
+ vector<string> cart = {"POS","Build_VacuumPipe"};
+ vector<string> sphe = {"R","Theta","Phi","Build_VacuumPipe"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
@@ -262,8 +195,6 @@ void SofTofW::ReadConfiguration(NPL::InputParser parser){
cout << endl << "//// SofTofW " << i+1 << endl;
G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
- m_Build_GLAD = blocks[i]->GetInt("Build_GLAD");
- m_Build_MagneticField = blocks[i]->GetInt("Build_MagneticField");
m_Build_VacuumPipe = blocks[i]->GetInt("Build_VacuumPipe");
AddDetector(Pos);
}
@@ -273,13 +204,7 @@ void SofTofW::ReadConfiguration(NPL::InputParser parser){
double R = blocks[i]->GetDouble("R","mm");
double Theta = blocks[i]->GetDouble("Theta","deg");
double Phi = blocks[i]->GetDouble("Phi","deg");
- m_GLAD_Yoffset = blocks[i]->GetDouble("Yoffset","mm");
- m_Build_GLAD = blocks[i]->GetInt("Build_GLAD");
- m_GLAD_TiltAngle = blocks[i]->GetDouble("GLAD_TiltAngle","deg");
m_Build_VacuumPipe = blocks[i]->GetInt("Build_VacuumPipe");
- m_Build_MagneticField = blocks[i]->GetInt("Build_MagneticField");
- m_GLAD_MagField = blocks[i]->GetDouble("GLAD_MagField","T");
- m_GLAD_DistanceFromTarget = blocks[i]->GetDouble("GLAD_DistanceFromTarget", "m");
m_VacuumPipeX = blocks[i]->GetDouble("VacuumPipeX","m");
m_VacuumPipeY = blocks[i]->GetDouble("VacuumPipeY","m");
m_VacuumPipeZ = blocks[i]->GetDouble("VacuumPipeZ","m");
@@ -320,30 +245,11 @@ void SofTofW::ConstructDetector(G4LogicalVolume* world){
u = u.unit();
G4RotationMatrix* Rot = new G4RotationMatrix(u,v,w);
- //G4ThreeVector Z_translation = G4ThreeVector(0,0,4*m);
- G4ThreeVector Z_translation = G4ThreeVector(0,0,m_GLAD_DistanceFromTarget);
- Det_pos += Z_translation;
BuildTOFDetector()->MakeImprint(world,Det_pos,Rot);
}
- G4ThreeVector GLAD_pos = G4ThreeVector(0,m_GLAD_Yoffset,m_GLAD_DistanceFromTarget);
- if(m_Build_MagneticField==1){
- BuildMagneticField()->MakeImprint(world,GLAD_pos,0);
- }
- if(m_Build_GLAD==1){
- G4RotationMatrix* RotGLAD = new G4RotationMatrix();
- RotGLAD->rotateX(90*deg);
- //RotGLAD->rotateZ(-(90-m_GLAD_TiltAngle)*deg);
- RotGLAD->rotateZ(-90*deg);
- RotGLAD->rotateZ(m_GLAD_TiltAngle);
- new G4PVPlacement(RotGLAD, GLAD_pos,
- BuildGLADFromSTL(),
- "GLAD",
- world, false, 0);
-
- }
-
+
if(m_Build_VacuumPipe==1){
G4ThreeVector Tube_Pos = G4ThreeVector(m_VacuumPipeX,m_VacuumPipeY,m_VacuumPipeZ);
BuildVacuumPipe()->MakeImprint(world,Tube_Pos,0);
diff --git a/Projects/Sofia/Analysis.cxx b/Projects/Sofia/Analysis.cxx
index bd7d89b37..f2cb83078 100644
--- a/Projects/Sofia/Analysis.cxx
+++ b/Projects/Sofia/Analysis.cxx
@@ -21,9 +21,11 @@
#include<iostream>
using namespace std;
-#include"Analysis.h"
-#include"NPAnalysisFactory.h"
-#include"NPDetectorManager.h"
+#include "Analysis.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "RootInput.h"
+
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -33,11 +35,118 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
- Sofia= (TSofiaPhysics*) m_DetectorManager->GetDetector("Sofia");
+ //Sofia= (TSofiaPhysics*) m_DetectorManager->GetDetector("Sofia");
+
+ InitialConditions = new TInitialConditions();
+ InteractionCoordinates = new TInteractionCoordinates();
+ FissionConditions = new TFissionConditions();
+
+ InitInputBranch();
+ InitOutputBranch();
+
+ m_GladField = new GladFieldMap();
+ m_GladField->SetCurrent(2185);
+ m_GladField->SetGladEntrance(0,0,-1.1135*m);
+ m_GladField->SetGladTiltAngle(-14*deg);
+ m_GladField->LoadMap("GladFieldMap_50mm.dat");
+ m_GladField->SetBin(50);
+ m_GladField->SetTimeStep(0.8);
+ m_GladField->SetCentralTheta(-20*deg);
+
+ double Z_MW3 = 4440 * cos(20*deg);
+ double X_MW3 = Z_MW3 * tan(20*deg);
+ m_GladField->Set_MWPC3_Position(X_MW3,0,Z_MW3);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
+
+ ReInit();
+
+ int mult = InteractionCoordinates->GetDetectedMultiplicity();
+ if(mult==2){
+ int A_CN = FissionConditions->GetA_CN();
+ int Z_CN = FissionConditions->GetZ_CN();
+ double Elab_CN = FissionConditions->GetELab_CN();
+ double z_start_beam = InitialConditions->GetIncidentPositionZ();
+ double z_target = -4093;
+ TVector3 v1 = TVector3(0,0,z_target-z_start_beam);
+ TVector3 v2 = TVector3(0,0,z_start_beam);
+ NPL::Particle* beam = new NPL::Particle(Z_CN,A_CN);
+ beam->SetKineticEnergy(Elab_CN);
+ double beta = beam->GetBeta();
+ double BeamTimeOffset = v1.Mag()*mm/(beta*NPUNITS::c_light);
+ for(int i=0; i<mult; i++){
+ double Time = ran.Gaus(InteractionCoordinates->GetTime(i) - BeamTimeOffset,0.02);
+ double XE = ran.Gaus(InteractionCoordinates->GetDetectedPositionX(i),0.1);
+ double YE = 0;//ran.Gaus(InteractionCoordinates->GetDetectedPositionY(i),0.1);
+ double ZE = ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(i),0.1);
+ TVector3 vE = TVector3(XE,YE,ZE);
+ TVector3 vA = TVector3(0,0,z_target);
+ int A = FissionConditions->GetFragmentA(i);
+ int Z = FissionConditions->GetFragmentZ(i);
+ double Theta = FissionConditions->GetFragmentTheta(i);
+ double Phi = FissionConditions->GetFragmentPhi(i);
+ double Brho = FissionConditions->GetFragmentBrho(i);
+
+ TVector3 dir = TVector3(sin(Theta*deg)*cos(Phi*deg), sin(Theta*deg)*sin(Phi*deg), cos(Theta*deg));
+ m_TOF.push_back(Time);
+ m_A.push_back(A);
+ m_Brho.push_back(Brho);
+
+ double Brho_calc = m_GladField->FindBrho(vA,dir,vE);
+ double PathLength = m_GladField->GetFlightPath(vA,Brho_calc,vA,dir);
+ double velocity = PathLength/Time;
+ double beta = velocity * m/ns / NPUNITS::c_light;
+ double gamma = 1. / sqrt(1 - beta*beta);
+ double AoQ = Brho_calc / (3.107 * beta * gamma);
+ double A_calc = AoQ*Z;
+
+ m_Brho_calc.push_back(Brho_calc);
+ m_FlightPath.push_back(PathLength);
+ m_A.push_back(A_calc);
+ }
+ }
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fIC_*",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions);
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("InteractionCoordinates",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fDetecetd_*",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&InteractionCoordinates);
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("FissionConditions",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fFC_*",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("FissionConditions",&FissionConditions);
+
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch(){
+ RootOutput::getInstance()->GetTree()->Branch("TOF",&m_TOF);
+ RootOutput::getInstance()->GetTree()->Branch("Brho",&m_Brho);
+ RootOutput::getInstance()->GetTree()->Branch("A",&m_A);
+ RootOutput::getInstance()->GetTree()->Branch("Brho_calc",&m_Brho_calc);
+ RootOutput::getInstance()->GetTree()->Branch("A_calc",&m_A_calc);
+ RootOutput::getInstance()->GetTree()->Branch("FlightPath",&m_FlightPath);
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInit(){
+ m_TOF.clear();
+ m_Brho.clear();
+ m_A.clear();
+
+ m_Brho_calc.clear();
+ m_A_calc.clear();
+ m_FlightPath.clear();
}
////////////////////////////////////////////////////////////////////////////////
@@ -56,13 +165,13 @@ NPL::VAnalysis* Analysis::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy{
- public:
- proxy(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
-};
+ class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
-proxy p;
+ proxy p;
}
diff --git a/Projects/Sofia/Analysis.h b/Projects/Sofia/Analysis.h
index a4f50f8fa..636524f5d 100644
--- a/Projects/Sofia/Analysis.h
+++ b/Projects/Sofia/Analysis.h
@@ -21,8 +21,13 @@
* *
*****************************************************************************/
-#include"NPVAnalysis.h"
-#include"TSofiaPhysics.h"
+#include "NPVAnalysis.h"
+#include "TInitialConditions.h"
+#include "TInteractionCoordinates.h"
+#include "TFissionConditions.h"
+#include "GladFieldMap.h"
+#include "TRandom3.h"
+
class Analysis: public NPL::VAnalysis{
public:
Analysis();
@@ -30,13 +35,28 @@ class Analysis: public NPL::VAnalysis{
public:
void Init();
+ void InitInputBranch();
+ void InitOutputBranch();
void TreatEvent();
+ void ReInit();
void End();
- static NPL::VAnalysis* Construct();
+ static NPL::VAnalysis* Construct();
+
+ private:
+ vector<double> m_TOF;
+ vector<double> m_Brho;
+ vector<int> m_A;
+ vector<double> m_Brho_calc;
+ vector<int> m_A_calc;
+ vector<double> m_FlightPath;
private:
- TSofiaPhysics* Sofia;
+ TInitialConditions* InitialConditions;
+ TInteractionCoordinates* InteractionCoordinates;
+ TFissionConditions* FissionConditions;
+ TRandom3 ran;
+ GladFieldMap* m_GladField;
};
#endif
diff --git a/Projects/Sofia/PhysicsListOption.txt b/Projects/Sofia/PhysicsListOption.txt
index ee6f9bf67..8edb4ba48 100644
--- a/Projects/Sofia/PhysicsListOption.txt
+++ b/Projects/Sofia/PhysicsListOption.txt
@@ -1,5 +1,5 @@
EmPhysicsList Option4
-DefaultCutOff 10000
+DefaultCutOff 100000
IonBinaryCascadePhysics 0
NPIonInelasticPhysics 0
EmExtraPhysics 0
diff --git a/Projects/Sofia/sofia.detector b/Projects/Sofia/sofia.detector
index 3cba4f8d2..6d08511d5 100644
--- a/Projects/Sofia/sofia.detector
+++ b/Projects/Sofia/sofia.detector
@@ -6,27 +6,29 @@ Target
ANGLE= 0 deg
X= 0 mm
Y= 0 mm
- Z= 0 m
+ Z= -4093 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SofTofW
- R= 4.5 m
- THETA= -18 deg
- PHI= 0 deg
- Yoffset= 0.0 mm
- Build_GLAD= 1
- GLAD_TiltAngle= 14 deg
- GLAD_DistanceFromTarget= 3.3 m
- Build_MagneticField= 1
- GLAD_MagField= 1.75 T
- Build_VacuumPipe= 1
- VacuumPipeX= 0 cm
- VacuumPipeY= 0 cm
- VacuumPipeZ= 0.9 m
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-SofTwim
- R= 2.00 m
+%SofTwim
+ R= -1850 mm
THETA= 0
PHI= 0
TwimGas= MixTwinMusic
Pressure= 1 bar
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+SofTofW
+ R= 4440 mm
+ THETA= -20 deg
+ PHI= 0 deg
+ Build_VacuumPipe= 0
+ VacuumPipeX= 0 cm
+ VacuumPipeY= 0 cm
+ VacuumPipeZ= -3250 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Glad
+ POS= 0 0 386.5 mm
+ TILT_ANGLE= 0 deg
+ FIELDMAP= GladFieldMap_50mm.dat
+ CURRENT= 2185
+ STEPSIZE= 50 mm
--
GitLab