Addapted Mike's geant program for ZDC beam test 2018, this includes adding...

Addapted Mike's geant program for ZDC beam test 2018, this includes adding magnetic field, zdc class, rpd class, new particle generator, rpdhit/sd class. Also, began to make code more readable.

CMakeLists.txt

@@ -74,4 +74,7 @@ install(TARGETS test_JZCaPA RUNTIME DESTINATION bin)
 install(FILES doxygen/JZCaPA_doxy.cnf DESTINATION doxygen)
 install(FILES Utils/ConfigFile2018.xml DESTINATION Utils)
 install(FILES Utils/Alignment_2018.xml DESTINATION Utils)
+install(FILES Utils/Survey_2018.xml DESTINATION Utils)
+install(FILES Utils/PurgMag3D.TABLE DESTINATION Utils)
+
 

MonteCarlo/CMakeLists.txt

@@ -29,6 +29,10 @@ include(${Geant4_USE_FILE})
 find_package (ROOT REQUIRED)
 include_directories (${ROOT_INCLUDE_DIR})
 
+#Xerces-C support
+find_package (Xerces REQUIRED)
+include_directories (${XERCESC_INCLUDE})
+
 # Add ROOTs header paths
 
 #message(${CMAKE_MODULE_PATH} )
@@ -63,7 +67,8 @@ if(useROOT)
 endif(useROOT)
 
 
-target_link_libraries(zdc ${Geant4_LIBRARIES} ${ROOT_LIBRARIES})
+target_link_libraries(zdc ${Geant4_LIBRARIES} ${ROOT_LIBRARIES} ${XERCESC_LIBRARY})
+
 
 #----------------------------------------------------------------------------
 # Copy all scripts to the build directory, i.e. the directory in which we

MonteCarlo/config/config.cfg

 #################################
+########## Survey Setup #########
+
+########        if using this code to simulate the 2018 ZDC beamtest set "ZDC_SETUP" = 1       #########
+########              set nModules = 3, mod1Type = 5, mod2Type = 5, mod3Type = 6               #########
+
+ZDC_SETUP						1
+
+RunNumber:						79
+
+
+
+
 ###### Module Attributes ########
 
 # currently 5 mods are maximum allowed
-# mod types: 1 (EM-A), 2 (H1), 3 (standardMod), 4 (custom)
+# mod types: 1 (EM-A), 2 (H1), 3 (standardMod), 4 (custom), 5 (zdc)
 # if mod type 1-3 chosen, geometry from ATHENA simulation used
 # if mod type 4 (custom) chosen, dimensions and materials must be set
 # if you want just a block of absorber (eg for a testbeam) use type 4 and set core diameter to 0
@@ -10,13 +22,88 @@
 # Width of radiator gap assumed to be 0.5 mm greater than fiber diameter
 # Units in mm
 
-nModules:                      1
+nModules:                      3
 
-mod1Type:                      4
-mod2Type:                      4
-mod3Type:                      4
+mod1Type:                      5
+mod2Type:                      5
+mod3Type:                      6
 mod4Type:                      1
 
