Project
Loading...
Searching...
No Matches
Detector.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
16
17#include <FairVolume.h>
18#include <TVirtualMC.h>
19#include <TVirtualMCStack.h>
20#include <TGeoManager.h>
21#include <TGeoVolume.h>
22#include <TGeoTube.h>
23#include <TGeoArb8.h>
24#include <TGeoTrd2.h>
25#include <Math/Point2D.h>
26#include <Math/Vector2D.h>
27#include <DetectorsBase/Stack.h>
28#include <ECalBase/Hit.h>
30#include <ECalBase/Geometry.h>
32
33using o2::ecal::Hit;
34
35namespace o2
36{
37namespace ecal
38{
40 : o2::base::DetImpl<Detector>("ECL", active),
41 mHits(o2::utils::createSimVector<o2::ecal::Hit>())
42{
43}
44
45//==============================================================================
47{
48 if (mHits) {
50 }
51}
52
53//==============================================================================
55{
56 createMaterials();
57 createGeometry();
58 defineSamplingFactor();
59}
60
61//==============================================================================
62void Detector::defineSamplingFactor()
63{
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.;
71 }
72}
73
74//==============================================================================
75void Detector::createMaterials()
76{
77 LOGP(info, "Creating materials for ECL");
78
79 // Air
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);
85
86 // Pb
87 Material(1, "Pb", 207.2, 82, 11.35, 0.56, 0., nullptr, 0);
88
89 // Polysterene scintillator (CH)
90 float aP[2] = {12.011, 1.00794};
91 float zP[2] = {6.0, 1.0};
92 float wP[2] = {1.0, 1.0};
93 float dP = 1.032;
94 Mixture(2, "Scintillator", aP, zP, dP, -2, wP);
95
96 // Al
97 Material(3, "Al", 26.98, 13., 2.7, 8.9, 999., nullptr, 0);
98
99 // PWO crystals
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};
103 float dX = 8.28;
104 Mixture(4, "PbWO4", aX, zX, dX, -3, wX);
105
106 int ifield = 2; // magnetic field flag
107 float fieldm = 10.; // maximum field value (in Kilogauss)
108 float deemax = 0.1; // maximum fractional energy loss in one step (0 < deemax <=1)
109 float tmaxfd = 10.; // maximum angle due to field permitted in one step (in degrees)
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);
116}
117
118//==============================================================================
120{
121 LOG(info) << "Initialize ECal O2Detector";
122 mGeometryTGeo = GeometryTGeo::Instance();
123}
124
125//==============================================================================
127{
128 if (!o2::utils::ShmManager::Instance().isOperational()) {
129 mHits->clear();
130 }
131 mSuperParentIndices.clear();
132 currentTrackId = -1;
133 superparentId = -1;
134}
135
136//==============================================================================
138{
139 // This will create a branch in the output tree called Hit, setting the last
140 // parameter to kFALSE means that this collection will not be written to the file,
141 // it will exist only during the simulation
142 LOGP(info, "Registering hits");
143 if (FairRootManager::Instance()) {
144 FairRootManager::Instance()->RegisterAny(addNameTo("Hit").data(), mHits, true);
145 }
146}
147
148//==============================================================================
149void Detector::createGeometry()
150{
151 LOGP(info, "Creating ECal geometry");
152
153 TGeoVolume* vALIC = gGeoManager->GetVolume("barrel");
154 if (!vALIC) {
155 LOGP(fatal, "Could not find barrel volume while constructing ECal geometry");
156 }
157 new TGeoVolumeAssembly(GeometryTGeo::getECalVolPattern());
158 TGeoVolume* vECal = gGeoManager->GetVolume(GeometryTGeo::getECalVolPattern());
159 vALIC->AddNode(vECal, 2, new TGeoTranslation(0, 30., 0));
160 vECal->SetTitle("ECalVol");
161
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");
167
168 // Get relevant parameters
169 auto& pars = ECalBaseParam::Instance();
170 auto& geo = Geometry::instance();
171
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);
180
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);
192
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);
212
213 double sectorR = rMin + sectorL / 2.;
214 for (int ism = 0; ism < pars.nSuperModules; ism++) {
215 // crystal
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));
223 }
224 // sampling
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));
234 }
235 }
236
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));
259 }
260
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));
282
283 // adding front aluminium plate into the volume
284 double fdx1 = dx1;
285 double fdx2 = dx1 + pars.frontPlateThickness * tanSamplingAlpha;
286 double fdy1 = dy1;
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);
295
296 // adding lead plates
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);
312 }
313 }
314
315 if (pars.enableFwdEndcap) {
316 // Build the ecal endcap
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));
322 }
323 // gGeoManager->CloseGeometry();
324 // gGeoManager->CheckOverlaps(0.0001);
325}
326
327//==============================================================================
328bool Detector::ProcessHits(FairVolume* vol)
329{
330 LOGP(debug, "Processing hits");
331 auto stack = (o2::data::Stack*)fMC->GetStack();
332 int trackId = stack->GetCurrentTrackNumber();
333 int parentId = stack->GetCurrentParentTrackNumber();
334
335 if (trackId != currentTrackId) {
336 auto superparentIndexIt = mSuperParentIndices.find(parentId);
337 if (superparentIndexIt != mSuperParentIndices.end()) {
338 superparentId = superparentIndexIt->second;
339 mSuperParentIndices[trackId] = superparentIndexIt->second;
340 } else {
341 // for new incoming tracks the superparent index is equal to the track ID (for recursion)
342 mSuperParentIndices[trackId] = trackId;
343 superparentId = trackId;
344 }
345 currentTrackId = trackId;
346 }
347
348 double eloss = fMC->Edep();
349 if (eloss < DBL_EPSILON) {
350 return false;
351 }
352
353 TString volName = vol->GetName();
354 bool isCrystal = volName.Contains("crystalmodule");
355 bool isSampling = volName.Contains("samplingmodule");
356
357 if (!isCrystal && !isSampling) {
358 return false;
359 }
360
361 if (isCrystal) {
362 LOGP(debug, "Processing crystal {}", volName.Data());
363 } else {
364 eloss *= mSamplingFactorTransportModel;
365 LOGP(debug, "Processing scintillator {}", volName.Data());
366 }
367 int sectorId, moduleId;
368 fMC->CurrentVolID(moduleId);
369 fMC->CurrentVolOffID(1, sectorId);
370 int cellID = Geometry::instance().getCellID(moduleId, sectorId, isCrystal);
371 LOGP(debug, "isCrystal={} sectorId={} moduleId={} cellID={} eloss={}", isCrystal, sectorId, moduleId, cellID, eloss);
372
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);
379 auto pos = math_utils::Point3D<float>(posX, posY, posZ);
380 auto mom = math_utils::Vector3D<float>(momX, momY, momZ);
381 float time = fMC->TrackTime() * 1e9; // time in ns
382 mHits->emplace_back(trackID, cellID, pos, mom, time, eloss);
383 stack->addHit(GetDetId());
384 } else {
385 hit->SetEnergyLoss(hit->GetEnergyLoss() + eloss);
386 }
387 return true;
388}
389
390} // namespace ecal
391} // namespace o2
392
Definition of the Stack class.
MC hit class to store energy loss per cell and per superparent.
Geometry parameters configurable via o2-sim –configKeyValues.
int16_t time
Definition RawEventData.h:4
int32_t i
uint16_t pos
Definition RawData.h:3
uint32_t stack
Definition RawData.h:1
Geometry helper class.
ECal geometry creation and hit processing.
ClassImp(o2::ecal::Detector)
std::ostringstream debug
V GetEnergyLoss() const
Definition BaseHits.h:103
void SetEnergyLoss(V val)
Definition BaseHits.h:104
Detector()
Default Constructor.
Definition Detector.cxx:36
void Mixture(Int_t imat, const char *name, Float_t *a, Float_t *z, Float_t dens, Int_t nlmat, Float_t *wmat)
Definition Detector.cxx:66
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)
Definition Detector.cxx:72
static void initFieldTrackingParams(int &mode, float &maxfield)
Definition Detector.cxx:143
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)
Definition Detector.cxx:59
std::string addNameTo(const char *ext) const
Definition Detector.h:150
void Register() override
Definition Detector.cxx:137
bool ProcessHits(FairVolume *v=nullptr) override
Definition Detector.cxx:328
void ConstructGeometry() override
Definition Detector.cxx:54
void InitializeO2Detector() override
Definition Detector.cxx:119
void Reset() override
Definition Detector.cxx:126
static GeometryTGeo * Instance()
static const char * getECalVolPattern()
static Geometry & instance()
Definition Geometry.h:31
int getCellID(int moduleId, int sectorId, bool isCrystal)
Definition Geometry.cxx:150
static ShmManager & Instance()
Definition ShmManager.h:61
GLint GLenum GLint x
Definition glcorearb.h:403
const GLfloat * m
Definition glcorearb.h:4066
GLboolean * data
Definition glcorearb.h:298
GLboolean r
Definition glcorearb.h:1233
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
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"
std::vector< int > row