// 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.

#include <iomanip>

#include <TGeoBBox.h>

#include <TGeoManager.h>

#include <TGeoMatrix.h>

#include <TList.h>

#include <FairLogger.h>

#include "EMCALBase/Geometry.h"

#include "EMCALBase/ShishKebabTrd1Module.h"

#include <boost/algorithm/string/predicate.hpp>

using namespace o2::emcal;

// these initialisations are needed for a singleton

Geometry* Geometry::sGeom = nullptr;

Geometry::Geometry(const Geometry& geo)

: mGeoName(geo.mGeoName),

mKey110DEG(geo.mKey110DEG),

mnSupModInDCAL(geo.mnSupModInDCAL),

mNCellsInSupMod(geo.mNCellsInSupMod),

mNETAdiv(geo.mNETAdiv),

mNPHIdiv(geo.mNPHIdiv),

mNCellsInModule(geo.mNCellsInModule),

mPhiBoundariesOfSM(geo.mPhiBoundariesOfSM),

mPhiCentersOfSM(geo.mPhiCentersOfSM),

mPhiCentersOfSMSec(geo.mPhiCentersOfSMSec),

mPhiCentersOfCells(geo.mPhiCentersOfCells),

mCentersOfCellsEtaDir(geo.mCentersOfCellsEtaDir),

mCentersOfCellsPhiDir(geo.mCentersOfCellsPhiDir),

mEtaCentersOfCells(geo.mEtaCentersOfCells),

mNCells(geo.mNCells),

mNPhi(geo.mNPhi),

mCentersOfCellsXDir(geo.mCentersOfCellsXDir),

mArm1EtaMin(geo.mArm1EtaMin),

mArm1EtaMax(geo.mArm1EtaMax),

mArm1PhiMin(geo.mArm1PhiMin),

mArm1PhiMax(geo.mArm1PhiMax),

mEtaMaxOfTRD1(geo.mEtaMaxOfTRD1),

mDCALPhiMin(geo.mDCALPhiMin),

mDCALPhiMax(geo.mDCALPhiMax),

mEMCALPhiMax(geo.mEMCALPhiMax),

mDCALStandardPhiMax(geo.mDCALStandardPhiMax),

mDCALInnerExtandedEta(geo.mDCALInnerExtandedEta),

mDCALInnerEdge(geo.mDCALInnerEdge),

mShishKebabTrd1Modules(geo.mShishKebabTrd1Modules),

mPhiModuleSize(geo.mPhiModuleSize),

mEtaModuleSize(geo.mEtaModuleSize),

mPhiTileSize(geo.mPhiTileSize),

mEtaTileSize(geo.mEtaTileSize),

mNZ(geo.mNZ),

mIPDistance(geo.mIPDistance),

mLongModuleSize(geo.mLongModuleSize),

mShellThickness(geo.mShellThickness),

mZLength(geo.mZLength),

mSampling(geo.mSampling),

mECPbRadThickness(geo.mECPbRadThickness),

mECScintThick(geo.mECScintThick),

mNECLayers(geo.mNECLayers),

mNumberOfSuperModules(geo.mNumberOfSuperModules),

mEMCSMSystem(geo.mEMCSMSystem),

mFrontSteelStrip(geo.mFrontSteelStrip),

mLateralSteelStrip(geo.mLateralSteelStrip),

mPassiveScintThick(geo.mPassiveScintThick),

mPhiSuperModule(geo.mPhiSuperModule),

mNPhiSuperModule(geo.mNPhiSuperModule),

mTrd1Angle(geo.mTrd1Angle),

m2Trd1Dx2(geo.m2Trd1Dx2),

mPhiGapForSM(geo.mPhiGapForSM),

mTrd1AlFrontThick(geo.mTrd1AlFrontThick),

mTrd1BondPaperThick(geo.mTrd1BondPaperThick),

mILOSS(geo.mILOSS),

mIHADR(geo.mIHADR),

mSteelFrontThick(geo.mSteelFrontThick), // obsolete data member?

mCellIndexLookup(geo.mCellIndexLookup)

{

memcpy(mEnvelop, geo.mEnvelop, sizeof(Float_t) * 3);

memcpy(mParSM, geo.mParSM, sizeof(Float_t) * 3);

memset(SMODULEMATRIX, 0, sizeof(TGeoHMatrix*) * EMCAL_MODULES);

}

Geometry::Geometry(const std::string_view name, const std::string_view mcname, const std::string_view mctitle)

: mGeoName(name),

mKey110DEG(0),

mnSupModInDCAL(0),

mNCellsInSupMod(0),

mNETAdiv(0),

mNPHIdiv(0),

mNCellsInModule(0),

mPhiBoundariesOfSM(),

mPhiCentersOfSM(),

mPhiCentersOfSMSec(),

mPhiCentersOfCells(),

mCentersOfCellsEtaDir(),

mCentersOfCellsPhiDir(),

mEtaCentersOfCells(),

mNCells(0),

mNPhi(0),

mCentersOfCellsXDir(),

mArm1EtaMin(0),

mArm1EtaMax(0),

mArm1PhiMin(0),

mArm1PhiMax(0),

mEtaMaxOfTRD1(0),

mDCALPhiMin(0),

mDCALPhiMax(0),

mEMCALPhiMax(0),

mDCALStandardPhiMax(0),

mDCALInnerExtandedEta(0),

