Addition of a rock envelope for background studies

CMakeLists.txt

@@ -90,14 +90,15 @@ endforeach()
 
 #----------------------------------------------------------------------------
 # generate dictionaries for the custom classes in the ROOT tree
-ROOT_GENERATE_DICTIONARY(DictOutput ${PROJECT_SOURCE_DIR}/include/FPFNeutrino.hh 
-                                    ${PROJECT_SOURCE_DIR}/include/FPFParticle.hh 
-				    LINKDEF ${PROJECT_SOURCE_DIR}/include/LinkDef.h)
+ROOT_GENERATE_DICTIONARY(FPFClassesDict ${PROJECT_SOURCE_DIR}/include/FPFNeutrino.hh 
+                                        ${PROJECT_SOURCE_DIR}/include/FPFParticle.hh
+                                        LINKDEF ${PROJECT_SOURCE_DIR}/include/LinkDef.h
+                                        MODULE FPFClasses)
 
 # create a shared library that includes the dictionary
 add_library(FPFClasses SHARED ${PROJECT_SOURCE_DIR}/src/FPFNeutrino.cc
                               ${PROJECT_SOURCE_DIR}/src/FPFParticle.cc
-                              DictOutput.cxx)
+                              FPFClassesDict.cxx)
 
 # include directories for headers
 target_include_directories(FPFClasses PUBLIC ${PROJECT_SOURCE_DIR}/include)