+
+
+
+#######################################ZDC_MODULE_DO_NOT_ALTER#######################################
+#####################################################################################################
+mod5CasingThickness:           7.94
+mod5NStripsPerGap:             52
+# W and Pb are two currently enabled
+mod5AbsorberMat:               W
+mod5AbsorberThickness:         10.2
+mod5AbsorberHeight:            182.8
+# AbsorberWidth is only used if NStripsPerGap = 0 -- ie absorber only mode
+mod5AbsorberWidth:             90.
+# CoreDiameter should never be 0 -- even in absorber only mode
+mod5CoreDiameter:              1.5
+# Number of absorber layers assumed 1 fewer than # of radiators, unless 0 Radiators is entered, in which case the entire "module" is assumed to be an absorber block (like in test beam case when a Pb block is inserted)
+mod5NRadiators:                12
+# make sure you account for core + 2*cladding thickness + any empty buffer space
+mod5RadiatorGapLength:         2.
+mod5CoreIndexRefraction:       1.46
+mod5Cladding:                  NO
+# only if cladding enabled
+mod5CladdingThickness:         0.
+mod5CladdingIndexRefraction:   1.43
+#####################################################################################################
+#######################################ZDC_MODULE_DO_NOT_ALTER#######################################
+
+
+
+####################################### RPD_MODULE #######################################
+##########################################################################################
+# X spacing between tiles in mm
+modGapX = 1.58
+modGapY = 0.5
+#1.58 mm =1/16in
+
+# dimensions of individual tile in mm
+tileXsize = 20
+tileYsize = 20
+tileZsize = 10
+
+mod6CasingThickness:           7.94
+mod6NStripsPerGap:             16
+# W and Pb are two currently enabled
+mod6AbsorberMat:               W
+mod6AbsorberThickness:         10.2
+mod6AbsorberHeight:            182.8
+# AbsorberWidth is only used if NStripsPerGap = 0 -- ie absorber only mode
+mod6AbsorberWidth:             90.
+# CoreDiameter should never be 0 -- even in absorber only mode
+mod6CoreDiameter:              1.5
+# Number of absorber layers assumed 1 fewer than # of radiators, unless 0 Radiators is entered, in which case the entire "module" is assumed to be an absorber block (like in test beam case when a Pb block is inserted)
+mod6NRadiators:                1
+# make sure you account for core + 2*cladding thickness + any empty buffer space
+mod6RadiatorGapLength:         2.
+mod6CoreIndexRefraction:       1.46
+mod6Cladding:                  NO
+# only if cladding enabled
+mod6CladdingThickness:         0.
+mod6CladdingIndexRefraction:   1.43
+
+####################################### RPD_MODULE #######################################
+##########################################################################################
+
+
+
+
+
+
+
+
+
+
+
+
 #Custom categorizes (use same syntax for modN -- w/ 5 being the max # of mods supported)
 mod1CasingThickness:           7.94
 mod1NStripsPerGap:             52

MonteCarlo/include/DetectorConstruction.hh

@@ -37,7 +37,10 @@
 #include "ModType2.hh"
 #include "ModType3.hh"
 #include "ModTypeCustom.hh"
-
+#include "ModTypeZDC.hh"
+#include "ModTypeRPD.hh"
+#include "G4Cache.hh"
+#include "G4MagneticField.hh"
 #include <vector>
 
 class G4Box;
@@ -47,6 +50,10 @@ class G4LogicalVolume;
 class G4Material;
 class SharedData;
 
+class G4UniformMagField;
+class PurgMagTabulatedField3D;
+
+
 /// Detector construction class to define materials and geometry.
 
 class DetectorConstruction : public G4VUserDetectorConstruction
@@ -69,6 +76,11 @@ protected:
   G4Box*               m_solidWorld;
   G4LogicalVolume*     m_logicWorld;
   G4VPhysicalVolume*   m_physWorld;
+  G4LogicalVolume*     logic_leadTarget;
+  G4LogicalVolume*     logic_leadBlock;
+
+private:
+	G4Cache<G4MagneticField*> fField;  //pointer to the thread-local fields
   
 };
 

MonteCarlo/include/PrimaryGeneratorAction.hh

@@ -36,14 +36,10 @@
 #include "G4IonTable.hh"
 #include "globals.hh"
 
-class G4ParticleGun;
+
+class G4GeneralParticleSource;
 class G4Event;
-class G4Box;
 
-/// The primary generator action class with particle gun.
-///
-/// The default kinematic is a 6 MeV gamma, randomly distribued 
-/// in front of the phantom across 80% of the (X,Y) phantom size.
 
 class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
 {
@@ -53,13 +49,10 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
 
     // method from the base class
     virtual void GeneratePrimaries(G4Event*);         
-  
-    // method to access particle gun
-    const G4ParticleGun* GetParticleGun() const { return m_particleGun; }
+   
   
   private:
-    G4ParticleGun*  m_particleGun; // pointer a to G4 gun class
-    G4Box* m_world;
+	G4GeneralParticleSource*   fParticleGun;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

MonteCarlo/include/RunAction.hh

@@ -77,6 +77,7 @@ class RunAction : public G4UserRunAction
         fVelocity_v.clear();
         fBeta_v.clear();}
   inline SharedData* GetSharedData() {return fSharedData;}