mDCALInnerEdge(0.),

mShishKebabTrd1Modules(),

mPhiModuleSize(0.),

mEtaModuleSize(0.),

mPhiTileSize(0.),

mEtaTileSize(0.),

mNZ(0),

mIPDistance(0.),

mLongModuleSize(0.),

mShellThickness(0.),

mZLength(0.),

mSampling(0.),

mECPbRadThickness(0.),

mECScintThick(0.),

mNECLayers(0),

mNumberOfSuperModules(0),

mEMCSMSystem(),

mFrontSteelStrip(0.),

mLateralSteelStrip(0.),

mPassiveScintThick(0.),

mPhiSuperModule(0),

mNPhiSuperModule(0),

mTrd1Angle(0.),

m2Trd1Dx2(0.),

mPhiGapForSM(0.),

mTrd1AlFrontThick(0.0),

mTrd1BondPaperThick(0.),

mILOSS(-1),

mIHADR(-1),

mSteelFrontThick(0.) // obsolete data member?

{

DefineEMC(mcname, mctitle);

mNCellsInModule = mNPHIdiv * mNETAdiv;

CreateListOfTrd1Modules();

mCellIndexLookup.resize(mNCells);

for (auto icell = 0; icell < mNCells; icell++) {

mCellIndexLookup[icell] = CalculateCellIndex(icell);

}

memset(SMODULEMATRIX, 0, sizeof(TGeoHMatrix*) * EMCAL_MODULES);

LOG(DEBUG) << "Name <<" << name << ">>";

}

Geometry& Geometry::operator=(const Geometry& /*rvalue*/)

{

LOG(FATAL) << "assignment operator, not implemented";

return *this;

}

Geometry::~Geometry()

{

if (this == sGeom) {

LOG(ERROR) << "Do not call delete on me";

return;

}

for (Int_t smod = 0; smod < mNumberOfSuperModules; smod++) {

if (SMODULEMATRIX[smod]) {

delete SMODULEMATRIX[smod];

}

}

}

Geometry* Geometry::GetInstance()

{

Geometry* rv = static_cast<Geometry*>(sGeom);

return rv;

}

Geometry* Geometry::GetInstance(const std::string_view name, const std::string_view mcname,

const std::string_view mctitle)

{

if (!sGeom) {

if (!name.length()) { // get default geometry

sGeom = new Geometry(DEFAULT_GEOMETRY, mcname, mctitle);

} else {

sGeom = new Geometry(name, mcname, mctitle);

} // end if strcmp(name,"")

return sGeom;

} else {

if (sGeom->GetName() != name) {

LOG(INFO) << "\n current geometry is " << sGeom->GetName() << " : you should not call " << name;

} // end

return sGeom;

} // end if sGeom

return nullptr;

}

Geometry* Geometry::GetInstanceFromRunNumber(Int_t runNumber, const std::string_view geoName,

const std::string_view mcname, const std::string_view mctitle)

{

using boost::algorithm::contains;

// printf("AliEMCALGeometry::GetInstanceFromRunNumber() - run %d, geoName <<%s>> \n",runNumber,geoName.Data());

if (runNumber >= 104064 && runNumber < 140000) {

// 2009-2010 runs

// First year geometry, 4 SM.

if (contains(geoName, "FIRSTYEARV1") && geoName != std::string("")) {

LOG(INFO) << "o2::emcal::Geometry::GetInstanceFromRunNumber() *** ATTENTION *** \n"

<< "\t Specified geometry name <<" << geoName << ">> for run " << runNumber

<< " is not considered! \n"

<< "\t In use <<EMCAL_FIRSTYEARV1>>, check run number and year";

} else {

LOG(INFO)

<< "o2::emcal::Geometry::GetInstanceFromRunNumber() - Initialized geometry with name <<EMCAL_FIRSTYEARV1>>";

}

return Geometry::GetInstance("EMCAL_FIRSTYEARV1", mcname, mctitle);

} else if (runNumber >= 140000 && runNumber <= 170593) {

// Almost complete EMCAL geometry, 10 SM. Year 2011 configuration

if (contains(geoName, "COMPLETEV1") && geoName != std::string("")) {

LOG(INFO) << "o2::emcal::Geometry::GetInstanceFromRunNumber() *** ATTENTION *** \n"

<< "\t Specified geometry name <<" << geoName << ">> for run " << runNumber

<< " is not considered! \n"

<< "\t In use <<EMCAL_COMPLETEV1>>, check run number and year";

} else {

LOG(INFO)

<< "o2::emcal::Geometry::GetInstanceFromRunNumber() - Initialized geometry with name <<EMCAL_COMPLETEV1>>";

}

return Geometry::GetInstance("EMCAL_COMPLETEV1", mcname, mctitle);

} else if (runNumber > 176000 && runNumber <= 197692) {

// Complete EMCAL geometry, 12 SM. Year 2012 and on

// The last 2 SM were not active, anyway they were there.

if (contains(geoName, "COMPLETE12SMV1") && geoName != std::string("")) {

LOG(INFO) << "o2::emcal::Geometry::GetInstanceFromRunNumber() *** ATTENTION *** \n"

<< "\t Specified geometry name <<" << geoName << " >> for run " << runNumber

<< " is not considered! \n"

<< "\t In use <<EMCAL_COMPLETE12SMV1>>, check run number and year";

} else {

LOG(INFO) << "o2::emcal::Geometry::GetInstanceFromRunNumber() - Initialized geometry with name "

"<<EMCAL_COMPLETE12SMV1>>";

}

return Geometry::GetInstance("EMCAL_COMPLETE12SMV1", mcname, mctitle);

} else // Run 2

{

// EMCAL + DCAL geometry, 20 SM. Year 2015 and on

if (contains(geoName, "DCAL_8SM") && geoName != std::string("")) {

LOG(INFO) << "o2::emcal::Geometry::GetInstanceFromRunNumber() *** ATTENTION *** \n"

<< "\t Specified geometry name <<" << geoName << ">> for run " << runNumber

<< " is not considered! \n"

<< "\t In use <<EMCAL_COMPLETE12SMV1_DCAL_8SM>>, check run number and year";

} else {

LOG(INFO) << "o2::emcal::Geometry::GetInstanceFromRunNumber() - Initialized geometry with name "

"<<EMCAL_COMPLETE12SMV1_DCAL_8SM>>";

}

return Geometry::GetInstance("EMCAL_COMPLETE12SMV1_DCAL_8SM", mcname, mctitle);

}

}

