Skip to content
Snippets Groups Projects
Commit 79b3aeef authored by Administrator's avatar Administrator Committed by Florian Uhlig
Browse files

Move missing geometry generation macros to the geometry repository, refs #2474

parent e21283f9
Branches
Tags
1 merge request!797Remove geometry generation macros
Pipeline #16877 passed
Stage: documentation
Stage: deploy
Expand all Show whitespace changes
Inline Side-by-side
macro/passive/create_platform_v22b.C deleted 100644 → 0
+ 0
− 150
View file @ e21283f9
/* 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 =====
// ============================================================================
macro/psd/create_psdgeo.C deleted 100644 → 0
+ 0
− 437
View file @ e21283f9
/* 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;
}
/** ======================================================================= **/
macro/psd/create_psdgeo_46modules.C deleted 100644 → 0
+ 0
− 576
View file @ e21283f9
This diff is collapsed.
macro/psd/create_psdgeo_52modules.C deleted 100644 → 0
+ 0
− 777
View file @ e21283f9
This diff is collapsed.
macro/psd/create_psdgeo_ideal.C deleted 100644 → 0
+ 0
− 312
View file @ e21283f9
/* 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;
}
macro/psd/create_psdgeo_with_hole.C deleted 100644 → 0
+ 0
− 838
View file @ e21283f9
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment