// Copyright 2019-2020 CERN and copyright holders of ALICE O2.

// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.

// All rights not expressly granted are reserved.

//

// This software is distributed under the terms of the GNU General Public

// License v3 (GPL Version 3), copied verbatim in the file "COPYING".

//

// In applying this license CERN does not waive the privileges and immunities

// granted to it by virtue of its status as an Intergovernmental Organization

// or submit itself to any jurisdiction.

/// \file Geometry.cxx

/// \brief Implementation of FV0 geometry.

///

/// \author Maciej Slupecki, University of Jyvaskyla, Finland

/// \author Andreas Molander, University of Helsinki, Finland

#include "FV0Base/Geometry.h"

#include <cmath>

#include <fairlogger/Logger.h>

#include <TGeoBBox.h>

#include <TGeoCompositeShape.h>

#include <TGeoCone.h>

#include <TGeoManager.h>

#include <TGeoMatrix.h>

#include <TGeoMedium.h>

#include <TGeoTube.h>

#include <TGeoVolume.h>

ClassImp(o2::fv0::Geometry);

using namespace o2::fv0;

Geometry::Geometry(EGeoType initType) : mGeometryType(initType)

{

initializeCellCenters();

initializeReadoutCenters();

if (initType != eUninitialized) {

initializeGeometry();

}

}

Geometry::Geometry(const Geometry& geometry) : mGeometryType(geometry.mGeometryType), mLeftTransformation(nullptr), mRightTransformation(nullptr)

{

this->mEnabledComponents = geometry.mEnabledComponents;

}

Geometry::~Geometry()

{

if (mRightTransformation) {

delete mRightTransformation;

}

if (mLeftTransformation) {

delete mLeftTransformation;

}

}

int Geometry::getCurrentCellId(const TVirtualMC* fMC) const

{

int detectorHalfID = -1;

int sectorID = -1;

int ringID = -1;

fMC->CurrentVolOffID(2, detectorHalfID);

fMC->CurrentVolOffID(1, sectorID);

fMC->CurrentVolOffID(0, ringID);

sectorID += detectorHalfID * sNumberOfCellSectors;

LOG(debug) << "FV0 Geometry::getCurrentCellId(): \n"

<< "Half id: " << detectorHalfID << "\n"

<< "Half name: " << fMC->CurrentVolOffName(2) << "\n"

<< "Sector id: " << sectorID << "\n"

<< "Sector name: " << fMC->CurrentVolOffName(1) << "\n"

<< "Ring id: " << ringID << "\n"

<< "Ring name: " << fMC->CurrentVolOffName(0) << "\n"

<< "Cell id : " << sectorID + 8 * ringID << "\n";

return sectorID + 8 * ringID;

}

bool Geometry::enableComponent(const EGeoComponent component, const bool enable)

{

if (mEnabledComponents.find(component) == mEnabledComponents.end()) {

LOG(debug) << "FV0 Geometry::enableComponent(): Component not initialized and cannot be enabled/disabled!";

return false;

}

return mEnabledComponents[component] = enable;

}

void Geometry::buildGeometry() const

{

TGeoVolume* vALIC = gGeoManager->GetVolume("barrel");

if (!vALIC) {

LOG(fatal) << "FV0: Could not find the top volume";

}

// Top volume of FV0 detector

TGeoVolumeAssembly* vFV0 = new TGeoVolumeAssembly(sDetectorName.c_str());

LOG(info) << "FV0: Building geometry. FV0 volume name is '" << vFV0->GetName() << "'";

TGeoVolumeAssembly* vFV0Right = new TGeoVolumeAssembly((sDetectorName + "RIGHT").c_str());

TGeoVolumeAssembly* vFV0Left = new TGeoVolumeAssembly((sDetectorName + "LEFT").c_str());

vFV0->AddNode(vFV0Right, 0, mRightTransformation);

vFV0->AddNode(vFV0Left, 1, mLeftTransformation);

assembleSensVols(vFV0Right, vFV0Left);

assembleNonSensVols(vFV0Right, vFV0Left);

vALIC->AddNode(vFV0, 1, new TGeoTranslation(sXGlobal, sYGlobal + 30., sZGlobal));

}

void Geometry::getGlobalPosition(float& x, float& y, float& z)

{

x = sXGlobal;

y = sYGlobal;

z = sZGlobal;

}

Point3Dsimple& Geometry::getCellCenter(UInt_t cellId)

{

return mCellCenter.at(cellId);

}

Point3Dsimple& Geometry::getReadoutCenter(UInt_t cellId)

{

return mReadoutCenter.at(cellId);

}

bool Geometry::isRing5(UInt_t cellId)

{

return cellId >= (sNumberOfCellRings - 1) * sNumberOfCellSectors * 2;

}

void Geometry::initializeGeometry()

{

initializeMaps();

initializeVectors();

initializeTransformations();

initializeSensVols();

initializeNonSensVols();

}

void Geometry::initializeMaps()

{

const bool isFull = mGeometryType == eFull;

const bool isRough = isFull || mGeometryType == eRough;

const bool hasScint = isRough || mGeometryType == eOnlySensitive;

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(eScintillator, hasScint));

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(ePlastics, isRough));

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(ePmts, isFull));

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(eFibers, isFull));

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(eScrews, isFull));

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(eRods, isFull));

mEnabledComponents.insert(std::pair<EGeoComponent, bool>(eContainer, isRough));

}

void Geometry::initializeVectors()

{

initializeCellRingRadii();

initializeSectorTransformations();

initializeFiberVolumeRadii();

initializeFiberMedium();

initializeScrewAndRodRadii();

initializeScrewTypeMedium();

initializeRodTypeMedium();

initializeScrewAndRodPositionsAndDimensions();

}

void Geometry::initializeCellRingRadii()

