create_rootgeom_HGND_run9.C

Go to the documentation of this file.

#include "../profiling/bmnloadlibs.C"
#include "TFile.h"
#include "TGeoManager.h"
#include "TGeoMedium.h"
#include "TGeoVolume.h"
 
class FairGeoMedia;
class FairGeoBuilder;
 
// ===========================================================================
// Detector geometry constants
// ===========================================================================
namespace NDetGeom
{
 
// Module grid: 1x1 (single detector element per layer)
constexpr Int_t kNRows = 1;
constexpr Int_t kNColumns = 1;
constexpr Int_t kNLayers = 2;
constexpr Int_t kNCellsPerLayer = kNRows * kNColumns;
 
// --- Full sizes, cm (NOT half-sizes) -----------------------------------
// TGeoBBox's constructor requires half-sizes -- that conversion is done
// exactly once, at each TGeoBBox(...) call below, via the Half() helper.
// Nowhere else in this file should a half-size ever appear.
constexpr Double_t kScintSizeX = 4.0, kScintSizeY = 4.0, kScintSizeZ = 2.5;   // scintillator cell
constexpr Double_t kG10SizeX = kNColumns * kScintSizeX;                       // PCB, same footprint as the cell
constexpr Double_t kG10SizeY = kNRows * kScintSizeY;
constexpr Double_t kG10SizeZ = 0.2;
 
constexpr Double_t kEnvelopeMargin = 0.0001;
 
// --- Layer positions, front-face-relative, cm ---------------------------
// z = 0 is the front face of the first layer. Both layers are identical
// (scintillator cell + PCB immediately behind it).
constexpr Double_t kScintFrontZ[kNLayers] = {0.0, 10};
 
constexpr Double_t kEnvelopeSizeX = kG10SizeX;
constexpr Double_t kEnvelopeSizeY = kG10SizeY;
constexpr Double_t kEnvelopeSizeZ = kScintFrontZ[kNLayers - 1] + kScintSizeZ + kG10SizeZ;
 
// ROOT's TGeoBBox is always centered on its own local origin -- there is no
// way around that. Half() converts a full size to the half-size the
// constructor wants; FrontToLocalZ() converts a front-face-relative z
// (0-based, always >= 0, the coordinate a physicist naturally writes down)
// into a z relative to the envelope's own centered local origin (what
// TGeoTranslation actually needs). Both conversions happen only inside
// these two helpers -- nowhere else in the file.
constexpr Double_t Half(Double_t fullSize)
{
    return fullSize / 2.0;
}
constexpr Double_t FrontToLocalZ(Double_t zFront)
{
    return zFront - Half(kEnvelopeSizeZ);
}
 
}   // namespace NDetGeom
 
// ===========================================================================
// Media
// ===========================================================================
TGeoMedium* pMedVacuum = nullptr;
TGeoMedium* pMedAir = nullptr;
TGeoMedium* pMedG10 = nullptr;
TGeoMedium* pMedScint = nullptr;
 
void DefineRequiredMedia(FairGeoMedia* geoMedia, FairGeoBuilder* geoBuild)
{
    auto loadMedium = [&](const char* name, TGeoMedium*& dest) {
        FairGeoMedium* m = geoMedia->getMedium(name);
        if (!m)
            Fatal("DefineRequiredMedia", "FairMedium '%s' not found", name);
        geoBuild->createMedium(m);
        dest = gGeoManager->GetMedium(name);
        if (!dest)
            Fatal("DefineRequiredMedia", "Medium '%s' not found after creation", name);
    };
 
    loadMedium("vacuum", pMedVacuum);
    loadMedium("air", pMedAir);
    loadMedium("G10", pMedG10);
    loadMedium("polystyrene", pMedScint);   // scintillator
}
 