+  
   inline  void SetRadNo(G4int rNo) {fRadNb_v.push_back(rNo);}
   inline  void SetRodNo(G4int rNo) {fRodNb_v.push_back(rNo);}  
   inline  void SetEdep(G4double edep) {fEdep_v.push_back(edep);}

MonteCarlo/include/SharedData.hh

@@ -16,10 +16,56 @@
 #include <TH1.h>
 #include <TTree.h>
 #include <TLorentzVector.h>
+#include "XMLSettingsReader.hh"
 
 #include <iostream>
 #include <string>
 
+class Survey {
+
+ public :
+  
+	/** Type of detector - ZDC or RPD **/
+    std::string detector;
+    /** x_pos for this survey */
+    double x_pos;
+	/** y_pos for this survey */
+    double y_pos;
+	/** z_pos for this survey */
+    double z_pos;
+	/** cos_x for this survey */
+    double cos_x;
+	/** cos_y for this survey */
+    double cos_y;
+	/** cos_z for this survey */
+    double cos_z;
+	
+};
+
+class Alignment {
+
+ public:
+
+    /** X position of the Desy Table **/
+    double x_table;
+    /** Y position of the Desy Table **/
+    double y_table;
+    /** First detector met by the beam **/
+    std::string upstream_Det;
+    /** Second detector met by the beam **/
+    std::string mid_Det;
+    /** Third detector met by the beam **/
+    std::string downstream_Det;
+    /** GOLIATH magnet status **/
+    bool magnet_On;
+    /** Target in **/
+    bool target_In;
+    /** Lead absorber in **/
+    bool lead_In;
+
+};
+
+
 class SharedData{
   
 public:
@@ -32,12 +78,11 @@ public:
   SharedData& operator= ( const SharedData& ) = delete ;
 
   // Only want MyRunManager to access these 
-public:
+
   void   Initialize();
   void   EndOfEvent();
   void   Finalize(); 
    
-public:
   template<class T> 
   void   AddOutputToTree    ( const std::string&, T*);
 
@@ -48,6 +93,15 @@ public:
   bool   DoPrint          ();
   inline TTree* GetTree          () {return m_tree;}
   inline int GetEventNo       () {return m_eventCounter;}
+  
+  Survey* GetSurvey(std::string name);
+  Alignment* GetAlignment();
+  
+  void LoadConfigurationFile(int runNum, std::string _inFile = std::getenv("JZCaPA") + std::string("/Utils/Survey_2018.xml"));
+  void LoadAlignmentFile(int runNum, std::string _inFile = std::getenv("JZCaPA") + std::string("/Utils/Alignment_2018.xml"));
+   
+   
+   
 private:
   TEnv*   m_config;
   TFile*  m_fout;
@@ -58,6 +112,16 @@ private:
 
   std::string  m_outputFileName;
   std::string  m_configFileName;
+  
+   //XML parser
+  XMLSettingsReader *m_XMLparser;
+
+  //Alignment information for the given run
+	Alignment* 	m_alignment;
+  //Survey information for the given run
+	Survey* m_survey;
+	std::vector < Survey* > surveyEntries;
+	
 };
 
 /** @brief Function to add an object to the tree

MonteCarlo/run1.mac

-# Macro file for example B1
+# Macro file for JZCaPA beam test 2018
 # 
 # To be run preferably in batch, without graphics:
-# % exampleB1 run2.mac
+# 
 #
-#/run/numberOfWorkers 4
 
-#
+# Initialize kernel
+/run/initialize
+
+
 /control/verbose 2
 /run/verbose 2
-# 
-# proton 210 MeV to the direction (0.,0.,-1.)
-# 1000 events
+/tracking/verbose 0
+/run/particle/verbose 1
+
+
+#/gps/particle ion
+#/gps/ion 82 208 0 0.0
+#/gps/energy 31.2 TeV
+
+# using alphas is useful if you would like to visualize the tracks (Pb ions require much more computing time)
+/gps/particle ion
+/gps/ion 2 4 0 0.0
+/gps/energy 600 GeV
+
+#/gps/particle proton
+#/gps/particle neutron
+
+# the beam energy is in gaussian profile
+#/gps/ene/type Gauss
+#/gps/ene/mono 16.8 GeV
+#/gps/ene/sigma 0.3 MeV
+
+
+# General particle source  
+/gps/pos/type Beam
+/gps/pos/shape Circle
+/gps/pos/centre 0. 0. -15999. mm
+/gps/pos/radius 1. mm
+/gps/pos/sigma_r 2. mm
+
 #
-#/testhadr/ionPhysics FTF
-/run/initialize
-/gun/particle proton
-#/gun/particle ion
-#/gun/ion Z A charge excitationEnergy
-#/gun/ion 82 208 0 0.0
-/gun/energy 3.5 TeV
 #
-/run/beamOn 1
+# the incident surface is in the x-y plane
+/gps/pos/rot1 1 0 0
+/gps/pos/rot2 0 1 0
+#
+#
+# the beam is travelling along the z-axis with any angular 
+#dispersion (angular despersion set to 0.057)
+/gps/ang/rot1 0 1 0
+/gps/ang/rot2 1 0 0 
+/gps/ang/type beam1d 
+/gps/ang/sigma_r 0.057 deg 
+# 1.0milliradians   0.057 deg 
+#
+
+#if using these commands in interactive mode to visualize, use this to allow more events to be displayed
+/vis/ogl/set/displayListLimit 1000000
+
+##############################################################
+##############################################################
+# number of events
+/run/beamOn 2
+##############################################################
+##############################################################
\ No newline at end of file

MonteCarlo/run2.mac

@@ -13,11 +13,11 @@
 # 1000 events
 #
 #/testhadr/ionPhysics FTF
-/run/initialize
-/gun/particle proton
+#/run/initialize
+#/gun/particle proton
 #/gun/particle ion
 #/gun/ion Z A charge excitationEnergy
 #/gun/ion 82 208 0 0.0
-/gun/energy 3.5 TeV
+#/gun/energy 3.5 TeV
 #
-/run/beamOn 1
+#/run/beamOn 1

MonteCarlo/src/CherenkovSD.cc

@@ -203,15 +203,17 @@ G4bool CherenkovSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+
 void CherenkovSD::EndOfEvent(G4HCofThisEvent* HCE)
 {
-  /*
+	
+  
   //  G4String name = HCE->GetHC()->GetName();
   //if(HCID<0) return;
    std::cout << " HCID end of event " << HCID << std::endl;
   CherenkovHitsCollection* HC = 0;
-  G4int nCollections =  HCE->GetNumberOfCollections();
-  G4int hitsCollID = 0;
+	__attribute__((unused)) G4int nCollections =  HCE->GetNumberOfCollections();
+  __attribute__((unused)) G4int hitsCollID = 0;
 
   
   if(HCE) {
@@ -260,8 +262,9 @@ void CherenkovSD::EndOfEvent(G4HCofThisEvent* HCE)
       // }
   }
   //  std::cout << " vec size " << v_fRadNb.size() << std::endl;
-  */
+  
 }
 
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 