{

// Index of mRAvgRing is NOT linked directly to any ring number

mRAvgRing.assign(sCellRingRadii, sCellRingRadii + sNumberOfCellRings + 1);

// Set real scintillator radii (reduced by paint thickness and separation gap)

for (int i = 0; i < mRAvgRing.size() - 1; ++i) {

mRMinScintillator.push_back(mRAvgRing[i] + sDrSeparationScint);

mRMaxScintillator.push_back(mRAvgRing[i + 1] - sDrSeparationScint);

}

// Now indices of mRMinScint and mRMaxScint correspond to the same ring

}

void Geometry::initializeSectorTransformations()

{

for (int iSector = 0; iSector < sNumberOfCellSectors; ++iSector) {

// iSector = 0 corresponds to the first sector clockwise from the y-axis

// iSector = 1 corresponds to the next sector in clockwise direction and so on

TGeoRotation* trans = createAndRegisterRot(sDetectorName + sSectorName + std::to_string(iSector) + "TRANS");

if (iSector == 2 || iSector == 3) {

// "a" and "b" mirrors.

// The reference to "a" and "b" can be understood with the CAD drawings of the detector.

trans->ReflectY(true);

}

mSectorTrans.push_back(trans);

}

}

void Geometry::initializeFiberVolumeRadii()

{

for (int i = 0; i < sNumberOfCellRings; i++) {

mRMinFiber.push_back(sCellRingRadii[i] + sEpsilon / 2);

mRMaxFiber.push_back(sCellRingRadii[i + 1] - sEpsilon / 2);

}

}

void Geometry::initializeFiberMedium()

{

TGeoMedium* medium;

TString mediumName;

// one fiber volume per ring

for (int i = 0; i < mRMinFiber.size(); i++) {

mediumName = Form("FV0_FiberRing%i$", i + 1);

medium = gGeoManager->GetMedium(mediumName);

if (!medium) {

LOG(warning) << Form("FV0 geometry: Fiber medium for ring no. %i (%s) not found!", i + 1, mediumName.Data());

}

mMediumFiberRings.push_back(medium);

}

// five fiber volumes in front of the PMTs, one from each scintillator cell (two volumes from cell 5 but they are identical)

for (int i = 0; i < sNumberOfPMTFiberVolumes; i++) {

mediumName = Form("FV0_FiberPMT%i$", i + 1);

medium = gGeoManager->GetMedium(mediumName);

if (!medium) {

LOG(warning) << Form("FV0 geometry: PMT fiber medium from cell no. %i (%s) not found!", i + 1, mediumName.Data());

}

mMediumFiberPMTs.push_back(medium);

}

}

void Geometry::initializeScrewAndRodRadii()

{

mRScrewAndRod.push_back(mRAvgRing[1]);

mRScrewAndRod.push_back(mRAvgRing[2]);

mRScrewAndRod.push_back(mRAvgRing[3]);

mRScrewAndRod.push_back(mRAvgRing[4]);

mRScrewAndRod.push_back((mRAvgRing[4] + mRAvgRing[5]) / 2);

mRScrewAndRod.push_back(mRAvgRing[5]);

}

void Geometry::initializeScrewTypeMedium()

{

// There are no further checks if the medium is actually found

TGeoMedium* medium = gGeoManager->GetMedium("FV0_Titanium$");

for (int i = 0; i < sNumberOfScrewTypes; ++i) {

mMediumScrewTypes.push_back(medium);

}

}

void Geometry::initializeRodTypeMedium()

{

// There are no further checks if the medium is actually found

TGeoMedium* medium = gGeoManager->GetMedium("FV0_Aluminium$");

for (int i = 0; i < sNumberOfRodTypes; ++i) {

mMediumRodTypes.push_back(medium);

}

}

void Geometry::addScrewProperties(const int screwTypeID, const int iRing, const float phi)

{

float r = mRScrewAndRod[iRing];

mScrewTypeIDs.push_back(screwTypeID);

mScrewPos.push_back(std::vector<float>{cosf(phi * M_PI / 180) * r,

sinf(phi * M_PI / 180) * r,

sZScintillator - sDzScintillator / 2 + sZShiftScrew + sDzMaxScrewTypes[screwTypeID] / 2});

mDrMinScrews.push_back(sDrMinScrewTypes[screwTypeID]);

mDrMaxScrews.push_back(sDrMaxScrewTypes[screwTypeID]);

mDzMaxScrews.push_back(sDzMaxScrewTypes[screwTypeID]);

mDzMinScrews.push_back(sDzMinScrewTypes[screwTypeID]);

}

void Geometry::addRodProperties(const int rodTypeID, const int iRing)

{

mRodTypeIDs.push_back(rodTypeID);

mRodPos.push_back(std::vector<float>{sDxMinRodTypes[rodTypeID] / 2,

mRScrewAndRod[iRing],

sZScintillator - sDzScintillator / 2 + sZShiftRod + sDzMaxRodTypes[rodTypeID] / 2});

mDxMinRods.push_back(sDxMinRodTypes[rodTypeID]);

mDzMaxRods.push_back(sDxMaxRodTypes[rodTypeID]);

mDyMinRods.push_back(sDyMinRodTypes[rodTypeID]);

mDyMaxRods.push_back(sDyMaxRodTypes[rodTypeID]);

mDzMaxRods.push_back(sDzMaxRodTypes[rodTypeID]);

mDzMinRods.push_back(sDzMinRodTypes[rodTypeID]);

mRodTypeIDs.push_back(rodTypeID);

mRodPos.push_back(std::vector<float>{sDxMinRodTypes[rodTypeID] / 2,

-mRScrewAndRod[iRing],

sZScintillator - sDzScintillator / 2 + sZShiftRod + sDzMaxRodTypes[rodTypeID] / 2});

mDxMinRods.push_back(sDxMinRodTypes[rodTypeID]);

mDzMaxRods.push_back(sDxMaxRodTypes[rodTypeID]);

mDyMinRods.push_back(sDyMinRodTypes[rodTypeID]);

mDyMaxRods.push_back(sDyMaxRodTypes[rodTypeID]);

mDzMaxRods.push_back(sDzMaxRodTypes[rodTypeID]);

mDzMinRods.push_back(sDzMinRodTypes[rodTypeID]);

}

