diff --git a/.github/workflows/build-run.yml b/.github/workflows/build-run.yml index 3e7a5f5b5b..77becc8510 100644 --- a/.github/workflows/build-run.yml +++ b/.github/workflows/build-run.yml @@ -64,8 +64,8 @@ jobs: matrix: vessel_option: [vacuums, helium] SND_option: [all] - muon_shield: [warm_opt, New_HA_Design] - target: [old, Jun25] + muon_shield: [TRY_2025] + target: [Jun25] steps: - name: Checkout diff --git a/CHANGELOG.md b/CHANGELOG.md index 3b279d268c..6cc495d2d4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,9 +14,13 @@ it in future. ### Added -### Changed +* New warm MS option TRY_2025 (Stellatryon v.2) added to the config and its field map +### Changed +* Change function SND in shipDet_conf.py +* Cahnge disable magnetic field in the MS and remove field map upload in run_fixedTarget.py +* Change the ShipMuonShield script to accept new magnet configuration: no fixed number of magnets and variable z-gap between them. * Change naming convention for simulation files to `{sim,geo,params}_{uuid4}.root`, with optional `--tag` parameter to specify custom identifier ### Fixed @@ -27,6 +31,9 @@ it in future. ### Removed +* New_HA_Design and warm_opt muon shield configuration are no longer supported +* CI check for old target version + ## 25.12 - 2025-12-22 ### Added diff --git a/files/TRY_2025.root b/files/TRY_2025.root new file mode 100644 index 0000000000..c94659c8e4 --- /dev/null +++ b/files/TRY_2025.root @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7d2e41500ff46749644e9936ff92cbeb807a8e50459c0aaccb5dc1cdb716c408 +size 21347465 diff --git a/macro/run_fixedTarget.py b/macro/run_fixedTarget.py index 760f7816be..97021441f4 100644 --- a/macro/run_fixedTarget.py +++ b/macro/run_fixedTarget.py @@ -69,7 +69,7 @@ def get_work_dir(run_number,tag=None): ap.add_argument('-o', '--output', type=str, help="output directory", dest='work_dir', default=None) ap.add_argument('-rs', '--seed', type=int, help="random seed; default value is 0, see TRrandom::SetSeed documentation", default=0) ap.add_argument('--DecayVolumeMedium', help='Set Decay Volume Medium. Choices are helium (default) or vacuums.', default='helium', choices=['helium', 'vacuums']) -ap.add_argument('--shieldName', help='Name of the shield in the database. New_HA_Design or warm_opt.', default='New_HA_Design', choices=['New_HA_Design', 'warm_opt']) +ap.add_argument('--shieldName', help='Name of the shield in the database. TRY_2025.', default='TRY_2025', choices=['TRY_2025']) ap.add_argument('--AddMuonShield', help='Whether or not to add the muon shield. Default set to False.', default=False, action=argparse.BooleanOptionalAction) ap.add_argument('--AddMuonShieldField', help='Whether or not to add the muon shield magnetic field. Default set to False.', default=False, action=argparse.BooleanOptionalAction) ap.add_argument('--AddHadronAbsorberOnly', help='Whether to only add the hadron absorber part of the muon shield. Default set to True.', default=True, action=argparse.BooleanOptionalAction) @@ -190,16 +190,14 @@ def get_work_dir(run_number,tag=None): if args.AddMuonShield or args.AddHadronAbsorberOnly: - n_magnets = 7 - n_params = 13 + n_params = 15 if not args.AddMuonShieldField: - for i in range(n_magnets): - ship_geo.muShield.params[n_magnets + i * n_params + 12] = 0 # set B field to 0 + for i in range(ship_geo.muShield.nMagnets): + ship_geo.muShield.params[i * n_params + 14] = 0 # set B field to 0 if args.AddHadronAbsorberOnly: - for i in range(1, n_magnets): - ship_geo.muShield.params[n_magnets + i * n_params] = 0 # set dXIn to 0 + ship_geo.muShield.params = ship_geo.muShield.params[:15] # set dXIn to 0 - MuonShield = ROOT.ShipMuonShield(in_params=list(ship_geo.muShield.params), z=ship_geo.muShield.z, WithConstShieldField=ship_geo.muShield.WithConstField, + MuonShield = ROOT.ShipMuonShield(in_params=list(ship_geo.muShield.params), z=ship_geo.muShield.z, WithConstShieldField=True, SC_key=ship_geo.SC_mag) # MuonShield.SetSupports(False) # otherwise overlap with sensitive Plane run.AddModule(MuonShield) # needs to be added because of magn hadron shield. diff --git a/macro/run_simScript.py b/macro/run_simScript.py index a83c2cfd46..eafd361e14 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -139,7 +139,7 @@ parser.add_argument("-t", "--test", dest="testFlag", help="quick test", action="store_true") parser.add_argument("--dry-run", dest="dryrun", help="stop after initialize", action="store_true") parser.add_argument("-D", "--display", dest="eventDisplay", help="store trajectories", action="store_true") -parser.add_argument("--shieldName", help="The name of the muon shield in the database to use.", default="New_HA_Design", choices=["New_HA_Design", "warm_opt"]) +parser.add_argument("--shieldName", help="The name of the muon shield in the database to use.", default="TRY_2025", choices=["TRY_2025"]) parser.add_argument("--MesonMother", dest="MM", help="Choose DP production meson source: pi0, eta, omega, eta1, eta11", default='pi0') parser.add_argument("--debug", help="Control FairLogger verbosity: 0=info (default), 1=+debug, 2=+debug1, 3=+debug2", default=0, type=int, choices=range(0,4)) parser.add_argument("--print-fields", help="Print VMC fields and weights information", action="store_true") diff --git a/passive/ShipMuonShield.cxx b/passive/ShipMuonShield.cxx index c6ced94d32..1fa85ffb37 100644 --- a/passive/ShipMuonShield.cxx +++ b/passive/ShipMuonShield.cxx @@ -35,19 +35,19 @@ ShipMuonShield::ShipMuonShield(std::vector in_params, Double_t z, const Bool_t WithConstShieldField, const Bool_t SC_key) : FairModule("MuonShield", "ShipMuonShield") { + LOG(debug) << " Function ShipMuonShield MS "; for (size_t i = 0; i < in_params.size(); i++) { shield_params.push_back(in_params[i]); } + nParams = 15; + nMagnets = in_params.size() / nParams; // integer division + + if (in_params.size() % 15 != 0) { + LOG(error) << "Warning: incomplete magnet data!"; + } + fWithConstShieldField = WithConstShieldField; fSC_mag = SC_key; - dZ1 = in_params[0]; - dZ2 = in_params[1]; - dZ3 = in_params[2]; - dZ4 = in_params[3]; - dZ5 = in_params[4]; - dZ6 = in_params[5]; - dZ7 = in_params[6]; - fMuonShieldHalfLength = dZ1 + dZ2 + dZ3 + dZ4 + dZ5 + dZ6 + dZ7; z_end_of_proximity_shielding = z; } @@ -83,9 +83,17 @@ void ShipMuonShield::CreateArb8(TString arbName, TGeoMedium* medium, Double_t z_translation) { TGeoVolume* magVol = nullptr; - if (snd_hole && - ((arbName == "Magn6_MiddleMagL" || arbName == "Magn6_MiddleMagR") || - (arbName == "Magn5_MiddleMagL" || arbName == "Magn5_MiddleMagR"))) { + LOG(debug) << " Create CreateArb8 of the MS "; + + TString magnLast = Form("Magn%zu", nMagnets - 1); + TString magnPrev = Form("Magn%zu", nMagnets - 2); + + bool snd_magnet = (arbName == magnLast + "_MiddleMagL") || + (arbName == magnLast + "_MiddleMagR") || + (arbName == magnPrev + "_MiddleMagL") || + (arbName == magnPrev + "_MiddleMagR"); + + if (snd_hole && snd_magnet) { // // 1) Raw Arb8 “shape” // @@ -120,6 +128,7 @@ void ShipMuonShield::CreateArb8(TString arbName, TGeoMedium* medium, // // 5) Wrap the composite in a volume // + LOG(debug) << " Create CreateArb8 of the MS 5"; magVol = new TGeoVolume(arbName, compShape, medium); } else { // original uncut magnet @@ -143,6 +152,7 @@ void ShipMuonShield::CreateMagnet( Double_t ratio_yoke_2, Double_t dY_yoke_1, Double_t dY_yoke_2, Double_t dZ, Double_t middleGap, Double_t middleGap2, Double_t gap, Double_t gap2, Double_t Z, Bool_t NotMagnet, Bool_t SC_key = false) { + LOG(debug) << " Create a magnet of the MS "; Double_t coil_gap, coil_gap2; Int_t color[4] = {45, 31, 30, 38}; gap = std::ceil(std::max(100. / dY, gap)); @@ -312,31 +322,31 @@ void ShipMuonShield::CreateMagnet( break; } } - -Int_t ShipMuonShield::Initialize( +void ShipMuonShield::Initialize( std::vector& magnetName, std::vector& fieldDirection, std::vector& dXIn, std::vector& dYIn, std::vector& dXOut, std::vector& dYOut, std::vector& ratio_yokesIn, std::vector& ratio_yokesOut, std::vector& dY_yokeIn, std::vector& dY_yokeOut, std::vector& dZ, - std::vector& midGapIn, std::vector& midGapOut, - std::vector& Bgoal, std::vector& gapIn, - std::vector& gapOut, std::vector& Z) { - const Int_t nMagnets = 7; - LOG(info) << " Initialize the MS "; + std::vector& Z_rel, std::vector& midGapIn, + std::vector& midGapOut, std::vector& Bgoal, + std::vector& gapIn, std::vector& gapOut, + std::vector& Z) { + LOG(debug) << " Initialize the MS "; magnetName.reserve(nMagnets); fieldDirection.reserve(nMagnets); - for (auto i : {&dXIn, &dXOut, &dYIn, &dYOut, &dZ, &midGapIn, &midGapOut, - &ratio_yokesIn, &ratio_yokesOut, &dY_yokeIn, &dY_yokeOut, - &Bgoal, &gapIn, &gapOut, &Z}) { + for (auto i : {&dXIn, &dXOut, &dYIn, &dYOut, &dZ, &Z_rel, &midGapIn, + &midGapOut, &ratio_yokesIn, &ratio_yokesOut, &dY_yokeIn, + &dY_yokeOut, &Bgoal, &gapIn, &gapOut, &Z}) { i->reserve(nMagnets); } - - Double_t z_gap = 10 * cm; // fixed distance between magnets in Z-axis - - magnetName = {"MagnAbsorb", "Magn1", "Magn2", "Magn3", - "Magn4", "Magn5", "Magn6"}; + magnetName.push_back("MagnAbsorb"); + for (size_t i = 1; i < nMagnets; ++i) { + TString name = Form("Magn%zu", i); + magnetName.push_back(name); + LOG(debug) << "magnet i: " << name; + } fieldDirection = {FieldDirection::up, FieldDirection::up, FieldDirection::up, FieldDirection::up, @@ -346,43 +356,49 @@ Int_t ShipMuonShield::Initialize( std::vector params; params = shield_params; - const int offset = nMagnets; - const int nParams = 13; - - for (Int_t i = 0; i < nMagnets; ++i) { - dXIn[i] = params[offset + i * nParams + 0]; - dXOut[i] = params[offset + i * nParams + 1]; - dYIn[i] = params[offset + i * nParams + 2]; - dYOut[i] = params[offset + i * nParams + 3]; - gapIn[i] = params[offset + i * nParams + 4]; - gapOut[i] = params[offset + i * nParams + 5]; - ratio_yokesIn[i] = params[offset + i * nParams + 6]; - ratio_yokesOut[i] = params[offset + i * nParams + 7]; - dY_yokeIn[i] = params[offset + i * nParams + 8]; - dY_yokeOut[i] = params[offset + i * nParams + 9]; - midGapIn[i] = params[offset + i * nParams + 10]; - midGapOut[i] = params[offset + i * nParams + 11]; - Bgoal[i] = params[offset + i * nParams + 12]; + for (size_t i = 0; i < nMagnets; ++i) { + // --- Load parameters for each magnet --- + + dZ[i] = params[i * nParams + 0]; + Z_rel[i] = params[i * nParams + 1]; + dXIn[i] = params[i * nParams + 2]; + dXOut[i] = params[i * nParams + 3]; + dYIn[i] = params[i * nParams + 4]; + dYOut[i] = params[i * nParams + 5]; + gapIn[i] = params[i * nParams + 6]; + gapOut[i] = params[i * nParams + 7]; + ratio_yokesIn[i] = params[i * nParams + 8]; + ratio_yokesOut[i] = params[i * nParams + 9]; + dY_yokeIn[i] = params[i * nParams + 10]; + dY_yokeOut[i] = params[i * nParams + 11]; + midGapIn[i] = params[i * nParams + 12]; + midGapOut[i] = params[i * nParams + 13]; + Bgoal[i] = params[i * nParams + 14]; + + // --- Compute Z position for each magnet --- + if (i == 0) { + // First magnet uses the initial offset + Z[i] = z_end_of_proximity_shielding + dZ[i] + Z_rel[i]; + } else { + // Subsequent magnets are placed relative to the previous one + Z[i] = Z[i - 1] + Z_rel[i - 1] + dZ[i] + Z_rel[i]; + } + // --- Print all values for this magnet --- + LOG(debug) << "Magnet " << i << ": dZ=" << dZ[i] << ", Z_rel=" << Z_rel[i] + << ", dXIn=" << dXIn[i] << ", dXOut=" << dXOut[i] + << ", dYIn=" << dYIn[i] << ", dYOut=" << dYOut[i] + << ", gapIn=" << gapIn[i] << ", gapOut=" << gapOut[i] + << ", ratio_yokesIn=" << ratio_yokesIn[i] + << ", ratio_yokesOut=" << ratio_yokesOut[i] + << ", dY_yokeIn=" << dY_yokeIn[i] + << ", dY_yokeOut=" << dY_yokeOut[i] + << ", midGapIn=" << midGapIn[i] << ", midGapOut=" << midGapOut[i] + << ", Bgoal=" << Bgoal[i] << ", Z=" << Z[i]; } - - dZ[0] = dZ1 - z_gap / 2; - Z[0] = z_end_of_proximity_shielding + dZ[0] + z_gap; - dZ[1] = dZ2 - z_gap / 2; - Z[1] = Z[0] + dZ[0] + dZ[1] + z_gap; - dZ[2] = dZ3 - z_gap / 2; - Z[2] = Z[1] + dZ[1] + dZ[2] + 2 * z_gap; - dZ[3] = dZ4 - z_gap / 2; - Z[3] = Z[2] + dZ[2] + dZ[3] + z_gap; - dZ[4] = dZ5 - z_gap / 2; - Z[4] = Z[3] + dZ[3] + dZ[4] + z_gap; - dZ[5] = dZ6 - z_gap / 2; - Z[5] = Z[4] + dZ[4] + dZ[5] + z_gap; - dZ[6] = dZ7 - z_gap / 2; - Z[6] = Z[5] + dZ[5] + dZ[6] + z_gap; - - return nMagnets; + LOG(debug) << " ending of Initialize the MS "; } void ShipMuonShield::ConstructGeometry() { + LOG(debug) << " Construct Geometry of the MS "; TGeoVolume* top = gGeoManager->GetTopVolume(); TGeoVolume* tShield = new TGeoVolumeAssembly("MuonShieldArea"); InitMedium("iron"); @@ -390,22 +406,18 @@ void ShipMuonShield::ConstructGeometry() { std::vector magnetName; std::vector fieldDirection; - std::vector dXIn, dYIn, dXOut, dYOut, dZf, midGapIn, midGapOut, - ratio_yokesIn, ratio_yokesOut, dY_yokeIn, dY_yokeOut, gapIn, gapOut, - Bgoal, Z; - const Int_t nMagnets = - Initialize(magnetName, fieldDirection, dXIn, dYIn, dXOut, dYOut, - ratio_yokesIn, ratio_yokesOut, dY_yokeIn, dY_yokeOut, dZf, - midGapIn, midGapOut, Bgoal, gapIn, gapOut, Z); - - // Create TCC8 tunnel around muon shield - - std::array fieldScale = {{1., 1., 1., 1., 1., 1., 1.}}; - for (Int_t nM = 0; nM < (nMagnets); nM++) { - if (dZf[nM] < 1e-5 || dXIn[nM] == 0) { + std::vector dXIn, dYIn, dXOut, dYOut, dZf, Z_relf, midGapIn, + midGapOut, ratio_yokesIn, ratio_yokesOut, dY_yokeIn, dY_yokeOut, gapIn, + gapOut, Bgoal, Z; + Initialize(magnetName, fieldDirection, dXIn, dYIn, dXOut, dYOut, + ratio_yokesIn, ratio_yokesOut, dY_yokeIn, dY_yokeOut, dZf, Z_relf, + midGapIn, midGapOut, Bgoal, gapIn, gapOut, Z); + + for (size_t nM = 0; nM < nMagnets; ++nM) { + if (Z_relf[nM] < 1e-5 || dXIn[nM] == 0) { continue; } - Double_t ironField_s = Bgoal[nM] * fieldScale[nM] * tesla; + Double_t ironField_s = Bgoal[nM] * tesla; TGeoUniformMagField* magFieldIron_s = new TGeoUniformMagField(0., ironField_s, 0.); TGeoUniformMagField* RetField_s = @@ -419,7 +431,7 @@ void ShipMuonShield::ConstructGeometry() { // Create the magnet CreateMagnet(magnetName[nM], iron, tShield, fields_s, fieldDirection[nM], dXIn[nM], dYIn[nM], dXOut[nM], dYOut[nM], ratio_yokesIn[nM], - ratio_yokesOut[nM], dY_yokeIn[nM], dY_yokeOut[nM], dZf[nM], + ratio_yokesOut[nM], dY_yokeIn[nM], dY_yokeOut[nM], Z_relf[nM], midGapIn[nM], midGapOut[nM], gapIn[nM], gapOut[nM], Z[nM], nM == 0, nM == 3 && fSC_mag); } @@ -427,9 +439,8 @@ void ShipMuonShield::ConstructGeometry() { // Place in origin of SHiP coordinate system as subnodes placed correctly top->AddNode(tShield, 1); - Double_t z_gap = 10 * cm; - Double_t absorber_offset = z_gap; - Double_t absorber_half_length = (dZf[0]); + Double_t absorber_offset = dZf[0]; + Double_t absorber_half_length = (Z_relf[0]); // Absorber diff --git a/passive/ShipMuonShield.h b/passive/ShipMuonShield.h index 761a017fec..afaedc0c6e 100644 --- a/passive/ShipMuonShield.h +++ b/passive/ShipMuonShield.h @@ -29,11 +29,11 @@ class ShipMuonShield : public FairModule { void SetSNDSpace(Bool_t hole, Double_t hole_dx, Double_t hole_dy); protected: - Double_t fMuonShieldHalfLength; // FIXME: HA_field to be removed in the next - // workshop meeting Double_t dZ0, dZ1, dZ2, dZ3, dZ4, dZ5, dZ6, dZ7, dXgap, z_end_of_proximity_shielding; + size_t nMagnets; Int_t InitMedium(TString name); + Int_t nParams; Bool_t fWithConstShieldField; Bool_t fSC_mag; std::vector shield_params; @@ -46,18 +46,18 @@ class ShipMuonShield : public FairModule { Double_t x_translation, Double_t y_translation, Double_t z_translation); - Int_t Initialize(std::vector& magnetName, - std::vector& fieldDirection, - std::vector& dXIn, std::vector& dYIn, - std::vector& dXOut, std::vector& dYOut, - std::vector& ratio_yokesIn, - std::vector& ratio_yokesOut, - std::vector& dY_yokeIn, - std::vector& dY_yokeOut, std::vector& dZ, - std::vector& midGapIn, - std::vector& midGapOut, - std::vector& Bgoal, std::vector& gapIn, - std::vector& gapOut, std::vector& Z); + void Initialize(std::vector& magnetName, + std::vector& fieldDirection, + std::vector& dXIn, std::vector& dYIn, + std::vector& dXOut, std::vector& dYOut, + std::vector& ratio_yokesIn, + std::vector& ratio_yokesOut, + std::vector& dY_yokeIn, + std::vector& dY_yokeOut, std::vector& dZ, + std::vector& Z_rel, std::vector& midGapIn, + std::vector& midGapOut, + std::vector& Bgoal, std::vector& gapIn, + std::vector& gapOut, std::vector& Z); void CreateMagnet(TString magnetName, TGeoMedium* medium, TGeoVolume* tShield, TGeoUniformMagField* fields[4], diff --git a/python/ShieldUtils.py b/python/ShieldUtils.py index 2b687606bd..9d86cc9b05 100644 --- a/python/ShieldUtils.py +++ b/python/ShieldUtils.py @@ -2,32 +2,4 @@ # SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration def find_shield_center(ship_geo): - zEndOfPassiveShield = ship_geo.muShield.z - dZ = [None] * 7 - Z = [None] * 7 - zgap = 10. - dZ[0] = ship_geo.muShield.dZ1 - zgap / 2 - Z[0] = zEndOfPassiveShield + dZ[0] + zgap - - dZ[1] = ship_geo.muShield.dZ2 - zgap / 2 - Z[1] = Z[0] + dZ[0] + dZ[1] + zgap - - dZ[2] = ship_geo.muShield.dZ3 - zgap / 2 - Z[2] = Z[1] + dZ[1] + dZ[2] + 2 * zgap - - dZ[3] = ship_geo.muShield.dZ4 - zgap / 2 - Z[3] = Z[2] + dZ[2] + dZ[3] + zgap - - dZ[4] = ship_geo.muShield.dZ5 - zgap / 2 - Z[4] = Z[3] + dZ[3] + dZ[4] + zgap - - dZ[5] = ship_geo.muShield.dZ6 - zgap / 2 - Z[5] = Z[4] + dZ[4] + dZ[5] + zgap - - dZ[6] = ship_geo.muShield.dZ7 - zgap / 2 - Z[6] = Z[5] + dZ[5] + dZ[6] + zgap - - - shield_center = (Z[1] + Z[6] + dZ[6] - dZ[1]) / 2 - shield_half_length = abs((Z[1] - dZ[1]) - (Z[6] + dZ[6])) / 2 - return shield_center, shield_half_length, Z + return ship_geo.muShield.z + ship_geo.muShield.length/2, ship_geo.muShield.length/2 diff --git a/python/geomGeant4.py b/python/geomGeant4.py index 30c9d6b5ef..fffeec13cc 100644 --- a/python/geomGeant4.py +++ b/python/geomGeant4.py @@ -6,6 +6,7 @@ import hepunit as G4Unit import ShieldUtils import ROOT +import os ROOT.gROOT.ProcessLine('#include "Geant4/G4TransportationManager.hh"') ROOT.gROOT.ProcessLine('#include "Geant4/G4FieldManager.hh"') ROOT.gROOT.ProcessLine('#include "Geant4/G4UIterminal.hh"') @@ -155,15 +156,18 @@ def addVMCFields(shipGeo, controlFile = '', verbose = False, withVirtualMC = Tru ROOT.TVector3(0.0, 0.0, shipGeo.Bfield.z)) fieldsList.append('MainSpecMap') - if not shipGeo.hadronAbsorber.WithConstField: + if not shipGeo.muShield.WithConstField: + offset = shipGeo.muShield.Entrance[0] + quadSymm = True + file_name = f"files/{shipGeo.shieldName}.root" + fieldMaker.defineFieldMap('muonShieldField', file_name, + ROOT.TVector3(0.0, 0.0, offset), ROOT.TVector3(0.0, 0.0, 0.0), quadSymm) + fieldsList.append('muonShieldField') + + elif not shipGeo.hadronAbsorber.WithConstField: fieldMaker.defineFieldMap('HadronAbsorberMap','files/FieldHadronStopper_raised_20190411.root', ROOT.TVector3(0.0,0.0,shipGeo.hadronAbsorber.z)) fieldsList.append('HadronAbsorberMap') - if not shipGeo.muShield.WithConstField: - field_center, _ = ShieldUtils.find_shield_center(shipGeo) - fieldMaker.defineFieldMap('muonShieldField', 'files/MuonShieldField.root', - ROOT.TVector3(0.0, 0.0, field_center), ROOT.TVector3(0.0, 0.0, 0.0), True) - fieldsList.append('muonShieldField') # Combine the fields to obtain the global field if len(fieldsList) > 1: fieldMaker.defineComposite('TotalField', *fieldsList) #fieldsList MUST have length <=4 diff --git a/python/geometry_config.py b/python/geometry_config.py index 9b375a11a9..cb26862264 100644 --- a/python/geometry_config.py +++ b/python/geometry_config.py @@ -17,212 +17,112 @@ # The first row is the length of the magnets # The other rows are the transverse dimensions of the magnets: dXIn[i], dXOut[i] , dYIn[i], dYOut[i], gapIn[i], gapOut[i]. shield_db = { - "warm_opt": { + "TRY_2025": { "hybrid": False, - "WithConstField": True, + "WithConstField": False, "params": [ - 231.0, - 208.0, - 207.0, - 281.0, - 172.82, - 212.54, - 168.64, - 50.0, - 50.0, - 119.0, - 119.0, - 2.0, - 2.0, - 1.0, - 1.0, - 50.0, - 50.0, - 0.0, - 0.0, - 1.6, - 72.0, - 51.0, - 29.0, - 46.0, - 10.0, - 7.0, - 1.0, - 1.0, - 72.0, - 51.0, - 0.0, - 0.0, - 1.7, - 54.0, - 38.0, - 46.0, - 122.0, - 14.0, - 9.0, - 1.0, - 1.0, - 54.0, - 38.0, - 0.0, - 0.0, - 1.7, - 10.0, - 31.0, - 35.0, - 31.0, - 51.0, - 11.0, - 1.0, - 1.0, - 0.0, - 31.0, - 0.0, - 0.0, - 1.7, - 3.0, - 32.0, - 54.0, - 24.0, - 8.0, - 8.0, - 3.0, - 1.0, - 1.0, - 32.0, - 0.0, - 0.0, - 1.7, - 22.0, - 32.0, - 209.0, - 35.0, - 8.0, - 13.0, - 1.0, - 1.0, - 22.0, - 32.0, - 0.0, - 0.0, - 1.7, - 33.0, - 77.0, - 85.0, - 241.0, - 9.0, - 26.0, - 1.0, - 1.0, - 33.0, - 77.0, - 0.0, - 0.0, - 1.7, - ], - }, - "New_HA_Design": { - "hybrid": False, - "WithConstField": True, - "params": [ - 231.0 / 2, - 208.0 + 231.0 / 2, - 207.0, - 281.0, - 172.82, - 212.54, - 168.64, - 50.0, - 50.0, - 119.0, - 119.0, - 2.0, - 2.0, - 1.0, - 1.0, - 50.0, - 50.0, - 0.0, - 0.0, - 1.9, - 72.0, - 51.0, - 29.0, - 46.0, - 10.0, - 7.0, - 1.0, - 1.0, - 72.0, - 51.0, - 0.0, - 0.0, - 1.7, - 54.0, - 38.0, - 46.0, - 122.0, - 14.0, - 9.0, - 1.0, - 1.0, - 54.0, - 38.0, - 0.0, - 0.0, - 1.7, - 10.0, - 31.0, - 35.0, - 31.0, - 51.0, - 11.0, - 1.0, - 1.0, - 0.0, - 31.0, - 0.0, - 0.0, - 1.7, - 3.0, - 32.0, - 54.0, - 24.0, - 8.0, - 8.0, - 3.0, - 1.0, - 1.0, - 32.0, - 0.0, - 0.0, - 1.7, - 22.0, - 32.0, - 209.0, - 35.0, - 8.0, - 13.0, - 1.0, - 1.0, - 22.0, - 32.0, - 0.0, - 0.0, - 1.7, - 33.0, - 77.0, - 85.0, - 241.0, - 9.0, - 26.0, - 1.0, - 1.0, - 33.0, - 77.0, - 0.0, - 0.0, - 1.7, + [ + 0, + 115.5, + 50.00, + 50.00, + 119.00, + 119.00, + 2.00, + 2.00, + 1.00, + 1.00, + 50.00, + 50.00, + 0.00, + 0.00, + 1.90, + ], + [ + 20, + 495.00, + 67.10, + 79.92, + 27.00, + 43.00, + 5.00, + 5.00, + 1.38, + 1.06, + 67.10, + 79.92, + 0.00, + 0.00, + 1.90, + ], + [ + 10, + 280.48, + 53.12, + 49.56, + 43.00, + 56.00, + 5.03, + 5.00, + 2.11, + 2.40, + 53.12, + 49.56, + 0.00, + 0.00, + 1.90, + ], + [ + 10, + 232.53, + 2.73, + 3.68, + 56.00, + 56.00, + 5.00, + 5.21, + 60.44, + 45.63, + 2.73, + 3.68, + 0.50, + 0.50, + -1.91, + ], + [ + 10, + 85.00, + 31.00, + 107.12, + 56.00, + 56.00, + 5.27, + 5.00, + 4.55, + 0.63, + 1.00, + 77.12, + 0.00, + 0.00, + -1.91, + ], + [ + 10, + 233.82, + 30.03, + 40.00, + 56.00, + 56.00, + 5.00, + 5.01, + 4.83, + 3.37, + 30.03, + 40.00, + 0.00, + 0.00, + -1.91, + ], ], }, } @@ -330,42 +230,6 @@ def create_config( c.target.z0 = 0 # Origin of SHiP coordinate system c.target.z = c.target.z0 + c.target.length / 2.0 - c.chambers = AttrDict() - magnetIncrease = 100.0 * u.cm - c.muShield = AttrDict() - c.muShield.Field = 1.7 # in units of Tesla expected by ShipMuonShield - c.muShield.LE = ( - 7 * u.m - ) # - 0.5 m air - Goliath: 4.5 m - 0.5 m air - nu-tau mu-det: 3 m - 0.5 m air. finally 10m asked by Giovanni - c.muShield.dZ0 = 1 * u.m - - # zGap to compensate automatic shortening of magnets - zGap = 0.05 * u.m # halflengh of gap - - params = shield_db[shieldName]["params"] - c.muShield.params = params - c.muShield.dZ1 = params[0] - c.muShield.dZ2 = params[1] - c.muShield.dZ3 = params[2] - c.muShield.dZ4 = params[3] - c.muShield.dZ5 = params[4] - c.muShield.dZ6 = params[5] - c.muShield.dZ7 = params[6] - c.muShield.dXgap = 0.0 * u.m - - c.muShield.length = ( - 2 - * ( - c.muShield.dZ1 - + c.muShield.dZ2 - + c.muShield.dZ3 - + c.muShield.dZ4 - + c.muShield.dZ5 - + c.muShield.dZ6 - + c.muShield.dZ7 - ) - + c.muShield.LE - ) c.hadronAbsorber = AttrDict() @@ -377,7 +241,42 @@ def create_config( + 207.5 * u.mm # Distance between hadron absorber and proximity shielding - 10 * u.cm # Remove spacing internal to hadron absorber ) + + # DEFINITION OF THE MUON SHIELD + c.muShield = AttrDict() + c.muShield.z = c.hadronAbsorber.z + + params = shield_db[shieldName]["params"] + c.muShield.params = params + + # MS length + c.muShield.length = sum(line[0] + line[1] * 2 for line in params) + c.muShield.nMagnets = len(params) + + c.muShield.Zgap = [] + c.muShield.half_length = [] + c.muShield.Entrance = [] + + for line in params: + c.muShield.Zgap.append(line[0]) + c.muShield.half_length.append(line[1]) + + # Compute Z position for each magnet + for i in range(len(c.muShield.Zgap)): + if i == 0: + # First magnet uses the initial offset + c.muShield.Entrance.append(c.muShield.z + c.muShield.Zgap[i]) + else: + # Subsequent magnets are placed relative to the previous one + c.muShield.Entrance.append( + c.muShield.Entrance[i - 1] + + c.muShield.half_length[i - 1] * 2 + + c.muShield.Zgap[i] + ) + # Flatten the params list + c.muShield.params = [item for sublist in params for item in sublist] + c.decayVolume = AttrDict() # target absorber muon shield setup, decayVolume.length = nominal EOI length, only kept to define z=0 @@ -390,6 +289,10 @@ def create_config( c.z - 31.450 * u.m ) # Relative position of decay vessel centre to spectrometer magnet c.decayVolume.z0 = c.decayVolume.z - c.decayVolume.length / 2.0 + + c.chambers = AttrDict() + magnetIncrease = 100.0 * u.cm + if strawDesign != 4 and strawDesign != 10: print( "this design ", strawDesign, " is not supported, use strawDesign = 4 or 10" @@ -513,7 +416,7 @@ def create_config( c.hadronAbsorber.WithConstField = shield_db[shieldName][ "WithConstField" - ] # TO BE CHECKED: NOT SURE IT IS NEEDED + ] c.muShield.WithConstField = shield_db[shieldName]["WithConstField"] # for the digitizing step diff --git a/python/shipDet_conf.py b/python/shipDet_conf.py index c79a7534b7..e51209e591 100644 --- a/python/shipDet_conf.py +++ b/python/shipDet_conf.py @@ -6,7 +6,6 @@ import os import shipunit as u from ShipGeoConfig import AttrDict -from ShieldUtils import find_shield_center from array import array import yaml @@ -133,9 +132,8 @@ def configure_snd_mtc(yaml_file, ship_geo): # Initialize detector if ship_geo.mtc_geo.zPosition == "auto": # Get the the center of the *last* magnet - ship_geo.mtc_geo.zPosition = find_shield_center(ship_geo)[2][-1] - print("MTC zPosition set to ", ship_geo.mtc_geo.zPosition) - mtc = ROOT.MTCDetector("MTC", ROOT.kTRUE) + mtc_total_length = (ship_geo.mtc_geo.ironThick + ship_geo.mtc_geo.sciFiThick + ship_geo.mtc_geo.scintThick) * ship_geo.mtc_geo.nLayers + ship_geo.mtc_geo.zPosition = ship_geo.muShield.Entrance[-1] + mtc_total_length / 2 mtc.SetMTCParameters( ship_geo.mtc_geo.width, ship_geo.mtc_geo.height, @@ -160,7 +158,8 @@ def configure_snd_siliconTarget(yaml_file, ship_geo): # Get the the center of the next to last magnet (temporary placement) # Offset placement of detector by 140 cm, magnet is 2* 212.54 cm, # 120 layers at 132 cm will fit, with 140 cm offset final layer within 10 cm of MTC. - ship_geo.SiliconTarget_geo.zPosition = find_shield_center(ship_geo)[2][-2] + 140 + SiliconTarget_total_length = ship_geo.SiliconTarget_geo.targetSpacing * ship_geo.SiliconTarget_geo.nLayers + ship_geo.SiliconTarget_geo.zPosition = ship_geo.muShield.Entrance[-1] - ship_geo.muShield.Zgap[-1] - SiliconTarget_total_length / 2 print("SiliconTarget zPosition set to ", ship_geo.SiliconTarget_geo.zPosition) SiliconTarget = ROOT.SiliconTarget("SiliconTarget", ROOT.kTRUE) SiliconTarget.SetSiliconTargetParameters(