MonteCarlo/src/DetectorConstruction.cc

@@ -52,6 +52,28 @@
 #include "G4VisAttributes.hh"
 #include "G4Colour.hh"
 
+#include "G4FieldManager.hh"
+#include "G4TransportationManager.hh"
+#include "G4EqMagElectricField.hh"
+#include "PurgMagTabulatedField3D.hh"
+#include "G4ChordFinder.hh"
+#include "G4Mag_UsualEqRhs.hh"
+#include "G4PVParameterised.hh"
+#include "G4ThreeVector.hh"
+#include "G4PVReplica.hh"
+#include "G4UniformMagField.hh"
+#include "G4ExplicitEuler.hh"
+#include "G4ImplicitEuler.hh"
+#include "G4SimpleRunge.hh"
+#include "G4SimpleHeum.hh"
+#include "G4ClassicalRK4.hh"
+#include "G4HelixExplicitEuler.hh"
+#include "G4HelixImplicitEuler.hh"
+#include "G4HelixSimpleRunge.hh"
+#include "G4CashKarpRKF45.hh"
+#include "G4RKG3_Stepper.hh"
+
+
 #include <iostream>
 #include <stdio.h>
 
@@ -98,20 +120,101 @@ void DetectorConstruction :: DefineBorderProperties()
 
 G4VPhysicalVolume* DetectorConstruction::ConstructDetector()
 {
-  // Get Config
-  TEnv* config = m_sd->GetConfig();
+	// Get Config
+	TEnv* config = m_sd->GetConfig();
+	int ZDC_SETUP = config->GetValue("ZDC_SETUP", 0);
+	int runNum = config->GetValue( "RunNumber", -1);
+	
+	m_sd->LoadConfigurationFile(runNum);
+	m_sd->LoadAlignmentFile(runNum);
+
+	// Create variables to be used in beamtest 2018 simulation
+	
+	G4ThreeVector zdc1Pos,zdc2Pos, rpdPos;
+	bool 		mag_on = false, targ_in = false, lead_in = false, 
+				ZDC1 = false, ZDC2 = false, RPD = false;
+	G4double 	mag_zOffset=0, density=0, fractionmass=0, leadblockZ=0, worldZoffset=16000*mm,
+				zdc1X=0,zdc1Y=0,zdc1Z=0,zdc2X=0,zdc2Y=0,zdc2Z=0,rpdX=0,rpdY=0,rpdZ=0,
+				tableX_shift=0, tableY_shift=0;
+	std::string detector[3];
+	G4int   	ncomponents;  
+	Survey		*srvy_zdc1 = m_sd->GetSurvey("ZDC1"), 
+				*srvy_zdc2 = m_sd->GetSurvey("ZDC2"), 
+				*srvy_rpd = m_sd->GetSurvey("RPD");;
+	Alignment	*align_run 	= m_sd->GetAlignment();
+
+	//################################ SURVEY/ALIGNMENT_SETUP
+	if(ZDC_SETUP == 1){
+	
+	//table(-2250,500) -> rpd/beam(0,0)	where 100=0.1cm in table coordinates
+	tableX_shift = (-2250.0 - (align_run->x_table)  )/100*mm ;
+	tableY_shift = (500.0   - (align_run->y_table)  )/100*mm ;
+		
+	rpdX  = (srvy_rpd ->x_pos)   *1000.0			 					*mm;	
+	zdc1X = (( (srvy_zdc1->x_pos)*1000.0 ) - rpdX + tableX_shift    )	*mm;
+	zdc2X = (( (srvy_zdc2->x_pos)*1000.0 ) - rpdX + tableX_shift	)	*mm;
+	rpdX  =  tableX_shift;
+		
+	rpdY  = (srvy_rpd ->y_pos)   *1000.0								*mm;	
+	zdc1Y = (( (srvy_zdc1->y_pos)*1000.0 ) - 320 - rpdY + tableY_shift)	*mm;
+	zdc2Y = (( (srvy_zdc2->y_pos)*1000.0 ) - 320 - rpdY + tableY_shift)	*mm;
+	rpdY  =  tableY_shift;
+	
+	rpdZ  = (( (srvy_rpd ->z_pos)*1000.0 ) -(worldZoffset) )*mm;
+	zdc1Z = (( (srvy_zdc1->z_pos)*1000.0 ) -(worldZoffset) )*mm;
+	zdc2Z = (( (srvy_zdc2->z_pos)*1000.0 ) -(worldZoffset) )*mm;
+	
+	 
+	//-320mm is offset to get from zdc mid to active area mid
+	zdc1Pos = G4ThreeVector( zdc1X, zdc1Y, zdc1Z); 
+	zdc2Pos = G4ThreeVector( zdc2X, zdc2Y, zdc2Z);
+	rpdPos 	= G4ThreeVector( rpdX , rpdY , rpdZ );
   
-  // Get nist material manager
+	
+	detector[0]=align_run->upstream_Det;
+	detector[1]=align_run->mid_Det;
+	detector[2]=align_run->downstream_Det;
+	
+	for(int i=0; i<3; i++){
+		if(detector[i]=="ZDC1") {
+			ZDC1 = true;
+		}
+		if(detector[i]=="ZDC2") {
+			ZDC2 = true;
+		}
+		if(detector[i]=="RPD") {
+			RPD = true;
+		}
+	}
+	
+	// Assign lead block position in mm
+	leadblockZ = ( zdc1Z- 250)*mm; //approximate position
+	
+	
+	targ_in =  	align_run->target_In;
+	mag_on 	= 	align_run->magnet_On;
+	lead_in = 	align_run->lead_In;
+	}
+	//################################ SURVEY/ALIGNMENT_END
+	
+	
   G4NistManager* nist = G4NistManager::Instance();
+		
+		G4Element* Pb = nist->FindOrBuildElement(82);
+		G4Material* Lead = new G4Material("Lead", density= 11.35*g/cm3, ncomponents=1);
+		Lead->AddElement(Pb, fractionmass=1.0);  
+	
+  // Get nist material manager
+  // G4NistManager* nist = G4NistManager::Instance();
 
   //----------------------------------------------     
   // Set Some Values
   //----------------------------------------------
   G4double modSizeX[5]; G4double modSizeY[5]; G4double modSizeZ[5]; G4double modCasingThickness[5]; G4double modWidth[5]; G4String modAbsorberMat[5]; G4double modAbsorberThickness[5]; G4double modAbsorberHeight[5]; G4double modAbsorberWidth[5]; G4double modCoreDiameter[5]; G4double modCladdingThickness[5]; G4int modNRadiators[5]; G4int modNAbsorbers[5]; G4int modType[5]; bool cladding[5]; G4double modRadiatorGapLength[5]; G4double stripPitch[5]; G4int modNStripsPerGap[5];
   
-  G4int    nModules            = config->GetValue( "nModules",4);
-  modType[0]            = config->GetValue( "mod1Type",4);
-  modType[1]            = config->GetValue( "mod2Type",2);
+  G4int    nModules            = config->GetValue( "nModules",2);
+  modType[0]            = config->GetValue( "mod1Type",5);
+  modType[1]            = config->GetValue( "mod2Type",5);
   modType[2]            = config->GetValue( "mod3Type",3);
   modType[3]            = config->GetValue( "mod4Type",3);
   modType[4]            = config->GetValue( "mod5Type",3);    
@@ -158,6 +261,46 @@ G4VPhysicalVolume* DetectorConstruction::ConstructDetector()
       modSizeZ[i] = 2*modCasingThickness[i]+modNRadiators[i]*modRadiatorGapLength[i] + modNAbsorbers[i]*modAbsorberThickness[i];
       
     }
