// 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 "DetectorsBase/MaterialManager.h"

#include "TPCSimulation/Detector.h"

#include "TPCSimulation/Point.h"

#include "TPCBase/ParameterGas.h"

#include "TPCBase/ParameterDetector.h"

#include "DetectorsBase/Stack.h"

#include "SimulationDataFormat/TrackReference.h"

#include "FairVolume.h" // for FairVolume

#include "TVirtualMC.h" // for TVirtualMC, gMC

#include <cstddef> // for NULL

#include "FairGeoVolume.h"

#include "FairGeoNode.h"

#include "FairGeoLoader.h"

#include "FairGeoInterface.h"

#include "FairRun.h"

#include "FairRuntimeDb.h"

#include <fairlogger/Logger.h>

#include "FairRootManager.h"

#include "TSystem.h"

#include "TVirtualMC.h"

#include "TFile.h"

// geo stuff

#include "TGeoManager.h"

#include "TGeoVolume.h"

#include "TGeoPcon.h"

#include "TGeoTube.h"

#include "TGeoCone.h"

#include "TGeoPgon.h"

#include "TGeoTrd1.h"

#include "TGeoCompositeShape.h"

#include "TGeoPara.h"

#include "TGeoPhysicalNode.h"

#include "TGeoHalfSpace.h"

#include "TGeoArb8.h"

#include "TGeoMatrix.h"

#include <iostream>

#include <cmath>

using std::cout;

using std::endl;

using std::ifstream;

using std::ios_base;

using namespace o2::tpc;

Detector::Detector(Bool_t active) : o2::base::DetImpl<Detector>("TPC", active), mGeoFileName()

{

for (int i = 0; i < Sector::MAXSECTOR; ++i) {

mHitsPerSectorCollection[i] = o2::utils::createSimVector<o2::tpc::HitGroup>(); // new std::vector<o2::tpc::HitGroup>;

}

}

// forward default constructor

Detector::Detector() : o2::tpc::Detector(kTRUE) {}

Detector::Detector(const Detector& rhs)

: o2::base::DetImpl<Detector>(rhs),

mHitCounter(0),

mElectronCounter(0),

mStepCounter(0),

mGeoFileName(rhs.mGeoFileName)

{

for (int i = 0; i < Sector::MAXSECTOR; ++i) {

mHitsPerSectorCollection[i] = o2::utils::createSimVector<o2::tpc::HitGroup>(); // new std::vector<o2::tpc::HitGroup>;new std::vector<o2::tpc::HitGroup>;

}

}

Detector::~Detector()

{

for (int i = 0; i < Sector::MAXSECTOR; ++i) {

// mHitsPerSectorCollection[i]->clear();

o2::utils::freeSimVector(mHitsPerSectorCollection[i]);

}

std::cout << "Produced hits " << mHitCounter << "\n";

std::cout << "Produced electrons " << mElectronCounter << "\n";

std::cout << "Stepping called " << mStepCounter << "\n";

}

void Detector::InitializeO2Detector()

{

// Define the list of sensitive volumes

defineSensitiveVolumes();

}

Bool_t Detector::ProcessHits(FairVolume* vol)