@@ -121,4 +122,4 @@ target_link_libraries(FPFSim
 #
 install(TARGETS FPFSim DESTINATION ${PROJECT_SOURCE_DIR}/bin)
 install(TARGETS FPFClasses DESTINATION ${PROJECT_SOURCE_DIR}/lib)
-install(FILES ${PROJECT_SOURCE_DIR}/build/libDictOutput_rdict.pcm DESTINATION ${PROJECT_SOURCE_DIR}/lib)
+install(FILES ${PROJECT_SOURCE_DIR}/build/libFPFClasses_rdict.pcm DESTINATION ${PROJECT_SOURCE_DIR}/lib)

README.md

@@ -120,7 +120,8 @@ Older versions of FORESEE output events in the HepMC2 format. To run over HepMC2
 |/det/addFLArE          | option for adding the FLArE detector, run before `/run/initialize`                   |`true`|
 |/det/addFORMOSA        | option for adding the FORMOSA detector, run before `/run/initialize`                 |`true`|
 |/det/addFASERnu2       | option for adding the FASERnu2 detector, run before `/run/initialize`                |`true`|
-|/det/faser/addFASER2   | option for adding the FASER2 detector, run before `/run/initialize`                  |`true`|
+|/det/addFASER2         | option for adding the FASER2 detector, run before `/run/initialize`                  |`true`|
+|/det/enableRock        | enable the rock envelope around the hall, run before `/run/initialize`               |`false`|
 |/det/flare/addFLArEPos | position of the FLArE detector, run before `/run/initialize`                         |`0 0 4300 mm`|
 |/det/flare/material    | option for detector material, choose `LAr` or `LKr`, run before `/run/initialize`    |`LAr`|
 |/det/flare/module      | option for tpc module option, choose `single` or `3x7`, run before `/run/initialize` |`single`|

include/DetectorConstruction.hh

@@ -36,6 +36,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction {
     void AddFASERnu2(G4bool i) { m_addFASERnu2 = i; }
     void AddFASER2(G4bool i) { m_addFASER2 = i; }
     void UseBabyMIND(G4bool i) { m_useBabyMIND = i; }
+    void EnableRockEnvelope(G4bool i) { m_enableRockEnvelope = i; }
     void UpdateGeometry();
 
   private:
@@ -54,6 +55,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction {
     G4bool m_addFORMOSA;
     G4bool m_addFASERnu2;
     G4bool m_addFASER2;
+    G4bool m_enableRockEnvelope;
 
     // FLArE
     G4LogicalVolume* TPCModuleLogical;

include/DetectorConstructionMessenger.hh

@@ -39,6 +39,7 @@ class DetectorConstructionMessenger: public G4UImessenger {
     G4UIcmdWithABool* detGdmlSaveCmd;
     G4UIcmdWithAString* detGdmlFileCmd;
     G4UIcmdWithABool* detCheckOverlapCmd;
+    G4UIcmdWithABool* detEnableRockCmd;
     G4UIcmdWithABool* detAddFLArECmd;
     G4UIcmdWithABool* detAddFORMOSACmd;
     G4UIcmdWithABool* detAddFASERnu2Cmd;

include/geometry/GeometricalParameters.hh

@@ -23,6 +23,13 @@ class GeometricalParameters  {
     G4double GetHallOffsetX() { return fHallOffsetX; }
     G4double GetHallOffsetY() { return fHallOffsetY; }
 
+    // rock envelope
+    G4bool GetEnableRockEnvelope() { return fEnableRockEnvelope; }
+    void SetEnableRockEnvelope(G4bool val) { fEnableRockEnvelope = val; }
+    G4double GetRockFrontThickness() { return fRockFrontThickness; }
+    G4double GetRockSideThickness() { return fRockSideThickness; }
+    G4double GetRockBackThickness() { return fRockBackThickness; }
+
     // FLArE TPC volume
     enum tpcMaterialOption { LiquidArgon, LiquidKrypton};
     tpcMaterialOption ConvertStringToTPCMaterialOption(G4String val);
@@ -223,6 +230,12 @@ class GeometricalParameters  {
     G4double fHallOffsetX; // x offset of hall center from the LOS
     G4double fHallOffsetY; // x offset of hall center from the LOS
 
+    // rock envelope
+    G4bool fEnableRockEnvelope;
+    G4double fRockFrontThickness;
+    G4double fRockSideThickness;
+    G4double fRockBackThickness;
+
     // FLArE TPC volume
     tpcMaterialOption fFLArETPCMaterialOption;
     tpcConfigOption   fFLArETPCConfigOption;

src/DetectorConstruction.cc

@@ -36,6 +36,7 @@
 #include <G4UniformMagField.hh>
 #include <G4FieldManager.hh>
 #include <G4GDMLParser.hh>
+#include <G4SubtractionSolid.hh>
 
 using namespace std;
 
@@ -47,8 +48,9 @@ G4ThreadLocal BabyMINDMagneticField* DetectorConstruction::babyMINDField = 0;
 G4ThreadLocal G4FieldManager* DetectorConstruction::babyMINDFieldMgr = 0;
 
 DetectorConstruction::DetectorConstruction()
-  : G4VUserDetectorConstruction(), 
-    m_addFLArE(true), m_addFORMOSA(true), m_addFASERnu2(true), m_addFASER2(true), m_useBabyMIND(false)
+  : G4VUserDetectorConstruction(),
+    m_addFLArE(true), m_addFORMOSA(true), m_addFASERnu2(true), m_addFASER2(true),
+    m_useBabyMIND(false), m_enableRockEnvelope(false)
 {
   DefineMaterial();
   messenger = new DetectorConstructionMessenger(this);
@@ -57,7 +59,7 @@ DetectorConstruction::DetectorConstruction()
   fCheckOverlap = false;
 }
 
-DetectorConstruction::~DetectorConstruction() 
+DetectorConstruction::~DetectorConstruction()
 {
   delete messenger;
 }
@@ -66,7 +68,7 @@ void DetectorConstruction::DefineMaterial() {
   //-----------------------------
   // construction of materials
   //-----------------------------
-  
+
   LArBoxMaterials = DetectorConstructionMaterial::GetInstance();
 }
 
@@ -79,27 +81,52 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
                                    worldLV,
                                    "worldPV",
                                    nullptr,
-                                   false, 
+                                   false,
                                    0);
 
   // FPF long paper: https://dx.doi.org/10.1088/1361-6471/ac865e
   // section 2.1.1, figure 5, and figure 7
   G4double hallSizeX  = 9.4 * m;
   G4double hallSizeY  = 7.6 * m;
   G4double hallSizeZ  = 64.6 * m;
-  
-  // this offset accounts for: 
-  // - distance between the entrance wall of the hall and the first detector, so the first detector 
+
+  // this offset accounts for:
+  // - distance between the entrance wall of the hall and the first detector, so the first detector
   //   starts at the center of the global coordinate
   // - position of the cavern center w.r.t. the line of sight, since it's not in the exact middle
-  G4ThreeVector hallOffset( GeometricalParameters::Get()->GetHallOffsetX(), 
-                            GeometricalParameters::Get()->GetHallOffsetY(), 
-                            hallSizeZ/2 - GeometricalParameters::Get()->GetHallHeadDistance()); 
-                                                           
+  G4ThreeVector hallOffset( GeometricalParameters::Get()->GetHallOffsetX(),
+                            GeometricalParameters::Get()->GetHallOffsetY(),
+                            hallSizeZ/2 - GeometricalParameters::Get()->GetHallHeadDistance());
+
   auto hallBox = new G4Box("hallBox", hallSizeX/2, hallSizeY/2, hallSizeZ/2);
   hallLV = new G4LogicalVolume(hallBox, LArBoxMaterials->Material("Air"), "hallLV");
   auto hallPV = new G4PVPlacement(nullptr, hallOffset, hallLV, "hallPV", worldLV, false, 0, fCheckOverlap);
 
+  //----------------------------------
+  // Rock envelope
+
+  if(m_enableRockEnvelope){
+
+    G4double rockSizeX = hallSizeX + 2*GeometricalParameters::Get()->GetRockSideThickness();
+    G4double rockSizeY = hallSizeY + 2*GeometricalParameters::Get()->GetRockSideThickness();
+    G4double rockSizeZ = hallSizeZ + GeometricalParameters::Get()->GetRockFrontThickness()
+                                   + GeometricalParameters::Get()->GetRockBackThickness();
+    G4ThreeVector rockOffset( 0, 0, (GeometricalParameters::Get()->GetRockFrontThickness()/2)-(GeometricalParameters::Get()->GetRockBackThickness()/2));
+
+    auto rockBox = new G4Box("rockBox",rockSizeX/2,rockSizeY/2,rockSizeZ/2);
+    auto rockEnvelopeSolid = new G4SubtractionSolid("rockEnvelopeSolid", rockBox, hallBox, 0, rockOffset);
+    auto rockEnvelope = new G4LogicalVolume(rockEnvelopeSolid, LArBoxMaterials->Material("Rock"), "rockEnvelope");
+    auto rockEnvelopePV = new G4PVPlacement(nullptr, hallOffset-rockOffset, rockEnvelope, "rockEnvelopePV", worldLV, false, 0, fCheckOverlap);
+
+    G4cout << "Rock envelope enabled: upstream " << GeometricalParameters::Get()->GetRockFrontThickness() << ","
+           << " downstream " << GeometricalParameters::Get()->GetRockBackThickness() << ","
+           << " side " << GeometricalParameters::Get()->GetRockSideThickness() << G4endl;
+
+    G4VisAttributes* rockVis = new G4VisAttributes(G4Colour(167./255, 168./255, 189./255));
+    rockVis->SetVisibility(true);
+    rockEnvelope->SetVisAttributes(rockVis);       
+  } else G4cout << "Rock envelope disabled" << G4endl;
+
   //-----------------------------------
   // FLArE TPC volume
 
@@ -112,64 +139,64 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
     TPCModuleLogical = FLArETPCAssembler->GetFLArETPCVolume();
 
     // positioning
-    G4double lengthFLArE = 2*TPCInsulationThickness + lArSizeZ; 
+    G4double lengthFLArE = 2*TPCInsulationThickness + lArSizeZ;
     G4ThreeVector FLArEPos = GeometricalParameters::Get()->GetFLArEPosition();
     FLArEPos -= hallOffset;
     new G4PVPlacement(nullptr, FLArEPos, FLArETPCAssembly, "FLArETPCPhysical", hallLV, false, 0, fCheckOverlap);
 
     G4cout << "Length of FLArE     : " << lengthFLArE << G4endl;
     G4cout << "Center of FLArE TPC : " << FLArEPos+hallOffset << G4endl; // w.r.t the global coordinate
-  
+
     //-----------------------------------
     // FLArE HadCal/MuonCatcher or BabyMIND
 
     if( m_useBabyMIND ){   /// use BabyMIND
-  
+
       BabyMINDDetectorConstruction *BabyMINDAssembler = new BabyMINDDetectorConstruction();
       G4LogicalVolume *BabyMINDAssembly = BabyMINDAssembler->GetBabyMINDAssembly();
-  
+
       BabyMINDMagnetPlateLogical = BabyMINDAssembler->GetMagnetPlate();
       BabyMINDVerticalBar = BabyMINDAssembler->GetVerticalBar();
       BabyMINDHorizontalBar = BabyMINDAssembler->GetHorizontalBar();
-    
+
       G4double babyMINDLengthZ  = GeometricalParameters::Get()->GetBabyMINDTotalSizeZ();
       G4ThreeVector babyMINDPos = GeometricalParameters::Get()->GetFLArEPosition() +
                                   G4ThreeVector(0.,0.,lArSizeZ/2.+TPCInsulationThickness) +
                                   G4ThreeVector(0.,0.,babyMINDLengthZ/2.);
       babyMINDPos -= hallOffset;
       new G4PVPlacement(nullptr, babyMINDPos, BabyMINDAssembly, "BabyMINDPhysical", hallLV, false, 0, fCheckOverlap);
-      
+
       G4cout << "Length of BabyMIND : " << babyMINDLengthZ << G4endl;
       G4cout << "Center of BabyMIND : " << babyMINDPos+hallOffset << G4endl; // w.r.t the global coordinate
-    
+
     }
     else{  //legacy HadCal/MuonCatcher
-    
+
       FLArEHadCatcherMuonFinderConstruction *HadCatMuonFindAssembler = new FLArEHadCatcherMuonFinderConstruction();
       G4double HadCatcherLength       = GeometricalParameters::Get()->GetHadCalLength();
       G4double MuonFinderLength       = GeometricalParameters::Get()->GetMuonCatcherLength();
-    
+
       G4LogicalVolume* HadCatMuonFindAssembly = HadCatMuonFindAssembler->GetHadCatcherMuonFinderAssembly();
       HadCalXCellLogical          = HadCatMuonFindAssembler->GetHadCalXVolume();
       HadCalYCellLogical          = HadCatMuonFindAssembler->GetHadCalYVolume();
       HadAbsorLayersLogical       = HadCatMuonFindAssembler->GetHadCalAbsorbVolume();
       MuonFinderXCellLogical      = HadCatMuonFindAssembler->GetMuonCatcherXVolume();
       MuonFinderYCellLogical      = HadCatMuonFindAssembler->GetMuonCatcherYVolume();
       MuonFinderAbsorLayersLogical= HadCatMuonFindAssembler->GetMuonCatcherAbsorbVolume();
-      
+
       G4double HadCatMuonFindLengthZ  = HadCatcherLength + MuonFinderLength;
       G4ThreeVector HadCatMuonFindPos = GeometricalParameters::Get()->GetFLArEPosition() +
 	                                      G4ThreeVector(0.,0.,lArSizeZ/2.+TPCInsulationThickness) +
 				                                G4ThreeVector(0.,0.,HadCatMuonFindLengthZ/2.);
-    
+
       HadCatMuonFindPos -= hallOffset;
       new G4PVPlacement(nullptr, HadCatMuonFindPos, HadCatMuonFindAssembly, "HadCatMuonFindPhysical", hallLV, false, 0, fCheckOverlap);
-      
+
       G4cout << "Length of HadCatherMuonFinder : " << HadCatMuonFindLengthZ << G4endl;
       G4cout << "Center of HadCatherMuonFinder : " << HadCatMuonFindPos+hallOffset << G4endl; // w.r.t the global coordinate
     }
   }
-  
+
   //-----------------------------------
   // FORMOSA
 
@@ -187,7 +214,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
     G4cout<<"Length of FORMOSA : "<<lengthFORMOSA<<G4endl;
     G4cout<<"Center of FORMOSA : "<<FORMOSAPos+hallOffset<<G4endl; // w.r.t the global coordinate
   }
-                         
+
   //-----------------------------------
   // FASERnu2 Emulsion Detector
 
@@ -196,7 +223,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
     FASERnu2EmulsionLogical = FASERnu2Assembler->GetEmulsionFilm();
     FASERnu2VetoInterfaceLogical = FASERnu2Assembler->GetVetoInterfaceDetector();
     G4LogicalVolume* FASERnu2Assembly = FASERnu2Assembler->GetFASERnu2Assembly();
-    
+
     // positioning
     G4double lengthFASERnu2 = GeometricalParameters::Get()->GetFASERnu2TotalSizeZ();
     G4ThreeVector FASERnu2Pos = GeometricalParameters::Get()->GetFASERnu2Position();
@@ -229,7 +256,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
     G4cout<<"Length of FASER2 Spectrometer : "<<lengthFASER2Assembly<<G4endl;
     G4cout<<"Center of FASER2 Spectrometer : "<<FASER2Pos+hallOffset<<G4endl; // w.r.t the global coordinate
   }
-  
+
   //-------------------------------------------------------------------
 
   // visualization
@@ -251,7 +278,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
 }
 
 void DetectorConstruction::ConstructSDandField() {
-  
+
   G4SDManager* sdManager = G4SDManager::GetSDMpointer();
   sdManager->SetVerboseLevel(2);
   int SDIdx = 0;
@@ -263,15 +290,15 @@ void DetectorConstruction::ConstructSDandField() {
     sdManager->AddNewDetector(TPCModuleSD);
     GeometricalParameters::Get()->AddSD2List(SDIdx, "lArBoxSD/lar_box");
     SDIdx++;
-    
+
     if (m_useBabyMIND) {
-    
+
       LArBoxSD* BabyMINDHorBarSD = new LArBoxSD("BabyMINDHorBarSD");
       BabyMINDHorizontalBar->SetSensitiveDetector(BabyMINDHorBarSD);
       sdManager->AddNewDetector(BabyMINDHorBarSD);
       GeometricalParameters::Get()->AddSD2List(SDIdx, "BabyMINDHorBarSD/lar_box");
       SDIdx++;
-      
+
       LArBoxSD* BabyMINDVerBarSD = new LArBoxSD("BabyMINDVerBarSD");
       BabyMINDVerticalBar->SetSensitiveDetector(BabyMINDVerBarSD);
       sdManager->AddNewDetector(BabyMINDVerBarSD);
@@ -284,9 +311,9 @@ void DetectorConstruction::ConstructSDandField() {
       babyMINDFieldMgr->SetDetectorField(babyMINDField);
       babyMINDFieldMgr->CreateChordFinder(babyMINDField);
       BabyMINDMagnetPlateLogical->SetFieldManager(babyMINDFieldMgr, true);
-    
+
     } else {
-    
+
       LArBoxSD* HadCalXSD = new LArBoxSD("HadCalXSD");
       HadCalXCellLogical->SetSensitiveDetector(HadCalXSD);
       sdManager->AddNewDetector(HadCalXSD);
@@ -322,7 +349,7 @@ void DetectorConstruction::ConstructSDandField() {
       sdManager->AddNewDetector(MuonFinderAbsorbSD);
       GeometricalParameters::Get()->AddSD2List(SDIdx, "MuonFinderAbsorbSD/lar_box");
       SDIdx++;
-    
+
       // magnetic field for HadCatcher + MuonFinder
       G4ThreeVector fieldValue = G4ThreeVector(0,fFieldValue, 0);
       magField = new G4UniformMagField(fieldValue);
@@ -352,7 +379,7 @@ void DetectorConstruction::ConstructSDandField() {
     sdManager->AddNewDetector(EmulsionFilmSD);
     GeometricalParameters::Get()->AddSD2List(SDIdx, "FASERnu2EmulsionSD/lar_box");
     SDIdx++;
-    
+
     LArBoxSD* VetoInterfaceSD = new LArBoxSD("FASERnu2VetoInterfaceSD");
     FASERnu2VetoInterfaceLogical->SetSensitiveDetector(VetoInterfaceSD);
     sdManager->AddNewDetector(VetoInterfaceSD);
@@ -406,4 +433,4 @@ void DetectorConstruction::ConstructSDandField() {
   //// define new one
   //G4RunManager::GetRunManager()->DefineWorldVolume(Construct());
 //  G4RunManager::GetRunManager()->GeometryHasBeenModified();
-//}
+//} 

src/DetectorConstructionMessenger.cc

@@ -32,6 +32,12 @@ DetectorConstructionMessenger::DetectorConstructionMessenger(DetectorConstructio
     detCheckOverlapCmd = new G4UIcmdWithABool("/det/checkOverlap", this);
     detCheckOverlapCmd->SetParameterName("checkOverlap", true);
     detCheckOverlapCmd->SetDefaultValue(false);
+
+    // enable rock envelope
+    detEnableRockCmd = new G4UIcmdWithABool("/det/enableRock", this);
+    detEnableRockCmd->SetParameterName("enableRock", true);
+    detEnableRockCmd->SetDefaultValue(false);
+
     // add FLARE volume
     detAddFLArECmd = new G4UIcmdWithABool("/det/addFLArE", this);
     detAddFLArECmd->SetParameterName("Add FLArE detector", true);
@@ -223,6 +229,7 @@ DetectorConstructionMessenger::~DetectorConstructionMessenger() {
   delete detGdmlSaveCmd;
   delete detGdmlFileCmd;
   delete detCheckOverlapCmd;
+  delete detEnableRockCmd;
   delete detAddFLArECmd;
   delete detAddFORMOSACmd;
   delete detAddFASER2Cmd;
@@ -295,6 +302,10 @@ void DetectorConstructionMessenger::SetNewValue(G4UIcommand* command, G4String n
   else if (command == detAddFORMOSACmd) det->AddFORMOSA(detAddFORMOSACmd->GetNewBoolValue(newValues));
   else if (command == detAddFASERnu2Cmd) det->AddFASERnu2(detAddFASERnu2Cmd->GetNewBoolValue(newValues));
   else if (command == detAddFASER2Cmd) det->AddFASER2(detAddFASER2Cmd->GetNewBoolValue(newValues));
+  else if (command == detEnableRockCmd){
+    det->EnableRockEnvelope(detEnableRockCmd->GetNewBoolValue(newValues));
+    GeometricalParameters::Get()->SetEnableRockEnvelope(detEnableRockCmd->GetNewBoolValue(newValues));
+  } 
 
   // FLARE COMMANDS
   else if (command == flarePosCmd)