diff --git a/CMakeLists.txt b/CMakeLists.txt index 4a24091..b7e5dd8 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/README.md b/README.md index b4c72e6..dee7a7a 100644 --- a/README.md +++ b/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`| diff --git a/include/DetectorConstruction.hh b/include/DetectorConstruction.hh index 13812da..3a87b3b 100644 --- a/include/DetectorConstruction.hh +++ b/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; diff --git a/include/DetectorConstructionMessenger.hh b/include/DetectorConstructionMessenger.hh index 0cc0930..29686ce 100644 --- a/include/DetectorConstructionMessenger.hh +++ b/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; diff --git a/include/geometry/GeometricalParameters.hh b/include/geometry/GeometricalParameters.hh index db30345..2fe0d65 100644 --- a/include/geometry/GeometricalParameters.hh +++ b/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; diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc index af28ede..c80de0d 100644 --- a/src/DetectorConstruction.cc +++ b/src/DetectorConstruction.cc @@ -36,6 +36,7 @@ #include #include #include +#include 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,7 +81,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct() worldLV, "worldPV", nullptr, - false, + false, 0); // FPF long paper: https://dx.doi.org/10.1088/1361-6471/ac865e @@ -87,19 +89,44 @@ G4VPhysicalVolume* DetectorConstruction::Construct() 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,43 +139,43 @@ 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(); @@ -156,20 +183,20 @@ G4VPhysicalVolume* DetectorConstruction::Construct() 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 : "<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 : "<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(); -//} +//} diff --git a/src/DetectorConstructionMessenger.cc b/src/DetectorConstructionMessenger.cc index b5f2c5e..f7ce3e2 100644 --- a/src/DetectorConstructionMessenger.cc +++ b/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) diff --git a/src/generators/BackgroundGenerator.cc b/src/generators/BackgroundGenerator.cc index f23ff72..29fe4da 100644 --- a/src/generators/BackgroundGenerator.cc +++ b/src/generators/BackgroundGenerator.cc @@ -110,7 +110,7 @@ int BackgroundGenerator::ExtractBackgroundParticles() const // use it to extract a realization... int Nparticles = 0; if(lambda < 100) Nparticles = int(G4Poisson(lambda) + 0.5); - else Nparticles = int(G4RandGauss::shoot(Nparticles, TMath::Sqrt(Nparticles))+0.5); + else Nparticles = int(G4RandGauss::shoot(lambda, TMath::Sqrt(lambda)+0.5); return Nparticles; } @@ -168,7 +168,9 @@ void BackgroundGenerator::GeneratePrimaries(G4Event* anEvent) // then extract the direction (directional cosines) from second histogram // TODO: you should actually be using a 5D histo for full correlations? - double z = -1.*GeometricalParameters::Get()->GetHallHeadDistance(); //entry wall z in mm + double z = (GeometricalParameters::Get()->GetEnableRockEnvelope()) ? + -1.*(GeometricalParameters::Get()->GetHallHeadDistance()+GeometricalParameters::Get()->GetRockFrontThickness()) + : -1.*GeometricalParameters::Get()->GetHallHeadDistance(); //entry wall z in mm double t = 0.; // TODO: imprint the bunch-crossing timing structure? double x, y, E; fhxyE->GetRandom3(x, y, E); //pos in cm, E in GeV diff --git a/src/geometry/GeometricalParameters.cc b/src/geometry/GeometricalParameters.cc index 61013c9..f935f09 100644 --- a/src/geometry/GeometricalParameters.cc +++ b/src/geometry/GeometricalParameters.cc @@ -12,6 +12,12 @@ GeometricalParameters::GeometricalParameters() fHallOffsetX = 1.44*m; fHallOffsetY = 2.21*m; + // rock envelope + fEnableRockEnvelope = false; + fRockFrontThickness = 10*m; + fRockSideThickness = 3*m; + fRockBackThickness = 3*m; + // FLArE TPC volume fFLArETPCMaterialOption = tpcMaterialOption::LiquidArgon; fFLArETPCConfigOption = tpcConfigOption::ThreeBySeven;