{

mStepCounter++;

auto& gasParam = ParameterGas::Instance();

auto& detParam = ParameterDetector::Instance();

const Int_t kMaxDistRef = 15; // maximal difference between 2 stored references - the parameter should be 15 cm as default

static Double_t lastReferenceR = 0; // keeps last reference point in radius (cm)

/* This method is called from the MC stepping for the sensitive volume only */

// LOG(info) << "tpc::ProcessHits";

const double trackCharge = fMC->TrackCharge();

if (static_cast<int>(trackCharge) == 0) {

// set a very large step size for neutral particles

fMC->SetMaxStep(1.e10);

return kFALSE; // take only charged particles

}

// ===| SET THE LENGTH OF THE NEXT ENERGY LOSS STEP |=========================

// (We do this first so we can skip out of the method in the following

// Eloss->Ionization part.)

//

// In all cases we have multiple collisions and we use 2mm (+ some

// random shift to avoid binning effects), which was tuned for GEANT4, see

// https://indico.cern.ch/event/316891/contributions/732168/

const double rnd = fMC->GetRandom()->Rndm();

fMC->SetMaxStep(0.2 + (2. * rnd - 1.) * 0.05); // 2 mm +- rndm*0.5mm step

// ===| check active sector |=================================================

//

// Get the sector ID and check if the sector is active

static thread_local TLorentzVector position;

fMC->TrackPosition(position);

// for processing reasons in the digitizer, the sectors are shifted by -10deg, so sector 0 will be

// from -10 - +10 deg instead of 0-20 and so on

const int sectorID = static_cast<int>(Sector::ToShiftedSector(position.X(), position.Y(), position.Z()));

// const int sectorID = static_cast<int>(Sector::ToSector(position.X(), position.Y(), position.Z()));

// TODO: Temporary hack to process only one sector

// if (sectorID != 0) return kFALSE;

// ---| momentum and beta gamma |---

static TLorentzVector momentum; // static to make avoid creation/deletion of this expensive object

fMC->TrackMomentum(momentum);

const float time = fMC->TrackTime() * 1.0e9;

const int trackID = fMC->GetStack()->GetCurrentTrackNumber();

const int detID = vol->getMCid();

o2::data::Stack* stack = (o2::data::Stack*)fMC->GetStack();

if (fMC->IsTrackEntering() || fMC->IsTrackExiting()) {

stack->addTrackReference(o2::TrackReference(position.X(), position.Y(), position.Z(), momentum.X(), momentum.Y(),

momentum.Z(), fMC->TrackLength(), time, trackID, GetDetId()));

lastReferenceR = fMC->TrackLength();

}

if (TMath::Abs(lastReferenceR - fMC->TrackLength()) > kMaxDistRef) { /// we can speedup

stack->addTrackReference(o2::TrackReference(position.X(), position.Y(), position.Z(), momentum.X(), momentum.Y(),

momentum.Z(), fMC->TrackLength(), time, trackID, GetDetId()));

lastReferenceR = fMC->TrackLength();

}

// ---| remove clusters between the IFC and the FC strips |---

// those should not enter the active readout area

// do coarse selection before, to limit number of transformations

if (detParam.ExcludeFCGap) {

const auto rCluster = std::sqrt(position.X() * position.X() + position.Y() * position.Y());

const float rodRin = 81.5 + 2.2; // radial position of the inner field cage rods + radial size of the field cage rods

const float rodRout = 254.25 + 2.2; // radial position of the outer field cage rods + radial size of the field cage rods

const float fcLxIn = 82.428409; // position of the inner FC strips in local x = cos(10 deg) * rodRin;

const float fcLxOut = 252.55395; // position of the outer FC strips in local x = cos(10 deg) * rodRin;

if (rCluster < rodRin || rCluster > fcLxOut) {

const int sectorIDnonShift = static_cast<int>(Sector::ToSector(position.X(), position.Y(), position.Z()));

const double alpha = TMath::DegToRad() * (10. + sectorIDnonShift * 20.);

const double cs = std::cos(-alpha), sn = std::sin(-alpha);

const auto localX = position.X() * cs - position.Y() * sn;

// fine cut

if (localX < fcLxIn || localX > fcLxOut) {

return kFALSE;

}

}

}

// ===| CONVERT THE ENERGY LOSS TO IONIZATION ELECTRONS |=====================

//

// The energy loss is implemented directly below and taken GEANT3,

// ILOSS model 5 (in gfluct.F), which gives

// the energy loss in a single collision (NA49 model).

// TODO: Add discussion about drawback

Int_t numberOfElectrons = 0;

// I.H. - the type expected in addHit is short

// ---| Stepsize in cm |---

const double stepSize = fMC->TrackStep();

double betaGamma = momentum.P() / fMC->TrackMass();

betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg

// ---| number of primary ionisations per cm |---

const double primaryElectronsPerCM =

gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],

gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],

gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);

// ---| mean number of collisions and random for this event |---

const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM;

const int nColl = static_cast<int>(fMC->GetRandom()->Poisson(meanNcoll));

// Variables needed to generate random powerlaw distributed energy loss

const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value

const double oneOverAlpha_p1 = 1. / alpha_p1;

const double eMin = gasParam.Ipot;

const double eMax = gasParam.Eend;