void Geometry::DefineSamplingFraction(const std::string_view mcname, const std::string_view mctitle)

{

// Jun 05,2006

// Look http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html

// Keep for compatibility

//

using boost::algorithm::contains;

// Sampling factor for G3

mSampling = 10.87; // Default value - Nov 25,2010

if (mNECLayers == 69) { // 10% layer reduction

mSampling = 12.55;

} else if (mNECLayers == 61) { // 20% layer reduction

mSampling = 12.80;

} else if (mNECLayers == 77) {

if (contains(mGeoName, "V1")) {

mSampling = 10.87; // Adding paper sheets and cover plate; Nov 25,2010

} else if (mECScintThick > 0.159 && mECScintThick < 0.161) { // original sampling fraction, equal layers

mSampling = 12.327; // fECScintThick = fECPbRadThickness = 0.160;

} else if (mECScintThick > 0.175 && mECScintThick < 0.177) { // 10% Pb thicknes reduction

mSampling = 10.5; // fECScintThick = 0.176, fECPbRadThickness=0.144;

} else if (mECScintThick > 0.191 && mECScintThick < 0.193) { // 20% Pb thicknes reduction

mSampling = 8.93; // fECScintThick = 0.192, fECPbRadThickness=0.128;

}

}

Float_t samplingFactorTranportModel = 1.;

if (contains(mcname, "Geant3")) {

samplingFactorTranportModel = 1.; // 0.988 // Do nothing

} else if (contains(mcname, "Fluka")) {

samplingFactorTranportModel = 1.; // To be set

} else if (contains(mcname, "Geant4")) {

if (contains(mctitle, "EMV-EMCAL")) {

samplingFactorTranportModel = 0.821; // EMC list but for EMCal, before 0.86

} else if (contains(mctitle, "EMV")) {

samplingFactorTranportModel = 1.096; // 0.906, 0.896 (OPT)

} else {

samplingFactorTranportModel = 0.821; // 1.15 (CHIPS), 1.149 (BERT), 1.147 (BERT_CHIPS)

}

}

LOG(INFO) << "MC modeler <" << mcname << ">, Title <" << mctitle << ">: Sampling " << std::setw(2)

<< std::setprecision(3) << mSampling << ", model fraction with respect to G3 "

<< samplingFactorTranportModel << ", final sampling " << mSampling * samplingFactorTranportModel;

mSampling *= samplingFactorTranportModel;

}

void Geometry::DefineEMC(std::string_view mcname, std::string_view mctitle)