+	if(ZDC_SETUP==1){
+		if (modType[i] == 5) {
+		char variable[256];
+		std::string modCladding;
+		sprintf(variable,"mod%dCasingThickness",5);
+		modCasingThickness[i] = config->GetValue(variable,.605);
+		sprintf(variable,"mod%dAbsorberThickness",5);
+		modAbsorberThickness[i] = config->GetValue(variable,10.);
+		sprintf(variable,"mod%dAbsorberHeight",5);
+		modAbsorberHeight[i] = config->GetValue(variable,180.);
+		sprintf(variable,"mod%dAbsorberWidth",5);
+		modAbsorberWidth[i] = config->GetValue(variable,100.);
+		sprintf(variable,"mod%dRadiatorGapLength",5);
+		modRadiatorGapLength[i] = config->GetValue(variable,2.);
+		sprintf(variable,"mod%dCoreDiameter",5);
+		modCoreDiameter[i] = config->GetValue(variable,1.5);
+		sprintf(variable,"mod%dCladding",5);
+		modCladding = config->GetValue(variable,"true");
+		sprintf(variable,"mod%dCladdingThickness",5);
+		modCladdingThickness[i] = config->GetValue(variable,0.605);
+		sprintf(variable,"mod%dNRadiators",5);
+		modNRadiators[i] = config->GetValue(variable,12);
+		sprintf(variable,"mod%dNStripsPerGap",5);
+		modNStripsPerGap[i] = config->GetValue(variable,52);
+		modNAbsorbers[i] = modNRadiators[i] - 1;
+		stripPitch[i] = modCoreDiameter[i] + 2*modCladdingThickness[i];