void Geometry::initializeScrewAndRodPositionsAndDimensions()

{

for (int iRing = 0; iRing < mRScrewAndRod.size(); ++iRing) {

switch (iRing) {

case 0:

addRodProperties(0, iRing);

for (float phi = 45; phi >= -45; phi -= 45) {

addScrewProperties(0, iRing, phi);

}

break;

case 1:

addRodProperties(0, iRing);

for (float phi = 45; phi >= -45; phi -= 45) {

addScrewProperties(1, iRing, phi);

}

break;

case 2:

addRodProperties(1, iRing);

for (float phi = 67.5; phi >= -67.5; phi -= 22.5) {

addScrewProperties(2, iRing, phi);

}

break;

case 3:

addRodProperties(2, iRing);

for (float phi = 67.5; phi >= -67.5; phi -= 22.5) {

addScrewProperties(3, iRing, phi);

}

break;

case 4:

addRodProperties(3, iRing);

for (float phi = 45; phi >= -45; phi -= 45) {

addScrewProperties(4, iRing, phi);

}

break;

case 5:

addRodProperties(3, iRing);

for (float phi = 67.5; phi >= -67.5; phi -= 22.5) {

addScrewProperties(5, iRing, phi);

}

break;

default:

break;

}

}

}

void Geometry::initializeTransformations()

{

TGeoTranslation leftTranslation(-sDxHalvesSeparation / 2, 0, 0);

TGeoRotation leftRotation;

leftRotation.ReflectX(true);

TGeoHMatrix leftTotalTransformation = leftTranslation * leftRotation;

mLeftTransformation = new TGeoHMatrix(leftTotalTransformation);

mRightTransformation = new TGeoTranslation(sDxHalvesSeparation / 2, sDyHalvesSeparation, sDzHalvesSeparation);

}

void Geometry::initializeSensVols()

{

initializeScintCells();

}

void Geometry::initializeNonSensVols()

{

initializeScrewHoles();

initializeRodHoles();

initializePlasticCells();

initializePmts();

initializeFibers();

initializeScrews();

initializeRods();

initializeMetalContainer();

}

void Geometry::initializeScrewHoles()

{

std::string boolFormula = "";

for (int i = 0; i < mScrewPos.size(); ++i) {

std::string holeShapeName = sDetectorName + sScrewName + "HOLE" + std::to_string(i);

std::string holeTransName = sDetectorName + sScrewName + "HOLETRANS" + std::to_string(i);

createScrewShape(holeShapeName, mScrewTypeIDs[i], sEpsilon, sEpsilon);

createAndRegisterTrans(holeTransName, mScrewPos[i][0] + sXShiftScrews, mScrewPos[i][1], mScrewPos[i][2]);

boolFormula += ((i != 0) ? "+" : "") + holeShapeName + ":" + holeTransName;

}

new TGeoCompositeShape(sScrewHolesCSName.c_str(), boolFormula.c_str());

}

void Geometry::initializeRodHoles()

{

std::string boolFormula = "";

for (int i = 0; i < mRodPos.size(); ++i) {

std::string holeShapeName = sDetectorName + sRodName + "HOLE" + std::to_string(i);

std::string holeTransName = sDetectorName + sRodName + "HOLETRANS" + std::to_string(i);

createRodShape(holeShapeName, mRodTypeIDs[i], sEpsilon, sEpsilon);

createAndRegisterTrans(holeTransName, mRodPos[i][0] + sXShiftScrews, mRodPos[i][1], mRodPos[i][2]);

boolFormula += ((i != 0) ? "+" : "") + holeShapeName + ":" + holeTransName;

}

new TGeoCompositeShape(sRodHolesCSName.c_str(), boolFormula.c_str());

}

void Geometry::initializeCells(const std::string& cellType, const float zThickness, const TGeoMedium* medium,

const bool isSensitive)