const double kMin = TMath::Power(eMin, alpha_p1);

const double kMax = TMath::Power(eMax, alpha_p1);

const double wIon = gasParam.Wion;

for (Int_t n = 0; n < nColl; n++) {

// Use GEANT3 / NA49 expression:

// P(eDep) ~ k * edep^-gasParam.getExp()

// eMin(~I) < eDep < eMax(300 electrons)

// k fixed so that Int_Emin^EMax P(Edep) = 1.

const double rndm = fMC->GetRandom()->Rndm();

const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1);

int nel_step = static_cast<int>(((eDep - eMin) / wIon) + 1);

nel_step = TMath::Min(nel_step, 300); // 300 electrons corresponds to 10 keV

numberOfElectrons += nel_step;

}

// LOG(info) << "tpc::AddHit" << FairLogger::endl << "Eloss: "

//<< fMC->Edep() << ", Nelectrons: "

//<< numberOfElectrons;

if (numberOfElectrons <= 0) { // Could maybe be smaller than 0 due to the Gamma function

return kFALSE;

}

// ADD HIT

static thread_local int oldTrackId = trackID;

static thread_local int oldDetId = detID;

static thread_local int groupCounter = 0;

static thread_local int oldSectorId = sectorID;

// a new group is starting -> put it into the container

static thread_local HitGroup* currentgroup = nullptr;

if (groupCounter == 0) {

mHitsPerSectorCollection[sectorID]->emplace_back(trackID);

currentgroup = &(mHitsPerSectorCollection[sectorID]->back());

}

if (trackID == oldTrackId && oldSectorId == sectorID) {

groupCounter++;

mHitCounter++;

mElectronCounter += numberOfElectrons;

currentgroup->addHit(position.X(), position.Y(), position.Z(), time, numberOfElectrons);

// add last buffered hit, which was not yet added to the currentgroup

if (mHitLast.GetEnergyLoss() >= 0) {

currentgroup->addHit(mHitLast.GetX(), mHitLast.GetY(), mHitLast.GetZ(), mHitLast.GetTime(), mHitLast.GetEnergyLoss());

mHitLast.mELoss = -1;

groupCounter++;

mHitCounter++;

mElectronCounter += mHitLast.GetEnergyLoss();

}

}

// finish group

else {

oldTrackId = trackID;

oldSectorId = sectorID;

groupCounter = 0;

// buffer this hit, otherwise it wouldnt be stored in the HitGroup

mHitLast = ElementalHit(position.X(), position.Y(), position.Z(), time, numberOfElectrons);

}

// LOG(info) << "tpc::AddHit" << FairLogger::endl

//<< " -- " << trackNumberID <<"," << volumeID << " " << vol->GetName()

//<< ", Pos: (" << position.X() << ", " << position.Y() <<", "<< position.Z()<< ", " << r << ") "

//<< ", Mom: (" << momentum.Px() << ", " << momentum.Py() << ", " << momentum.Pz() << ") "

//<< " Time: "<< time <<", Len: " << length << ", Nelectrons: " <<

// numberOfElectrons;

// I.H. - the code above does not compile if uncommented

// Increment number of Detector det points in TParticle

stack->addHit(GetDetId());

return kTRUE;

}

void Detector::EndOfEvent()

{

if (!o2::utils::ShmManager::Instance().isOperational()) {

for (int i = 0; i < Sector::MAXSECTOR; ++i) {

mHitsPerSectorCollection[i]->clear();

}

}

}

void Detector::Register()

{

/** This will create a branch in the output tree called

DetectorPoint, setting the last parameter to kFALSE means:

this collection will not be written to the file, it will exist

only during the simulation.

*/

auto* mgr = FairRootManager::Instance();

for (int i = 0; i < Sector::MAXSECTOR; ++i) {

mgr->RegisterAny(getHitBranchNames(i).c_str(), mHitsPerSectorCollection[i], kTRUE);

}

}

void Detector::Reset()

{

for (int i = 0; i < Sector::MAXSECTOR; ++i) {

// cout << "CLEARED (reset) hitsCollection " << i << " " << mHitsPerSectorCollection[i] << endl;

mHitsPerSectorCollection[i]->clear();

}

}

