17#include <FairVolume.h>
18#include <TVirtualMC.h>
19#include <TVirtualMCStack.h>
20#include <TGeoManager.h>
21#include <TGeoVolume.h>
25#include <Math/Point2D.h>
26#include <Math/Vector2D.h>
58 defineSamplingFactor();
62void Detector::defineSamplingFactor()
64 TString mcname = TVirtualMC::GetMC()->GetName();
65 TString mctitle = TVirtualMC::GetMC()->GetTitle();
66 LOGP(info,
"Defining sampling factor for mc={}' and title='{}'", mcname.Data(), mctitle.Data());
67 if (mcname.Contains(
"Geant3")) {
68 mSamplingFactorTransportModel = 0.983;
69 }
else if (mcname.Contains(
"Geant4")) {
70 mSamplingFactorTransportModel = 1.;
75void Detector::createMaterials()
77 LOGP(info,
"Creating materials for ECL");
80 float aAir[4] = {12.0107, 14.0067, 15.9994, 39.948};
81 float zAir[4] = {6., 7., 8., 18.};
82 float wAir[4] = {0.000124, 0.755267, 0.231781, 0.012827};
83 float dAir = 1.20479E-3;
84 Mixture(0,
"Air", aAir, zAir, dAir, 4, wAir);
87 Material(1,
"Pb", 207.2, 82, 11.35, 0.56, 0.,
nullptr, 0);
90 float aP[2] = {12.011, 1.00794};
91 float zP[2] = {6.0, 1.0};
92 float wP[2] = {1.0, 1.0};
94 Mixture(2,
"Scintillator", aP, zP, dP, -2, wP);
97 Material(3,
"Al", 26.98, 13., 2.7, 8.9, 999.,
nullptr, 0);
100 float aX[3] = {207.19, 183.85, 16.0};
101 float zX[3] = {82.0, 74.0, 8.0};
102 float wX[3] = {1.0, 1.0, 4.0};
104 Mixture(4,
"PbWO4", aX, zX, dX, -3, wX);
111 Medium(0,
"Air", 0, 0, ifield, fieldm, tmaxfd, 1.0, deemax, 0.1, 10.0,
nullptr, 0);
112 Medium(1,
"Pb", 1, 0, ifield, fieldm, tmaxfd, 0.1, deemax, 0.1, 0.1,
nullptr, 0);
113 Medium(2,
"Scintillator", 2, 1, ifield, fieldm, tmaxfd, 0.001, deemax, 0.001, 0.001,
nullptr, 0);
114 Medium(3,
"Al", 3, 0, ifield, fieldm, tmaxfd, 0.1, deemax, 0.001, 0.001,
nullptr, 0);
115 Medium(4,
"Crystal", 4, 1, ifield, fieldm, tmaxfd, 0.1, deemax, 0.1, 0.1,
nullptr, 0);
121 LOG(info) <<
"Initialize ECal O2Detector";
131 mSuperParentIndices.clear();
142 LOGP(info,
"Registering hits");
143 if (FairRootManager::Instance()) {
144 FairRootManager::Instance()->RegisterAny(
addNameTo(
"Hit").
data(), mHits,
true);
149void Detector::createGeometry()
151 LOGP(info,
"Creating ECal geometry");
153 TGeoVolume* vALIC = gGeoManager->GetVolume(
"barrel");
155 LOGP(fatal,
"Could not find barrel volume while constructing ECal geometry");
159 vALIC->AddNode(vECal, 2,
new TGeoTranslation(0, 30., 0));
160 vECal->SetTitle(
"ECalVol");
162 TGeoMedium* medAir = gGeoManager->GetMedium(
"ECL_Air");
163 TGeoMedium* medPb = gGeoManager->GetMedium(
"ECL_Pb");
164 TGeoMedium* medAl = gGeoManager->GetMedium(
"ECL_Al");
165 TGeoMedium* medSc = gGeoManager->GetMedium(
"ECL_Scintillator");
166 TGeoMedium* medPWO = gGeoManager->GetMedium(
"ECL_Crystal");
172 double rMin = pars.rMin;
173 double rMax = pars.rMax;
174 double layerThickness = pars.pbLayerThickness + pars.scLayerThickness;
175 double samplingModL = pars.frontPlateThickness + layerThickness * pars.nSamplingLayers - pars.pbLayerThickness;
176 double crystalAlpha = geo.getCrystalAlpha();
177 double samplingAlpha = geo.getSamplingAlpha();
178 double tanCrystalAlpha = std::tan(crystalAlpha);
179 double tanSamplingAlpha = std::tan(samplingAlpha);
181 double sectorL = rMax - rMin;
182 double crystalThetaMin = geo.getFrontFaceCenterTheta(pars.nCrystalModulesZ - 1) - crystalAlpha;
183 double crystalHlMin = geo.getFrontFaceZatMinR(pars.nCrystalModulesZ - 1);
184 double crystalHlMax = crystalHlMin + sectorL / std::tan(crystalThetaMin);
185 double crystalHwMin = geo.getCrystalModW() / 2.;
186 double crystalHwMax = crystalHwMin * rMax / rMin;
187 auto crystalSectorShape =
new TGeoTrap(sectorL / 2., 0, 0, crystalHlMin, crystalHwMin, crystalHwMin, 0, crystalHlMax, crystalHwMax, crystalHwMax, 0);
188 auto crystalSectorVolume =
new TGeoVolume(
"crystalSectorVolume", crystalSectorShape, medAir);
189 AddSensitiveVolume(crystalSectorVolume);
190 crystalSectorVolume->SetLineColor(kCyan + 1);
191 crystalSectorVolume->SetTransparency(0);
193 double samplingThetaAtMinZ = geo.getFrontFaceCenterTheta(pars.nCrystalModulesZ) + samplingAlpha;
194 double samplingThetaAtMaxZ = geo.getFrontFaceCenterTheta(pars.nCrystalModulesZ + pars.nSamplingModulesZ - 1) - samplingAlpha;
195 double samplingMinZatMinR = geo.getFrontFaceZatMinR(pars.nCrystalModulesZ) - geo.getSamplingModW() / std::sin(samplingThetaAtMinZ);
196 double samplingMaxZatMinR = geo.getFrontFaceZatMinR(pars.nCrystalModulesZ + pars.nSamplingModulesZ - 1);
197 double samplingMinZatMaxR = samplingMinZatMinR + sectorL / std::tan(samplingThetaAtMinZ);
198 double samplingMaxZatMaxR = samplingMaxZatMinR + sectorL / std::tan(samplingThetaAtMaxZ);
199 double hlMin = (samplingMaxZatMinR - samplingMinZatMinR) / 2.;
200 double hlMax = (samplingMaxZatMaxR - samplingMinZatMaxR) / 2.;
201 double zCenterMin = (samplingMaxZatMinR + samplingMinZatMinR) / 2.;
202 double zCenterMax = (samplingMaxZatMaxR + samplingMinZatMaxR) / 2.;
203 double zCenter = (zCenterMax + zCenterMin) / 2.;
204 double thetaCenter = std::atan((zCenterMax - zCenterMin) / sectorL) * TMath::RadToDeg();
205 double samplingHwMin = geo.getSamplingModW() / 2.;
206 double samplingHwMax = samplingHwMin * rMax / rMin;
207 auto samplingSectorShape =
new TGeoTrap(sectorL / 2., thetaCenter, 90, hlMin, samplingHwMin, samplingHwMin, 0, hlMax, samplingHwMax, samplingHwMax, 0);
208 auto samplingSectorVolume =
new TGeoVolume(
"samplingSectorVolume", samplingSectorShape, medAir);
209 AddSensitiveVolume(samplingSectorVolume);
210 samplingSectorVolume->SetLineColor(kBlue + 1);
211 samplingSectorVolume->SetTransparency(0);
213 double sectorR = rMin + sectorL / 2.;
214 for (
int ism = 0; ism < pars.nSuperModules; ism++) {
216 for (
int i = 0;
i < pars.nCrystalModulesPhi;
i++) {
217 int row = ism * pars.nCrystalModulesPhi +
i;
218 double phi0 = geo.getFrontFaceCenterCrystalPhi(
row);
219 double x = sectorR * std::cos(phi0);
220 double y = sectorR * std::sin(phi0);
221 auto rot =
new TGeoRotation(Form(
"ecalcrystalsecrot%d",
row), 90 + phi0 * TMath::RadToDeg(), 90, 0);
222 vECal->AddNode(crystalSectorVolume,
row,
new TGeoCombiTrans(
x,
y, 0., rot));
225 for (
int i = 0;
i < pars.nSamplingModulesPhi;
i++) {
226 int row = ism * pars.nSamplingModulesPhi +
i;
227 double phi0 = geo.getFrontFaceCenterSamplingPhi(
row);
228 double x = sectorR * std::cos(phi0);
229 double y = sectorR * std::sin(phi0);
230 auto rot1 =
new TGeoRotation(Form(
"ecalsamplingsec1rot%d",
row), 90 + phi0 * TMath::RadToDeg(), 90, 0.);
231 auto rot2 =
new TGeoRotation(Form(
"ecalsamplingsec2rot%d",
row), 90 + phi0 * TMath::RadToDeg(), 90, 180);
232 vECal->AddNode(samplingSectorVolume, 2 *
row + 0,
new TGeoCombiTrans(
x,
y, zCenter, rot1));
233 vECal->AddNode(samplingSectorVolume, 2 *
row + 1,
new TGeoCombiTrans(
x,
y, -zCenter, rot2));
237 for (
int m = 0;
m < pars.nCrystalModulesZ;
m++) {
238 double tanBeta = geo.getTanBeta(
m);
239 double dx1 = crystalHwMin;
240 double dx2 = crystalHwMin + pars.crystalModuleLength * tanCrystalAlpha;
241 double dy1 = crystalHwMin;
242 double dy2 = crystalHwMin + pars.crystalModuleLength * tanBeta;
243 double dz = pars.crystalModuleLength / 2.;
244 auto crystalModuleShape =
new TGeoTrd2(dx1, dx2, dy1, dy2, dz);
245 auto crystalModuleVolume =
new TGeoVolume(Form(
"crystalmodule%d",
m), crystalModuleShape, medPWO);
246 AddSensitiveVolume(crystalModuleVolume);
247 crystalModuleVolume->SetLineColor(kCyan + 1);
248 crystalModuleVolume->SetTransparency(0);
249 double theta = geo.getFrontFaceCenterTheta(
m);
250 double r = geo.getFrontFaceCenterR(
m);
251 double z = geo.getFrontFaceCenterZ(
m);
252 ROOT::Math::XYPoint pFrontFace(
z,
r - sectorR);
253 ROOT::Math::Polar2DVector vFrontFaceToCenter(dz, theta);
254 ROOT::Math::XYPoint pc = pFrontFace + vFrontFaceToCenter;
255 auto rot1 =
new TGeoRotation(Form(
"ecalcrystalrot%d", 2 *
m), 0, 270 + theta * TMath::RadToDeg(), 90);
256 crystalSectorVolume->AddNode(crystalModuleVolume, 2 *
m,
new TGeoCombiTrans(0, pc.x(), pc.y(), rot1));
257 auto rot2 =
new TGeoRotation(Form(
"ecalcrystalrot%d", 2 *
m + 1), 0, 90 - theta * TMath::RadToDeg(), 90);
258 crystalSectorVolume->AddNode(crystalModuleVolume, 2 *
m + 1,
new TGeoCombiTrans(0, -pc.x(), pc.y(), rot2));
261 for (
int m = 0;
m < pars.nSamplingModulesZ;
m++) {
262 int k = pars.nCrystalModulesZ +
m;
263 double tanBeta = geo.getTanBeta(k);
264 double dx1 = samplingHwMin;
265 double dx2 = samplingHwMin + samplingModL * tanSamplingAlpha;
266 double dy1 = samplingHwMin;
267 double dy2 = samplingHwMin + samplingModL * tanBeta;
268 double dz = samplingModL / 2.;
269 auto samplingModuleShape =
new TGeoTrd2(dx1, dx2, dy1, dy2, dz);
270 auto samplingModuleVolume =
new TGeoVolume(Form(
"samplingmodule%d",
m), samplingModuleShape, medSc);
271 AddSensitiveVolume(samplingModuleVolume);
272 samplingModuleVolume->SetLineColor(kAzure - 9);
273 samplingModuleVolume->SetTransparency(0);
274 double theta = geo.getFrontFaceCenterTheta(k);
275 double r = geo.getFrontFaceCenterR(k);
276 double z = geo.getFrontFaceCenterZ(k);
277 ROOT::Math::XYPoint pFrontFace(
z - zCenter,
r - sectorR);
278 ROOT::Math::Polar2DVector vFrontFaceToCenter(dz, theta);
279 ROOT::Math::XYPoint pc = pFrontFace + vFrontFaceToCenter;
280 auto rot1 =
new TGeoRotation(Form(
"ecalsamplingrot%d",
m), 0, 270 + theta * TMath::RadToDeg(), 90);
281 samplingSectorVolume->AddNode(samplingModuleVolume,
m,
new TGeoCombiTrans(0, pc.x(), pc.y(), rot1));
285 double fdx2 = dx1 + pars.frontPlateThickness * tanSamplingAlpha;
287 double fdy2 = fdy1 + pars.frontPlateThickness * tanBeta;
288 double fdz = pars.frontPlateThickness / 2.;
289 auto frontShape =
new TGeoTrd2(fdx1, fdx2, fdy1, fdy2, fdz);
290 auto frontVolume =
new TGeoVolume(Form(
"front%d",
m), frontShape, medAl);
291 samplingModuleVolume->AddNode(frontVolume, 0,
new TGeoTranslation(0., 0., -dz + pars.frontPlateThickness / 2.));
292 AddSensitiveVolume(frontVolume);
293 frontVolume->SetLineColor(kAzure - 7);
294 frontVolume->SetTransparency(0);
297 for (
int i = 0;
i < pars.nSamplingLayers - 1;
i++) {
298 double lz1 = pars.frontPlateThickness + pars.scLayerThickness + layerThickness *
i;
299 double lz2 = lz1 + pars.pbLayerThickness;
300 double lzc = -dz + (lz1 + lz2) / 2.;
301 double ldx1 = dx1 + lz1 * tanSamplingAlpha;
302 double ldx2 = dx1 + lz2 * tanSamplingAlpha;
303 double ldy1 = dy1 + lz1 * tanBeta;
304 double ldy2 = dy1 + lz2 * tanBeta;
305 double ldz = pars.pbLayerThickness / 2.;
306 auto leadShape =
new TGeoTrd2(ldx1, ldx2, ldy1, ldy2, ldz);
307 auto leadVolume =
new TGeoVolume(Form(
"lead%d_%d",
m,
i), leadShape, medPb);
308 samplingModuleVolume->AddNode(leadVolume,
i,
new TGeoTranslation(0., 0., lzc));
309 AddSensitiveVolume(leadVolume);
310 leadVolume->SetLineColor(kAzure - 7);
311 leadVolume->SetTransparency(0);
315 if (pars.enableFwdEndcap) {
317 TGeoTube* ecalEndcapShape =
new TGeoTube(
"ECLECsh", 15.f, 160.f, 0.5 * (rMax - rMin));
318 TGeoVolume* ecalEndcapVol =
new TGeoVolume(
"ECLEC", ecalEndcapShape, medPb);
319 ecalEndcapVol->SetLineColor(kAzure - 9);
320 ecalEndcapVol->SetTransparency(0);
321 vECal->AddNode(ecalEndcapVol, 1,
new TGeoTranslation(0, 0, -450.f));
330 LOGP(
debug,
"Processing hits");
332 int trackId =
stack->GetCurrentTrackNumber();
333 int parentId =
stack->GetCurrentParentTrackNumber();
335 if (trackId != currentTrackId) {
336 auto superparentIndexIt = mSuperParentIndices.find(parentId);
337 if (superparentIndexIt != mSuperParentIndices.end()) {
338 superparentId = superparentIndexIt->second;
339 mSuperParentIndices[trackId] = superparentIndexIt->second;
342 mSuperParentIndices[trackId] = trackId;
343 superparentId = trackId;
345 currentTrackId = trackId;
348 double eloss = fMC->Edep();
349 if (eloss < DBL_EPSILON) {
353 TString volName = vol->GetName();
354 bool isCrystal = volName.Contains(
"crystalmodule");
355 bool isSampling = volName.Contains(
"samplingmodule");
357 if (!isCrystal && !isSampling) {
362 LOGP(
debug,
"Processing crystal {}", volName.Data());
364 eloss *= mSamplingFactorTransportModel;
365 LOGP(
debug,
"Processing scintillator {}", volName.Data());
367 int sectorId, moduleId;
368 fMC->CurrentVolID(moduleId);
369 fMC->CurrentVolOffID(1, sectorId);
371 LOGP(
debug,
"isCrystal={} sectorId={} moduleId={} cellID={} eloss={}", isCrystal, sectorId, moduleId, cellID, eloss);
373 int trackID = superparentId;
374 auto hit = std::find_if(mHits->begin(), mHits->end(), [cellID, trackID](
const Hit& hit) { return hit.GetTrackID() == trackID && hit.GetCellID() == cellID; });
375 if (hit == mHits->end()) {
376 float posX, posY, posZ, momX, momY, momZ, energy;
377 fMC->TrackPosition(posX, posY, posZ);
378 fMC->TrackMomentum(momX, momY, momZ, energy);
381 float time = fMC->TrackTime() * 1e9;
382 mHits->emplace_back(trackID, cellID,
pos, mom,
time, eloss);
383 stack->addHit(GetDetId());
Definition of the Stack class.
MC hit class to store energy loss per cell and per superparent.
Geometry parameters configurable via o2-sim –configKeyValues.
ECal geometry creation and hit processing.
ClassImp(o2::ecal::Detector)
void SetEnergyLoss(V val)
Detector()
Default Constructor.
void Mixture(Int_t imat, const char *name, Float_t *a, Float_t *z, Float_t dens, Int_t nlmat, Float_t *wmat)
void Medium(Int_t numed, const char *name, Int_t nmat, Int_t isvol, Int_t ifield, Float_t fieldm, Float_t tmaxfd, Float_t stemax, Float_t deemax, Float_t epsil, Float_t stmin, Float_t *ubuf=nullptr, Int_t nbuf=0)
static void initFieldTrackingParams(int &mode, float &maxfield)
void Material(Int_t imat, const char *name, Float_t a, Float_t z, Float_t dens, Float_t radl, Float_t absl, Float_t *buf=nullptr, Int_t nwbuf=0)
std::string addNameTo(const char *ext) const
static const ECalBaseParam & Instance()
bool ProcessHits(FairVolume *v=nullptr) override
void ConstructGeometry() override
void InitializeO2Detector() override
static GeometryTGeo * Instance()
static const char * getECalVolPattern()
static Geometry & instance()
int getCellID(int moduleId, int sectorId, bool isCrystal)
static ShmManager & Instance()
GLdouble GLdouble GLdouble z
void freeSimVector(std::vector< T > *ptr)
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Common utility functions.
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"