{

// Creating the two types of cells, "a" and "b", for each ring.

// All sectors can be assembled with these cells.

//

// The reference to "a" and "b" can be understood with the CAD drawings of the detector.

const float dxHoleCut = sDxHoleExtensionScintillator; // width of extension of hole 1, 2 and 7 in the "a" cell

const float xHole = sDrSeparationScint + dxHoleCut; // x-placement of holes 1, 2 and 7 in the "a" cell

// Sector separation gap shape

const std::string secSepShapeName = "FV0_" + cellType + "SectorSeparation";

new TGeoBBox(secSepShapeName.c_str(), mRMaxScintillator.back() + sEpsilon, sDrSeparationScint, zThickness / 2);

// Sector separation gap rotations

const std::string secSepRot45Name = "FV0_" + cellType + "SecSepRot45";

const std::string secSepRot90Name = "FV0_" + cellType + "SecSepRot90";

createAndRegisterRot(secSepRot45Name, 45, 0, 0);

createAndRegisterRot(secSepRot90Name, 90, 0, 0);

// Hole shapes

const std::string holeSmallName = "FV0_" + cellType + "HoleSmall";

const std::string holeLargeName = "FV0_" + cellType + "HoleLarge";

const std::string holeSmallCutName = "FV0_" + cellType + "HoleSmallCut";

const std::string holeLargeCutName = "FV0_" + cellType + "HoleLargeCut";

new TGeoTube(holeSmallName.c_str(), 0, sDrHoleSmallScintillator, zThickness / 2);

new TGeoTube(holeLargeName.c_str(), 0, sDrHoleLargeScintillator, zThickness / 2);

new TGeoBBox(holeSmallCutName.c_str(), dxHoleCut, sDrHoleSmallScintillator, zThickness / 2);

new TGeoBBox(holeLargeCutName.c_str(), dxHoleCut, sDrHoleLargeScintillator, zThickness / 2);

for (int ir = 0; ir < sNumberOfCellRings; ++ir) {

// Radii without separation

const float rMin = mRAvgRing[ir];

const float rMax = mRAvgRing[ir + 1];

const float rMid = rMin + (rMax - rMin) / 2;

// "a"-type cell

//

// Initial placement:

//

// y

// ^

// | 1******

// | ************5

// | 7*****************

// | *********************3

// | *******************

// | 2**************8

// | 6********

// | **4

// |

// |

// | O

// ------------------------> x

//

// * = cell volume

// numbers = hole numbers (as numbered in the code below)

// O = beam pipe

const std::string aCellName = createVolumeName(cellType + sCellName + "a", ir);

// Base shape

const std::string aCellShapeName = aCellName + "Shape";

// The cells in the innermost ring have a slightly shifted inner radius origin.

if (ir == 0) {

// The innermost "a"-type cell

const std::string a1CellShapeFullName = aCellShapeName + "Full";

const std::string a1CellShapeHoleCutName = aCellShapeName + "HoleCut";

const std::string a1CellShapeHoleCutTransName = a1CellShapeHoleCutName + "Trans";

new TGeoTubeSeg(a1CellShapeFullName.c_str(), 0, mRMaxScintillator[ir], zThickness / 2 - sEpsilon, 45, 90);

new TGeoTube(a1CellShapeHoleCutName.c_str(), 0, mRMinScintillator[ir], zThickness);

createAndRegisterTrans(a1CellShapeHoleCutTransName, sXShiftInnerRadiusScintillator, 0, 0);

const std::string a1BoolFormula = a1CellShapeFullName + "-" + a1CellShapeHoleCutName + ":" + a1CellShapeHoleCutTransName;

new TGeoCompositeShape(aCellShapeName.c_str(), a1BoolFormula.c_str());

} else {

// The rest of the "a"-type cells

new TGeoTubeSeg(aCellShapeName.c_str(), mRMinScintillator[ir], mRMaxScintillator[ir], zThickness / 2, 45, 90);

}

// Translations for screw holes (inner = rmin, half-length = rmid, outer = rmax)

//

// 1 = outer left

// 2 = inner left

// 3 = outer right

// 4 = inner right

// 5 = outer middle

// 6 = inner middle

// 7 = half-length left

// 8 = half-length right

//

// holes 1, 2 and 7 are slightly shifted along the rim of the cell

const std::string aHole1TransName = aCellName + "Hole1Trans";

const std::string aHole2TransName = aCellName + "Hole2Trans";

const std::string aHole3TransName = aCellName + "Hole3Trans";

const std::string aHole4TransName = aCellName + "Hole4Trans";

const std::string aHole5TransName = aCellName + "Hole5Trans";

const std::string aHole6TransName = aCellName + "Hole6Trans";

const std::string aHole7TransName = aCellName + "Hole7Trans";

const std::string aHole8TransName = aCellName + "Hole8Trans";

const std::string aHole1CutTransName = aCellName + "Hole1CutTrans";

const std::string aHole2CutTransName = aCellName + "Hole2CutTrans";

const std::string aHole7CutTransName = aCellName + "Hole7CutTrans";

createAndRegisterTrans(aHole1TransName, xHole, cos(asin(xHole / rMax)) * rMax, 0);

createAndRegisterTrans(aHole2TransName, xHole, cos(asin(xHole / rMin)) * rMin, 0);

createAndRegisterTrans(aHole3TransName, sin(45 * M_PI / 180) * rMax, cos(45 * M_PI / 180) * rMax, 0);

createAndRegisterTrans(aHole4TransName, sin(45 * M_PI / 180) * rMin, cos(45 * M_PI / 180) * rMin, 0);

createAndRegisterTrans(aHole5TransName, sin(22.5 * M_PI / 180) * rMax, cos(22.5 * M_PI / 180) * rMax, 0);

createAndRegisterTrans(aHole6TransName, sin(22.5 * M_PI / 180) * rMin, cos(22.5 * M_PI / 180) * rMin, 0);

createAndRegisterTrans(aHole7TransName, xHole, cos(asin(xHole / rMid)) * rMid, 0);

createAndRegisterTrans(aHole8TransName, sin(45 * M_PI / 180) * rMid, cos(45 * M_PI / 180) * rMid, 0);

createAndRegisterTrans(aHole1CutTransName, 0, cos(asin(xHole / rMax)) * rMax, 0);

createAndRegisterTrans(aHole2CutTransName, 0, cos(asin(xHole / rMin)) * rMin, 0);

createAndRegisterTrans(aHole7CutTransName, 0, cos(asin(xHole / rMid)) * rMid, 0);

// Composite shape

std::string aBoolFormula = aCellShapeName;

// sector separation

aBoolFormula += "-" + secSepShapeName + ":" + secSepRot45Name;

aBoolFormula += "-" + secSepShapeName + ":" + secSepRot90Name;

// outer holes

aBoolFormula += "-" + ((ir < 2) ? holeSmallName : holeLargeName) + ":" + aHole1TransName;

aBoolFormula += "-" + ((ir < 2) ? holeSmallCutName : holeLargeCutName) + ":" + aHole1CutTransName;

aBoolFormula += "-" + ((ir < 2) ? holeSmallName : holeLargeName) + ":" + aHole3TransName;

// inner holes

if (ir > 0) {

const std::string screwHoleName = (ir < 3) ? holeSmallName : holeLargeName;

const std::string screwHoleCutName = (ir < 3) ? holeSmallCutName : holeLargeCutName;

aBoolFormula += "-" + screwHoleName + ":" + aHole2TransName;

aBoolFormula += "-" + screwHoleCutName + ":" + aHole2CutTransName;

aBoolFormula += "-" + screwHoleName + ":" + aHole4TransName;

}

// outer middle hole

if (ir > 1) {

aBoolFormula += "-" + holeLargeName + ":" + aHole5TransName;

}

// inner middle hole

if (ir > 2) {

aBoolFormula += "-" + holeLargeName + ":" + aHole6TransName;

}

// half-length holes

if (ir == 4) {

aBoolFormula += "-" + holeLargeName + ":" + aHole7TransName;

aBoolFormula += "-" + holeLargeCutName + ":" + aHole7CutTransName;

aBoolFormula += "-" + holeLargeName + ":" + aHole8TransName;

}

const std::string aCellCSName = aCellName + "CS";

const TGeoCompositeShape* aCellCs = new TGeoCompositeShape(aCellCSName.c_str(), aBoolFormula.c_str());

// Cell volume

const TGeoVolume* aCell = new TGeoVolume(aCellName.c_str(), aCellCs, medium);

// "b"-type cell

//

// Initial placement:

//

// y

// ^

// | 1

// | *****

// | *********

// | 7*************

// | ****************5

// | 2*******************

// | 6*****************

// | *****************

// | O 4*******8******3

// |

// ------------------------------> x

//

// * = cell volume

// numbers = number of holes (as numbered in the code below)

// O = beam pipe

const std::string bCellName = createVolumeName(cellType + sCellName + "b", ir);

// Base shape

const std::string bCellShapeName = bCellName + "Shape";

// The cells in the innermost ring are slightly different than the rest

if (ir == 0) {

// The innermost "b"-type cell

const std::string b1CellShapeFullName = bCellShapeName + "Full";

const std::string b1CellShapeHoleCutName = bCellShapeName + "Cut";

const std::string b1CellShapeHoleCutTransName = b1CellShapeHoleCutName + "Trans";

new TGeoTubeSeg(b1CellShapeFullName.c_str(), 0, mRMaxScintillator[ir], zThickness / 2 - sEpsilon, 0, 45);

new TGeoTube(b1CellShapeHoleCutName.c_str(), 0, mRMinScintillator[ir], zThickness);

createAndRegisterTrans(b1CellShapeHoleCutTransName, sXShiftInnerRadiusScintillator, 0, 0);

const std::string b1BoolFormula = b1CellShapeFullName + "-" + b1CellShapeHoleCutName + ":" + b1CellShapeHoleCutTransName;

new TGeoCompositeShape(bCellShapeName.c_str(), b1BoolFormula.c_str());

} else {

// The rest of the "b"-type cells

new TGeoTubeSeg(bCellShapeName.c_str(), mRMinScintillator[ir], mRMaxScintillator[ir], zThickness / 2, 0, 45);

}

// Translations for holes

//

// 1 = outer left

// 2 = inner left

// 3 = outer right

// 4 = inner right

// 5 = outer middle

// 6 = inner middle

// 7 = half-lenght left

// 8 = half-length right

const std::string bHole1TransName = bCellName + "Hole1Trans";

const std::string bHole2TransName = bCellName + "Hole2Trans";

const std::string bHole3TransName = bCellName + "Hole3Trans";

const std::string bHole4TransName = bCellName + "Hole4Trans";

const std::string bHole5TransName = bCellName + "Hole5Trans";

const std::string bHole6TransName = bCellName + "Hole6Trans";

const std::string bHole7TransName = bCellName + "Hole7Trans";

const std::string bHole8TransName = bCellName + "Hole8Trans";

createAndRegisterTrans(bHole1TransName, sin(45 * M_PI / 180) * rMax, cos(45 * M_PI / 180) * rMax, 0);

createAndRegisterTrans(bHole2TransName, sin(45 * M_PI / 180) * rMin, cos(45 * M_PI / 180) * rMin, 0);

createAndRegisterTrans(bHole3TransName, rMax, 0, 0);

createAndRegisterTrans(bHole4TransName, rMin, 0, 0);

createAndRegisterTrans(bHole5TransName, cos(22.5 * M_PI / 180) * rMax, sin(22.5 * M_PI / 180) * rMax, 0);

createAndRegisterTrans(bHole6TransName, cos(22.5 * M_PI / 180) * rMin, sin(22.5 * M_PI / 180) * rMin, 0);

createAndRegisterTrans(bHole7TransName, sin(45 * M_PI / 180) * rMid, cos(45 * M_PI / 180) * rMid, 0);

createAndRegisterTrans(bHole8TransName, rMid, 0, 0);

// Composite shape

std::string bBoolFormula = bCellShapeName;

// sector separation

bBoolFormula += "-" + secSepShapeName;

bBoolFormula += "-" + secSepShapeName + ":" + secSepRot45Name;

// outer holes

bBoolFormula += "-" + ((ir < 2) ? holeSmallName : holeLargeName) + ":" + bHole1TransName;

bBoolFormula += "-" + ((ir < 2) ? holeSmallName : holeLargeName) + ":" + bHole3TransName;

// inner holes

if (ir > 0) {

const std::string holeName = (ir < 3) ? holeSmallName : holeLargeName;

bBoolFormula += "-" + holeName + ":" + bHole2TransName;

bBoolFormula += "-" + holeName + ":" + bHole4TransName;

}

// outer middle hole

if (ir > 1) {

bBoolFormula += "-" + holeLargeName + ":" + bHole5TransName;

}

// inner middle hole

if (ir > 2) {

bBoolFormula += "-" + holeLargeName + ":" + bHole6TransName;

}

// half-lenght holes

if (ir == 4) {

bBoolFormula += "-" + holeLargeName + ":" + bHole7TransName;

bBoolFormula += "-" + holeLargeName + ":" + bHole8TransName;

}

const std::string bCellCSName = bCellName + "CS";

const TGeoCompositeShape* bCellCs = new TGeoCompositeShape(bCellCSName.c_str(), bBoolFormula.c_str());

// Cell volume

const TGeoVolume* bCell = new TGeoVolume(bCellName.c_str(), bCellCs, medium);

if (isSensitive) {

mSensitiveVolumeNames.push_back(aCell->GetName());

mSensitiveVolumeNames.push_back(bCell->GetName());

}

}

}