void Detector::ConstructGeometry()

{

// Create the detector materials

CreateMaterials();

// Load geometry

// LoadGeometryFromFile();

ConstructTPCGeometry();

}

void Detector::CreateMaterials()

{

//-----------------------------------------------

// Create Materials for for TPC simulations

//-----------------------------------------------

//-----------------------------------------------------------------

// Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl

//-----------------------------------------------------------------

const auto& gasParam = ParameterGas::Instance();

Int_t iSXFLD = 2;

Float_t sXMGMX = 10.0;

// init the field tracking params

o2::base::Detector::initFieldTrackingParams(iSXFLD, sXMGMX);

Float_t amat[7]; // atomic numbers

Float_t zmat[7]; // z

Float_t wmat[7]; // proportions

Float_t density;

// TODO: load pressure and temperature values from CCDB

const Double_t pressure = gasParam.Pressure; // in mbar

const Double_t temperature = gasParam.Temperature + 273.15; // in K

// densities were taken for these values

const Double_t t1 = 293.15; // 20°C in K

const Double_t p1 = 1013.25; // 1 atm in mbars

// sanity check - temperature between 10 and 30 deg, pressure between 800 and 1200 mbar

Double_t ptCorr = 1.;

if (TMath::Abs(temperature - 293.15) > 10. || TMath::Abs(pressure - 1000.) > 200.) {

ptCorr = 1.;

} else {

ptCorr = (pressure * t1) / (p1 * temperature);

}

LOG(info) << "Setting gas density correction to: " << ptCorr;

//***************** Gases *************************

//--------------------------------------------------------------

// gases - air and CO2

//--------------------------------------------------------------

// CO2

amat[0] = 12.011;

amat[1] = 15.9994;

zmat[0] = 6.;

zmat[1] = 8.;

wmat[0] = 0.2729;

wmat[1] = 0.7271;

density = 1.842e-3;

o2::base::Detector::Mixture(10, "CO2", amat, zmat, density * ptCorr, 2, wmat);

//

// Air

//

amat[0] = 15.9994;

amat[1] = 14.007;

//

zmat[0] = 8.;

zmat[1] = 7.;

//

wmat[0] = 0.233;

wmat[1] = 0.767;

//

density = 0.001205;

o2::base::Detector::Mixture(11, "Air", amat, zmat, density * ptCorr, 2, wmat);

//----------------------------------------------------------------

// drift gases 5 mixtures, 5 materials

//----------------------------------------------------------------

//

// Drift gases 1 - nonsensitive, 2 - sensitive, 3 - for Kr

// Composition by % of volume, values at 20deg and 1 atm.

//

// get the geometry title - defined in Config.C

//

//--------------------------------------------------------------

// predefined gases, composition taken from param file

//--------------------------------------------------------------

TString names[6] = {"Ne", "Ar", "CO2", "N", "CF4", "CH4"};

TString gname;

/// @todo: Gas mixture is hard coded here, this should be moved to some kind of parameter

// container in the future

Float_t comp[6] = {90. / 105., 0., 10. / 105., 5. / 105., 0., 0.};

// indices:

// 0-Ne, 1-Ar, 2-CO2, 3-N, 4-CF4, 5-CH4

//

// elements' masses

//

amat[0] = 20.18; // Ne

amat[1] = 39.95; // Ar

amat[2] = 12.011; // C

amat[3] = 15.9994; // O

amat[4] = 14.007; // N

amat[5] = 18.998; // F

amat[6] = 1.; // H

//

// elements' atomic numbers

//

//

zmat[0] = 10.; // Ne

zmat[1] = 18.; // Ar

zmat[2] = 6.; // C

zmat[3] = 8.; // O

zmat[4] = 7.; // N

zmat[5] = 9.; // F

zmat[6] = 1.; // H

//

// Mol masses

//

Float_t wmol[6];

wmol[0] = 20.18; // Ne

wmol[1] = 39.948; // Ar

wmol[2] = 44.0098; // CO2

wmol[3] = 2. * 14.0067; // N2

wmol[4] = 88.0046; // CF4

wmol[5] = 16.011; // CH4

//

Float_t wtot = 0.; // total mass of the mixture

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

wtot += *(comp + i) * wmol[i];

}