// ===========================================================================
// Main geometry construction macro
// ===========================================================================
/**
 * HGND Run9 geometry creation macro
 * Two identical layers (scintillator cell + PCB),
 * no absorber between them -- a bare test-stand / cosmic-test geometry.
 *
 * \param distToDetector  Distance from target to the front face of detector, cm.
 * \param yShift          Vertical shift of the detector center, cm.
 * \param rotYdeg         Rotation of the detector around the Y axis, degrees.
 * \param rotXdeg         Rotation of the detector around the X axis, degrees.
 * \param geoTag          Optional short tag appended to the output geometry
 */
void create_rootgeom_HGND_run9(Double_t distToDetector = 699.7,
                               Double_t yShift = 0.0,
                               Double_t rotYdeg = -10.0,
                               Double_t rotXdeg = 0.0,
                               const char* geoTag = "")
{
    using namespace NDetGeom;
 
    bmnloadlibs();   // load libraries
 
    TString workDir = gSystem->Getenv("VMCWORKDIR");
    const TString geoDetectorName = "HGND";
 
    TString geoFileName = Form("%s/geometry/%s_dist%.3fcm_yshift%.3fcm_rotY%.2fdeg_rotX%.2fdeg%s%s.root",
                               workDir.Data(), geoDetectorName.Data(), distToDetector, yShift, rotYdeg, rotXdeg,
                               (geoTag && geoTag[0]) ? "_" : "", geoTag);
 
    // ----  global geometry parameters  ---------------------------------------
    FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
    FairGeoInterface* geoFace = geoLoad->getGeoInterface();
 
    // -------   Load media from media file   ----------------------------------
    TString medFile = workDir + "/geometry/media.geo";
    geoFace->setMediaFile(medFile);
    geoFace->readMedia();
    FairGeoMedia* geoMedia = geoFace->getMedia();
    FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
    DefineRequiredMedia(geoMedia, geoBuild);
 
    // --------------   Create geometry and global top volume  --------------------
    gGeoManager = (TGeoManager*)gROOT->FindObject("FAIRGeom");
    gGeoManager->SetName(geoDetectorName + "_geom");
 
    TGeoVolume* top = new TGeoVolumeAssembly(BmnNdetGeoPar::fTopVolumeName);
    top->SetMedium(pMedAir);
    gGeoManager->SetTopVolume(top);
 
    TGeoVolume* ndetCommon = new TGeoVolumeAssembly(BmnNdetGeoPar::fCommonVolumeName);
    ndetCommon->SetMedium(pMedAir);
 
    TGeoBBox* envShape = new TGeoBBox(BmnNdetGeoPar::fEnvelopeVolumeName, Half(kEnvelopeSizeX) + kEnvelopeMargin,
                                      Half(kEnvelopeSizeY) + kEnvelopeMargin, Half(kEnvelopeSizeZ) + kEnvelopeMargin);
    TGeoVolume* ndetEnvelope = new TGeoVolume(BmnNdetGeoPar::fEnvelopeVolumeName, envShape, pMedVacuum);
 
    TGeoBBox* detShape = new TGeoBBox(BmnNdetGeoPar::fDetectorVolumeName, Half(kEnvelopeSizeX), Half(kEnvelopeSizeY),
                                      Half(kEnvelopeSizeZ));
    TGeoVolume* detVolume = new TGeoVolume(BmnNdetGeoPar::fDetectorVolumeName, detShape, pMedAir);
    ndetEnvelope->AddNode(detVolume, 0, new TGeoTranslation(0., 0., 0.));
 
    // ----  material/slice volume prototypes (shared across all layers)  --------
    auto makeSliceVolume = [](const char* name, Double_t sizeX, Double_t sizeY, Double_t sizeZ, TGeoMedium* medium,
                              Color_t color) {
        TGeoShape* shape = new TGeoBBox(name, Half(sizeX), Half(sizeY), Half(sizeZ));
        TGeoVolume* vol = new TGeoVolume(name, shape);
        vol->SetLineColor(color);
        vol->SetFillColor(color);
        if (!medium)
            Fatal("create_ndet_geometry", "Invalid medium for volume '%s'", name);
        vol->SetMedium(medium);
        return vol;
    };
 
    TGeoVolume* volG10 =
        makeSliceVolume(BmnNdetGeoPar::G10_name_NICA, kG10SizeX, kG10SizeY, kG10SizeZ, pMedG10, kGreen);
    TGeoVolume* volScint = makeSliceVolume(BmnNdetGeoPar::SensitiveVolume_name_NICA, kScintSizeX, kScintSizeY,
                                           kScintSizeZ, pMedScint, kBlue);
 
    // ----  helper: place the single cell of one layer  --------------------------
    // Module numbering: modNb = row*NColumns + column + layerIndex*NCellsPerLayer + 1
    // (kNRows=kNColumns=1 here, so this is just layerIndex + 1, but the general
    // formula is kept so this still works if the grid is ever made bigger again).
    auto placeLayer = [&](Int_t layerIndex, Double_t zFrontFace) {
        const Double_t zScintCenter = FrontToLocalZ(zFrontFace + Half(kScintSizeZ));
        const Double_t zG10Center = FrontToLocalZ(zFrontFace + kScintSizeZ + Half(kG10SizeZ));
 
        for (Int_t iy = 0; iy < kNRows; ++iy) {
            for (Int_t ix = 0; ix < kNColumns; ++ix) {
                const Int_t modNb = iy * kNColumns + ix + layerIndex * kNCellsPerLayer + 1;
                const Int_t rowNb = iy + 1;
                const Int_t columnNb = ix + 1;
 
                TGeoBBox* moduleShape = new TGeoBBox(BmnNdetGeoPar::fModuleVolumeName, Half(kScintSizeX),
                                                     Half(kScintSizeY), Half(kScintSizeZ));
                TGeoVolume* moduleVolume = new TGeoVolume(BmnNdetGeoPar::fModuleVolumeName, moduleShape);
                moduleVolume->SetMedium(pMedAir);
                moduleVolume->SetName(
                    Form("%s_%d_%d_%d", BmnNdetGeoPar::fModuleVolumeName.Data(), rowNb, columnNb, layerIndex + 1));
                moduleVolume->AddNode(volScint, modNb, new TGeoTranslation(0, 0, 0));
 
                detVolume->AddNode(moduleVolume, modNb, new TGeoTranslation(0, 0, zScintCenter));
            }
        }
 
        detVolume->AddNode(volG10, layerIndex + 1, new TGeoTranslation(0, 0, zG10Center));
    };
 
    for (Int_t layer = 0; layer < kNLayers; ++layer) {
        placeLayer(layer, kScintFrontZ[layer]);
    }
 
    // ----  position and rotate the whole detector  ------------------------------
    TGeoRotation* rot = new TGeoRotation("rot");
    rot->RotateX(rotXdeg);
    rot->RotateY(rotYdeg);
 
    Double_t local[3] = {0., 0., distToDetector + Half(kEnvelopeSizeZ)};
    Double_t master[3] = {0., 0., 0.};
    rot->LocalToMaster(local, master);
    master[1] += yShift;
 
    printf("NDet center: x=%.3f y=%.3f z=%.3f cm\n", master[0], master[1], master[2]);
 
    TGeoCombiTrans* ndetTransform = new TGeoCombiTrans("ndetTransform", master[0], master[1], master[2], rot);
    ndetCommon->AddNode(ndetEnvelope, 1, ndetTransform);
    top->AddNode(ndetCommon, 0);
    top->SetVisContainers(kTRUE);
 
    // ---------------   Finish   ----------------------------------------------
    gGeoManager->CloseGeometry();
    gGeoManager->CheckOverlaps(0.001);
    gGeoManager->PrintOverlaps();
    gGeoManager->Test();
 
    TFile* geoFile = new TFile(geoFileName, "RECREATE");
    top->Write();
    geoFile->Close();
 
    printf("Geometry written to: %s\n", geoFileName.Data());
}