{

using boost::algorithm::contains;

// geometry

std::transform(mGeoName.begin(), mGeoName.end(), mGeoName.begin(), ::toupper);

// Convert old geometry names to new ones

if (contains(mGeoName, "SHISH_77_TRD1_2X2_FINAL_110DEG")) {

if (contains(mGeoName, "PBTH=0.144") && contains(mGeoName, "SCTH=0.176")) {

mGeoName = "EMCAL_COMPLETE";

} else {

mGeoName = "EMCAL_PDC06";

}

}

if (contains(mGeoName, "WSUC")) {

mGeoName = "EMCAL_WSUC";

}

// check that we have a valid geometry name

if (!(contains(mGeoName, "EMCAL_PDC06") || contains(mGeoName, "EMCAL_WSUC") || contains(mGeoName, "EMCAL_COMPLETE") ||

contains(mGeoName, "EMCAL_COMPLETEV1") || contains(mGeoName, "EMCAL_COMPLETE12SMV1") ||

contains(mGeoName, "EMCAL_FIRSTYEAR") || contains(mGeoName, "EMCAL_FIRSTYEARV1"))) {

LOG(FATAL) << "Init, " << mGeoName << " is an undefined geometry!\n";

}

// Option to know whether we have the "half" supermodule(s) or not

mKey110DEG = 0;

if (contains(mGeoName, "COMPLETE") || contains(mGeoName, "PDC06") || contains(mGeoName, "12SM")) {

mKey110DEG = 1; // for GetAbsCellId

}

if (contains(mGeoName, "COMPLETEV1")) {

mKey110DEG = 0;

}

mnSupModInDCAL = 0;

if (contains(mGeoName, "DCAL_DEV")) {

mnSupModInDCAL = 10;

} else if (contains(mGeoName, "DCAL_8SM")) {

mnSupModInDCAL = 8;

} else if (contains(mGeoName, "DCAL")) {

mnSupModInDCAL = 6;

}

// JLK 13-Apr-2008

// default parameters are those of EMCAL_COMPLETE geometry

// all others render variations from these at the end of

// geometry-name specific options

mNumberOfSuperModules = 12; // 12 = 6 * 2 (6 in phi, 2 in Z)

mNPhi = 12; // module granularity in phi within smod (azimuth)

mNZ = 24; // module granularity along Z within smod (eta)

mNPHIdiv = mNETAdiv = 2; // tower granularity within module

mArm1PhiMin = 80.0; // degrees, Starting EMCAL Phi position

mArm1PhiMax = 200.0; // degrees, Ending EMCAL Phi position

mArm1EtaMin = -0.7; // pseudorapidity, Starting EMCAL Eta position

mArm1EtaMax = +0.7; // pseudorapidity, Ending EMCAL Eta position

mIPDistance = 428.0; // cm, radial distance to front face from nominal vertex point

mPhiGapForSM = 2.; // cm, only for final TRD1 geometry

mFrontSteelStrip = 0.025; // 0.025cm = 0.25mm (13-may-05 from V.Petrov)

mPassiveScintThick = 0.8; // 0.8cm = 8mm (13-may-05 from V.Petrov)

mLateralSteelStrip = 0.01; // 0.01cm = 0.1mm (13-may-05 from V.Petrov) - was 0.025

mTrd1Angle = 1.5; // in degrees

mSampling = 1.; // should be calculated with call to DefineSamplingFraction()

mNECLayers = 77; // (13-may-05 from V.Petrov) - can be changed with additional options

mECScintThick = 0.176; // scintillator layer thickness

mECPbRadThickness = 0.144; // lead layer thickness

mPhiModuleSize = 12.26 - mPhiGapForSM / Float_t(mNPhi); // first assumption

mEtaModuleSize = mPhiModuleSize;

mZLength = 700.; // Z coverage (cm)

mPhiSuperModule = 20.; // phi in degree

mDCALInnerEdge = mIPDistance * TMath::Tan(mTrd1Angle * 8. * TMath::DegToRad());

// modifications to the above for PDC06 geometry

if (contains(mGeoName, "PDC06")) { // 18-may-05 - about common structure

mECScintThick = mECPbRadThickness = 0.16; // (13-may-05 from V.Petrov)

}

// modifications to the above for WSUC geometry

if (contains(mGeoName, "WSUC")) { // 18-may-05 - about common structure

mNumberOfSuperModules = 2; // 27-may-05; Nov 24,2010 for TB

mNPhi = mNZ = 4;

mTrd1AlFrontThick = 1.0; // one cm

// Bond paper - two sheets around Sc tile

mTrd1BondPaperThick = 0.01; // 0.01cm = 0.1 mm

mPhiModuleSize = 12.0;

mEtaModuleSize = mPhiModuleSize;

mLateralSteelStrip = 0.015; // 0.015cm = 0.15mm

}

// In 2009-2010 data taking runs only 4 SM, in the upper position.

if (contains(mGeoName, "FIRSTYEAR")) {

mNumberOfSuperModules = 4;

mArm1PhiMax = 120.0;

}

if (contains(mGeoName, "FIRSTYEARV1") || contains(mGeoName, "COMPLETEV1") || contains(mGeoName, "COMPLETE12SMV1")) {

// Oct 26,2010 : First module has tilt = 0.75 degree :

// look to AliEMCALShishKebabTrd1Module::DefineFirstModule(key)

// New sizes from production drawing, added Al front plate.

// The thickness of sampling is change due to existing two sheets of paper.

// Will replace fFrontSteelStrip

mTrd1AlFrontThick = 1.0; // one cm

// Bond paper - two sheets around Sc tile

mTrd1BondPaperThick = 0.01; // 0.01cm = 0.1 mm

mPhiModuleSize = 12.0;

mEtaModuleSize = mPhiModuleSize;

mLateralSteelStrip = 0.015; // 0.015cm = 0.15mm

if (contains(mGeoName, "COMPLETEV1")) {

mNumberOfSuperModules = 10;

mArm1PhiMax = 180.0;

} else if (contains(mGeoName, "COMPLETE12SMV1")) {

mNumberOfSuperModules = 12;

mArm1PhiMax = 200.0;

}

if (contains(mGeoName, "DCAL")) {

mNumberOfSuperModules = 12 + mnSupModInDCAL;

mArm1PhiMax = 320.0;

if (contains(mGeoName, "DCAL_8SM")) {

mArm1PhiMax = 340.0; // degrees, End of DCAL Phi position

} else if (contains(mGeoName, "DCAL_DEV")) {

mArm1PhiMin = 40.0; // degrees, Starting EMCAL(shifted) Phi position

}

mDCALPhiMin = mArm1PhiMax - 10. * mnSupModInDCAL;

}

}

//

// Init EMCal/DCal SMs type array

mEMCSMSystem.clear();

mEMCSMSystem.resize(mNumberOfSuperModules);

for (Int_t i = 0; i < mNumberOfSuperModules; i++) {

mEMCSMSystem[i] = NOT_EXISTENT;

}

Int_t iSM = 0;

//

// BASIC EMCAL SM

if (contains(mGeoName, "WSUC")) {

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

mEMCSMSystem[iSM] = EMCAL_STANDARD;

iSM++;

}

} else if (contains(mGeoName, "FIRSTYEAR")) {

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

mEMCSMSystem[iSM] = EMCAL_STANDARD;

iSM++;

}

} else if (contains(mGeoName, "PDC06") || contains(mGeoName, "COMPLETE")) {

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

mEMCSMSystem[iSM] = EMCAL_STANDARD;

iSM++;

}

}

