From 79b3aeef700a7f49c860f51ead7a9f91f9b6d997 Mon Sep 17 00:00:00 2001
From: Florian Uhlig <f.uhlig@gsi.de>
Date: Wed, 6 Apr 2022 15:31:17 +0200
Subject: [PATCH] Move missing geometry generation macros to the geometry
repository, refs #2474
---
macro/passive/create_platform_v22b.C | 150 -----
macro/psd/create_psdgeo.C | 437 --------------
macro/psd/create_psdgeo_46modules.C | 576 ------------------
macro/psd/create_psdgeo_52modules.C | 777 -------------------------
macro/psd/create_psdgeo_ideal.C | 312 ----------
macro/psd/create_psdgeo_with_hole.C | 838 ---------------------------
6 files changed, 3090 deletions(-)
delete mode 100644 macro/passive/create_platform_v22b.C
delete mode 100644 macro/psd/create_psdgeo.C
delete mode 100644 macro/psd/create_psdgeo_46modules.C
delete mode 100644 macro/psd/create_psdgeo_52modules.C
delete mode 100644 macro/psd/create_psdgeo_ideal.C
delete mode 100644 macro/psd/create_psdgeo_with_hole.C
diff --git a/macro/passive/create_platform_v22b.C b/macro/passive/create_platform_v22b.C
deleted file mode 100644
index 6432ba4ebb..0000000000
--- a/macro/passive/create_platform_v22b.C
+++ /dev/null
@@ -1,150 +0,0 @@
-/* Copyright (C) 2022 Facility of Antiproton and Ion Research in Europe, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Eoin Clerkin [committer] */
-
-/******************************************************************************
- ** based on create_platform_v16.C
- **
- **/
-
-
-#include "TGeoManager.h"
-
-#include <iomanip>
-#include <iostream>
-
-
-using namespace std;
-
-// ------------- Other global variables -----------------------------------
-// ---> TGeoManager (too lazy to write out 'Manager' all the time
-TGeoManager* gGeoMan = NULL; // will be set later
-// ----------------------------------------------------------------------------
-
-// ============================================================================
-// ====== Main function =====
-// ============================================================================
-
-void create_platform_v22b()
-{
-
- // ----- Define platform parts ----------------------------------------------
-
- /** For v22b (SIS 100) **/
- TString geoTag = "v22b";
- Double_t sizeX = 725.0; // symmetric in x
- Double_t sizeY = 234.0; // without rails
- Double_t sizeZ = 297.0; // short version
- Double_t endZ = 410.0; // ends at z = 450.0
-
- // /** For v16b (SIS 300) **/
- // TString geoTag = "v16b";
- // Double_t sizeX = 725.0; // symmetric in x
- // Double_t sizeY = 234.0; // without rails
- // Double_t sizeZ = 455.0; // long version
- // Double_t endZ = 540.0; // ends at z = 450.0
-
- Double_t beamY = 570.0; // nominal beam height
- Double_t posX = 0.0;
- Double_t posY = -beamY + sizeY / 2.; // rest on the floor at -beamY
- Double_t posZ = endZ - sizeZ / 2.; // end at MUCH z-end
-
- // --------------------------------------------------------------------------
-
-
- // ------- Geometry file name (output) ----------------------------------
- TString geoFileName = "platform_";
- geoFileName = geoFileName + geoTag + ".geo.root";
- // --------------------------------------------------------------------------
-
-
- // ------- Open info file -----------------------------------------------
- TString infoFileName = geoFileName;
- infoFileName.ReplaceAll("root", "info");
- fstream infoFile;
- infoFile.open(infoFileName.Data(), fstream::out);
- infoFile << "Platform geometry created with create_platform_v16.C" << std::endl << std::endl;
- // --------------------------------------------------------------------------
-
-
- // ------- Load media from media file -----------------------------------
- FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- TString geoPath = gSystem->Getenv("VMCWORKDIR");
- TString medFile = geoPath + "/geometry/media.geo";
- geoFace->setMediaFile(medFile);
- geoFace->readMedia();
- gGeoMan = gGeoManager;
- // --------------------------------------------------------------------------
-
-
- // ----------------- Get and create the required media -----------------
- FairGeoMedia* geoMedia = geoFace->getMedia();
- FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
-
- // ---> aluminium
- FairGeoMedium* mAluminium = geoMedia->getMedium("aluminium");
- if (!mAluminium) Fatal("Main", "FairMedium aluminium not found");
- geoBuild->createMedium(mAluminium);
- TGeoMedium* aluminium = gGeoMan->GetMedium("aluminium");
- if (!aluminium) Fatal("Main", "Medium aluminium not found");
-
- // // ---> air
- // FairGeoMedium* mAir = geoMedia->getMedium("air");
- // if ( ! mAir ) Fatal("Main", "FairMedium air not found");
- // geoBuild->createMedium(mAir);
- // TGeoMedium* air = gGeoMan->GetMedium("air");
- // if ( ! air ) Fatal("Main", "Medium air not found");
- // --------------------------------------------------------------------------
-
-
- // -------------- Create geometry and top volume -------------------------
- gGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
- gGeoMan->SetName("PLATFORMgeom");
- TGeoVolume* top = new TGeoVolumeAssembly("top");
- gGeoMan->SetTopVolume(top);
- TString platformName = "platform_";
- platformName += geoTag;
- TGeoVolume* platform = new TGeoVolumeAssembly(platformName.Data());
- // --------------------------------------------------------------------------
-
-
- // ----- Create ---------------------------------------------------------
- TGeoBBox* platform_base = new TGeoBBox("", sizeX / 2., sizeY / 2., sizeZ / 2.);
- TGeoVolume* platform_vol = new TGeoVolume("platform", platform_base, aluminium);
- platform_vol->SetLineColor(kBlue);
- platform_vol->SetTransparency(70);
-
- TGeoTranslation* platform_trans = new TGeoTranslation("", posX, posY, posZ);
- platform->AddNode(platform_vol, 1, platform_trans);
-
- infoFile << "sizeX: " << setprecision(2) << sizeX << "sizeY: " << setprecision(2) << sizeY
- << "sizeZ: " << setprecision(2) << sizeZ << endl;
- infoFile << "posX : " << setprecision(2) << posX << "posY : " << setprecision(2) << posY
- << "posZ : " << setprecision(2) << posZ << endl;
-
- // --------------- Finish -----------------------------------------------
- top->AddNode(platform, 1);
- cout << endl << endl;
- gGeoMan->CloseGeometry();
- gGeoMan->CheckOverlaps(0.001);
- gGeoMan->PrintOverlaps();
- gGeoMan->Test();
-
- platform->Export(geoFileName); // an alternative way of writing the platform volume
-
- TFile* geoFile = new TFile(geoFileName, "UPDATE");
- TGeoTranslation* platform_placement = new TGeoTranslation("platform_trans", 0., 0., 0.);
- platform_placement->Write();
- geoFile->Close();
-
- cout << endl;
- cout << "Geometry " << top->GetName() << " written to " << geoFileName << endl;
-
- top->Draw("ogl");
-
- infoFile.close();
-}
-// ============================================================================
-// ====== End of main function =====
-// ============================================================================
diff --git a/macro/psd/create_psdgeo.C b/macro/psd/create_psdgeo.C
deleted file mode 100644
index c2c736b6bb..0000000000
--- a/macro/psd/create_psdgeo.C
+++ /dev/null
@@ -1,437 +0,0 @@
-/* Copyright (C) 2017 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Volker Friese [committer] */
-
-/** @file create_psdgeo.C
- ** @author Volker Friese <v.friese@gsi.de>
- ** @date 16 September 2017
- ** @version 1.0
- **
- ** This macro creates a PSD geometry in ROOT format to be used as input
- ** for the transport simulation. It allows to create module stacks of
- ** varying sizes; the numbers of rows and columns must both be uneven.
- ** The centre module as well as the edge modules are omitted.
- ** The module definition is implemented in the function ConstructModule.
- **
- ** The user can adjust the steering variables to the specific needs.
- **/
-
-
-using std::cout;
-using std::endl;
-using std::fstream;
-using std::right;
-using std::setw;
-
-// Forward declarations
-TGeoVolume* ConstructModule(const char* name, Double_t moduleSize, Int_t nLayers, fstream* info = 0);
-
-
-// ============================================================================
-// ====== Main function =====
-// ============================================================================
-
-void create_psdgeo()
-{
-
- // ----- Steering variables ---------------------------------------------
- const char* geoTag = "v17a"; // Geometry tag
-
- Double_t psdX = 11.; // x position of PSD in cave
- Double_t psdY = 0.; // y position of PSD in cave
- Double_t psdZ = 800.; // z position of PSD in cave (back side)
- Double_t psdRotY = 0.; // Rotation of PSD around y axis [degrees]
- Double_t moduleSize = 20.; // Module size [cm]
- Int_t nModulesX = 7; // Number of modules in a row (x direction)
- Int_t nModulesY = 7; // Number of modules in a row (x direction)
- // --------------------------------------------------------------------------
-
-
- // ---- Number of modules per row/column must be uneven ---------------------
- // Otherwise, the central one (to be skipped) is not defined.
- if (nModulesX % 2 != 1 || nModulesY % 2 != 1) {
- cout << "Error: number of modules per row and column must be uneven! " << nModulesX << " " << nModulesY << endl;
- return;
- }
- // --------------------------------------------------------------------------
-
-
- // ------- Open info file -----------------------------------------------
- TString infoFileName = "psd_";
- infoFileName = infoFileName + geoTag + ".geo.info";
- fstream infoFile;
- infoFile.open(infoFileName.Data(), fstream::out);
- infoFile << "PSD geometry " << geoTag << " created with create_psdgeo.C" << endl << endl;
- infoFile << "Number of modules: " << nModulesX << " x " << nModulesY << endl;
- infoFile << "Module size: " << moduleSize << " cm x " << moduleSize << " cm" << endl;
- infoFile << "PSD translation in cave: (" << psdX << ", " << psdY << ", " << psdZ << ") cm" << endl;
- infoFile << "PSD rotation around y axis: " << psdRotY << " degrees" << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ------- Load media from media file -----------------------------------
- cout << endl << "==> Reading media..." << endl;
- FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- TString geoPath = gSystem->Getenv("VMCWORKDIR");
- TString medFile = geoPath + "/geometry/media.geo";
- geoFace->setMediaFile(medFile);
- geoFace->readMedia();
- // --------------------------------------------------------------------------
-
-
- // ------ Create the required media from the media file ----------------
- FairGeoMedia* geoMedia = geoFace->getMedia();
- FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
-
- // ---> Air
- FairGeoMedium* mAir = geoMedia->getMedium("air");
- if (!mAir) Fatal("Main", "FairMedium air not found");
- geoBuild->createMedium(mAir);
- TGeoMedium* air = gGeoManager->GetMedium("air");
- if (!air) Fatal("Main", "Medium air not found");
- cout << "Created medium: air" << endl;
-
- // ---> Iron
- FairGeoMedium* mIron = geoMedia->getMedium("iron");
- if (!mIron) Fatal("Main", "FairMedium iron not found");
- geoBuild->createMedium(mIron);
- TGeoMedium* iron = gGeoManager->GetMedium("iron");
- if (!iron) Fatal("Main", "Medium iron not found");
- cout << "Created medium: iron" << endl;
-
- // ---> Lead
- FairGeoMedium* mLead = geoMedia->getMedium("lead");
- if (!mIron) Fatal("Main", "FairMedium lead not found");
- geoBuild->createMedium(mLead);
- TGeoMedium* lead = gGeoManager->GetMedium("lead");
- if (!lead) Fatal("Main", "Medium iron not found");
- cout << "Created medium: lead" << endl;
-
- // ---> Tyvek
- FairGeoMedium* mTyvek = geoMedia->getMedium("PsdTyvek");
- if (!mTyvek) Fatal("Main", "FairMedium PsdTyvek not found");
- geoBuild->createMedium(mTyvek);
- TGeoMedium* tyvek = gGeoManager->GetMedium("PsdTyvek");
- if (!tyvek) Fatal("Main", "Medium PsdTyvek not found");
- cout << "Created medium: PsdTyvek" << endl;
-
- // ---> Scintillator
- FairGeoMedium* mScint = geoMedia->getMedium("PsdScint");
- if (!mScint) Fatal("Main", "FairMedium PsdScint not found");
- geoBuild->createMedium(mScint);
- TGeoMedium* scint = gGeoManager->GetMedium("PsdScint");
- if (!scint) Fatal("Main", "Medium PsdScint not found");
- cout << "Created medium: PsdScint" << endl;
-
- // ---> Fibres
- FairGeoMedium* mFibre = geoMedia->getMedium("PsdFibre");
- if (!mFibre) Fatal("Main", "FairMedium PsdFibre not found");
- geoBuild->createMedium(mFibre);
- TGeoMedium* fibre = gGeoManager->GetMedium("PsdFibre");
- if (!fibre) Fatal("Main", "Medium PsdFibre not found");
- cout << "Created medium: PsdFibre" << endl;
-
- // ---> Medium for PSD volume
- TGeoMedium* psdMedium = air;
- // --------------------------------------------------------------------------
-
-
- // -------------- Create geometry and top volume -------------------------
- cout << endl << "==> Set top volume..." << endl;
- gGeoManager->SetName("PSDgeom");
- TGeoVolume* top = new TGeoVolumeAssembly("TOP");
- gGeoManager->SetTopVolume(top);
- // --------------------------------------------------------------------------
-
-
- // ----- Create the PSD modules -----------------------------------------
- cout << endl;
- TGeoVolume* module2060 = ConstructModule("module2060", 20., 60, &infoFile);
- // --------------------------------------------------------------------------
-
-
- // --------------- Create PSD volume ------------------------------------
- cout << endl;
- cout << "===> Creating PSD...." << endl;
-
- // --- Determine size of PSD box (half-lengths)
- Double_t psdSizeX = 0.5 * nModulesX * moduleSize;
- Double_t psdSizeY = 0.5 * nModulesY * moduleSize;
- TGeoBBox* shape = dynamic_cast<TGeoBBox*>(module2060->GetShape());
- Double_t psdSizeZ = shape->GetDZ();
-
- TString psdName = "psd_";
- psdName += geoTag;
- TGeoVolume* psd = gGeoManager->MakeBox(psdName, psdMedium, psdSizeX, psdSizeY, psdSizeZ);
- cout << "Module array is " << nModulesX << " x " << nModulesY << endl;
- cout << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm" << endl;
- infoFile << endl
- << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm"
- << endl;
- // --------------------------------------------------------------------------
-
-
- // ----- Place modules into the system volume ---------------------------
- Double_t xModCurrent = -0.5 * (nModulesX + 1) * moduleSize;
-
- Int_t nModules = 0;
- for (Int_t iModX = 0; iModX < nModulesX; iModX++) {
- xModCurrent += moduleSize;
- Double_t yModCurrent = -0.5 * (nModulesY + 1) * moduleSize;
- for (Int_t iModY = 0; iModY < nModulesY; iModY++) {
- yModCurrent += moduleSize;
-
- // Skip edge modules
- if ((iModX == 0 || iModX == nModulesX - 1) && (iModY == 0 || iModY == nModulesY - 1)) continue;
-
- // Skip centre module (only for uneven number of modules)
- if (iModX == (nModulesX - 1) / 2 && iModY == (nModulesY - 1) / 2) continue;
-
- // Node number and name (convention)
- Int_t iNode = 100 * iModX + iModY;
-
- // Position of module inside PSD
- TGeoTranslation* trans = new TGeoTranslation(xModCurrent, yModCurrent, 0.);
-
- psd->AddNode(module2060, iNode, trans);
- cout << "Adding module " << setw(5) << right << iNode << " at x = " << setw(6) << right << xModCurrent
- << " cm , y = " << setw(6) << right << yModCurrent << " cm" << endl;
- nModules++;
-
- } //# modules in y direction
- } //# modules in x direction
- cout << "PSD contains " << nModules << " modules." << endl;
- infoFile << "PSD contains " << nModules << " modules." << endl;
- // --------------------------------------------------------------------------
-
-
- // ----- Place PSD in top node (cave) -------------------------------------
- TGeoRotation* psdRot = new TGeoRotation();
- psdRot->RotateY(psdRotY);
- TGeoCombiTrans* psdTrans = new TGeoCombiTrans(psdX, psdY, psdZ + psdSizeZ, psdRot);
- top->AddNode(psd, 0, psdTrans);
- cout << endl << "==> PSD position in cave: " << endl;
- psdTrans->Print();
- // --------------------------------------------------------------------------
-
-
- // ----- Close and check geometry ---------------------------------------
- cout << endl << "==> Closing geometry..." << endl;
- gGeoManager->CloseGeometry();
- gGeoManager->CheckGeometryFull();
- cout << endl;
- gGeoManager->CheckOverlaps(0.0001, "s");
- // --------------------------------------------------------------------------
-
-
- // ----- Write PSD volume and placement matrix to geometry file ---------
- cout << endl;
- TString geoFileName = "psd_";
- geoFileName = geoFileName + geoTag + ".geo.root";
- psd->Export(geoFileName);
- TFile* geoFile = new TFile(geoFileName, "UPDATE");
- psdTrans->Write();
- cout << endl;
- cout << "==> Geometry " << psd->GetName() << " written to " << geoFileName << endl;
- cout << "==> Info written to " << infoFileName << endl;
- geoFile->Close();
- infoFile.close();
- // --------------------------------------------------------------------------
-
-
- // ----- Write entire TGeoManager to file -------------------------------
- TString geoManFileName = "psd_";
- geoManFileName = geoManFileName + geoTag + ".geoman.root";
- TFile* geoManFile = new TFile(geoManFileName, "RECREATE");
- gGeoManager->Write();
- geoManFile->Close();
- cout << "==> TGeoManager " << gGeoManager->GetName() << " written to " << geoManFileName << endl << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ---- Display geometry ------------------------------------------------
- gGeoManager->SetVisLevel(5); // Scintillator level
- top->Draw("ogl");
- // --------------------------------------------------------------------------
-}
-// End of main function
-/** ======================================================================= **/
-
-
-/** ======================================================================= **/
-TGeoVolume* ConstructModule(const char* name, Double_t sizeXY, Int_t nLayers, fstream* info)
-{
-
- // The module consists of nLayers of scintillators covered with Tyvek.
- // Between each two scintillators, there is a lead layer (total nLayers -1).
- // At the top, there is a gap to stretch the light fibres. This is modelled
- // by a volume made of the same material as the scintillator (but inactive).
- // The arrangement is surrounded by iron. At the front and at the back,
- // there is a thick iron layer (the back one acting as first absorber)
- // for constructional reasons.
-
- cout << "===> Creating Module " << name << ", size " << sizeXY << " cm with " << nLayers << " layers...." << endl;
-
- // ----- Constructional parameters ----------------------------
- Double_t leadD = 1.60; // Thickness of lead layer
- Double_t scintD = 0.40; // Thickness of scintillator layer
- Double_t tyvekD = 0.02; // Thickness of Tyvek wrapper around scintillator
- Double_t boxDx = 0.15; // Thickness of iron box lateral
- Double_t boxDy = 0.05; // Thickness of iron box top/bottom
- Double_t boxDz = 2.00; // Thickness of iron box front/back
- Double_t chanDy = 0.20; // Height of fibre channel
- Double_t chanDistL = 0.50; // Distance of channel from left edge of lead
- Double_t chanDistR = 2.00; // Distance of channel from right edge of lead
- // ------------------------------------------------------------------------
-
-
- // ----- Info to file -------------------------------------------------
- if (info) {
- (*info) << "Parameters of module " << name << ": " << endl;
- (*info) << "Size: " << sizeXY << " cm x " << sizeXY << " cm" << endl;
- (*info) << "Number of layers: " << nLayers << endl;
- (*info) << "Thickness of lead layers: " << leadD << " cm" << endl;
- (*info) << "Thickness of scintillators: " << scintD << " cm" << endl;
- (*info) << "Thickness of Tyvek wrap: " << tyvekD << " cm" << endl;
- (*info) << "Thickness of iron box: (" << boxDx << " / " << boxDy << " / " << boxDz << ") cm" << endl;
- (*info) << "Height of fibre channel: " << chanDy << " cm" << endl;
- (*info) << "Distance of channel from edges: left " << chanDistL << " cm, right " << chanDistR << " cm" << endl;
- }
- // ------------------------------------------------------------------------
-
-
- // ----- Get required media -------------------------------------------
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
- TGeoMedium* medLead = gGeoManager->GetMedium("lead"); // Lead layers
- if (!medLead) Fatal("ConstructModule", "Medium lead not found");
- TGeoMedium* medTyvek = gGeoManager->GetMedium("PsdTyvek"); // Tyvek layers
- if (!medTyvek) Fatal("ConstructModule", "Medium PsdTyvek not found");
- TGeoMedium* medScint = gGeoManager->GetMedium("PsdScint"); // Scintillator
- if (!medScint) Fatal("ConstructModule", "Medium PsdScint not found");
- TGeoMedium* medFibre = gGeoManager->GetMedium("PsdFibre"); // Fibres
- if (!medFibre) Fatal("ConstructModule", "Medium PsdFibre not found");
- // ------------------------------------------------------------------------
-
-
- // ----- Create module volume -----------------------------------------
- // Module length: nLayers of scintillators, nLayers - 1 of lead
- // plus the iron box front and back.
- // Material is iron.
- Double_t moduleLength = 2. * boxDz + nLayers * (scintD + 2. * tyvekD) + (nLayers - 1) * leadD;
- TGeoVolume* module = gGeoManager->MakeBox(name, medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * moduleLength);
- module->SetLineColor(kRed);
- module->Print();
- cout << endl;
- // ------------------------------------------------------------------------
-
-
- // ----- Lead filler --------------------------------------------------
- // The lead filler represents all lead absorber layers.
- // The Tyvek/scintillator layers and the fibre channel will be placed
- // inside.
- // Dimensions are half-lengths.
- Double_t leadLx = 0.5 * sizeXY - boxDx;
- Double_t leadLy = 0.5 * sizeXY - boxDy;
- Double_t leadLz = 0.5 * moduleLength - boxDz;
- TGeoVolume* lead = gGeoManager->MakeBox("lead", medLead, leadLx, leadLy, leadLz);
- lead->SetLineColor(kBlue);
- // ------------------------------------------------------------------------
-
-
- // ----- Fibre channel ------------------------------------------------
- // This volume represents the air gap at the top of the volume hosting the
- // light guide fibres. The material is the same as for the scintillators
- // (but inactive).
- // The (half-)width of the channel is defined by the distances from the
- // lead edges.
- Double_t chanLx = leadLx - 0.5 * chanDistL - 0.5 * chanDistR;
- if (chanLx < 0.) {
- cout << "Error: lead volume is not large enough to host fibre channel!" << endl;
- cout << "Lead width: " << 2. * leadLx << ", distance from left edge " << chanDistL << ", distance from right edge "
- << chanDistR << endl;
- return 0;
- }
- Double_t chanLy = 0.5 * chanDy;
- Double_t chanLz = leadLz;
- TGeoVolume* channel = gGeoManager->MakeBox("channel", medFibre, chanLx, chanLy, chanLz);
- channel->SetLineColor(kMagenta);
- // ------------------------------------------------------------------------
-
-
- // ---- Positioning of the fibre channel ------------------------------
- // The channel extends from chanDistL from the left edge of the lead filler
- // to chanDistR from its right edge (seen from the front). It is top-aligned
- // with the lead filler.
- Double_t chanShiftX = 0.5 * (chanDistL - chanDistR);
- Double_t chanShiftY = leadLy - 0.5 * chanDy;
- // ------------------------------------------------------------------------
-
-
- // ----- Tyvek layer --------------------------------------------------
- // The scintillator will be placed inside this, leaving only the thin
- // Tyvek as a wrapper. Since these layers will be placed in the lead filler,
- // there has to be a cut-out for the fibre channel.
- Double_t tyvekLz = 0.5 * scintD + tyvekD; // half-length
- TGeoBBox* tyv1 = new TGeoBBox("psdTyv1", leadLx, leadLy, tyvekLz);
- TGeoBBox* tyv2 = new TGeoBBox("psdTyv2", chanLx, chanLy, tyvekLz);
- TGeoTranslation* trans1 = new TGeoTranslation("tPsd1", chanShiftX, chanShiftY, 0.);
- trans1->RegisterYourself();
- TGeoCompositeShape* tyvekShape = new TGeoCompositeShape("psdTyv1-psdTyv2:tPsd1");
- TGeoVolume* tyvek = new TGeoVolume("tyvek", tyvekShape, medTyvek);
- tyvek->SetLineColor(kGreen);
- // ------------------------------------------------------------------------
-
-
- // ----- Scintillator layer -------------------------------------------
- // The scintillator layer is embedded (as daughter) into the Tyvek layer.
- // It is also a box minus the cut-out for the fibres. The cut-out is
- // slightly smaller than for the Tyvek volume.
- Double_t scintLx = leadLx - tyvekD;
- Double_t scintLy = leadLy - tyvekD;
- Double_t scintLz = 0.5 * scintD;
- TGeoBBox* sci1 = new TGeoBBox("sci1", scintLx, scintLy, scintLz);
- Double_t cutLx = chanLx;
- Double_t cutLy = chanLy - tyvekD;
- Double_t cutLz = scintLz;
- TGeoBBox* sci2 = new TGeoBBox("sci2", cutLx, cutLy, cutLz);
- Double_t cutShiftX = chanShiftX;
- Double_t cutShiftY = scintLy - cutLy;
- TGeoTranslation* trans2 = new TGeoTranslation("tPsd2", cutShiftX, cutShiftY, 0.);
- trans2->RegisterYourself();
- TGeoCompositeShape* scintShape = new TGeoCompositeShape("scintShape", "sci1-sci2:tPsd2");
- TGeoVolume* scint = new TGeoVolume("scint", scintShape, medScint);
- scint->SetLineColor(kBlack);
- // ------------------------------------------------------------------------
-
-
- // ----- Place volumes ------------------------------------------------
-
- // ---> Scintillator into Tyvek
- tyvek->AddNode(scint, 0);
-
- // ---> Fibre channel into lead filler
- lead->AddNode(channel, 0, trans1);
-
- // ---> Tyvek layers into lead
- Double_t zFirst = 0.;
- Double_t zLast = 0.;
- for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
- Double_t zPos = -1. * leadLz + iLayer * leadD + (2. * iLayer + 1.) * tyvekLz;
- if (iLayer == 0) zFirst = zPos;
- if (iLayer == nLayers - 1) zLast = zPos;
- lead->AddNode(tyvek, iLayer, new TGeoTranslation(0., 0., zPos));
- }
- cout << module->GetName() << ": Positioned " << nLayers << " Tyvek layers; first at z = " << zFirst
- << ", last at z = " << zLast << endl;
-
- // ---> Lead into module
- module->AddNode(lead, 0);
- // ------------------------------------------------------------------------
-
- return module;
-}
-/** ======================================================================= **/
diff --git a/macro/psd/create_psdgeo_46modules.C b/macro/psd/create_psdgeo_46modules.C
deleted file mode 100644
index 41494145d2..0000000000
--- a/macro/psd/create_psdgeo_46modules.C
+++ /dev/null
@@ -1,576 +0,0 @@
-/* Copyright (C) 2017-2018 National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Moscow
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Volker Friese, Oleg Golosov [committer] */
-
-/** @file create_psdgeo.C
- ** @author Volker Friese <v.friese@gsi.de>
- ** @date 16 September 2017
- ** @version 1.0
- **
- ** This macro creates a PSD geometry in ROOT format to be used as input
- ** for the transport simulation. It allows to create module stacks of
- ** varying sizes; the numbers of rows and columns must both be uneven.
- ** The centre module as well as the edge modules are omitted.
- ** The module definition is implemented in the function ConstructModule.
- **
- ** The user can adjust the steering variables to the specific needs.
- **/
-
-
-using std::cout;
-using std::endl;
-using std::fstream;
-using std::right;
-using std::setw;
-
-// Forward declarations
-TGeoVolume* ConstructModule(const char* name, Double_t moduleSize, Int_t nLayers, fstream* info = 0);
-
-TGeoVolume* ConstructShield(const char* name, Double_t sizeXY, Double_t holesize, Int_t hole_pos, fstream* info);
-
-// ============================================================================
-// ====== Main function =====
-// ============================================================================
-
-void create_psdgeo_46modules()
-{
- // ----- Steering variables ---------------------------------------------
- const Double_t psdX = 9.65; // x position of PSD in cave (front plane center)
- const Double_t psdY = 0.; // y position of PSD in cave (front plane center)
- const Double_t psdZ = 800.; // z position of PSD in cave (front plane center)
- const Double_t psdRotY = 0.01321; // Rotation of PSD around y axis [rad]
- const Double_t holeSize = 0.; // diagonal size of the square shaped hole
-
- Double_t moduleSize = 20.; // Module size [cm]
- const Double_t shieldWidth = 16.;
-
- TString geoTag = Form("46mod_hole20cm_xshift%4.2fcm", psdX); // Geometry tag
- // --------------------------------------------------------------------------
-
-
- // ------- Open info file -----------------------------------------------
- TString infoFileName = "psd_";
- infoFileName = infoFileName + geoTag + ".geo.info";
- fstream infoFile;
- infoFile.open(infoFileName.Data(), fstream::out);
- infoFile << "PSD geometry " << geoTag << " created with create_psdgeo_46modules.C" << endl << endl;
- infoFile << "PSD front plane center coordinates: (" << psdX << ", " << psdY << ", " << psdZ << ") cm" << endl;
- infoFile << "PSD rotation around y axis: " << psdRotY * TMath::RadToDeg() << " degrees" << endl;
- // --------------------------------------------------------------------------
-
-
- // ------- Load media from media file -----------------------------------
- cout << endl << "==> Reading media..." << endl;
- FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- TString geoPath = gSystem->Getenv("VMCWORKDIR");
- TString medFile = geoPath + "/geometry/media.geo";
- geoFace->setMediaFile(medFile);
- geoFace->readMedia();
- // --------------------------------------------------------------------------
-
-
- // ------ Create the required media from the media file ----------------
- FairGeoMedia* geoMedia = geoFace->getMedia();
- FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
-
- // ---> Air
- FairGeoMedium* mAir = geoMedia->getMedium("air");
- if (!mAir) Fatal("Main", "FairMedium air not found");
- geoBuild->createMedium(mAir);
- TGeoMedium* air = gGeoManager->GetMedium("air");
- if (!air) Fatal("Main", "Medium air not found");
- cout << "Created medium: air" << endl;
-
- // ---> Iron
- FairGeoMedium* mIron = geoMedia->getMedium("iron");
- if (!mIron) Fatal("Main", "FairMedium iron not found");
- geoBuild->createMedium(mIron);
- TGeoMedium* iron = gGeoManager->GetMedium("iron");
- if (!iron) Fatal("Main", "Medium iron not found");
- cout << "Created medium: iron" << endl;
-
- // ---> Lead
- FairGeoMedium* mLead = geoMedia->getMedium("lead");
- if (!mIron) Fatal("Main", "FairMedium lead not found");
- geoBuild->createMedium(mLead);
- TGeoMedium* lead = gGeoManager->GetMedium("lead");
- if (!lead) Fatal("Main", "Medium iron not found");
- cout << "Created medium: lead" << endl;
-
- // ---> Tyvek
- FairGeoMedium* mTyvek = geoMedia->getMedium("PsdTyvek");
- if (!mTyvek) Fatal("Main", "FairMedium PsdTyvek not found");
- geoBuild->createMedium(mTyvek);
- TGeoMedium* tyvek = gGeoManager->GetMedium("PsdTyvek");
- if (!tyvek) Fatal("Main", "Medium PsdTyvek not found");
- cout << "Created medium: PsdTyvek" << endl;
-
- // ---> Scintillator
- FairGeoMedium* mScint = geoMedia->getMedium("PsdScint");
- if (!mScint) Fatal("Main", "FairMedium PsdScint not found");
- geoBuild->createMedium(mScint);
- TGeoMedium* scint = gGeoManager->GetMedium("PsdScint");
- if (!scint) Fatal("Main", "Medium PsdScint not found");
- cout << "Created medium: PsdScint" << endl;
-
- // ---> Fibres
- FairGeoMedium* mFibre = geoMedia->getMedium("PsdFibre");
- if (!mFibre) Fatal("Main", "FairMedium PsdFibre not found");
- geoBuild->createMedium(mFibre);
- TGeoMedium* fibre = gGeoManager->GetMedium("PsdFibre");
- if (!fibre) Fatal("Main", "Medium PsdFibre not found");
- cout << "Created medium: PsdFibre" << endl;
-
- // ---> Shield Polyethylene
- FairGeoMedium* mShieldPoly = geoMedia->getMedium("PsdPolyethylene");
- if (!mShieldPoly) Fatal("Main", "FairMedium PsdPolyethylene not found");
- geoBuild->createMedium(mShieldPoly);
- TGeoMedium* ShieldPoly = gGeoManager->GetMedium("PsdPolyethylene");
- if (!ShieldPoly) Fatal("Main", "Medium PsdPolyethylene not found");
- cout << "Created medium: PsdPolyethylene" << endl;
-
-
- // ---> Medium for PSD volume
- TGeoMedium* psdMedium = air;
- // --------------------------------------------------------------------------
-
-
- // -------------- Create geometry and top volume -------------------------
- cout << endl << "==> Set top volume..." << endl;
- gGeoManager->SetName("PSDgeom");
- TGeoVolume* top = new TGeoVolumeAssembly("TOP");
- gGeoManager->SetTopVolume(top);
- // --------------------------------------------------------------------------
-
-
- // ----- Create the PSD modules -----------------------------------------
- cout << endl;
- TGeoVolume* module2060 = ConstructModule("module2060", 20., 60, &infoFile);
- TGeoVolume* module_shield = ConstructShield("shield20", moduleSize, 0., 0, &infoFile);
- // --------------------------------------------------------------------------
-
-
- // --------------- Create PSD volume ------------------------------------
- cout << endl;
- cout << "===> Creating PSD...." << endl;
-
- // --- Determine size of PSD box (half-lengths)
- Double_t psdSizeX = 180;
- Double_t psdSizeY = 140;
- TGeoBBox* shape = dynamic_cast<TGeoBBox*>(module2060->GetShape());
- Double_t psdSizeZ = shape->GetDZ();
-
- TString psdName = "psd_";
- psdName += geoTag;
- TGeoVolume* psd = gGeoManager->MakeBox(psdName, psdMedium, psdSizeX, psdSizeY, psdSizeZ + 0.5 * shieldWidth);
- cout << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm" << endl;
- infoFile << endl
- << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm"
- << endl;
- // --------------------------------------------------------------------------
-
- struct PSDStatic_t {
- Double_t posX;
- Double_t posY;
- };
-
- std::vector<PSDStatic_t> module2060Positions = {
- // inner part
- {.0, 20.0},
- {20.0, 20.0},
- {20.0, .0},
- {20.0, -20.0},
- {.0, -20.0},
- {-20.0, -20.0},
- {-20.0, .0},
- {-20.0, 20.0},
- // mid-central part
- {.0, 40.0},
- {20.0, 40.0},
- {40.0, 40.0},
- {40.0, 20.0},
- {40.0, .0},
- {40.0, -20.0},
- {40.0, -40.0},
- {20.0, -40.0},
- {.0, -40.0},
- {-20.0, -40.0},
- {-40.0, -40.0},
- {-40.0, -20.0},
- {-40.0, .0},
- {-40.0, 20.0},
- {-40.0, 40.0},
- {-20.0, 40.0},
-
- {-60.0, -40.0},
- {-60.0, -20.0},
- {-60.0, .0},
- {-60.0, 20.0},
- {-60.0, 40.0},
-
- {60.0, -40.0},
- {60.0, -20.0},
- {60.0, .0},
- {60.0, 20.0},
- {60.0, 40.0},
- // outer-top part
- {.0, 60.0},
- {-20.0, 60.0},
- {20.0, 60.0},
- // outer-bottom part
- {.0, -60.0},
- {-20.0, -60.0},
- {20.0, -60.0},
- // outer-left part
- {-80.0, .0},
- {-80.0, -20.0},
- {-80.0, 20.0},
- // outer-right part
- {80.0, .0},
- {80.0, -20.0},
- {80.0, 20.0},
- };
-
-
- // ----- Place modules into the system volume ---------------------------
- Double_t xModCurrent;
- Double_t yModCurrent;
-
- Int_t nModules = 0;
- for (auto mm : module2060Positions) {
- xModCurrent = mm.posX;
- yModCurrent = mm.posY;
- // Node number and name (convention)
- Int_t iNode = nModules;
-
- // Position of module inside PSD
- TGeoTranslation* trans = new TGeoTranslation(xModCurrent, yModCurrent, -0.5 * shieldWidth);
- TGeoTranslation* trans_shield = new TGeoTranslation(xModCurrent, yModCurrent, psdSizeZ);
-
- psd->AddNode(module2060, iNode, trans);
- psd->AddNode(module_shield, iNode, trans_shield);
-
- cout << "Adding module " << setw(5) << right << iNode << " at x = " << setw(6) << right << xModCurrent
- << " cm , y = " << setw(6) << right << yModCurrent << " cm" << endl;
- nModules++;
-
- } //# modules
-
- cout << endl;
-
- cout << "PSD contains " << nModules << " modules." << endl;
- infoFile << "PSD contains " << nModules << " modules." << endl;
- // --------------------------------------------------------------------------
-
-
- // ----- Place PSD in top node (cave) -------------------------------------
- Double_t psdHalfLength = psdSizeZ + 0.5 * shieldWidth;
- Double_t psdVolCenterX = psdX + psdHalfLength * sin(psdRotY);
- Double_t psdVolCenterY = psdY;
- Double_t psdVolCenterZ = psdZ + psdHalfLength * cos(psdRotY);
- infoFile << "PSD volume center coordinates: (" << psdVolCenterX << ", " << psdVolCenterY << ", " << psdVolCenterZ
- << ") cm" << endl;
- TGeoRotation* psdRot = new TGeoRotation();
- psdRot->RotateY(psdRotY * TMath::RadToDeg());
- TGeoCombiTrans* psdTrans = new TGeoCombiTrans(psdVolCenterX, psdVolCenterY, psdVolCenterZ, psdRot);
- top->AddNode(psd, 0, psdTrans);
- cout << endl << "==> PSD position in cave: " << endl;
- psdTrans->Print();
- // --------------------------------------------------------------------------
-
-
- // ----- Close and check geometry ---------------------------------------
- cout << endl << "==> Closing geometry..." << endl;
- gGeoManager->CloseGeometry();
- gGeoManager->CheckGeometryFull();
- cout << endl;
- gGeoManager->CheckOverlaps(0.0001, "s");
- // --------------------------------------------------------------------------
-
-
- // ----- Write PSD volume and placement matrix to geometry file ---------
- cout << endl;
- TString geoFileName = "psd_";
- geoFileName = geoFileName + geoTag + ".geo.root";
- psd->Export(geoFileName);
- TFile* geoFile = new TFile(geoFileName, "UPDATE");
- psdTrans->Write();
- cout << endl;
- cout << "==> Geometry " << psd->GetName() << " written to " << geoFileName << endl;
- cout << "==> Info written to " << infoFileName << endl;
- geoFile->Close();
- infoFile.close();
- // --------------------------------------------------------------------------
-
-
- // ----- Write entire TGeoManager to file -------------------------------
- TString geoManFileName = "psd_";
- geoManFileName = geoManFileName + geoTag + ".geoman.root";
- TFile* geoManFile = new TFile(geoManFileName, "RECREATE");
- gGeoManager->Write();
- geoManFile->Close();
- cout << "==> TGeoManager " << gGeoManager->GetName() << " written to " << geoManFileName << endl << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ---- Display geometry ------------------------------------------------
- gGeoManager->SetVisLevel(5); // Scintillator level
- top->Draw("ogl");
- // --------------------------------------------------------------------------
-}
-// End of main function
-/** ======================================================================= **/
-
-TGeoVolume* ConstructShield(const char* name, Double_t sizeXY, Double_t holesize, Int_t hole_pos, fstream* info)
-{
-
- const TString suffix = Form("_shield_%d_%d_%d", int(sizeXY), int(holesize), int(hole_pos));
-
- const Double_t psd2iron = 2.;
- const Double_t ironD = 4.;
- const Double_t iron2poly = 2.;
- const Double_t polyD = 8.;
- const Double_t shieldLength = psd2iron + ironD + iron2poly + polyD;
-
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
-
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
-
- TGeoMedium* medShieldPoly = gGeoManager->GetMedium("PsdPolyethylene");
- if (!medShieldPoly) Fatal("Main", "Medium PsdPolyethylene not found");
-
- TGeoVolume* shield = gGeoManager->MakeBox(name, medAir, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * shieldLength);
- TGeoVolume* iron {nullptr};
- TGeoVolume* poly {nullptr};
-
- if (holesize > 0) {
- Int_t sx, sy;
- switch (hole_pos) {
- case 0:
- sx = 1;
- sy = 1;
- break; // defines hole position in x-y plane
- case 1:
- sx = 1;
- sy = -1;
- break;
- case 2:
- sx = -1;
- sy = 1;
- break;
- case 3:
- sx = -1;
- sy = -1;
- break;
- }
- TGeoCombiTrans* rot45 = new TGeoCombiTrans(Form("rot45%s", suffix.Data()), sx * sizeXY / 2., sy * sizeXY / 2., 0.,
- new TGeoRotation("rot", 45., 0., 0.));
- rot45->RegisterYourself();
-
- TGeoBBox* iron1 = new TGeoBBox(Form("iron1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * ironD);
- TGeoBBox* iron2 = new TGeoBBox(Form("iron2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * ironD); //hole
-
- TGeoCompositeShape* ironShape =
- new TGeoCompositeShape(Form("iron1%s-iron2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- iron = new TGeoVolume(Form("iron%s", suffix.Data()), ironShape, medIron);
-
- TGeoBBox* poly1 = new TGeoBBox(Form("poly1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * polyD);
- TGeoBBox* poly2 = new TGeoBBox(Form("poly2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * polyD); //hole
-
- TGeoCompositeShape* polyShape =
- new TGeoCompositeShape(Form("poly1%s-poly2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- poly = new TGeoVolume(Form("poly%s", suffix.Data()), polyShape, medShieldPoly);
- }
- else {
- iron = gGeoManager->MakeBox("shield_iron", medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * ironD);
- poly = gGeoManager->MakeBox("shield_poly", medShieldPoly, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * polyD);
- }
-
- shield->AddNode(iron, 0, new TGeoTranslation(0., 0., -shieldLength / 2. + psd2iron + ironD / 2.));
- shield->AddNode(poly, 0, new TGeoTranslation(0., 0., -shieldLength / 2. + psd2iron + ironD + iron2poly + polyD / 2.));
-
- return shield;
-}
-
-
-/** ======================================================================= **/
-TGeoVolume* ConstructModule(const char* name, Double_t sizeXY, Int_t nLayers, fstream* info)
-{
-
- // The module consists of nLayers of scintillators covered with Tyvek.
- // Between each two scintillators, there is a lead layer (total nLayers -1).
- // At the top, there is a gap to stretch the light fibres. This is modelled
- // by a volume made of the same material as the scintillator (but inactive).
- // The arrangement is surrounded by iron. At the front and at the back,
- // there is a thick iron layer (the back one acting as first absorber)
- // for constructional reasons.
-
- cout << "===> Creating Module " << name << ", size " << sizeXY << " cm with " << nLayers << " layers...." << endl;
-
- // ----- Constructional parameters ----------------------------
- Double_t leadD = 1.60; // Thickness of lead layer
- Double_t scintD = 0.40; // Thickness of scintillator layer
- Double_t tyvekD = 0.02; // Thickness of Tyvek wrapper around scintillator
- Double_t boxDx = 0.15; // Thickness of iron box lateral
- Double_t boxDy = 0.05; // Thickness of iron box top/bottom
- Double_t boxDz = 2.00; // Thickness of iron box front/back
- Double_t chanDy = 0.20; // Height of fibre channel
- Double_t chanDistL = 0.50; // Distance of channel from left edge of lead
- Double_t chanDistR = 2.00; // Distance of channel from right edge of lead
- // ------------------------------------------------------------------------
-
-
- // ----- Info to file -------------------------------------------------
- if (info) {
- (*info) << "Parameters of module " << name << ": " << endl;
- (*info) << "Size: " << sizeXY << " cm x " << sizeXY << " cm" << endl;
- (*info) << "Number of layers: " << nLayers << endl;
- (*info) << "Thickness of lead layers: " << leadD << " cm" << endl;
- (*info) << "Thickness of scintillators: " << scintD << " cm" << endl;
- (*info) << "Thickness of Tyvek wrap: " << tyvekD << " cm" << endl;
- (*info) << "Thickness of iron box: (" << boxDx << " / " << boxDy << " / " << boxDz << ") cm" << endl;
- (*info) << "Height of fibre channel: " << chanDy << " cm" << endl;
- (*info) << "Distance of channel from edges: left " << chanDistL << " cm, right " << chanDistR << " cm" << endl;
- }
- // ------------------------------------------------------------------------
-
-
- // ----- Get required media -------------------------------------------
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
- TGeoMedium* medLead = gGeoManager->GetMedium("lead"); // Lead layers
- if (!medLead) Fatal("ConstructModule", "Medium lead not found");
- TGeoMedium* medTyvek = gGeoManager->GetMedium("PsdTyvek"); // Tyvek layers
- if (!medTyvek) Fatal("ConstructModule", "Medium PsdTyvek not found");
- TGeoMedium* medScint = gGeoManager->GetMedium("PsdScint"); // Scintillator
- if (!medScint) Fatal("ConstructModule", "Medium PsdScint not found");
- TGeoMedium* medFibre = gGeoManager->GetMedium("PsdFibre"); // Fibres
- if (!medFibre) Fatal("ConstructModule", "Medium PsdFibre not found");
- // ------------------------------------------------------------------------
-
-
- // ----- Create module volume -----------------------------------------
- // Module length: nLayers of scintillators, nLayers - 1 of lead
- // plus the iron box front and back.
- // Material is iron.
- Double_t moduleLength = 2. * boxDz + nLayers * (scintD + 2. * tyvekD) + (nLayers - 1) * leadD;
- TGeoVolume* module = gGeoManager->MakeBox(name, medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * moduleLength);
- module->SetLineColor(kRed);
- module->Print();
- cout << endl;
- // ------------------------------------------------------------------------
-
-
- // ----- Lead filler --------------------------------------------------
- // The lead filler represents all lead absorber layers.
- // The Tyvek/scintillator layers and the fibre channel will be placed
- // inside.
- // Dimensions are half-lengths.
- Double_t leadLx = 0.5 * sizeXY - boxDx;
- Double_t leadLy = 0.5 * sizeXY - boxDy;
- Double_t leadLz = 0.5 * moduleLength - boxDz;
- TGeoVolume* lead = gGeoManager->MakeBox("lead", medLead, leadLx, leadLy, leadLz);
- lead->SetLineColor(kBlue);
- // ------------------------------------------------------------------------
-
-
- // ----- Fibre channel ------------------------------------------------
- // This volume represents the air gap at the top of the volume hosting the
- // light guide fibres. The material is the same as for the scintillators
- // (but inactive).
- // The (half-)width of the channel is defined by the distances from the
- // lead edges.
- Double_t chanLx = leadLx - 0.5 * chanDistL - 0.5 * chanDistR;
- if (chanLx < 0.) {
- cout << "Error: lead volume is not large enough to host fibre channel!" << endl;
- cout << "Lead width: " << 2. * leadLx << ", distance from left edge " << chanDistL << ", distance from right edge "
- << chanDistR << endl;
- return 0;
- }
- Double_t chanLy = 0.5 * chanDy;
- Double_t chanLz = leadLz;
- TGeoVolume* channel = gGeoManager->MakeBox("channel", medFibre, chanLx, chanLy, chanLz);
- channel->SetLineColor(kMagenta);
- // ------------------------------------------------------------------------
-
-
- // ---- Positioning of the fibre channel ------------------------------
- // The channel extends from chanDistL from the left edge of the lead filler
- // to chanDistR from its right edge (seen from the front). It is top-aligned
- // with the lead filler.
- Double_t chanShiftX = 0.5 * (chanDistL - chanDistR);
- Double_t chanShiftY = leadLy - 0.5 * chanDy;
- // ------------------------------------------------------------------------
-
-
- // ----- Tyvek layer --------------------------------------------------
- // The scintillator will be placed inside this, leaving only the thin
- // Tyvek as a wrapper. Since these layers will be placed in the lead filler,
- // there has to be a cut-out for the fibre channel.
- Double_t tyvekLz = 0.5 * scintD + tyvekD; // half-length
- TGeoBBox* tyv1 = new TGeoBBox("psdTyv1", leadLx, leadLy, tyvekLz);
- TGeoBBox* tyv2 = new TGeoBBox("psdTyv2", chanLx, chanLy, tyvekLz);
- TGeoTranslation* trans1 = new TGeoTranslation("tPsd1", chanShiftX, chanShiftY, 0.);
- trans1->RegisterYourself();
- TGeoCompositeShape* tyvekShape = new TGeoCompositeShape("psdTyv1-psdTyv2:tPsd1");
- TGeoVolume* tyvek = new TGeoVolume("tyvek", tyvekShape, medTyvek);
- tyvek->SetLineColor(kGreen);
- // ------------------------------------------------------------------------
-
-
- // ----- Scintillator layer -------------------------------------------
- // The scintillator layer is embedded (as daughter) into the Tyvek layer.
- // It is also a box minus the cut-out for the fibres. The cut-out is
- // slightly smaller than for the Tyvek volume.
- Double_t scintLx = leadLx - tyvekD;
- Double_t scintLy = leadLy - tyvekD;
- Double_t scintLz = 0.5 * scintD;
- TGeoBBox* sci1 = new TGeoBBox("sci1", scintLx, scintLy, scintLz);
- Double_t cutLx = chanLx;
- Double_t cutLy = chanLy - tyvekD;
- Double_t cutLz = scintLz;
- TGeoBBox* sci2 = new TGeoBBox("sci2", cutLx, cutLy, cutLz);
- Double_t cutShiftX = chanShiftX;
- Double_t cutShiftY = scintLy - cutLy;
- TGeoTranslation* trans2 = new TGeoTranslation("tPsd2", cutShiftX, cutShiftY, 0.);
- trans2->RegisterYourself();
- TGeoCompositeShape* scintShape = new TGeoCompositeShape("scintShape", "sci1-sci2:tPsd2");
- TGeoVolume* scint = new TGeoVolume("scint", scintShape, medScint);
- scint->SetLineColor(kBlack);
- // ------------------------------------------------------------------------
-
-
- // ----- Place volumes ------------------------------------------------
-
- // ---> Scintillator into Tyvek
- tyvek->AddNode(scint, 0);
-
- // ---> Fibre channel into lead filler
- lead->AddNode(channel, 0, trans1);
-
- // ---> Tyvek layers into lead
- Double_t zFirst = 0.;
- Double_t zLast = 0.;
- for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
- Double_t zPos = -1. * leadLz + iLayer * leadD + (2. * iLayer + 1.) * tyvekLz;
- if (iLayer == 0) zFirst = zPos;
- if (iLayer == nLayers - 1) zLast = zPos;
- lead->AddNode(tyvek, iLayer, new TGeoTranslation(0., 0., zPos));
- }
- cout << module->GetName() << ": Positioned " << nLayers << " Tyvek layers; first at z = " << zFirst
- << ", last at z = " << zLast << endl;
-
- // ---> Lead into module
- module->AddNode(lead, 0);
- // ------------------------------------------------------------------------
-
- return module;
-}
-/** ======================================================================= **/
diff --git a/macro/psd/create_psdgeo_52modules.C b/macro/psd/create_psdgeo_52modules.C
deleted file mode 100644
index fc6e63b7c4..0000000000
--- a/macro/psd/create_psdgeo_52modules.C
+++ /dev/null
@@ -1,777 +0,0 @@
-/* Copyright (C) 2017-2018 National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Moscow
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Volker Friese, Evgeny Kashirin [committer] */
-
-/** @file create_psdgeo_with_hole.C
- ** @author Volker Friese <v.friese@gsi.de>
- ** @date 16 September 2017
- ** @version 1.0
- **
- ** This macro creates a PSD geometry in ROOT format to be used as input
- ** for the transport simulation. It allows to create module stacks of
- ** varying sizes; the numbers of rows and columns must both be uneven.
- ** The centre module as well as the edge modules are omitted.
- ** The module definition is implemented in the function ConstructModule.
- **
- ** The user can adjust the steering variables to the specific needs.
- **/
-
-#include <iostream>
-
-
-using std::cout;
-using std::endl;
-using std::fstream;
-using std::right;
-using std::setw;
-
-// Forward declarations
-TGeoVolume* ConstructModule(const char* name, Double_t moduleSize, Int_t nLayers, fstream* info = 0);
-
-TGeoVolume* ConstructModuleWithHole(const char* name, Double_t moduleSize, Int_t nLayers, fstream* info = 0,
- float holesize = 20, Int_t hole_pos = 0);
-
-TGeoVolume* ConstructShield(const char* name, Double_t sizeXY, Double_t holesize, Int_t hole_pos, fstream* info);
-
-// ============================================================================
-// ====== Main function =====
-// ============================================================================
-
-void create_psdgeo_52modules()
-{
- // ----- Steering variables ---------------------------------------------
- const Double_t psdX = 9.65; // x position of PSD in cave (front plane center)
- const Double_t psdY = 0.; // y position of PSD in cave (front plane center)
- const Double_t psdZ = 800.; // z position of PSD in cave (front plane center)
- const Double_t psdRotY = 0.01321; // Rotation of PSD around y axis [rad]
-
- const Double_t bigModuleSize = 20.; // Module size [cm]
- const Double_t smallModuleSize = 10.; // Module size [cm]
- const Int_t nModulesX = 8; // Number of modules in a row (x direction)
- const Int_t nModulesY = 6; // Number of modules in a row (x direction)
- const Int_t nLayers = 60; // Number of sections in a module (z direction)
-
- TString geoTag = Form("52mod_hole20cm_xshift%4.2fcm", psdX); // Geometry tag
-
- const Double_t shieldWidth = 16.;
- // --------------------------------------------------------------------------
-
- // ---- Number of modules per row/column must be uneven ---------------------
- // Otherwise, the central one (to be skipped) is not defined.
- // if ( nModulesX % 2 != 1 || nModulesY % 2 != 1 ) {
- cout << "Number of modules per row and column = " << nModulesX << " " << nModulesY << endl;
- // return;
- // }
- // --------------------------------------------------------------------------
-
-
- // ------- Open info file -----------------------------------------------
- TString infoFileName = "psd_";
- infoFileName = infoFileName + geoTag + ".geo.info";
- fstream infoFile;
- infoFile.open(infoFileName.Data(), fstream::out);
- infoFile << "PSD geometry " << geoTag << " created with create_psdgeo_52modules.C" << endl << endl;
- infoFile << "Number of modules: " << nModulesX << " x " << nModulesY << endl;
- infoFile << "Small module size: " << smallModuleSize << " cm x " << smallModuleSize << " cm" << endl;
- infoFile << "Big module size: " << bigModuleSize << " cm x " << bigModuleSize << " cm" << endl;
- infoFile << "PSD front plane center coordinates: (" << psdX << ", " << psdY << ", " << psdZ << ") cm" << endl;
- infoFile << "PSD rotation around y axis: " << psdRotY * TMath::RadToDeg() << " degrees" << endl;
- infoFile << "Size of the square shaped hole in PSD center: 20 cm" << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ------- Load media from media file -----------------------------------
- cout << endl << "==> Reading media..." << endl;
- FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- TString geoPath = gSystem->Getenv("VMCWORKDIR");
- TString medFile = geoPath + "/geometry/media.geo";
- geoFace->setMediaFile(medFile);
- geoFace->readMedia();
- // --------------------------------------------------------------------------
-
-
- // ------ Create the required media from the media file ----------------
- FairGeoMedia* geoMedia = geoFace->getMedia();
- FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
-
- // ---> Air
- FairGeoMedium* mAir = geoMedia->getMedium("air");
- if (!mAir) Fatal("Main", "FairMedium air not found");
- geoBuild->createMedium(mAir);
- TGeoMedium* air = gGeoManager->GetMedium("air");
- if (!air) Fatal("Main", "Medium air not found");
- cout << "Created medium: air" << endl;
-
- // ---> Iron
- FairGeoMedium* mIron = geoMedia->getMedium("iron");
- if (!mIron) Fatal("Main", "FairMedium iron not found");
- geoBuild->createMedium(mIron);
- TGeoMedium* iron = gGeoManager->GetMedium("iron");
- if (!iron) Fatal("Main", "Medium iron not found");
- cout << "Created medium: iron" << endl;
-
- // ---> Lead
- FairGeoMedium* mLead = geoMedia->getMedium("lead");
- if (!mIron) Fatal("Main", "FairMedium lead not found");
- geoBuild->createMedium(mLead);
- TGeoMedium* lead = gGeoManager->GetMedium("lead");
- if (!lead) Fatal("Main", "Medium iron not found");
- cout << "Created medium: lead" << endl;
-
- // ---> Tyvek
- FairGeoMedium* mTyvek = geoMedia->getMedium("PsdTyvek");
- if (!mTyvek) Fatal("Main", "FairMedium PsdTyvek not found");
- geoBuild->createMedium(mTyvek);
- TGeoMedium* tyvek = gGeoManager->GetMedium("PsdTyvek");
- if (!tyvek) Fatal("Main", "Medium PsdTyvek not found");
- cout << "Created medium: PsdTyvek" << endl;
-
- // ---> Scintillator
- FairGeoMedium* mScint = geoMedia->getMedium("PsdScint");
- if (!mScint) Fatal("Main", "FairMedium PsdScint not found");
- geoBuild->createMedium(mScint);
- TGeoMedium* scint = gGeoManager->GetMedium("PsdScint");
- if (!scint) Fatal("Main", "Medium PsdScint not found");
- cout << "Created medium: PsdScint" << endl;
-
- // ---> Fibres
- FairGeoMedium* mFibre = geoMedia->getMedium("PsdFibre");
- if (!mFibre) Fatal("Main", "FairMedium PsdFibre not found");
- geoBuild->createMedium(mFibre);
- TGeoMedium* fibre = gGeoManager->GetMedium("PsdFibre");
- if (!fibre) Fatal("Main", "Medium PsdFibre not found");
- cout << "Created medium: PsdFibre" << endl;
-
- // ---> Shield Polyethylene
- FairGeoMedium* mShieldPoly = geoMedia->getMedium("PsdPolyethylene");
- if (!mShieldPoly) Fatal("Main", "FairMedium PsdPolyethylene not found");
- geoBuild->createMedium(mShieldPoly);
- TGeoMedium* ShieldPoly = gGeoManager->GetMedium("PsdPolyethylene");
- if (!ShieldPoly) Fatal("Main", "Medium PsdPolyethylene not found");
- cout << "Created medium: PsdPolyethylene" << endl;
-
-
- // ---> Medium for PSD volume
- TGeoMedium* psdMedium = air;
- // --------------------------------------------------------------------------
-
-
- // -------------- Create geometry and top volume -------------------------
- cout << endl << "==> Set top volume..." << endl;
- gGeoManager->SetName("PSDgeom");
- TGeoVolume* top = new TGeoVolumeAssembly("TOP");
- gGeoManager->SetTopVolume(top);
- // --------------------------------------------------------------------------
-
- // ----- Create the PSD modules -----------------------------------------
- cout << endl;
- TGeoVolume* module2060 = ConstructModule("module2060", bigModuleSize, nLayers, &infoFile);
- TGeoVolume* module1060 = ConstructModule("module1060", smallModuleSize, nLayers, &infoFile);
- TGeoVolume* module_shield20 = ConstructShield("shield20", bigModuleSize, 0., 0, &infoFile);
- TGeoVolume* module_shield10 = ConstructShield("shield10", smallModuleSize, 0., 0, &infoFile);
-
- // TGeoVolume* module2060_hole = ConstructModuleWithHole("module2060_hole", bigModuleSize, 60, &infoFile);
- // --------------------------------------------------------------------------
-
-
- // --------------- Create PSD volume ------------------------------------
- cout << endl;
- cout << "===> Creating PSD...." << endl;
-
- // --- Determine size of PSD box (half-lengths)
- const Double_t psdSizeX = 0.5 * nModulesX * bigModuleSize;
- const Double_t psdSizeY = 0.5 * nModulesY * bigModuleSize;
- TGeoBBox* shape = dynamic_cast<TGeoBBox*>(module2060->GetShape());
- const Double_t psdSizeZ = shape->GetDZ();
-
- TString psdName = "psd_";
- psdName += geoTag;
- TGeoVolume* psd = gGeoManager->MakeBox(psdName, psdMedium, psdSizeX, psdSizeY, psdSizeZ + 0.5 * shieldWidth);
- cout << "Module array is " << nModulesX << " x " << nModulesY << endl;
- cout << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm" << endl;
- infoFile << endl
- << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ + shieldWidth
- << " cm" << endl;
- // --------------------------------------------------------------------------
-
-
- // ----- Place modules into the system volume ---------------------------
-
- Int_t nModules {0};
- Int_t iModule {0};
-
- const Int_t nSmallModulesX = 4;
- const Int_t nSmallModulesY = 4;
- Double_t xModCurrent = -0.5 * (nSmallModulesX + 1) * smallModuleSize;
- for (Int_t iModX = 0; iModX < nSmallModulesX; iModX++) {
- xModCurrent += smallModuleSize;
- Double_t yModCurrent = -0.5 * (nSmallModulesY + 1) * smallModuleSize;
- for (Int_t iModY = 0; iModY < nSmallModulesY; iModY++) {
- yModCurrent += smallModuleSize;
-
- // remove 4 central modules
- if ((iModX == nSmallModulesX / 2 || iModX == nSmallModulesX / 2 - 1)
- && (iModY == nSmallModulesY / 2 || iModY == nSmallModulesY / 2 - 1))
- continue;
-
- // Position of module inside PSD
- TGeoTranslation* trans = new TGeoTranslation(xModCurrent, yModCurrent, -0.5 * shieldWidth);
- TGeoTranslation* trans_shield = new TGeoTranslation(xModCurrent, yModCurrent, psdSizeZ);
-
- psd->AddNode(module1060, iModule, trans);
- psd->AddNode(module_shield10, iModule, trans_shield);
-
- cout << "Adding module " << setw(5) << right << iModule << " at x = " << setw(6) << right << xModCurrent
- << " cm , y = " << setw(6) << right << yModCurrent << " cm" << endl;
- iModule++;
- nModules++;
-
- } //# modules in y direction
- } //# modules in x direction
-
-
- xModCurrent = -0.5 * (nModulesX + 1) * bigModuleSize;
- for (Int_t iModX = 0; iModX < nModulesX; iModX++) {
- xModCurrent += bigModuleSize;
- Double_t yModCurrent = -0.5 * (nModulesY + 1) * bigModuleSize;
- for (Int_t iModY = 0; iModY < nModulesY; iModY++) {
- yModCurrent += bigModuleSize;
-
- // Skip edge modules
- if ((iModX == 0 || iModX == nModulesX - 1) && (iModY == 0 || iModY == nModulesY - 1)) continue;
-
- // Replace 4 central modules with 12 (10 cm x 10 cm) and a (20 cm x 20 cm) hole
- if ((iModX == nModulesX / 2 || iModX == nModulesX / 2 - 1)
- && (iModY == nModulesY / 2 || iModY == nModulesY / 2 - 1))
- continue;
-
- // Position of module inside PSD
- TGeoTranslation* trans = new TGeoTranslation(xModCurrent, yModCurrent, -0.5 * shieldWidth);
- TGeoTranslation* trans_shield = new TGeoTranslation(xModCurrent, yModCurrent, psdSizeZ);
-
- psd->AddNode(module2060, iModule, trans);
- psd->AddNode(module_shield20, iModule, trans_shield);
-
- cout << "Adding module " << setw(5) << right << iModule << " at x = " << setw(6) << right << xModCurrent
- << " cm , y = " << setw(6) << right << yModCurrent << " cm" << endl;
- iModule++;
- nModules++;
-
- } //# modules in y direction
- } //# modules in x direction
-
-
- cout << "PSD contains " << nModules << " modules." << endl;
- infoFile << "PSD contains " << nModules << " modules." << endl;
- // --------------------------------------------------------------------------
-
- // ----- Place PSD in top node (cave) -------------------------------------
- Double_t psdHalfLength = psdSizeZ + 0.5 * shieldWidth;
- Double_t psdVolCenterX = psdX + psdHalfLength * sin(psdRotY);
- Double_t psdVolCenterY = psdY;
- Double_t psdVolCenterZ = psdZ + psdHalfLength * cos(psdRotY);
- infoFile << "PSD volume center coordinates: (" << psdVolCenterX << ", " << psdVolCenterY << ", " << psdVolCenterZ
- << ") cm" << endl;
- TGeoRotation* psdRot = new TGeoRotation();
- psdRot->RotateY(psdRotY * TMath::RadToDeg());
- TGeoCombiTrans* psdTrans = new TGeoCombiTrans(psdVolCenterX, psdVolCenterY, psdVolCenterZ, psdRot);
- top->AddNode(psd, 0, psdTrans);
- cout << endl << "==> PSD position in cave: " << endl;
- psdTrans->Print();
- // --------------------------------------------------------------------------
-
-
- // ----- Close and check geometry ---------------------------------------
- cout << endl << "==> Closing geometry..." << endl;
- gGeoManager->CloseGeometry();
- gGeoManager->CheckGeometryFull();
- cout << endl;
- gGeoManager->CheckOverlaps(0.0001, "s");
- // --------------------------------------------------------------------------
-
-
- // ----- Write PSD volume and placement matrix to geometry file ---------
- cout << endl;
- TString geoFileName = "psd_";
- geoFileName = geoFileName + geoTag + ".geo.root";
- psd->Export(geoFileName);
- TFile* geoFile = new TFile(geoFileName, "UPDATE");
- psdTrans->Write();
- cout << endl;
- cout << "==> Geometry " << psd->GetName() << " written to " << geoFileName << endl;
- cout << "==> Info written to " << infoFileName << endl;
- geoFile->Close();
- infoFile.close();
- // --------------------------------------------------------------------------
-
-
- // ----- Write entire TGeoManager to file -------------------------------
- TString geoManFileName = "psd_";
- geoManFileName = geoManFileName + geoTag + ".geoman.root";
- TFile* geoManFile = new TFile(geoManFileName, "RECREATE");
- gGeoManager->Write();
- geoManFile->Close();
- cout << "==> TGeoManager " << gGeoManager->GetName() << " written to " << geoManFileName << endl << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ---- Display geometry ------------------------------------------------
- // gGeoManager->SetVisLevel(5); // Scintillator level
- top->Draw("ogl");
- // --------------------------------------------------------------------------
- // top->Raytrace();
-}
-// End of main function
-/** ======================================================================= **/
-
-TGeoVolume* ConstructShield(const char* name, Double_t sizeXY, Double_t holesize, Int_t hole_pos, fstream* info)
-{
-
- const TString suffix = Form("_shield_%d_%d_%d", int(sizeXY), int(holesize), int(hole_pos));
-
- const Double_t psd2iron = 2.;
- const Double_t ironD = 4.;
- const Double_t iron2poly = 2.;
- const Double_t polyD = 8.;
- const Double_t shieldLength = psd2iron + ironD + iron2poly + polyD;
-
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
-
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
-
- TGeoMedium* medShieldPoly = gGeoManager->GetMedium("PsdPolyethylene");
- if (!medShieldPoly) Fatal("Main", "Medium PsdPolyethylene not found");
-
- TGeoVolume* shield = gGeoManager->MakeBox(name, medAir, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * shieldLength);
- TGeoVolume* iron {nullptr};
- TGeoVolume* poly {nullptr};
-
- if (holesize > 0) {
- Int_t sx, sy;
- switch (hole_pos) {
- case 0:
- sx = 1;
- sy = 1;
- break; // defines hole position in x-y plane
- case 1:
- sx = 1;
- sy = -1;
- break;
- case 2:
- sx = -1;
- sy = 1;
- break;
- case 3:
- sx = -1;
- sy = -1;
- break;
- }
- TGeoCombiTrans* rot45 = new TGeoCombiTrans(Form("rot45%s", suffix.Data()), sx * sizeXY / 2., sy * sizeXY / 2., 0.,
- new TGeoRotation("rot", 45., 0., 0.));
- rot45->RegisterYourself();
-
- TGeoBBox* iron1 = new TGeoBBox(Form("iron1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * ironD);
- TGeoBBox* iron2 = new TGeoBBox(Form("iron2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * ironD); //hole
-
- TGeoCompositeShape* ironShape =
- new TGeoCompositeShape(Form("iron1%s-iron2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- iron = new TGeoVolume(Form("iron%s", suffix.Data()), ironShape, medIron);
-
- TGeoBBox* poly1 = new TGeoBBox(Form("poly1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * polyD);
- TGeoBBox* poly2 = new TGeoBBox(Form("poly2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * polyD); //hole
-
- TGeoCompositeShape* polyShape =
- new TGeoCompositeShape(Form("poly1%s-poly2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- poly = new TGeoVolume(Form("poly%s", suffix.Data()), polyShape, medShieldPoly);
- }
- else {
- iron = gGeoManager->MakeBox("shield_iron", medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * ironD);
- poly = gGeoManager->MakeBox("shield_poly", medShieldPoly, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * polyD);
- }
-
- shield->AddNode(iron, 0, new TGeoTranslation(0., 0., -shieldLength / 2. + psd2iron + ironD / 2.));
- shield->AddNode(poly, 0, new TGeoTranslation(0., 0., -shieldLength / 2. + psd2iron + ironD + iron2poly + polyD / 2.));
-
- return shield;
-}
-
-/** ======================================================================= **/
-TGeoVolume* ConstructModule(const char* name, Double_t sizeXY, Int_t nLayers, fstream* info)
-{
-
- // The module consists of nLayers of scintillators covered with Tyvek.
- // Between each two scintillators, there is a lead layer (total nLayers -1).
- // At the top, there is a gap to stretch the light fibres. This is modelled
- // by a volume made of the same material as the scintillator (but inactive).
- // The arrangement is surrounded by iron. At the front and at the back,
- // there is a thick iron layer (the back one acting as first absorber)
- // for constructional reasons.
-
- cout << "===> Creating Module " << name << ", size " << sizeXY << " cm with " << nLayers << " layers...." << endl;
-
- const TString suffix = Form("_%d_%d", int(sizeXY), nLayers); // unique suffix for volume names
-
- // ----- Constructional parameters ----------------------------
- const Double_t leadD = 1.60; // Thickness of lead layer
- const Double_t scintD = 0.40; // Thickness of scintillator layer
- const Double_t tyvekD = 0.02; // Thickness of Tyvek wrapper around scintillator
- const Double_t boxDx = 0.15; // Thickness of iron box lateral
- const Double_t boxDy = 0.05; // Thickness of iron box top/bottom
- const Double_t boxDz = 2.00; // Thickness of iron box front/back
- const Double_t chanDy = 0.20; // Height of fibre channel
- const Double_t chanDistL = 0.50; // Distance of channel from left edge of lead
- const Double_t chanDistR = 2.00; // Distance of channel from right edge of lead
- // ------------------------------------------------------------------------
-
-
- // ----- Info to file -------------------------------------------------
- if (info) {
- (*info) << "Parameters of module " << name << ": " << endl;
- (*info) << "Size: " << sizeXY << " cm x " << sizeXY << " cm" << endl;
- (*info) << "Number of layers: " << nLayers << endl;
- (*info) << "Thickness of lead layers: " << leadD << " cm" << endl;
- (*info) << "Thickness of scintillators: " << scintD << " cm" << endl;
- (*info) << "Thickness of Tyvek wrap: " << tyvekD << " cm" << endl;
- (*info) << "Thickness of iron box: (" << boxDx << " / " << boxDy << " / " << boxDz << ") cm" << endl;
- (*info) << "Height of fibre channel: " << chanDy << " cm" << endl;
- (*info) << "Distance of channel from edges: left " << chanDistL << " cm, right " << chanDistR << " cm" << endl;
- }
- // ------------------------------------------------------------------------
-
-
- // ----- Get required media -------------------------------------------
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
- TGeoMedium* medLead = gGeoManager->GetMedium("lead"); // Lead layers
- if (!medLead) Fatal("ConstructModule", "Medium lead not found");
- TGeoMedium* medTyvek = gGeoManager->GetMedium("PsdTyvek"); // Tyvek layers
- if (!medTyvek) Fatal("ConstructModule", "Medium PsdTyvek not found");
- TGeoMedium* medScint = gGeoManager->GetMedium("PsdScint"); // Scintillator
- if (!medScint) Fatal("ConstructModule", "Medium PsdScint not found");
- TGeoMedium* medFibre = gGeoManager->GetMedium("PsdFibre"); // Fibres
- if (!medFibre) Fatal("ConstructModule", "Medium PsdFibre not found");
- // ------------------------------------------------------------------------
-
-
- // ----- Create module volume -----------------------------------------
- // Module length: nLayers of scintillators, nLayers - 1 of lead
- // plus the iron box front and back.
- // Material is iron.
- const Double_t moduleLength = 2. * boxDz + nLayers * (scintD + 2. * tyvekD) + (nLayers - 1) * leadD;
- TGeoVolume* module = gGeoManager->MakeBox(name, medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * moduleLength);
- module->SetLineColor(kRed);
- module->Print();
- cout << endl;
- // ------------------------------------------------------------------------
-
-
- // ----- Lead filler --------------------------------------------------
- // The lead filler represents all lead absorber layers.
- // The Tyvek/scintillator layers and the fibre channel will be placed
- // inside.
- // Dimensions are half-lengths.
- const Double_t leadLx = 0.5 * sizeXY - boxDx;
- const Double_t leadLy = 0.5 * sizeXY - boxDy;
- const Double_t leadLz = 0.5 * moduleLength - boxDz;
- TGeoVolume* lead = gGeoManager->MakeBox(Form("lead%s", suffix.Data()), medLead, leadLx, leadLy, leadLz);
- lead->SetLineColor(kBlue);
- // ------------------------------------------------------------------------
-
-
- // ----- Fibre channel ------------------------------------------------
- // This volume represents the air gap at the top of the volume hosting the
- // light guide fibres. The material is the same as for the scintillators
- // (but inactive).
- // The (half-)width of the channel is defined by the distances from the
- // lead edges.
- const Double_t chanLx = leadLx - 0.5 * chanDistL - 0.5 * chanDistR;
- if (chanLx < 0.) {
- cout << "Error: lead volume is not large enough to host fibre channel!" << endl;
- cout << "Lead width: " << 2. * leadLx << ", distance from left edge " << chanDistL << ", distance from right edge "
- << chanDistR << endl;
- return 0;
- }
- const Double_t chanLy = 0.5 * chanDy;
- const Double_t chanLz = leadLz;
- TGeoVolume* channel = gGeoManager->MakeBox(Form("channel%s", suffix.Data()), medFibre, chanLx, chanLy, chanLz);
- channel->SetLineColor(kMagenta);
- // ------------------------------------------------------------------------
-
-
- // ---- Positioning of the fibre channel ------------------------------
- // The channel extends from chanDistL from the left edge of the lead filler
- // to chanDistR from its right edge (seen from the front). It is top-aligned
- // with the lead filler.
- const Double_t chanShiftX = 0.5 * (chanDistL - chanDistR);
- const Double_t chanShiftY = leadLy - 0.5 * chanDy;
- // ------------------------------------------------------------------------
- // cout << " leadLx = " << leadLx << " leadLy = " << leadLy << endl;
- // cout << " chanShiftX = " << chanShiftX << " chanShiftY = " << chanShiftY << endl;
-
- // ----- Tyvek layer --------------------------------------------------
- // The scintillator will be placed inside this, leaving only the thin
- // Tyvek as a wrapper. Since these layers will be placed in the lead filler,
- // there has to be a cut-out for the fibre channel.
- const Double_t tyvekLz = 0.5 * scintD + tyvekD; // half-length
- TGeoBBox* tyv1 = new TGeoBBox(Form("psdTyv1%s", suffix.Data()), leadLx, leadLy, tyvekLz);
- TGeoBBox* tyv2 = new TGeoBBox(Form("psdTyv2%s", suffix.Data()), chanLx, chanLy, tyvekLz);
- TGeoTranslation* trans1 = new TGeoTranslation(Form("tPsd1%s", suffix.Data()), chanShiftX, chanShiftY, 0.);
-
- trans1->RegisterYourself();
- TGeoCompositeShape* tyvekShape =
- new TGeoCompositeShape(Form("psdTyv1%s-psdTyv2%s:tPsd1%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* tyvek = new TGeoVolume(Form("tyvek%s", suffix.Data()), tyvekShape, medTyvek);
- tyvek->SetLineColor(kGreen);
- // ------------------------------------------------------------------------
-
-
- // ----- Scintillator layer -------------------------------------------
- // The scintillator layer is embedded (as daughter) into the Tyvek layer.
- // It is also a box minus the cut-out for the fibres. The cut-out is
- // slightly smaller than for the Tyvek volume.
- const Double_t scintLx = leadLx - tyvekD;
- const Double_t scintLy = leadLy - tyvekD;
- const Double_t scintLz = 0.5 * scintD;
- TGeoBBox* sci1 = new TGeoBBox(Form("sci1%s", suffix.Data()), scintLx, scintLy, scintLz);
- const Double_t cutLx = chanLx;
- const Double_t cutLy = chanLy - tyvekD;
- const Double_t cutLz = scintLz;
- TGeoBBox* sci2 = new TGeoBBox(Form("sci2%s", suffix.Data()), cutLx, cutLy, cutLz);
- const Double_t cutShiftX = chanShiftX;
- const Double_t cutShiftY = scintLy - cutLy;
- TGeoTranslation* trans2 = new TGeoTranslation(Form("tPsd2%s", suffix.Data()), cutShiftX, cutShiftY, 0.);
- trans2->RegisterYourself();
- TGeoCompositeShape* scintShape = new TGeoCompositeShape(
- Form("scintShape%s", suffix.Data()), Form("sci1%s-sci2%s:tPsd2%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* scint = new TGeoVolume(Form("scint%s", suffix.Data()), scintShape, medScint);
- scint->SetLineColor(kBlack);
-
- // ------------------------------------------------------------------------
- // cout << " scintLx = " << scintLx << " scintLy = " << scintLy << endl;
- // cout << " cutLx = " << cutLx << " cutLy = " << cutLy << endl;
- // cout << " cutShiftX = " << cutShiftX << " cutShiftY = " << cutShiftY << endl;
-
- // ----- Place volumes ------------------------------------------------
-
- // ---> Scintillator into Tyvek
- tyvek->AddNode(scint, 0);
-
- // ---> Fibre channel into lead filler
- lead->AddNode(channel, 0, trans1);
-
- // ---> Tyvek layers into lead
- Double_t zFirst {0.};
- Double_t zLast {0.};
- for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
- const Double_t zPos = -1. * leadLz + iLayer * leadD + (2. * iLayer + 1.) * tyvekLz;
- if (iLayer == 0) zFirst = zPos;
- if (iLayer == nLayers - 1) zLast = zPos;
- lead->AddNode(tyvek, iLayer, new TGeoTranslation(0., 0., zPos));
- }
- cout << module->GetName() << ": Positioned " << nLayers << " Tyvek layers; first at z = " << zFirst
- << ", last at z = " << zLast << endl;
-
- // ---> Lead into module
- module->AddNode(lead, 0);
- // ------------------------------------------------------------------------
-
- return module;
-}
-
-
-TGeoVolume* ConstructModuleWithHole(const char* name, Double_t sizeXY, Int_t nLayers, fstream* info, float holesize,
- Int_t hole_pos)
-{
-
- // The module consists of nLayers of scintillators covered with Tyvek.
- // Between each two scintillators, there is a lead layer (total nLayers -1).
- // At the top, there is a gap to stretch the light fibres. This is modelled
- // by a volume made of the same material as the scintillator (but inactive).
- // The arrangement is surrounded by iron. At the front and at the back,
- // there is a thick iron layer (the back one acting as first absorber)
- // for constructional reasons.
-
-
- if (holesize == 0.) return ConstructModule(name, sizeXY, nLayers, info);
-
- const TString suffix = Form("_%d_%d_%d", int(sizeXY), int(holesize), int(hole_pos));
-
- cout << "===> Creating Module " << name << ", size " << sizeXY << " cm with " << nLayers << " layers...." << endl;
-
- // ----- Constructional parameters ----------------------------
- const Double_t leadD = 1.60; // Thickness of lead layer
- const Double_t scintD = 0.40; // Thickness of scintillator layer
- const Double_t tyvekD = 0.02; // Thickness of Tyvek wrapper around scintillator
-
- const Double_t boxDx = 0.15; // Thickness of iron box lateral
- const Double_t boxDy = 0.05; // Thickness of iron box top/bottom
- const Double_t boxDz = 2.00; // Thickness of iron box front/back
-
- const Double_t chanDy = 0.20; // Height of fibre channel
- const Double_t chanDistL = 0.50; // Distance of channel from left edge of lead
- const Double_t chanDistR = 2.00; // Distance of channel from right edge of lead
-
- const Double_t moduleLength = 2. * boxDz + nLayers * (scintD + 2. * tyvekD) + (nLayers - 1) * leadD;
-
- // Dimensions are half-lengths.
- const Double_t leadLx = 0.5 * sizeXY - boxDx;
- const Double_t leadLy = 0.5 * sizeXY - boxDy;
- const Double_t leadLz = 0.5 * moduleLength - boxDz;
- const Double_t tyvekLz = 0.5 * scintD + tyvekD; // half-length
-
- const Double_t scintLx = leadLx - tyvekD; // Scintillator size x
- const Double_t scintLy = leadLy - tyvekD; // Scintillator size y
- const Double_t scintLz = 0.5 * scintD; // Scintillator size z
-
-
- // ------------------------------------------------------------------------
- Int_t sx, sy;
- switch (hole_pos) {
- case 0:
- sx = 1;
- sy = 1;
- break; // defines hole position in x-y plane
- case 1:
- sx = 1;
- sy = -1;
- break;
- case 2:
- sx = -1;
- sy = 1;
- break;
- case 3:
- sx = -1;
- sy = -1;
- break;
- }
- // ----- Info to file -------------------------------------------------
- if (info) {
- (*info) << "Parameters of module " << name << ": " << endl;
- (*info) << "Size: " << sizeXY << " cm x " << sizeXY << " cm" << endl;
- (*info) << "Number of layers: " << nLayers << endl;
- (*info) << "Thickness of lead layers: " << leadD << " cm" << endl;
- (*info) << "Thickness of scintillators: " << scintD << " cm" << endl;
- (*info) << "Thickness of Tyvek wrap: " << tyvekD << " cm" << endl;
- (*info) << "Thickness of iron box: (" << boxDx << " / " << boxDy << " / " << boxDz << ") cm" << endl;
- (*info) << "Height of fibre channel: " << chanDy << " cm" << endl;
- (*info) << "Distance of channel from edges: left " << chanDistL << " cm, right " << chanDistR << " cm" << endl;
- }
- // ------------------------------------------------------------------------
-
- // ----- Get required media -------------------------------------------
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
- TGeoMedium* medLead = gGeoManager->GetMedium("lead"); // Lead layers
- if (!medLead) Fatal("ConstructModule", "Medium lead not found");
- TGeoMedium* medTyvek = gGeoManager->GetMedium("PsdTyvek"); // Tyvek layers
- if (!medTyvek) Fatal("ConstructModule", "Medium PsdTyvek not found");
- TGeoMedium* medScint = gGeoManager->GetMedium("PsdScint"); // Scintillator
- if (!medScint) Fatal("ConstructModule", "Medium PsdScint not found");
- TGeoMedium* medFibre = gGeoManager->GetMedium("PsdFibre"); // Fibres
- if (!medFibre) Fatal("ConstructModule", "Medium PsdFibre not found");
- // ------------------------------------------------------------------------
-
- TGeoCombiTrans* rot45 = new TGeoCombiTrans(Form("rot45%s", suffix.Data()), sx * sizeXY / 2., sy * sizeXY / 2., 0.,
- new TGeoRotation("rot", 45., 0., 0.));
- rot45->RegisterYourself();
-
- // ----- Create module volume -----------------------------------------
- // Module length: nLayers of scintillators, nLayers - 1 of lead
- // plus the iron box front and back.
- // Material is iron.
- // TGeoVolume* module = gGeoManager->MakeBox(name, medIron, 0.5 * sizeXY,
- // 0.5 * sizeXY, 0.5 * moduleLength);
-
- TGeoBBox* module1 = new TGeoBBox(Form("module1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * moduleLength);
- TGeoBBox* module2 = new TGeoBBox(Form("module2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * moduleLength); //hole
- TGeoCompositeShape* moduleShape =
- new TGeoCompositeShape(Form("module1%s-module2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* module = new TGeoVolume(name, moduleShape, medIron);
-
- module->SetLineColor(kRed);
- module->Print();
- cout << endl;
- // ------------------------------------------------------------------------
-
-
- // ----- Lead filler --------------------------------------------------
- // The lead filler represents all lead absorber layers.
- // The Tyvek/scintillator layers and the fibre channel will be placed
- // inside.
- // TGeoVolume* lead = gGeoManager->MakeBox(Form("lead%s", suffix.Data()), medLead, leadLx, leadLy, leadLz);
-
- TGeoBBox* lead1 = new TGeoBBox(Form("lead1%s", suffix.Data()), leadLx, leadLy, leadLz);
- TGeoBBox* lead2 = new TGeoBBox(Form("lead2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- leadLz); //hole
- TGeoCompositeShape* leadShape =
- new TGeoCompositeShape(Form("lead1%s-lead2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* lead = new TGeoVolume(Form("lead%s", suffix.Data()), leadShape, medLead);
-
- lead->SetLineColor(kBlue);
- // ------------------------------------------------------------------------
-
- // ----- Tyvek layer --------------------------------------------------
- // The scintillator will be placed inside this, leaving only the thin
- // Tyvek as a wrapper. Since these layers will be placed in the lead filler,
- // there has to be a cut-out for the fibre channel.
-
-
- TGeoBBox* tyv1 = new TGeoBBox(Form("psdTyv1%s", suffix.Data()), leadLx, leadLy, tyvekLz);
- TGeoBBox* tyv2 = new TGeoBBox(Form("psdTyv2%s", suffix.Data()), holesize / 2., holesize / 2.,
- tyvekLz); //hole
-
- TGeoCompositeShape* tyvekShape =
- new TGeoCompositeShape(Form("psdTyv1%s-psdTyv2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* tyvek = new TGeoVolume(Form("tyvek%s", suffix.Data()), tyvekShape, medTyvek);
- tyvek->SetLineColor(kGreen);
- // ------------------------------------------------------------------------
-
-
- // ----- Scintillator layer -------------------------------------------
- // The scintillator layer is embedded (as daughter) into the Tyvek layer.
- // It is also a box minus the cut-out for the fibres. The cut-out is
- // slightly smaller than for the Tyvek volume.
- TGeoBBox* sci1 = new TGeoBBox(Form("sci1%s", suffix.Data()), scintLx, scintLy, scintLz);
- TGeoBBox* sci2 = new TGeoBBox(Form("sci2%s", suffix.Data()), holesize / 2., holesize / 2.,
- scintLz); //hole
- TGeoCompositeShape* scintShape = new TGeoCompositeShape(
- Form("scintShape%s", suffix.Data()), Form("sci1%s-sci2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* scint = new TGeoVolume(Form("scint%s", suffix.Data()), scintShape, medScint);
- scint->SetLineColor(kBlack);
-
- // ----- Place volumes ------------------------------------------------
-
- // ---> Scintillator into Tyvek
- tyvek->AddNode(scint, 0);
-
- // ---> Fibre channel into lead filler
- // lead->AddNode(channel, 0, trans1); //Needs addidional size check
-
- // ---> Tyvek layers into lead
- Double_t zFirst {0.};
- Double_t zLast {0.};
- for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
- const Double_t zPos = -1. * leadLz + iLayer * leadD + (2. * iLayer + 1.) * tyvekLz;
- if (iLayer == 0) zFirst = zPos;
- if (iLayer == nLayers - 1) zLast = zPos;
- lead->AddNode(tyvek, iLayer, new TGeoTranslation(0., 0., zPos));
- }
- cout << module->GetName() << ": Positioned " << nLayers << " Tyvek layers; first at z = " << zFirst
- << ", last at z = " << zLast << endl;
-
- // ---> Lead into module
- module->AddNode(lead, 0);
- // ------------------------------------------------------------------------
-
- return module;
-}
diff --git a/macro/psd/create_psdgeo_ideal.C b/macro/psd/create_psdgeo_ideal.C
deleted file mode 100644
index 63eb66745a..0000000000
--- a/macro/psd/create_psdgeo_ideal.C
+++ /dev/null
@@ -1,312 +0,0 @@
-/* Copyright (C) 2017-2018 National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Moscow
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Volker Friese, Evgeny Kashirin [committer] */
-
-/** @file create_psdgeo_with_hole.C
- ** @author Volker Friese <v.friese@gsi.de>
- ** @date 16 September 2017
- ** @version 1.0
- **
- ** This macro creates a PSD geometry in ROOT format to be used as input
- ** for the transport simulation. It allows to create module stacks of
- ** varying sizes; the numbers of rows and columns must both be uneven.
- ** The centre module as well as the edge modules are omitted.
- ** The module definition is implemented in the function ConstructModule.
- **
- ** The user can adjust the steering variables to the specific needs.
- **/
-
-#include <iostream>
-
-
-using std::cout;
-using std::endl;
-using std::fstream;
-using std::right;
-using std::setw;
-
-// Forward declarations
-
-TGeoVolume* ConstructIdealPsd(const char* name, Double_t sizeX, Double_t sizeY, Double_t sizeZ, Double_t modSize,
- Double_t holesize, fstream* info);
-
-// ============================================================================
-// ====== Main function =====
-// ============================================================================
-
-void create_psdgeo_ideal(TString geoTag = "v22a")
-{
-
- // ----- Steering variables ---------------------------------------------
- Double_t psdX; // x position (cm) of PSD in cave (front plane center)
- Double_t psdY; // y position (cm) of PSD in cave (front plane center)
- Double_t psdZ; // z position (cm) of PSD in cave (front plane center)
- Double_t psdRotY; // Rotation of PSD around y axis (rad)
- Double_t holeSize; // side length of the square shaped hole (cm)
- TString comment; // short description
-
- if (geoTag == "v22a") {
- psdX = 13.1;
- psdY = 0.;
- psdZ = 1010;
- psdRotY = 0.01335;
- holeSize = 0.;
- comment = "Position for Au beam at 12A GeV/c and 100% magnetic field strength";
- }
- else if (geoTag == "v22b") {
- psdX = 28.56;
- psdY = 0.;
- psdZ = 1010;
- psdRotY = 0.02906;
- holeSize = 0.;
- comment = "Position for Au beam at 3.3A GeV/c and 60% magnetic field strength (to fit beam dump)";
- }
- else if (geoTag == "v20c") {
- psdX = 12.95;
- psdY = 0.;
- psdZ = 1010;
- psdRotY = 0.0132;
- holeSize = 20.;
- }
- else if (geoTag == "v20a") {
- psdX = 12.95;
- psdY = 0.;
- psdZ = 1050;
- psdRotY = 0.0132;
- holeSize = 20.;
- }
- else if (geoTag == "v20b") {
- psdX = 13.75;
- psdY = 0.;
- psdZ = 1100;
- psdRotY = 0.0132;
- holeSize = 20.;
- }
- else if (geoTag == "v18e") {
- psdX = 9.65;
- psdY = 0.;
- psdZ = 800.;
- psdRotY = 0.01321;
- holeSize = 20.;
- }
- else if (geoTag == "v18f") {
- psdX = 9.65;
- psdY = 0.;
- psdZ = 800.;
- psdRotY = 0.01321;
- holeSize = 6.;
- }
- else if (geoTag == "v18g") {
- psdX = 9.65;
- psdY = 0.;
- psdZ = 800.;
- psdRotY = 0.01321;
- holeSize = 0.;
- }
- else if (geoTag == "v18m") {
- psdX = 8.9;
- psdY = 0.;
- psdZ = 1500.;
- psdRotY = 0.01321;
- holeSize = 20.;
- }
- else {
- cout << "\n\n\nERROR: Unrecognized geoTag! Exiting!\n\n\n";
- return;
- }
-
- geoTag += "_ideal";
- const Double_t bigModuleSize = 20.; // Module size [cm]
- const Int_t nModulesX = 8; // Number of modules in a row (x direction)
- const Int_t nModulesY = 6; // Number of modules in a row (x direction)
-
- // --------------------------------------------------------------------------
- cout << "Number of modules per row and column = " << nModulesX << " " << nModulesY << endl;
-
- // ------- Open info file -----------------------------------------------
- TString infoFileName = "psd_";
- infoFileName = infoFileName + geoTag + ".geo.info";
- fstream infoFile;
- infoFile.open(infoFileName.Data(), fstream::out);
- infoFile << "PSD geometry " << geoTag << " created with create_psdgeo_ideal.C" << endl << endl;
- infoFile << "Ideal PSD geometry for performance studies - whole volume as a "
- "single module"
- << endl;
- infoFile << comment << endl << endl;
- infoFile << "Number of modules: " << nModulesX << " x " << nModulesY << endl;
- infoFile << "Big module size: " << bigModuleSize << " cm x " << bigModuleSize << " cm" << endl;
- infoFile << "PSD front plane center coordinates: (" << psdX << ", " << psdY << ", " << psdZ << ") cm" << endl;
- infoFile << "PSD rotation around y axis: " << psdRotY * TMath::RadToDeg() << " degrees" << endl;
- infoFile << "Side length of the square shaped hole in PSD center: " << holeSize << " cm" << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ------- Load media from media file -----------------------------------
- cout << endl << "==> Reading media..." << endl;
- FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- TString geoPath = gSystem->Getenv("VMCWORKDIR");
- TString medFile = geoPath + "/geometry/media.geo";
- geoFace->setMediaFile(medFile);
- geoFace->readMedia();
- // --------------------------------------------------------------------------
-
-
- // ------ Create the required media from the media file ----------------
- FairGeoMedia* geoMedia = geoFace->getMedia();
- FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
-
- // ---> Air
- FairGeoMedium* mAir = geoMedia->getMedium("air");
- if (!mAir) Fatal("Main", "FairMedium air not found");
- geoBuild->createMedium(mAir);
- TGeoMedium* air = gGeoManager->GetMedium("air");
- if (!air) Fatal("Main", "Medium air not found");
- cout << "Created medium: air" << endl;
-
- // ---> Vacuum
- FairGeoMedium* mVacuum = geoMedia->getMedium("vacuum");
- if (!mVacuum) Fatal("Main", "FairMedium vacuum not found");
- geoBuild->createMedium(mVacuum);
- TGeoMedium* vacuum = gGeoManager->GetMedium("vacuum");
- if (!vacuum) Fatal("Main", "Medium vacuum not found");
- cout << "Created medium: vacuum" << endl;
-
- // ---> Medium for PSD volume
- TGeoMedium* psdMedium = vacuum;
- // --------------------------------------------------------------------------
-
- // -------------- Create geometry and top volume -------------------------
- cout << endl << "==> Set top volume..." << endl;
- gGeoManager->SetName("PSDgeom");
- TGeoVolume* top = new TGeoVolumeAssembly("TOP");
- gGeoManager->SetTopVolume(top);
- // --------------------------------------------------------------------------
-
-
- // --------------- Create PSD volume ------------------------------------
- cout << endl;
- cout << "===> Creating PSD...." << endl;
-
- // --- Determine size of PSD box (half-lengths)
- const Double_t psdSizeX = 0.5 * nModulesX * bigModuleSize;
- const Double_t psdSizeY = 0.5 * nModulesY * bigModuleSize;
- const Double_t psdSizeZ = 0.5 * 120.;
-
- TString psdName = "psd_";
- psdName += geoTag;
- TGeoVolume* psd = gGeoManager->MakeBox(psdName, psdMedium, psdSizeX, psdSizeY, psdSizeZ);
- cout << "Module array is " << nModulesX << " x " << nModulesY << endl;
- cout << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm" << endl;
- infoFile << endl
- << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm"
- << endl;
- infoFile << "PSD front plane center coordinates: (" << psdX << ", " << psdY << ", " << psdZ << ") cm" << endl;
- infoFile << "PSD rotation around y axis: " << psdRotY * TMath::RadToDeg() << " degrees" << endl;
- infoFile << "Diagonal size of the square shaped hole in PSD center: " << holeSize << " cm" << endl << endl;
- // --------------------------------------------------------------------------
-
- psd->AddNode(
- ConstructIdealPsd("ideal_psd", psdSizeX, psdSizeY, psdSizeZ, bigModuleSize / 2., holeSize / 2., &infoFile), 0);
-
- // --------------------------------------------------------------------------
-
-
- // ----- Place PSD in top node (cave) -------------------------------------
- Double_t psdHalfLength = psdSizeZ;
- Double_t psdVolCenterX = psdX + psdHalfLength * sin(psdRotY);
- Double_t psdVolCenterY = psdY;
- Double_t psdVolCenterZ = psdZ + psdHalfLength * cos(psdRotY);
- infoFile << "PSD volume center coordinates: (" << psdVolCenterX << ", " << psdVolCenterY << ", " << psdVolCenterZ
- << ") cm" << endl;
- TGeoRotation* psdRot = new TGeoRotation();
- psdRot->RotateY(psdRotY * TMath::RadToDeg());
- TGeoCombiTrans* psdTrans = new TGeoCombiTrans(psdVolCenterX, psdVolCenterY, psdVolCenterZ, psdRot);
- top->AddNode(psd, 0, psdTrans);
- cout << endl << "==> PSD position in cave: " << endl;
- psdTrans->Print();
- // --------------------------------------------------------------------------
-
-
- // ----- Close and check geometry ---------------------------------------
- cout << endl << "==> Closing geometry..." << endl;
- gGeoManager->CloseGeometry();
- gGeoManager->CheckGeometryFull();
- cout << endl;
- gGeoManager->CheckOverlaps(0.0001, "s");
- // --------------------------------------------------------------------------
-
-
- // ----- Write PSD volume and placement matrix to geometry file ---------
- cout << endl;
- TString geoFileName = "psd_";
- geoFileName = geoFileName + geoTag + ".geo.root";
- psd->Export(geoFileName);
- TFile* geoFile = new TFile(geoFileName, "UPDATE");
- psdTrans->Write();
- cout << endl;
- cout << "==> Geometry " << psd->GetName() << " written to " << geoFileName << endl;
- cout << "==> Info written to " << infoFileName << endl;
- geoFile->Close();
- infoFile.close();
- // --------------------------------------------------------------------------
-
-
- // ----- Write entire TGeoManager to file -------------------------------
- TString geoManFileName = "psd_";
- geoManFileName = geoManFileName + geoTag + ".geoman.root";
- TFile* geoManFile = new TFile(geoManFileName, "RECREATE");
- gGeoManager->Write();
- geoManFile->Close();
- cout << "==> TGeoManager " << gGeoManager->GetName() << " written to " << geoManFileName << endl << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ---- Display geometry ------------------------------------------------
- // gGeoManager->SetVisLevel(5); // Scintillator level
- top->Draw("ogl");
- // --------------------------------------------------------------------------
- // top->Raytrace();
-}
-// End of main function
-/** ======================================================================= **/
-
-TGeoVolume* ConstructIdealPsd(const char* name, Double_t sizeX, Double_t sizeY, Double_t sizeZ, Double_t modSize,
- Double_t holesize, fstream* info)
-{
-
- TGeoMedium* medVacuum = gGeoManager->GetMedium("vacuum"); // PSD
- if (!medVacuum) Fatal("ConstructModule", "Medium vacuum not found");
-
- TGeoVolume* psd = gGeoManager->MakeBox(name, medVacuum, sizeX, sizeY, sizeZ);
- TGeoVolume* scint {nullptr};
-
- TGeoCombiTrans* rot45 = new TGeoCombiTrans("rot45", 0., 0., 0., new TGeoRotation("rot", 45., 0., 0.));
- rot45->RegisterYourself();
-
- TGeoBBox* moduleBox = new TGeoBBox("ideal_scint", sizeX, sizeY, sizeZ);
-
- TGeoBBox* Hole = new TGeoBBox("hole", holesize, holesize, sizeZ); //hole
- TGeoBBox* Module = new TGeoBBox("module", modSize, modSize, sizeZ); //hole
-
- TGeoTranslation* trans1 = new TGeoTranslation("transl1", (sizeX - modSize), (sizeY - modSize), 0.);
- TGeoTranslation* trans2 = new TGeoTranslation("transl2", (sizeX - modSize), -(sizeY - modSize), 0.);
- TGeoTranslation* trans3 = new TGeoTranslation("transl3", -(sizeX - modSize), (sizeY - modSize), 0.);
- TGeoTranslation* trans4 = new TGeoTranslation("transl4", -(sizeX - modSize), -(sizeY - modSize), 0.);
-
- trans1->RegisterYourself();
- trans2->RegisterYourself();
- trans3->RegisterYourself();
- trans4->RegisterYourself();
-
- TString shape = "ideal_scint-(module:transl1)-(module:transl2)-(module:"
- "transl3)-(module:transl4)";
- if (holesize > 0) shape += "-(hole:rot45)";
-
- TGeoCompositeShape* modShape = new TGeoCompositeShape("sens_area", shape);
- scint = new TGeoVolume("ideal_scint", modShape, medVacuum);
-
- psd->AddNode(scint, 0);
-
- return psd;
-}
diff --git a/macro/psd/create_psdgeo_with_hole.C b/macro/psd/create_psdgeo_with_hole.C
deleted file mode 100644
index 41f70c4438..0000000000
--- a/macro/psd/create_psdgeo_with_hole.C
+++ /dev/null
@@ -1,838 +0,0 @@
-/* Copyright (C) 2017-2020 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tübingen
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Volker Friese, Viktor Klochkov [committer] */
-
-/** @file create_psdgeo_with_hole.C
- ** @author Volker Friese <v.friese@gsi.de>
- ** @date 16 September 2017
- ** @version 1.0
- **
- ** This macro creates a PSD geometry in ROOT format to be used as input
- ** for the transport simulation. It allows to create module stacks of
- ** varying sizes; the numbers of rows and columns must both be uneven.
- ** The centre module as well as the edge modules are omitted.
- ** The module definition is implemented in the function ConstructModule.
- **
- ** The user can adjust the steering variables to the specific needs.
- **/
-
-#include <iostream>
-
-
-using std::cout;
-using std::endl;
-using std::fstream;
-using std::right;
-using std::setw;
-
-// Forward declarations
-TGeoVolume* ConstructModule(const char* name, Double_t moduleSize, Int_t nLayers, fstream* info = 0);
-
-TGeoVolume* ConstructModuleWithHole(const char* name, Double_t moduleSize, Int_t nLayers, fstream* info = 0,
- float holesize = 20, Int_t hole_pos = 0);
-
-TGeoVolume* ConstructShield(const char* name, Double_t sizeXY, Double_t holesize, Int_t hole_pos, fstream* info);
-
-// ============================================================================
-// ====== Main function =====
-// ============================================================================
-
-void create_psdgeo_with_hole(TString geoTag = "v22c")
-{
-
- // ----- Steering variables ---------------------------------------------
- Double_t psdX; // x position (cm) of PSD in cave (front plane center)
- Double_t psdY; // y position (cm) of PSD in cave (front plane center)
- Double_t psdZ; // z position (cm) of PSD in cave (front plane center)
- Double_t psdRotY; // Rotation of PSD around y axis (rad)
- Double_t holeSize; // side length of the square shaped hole (cm)
- TString comment; // short description
-
- if (geoTag == "v22a") {
- psdX = 13.1;
- psdY = 0.;
- psdZ = 1010;
- psdRotY = 0.01335;
- holeSize = 20.;
- comment = "Position for Au beam at 12A GeV/c and 100% magnetic field strength";
- }
- else if (geoTag == "v22b") {
- psdX = 28.56;
- psdY = 0.;
- psdZ = 1010;
- psdRotY = 0.02906;
- holeSize = 20.;
- comment = "Position for Au beam at 3.3A GeV/c and 60% magnetic field strength (to fit beam dump)";
- }
- else if (geoTag == "v22c") {
- psdX = -90;
- psdY = -70;
- psdZ = 1756;
- psdRotY = 0.0;
- holeSize = 20.;
- comment = "This is the PSD parking position for an out-of-the-way position to be used in hadron and muon setuos.";
- }
- else if (geoTag == "v20c") {
- psdX = 12.95;
- psdY = 0.;
- psdZ = 1010;
- psdRotY = 0.0132;
- holeSize = 20.;
- comment = "This is a coordinate shifted version of v20a";
- }
- else if (geoTag == "v20a") {
- psdX = 12.95;
- psdY = 0.;
- psdZ = 1050;
- psdRotY = 0.0132;
- holeSize = 20.;
- }
- else if (geoTag == "v20b") {
- psdX = 13.75;
- psdY = 0.;
- psdZ = 1100;
- psdRotY = 0.0132;
- holeSize = 20.;
- }
- else if (geoTag == "v18e") {
- psdX = 9.65;
- psdY = 0.;
- psdZ = 800.;
- psdRotY = 0.01321;
- holeSize = 20.;
- }
- else if (geoTag == "v18f") {
- psdX = 9.65;
- psdY = 0.;
- psdZ = 800.;
- psdRotY = 0.01321;
- holeSize = 6.;
- }
- else if (geoTag == "v18g") {
- psdX = 9.65;
- psdY = 0.;
- psdZ = 800.;
- psdRotY = 0.01321;
- holeSize = 0.;
- }
- else if (geoTag == "v18m") {
- psdX = 8.9;
- psdY = 0.;
- psdZ = 1500.;
- psdRotY = 0.01321;
- holeSize = 20.;
- }
- else {
- cout << "\n\n\nERROR: Unrecognized geoTag! Exiting!\n\n\n";
- return;
- }
-
- const Double_t bigModuleSize = 20.; // Module size (cm)
- const Int_t nModulesX = 8; // Number of modules in a row (x direction)
- const Int_t nModulesY = 6; // Number of modules in a row (x direction)
- const Int_t nLayers = 60; // Number of sections in a module (z direction)
- const Double_t shieldWidth = 16.;
-
- // --------------------------------------------------------------------------
-
- // ---- Number of modules per row/column must be uneven ---------------------
- // Otherwise, the central one (to be skipped) is not defined.
- // if ( nModulesX % 2 != 1 || nModulesY % 2 != 1 ) {
- cout << "Number of modules per row and column = " << nModulesX << " " << nModulesY << endl;
- // return;
- // }
- // --------------------------------------------------------------------------
-
-
- // ------- Open info file -----------------------------------------------
- TString infoFileName = "psd_";
- infoFileName = infoFileName + geoTag + ".geo.info";
- fstream infoFile;
- infoFile.open(infoFileName.Data(), fstream::out);
- infoFile << "PSD geometry " << geoTag << " created with create_psdgeo_with_hole.C" << endl << endl;
- infoFile << comment << endl << endl;
- infoFile << "Number of modules: " << nModulesX << " x " << nModulesY << endl;
- infoFile << "Big module size: " << bigModuleSize << " cm x " << bigModuleSize << " cm" << endl;
- infoFile << "PSD front plane center coordinates: (" << psdX << ", " << psdY << ", " << psdZ << ") cm" << endl;
- infoFile << "PSD rotation around y axis: " << psdRotY << " rad" << endl;
- infoFile << "Side length of the square shaped hole in PSD center: " << holeSize << " cm" << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ------- Load media from media file -----------------------------------
- cout << endl << "==> Reading media..." << endl;
- FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- TString geoPath = gSystem->Getenv("VMCWORKDIR");
- TString medFile = geoPath + "/geometry/media.geo";
- geoFace->setMediaFile(medFile);
- geoFace->readMedia();
- // --------------------------------------------------------------------------
-
-
- // ------ Create the required media from the media file ----------------
- FairGeoMedia* geoMedia = geoFace->getMedia();
- FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
-
- // ---> Air
- FairGeoMedium* mAir = geoMedia->getMedium("air");
- if (!mAir) Fatal("Main", "FairMedium air not found");
- geoBuild->createMedium(mAir);
- TGeoMedium* air = gGeoManager->GetMedium("air");
- if (!air) Fatal("Main", "Medium air not found");
- cout << "Created medium: air" << endl;
-
- // ---> Iron
- FairGeoMedium* mIron = geoMedia->getMedium("iron");
- if (!mIron) Fatal("Main", "FairMedium iron not found");
- geoBuild->createMedium(mIron);
- TGeoMedium* iron = gGeoManager->GetMedium("iron");
- if (!iron) Fatal("Main", "Medium iron not found");
- cout << "Created medium: iron" << endl;
-
- // ---> Lead
- FairGeoMedium* mLead = geoMedia->getMedium("lead");
- if (!mIron) Fatal("Main", "FairMedium lead not found");
- geoBuild->createMedium(mLead);
- TGeoMedium* lead = gGeoManager->GetMedium("lead");
- if (!lead) Fatal("Main", "Medium iron not found");
- cout << "Created medium: lead" << endl;
-
- // ---> Tyvek
- FairGeoMedium* mTyvek = geoMedia->getMedium("PsdTyvek");
- if (!mTyvek) Fatal("Main", "FairMedium PsdTyvek not found");
- geoBuild->createMedium(mTyvek);
- TGeoMedium* tyvek = gGeoManager->GetMedium("PsdTyvek");
- if (!tyvek) Fatal("Main", "Medium PsdTyvek not found");
- cout << "Created medium: PsdTyvek" << endl;
-
- // ---> Scintillator
- FairGeoMedium* mScint = geoMedia->getMedium("PsdScint");
- if (!mScint) Fatal("Main", "FairMedium PsdScint not found");
- geoBuild->createMedium(mScint);
- TGeoMedium* scint = gGeoManager->GetMedium("PsdScint");
- if (!scint) Fatal("Main", "Medium PsdScint not found");
- cout << "Created medium: PsdScint" << endl;
-
- // ---> Fibres
- FairGeoMedium* mFibre = geoMedia->getMedium("PsdFibre");
- if (!mFibre) Fatal("Main", "FairMedium PsdFibre not found");
- geoBuild->createMedium(mFibre);
- TGeoMedium* fibre = gGeoManager->GetMedium("PsdFibre");
- if (!fibre) Fatal("Main", "Medium PsdFibre not found");
- cout << "Created medium: PsdFibre" << endl;
-
- // ---> Shield Polyethylene
- FairGeoMedium* mShieldPoly = geoMedia->getMedium("PsdPolyethylene");
- if (!mShieldPoly) Fatal("Main", "FairMedium PsdPolyethylene not found");
- geoBuild->createMedium(mShieldPoly);
- TGeoMedium* ShieldPoly = gGeoManager->GetMedium("PsdPolyethylene");
- if (!ShieldPoly) Fatal("Main", "Medium PsdPolyethylene not found");
- cout << "Created medium: PsdPolyethylene" << endl;
-
-
- // ---> Medium for PSD volume
- TGeoMedium* psdMedium = air;
- // --------------------------------------------------------------------------
-
-
- // -------------- Create geometry and top volume -------------------------
- cout << endl << "==> Set top volume..." << endl;
- gGeoManager->SetName("PSDgeom");
- TGeoVolume* top = new TGeoVolumeAssembly("TOP");
- gGeoManager->SetTopVolume(top);
- // --------------------------------------------------------------------------
-
- // ----- Create the PSD modules -----------------------------------------
- cout << endl;
- TGeoVolume* module2060 = ConstructModule("module2060", bigModuleSize, nLayers, &infoFile);
- TGeoVolume* module_shield = ConstructShield("shield20", bigModuleSize, 0., 0, &infoFile);
-
- // TGeoVolume* module2060_hole = ConstructModuleWithHole("module2060_hole", bigModuleSize, 60, &infoFile);
- // --------------------------------------------------------------------------
-
-
- // --------------- Create PSD volume ------------------------------------
- cout << endl;
- cout << "===> Creating PSD...." << endl;
-
- // --- Determine size of PSD box (half-lengths)
- const Double_t psdSizeX = 0.5 * nModulesX * bigModuleSize;
- const Double_t psdSizeY = 0.5 * nModulesY * bigModuleSize;
- TGeoBBox* shape = dynamic_cast<TGeoBBox*>(module2060->GetShape());
- const Double_t psdSizeZ = shape->GetDZ();
-
- TString psdName = "psd_";
- psdName += geoTag;
- TGeoVolume* psd = gGeoManager->MakeBox(psdName, psdMedium, psdSizeX, psdSizeY, psdSizeZ + 0.5 * shieldWidth);
- cout << "Module array is " << nModulesX << " x " << nModulesY << endl;
- cout << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ << " cm" << endl;
- infoFile << endl
- << "PSD size is " << 2. * psdSizeX << " cm x " << 2. * psdSizeY << " cm x " << 2. * psdSizeZ + shieldWidth
- << " cm" << endl;
- // --------------------------------------------------------------------------
-
-
- // ----- Place modules into the system volume ---------------------------
- Double_t xModCurrent = -0.5 * (nModulesX + 1) * bigModuleSize;
-
- Int_t nModules {0};
- Int_t iHole {0};
- Int_t iModule {0};
-
- for (Int_t iModX = 0; iModX < nModulesX; iModX++) {
- xModCurrent += bigModuleSize;
- Double_t yModCurrent = -0.5 * (nModulesY + 1) * bigModuleSize;
- for (Int_t iModY = 0; iModY < nModulesY; iModY++) {
- yModCurrent += bigModuleSize;
-
- // Skip edge modules
- if ((iModX == 0 || iModX == nModulesX - 1) && (iModY == 0 || iModY == nModulesY - 1)) continue;
-
- // Skip centre module (only for uneven number of modules)
- // Replace 4 central modules with 12 (10 cm x 10 cm) and a (20 cm x 20 cm) hole
- // Node number and name (convention)
- Int_t iNode = 100 * iModX + iModY;
-
- // Position of module inside PSD
- TGeoTranslation* trans = new TGeoTranslation(xModCurrent, yModCurrent, -0.5 * shieldWidth);
- TGeoTranslation* trans_shield = new TGeoTranslation(xModCurrent, yModCurrent, psdSizeZ);
-
- // Make hole in 4 central modules
- if ((iModX == nModulesX / 2 || iModX == nModulesX / 2 - 1)
- && (iModY == nModulesY / 2 || iModY == nModulesY / 2 - 1)) {
- psd->AddNode(ConstructModuleWithHole(Form("module2060_hole_%d", iHole), bigModuleSize, nLayers, &infoFile,
- holeSize, iHole),
- iModule, trans);
- psd->AddNode(ConstructShield(Form("shield20_hole_%d", iHole), bigModuleSize, holeSize, iHole, &infoFile),
- iModule, trans_shield);
- iHole++;
- }
- else {
- psd->AddNode(module2060, iModule, trans);
- psd->AddNode(module_shield, iModule, trans_shield);
- }
-
-
- cout << "Adding module " << setw(5) << right << iModule << " at x = " << setw(6) << right << xModCurrent
- << " cm , y = " << setw(6) << right << yModCurrent << " cm" << endl;
- iModule++;
- nModules++;
-
- } //# modules in y direction
- } //# modules in x direction
-
-
- cout << "PSD contains " << nModules << " modules." << endl;
- infoFile << "PSD contains " << nModules << " modules." << endl;
- // --------------------------------------------------------------------------
-
- // ----- Place PSD in top node (cave) -------------------------------------
- Double_t psdHalfLength = psdSizeZ + 0.5 * shieldWidth;
- Double_t psdVolCenterX = psdX + psdHalfLength * sin(psdRotY);
- Double_t psdVolCenterY = psdY;
- Double_t psdVolCenterZ = psdZ + psdHalfLength * cos(psdRotY);
- infoFile << "PSD volume center coordinates: (" << psdVolCenterX << ", " << psdVolCenterY << ", " << psdVolCenterZ
- << ") cm" << endl;
- TGeoRotation* psdRot = new TGeoRotation();
- psdRot->RotateY(psdRotY * TMath::RadToDeg());
- TGeoCombiTrans* psdTrans = new TGeoCombiTrans(psdVolCenterX, psdVolCenterY, psdVolCenterZ, psdRot);
- top->AddNode(psd, 0, psdTrans);
- cout << endl << "==> PSD position in cave: " << endl;
- psdTrans->Print();
- // --------------------------------------------------------------------------
-
-
- // ----- Close and check geometry ---------------------------------------
- cout << endl << "==> Closing geometry..." << endl;
- gGeoManager->CloseGeometry();
- gGeoManager->CheckGeometryFull();
- cout << endl;
- gGeoManager->CheckOverlaps(0.0001, "s");
- // --------------------------------------------------------------------------
-
-
- // ----- Write PSD volume and placement matrix to geometry file ---------
- cout << endl;
- TString geoFileName = "psd_";
- geoFileName = geoFileName + geoTag + ".geo.root";
- psd->Export(geoFileName);
- TFile* geoFile = new TFile(geoFileName, "UPDATE");
- psdTrans->Write();
- cout << endl;
- cout << "==> Geometry " << psd->GetName() << " written to " << geoFileName << endl;
- cout << "==> Info written to " << infoFileName << endl;
- geoFile->Close();
- infoFile.close();
- // --------------------------------------------------------------------------
-
-
- // ----- Write entire TGeoManager to file -------------------------------
- TString geoManFileName = "psd_";
- geoManFileName = geoManFileName + geoTag + ".geoman.root";
- TFile* geoManFile = new TFile(geoManFileName, "RECREATE");
- gGeoManager->Write();
- geoManFile->Close();
- cout << "==> TGeoManager " << gGeoManager->GetName() << " written to " << geoManFileName << endl << endl << endl;
- // --------------------------------------------------------------------------
-
-
- // ---- Display geometry ------------------------------------------------
- // gGeoManager->SetVisLevel(5); // Scintillator level
- top->Draw("ogl");
- // --------------------------------------------------------------------------
- // top->Raytrace();
-}
-// End of main function
-/** ======================================================================= **/
-
-TGeoVolume* ConstructShield(const char* name, Double_t sizeXY, Double_t holesize, Int_t hole_pos, fstream* info)
-{
-
- const TString suffix = Form("_shield_%d_%d_%d", int(sizeXY), int(holesize), int(hole_pos));
-
- const Double_t psd2iron = 2.;
- const Double_t ironD = 4.;
- const Double_t iron2poly = 2.;
- const Double_t polyD = 8.;
- const Double_t shieldLength = psd2iron + ironD + iron2poly + polyD;
-
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
-
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
-
- TGeoMedium* medShieldPoly = gGeoManager->GetMedium("PsdPolyethylene");
- if (!medShieldPoly) Fatal("Main", "Medium PsdPolyethylene not found");
-
- TGeoVolume* shield = gGeoManager->MakeBox(name, medAir, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * shieldLength);
- TGeoVolume* iron {nullptr};
- TGeoVolume* poly {nullptr};
-
- if (holesize > 0) {
- Int_t sx, sy;
- switch (hole_pos) {
- case 0:
- sx = 1;
- sy = 1;
- break; // defines hole position in x-y plane
- case 1:
- sx = 1;
- sy = -1;
- break;
- case 2:
- sx = -1;
- sy = 1;
- break;
- case 3:
- sx = -1;
- sy = -1;
- break;
- }
- TGeoCombiTrans* rot45 = new TGeoCombiTrans(Form("rot45%s", suffix.Data()), sx * sizeXY / 2., sy * sizeXY / 2., 0.,
- new TGeoRotation("rot", 45., 0., 0.));
- rot45->RegisterYourself();
-
- TGeoBBox* iron1 = new TGeoBBox(Form("iron1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * ironD);
- TGeoBBox* iron2 = new TGeoBBox(Form("iron2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * ironD); //hole
-
- TGeoCompositeShape* ironShape =
- new TGeoCompositeShape(Form("iron1%s-iron2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- iron = new TGeoVolume(Form("iron%s", suffix.Data()), ironShape, medIron);
-
- TGeoBBox* poly1 = new TGeoBBox(Form("poly1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * polyD);
- TGeoBBox* poly2 = new TGeoBBox(Form("poly2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * polyD); //hole
-
- TGeoCompositeShape* polyShape =
- new TGeoCompositeShape(Form("poly1%s-poly2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- poly = new TGeoVolume(Form("poly%s", suffix.Data()), polyShape, medShieldPoly);
- }
- else {
- iron = gGeoManager->MakeBox("shield_iron", medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * ironD);
- poly = gGeoManager->MakeBox("shield_poly", medShieldPoly, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * polyD);
- }
-
- shield->AddNode(iron, 0, new TGeoTranslation(0., 0., -shieldLength / 2. + psd2iron + ironD / 2.));
- shield->AddNode(poly, 0, new TGeoTranslation(0., 0., -shieldLength / 2. + psd2iron + ironD + iron2poly + polyD / 2.));
-
- return shield;
-}
-
-/** ======================================================================= **/
-TGeoVolume* ConstructModule(const char* name, Double_t sizeXY, Int_t nLayers, fstream* info)
-{
-
- // The module consists of nLayers of scintillators covered with Tyvek.
- // Between each two scintillators, there is a lead layer (total nLayers -1).
- // At the top, there is a gap to stretch the light fibres. This is modelled
- // by a volume made of the same material as the scintillator (but inactive).
- // The arrangement is surrounded by iron. At the front and at the back,
- // there is a thick iron layer (the back one acting as first absorber)
- // for constructional reasons.
-
- cout << "===> Creating Module " << name << ", size " << sizeXY << " cm with " << nLayers << " layers...." << endl;
-
- const TString suffix = Form("_%d_%d", int(sizeXY), nLayers); // unique suffix for volume names
-
- // ----- Constructional parameters ----------------------------
- const Double_t leadD = 1.60; // Thickness of lead layer
- const Double_t scintD = 0.40; // Thickness of scintillator layer
- const Double_t tyvekD = 0.02; // Thickness of Tyvek wrapper around scintillator
- const Double_t boxDx = 0.15; // Thickness of iron box lateral
- const Double_t boxDy = 0.05; // Thickness of iron box top/bottom
- const Double_t boxDz = 2.00; // Thickness of iron box front/back
- const Double_t chanDy = 0.20; // Height of fibre channel
- const Double_t chanDistL = 0.50; // Distance of channel from left edge of lead
- const Double_t chanDistR = 2.00; // Distance of channel from right edge of lead
- // ------------------------------------------------------------------------
-
-
- // ----- Info to file -------------------------------------------------
- if (info) {
- (*info) << "Parameters of module " << name << ": " << endl;
- (*info) << "Size: " << sizeXY << " cm x " << sizeXY << " cm" << endl;
- (*info) << "Number of layers: " << nLayers << endl;
- (*info) << "Thickness of lead layers: " << leadD << " cm" << endl;
- (*info) << "Thickness of scintillators: " << scintD << " cm" << endl;
- (*info) << "Thickness of Tyvek wrap: " << tyvekD << " cm" << endl;
- (*info) << "Thickness of iron box: (" << boxDx << " / " << boxDy << " / " << boxDz << ") cm" << endl;
- (*info) << "Height of fibre channel: " << chanDy << " cm" << endl;
- (*info) << "Distance of channel from edges: left " << chanDistL << " cm, right " << chanDistR << " cm" << endl;
- }
- // ------------------------------------------------------------------------
-
-
- // ----- Get required media -------------------------------------------
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
- TGeoMedium* medLead = gGeoManager->GetMedium("lead"); // Lead layers
- if (!medLead) Fatal("ConstructModule", "Medium lead not found");
- TGeoMedium* medTyvek = gGeoManager->GetMedium("PsdTyvek"); // Tyvek layers
- if (!medTyvek) Fatal("ConstructModule", "Medium PsdTyvek not found");
- TGeoMedium* medScint = gGeoManager->GetMedium("PsdScint"); // Scintillator
- if (!medScint) Fatal("ConstructModule", "Medium PsdScint not found");
- TGeoMedium* medFibre = gGeoManager->GetMedium("PsdFibre"); // Fibres
- if (!medFibre) Fatal("ConstructModule", "Medium PsdFibre not found");
- // ------------------------------------------------------------------------
-
-
- // ----- Create module volume -----------------------------------------
- // Module length: nLayers of scintillators, nLayers - 1 of lead
- // plus the iron box front and back.
- // Material is iron.
- const Double_t moduleLength = 2. * boxDz + nLayers * (scintD + 2. * tyvekD) + (nLayers - 1) * leadD;
- TGeoVolume* module = gGeoManager->MakeBox(name, medIron, 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * moduleLength);
- module->SetLineColor(kRed);
- module->Print();
- cout << endl;
- // ------------------------------------------------------------------------
-
-
- // ----- Lead filler --------------------------------------------------
- // The lead filler represents all lead absorber layers.
- // The Tyvek/scintillator layers and the fibre channel will be placed
- // inside.
- // Dimensions are half-lengths.
- const Double_t leadLx = 0.5 * sizeXY - boxDx;
- const Double_t leadLy = 0.5 * sizeXY - boxDy;
- const Double_t leadLz = 0.5 * moduleLength - boxDz;
- TGeoVolume* lead = gGeoManager->MakeBox(Form("lead%s", suffix.Data()), medLead, leadLx, leadLy, leadLz);
- lead->SetLineColor(kBlue);
- // ------------------------------------------------------------------------
-
-
- // ----- Fibre channel ------------------------------------------------
- // This volume represents the air gap at the top of the volume hosting the
- // light guide fibres. The material is the same as for the scintillators
- // (but inactive).
- // The (half-)width of the channel is defined by the distances from the
- // lead edges.
- const Double_t chanLx = leadLx - 0.5 * chanDistL - 0.5 * chanDistR;
- if (chanLx < 0.) {
- cout << "Error: lead volume is not large enough to host fibre channel!" << endl;
- cout << "Lead width: " << 2. * leadLx << ", distance from left edge " << chanDistL << ", distance from right edge "
- << chanDistR << endl;
- return 0;
- }
- const Double_t chanLy = 0.5 * chanDy;
- const Double_t chanLz = leadLz;
- TGeoVolume* channel = gGeoManager->MakeBox(Form("channel%s", suffix.Data()), medFibre, chanLx, chanLy, chanLz);
- channel->SetLineColor(kMagenta);
- // ------------------------------------------------------------------------
-
-
- // ---- Positioning of the fibre channel ------------------------------
- // The channel extends from chanDistL from the left edge of the lead filler
- // to chanDistR from its right edge (seen from the front). It is top-aligned
- // with the lead filler.
- const Double_t chanShiftX = 0.5 * (chanDistL - chanDistR);
- const Double_t chanShiftY = leadLy - 0.5 * chanDy;
- // ------------------------------------------------------------------------
- // cout << " leadLx = " << leadLx << " leadLy = " << leadLy << endl;
- // cout << " chanShiftX = " << chanShiftX << " chanShiftY = " << chanShiftY << endl;
-
- // ----- Tyvek layer --------------------------------------------------
- // The scintillator will be placed inside this, leaving only the thin
- // Tyvek as a wrapper. Since these layers will be placed in the lead filler,
- // there has to be a cut-out for the fibre channel.
- const Double_t tyvekLz = 0.5 * scintD + tyvekD; // half-length
- TGeoBBox* tyv1 = new TGeoBBox(Form("psdTyv1%s", suffix.Data()), leadLx, leadLy, tyvekLz);
- TGeoBBox* tyv2 = new TGeoBBox(Form("psdTyv2%s", suffix.Data()), chanLx, chanLy, tyvekLz);
- TGeoTranslation* trans1 = new TGeoTranslation(Form("tPsd1%s", suffix.Data()), chanShiftX, chanShiftY, 0.);
-
- trans1->RegisterYourself();
- TGeoCompositeShape* tyvekShape =
- new TGeoCompositeShape(Form("psdTyv1%s-psdTyv2%s:tPsd1%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* tyvek = new TGeoVolume(Form("tyvek%s", suffix.Data()), tyvekShape, medTyvek);
- tyvek->SetLineColor(kGreen);
- // ------------------------------------------------------------------------
-
-
- // ----- Scintillator layer -------------------------------------------
- // The scintillator layer is embedded (as daughter) into the Tyvek layer.
- // It is also a box minus the cut-out for the fibres. The cut-out is
- // slightly smaller than for the Tyvek volume.
- const Double_t scintLx = leadLx - tyvekD;
- const Double_t scintLy = leadLy - tyvekD;
- const Double_t scintLz = 0.5 * scintD;
- TGeoBBox* sci1 = new TGeoBBox(Form("sci1%s", suffix.Data()), scintLx, scintLy, scintLz);
- const Double_t cutLx = chanLx;
- const Double_t cutLy = chanLy - tyvekD;
- const Double_t cutLz = scintLz;
- TGeoBBox* sci2 = new TGeoBBox(Form("sci2%s", suffix.Data()), cutLx, cutLy, cutLz);
- const Double_t cutShiftX = chanShiftX;
- const Double_t cutShiftY = scintLy - cutLy;
- TGeoTranslation* trans2 = new TGeoTranslation(Form("tPsd2%s", suffix.Data()), cutShiftX, cutShiftY, 0.);
- trans2->RegisterYourself();
- TGeoCompositeShape* scintShape = new TGeoCompositeShape(
- Form("scintShape%s", suffix.Data()), Form("sci1%s-sci2%s:tPsd2%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* scint = new TGeoVolume(Form("scint%s", suffix.Data()), scintShape, medScint);
- scint->SetLineColor(kBlack);
-
- // ------------------------------------------------------------------------
- // cout << " scintLx = " << scintLx << " scintLy = " << scintLy << endl;
- // cout << " cutLx = " << cutLx << " cutLy = " << cutLy << endl;
- // cout << " cutShiftX = " << cutShiftX << " cutShiftY = " << cutShiftY << endl;
-
- // ----- Place volumes ------------------------------------------------
-
- // ---> Scintillator into Tyvek
- tyvek->AddNode(scint, 0);
-
- // ---> Fibre channel into lead filler
- lead->AddNode(channel, 0, trans1);
-
- // ---> Tyvek layers into lead
- Double_t zFirst {0.};
- Double_t zLast {0.};
- for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
- const Double_t zPos = -1. * leadLz + iLayer * leadD + (2. * iLayer + 1.) * tyvekLz;
- if (iLayer == 0) zFirst = zPos;
- if (iLayer == nLayers - 1) zLast = zPos;
- lead->AddNode(tyvek, iLayer, new TGeoTranslation(0., 0., zPos));
- }
- cout << module->GetName() << ": Positioned " << nLayers << " Tyvek layers; first at z = " << zFirst
- << ", last at z = " << zLast << endl;
-
- // ---> Lead into module
- module->AddNode(lead, 0);
- // ------------------------------------------------------------------------
-
- return module;
-}
-
-
-TGeoVolume* ConstructModuleWithHole(const char* name, Double_t sizeXY, Int_t nLayers, fstream* info, float holesize,
- Int_t hole_pos)
-{
-
- // The module consists of nLayers of scintillators covered with Tyvek.
- // Between each two scintillators, there is a lead layer (total nLayers -1).
- // At the top, there is a gap to stretch the light fibres. This is modelled
- // by a volume made of the same material as the scintillator (but inactive).
- // The arrangement is surrounded by iron. At the front and at the back,
- // there is a thick iron layer (the back one acting as first absorber)
- // for constructional reasons.
-
-
- if (holesize == 0.) return ConstructModule(name, sizeXY, nLayers, info);
-
- const TString suffix = Form("_%d_%d_%d", int(sizeXY), int(holesize), int(hole_pos));
-
- cout << "===> Creating Module " << name << ", size " << sizeXY << " cm with " << nLayers << " layers...." << endl;
-
- // ----- Constructional parameters ----------------------------
- const Double_t leadD = 1.60; // Thickness of lead layer
- const Double_t scintD = 0.40; // Thickness of scintillator layer
- const Double_t tyvekD = 0.02; // Thickness of Tyvek wrapper around scintillator
-
- const Double_t boxDx = 0.15; // Thickness of iron box lateral
- const Double_t boxDy = 0.05; // Thickness of iron box top/bottom
- const Double_t boxDz = 2.00; // Thickness of iron box front/back
-
- const Double_t chanDy = 0.20; // Height of fibre channel
- const Double_t chanDistL = 0.50; // Distance of channel from left edge of lead
- const Double_t chanDistR = 2.00; // Distance of channel from right edge of lead
-
- const Double_t moduleLength = 2. * boxDz + nLayers * (scintD + 2. * tyvekD) + (nLayers - 1) * leadD;
-
- // Dimensions are half-lengths.
- const Double_t leadLx = 0.5 * sizeXY - boxDx;
- const Double_t leadLy = 0.5 * sizeXY - boxDy;
- const Double_t leadLz = 0.5 * moduleLength - boxDz;
- const Double_t tyvekLz = 0.5 * scintD + tyvekD; // half-length
-
- const Double_t scintLx = leadLx - tyvekD; // Scintillator size x
- const Double_t scintLy = leadLy - tyvekD; // Scintillator size y
- const Double_t scintLz = 0.5 * scintD; // Scintillator size z
-
-
- // ------------------------------------------------------------------------
- Int_t sx, sy;
- switch (hole_pos) {
- case 0:
- sx = 1;
- sy = 1;
- break; // defines hole position in x-y plane
- case 1:
- sx = 1;
- sy = -1;
- break;
- case 2:
- sx = -1;
- sy = 1;
- break;
- case 3:
- sx = -1;
- sy = -1;
- break;
- }
- // ----- Info to file -------------------------------------------------
- if (info) {
- (*info) << "Parameters of module " << name << ": " << endl;
- (*info) << "Size: " << sizeXY << " cm x " << sizeXY << " cm" << endl;
- (*info) << "Number of layers: " << nLayers << endl;
- (*info) << "Thickness of lead layers: " << leadD << " cm" << endl;
- (*info) << "Thickness of scintillators: " << scintD << " cm" << endl;
- (*info) << "Thickness of Tyvek wrap: " << tyvekD << " cm" << endl;
- (*info) << "Thickness of iron box: (" << boxDx << " / " << boxDy << " / " << boxDz << ") cm" << endl;
- (*info) << "Height of fibre channel: " << chanDy << " cm" << endl;
- (*info) << "Distance of channel from edges: left " << chanDistL << " cm, right " << chanDistR << " cm" << endl;
- }
- // ------------------------------------------------------------------------
-
- // ----- Get required media -------------------------------------------
- TGeoMedium* medAir = gGeoManager->GetMedium("air"); // PSD
- if (!medAir) Fatal("ConstructModule", "Medium air not found");
- TGeoMedium* medIron = gGeoManager->GetMedium("iron"); // Box
- if (!medIron) Fatal("ConstructModule", "Medium iron not found");
- TGeoMedium* medLead = gGeoManager->GetMedium("lead"); // Lead layers
- if (!medLead) Fatal("ConstructModule", "Medium lead not found");
- TGeoMedium* medTyvek = gGeoManager->GetMedium("PsdTyvek"); // Tyvek layers
- if (!medTyvek) Fatal("ConstructModule", "Medium PsdTyvek not found");
- TGeoMedium* medScint = gGeoManager->GetMedium("PsdScint"); // Scintillator
- if (!medScint) Fatal("ConstructModule", "Medium PsdScint not found");
- TGeoMedium* medFibre = gGeoManager->GetMedium("PsdFibre"); // Fibres
- if (!medFibre) Fatal("ConstructModule", "Medium PsdFibre not found");
- // ------------------------------------------------------------------------
-
- TGeoCombiTrans* rot45 = new TGeoCombiTrans(Form("rot45%s", suffix.Data()), sx * sizeXY / 2., sy * sizeXY / 2., 0.,
- new TGeoRotation("rot", 45., 0., 0.));
- rot45->RegisterYourself();
-
- // ----- Create module volume -----------------------------------------
- // Module length: nLayers of scintillators, nLayers - 1 of lead
- // plus the iron box front and back.
- // Material is iron.
- // TGeoVolume* module = gGeoManager->MakeBox(name, medIron, 0.5 * sizeXY,
- // 0.5 * sizeXY, 0.5 * moduleLength);
-
- TGeoBBox* module1 = new TGeoBBox(Form("module1%s", suffix.Data()), 0.5 * sizeXY, 0.5 * sizeXY, 0.5 * moduleLength);
- TGeoBBox* module2 = new TGeoBBox(Form("module2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- 0.5 * moduleLength); //hole
- TGeoCompositeShape* moduleShape =
- new TGeoCompositeShape(Form("module1%s-module2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* module = new TGeoVolume(name, moduleShape, medIron);
-
- module->SetLineColor(kRed);
- module->Print();
- cout << endl;
- // ------------------------------------------------------------------------
-
-
- // ----- Lead filler --------------------------------------------------
- // The lead filler represents all lead absorber layers.
- // The Tyvek/scintillator layers and the fibre channel will be placed
- // inside.
- // TGeoVolume* lead = gGeoManager->MakeBox(Form("lead%s", suffix.Data()), medLead, leadLx, leadLy, leadLz);
-
- TGeoBBox* lead1 = new TGeoBBox(Form("lead1%s", suffix.Data()), leadLx, leadLy, leadLz);
- TGeoBBox* lead2 = new TGeoBBox(Form("lead2%s", suffix.Data()), 0.5 * holesize, 0.5 * holesize,
- leadLz); //hole
- TGeoCompositeShape* leadShape =
- new TGeoCompositeShape(Form("lead1%s-lead2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* lead = new TGeoVolume(Form("lead%s", suffix.Data()), leadShape, medLead);
-
- lead->SetLineColor(kBlue);
- // ------------------------------------------------------------------------
-
- // ----- Tyvek layer --------------------------------------------------
- // The scintillator will be placed inside this, leaving only the thin
- // Tyvek as a wrapper. Since these layers will be placed in the lead filler,
- // there has to be a cut-out for the fibre channel.
-
-
- TGeoBBox* tyv1 = new TGeoBBox(Form("psdTyv1%s", suffix.Data()), leadLx, leadLy, tyvekLz);
- TGeoBBox* tyv2 = new TGeoBBox(Form("psdTyv2%s", suffix.Data()), holesize / 2., holesize / 2.,
- tyvekLz); //hole
-
- TGeoCompositeShape* tyvekShape =
- new TGeoCompositeShape(Form("psdTyv1%s-psdTyv2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* tyvek = new TGeoVolume(Form("tyvek%s", suffix.Data()), tyvekShape, medTyvek);
- tyvek->SetLineColor(kGreen);
- // ------------------------------------------------------------------------
-
-
- // ----- Scintillator layer -------------------------------------------
- // The scintillator layer is embedded (as daughter) into the Tyvek layer.
- // It is also a box minus the cut-out for the fibres. The cut-out is
- // slightly smaller than for the Tyvek volume.
- TGeoBBox* sci1 = new TGeoBBox(Form("sci1%s", suffix.Data()), scintLx, scintLy, scintLz);
- TGeoBBox* sci2 = new TGeoBBox(Form("sci2%s", suffix.Data()), holesize / 2., holesize / 2.,
- scintLz); //hole
- TGeoCompositeShape* scintShape = new TGeoCompositeShape(
- Form("scintShape%s", suffix.Data()), Form("sci1%s-sci2%s:rot45%s", suffix.Data(), suffix.Data(), suffix.Data()));
- TGeoVolume* scint = new TGeoVolume(Form("scint%s", suffix.Data()), scintShape, medScint);
- scint->SetLineColor(kBlack);
-
- // ----- Place volumes ------------------------------------------------
-
- // ---> Scintillator into Tyvek
- tyvek->AddNode(scint, 0);
-
- // ---> Fibre channel into lead filler
- // lead->AddNode(channel, 0, trans1); //Needs addidional size check
-
- // ---> Tyvek layers into lead
- Double_t zFirst {0.};
- Double_t zLast {0.};
- for (Int_t iLayer = 0; iLayer < nLayers; iLayer++) {
- const Double_t zPos = -1. * leadLz + iLayer * leadD + (2. * iLayer + 1.) * tyvekLz;
- if (iLayer == 0) zFirst = zPos;
- if (iLayer == nLayers - 1) zLast = zPos;
- lead->AddNode(tyvek, iLayer, new TGeoTranslation(0., 0., zPos));
- }
- cout << module->GetName() << ": Positioned " << nLayers << " Tyvek layers; first at z = " << zFirst
- << ", last at z = " << zLast << endl;
-
- // ---> Lead into module
- module->AddNode(lead, 0);
- // ------------------------------------------------------------------------
-
- return module;
-}
--
GitLab