wmat[0] = comp[0] * amat[0] / wtot; // Ne

wmat[1] = comp[1] * amat[1] / wtot; // Ar

wmat[2] = (comp[2] * amat[2] + comp[4] * amat[2] + comp[5] * amat[2]) / wtot; // C

wmat[3] = comp[2] * amat[3] * 2. / wtot; // O

wmat[4] = comp[3] * amat[4] * 2. / wtot; // N

wmat[5] = comp[4] * amat[5] * 4. / wtot; // F

wmat[6] = comp[5] * amat[6] * 4. / wtot; // H

//

// densities (NTP)

//

Float_t dens[6] = {0.839e-3, 1.661e-3, 1.842e-3, 1.165e-3, 3.466e-3, 0.668e-3};

//

density = 0.;

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

density += comp[i] * dens[i];

}

//

// names

//

Int_t cnt = 0;

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

if (comp[i]) {

if (cnt) {

gname += "-";

}

gname += names[i];

cnt++;

}

}

TString gname1, gname2, gname3;

gname1 = gname + "-1";

gname2 = gname + "-2";

gname3 = gname + "-3";

//

// take only elements with nonzero weights

//

Float_t amat1[6], zmat1[6], wmat1[6];

cnt = 0;

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

if (wmat[i]) {

zmat1[cnt] = zmat[i];

amat1[cnt] = amat[i];

wmat1[cnt] = wmat[i];

cnt++;

}

}

//

o2::base::Detector::Mixture(12, gname1.Data(), amat1, zmat1, density * ptCorr, cnt, wmat1); // nonsensitive

o2::base::Detector::Mixture(13, gname2.Data(), amat1, zmat1, density * ptCorr, cnt, wmat1); // sensitive

o2::base::Detector::Mixture(40, gname3.Data(), amat1, zmat1, density * ptCorr, cnt, wmat1); // sensitive Kr

//----------------------------------------------------------------------

// solid materials

//----------------------------------------------------------------------

// Kevlar C14H22O2N2

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.999;

amat[3] = 14.006;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

zmat[3] = 7.;

wmat[0] = 14.;

wmat[1] = 22.;

wmat[2] = 2.;

wmat[3] = 2.;

density = 1.45;

o2::base::Detector::Mixture(14, "Kevlar", amat, zmat, density, -4, wmat);

// NOMEX

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.999;

amat[3] = 14.006;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

zmat[3] = 7.;

wmat[0] = 14.;

wmat[1] = 22.;

wmat[2] = 2.;

wmat[3] = 2.;

density = 0.029;

o2::base::Detector::Mixture(15, "NOMEX", amat, zmat, density, -4, wmat);

// Makrolon C16H18O3

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.999;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 16.;

wmat[1] = 18.;

wmat[2] = 3.;

density = 1.2;

o2::base::Detector::Mixture(16, "Makrolon", amat, zmat, density, -3, wmat);

// Tedlar C2H3F

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 18.998;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 9.;

wmat[0] = 2.;

wmat[1] = 3.;

wmat[2] = 1.;

density = 1.71;

o2::base::Detector::Mixture(17, "Tedlar", amat, zmat, density, -3, wmat);

// Mylar C5H4O2

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.9994;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 5.;

wmat[1] = 4.;

wmat[2] = 2.;

density = 1.39;

o2::base::Detector::Mixture(18, "Mylar", amat, zmat, density, -3, wmat);

// material for "prepregs"

// Epoxy - C14 H20 O3

// Quartz SiO2

// Carbon C

// prepreg1 60% C-fiber, 40% epoxy (vol)

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.994;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 0.923;

wmat[1] = 0.023;

wmat[2] = 0.054;

density = 1.859;

o2::base::Detector::Mixture(19, "Prepreg1", amat, zmat, density, 3, wmat);

// prepreg2 60% glass-fiber, 40% epoxy

amat[0] = 12.01;

amat[1] = 1.;

amat[2] = 15.994;

amat[3] = 28.086;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

zmat[3] = 14.;

wmat[0] = 0.194;

wmat[1] = 0.023;

wmat[2] = 0.443;

wmat[3] = 0.34;

density = 1.82;