//

// EMCAL 110SM

if (mKey110DEG && contains(mGeoName, "12SM")) {

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

mEMCSMSystem[iSM] = EMCAL_HALF;

if (contains(mGeoName, "12SMV1")) {

mEMCSMSystem[iSM] = EMCAL_THIRD;

}

iSM++;

}

}

//

// DCAL SM

if (mnSupModInDCAL && contains(mGeoName, "DCAL")) {

if (contains(mGeoName, "8SM")) {

for (int i = 0; i < mnSupModInDCAL - 2; i++) {

mEMCSMSystem[iSM] = DCAL_STANDARD;

iSM++;

}

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

mEMCSMSystem[iSM] = DCAL_EXT;

iSM++;

}

} else {

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

mEMCSMSystem[iSM] = DCAL_STANDARD;

iSM++;

}

}

}

// constant for transition absid <--> indexes

mNCellsInModule = mNPHIdiv * mNETAdiv;

mNCellsInSupMod = mNCellsInModule * mNPhi * mNZ;

mNCells = 0;

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

if (GetSMType(i) == EMCAL_STANDARD) {

mNCells += mNCellsInSupMod;

} else if (GetSMType(i) == EMCAL_HALF) {

mNCells += mNCellsInSupMod / 2;

} else if (GetSMType(i) == EMCAL_THIRD) {

mNCells += mNCellsInSupMod / 3;

} else if (GetSMType(i) == DCAL_STANDARD) {

mNCells += 2 * mNCellsInSupMod / 3;

} else if (GetSMType(i) == DCAL_EXT) {

mNCells += mNCellsInSupMod / 3;

} else {

LOG(ERROR) << "Uknown SuperModule Type !!\n";

}

}

mNPhiSuperModule = mNumberOfSuperModules / 2;

if (mNPhiSuperModule < 1) {

mNPhiSuperModule = 1;

}

mPhiTileSize = mPhiModuleSize / double(mNPHIdiv) - mLateralSteelStrip; // 13-may-05

mEtaTileSize = mEtaModuleSize / double(mNETAdiv) - mLateralSteelStrip; // 13-may-05

mLongModuleSize = mNECLayers * (mECScintThick + mECPbRadThickness);

if (contains(mGeoName, "V1")) {

Double_t ws = mECScintThick + mECPbRadThickness + 2. * mTrd1BondPaperThick; // sampling width

// Number of Pb tiles = Number of Sc tiles - 1

mLongModuleSize = mTrd1AlFrontThick + (ws * mNECLayers - mECPbRadThickness);

}

m2Trd1Dx2 = mEtaModuleSize + 2. * mLongModuleSize * TMath::Tan(mTrd1Angle * TMath::DegToRad() / 2.);

if (!contains(mGeoName, "WSUC")) {

mShellThickness = TMath::Sqrt(mLongModuleSize * mLongModuleSize + m2Trd1Dx2 * m2Trd1Dx2);

}

// These parameters are used to create the mother volume to hold the supermodules

// 2cm padding added to allow for misalignments - JLK 30-May-2008

mEnvelop[0] = mIPDistance - 1.; // mother volume inner radius

mEnvelop[1] = mIPDistance + mShellThickness + 1.; // mother volume outer r.

mEnvelop[2] = mZLength + 2.; // mother volume length

// Local coordinates

mParSM[0] = GetShellThickness() / 2.;

mParSM[1] = GetPhiModuleSize() * GetNPhi() / 2.;

mParSM[2] = mZLength / 4.; // divide by 4 to get half-length of SM

// SM phi boundaries - (0,1),(2,3) ... - has the same boundaries;

mPhiBoundariesOfSM.resize(mNumberOfSuperModules);

mPhiCentersOfSM.resize(mNumberOfSuperModules / 2);

mPhiCentersOfSMSec.resize(mNumberOfSuperModules / 2);

Double_t kfSupermodulePhiWidth = mPhiSuperModule * TMath::DegToRad();

mPhiCentersOfSM[0] = (mArm1PhiMin + mPhiSuperModule / 2.) * TMath::DegToRad(); // Define from First SM

mPhiCentersOfSMSec[0] = mPhiCentersOfSM[0]; // the same in the First SM

mPhiBoundariesOfSM[0] = mPhiCentersOfSM[0] - TMath::ATan2(mParSM[1], mIPDistance); // 1th and 2th modules)

mPhiBoundariesOfSM[1] = mPhiCentersOfSM[0] + TMath::ATan2(mParSM[1], mIPDistance);