void Geometry::initializeScintCells()

{

const TGeoMedium* medium = gGeoManager->GetMedium("FV0_Scintillator$");

initializeCells(sScintillatorName, sDzScintillator, medium, true);

}

void Geometry::initializePlasticCells()

{

const TGeoMedium* medium = gGeoManager->GetMedium("FV0_Plastic$");

initializeCells(sPlasticName, sDzPlastic, medium, false);

}

void Geometry::initializePmts()

{

const TGeoMedium* medium = gGeoManager->GetMedium("FV0_PMT$");

new TGeoVolume(createVolumeName(sPmtName).c_str(), new TGeoTube(createVolumeName(sPmtName + "Shape").c_str(), 0, sDrPmt, sDzPmt / 2), medium);

}

void Geometry::initializeFibers()

{

// depth of the fiber volumes

const float dzFibers = sDzContainer - sDzContainerBack - sDzContainerFront - sDzScintillator - sDzPlastic - 2 * sEpsilon;

const std::string fiberName = sDetectorName + "Fibers";

const std::string fiberSepCutName = fiberName + "SepCut";

const std::string fiberConeCutName = fiberName + "ConeCut";

const std::string fiberHoleCutName = fiberName + "HoleCut";

const std::string fiberTransName = fiberName + "Trans";

const std::string fiberConeCutTransName = fiberConeCutName + "Trans";

const std::string fiberHoleCutTransName = fiberHoleCutName + "Trans";

new TGeoBBox(fiberSepCutName.c_str(), sDrSeparationScint, mRMaxFiber.back() + sEpsilon, dzFibers / 2 + sEpsilon);

new TGeoConeSeg(fiberConeCutName.c_str(), sDzContainerCone / 2 + sEpsilon, 0,

sDrMinContainerCone + sXYThicknessContainerCone + sEpsilon, 0, sDrMinContainerFront + sEpsilon, -90, 90);

new TGeoTube(fiberHoleCutName.c_str(), 0, mRMinScintillator.front(), dzFibers / 2 + sEpsilon);

createAndRegisterTrans(fiberTransName, sXScintillator, 0, sZFiber);

createAndRegisterTrans(fiberConeCutTransName, sXScintillator, 0, sZCone);

createAndRegisterTrans(fiberHoleCutTransName, sXScintillator + sXShiftInnerRadiusScintillator, 0, sZFiber);

for (int i = 0; i < mRMinFiber.size(); ++i) {

const std::string fiberShapeName = fiberName + std::to_string(i + 1);

new TGeoTubeSeg(fiberShapeName.c_str(), mRMinFiber[i], mRMaxFiber[i] - sEpsilon, dzFibers / 2, -90, 90);

// Composite shape

std::string boolFormula = "";

boolFormula += fiberShapeName + ":" + fiberTransName;

boolFormula += "-" + fiberSepCutName + ":" + fiberTransName;

boolFormula += "-" + fiberConeCutName + ":" + fiberConeCutTransName;

if (i == 0) {

// Cut out the hole in the innermost fiber volume

boolFormula += "-" + fiberHoleCutName + ":" + fiberHoleCutTransName;

}

// Remove holes for screws and rods

boolFormula += "-" + sScrewHolesCSName;

boolFormula += "-" + sRodHolesCSName;

const TGeoCompositeShape* fiberCS = new TGeoCompositeShape((fiberShapeName + "CS").c_str(), boolFormula.c_str());

new TGeoVolume(createVolumeName(sFiberName, i + 1).c_str(), fiberCS, mMediumFiberRings[i]);

}

// Volume for fibers in front of PMTs

for (int i = 0; i < sNumberOfPMTFiberVolumes; i++) {

new TGeoVolume(createVolumeName(sFiberName + sPmtName, i + 1).c_str(),

new TGeoTube(createVolumeName(sFiberName + sPmtName + "Shape", i + 1).c_str(), 0, sDrPmt / 2, sDzPmt / 2),

mMediumFiberPMTs[i]);

}

}