o2::base::Detector::Mixture(20, "Prepreg2", amat, zmat, density, 4, wmat);

// prepreg3 50% glass-fiber, 50% epoxy

amat[0] = 12.01;

amat[1] = 1.;

amat[2] = 15.994;

amat[3] = 28.086;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

zmat[3] = 14.;

wmat[0] = 0.257;

wmat[1] = 0.03;

wmat[2] = 0.412;

wmat[3] = 0.3;

density = 1.725;

o2::base::Detector::Mixture(21, "Prepreg3", amat, zmat, density, 4, wmat);

// G10 60% SiO2 40% epoxy

amat[0] = 12.01;

amat[1] = 1.;

amat[2] = 15.994;

amat[3] = 28.086;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

zmat[3] = 14.;

wmat[0] = 0.194;

wmat[1] = 0.023;

wmat[2] = 0.443;

wmat[3] = 0.340;

density = 1.7;

o2::base::Detector::Mixture(22, "G10", amat, zmat, density, 4, wmat);

// Al

amat[0] = 26.98;

zmat[0] = 13.;

density = 2.7;

o2::base::Detector::Material(23, "Al", amat[0], zmat[0], density, 999., 999.);

// Si (for electronics

amat[0] = 28.086;

zmat[0] = 14.;

density = 2.33;

o2::base::Detector::Material(24, "Si", amat[0], zmat[0], density, 999., 999.);

// Cu

amat[0] = 63.546;

zmat[0] = 29.;

density = 8.96;

o2::base::Detector::Material(25, "Cu", amat[0], zmat[0], density, 999., 999.);

// brass

amat[0] = 63.546;

zmat[0] = 29.;

//

amat[1] = 65.409;

zmat[1] = 30.;

//

wmat[0] = 0.6;

wmat[1] = 0.4;

//

density = 8.23;

//

o2::base::Detector::Mixture(33, "Brass", amat, zmat, density, 2, wmat);

// Epoxy - C14 H20 O3

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.9994;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 14.;

wmat[1] = 20.;

wmat[2] = 3.;

density = 1.25;

o2::base::Detector::Mixture(26, "Epoxy", amat, zmat, density, -3, wmat);

// Epoxy - C14 H20 O3 for glue

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.9994;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 14.;

wmat[1] = 20.;

wmat[2] = 3.;

density = 1.25;

density *= 1.25;

o2::base::Detector::Mixture(35, "Epoxy1", amat, zmat, density, -3, wmat);

//

// epoxy film - 90% epoxy, 10% glass fiber

//

amat[0] = 12.01;

amat[1] = 1.;

amat[2] = 15.994;

amat[3] = 28.086;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

zmat[3] = 14.;

wmat[0] = 0.596;

wmat[1] = 0.071;

wmat[2] = 0.257;

wmat[3] = 0.076;

density = 1.345;

o2::base::Detector::Mixture(34, "Epoxy-film", amat, zmat, density, 4, wmat);

// Plexiglas C5H8O2

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.9994;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 5.;

wmat[1] = 8.;

wmat[2] = 2.;

density = 1.18;

o2::base::Detector::Mixture(27, "Plexiglas", amat, zmat, density, -3, wmat);

// Carbon

amat[0] = 12.011;

zmat[0] = 6.;

density = 2.265;

o2::base::Detector::Material(28, "C", amat[0], zmat[0], density, 999., 999.);

// Fe (steel for the inner heat screen)

amat[0] = 55.845;

zmat[0] = 26.;

density = 7.87;

o2::base::Detector::Material(29, "Fe", amat[0], zmat[0], density, 999., 999.);

//

// Peek - (C6H4-O-OC6H4-O-C6H4-CO)n

amat[0] = 12.011;

amat[1] = 1.;

amat[2] = 15.9994;

zmat[0] = 6.;

zmat[1] = 1.;

zmat[2] = 8.;

wmat[0] = 19.;

wmat[1] = 12.;

wmat[2] = 3.;

//

density = 1.3;

//

o2::base::Detector::Mixture(30, "Peek", amat, zmat, density, -3, wmat);

//

// Ceramics - Al2O3

//

amat[0] = 26.98;

amat[1] = 15.9994;

zmat[0] = 13.;

zmat[1] = 8.;