if (mNumberOfSuperModules > 2) { // 2 to Max

Int_t tmpSMType = GetSMType(2);

for (int i = 1; i < mNPhiSuperModule; i++) {

mPhiBoundariesOfSM[2 * i] += mPhiBoundariesOfSM[2 * i - 2] + kfSupermodulePhiWidth;

if (tmpSMType == GetSMType(2 * i)) {

mPhiBoundariesOfSM[2 * i + 1] += mPhiBoundariesOfSM[2 * i - 1] + kfSupermodulePhiWidth;

} else {

// changed SM Type, redefine the [2*i+1] Boundaries

tmpSMType = GetSMType(2 * i);

switch (GetSMType(2 * i)) {

case EMCAL_STANDARD:

mPhiBoundariesOfSM[2 * i + 1] = mPhiBoundariesOfSM[2 * i] + kfSupermodulePhiWidth;

break;

case EMCAL_HALF:

mPhiBoundariesOfSM[2 * i + 1] = mPhiBoundariesOfSM[2 * i] + 2. * TMath::ATan2((mParSM[1]) / 2, mIPDistance);

break;

case EMCAL_THIRD:

mPhiBoundariesOfSM[2 * i + 1] = mPhiBoundariesOfSM[2 * i] + 2. * TMath::ATan2((mParSM[1]) / 3, mIPDistance);

break;

case DCAL_STANDARD:

mPhiBoundariesOfSM[2 * i] = (mDCALPhiMin - mArm1PhiMin) * TMath::DegToRad() + mPhiBoundariesOfSM[0];

mPhiBoundariesOfSM[2 * i + 1] = (mDCALPhiMin - mArm1PhiMin) * TMath::DegToRad() + mPhiBoundariesOfSM[1];

break;

case DCAL_EXT:

mPhiBoundariesOfSM[2 * i + 1] = mPhiBoundariesOfSM[2 * i] + 2. * TMath::ATan2((mParSM[1]) / 3, mIPDistance);

break;

default:

break;

};

}

mPhiCentersOfSM[i] = (mPhiBoundariesOfSM[2 * i] + mPhiBoundariesOfSM[2 * i + 1]) / 2.;

mPhiCentersOfSMSec[i] = mPhiBoundariesOfSM[2 * i] + TMath::ATan2(mParSM[1], mIPDistance);

}

}

// inner extend in eta (same as outer part) for DCal (0.189917), //calculated from the smallest gap (1# cell to the

// 80-degree-edge),

const double INNNER_EXTENDED_PHI =

1.102840997; // calculated from the smallest gap (1# cell to the 80-degree-edge), too complicatd to explain...

mDCALInnerExtandedEta = -TMath::Log(

TMath::Tan((TMath::Pi() / 2. - 8 * mTrd1Angle * TMath::DegToRad() +

(TMath::Pi() / 2 - mNZ * mTrd1Angle * TMath::DegToRad() - TMath::ATan(TMath::Exp(mArm1EtaMin)) * 2)) /

2.));

mEMCALPhiMax = mArm1PhiMin;

mDCALPhiMax = mDCALPhiMin; // DCAl extention will not be included

for (Int_t i = 0; i < mNumberOfSuperModules; i += 2) {

switch (GetSMType(i)) {

case EMCAL_STANDARD:

mEMCALPhiMax += 20.;

break;

case EMCAL_HALF:

mEMCALPhiMax += mPhiSuperModule / 2. + INNNER_EXTENDED_PHI;

break;

case EMCAL_THIRD:

mEMCALPhiMax += mPhiSuperModule / 3. + 4.0 * INNNER_EXTENDED_PHI / 3.0;

break;

case DCAL_STANDARD:

mDCALPhiMax += 20.;

mDCALStandardPhiMax = mDCALPhiMax;

break;

case DCAL_EXT:

mDCALPhiMax += mPhiSuperModule / 3. + 4.0 * INNNER_EXTENDED_PHI / 3.0;

break;

default:

LOG(ERROR) << "Unkown SM Type!!\n";

break;

};

}

// for compatible reason

// if(fNumberOfSuperModules == 4) {fEMCALPhiMax = fArm1PhiMax ;}

if (mNumberOfSuperModules == 12) {

mEMCALPhiMax = mArm1PhiMax;

}

// called after setting of scintillator and lead layer parameters

// called now in AliEMCALv0::CreateGeometry() - 15/03/16

// DefineSamplingFraction(mcname,mctitle);

}

void Geometry::GetGlobal(const Double_t* loc, Double_t* glob, int iSM) const

{

const TGeoHMatrix* m = GetMatrixForSuperModule(iSM);

if (m) {

m->LocalToMaster(loc, glob);

} else {

LOG(FATAL) << "Geo matrixes are not loaded \n";

}

}

void Geometry::GetGlobal(const TVector3& vloc, TVector3& vglob, int iSM) const

{

Double_t tglob[3], tloc[3];

vloc.GetXYZ(tloc);

GetGlobal(tloc, tglob, iSM);

vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);

}

void Geometry::GetGlobal(Int_t absId, Double_t glob[3]) const

{

double loc[3];

memset(glob, 0, sizeof(Double_t) * 3);

try {

auto cellpos = RelPosCellInSModule(absId);

loc[0] = cellpos.X();

loc[1] = cellpos.Y();

loc[2] = cellpos.Z();

} catch (InvalidCellIDException& e) {

LOG(ERROR) << e.what();

return;

}

Int_t nSupMod = std::get<0>(GetCellIndex(absId));

const TGeoHMatrix* m = GetMatrixForSuperModule(nSupMod);

if (m) {

m->LocalToMaster(loc, glob);

} else {

LOG(FATAL) << "Geo matrixes are not loaded \n";

}

}

void Geometry::GetGlobal(Int_t absId, TVector3& vglob) const