void Geometry::initializeScrews()

{

for (int i = 0; i < sNumberOfScrewTypes; ++i) {

const std::string screwName = createVolumeName(sScrewName, i);

const TGeoShape* screwShape = createScrewShape(screwName + "Shape", i, 0, 0, 0);

// If modifying materials, make sure the appropriate initialization is done in initializeXxxMedium() methods

new TGeoVolume(screwName.c_str(), screwShape, mMediumScrewTypes[i]);

}

}

void Geometry::initializeRods()

{

for (int i = 0; i < sNumberOfRodTypes; ++i) {

const std::string rodName = createVolumeName(sRodName, i);

const TGeoShape* rodShape = createRodShape(rodName + "Shape", i, -sEpsilon, -sEpsilon);

// If modifying materials, make sure the appropriate initialization is done in initializeXxxMedium() methods

new TGeoVolume(rodName.c_str(), rodShape, mMediumRodTypes[i]);

}

}

void Geometry::initializeMetalContainer()

{

// The metal container is constructed starting from the backplate. The backplate is positioned first, relative to

// the scintillator cells. The rest of the container parts are positioned relative to the backplate.

//

// Caution: some position variables are in global coords, and some are relative to some other part of the container

// Backplate

const std::string backPlateName = "FV0_BackPlate"; // the full backplate

const std::string backPlateStandName = backPlateName + "Stand"; // the stand part of the backplate

const std::string backPlateHoleName = backPlateName + "Hole"; // the hole in the middle of the backplate

const std::string backPlateHoleCutName = backPlateHoleName + "Cut"; // extension of the hole

const std::string backPlateStandTransName = backPlateStandName + "Trans"; // shift of the backplate stand

const std::string backPlateHoleTransName = backPlateHoleName + "Trans"; // shift of the backplate inner radius

new TGeoTubeSeg(backPlateName.c_str(), 0, sDrMaxContainerBack, sDzContainerBack / 2, -90, 90);

new TGeoBBox(backPlateStandName.c_str(), sDxContainerStand / 2, (sDrMaxContainerBack + sDyContainerStand) / 2,

sDzContainerBack / 2);

new TGeoTubeSeg(backPlateHoleName.c_str(), 0, sDrContainerHole, sDzContainerBack / 2, -90, 90);

new TGeoBBox(backPlateHoleCutName.c_str(), -sXShiftContainerHole, sDrContainerHole, sDzContainerBack);

createAndRegisterTrans(backPlateStandTransName, sDxContainerStand / 2,

-(sDrMaxContainerBack + sDyContainerStand) / 2, 0);

createAndRegisterTrans(backPlateHoleTransName, sXShiftContainerHole, 0, 0);

// Backplate composite shape

std::string backPlateBoolFormula = "";

backPlateBoolFormula += backPlateName;

backPlateBoolFormula += "+" + backPlateStandName + ":" + backPlateStandTransName;

backPlateBoolFormula += "-" + backPlateHoleName + ":" + backPlateHoleTransName;

backPlateBoolFormula += "-" + backPlateHoleCutName;

const std::string backPlateCSName = backPlateName + "CS";

const std::string backPlateCSTransName = backPlateCSName + "Trans";

new TGeoCompositeShape(backPlateCSName.c_str(), backPlateBoolFormula.c_str());

createAndRegisterTrans(backPlateCSTransName, 0, 0, sZContainerBack);

// Frontplate

// the z-position o the frontplate

const float zPosFrontPlate = sZContainerFront;

// the height of the total stand overlapping with the rest of the plate

const float dyFrontPlateStand = sDyContainerStand + (sDrMaxContainerFront - sDrMinContainerFront) / 2;

// the y-position of the total stand

const float yPosFrontPlateStand = -sDrMaxContainerFront - sDyContainerStand + dyFrontPlateStand / 2;

const std::string frontPlateName = "FV0_FrontPlate";

const std::string frontPlateStandName = frontPlateName + "Stand";

const std::string frontPlateTransName = frontPlateName + "Trans";

const std::string frontPlateStandTransName = frontPlateStandName + "Trans";

new TGeoTubeSeg(frontPlateName.c_str(), sDrMinContainerFront, sDrMaxContainerFront, sDzContainerFront / 2, -90, 90);

new TGeoBBox(frontPlateStandName.c_str(), sDxContainerStand / 2, dyFrontPlateStand / 2, sDzContainerBack / 2);

createAndRegisterTrans(frontPlateTransName, 0, 0, zPosFrontPlate);

createAndRegisterTrans(frontPlateStandTransName, sDxContainerStand / 2, yPosFrontPlateStand, 0);

// Frontplate cone composite shape

std::string frontPlateBoolFormula = "";

frontPlateBoolFormula += frontPlateName;

frontPlateBoolFormula += "+" + frontPlateStandName + ":" + frontPlateStandTransName;

const std::string frontPlateCSName = frontPlateName + "CS";

new TGeoCompositeShape(frontPlateCSName.c_str(), frontPlateBoolFormula.c_str());

// Frontplate cone

// radial thickness of frontplate cone in the xy-plane

const float thicknessFrontPlateCone = sXYThicknessContainerCone;

// z-position of the frontplate cone relative to the frontplate

const float zPosCone = sDzContainerFront / 2 - sDzContainerCone / 2;

const std::string frontPlateConeName = "FV0_FrontPlateCone"; // no volume with this name

const std::string frontPlateConeShieldName = frontPlateConeName + "Shield"; // the "sides" of the cone

const std::string frontPlateConeShieldTransName = frontPlateConeShieldName + "Trans";

new TGeoConeSeg(frontPlateConeShieldName.c_str(), sDzContainerCone / 2, sDrMinContainerCone,

sDrMinContainerCone + thicknessFrontPlateCone, sDrMinContainerFront - thicknessFrontPlateCone,

sDrMinContainerFront,

-90, 90);

createAndRegisterTrans(frontPlateConeShieldTransName, 0, 0, zPosCone);

// Frontplate cone "bottom"

// z-position of the cone bottom relative to the frontplate

const float zPosConePlate = sDzContainerFront / 2 - sDzContainerCone + thicknessFrontPlateCone / 2;

// the bottom of the cone

const std::string frontPlateConePlateName = frontPlateConeName + "Plate";

new TGeoTubeSeg(frontPlateConePlateName.c_str(), 0, sDrMinContainerCone + thicknessFrontPlateCone,

thicknessFrontPlateCone / 2, -90, 90);

// Frontplate cone bottom composite shape

std::string frontPlateConePlateCSBoolFormula;

frontPlateConePlateCSBoolFormula += frontPlateConePlateName;

frontPlateConePlateCSBoolFormula += "-" + backPlateHoleName + ":" + backPlateHoleTransName;

const std::string frontPlateConePlateCSName = frontPlateConePlateName + "CS";

const std::string frontPlateConePlateCSTransName = frontPlateConePlateCSName + "Trans";

new TGeoCompositeShape(frontPlateConePlateCSName.c_str(), frontPlateConePlateCSBoolFormula.c_str());

createAndRegisterTrans(frontPlateConePlateCSTransName, 0, 0, zPosConePlate);

// Frontplate cone composite shape

std::string frontPlateConeCSBoolFormula = "";

frontPlateConeCSBoolFormula += frontPlateConeShieldName + ":" + frontPlateConeShieldTransName;

frontPlateConeCSBoolFormula += "+" + frontPlateConePlateCSName + ":" + frontPlateConePlateCSTransName;

const std::string frontPlateConeCSName = frontPlateConeName + "CS";

new TGeoCompositeShape(frontPlateConeCSName.c_str(), frontPlateConeCSBoolFormula.c_str());

// Shields

const float dzShieldGap = 0.7; // z-distance between the shields and the front- and backplate outer edges (in z-direction)

const float dzShield = sDzContainer - 2 * dzShieldGap; // depth of the shields

// Outer shield

const float zPosOuterShield = (sZContainerBack + sZContainerFront) / 2; // z-position of the outer shield

const std::string outerShieldName = "FV0_OuterShield";

const std::string outerShieldTransName = outerShieldName + "Trans";

new TGeoTubeSeg(outerShieldName.c_str(), sDrMinContainerOuterShield, sDrMaxContainerOuterShield, dzShield / 2, -90,

90);

createAndRegisterTrans(outerShieldTransName, 0, 0, zPosOuterShield);

// Inner shield

const float dzInnerShield = sDzContainer - sDzContainerCone - dzShieldGap; // depth of the inner shield

const float zPosInnerShield = sZContainerBack - sDzContainerBack / 2 + dzShieldGap + dzInnerShield / 2; // z-position of the inner shield relative to the backplate

const std::string innerShieldName = "FV0_InnerShield";

const std::string innerShieldCutName = innerShieldName + "Cut";

new TGeoTubeSeg(innerShieldName.c_str(), sDrMinContainerInnerShield, sDrMaxContainerInnerShield, dzInnerShield / 2, -90, 90);

new TGeoBBox(innerShieldCutName.c_str(), fabs(sXShiftContainerHole), sDrMaxContainerInnerShield, dzInnerShield / 2);

// Inner shield composite shape

std::string innerShieldCSBoolFormula;

innerShieldCSBoolFormula = innerShieldName;

innerShieldCSBoolFormula += "-" + innerShieldCutName;

const std::string innerShieldCSName = innerShieldName + "CS";

const std::string innerShieldCSTransName = innerShieldCSName + "Trans";

new TGeoCompositeShape(innerShieldCSName.c_str(), innerShieldCSBoolFormula.c_str());

createAndRegisterTrans(innerShieldCSTransName, sXShiftContainerHole, 0, zPosInnerShield);

// Cover

const float dzCover = sDzContainer; // Depth of the covers

const float zPosCoverConeCut = zPosFrontPlate + zPosCone; // Set the cone cut relative to the frontplate so that the exact position of the aluminium cone part can be used.

const std::string coverName = "FV0_Cover";

const std::string coverConeCutName = coverName + "ConeCut";

const std::string coverHoleCutName = coverName + "HoleCut";

new TGeoBBox(coverName.c_str(), sDxContainerCover / 2, sDrMaxContainerOuterShield, dzCover / 2);

new TGeoCone(coverConeCutName.c_str(), sDzContainerCone / 2, 0, sDrMinContainerCone + thicknessFrontPlateCone, 0,

sDrMinContainerFront);

new TGeoTubeSeg(coverHoleCutName.c_str(), 0, sDrMinContainerInnerShield, dzCover / 2, 0, 360);

const std::string coverTransName = coverName + "Trans";

const std::string coverConeCutTransName = coverConeCutName + "Trans";

const std::string coverHoleCutTransName = coverHoleCutName + "Trans";

createAndRegisterTrans(coverTransName, sDxContainerCover / 2, 0, zPosOuterShield);

createAndRegisterTrans(coverConeCutTransName, 0, 0, zPosCoverConeCut);

createAndRegisterTrans(coverHoleCutTransName.c_str(), sXShiftContainerHole, 0, zPosOuterShield);

// Cover composite shape

std::string coverCSBoolFormula = "";

coverCSBoolFormula += coverName + ":" + coverTransName;

coverCSBoolFormula += "-" + coverConeCutName + ":" + coverConeCutTransName;

coverCSBoolFormula += "-" + coverHoleCutName + ":" + coverHoleCutTransName;

const std::string coverCSName = coverName + "CS";

new TGeoCompositeShape(coverCSName.c_str(), coverCSBoolFormula.c_str());

// Stand bottom

const float dzStandBottom = sDzContainer - sDzContainerBack - sDzContainerFront;

const float dyStandBottomGap = 0.5; // This bottom part is not vertically aligned with the "front and backplate stands"

const float dxStandBottomHole = 9.4;

const float dzStandBottomHole = 20.4;

const float dxStandBottomHoleSpacing = 3.1;

const std::string standName = "FV0_StandBottom";

const std::string standHoleName = standName + "Hole";

new TGeoBBox(standName.c_str(), sDxContainerStandBottom / 2, sDyContainerStandBottom / 2, dzStandBottom / 2);

new TGeoBBox(standHoleName.c_str(), dxStandBottomHole / 2, sDyContainerStandBottom / 2 + sEpsilon,

dzStandBottomHole / 2);

const std::string standHoleTrans1Name = standHoleName + "Trans1";

const std::string standHoleTrans2Name = standHoleName + "Trans2";

const std::string standHoleTrans3Name = standHoleName + "Trans3";

createAndRegisterTrans(standHoleTrans1Name, -dxStandBottomHoleSpacing - dxStandBottomHole, 0, 0);

createAndRegisterTrans(standHoleTrans2Name, 0, 0, 0);

createAndRegisterTrans(standHoleTrans3Name, dxStandBottomHoleSpacing + dxStandBottomHole, 0, 0);

// Stand bottom composite shape