wmat[0] = 2.;

wmat[1] = 3.;

density = 3.97;

o2::base::Detector::Mixture(31, "Alumina", amat, zmat, density, -2, wmat);

//

// Ceramics for resistors

//

amat[0] = 26.98;

amat[1] = 15.9994;

zmat[0] = 13.;

zmat[1] = 8.;

wmat[0] = 2.;

wmat[1] = 3.;

density = 3.97;

//

density *= 1.25;

o2::base::Detector::Mixture(36, "Alumina1", amat, zmat, density, -2, wmat);

//

// liquids

//

// water

amat[0] = 1.;

amat[1] = 15.9994;

zmat[0] = 1.;

zmat[1] = 8.;

wmat[0] = 2.;

wmat[1] = 1.;

density = 1.;

o2::base::Detector::Mixture(32, "Water", amat, zmat, density, -2, wmat);

//----------------------------------------------------------

// tracking media for gases

//----------------------------------------------------------

o2::base::Detector::Medium(kAir, "Air", 11, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);

o2::base::Detector::Medium(kDriftGas1, "DriftGas1", 12, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kDriftGas2, "DriftGas2", 13, 1, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kCO2, "CO2", 10, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kDriftGas3, "DriftGas3", 40, 1, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

//-----------------------------------------------------------

// tracking media for solids

//-----------------------------------------------------------

o2::base::Detector::Medium(kAl, "Al", 23, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kKevlar, "Kevlar", 14, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kNomex, "Nomex", 15, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kMakrolon, "Makrolon", 16, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kMylar, "Mylar", 18, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kTedlar, "Tedlar", 17, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

//

o2::base::Detector::Medium(kPrepreg1, "Prepreg1", 19, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kPrepreg2, "Prepreg2", 20, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kPrepreg3, "Prepreg3", 21, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kEpoxy, "Epoxy", 26, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kCu, "Cu", 25, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kSi, "Si", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kG10, "G10", 22, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kPlexiglas, "Plexiglas", 27, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kSteel, "Steel", 29, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kPeek, "Peek", 30, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kAlumina, "Alumina", 31, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kWater, "Water", 32, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kBrass, "Brass", 33, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

o2::base::Detector::Medium(kEpoxyfm, "Epoxyfm", 34, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kEpoxy1, "Epoxy1", 35, 0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);

o2::base::Detector::Medium(kAlumina1, "Alumina1", 36, 0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);

}

void Detector::ConstructTPCGeometry()

{

//

// Create the geometry of Time Projection Chamber version 2

//

// Begin_Html

/*

* <img src="picts/AliTPC.gif">

*/

// End_Html

// Begin_Html

/*

* <img src="picts/AliTPCv2Tree.gif">

*/

// End_Html

//----------------------------------------------------------

// This geometry is written using TGeo class

// Firstly the shapes are defined, and only then the volumes

// What is recognized by the MC are volumes

//----------------------------------------------------------

//

// tpc - this will be the mother volume

//

// if (!mParam) {

// LOG(error) << "TPC Parameters not available, cannot create Geometry";

// return;

// }

//

// here I define a volume TPC

// retrive the medium name with "TPC_" as a leading string

//

auto* tpc = new TGeoPcon(0., 360., 30); // 30 sections

//

tpc->DefineSection(0, -289.6, 77., 278.);

tpc->DefineSection(1, -262.1, 77., 278.);

//

tpc->DefineSection(2, -262.1, 83.1, 278.);

tpc->DefineSection(3, -260., 83.1, 278.);

//

tpc->DefineSection(4, -260., 70., 278.);

tpc->DefineSection(5, -259.6, 70., 278.);

//

tpc->DefineSection(6, -259.6, 68.1, 278.);

tpc->DefineSection(7, -253.6, 68.1, 278.);

//

tpc->DefineSection(8, -253.6, 67.88, 278.); // hs

tpc->DefineSection(9, -74.0, 60.68, 278.); // hs

//

tpc->DefineSection(10, -74.0, 60.1, 278.);

tpc->DefineSection(11, -73.3, 60.1, 278.);

//

tpc->DefineSection(12, -73.3, 56.9, 278.);

tpc->DefineSection(13, -68.5, 56.9, 278.);