{

Double_t glob[3];

GetGlobal(absId, glob);

vglob.SetXYZ(glob[0], glob[1], glob[2]);

}

std::tuple<double, double> Geometry::EtaPhiFromIndex(Int_t absId) const

{

TVector3 vglob;

GetGlobal(absId, vglob);

return std::make_tuple(vglob.Eta(), vglob.Phi());

}

int Geometry::GetAbsCellId(int supermoduleID, int moduleID, int phiInModule, int etaInModule) const

{

// 0 <= nSupMod < fNumberOfSuperModules

// 0 <= nModule < fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)

// 0 <= nIphi < fNPHIdiv

// 0 <= nIeta < fNETAdiv

// 0 <= absid < fNCells

int cellid = 0; // have to change from 0 to fNCells-1

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

switch (GetSMType(i)) {

case EMCAL_STANDARD:

cellid += mNCellsInSupMod;

break;

case EMCAL_HALF:

cellid += mNCellsInSupMod / 2;

break;

case EMCAL_THIRD:

cellid += mNCellsInSupMod / 3;

break;

case DCAL_STANDARD:

cellid += 2 * mNCellsInSupMod / 3;

break;

case DCAL_EXT:

cellid += mNCellsInSupMod / 3;

break;

default:

throw InvalidSupermoduleTypeException();

};

}

cellid += mNCellsInModule * moduleID;

cellid += mNPHIdiv * phiInModule;

cellid += etaInModule;

if (!CheckAbsCellId(cellid)) {

throw InvalidCellIDException(cellid);

}

return cellid;

}

std::tuple<int, int, int> Geometry::GetModuleIndexesFromCellIndexesInSModule(int supermoduleID, int phiInSupermodule, int etaInSupermodule) const

{

int nModulesInSMPhi = GetNumberOfModuleInPhiDirection(supermoduleID);

int moduleEta = etaInSupermodule / mNETAdiv,

modulePhi = phiInSupermodule / mNPHIdiv,

moduleID = moduleEta * nModulesInSMPhi + modulePhi;

int etaInModule = etaInSupermodule % mNETAdiv,

phiInModule = phiInSupermodule % mNPHIdiv;

//return std::make_tuple(modulePhi, moduleEta, moduleID);

return std::make_tuple(phiInModule, etaInModule, moduleID);

}

Int_t Geometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const

{

// Check if the indeces correspond to existing SM or tower indeces

if (iphi < 0 || iphi >= EMCAL_ROWS || ieta < 0 || ieta >= EMCAL_COLS || nSupMod < 0 ||

nSupMod >= GetNumberOfSuperModules()) {

LOG(DEBUG) << "Wrong cell indexes : SM " << nSupMod << ", column (eta) " << ieta << ", row (phi) " << iphi;

return -1;

}

auto indexmod = GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta);

Int_t nIeta = ieta % mNETAdiv,

nIphi = iphi % mNPHIdiv;

nIeta = mNETAdiv - 1 - nIeta;

return GetAbsCellId(nSupMod, std::get<2>(indexmod), nIphi, nIeta);

}

std::tuple<int, int> Geometry::GlobalRowColFromIndex(int cellID) const

{

if (!CheckAbsCellId(cellID)) {

throw InvalidCellIDException(cellID);

}

auto [supermodule, module, phiInModule, etaInModule] = GetCellIndex(cellID);

auto [row, col] = GetCellPhiEtaIndexInSModule(supermodule, module, phiInModule, etaInModule);

// add offsets (row / col per supermodule)

if (supermodule == 13 || supermodule == 15 || supermodule == 17) {

// DCal odd SMs need shift of the col. index in oder to get the global col. index

col += 16;

}

if (supermodule % 2) {

col += mNZ * 2;

}

int sector = supermodule / 2;

if (sector > 0) {

for (int isec = 0; isec < sector; isec++) {

auto smtype = GetSMType(isec * 2);

auto nphism = (smtype == EMCAL_THIRD || smtype == DCAL_EXT) ? GetNPhi() / 3 : GetNPhi();

row += 2 * nphism;

}

}

return std::make_tuple(row, col);

}

std::tuple<int, int, int> Geometry::GetPositionInSupermoduleFromGlobalRowCol(int row, int col) const

{

if (col < 0 || col >= 4 * GetNEta()) {

throw RowColException(row, col);

}

int side = col < GetNEta() * 2 ? 0 : 1,

colSM = col % (GetNEta() * 2);

int sector = -1,

rowSM = row;

for (int isec = 0; isec < GetNPhiSuperModule(); isec++) {

auto smtype = GetSMType(isec * 2);

auto nphism = GetNPhi() * 2;

if (smtype == EMCAL_THIRD || smtype == DCAL_EXT) {

nphism /= 3;

}

if (rowSM < nphism) {

sector = isec;

break;

}

rowSM -= nphism;

}

if (sector < 0) {

throw RowColException(row, col);

}

int supermodule = sector * 2 + side;

if (supermodule == 13 || supermodule == 15 || supermodule == 17) {

// DCal odd SMs need shift of the col. index as global col index includes PHOS hole

colSM -= 16;

if (colSM < 0) {

throw RowColException(row, col); // Position inside PHOS hole specified

}

}

if (supermodule == 12 || supermodule == 14 || supermodule == 16) {

if (colSM > 32) {

throw RowColException(row, col); // Position inside PHOS hole specified

}

}

return std::make_tuple(supermodule, rowSM, colSM);

}

int Geometry::GetCellAbsIDFromGlobalRowCol(int row, int col) const

{

auto [supermodule, rowSM, colSM] = GetPositionInSupermoduleFromGlobalRowCol(row, col);

return GetAbsCellIdFromCellIndexes(supermodule, rowSM, colSM);

}

std::tuple<int, int, int, int> Geometry::GetCellIndexFromGlobalRowCol(int row, int col) const

{

auto [supermodule, rowSM, colSM] = GetPositionInSupermoduleFromGlobalRowCol(row, col);

auto indexmod = GetModuleIndexesFromCellIndexesInSModule(supermodule, rowSM, colSM);

Int_t colInModule = colSM % mNETAdiv,

rowInMOdule = rowSM % mNPHIdiv;

colInModule = mNETAdiv - 1 - colInModule;

return std::make_tuple(supermodule, std::get<2>(indexmod), rowInMOdule, colInModule);

}

int Geometry::GlobalCol(int cellID) const

{

return std::get<1>(GlobalRowColFromIndex(cellID));

}

int Geometry::GlobalRow(int cellID) const

{

return std::get<0>(GlobalRowColFromIndex(cellID));

}

Int_t Geometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi) const

{

if (TMath::Abs(eta) > mEtaMaxOfTRD1) {

throw InvalidPositionException(eta, phi);

}

phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries

Int_t nphism = mNumberOfSuperModules / 2;

Int_t nSupMod = 0;

for (Int_t i = 0; i < nphism; i++) {

LOG(DEBUG) << "Sec " << i << ": Min " << mPhiBoundariesOfSM[2 * i] << ", Max " << mPhiBoundariesOfSM[2 * i + 1];

if (phi >= mPhiBoundariesOfSM[2 * i] && phi <= mPhiBoundariesOfSM[2 * i + 1]) {

nSupMod = 2 * i;

if (eta < 0.0) {

nSupMod++;

}

if (GetSMType(nSupMod) == DCAL_STANDARD) { // Gap between DCAL

if (TMath::Abs(eta) < GetNEta() / 3 * mTrd1Angle * TMath::DegToRad()) {

throw InvalidPositionException(eta, phi);

}

}

LOG(DEBUG) << "eta " << eta << " phi " << phi << " (" << std::setw(5) << std::setprecision(2)

<< phi * TMath::RadToDeg() << ") : nSupMod " << nSupMod << ": #bound " << i;

return nSupMod;

}

}

throw InvalidPositionException(eta, phi);

}

Int_t Geometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi) const

{

Int_t nSupMod = SuperModuleNumberFromEtaPhi(eta, phi);

// phi index first

phi = TVector2::Phi_0_2pi(phi);

Double_t phiLoc = phi - mPhiCentersOfSMSec[nSupMod / 2];

Int_t nphi = mPhiCentersOfCells.size();

switch (GetSMType(nSupMod)) {

case EMCAL_HALF:

nphi /= 2;

case EMCAL_THIRD:

case DCAL_EXT:

nphi /= 3;

break;

default:

// All other supermodules have full number of cells in phi

break;

};

Double_t dmin = TMath::Abs(mPhiCentersOfCells[0] - phiLoc),

d = 0.;

Int_t iphi = 0;

for (Int_t i = 1; i < nphi; i++) {

d = TMath::Abs(mPhiCentersOfCells[i] - phiLoc);

if (d < dmin) {

dmin = d;

iphi = i;

}

// printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);

}

// odd SM are turned with respect of even SM - reverse indexes

LOG(DEBUG2) << " iphi " << iphi << " : dmin " << dmin << " (phi " << phi << ", phiLoc " << phiLoc << ")\n";

// eta index

Double_t absEta = TMath::Abs(eta);

Int_t neta = mCentersOfCellsEtaDir.size(),

etaShift = iphi * neta,

ieta = 0;

if (GetSMType(nSupMod) == DCAL_STANDARD) {

ieta += 16; // jump 16 cells for DCSM

}

dmin = TMath::Abs(mEtaCentersOfCells[etaShift + ieta] - absEta);

for (Int_t i = ieta + 1; i < neta; i++) {

d = TMath::Abs(mEtaCentersOfCells[i + etaShift] - absEta);

if (d < dmin) {

dmin = d;

ieta = i;

}

}

if (GetSMType(nSupMod) == DCAL_STANDARD) {

ieta -= 16; // jump 16 cells for DCSM

}

LOG(DEBUG2) << " ieta " << ieta << " : dmin " << dmin << " (eta=" << eta << ") : nSupMod " << nSupMod;

// patch for mapping following alice convention

if (nSupMod % 2 ==

0) { // 47 + 16 -ieta for DCSM, 47 - ieta for others, revert the ordering on A side in order to keep convention.

ieta = (neta - 1) - ieta;

if (GetSMType(nSupMod) == DCAL_STANDARD) {

ieta -= 16; // recover cells for DCSM

}

}

return GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);

}

std::tuple<int, int, int, int> Geometry::CalculateCellIndex(Int_t absId) const

{

if (!CheckAbsCellId(absId)) {

throw InvalidCellIDException(absId);

}

Int_t tmp = absId;

Int_t test = absId;

Int_t nSupMod;

for (nSupMod = -1; test >= 0;) {

nSupMod++;

tmp = test;