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
12#include <TParticle.h>
13#include <TVirtualMC.h>
14#include <TGeoVolume.h>
15#include <TGeoManager.h>
16#include <TGeoBBox.h>
17#include <TGeoCompositeShape.h>
18
19#include <FairVolume.h>
20
21#include "DetectorsBase/Stack.h"
23#include "FOCALBase/Geometry.h"
24#include "FOCALBase/Hit.h"
25
26using namespace o2::focal;
27
28Detector::Detector(bool active, std::string geofilename)
29 : o2::base::DetImpl<Detector>("FOC", active),
30 mHits(o2::utils::createSimVector<Hit>()),
31 mHitIndexMapping(),
32 mGeometry(nullptr),
33 mMedSensHCal(-1),
34 mMedSensECalPad(-1),
35 mMedSensECalPix(-1),
36 mGeoCompositions(),
37 mSuperParentsIndices(),
38 mSuperParents(),
39 mCurrentSuperparent(nullptr),
40 mCurrentTrack(-1),
41 mCurrentPrimaryID(-1),
42 mCurrentParentID(-1),
43 mVolumeIDScintillator(-1)
44{
45 mGeometry = getGeometry(geofilename);
46 if (!mGeometry) {
47 LOG(fatal) << "Geometry is nullptr";
48 }
49}
50
52 : o2::base::DetImpl<Detector>(rhs)
53{
54 mGeometry = rhs.mGeometry;
55 mMedSensHCal = rhs.mMedSensHCal;
56 mMedSensECalPad = rhs.mMedSensECalPad;
57 mMedSensECalPix = rhs.mMedSensECalPix;
58 mSensitive = rhs.mSensitive;
59 // mSensitiveHCAL = rhs.mSensitiveHCAL;
60 // mSensitiveECALPad = rhs.mSensitiveECALPad;
61 // mSensitiveECALPix = rhs.mSensitiveECALPix;
62 mVolumeIDScintillator = rhs.mVolumeIDScintillator;
63}
64
69
71{
72 if (!mGeometry) {
73 mGeometry = Geometry::getInstance(name);
74 }
75 if (!mGeometry) {
76 LOG(error) << "Failure accessing geometry";
77 }
78 return mGeometry;
79}
80
82{
83
84 LOG(info) << "Intializing FOCAL detector";
85
86 // All FOCAL volumes must be declared as sensitive, otherwise
87 // the decay chains are broken by volumes not processed in ProceeHits
88 for (const auto& child : mSensitive) {
89 LOG(debug1) << "Adding sensitive volume " << child;
90 auto svolID = registerSensitiveVolumeAndGetVolID(child);
91 // HCAL
92 if (child == "ScintFiber" || child == "HScint") {
93 LOG(debug1) << "Adding ScintFiber/HScint volume as sensitive volume with ID " << svolID;
94 mVolumeIDScintillator = svolID;
95 }
96 // ECAL Pads
97 else if (child == "EMSC1" || child == "EMSC2") {
98 LOG(debug1) << "Adding EMC SILICON volume as sensitive volume with ID " << svolID;
99 mVolumeIDScintillator = svolID;
100 }
101 }
102
103 mMedSensHCal = getMediumID(ID_SC);
104 mMedSensECalPad = getMediumID(ID_SIPAD);
105 mMedSensECalPix = getMediumID(ID_SIPIX);
106}
107
108Bool_t Detector::ProcessHits(FairVolume* v)
109{
110 int track = fMC->GetStack()->GetCurrentTrackNumber(),
111 directparent = fMC->GetStack()->GetCurrentParentTrackNumber();
112 // Like other calorimeters FOCAL will create a huge amount of shower particles during tracking
113 // Instead, the hits should be assigned to the incoming particle in FOCAL.
114 // Implementation of the incoming particle search taken from implementation in EMCAL.
115 if (track != mCurrentTrack) {
116 LOG(debug4) << "Doing new track " << track << " current (" << mCurrentTrack << "), direct parent (" << directparent << ")";
117 // new current track - check parentage
118 auto hasSuperParent = mSuperParentsIndices.find(directparent);
119 if (hasSuperParent != mSuperParentsIndices.end()) {
120 // same superparent as direct parent
121 mCurrentParentID = hasSuperParent->second;
122 mSuperParentsIndices[track] = hasSuperParent->second;
123 auto superparent = mSuperParents.find(mCurrentParentID);
124 if (superparent != mSuperParents.end()) {
125 mCurrentSuperparent = &(superparent->second);
126 } else {
127 LOG(error) << "Attention: No superparent object found (parent " << mCurrentParentID << ")";
128 mCurrentSuperparent = nullptr;
129 }
130 LOG(debug4) << "Found superparent " << mCurrentParentID;
131 } else {
132 // start of new chain
133 // for new incoming tracks the super parent index is equal to the track ID (for recursion)
134 mSuperParentsIndices[track] = track;
135 mCurrentSuperparent = AddSuperparent(track, fMC->TrackPid(), fMC->Etot());
136 mCurrentParentID = track;
137 }
138 mCurrentTrack = track;
139 }
140
141 // Processing HCAL hits
142 bool flagHCAL = true;
143 if (fMC->CurrentMedium() == mMedSensHCal) {
144 flagHCAL = ProcessHitsHCAL(v);
145 }
146
147 // Processing ECAL Pad hits
148 bool flagECALPad = true;
149 if (TVirtualMC::GetMC()->CurrentMedium() == mMedSensECalPad) {
150 flagECALPad = ProcessHitsEPad(v);
151 }
152
153 // Processing ECAL Pixel hits
154 bool flagECALPix = true;
155 if (fMC->CurrentMedium() == mMedSensECalPix) {
156 flagECALPix = ProcessHitsEPix(v);
157 }
158
159 return (flagHCAL || flagECALPad || flagECALPix);
160 // return true;
161}
162
163Hit* Detector::AddHit(int trackID, int primary, double initialEnergy, int detID, o2::focal::Hit::Subsystem_t subsystem,
164 const math_utils::Point3D<float>& pos, double time, double eLoss)
165{
166 LOG(debug3) << "Adding hit for track " << trackID << " with position (" << pos.X() << ", "
167 << pos.Y() << ", " << pos.Z() << ") with energy " << initialEnergy << " loosing " << eLoss;
168 mHits->emplace_back(primary, trackID, detID, subsystem, initialEnergy, pos, time, eLoss);
169
170 auto [isin, col, row, layer, segment] = mGeometry->getVirtualInfo(pos.X(), pos.Y(), pos.Z());
171 mHitIndexMapping.insert(std::pair<Hit::HitID, unsigned int>({trackID, uint8_t(row), uint8_t(col), uint8_t(layer)}, static_cast<unsigned int>(mHits->size() - 1)));
172
173 return &(mHits->back());
174}
175
176Hit* Detector::FindHit(int parentID, int col, int row, int layer)
177{
178 Hit::HitID hitToFind{parentID, uint8_t(row), uint8_t(col), uint8_t(layer)};
179 auto found = mHitIndexMapping.find(hitToFind);
180
181 if (found == mHitIndexMapping.end()) {
182 return nullptr;
183 }
184
185 return &((*mHits)[found->second]);
186}
187
188Parent* Detector::AddSuperparent(int trackID, int pdg, double energy)
189{
190 LOG(debug3) << "Adding superparent for track " << trackID << " with PID " << pdg << " and energy " << energy;
191 auto entry = mSuperParents.insert({trackID, {pdg, energy, false}});
192 return &(entry.first->second);
193}
194
196
198{
199 FairRootManager::Instance()->RegisterAny(addNameTo("Hit").data(), mHits, kTRUE);
200}
201
203{
204 LOG(debug) << "Cleaning FOCAL hits ...";
205 if (!o2::utils::ShmManager::Instance().isOperational()) {
206 mHits->clear();
207 }
208 mHitIndexMapping.clear();
209
210 mSuperParentsIndices.clear();
211 mSuperParents.clear();
212 mCurrentTrack = -1;
213 mCurrentParentID = -1;
214}
215
217{
218
219 // --- Define the various materials for GEANT ---
220
222 float aSi = 28.09;
223 float zSi = 14.0;
224 float dSi = 2.33;
225 float x0Si = 9.36;
226 Material(1, "Si $", aSi, zSi, dSi, x0Si, 18.5);
227
229 float aW = 183.84;
230 float zW = 74.0;
231 float dW = 19.3;
232 float x0W = 0.35;
233 Material(0, "W $", aW, zW, dW, x0W, 17.1);
234
235 // Cu
236 Material(3, "Cu $", 63.54, 29., 8.96, 1.43, 15.);
237
238 // Al
239 Material(9, "Al$", 26.98, 13.0, 2.7, 8.9, 37.2);
240
242 Material(10, "Pb $", 207.19, 82., 11.35, .56, 18.5);
243
245 // --- The polysterene scintillator (CH) ---
246 float aP[2] = {12.011, 1.00794};
247 float zP[2] = {6.0, 1.0};
248 float wP[2] = {1.0, 1.0};
249 float dP = 1.032;
250 Mixture(11, "Polystyrene$", aP, zP, dP, -2, wP);
251
252 // G10
253
254 float aG10[4] = {1., 12.011, 15.9994, 28.086};
255 float zG10[4] = {1., 6., 8., 14.};
256 // PH float wG10[4]={0.148648649,0.104054054,0.483499056,0.241666667};
257 float wG10[4] = {0.15201, 0.10641, 0.49444, 0.24714};
258 Mixture(2, "G10 $", aG10, zG10, 1.7, 4, wG10);
259
261 float aAlloy[3] = {183.84, 58.6934, 63.54};
262 float zAlloy[3] = {74.0, 28, 29};
263 float wAlloy[3] = {0.94, 0.04, 0.02};
264 float dAlloy = wAlloy[0] * 19.3 + wAlloy[1] * 8.908 + wAlloy[2] * 8.96;
265 Mixture(5, "Alloy $", aAlloy, zAlloy, dAlloy, 3, wAlloy);
266
267 // Steel
268 float aSteel[4] = {55.847, 51.9961, 58.6934, 28.0855};
269 float zSteel[4] = {26., 24., 28., 14.};
270 float wSteel[4] = {.715, .18, .1, .005};
271 float dSteel = 7.88;
272 Mixture(4, "STAINLESS STEEL$", aSteel, zSteel, dSteel, 4, wSteel);
273
274 // Air
275 float aAir[4] = {12.0107, 14.0067, 15.9994, 39.948};
276 float zAir[4] = {6., 7., 8., 18.};
277 float wAir[4] = {0.000124, 0.755268, 0.231781, 0.012827};
278 float dAir1 = 1.20479E-10;
279 float dAir = 1.20479E-3;
280 Mixture(98, "Vacum$", aAir, zAir, dAir1, 4, wAir);
281 Mixture(99, "Air $", aAir, zAir, dAir, 4, wAir);
282
283 // Ceramic
284 // Ceramic 97.2% Al2O3 , 2.8% SiO2
285 // float wcer[2]={0.972,0.028}; // Not used
286 float aal2o3[2] = {26.981539, 15.9994};
287 float zal2o3[2] = {13., 8.};
288 float wal2o3[2] = {2., 3.};
289 float denscer = 3.6;
290 // SiO2
291 float aglass[2] = {28.0855, 15.9994};
292 float zglass[2] = {14., 8.};
293 float wglass[2] = {1., 2.};
294 float dglass = 2.65;
295 Mixture(6, "Al2O3 $", aal2o3, zal2o3, denscer, -2, wal2o3);
296 Mixture(7, "glass $", aglass, zglass, dglass, -2, wglass);
297
298 // Ceramic is a mixtur of glass and Al2O3 ?
299 // Not clear how to do this with AliMixture
300 // Not needed; so skip for now
301
302 /*
303float acer[2],zcer[2];
304char namate[21]="";
305float a,z,d,radl,absl,buf[1];
306Int_t nbuf;
307fMC->Gfmate((*fIdmate)[6], namate, a, z, d, radl, absl, buf, nbuf);
308acer[0]=a;
309zcer[0]=z;
310fMC->Gfmate((*fIdmate)[7], namate, a, z, d, radl, absl, buf, nbuf);
311acer[1]=a;
312zcer[1]=z;
313
314AliMixture( 8, "Ceramic $", acer, zcer, denscer, 2, wcer);
315*/
316
317 // Use Al2O3 instead:
318
319 Mixture(8, "Ceramic $", aal2o3, zal2o3, denscer, -2, wal2o3);
320
321 // Define tracking media
322 // format
323
324 float tmaxfdSi = 10.0; // 0.1; // .10000E+01; // Degree
325 float stemaxSi = 0.1; // .10000E+01; // cm
326 float deemaxSi = 0.1; // 0.30000E-02; // Fraction of particle's energy 0<deemax<=1
327 // float epsilSi = 1.e-3;//1e-3;//1.0E-4;// .10000E+01;
328 float epsilSi = 1.e-3; // 1.0E-4;// .10000E+01; // This drives the step size ? 1e-4 makes multiple steps even in pixels?
329 float stminSi = 0.001; // cm "Default value used"
330
331 float epsil = 0.001;
332 // MvL: need to look up itdmed dynamically?
333 // or move to TGeo: uses pointers for medium
334
335 int isxfld = 2;
336 float sxmgmx = 10.0;
338
340 Medium(ID_TUNGSTEN, "W conv.$", 0, 0,
341 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, nullptr, 0);
343 Medium(ID_SIPAD, "Si sens pad$", 1, 0,
344 isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi, nullptr, 0);
346 Medium(ID_SIPIX, "Si sens pix$", 1, 0,
347 isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi, nullptr, 0);
348
350 Medium(ID_G10, "G10 plate$", 2, 0,
351 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.01, nullptr, 0);
352
354 Medium(ID_COPPER, "Cu$", 3, 0,
355 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.0001, nullptr, 0);
356
358 Medium(ID_STEEL, "S steel$", 4, 0,
359 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.0001, nullptr, 0);
360
362 Medium(ID_ALLOY, "Alloy conv.$", 5, 0,
363 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, nullptr, 0);
364
366 Medium(ID_CERAMIC, "Ceramic$", 8, 0,
367 isxfld, sxmgmx, 10.0, 0.01, 0.1, 0.003, 0.003, nullptr, 0);
368
369 // HCAL materials // Need to double-check tracking pars for this
371 Medium(ID_PB, "Pb // The Scintillator must be first in order in vector for Rin to be set$", 10, 0,
372 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, nullptr, 0);
374 Medium(ID_SC, "Scint$", 11, 0,
375 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, nullptr, 0);
376
378 Medium(ID_SIINSENS, "Si insens$", 1, 0,
379 isxfld, sxmgmx, 10.0, 0.1, 0.1, epsil, 0.001, nullptr, 0);
380
381 // Al for the cold plates
382 Medium(ID_ALUMINIUM, "Aluminium$", 9, 0,
383 isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.1, 0.1, nullptr, 0);
384
386 Medium(ID_VAC, "Vacuum $", 98, 0,
387 isxfld, sxmgmx, 10.0, 1.0, 0.1, 0.1, 1.0, nullptr, 0);
388
390 Medium(ID_AIR, "Air gaps$", 99, 0,
391 isxfld, sxmgmx, 10.0, 1.0, 0.1, epsil, 0.001, nullptr, 0);
392}
393
394//____________________________________________________________________________
396{
397 // Create entries for alignable volumes associating the symbolic volume
398 // name with the corresponding volume path. Needs to be syncronized with
399 // eventual changes in the geometry
400 // Alignable volumes are:
401
404}
405
406//____________________________________________________________________________
408{
409 const std::string vpsector = "/cave_1/barrel_1/FOCAL_1/HCAL_1";
410 const std::string snsector = "FOCAL/HCAL";
411
412 if (!gGeoManager->SetAlignableEntry(snsector.c_str(), vpsector.c_str())) {
413 LOG(fatal) << "Alignable entry " << snsector << " not created. Volume path " << vpsector << "not valid.";
414 }
415}
416
417//____________________________________________________________________________
419{
420 const std::string vpsector = "/cave_1/barrel_1/FOCAL_1/ECAL_1";
421 const std::string snsector = "FOCAL/ECAL";
422
423 if (!gGeoManager->SetAlignableEntry(snsector.c_str(), vpsector.c_str())) {
424 LOG(fatal) << "Alignable entry " << snsector << " not created. Volume path " << vpsector << "not valid.";
425 }
426}
427
429{
430
441
442 LOG(debug) << "Creating FOCAL geometry\n";
443
445
447 mGeoCompositions = mGeometry->getFOCALMicroModule(-1);
448 if (!mGeoCompositions.size()) {
449 LOG(error) << "FOCAL compositions not found!!";
450 return;
451 }
452
453 float pars[4];
454 pars[0] = (mGeometry->getFOCALSizeX() + 2 * mGeometry->getMiddleTowerOffset() + mGeometry->getDetectorOpeningRight() + mGeometry->getDetectorOpeningLeft()) / 2;
455 pars[1] = mGeometry->getFOCALSizeY() / 2;
456 pars[2] = mGeometry->getFOCALSizeZ() / 2;
457 // Add space to place 2 SiPad layers in front of ECAL
458 // The global position of ECAL and HCAL remains the same, but the FOCAL box needs to be slightly larger to accomodate
459 // the 2 SiPad layers which will sit at z=698 and 699cm (2 and 1 cm in front of ECAL)
460 if (mGeometry->getInsertFrontPadLayers()) {
461 pars[2] += 1.0;
462 }
463 if (mGeometry->getInsertHCalReadoutMaterial()) {
464 pars[2] += (1.0 + 0.5); // place Aluminium 1cm thick box (0.5 means half) at 2cm behind HCal to simulate SiPM readout material
465 pars[1] += (10.0 + 1.0); // place Aluminium 1cm thick box at 10 cm below FOCAL
466 }
467 pars[3] = 0;
468
469 LOG(info) << "Creating FOCAL with dimensions X: " << (mGeometry->getFOCALSizeX() + 2 * mGeometry->getMiddleTowerOffset()) << ", Y: "
470 << mGeometry->getFOCALSizeY() << ", Z: " << mGeometry->getFOCALSizeZ() + (mGeometry->getInsertFrontPadLayers() ? 2.0 : 0.0) + (mGeometry->getInsertHCalReadoutMaterial() ? 3.0 : 0.0);
471
472 TVirtualMC::GetMC()->Gsvolu("FOCAL", "BOX", getMediumID(ID_AIR), pars, 4);
473 mSensitive.push_back("FOCAL");
474 // mSensitiveHCAL.push_back("FOCAL");
475 // mSensitiveECALPad.push_back("FOCAL");
476 // mSensitiveECALPix.push_back("FOCAL");
477
478 // ECAL part
479 LOG(debug2) << "ECAL geometry : " << GetTitle();
481
482 // HCAL part
483 switch (mGeometry->getHCALDesign()) {
486 break;
487
490 break;
491
494 break;
495
496 default:
497 break;
498 }
499
500 // const float z0 = 1312.5; // center of barrel mother volume
501 TVirtualMC::GetMC()->Gspos("FOCAL", 1, "barrel", 0, 30., mGeometry->getFOCALZ0() - (mGeometry->getInsertFrontPadLayers() ? 2.0 : 0.0) + (mGeometry->getInsertHCalReadoutMaterial() ? 1.5 : 0.0), 0, "ONLY");
502}
503
505{
506 TGeoVolumeAssembly* volHCAL = new TGeoVolumeAssembly("HCAL");
507
508 TGeoVolumeAssembly* HcalTube = gGeoManager->MakeVolumeAssembly("ScintCuTubes");
509
510 TGeoVolume* volCuTube;
511 TGeoVolume* volSciFi;
512
513 float RScint = 0.;
514 float Rin = 0.;
515 float Rout = 0.;
516 float Length = 0.;
517
518 for (auto& icomp : mGeoCompositions) {
519 Length = icomp->sizeZ() / 2;
520
521 if (icomp->material() == "Pb") {
522 Rout = icomp->sizeX() / 2;
523 TGeoMedium* medium = gGeoManager->GetMedium(getMediumID(ID_PB));
524 volCuTube = gGeoManager->MakeTube("Tube", medium, Rin, Rout, Length); // The Scintillator must be first in order in vector for Rin to be set
525 volCuTube->SetLineWidth(2);
526 volCuTube->SetLineColor(kRed);
527 mSensitive.push_back(volCuTube->GetName());
528 // mSensitiveHCAL.push_back(volCuTube->GetName());
529 HcalTube->AddNode(volCuTube, 1, nullptr);
530 }
531 if (icomp->material() == "Scint") {
532 RScint = icomp->sizeX() / 2;
533 Rin = RScint + 0.005;
534 TGeoMedium* medium = gGeoManager->GetMedium(getMediumID(ID_SC));
535 volSciFi = gGeoManager->MakeTube("ScintFiber", medium, 0., RScint, Length);
536 volSciFi->SetLineWidth(2);
537 volSciFi->SetLineColor(kBlue);
538 // mSensitiveHCAL.push_back(volSciFi->GetName());
539 mSensitive.push_back(volSciFi->GetName());
540 HcalTube->AddNode(volSciFi, 1, nullptr);
541 }
542 if (icomp->material() == "CuHCAL") {
543 Rout = icomp->sizeX() / 2;
544 TGeoMedium* medium = gGeoManager->GetMedium(getMediumID(ID_COPPER));
545 volCuTube = gGeoManager->MakeTube("Tube", medium, Rin, Rout, Length); // The Scintillator must be first in order in vector for Rin to be set
546 volCuTube->SetLineWidth(2);
547 volCuTube->SetLineColor(kRed);
548 // mSensitiveHCAL.push_back(volCuTube->GetName());
549 mSensitive.push_back(volCuTube->GetName());
550 HcalTube->AddNode(volCuTube, 1, nullptr);
551 }
552 }
553
554 bool splitDet = mGeometry->getDetectorOpeningRight() > 0.0 || mGeometry->getDetectorOpeningLeft() > 0.0;
555
556 double TowerSize = mGeometry->getHCALTowerSize();
557 double CuBoxThickness = 0.3; // Thickness of the Cu box carrying capillary tubes
558
559 TGeoBBox* ODBox = new TGeoBBox("TowerOD", TowerSize / 2, TowerSize / 2, Length);
560 TGeoBBox* IDBox = new TGeoBBox("TowerID", (TowerSize - CuBoxThickness) / 2, (TowerSize - CuBoxThickness) / 2, Length + 0.01);
561 TGeoCompositeShape* TowerHCAL = new TGeoCompositeShape("TowerHCAL", "TowerOD - TowerID");
562 TGeoVolume* volTower = new TGeoVolume("volTower", TowerHCAL, gGeoManager->GetMedium(getMediumID(ID_COPPER)));
563 volTower->SetLineWidth(2);
564 volTower->SetLineColor(42);
565 // mSensitiveHCAL.push_back(volTower->GetName());
566 mSensitive.push_back(volTower->GetName());
567
568 TGeoVolumeAssembly* volTowerHCAL = new TGeoVolumeAssembly("volTowerHCAL");
569 volTowerHCAL->AddNode(volTower, 1, nullptr);
570
571 int Rows = 0;
572 float RowPos = 0.;
573 int Columns = 0;
574 int NumTubes = 1;
575
576 // Packing circles in Hexagonal shape
577 while (RowPos + CuBoxThickness / 2 + Rout + 2 * Rout < TowerSize) {
578 Columns = 0;
579 float ColumnPos = (Rows % 2 == 0) ? 0. : Rout;
580 while (ColumnPos + CuBoxThickness / 2 + Rout + 2 * Rout < TowerSize) {
581
582 TGeoTranslation* trans = new TGeoTranslation(ColumnPos - TowerSize / 2 + CuBoxThickness / 2 + Rout, RowPos - TowerSize / 2 + CuBoxThickness / 2 + Rout, 0.);
583
584 trans->SetName(Form("trans_Num_%d", NumTubes));
585 trans->RegisterYourself();
586
587 volTowerHCAL->AddNode(HcalTube, NumTubes, trans);
588 // volTowerHCAL->AddNode(volCuTube, NumTubes, trans);
589 // volTowerHCAL->AddNode(volSciFi, NumTubes, trans);
590
591 Columns++;
592 ColumnPos = Columns * 2 * Rout + ((Rows % 2 == 0) ? 0. : Rout);
593 NumTubes++;
594 }
595
596 Rows++;
597 RowPos = Rows * 2 * Rout * TMath::Sin(TMath::Pi() / 3);
598 }
599
600 // Define the distance from the beam pipe in which towers will ommitted
601 Double_t BeamPipeRadius = 3.0; // in cm To be changed later
602 Double_t TowerHalfDiag = TMath::Sqrt2() * 0.5 * TowerSize; // tower half diagonal
603 Double_t MinRadius = BeamPipeRadius + TowerSize / 2;
604
605 float SizeXHCAL = mGeometry->getHCALTowersInX() * TowerSize;
606 float SizeYHCAL = mGeometry->getHCALTowersInY() * TowerSize;
607
608 int nTowersX = mGeometry->getHCALTowersInX();
609 int nTowersY = mGeometry->getHCALTowersInY();
610
611 Rows = 0;
612 Columns = 0;
613 RowPos = 0.;
614 Int_t NumTowers = 1;
615
616 if (splitDet) {
617 SizeXHCAL = SizeXHCAL / 2;
618
619 TGeoVolumeAssembly* volHalfHCAL = new TGeoVolumeAssembly("HalfHCAL");
620
621 for (Rows = 0; Rows < nTowersY; Rows++) {
622
623 float ColumnPos = 0.;
624 RowPos = Rows * TowerSize;
625 for (Columns = 0; Columns < nTowersX / 2; Columns++) {
626 ColumnPos = Columns * TowerSize;
627 TGeoTranslation* trans = new TGeoTranslation(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, RowPos - SizeYHCAL / 2 + TowerSize / 2, 0.);
628
629 // Shit the beampipe towers by TowerSize/2
630 if (Rows == nTowersY / 2) {
631 trans->SetDx(trans->GetTranslation()[0] + TowerSize / 2);
632 }
633
634 // Adding the Tower to the HCAL
635 volHalfHCAL->AddNode(volTowerHCAL, NumTowers, trans);
636
637 NumTowers++;
638 }
639 volHCAL->AddNode(volHalfHCAL, 1, new TGeoTranslation(SizeXHCAL / 2 + mGeometry->getDetectorOpeningRight(), 0, 0));
640 TGeoRotation* rotFlipZ = new TGeoRotation();
641 rotFlipZ->RotateY(180); // Flip around Y to reverse Z
642 TGeoCombiTrans* combHalf = new TGeoCombiTrans(-SizeXHCAL / 2 - mGeometry->getDetectorOpeningLeft(), 0., 0., rotFlipZ);
643 volHCAL->AddNode(volHalfHCAL, 2, combHalf);
644 }
645 } else {
646 for (Rows = 0; Rows < nTowersY; Rows++) {
647
648 float ColumnPos = 0.;
649 RowPos = Rows * TowerSize;
650 for (Columns = 0; Columns < nTowersX; Columns++) {
651 ColumnPos = Columns * TowerSize;
652 TGeoTranslation* trans = new TGeoTranslation(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, RowPos - SizeYHCAL / 2 + TowerSize / 2, 0.);
653
654 // Remove the Towers that overlaps with the beam pipe
655 Double_t RadialDistance = TMath::Power(trans->GetTranslation()[0], 2) + TMath::Power(trans->GetTranslation()[1], 2);
656
657 if (RadialDistance < MinRadius * MinRadius || TMath::Abs(trans->GetTranslation()[0]) > SizeXHCAL / 2) {
658 continue;
659 }
660
661 // Adding the Tower to the HCAL
662 volHCAL->AddNode(volTowerHCAL, NumTowers, trans);
663
664 NumTowers++;
665 }
666 }
667 }
668
669 LOG(info) << "Number of Towers is: " << (NumTowers - 1);
670 LOG(info) << "Number of tubes is: " << (NumTubes - 1) * (NumTowers - 1);
671
672 // Create Aluminium plate at the back of HCal to simulate the electronics readout material
673 // Hardcoded thickness of 1cm and placement at 2 cm behind HCAL
674 TGeoBBox* alHcalBox = new TGeoBBox("AlHCalBox", SizeXHCAL / 2.0, SizeYHCAL / 2.0, 0.5 / 2.0);
675 TGeoVolume* volumeAlHcalBox = new TGeoVolume("volAlHcalBox", alHcalBox, gGeoManager->GetMedium(getMediumID(ID_ALUMINIUM)));
676 volumeAlHcalBox->SetLineColor(kOrange);
677 if (mGeometry->getInsertHCalReadoutMaterial()) {
678 TVirtualMC::GetMC()->Gspos("volAlHcalBox", 9999, "FOCAL", 0.0, 0.0, +1.0 * mGeometry->getFOCALSizeZ() / 2.0 + 1.0, 0, "ONLY");
679 // mSensitiveHCAL.push_back("volAlHcalBox");
680 mSensitive.push_back("volAlHcalBox");
681 }
682 TGeoBBox* alUnderBox = new TGeoBBox("AlUnderBox", SizeXHCAL / 2.0, 0.5, mGeometry->getFOCALSizeZ() / 2.0 + 1.5);
683 TGeoVolume* volumeAlUnderBox = new TGeoVolume("volAlUnderBox", alUnderBox, gGeoManager->GetMedium(getMediumID(ID_ALUMINIUM)));
684 volumeAlUnderBox->SetLineColor(kOrange);
685 if (mGeometry->getInsertHCalReadoutMaterial()) {
686 TVirtualMC::GetMC()->Gspos("volAlUnderBox", 9999, "FOCAL", 0.0, -1.0 * mGeometry->getFOCALSizeY() / 2 - 10.5, 0.0, 0, "ONLY");
687 // mSensitiveHCAL.push_back("volAlUnderBox");
688 mSensitive.push_back("volAlUnderBox");
689 }
690
691 volHCAL->SetVisibility();
692 volHCAL->SetVisDaughters();
693 TVirtualMC::GetMC()->Gspos("HCAL", 1, "FOCAL", 0, 0, mGeometry->getHCALCenterZ() - mGeometry->getFOCALSizeZ() / 2 + 0.01 + (mGeometry->getInsertFrontPadLayers() ? 2.0 : 0.0) - (mGeometry->getInsertHCalReadoutMaterial() ? 1.5 : 0.0), 0, "ONLY");
694}
695
696//_____________________________________________________________________________
697TGeoVolumeAssembly* Detector::CreatePitchAssembly(double Lx,
698 double Ly1,
699 double Ly2,
700 double Lz,
701 double hole_diameter,
702 double hole_spacing,
703 int nholes,
704 double fiber_radius,
705 std::string suffix)
706{
707
708 // Z-alignment doesn't change
709 double zpos = 0;
710
711 TGeoMedium* copper = gGeoManager->GetMedium(getMediumID(ID_COPPER));
712 TGeoMedium* scint = gGeoManager->GetMedium(getMediumID(ID_SC));
713
714 TGeoVolumeAssembly* pitchAssembly = new TGeoVolumeAssembly("pitchAssembly");
715
716 // Hardcoded values for hole placement, to be set from outside
717 float holeStart = 0.15; // cm
718 float holeEnd = 0.35; // cm
719
720 TGeoVolumeAssembly* volLowerSheetwHoles = new TGeoVolumeAssembly(Form("volLowerSheetwHoles_%s", suffix.c_str()));
721 TGeoVolume* cuSheet = gGeoManager->MakeBox("cuSheet", copper, Lx / 2, (Ly1 - fiber_radius * 2) / 2, Lz / 2);
722 cuSheet->SetLineColor(kOrange + 2);
723 mSensitive.push_back(cuSheet->GetName());
724 TGeoVolume* boxbegin = gGeoManager->MakeBox("BoxBegin", copper, holeStart / 2, fiber_radius, Lz / 2);
725 boxbegin->SetLineColor(kOrange + 2);
726 mSensitive.push_back(boxbegin->GetName());
727 TGeoVolume* boxMiddle = gGeoManager->MakeBox("BoxMiddle", copper, (hole_spacing - hole_diameter) / 2, fiber_radius, Lz / 2);
728 boxMiddle->SetLineColor(kOrange + 2);
729 mSensitive.push_back(boxMiddle->GetName());
730 TGeoVolume* boxEnd = gGeoManager->MakeBox("BoxEnd", copper, holeEnd / 2, fiber_radius, Lz / 2);
731 boxEnd->SetLineColor(kOrange + 2);
732 mSensitive.push_back(boxEnd->GetName());
733
734 double yPlacement = Ly1 / 2 - fiber_radius;
735
736 // -----------------
737 // Layer 1: Lower sheet with holes (y = 0)
738 // -----------------
739
740 volLowerSheetwHoles->AddNode(cuSheet, 0, new TGeoTranslation(0, -Ly1 / 2 + (Ly1 - fiber_radius * 2) / 2, zpos));
741
742 // Add holes starting at x = 1.5 mm
743 float start_x = -Lx / 2 + holeStart;
744
745 for (int ihole = 0; ihole < nholes; ++ihole) {
746 float holePlacement = start_x + ihole * hole_spacing + hole_diameter / 2;
747 if (ihole == 0) {
748 volLowerSheetwHoles->AddNode(boxbegin, ihole, new TGeoTranslation(holePlacement - holeStart / 2 - hole_diameter / 2, yPlacement, zpos));
749 volLowerSheetwHoles->AddNode(boxMiddle, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + (hole_spacing - hole_diameter) / 2, yPlacement, zpos));
750 } else if (ihole == nholes - 1) {
751 if ((holePlacement + hole_diameter / 2 + holeStart) < Lx / 2 - 0.005) {
752 volLowerSheetwHoles->AddNode(boxEnd, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + holeEnd / 2, yPlacement, zpos));
753 } else {
754 volLowerSheetwHoles->AddNode(boxbegin, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + holeStart / 2, yPlacement, zpos));
755 }
756 } else {
757 volLowerSheetwHoles->AddNode(boxMiddle, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + (hole_spacing - hole_diameter) / 2, yPlacement, zpos));
758 }
759 }
760
761 pitchAssembly->AddNode(volLowerSheetwHoles, 0, new TGeoTranslation(0, Ly1 / 2, zpos)); // Add Ly1 / 2 so the lower edge of the sheets start y=0
762
763 // -----------------
764 // Layer 2: Full copper sheet
765 // -----------------
766 TGeoVolume* fullSheet1 = gGeoManager->MakeBox("FullSheet1", copper, Lx / 2, Ly2 / 2, Lz / 2);
767 fullSheet1->SetLineColor(kOrange + 2);
768 mSensitive.push_back(fullSheet1->GetName());
769 pitchAssembly->AddNode(fullSheet1, 0, new TGeoTranslation(0, Ly1 / 2 + Ly2 / 2 + Ly1 / 2, zpos)); // Add Ly1 / 2 so the lower edge of the sheets start y=0
770
771 // -----------------
772 // Layer 3: Upper sheet with holes (shifted)
773 // -----------------
774
775 TGeoVolumeAssembly* volUpperSheetwHoles = new TGeoVolumeAssembly(Form("volUpperSheetwHoles_%s", suffix.c_str()));
776
777 volUpperSheetwHoles->AddNode(cuSheet, 0, new TGeoTranslation(0, -Ly1 / 2 + (Ly1 - fiber_radius * 2) / 2, zpos));
778
779 // Add holes starting at x = 3.5 mm
780 float start_x2 = -Lx / 2 + holeEnd;
781
782 for (int ihole = 0; ihole < nholes; ++ihole) {
783 float holePlacement = start_x2 + ihole * hole_spacing + hole_diameter / 2;
784 if (ihole == 0) {
785 volUpperSheetwHoles->AddNode(boxEnd, ihole, new TGeoTranslation(holePlacement - hole_diameter / 2 - holeEnd / 2, yPlacement, zpos));
786 volUpperSheetwHoles->AddNode(boxMiddle, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + (hole_spacing - hole_diameter) / 2, yPlacement, zpos));
787 } else if (ihole == nholes - 1) {
788 volUpperSheetwHoles->AddNode(boxbegin, ihole, new TGeoTranslation(holePlacement + holeStart / 2 + hole_diameter / 2, yPlacement, zpos));
789 } else {
790 if ((holePlacement + hole_spacing + hole_diameter / 2) < Lx / 2 - 0.005) {
791 volUpperSheetwHoles->AddNode(boxMiddle, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + (hole_spacing - hole_diameter) / 2, yPlacement, zpos));
792 } else {
793 volUpperSheetwHoles->AddNode(boxEnd, ihole, new TGeoTranslation(holePlacement + hole_diameter / 2 + holeEnd / 2, yPlacement, zpos));
794 break;
795 }
796 }
797 }
798
799 pitchAssembly->AddNode(volUpperSheetwHoles, 0, new TGeoTranslation(0, Ly1 / 2 + Ly2 + Ly1 / 2 + Ly1 / 2, zpos)); // Add Ly1 / 2 so the lower edge of the sheets start y=0
800
801 // -----------------
802 // Layer 4: Full copper sheet
803 // -----------------
804 pitchAssembly->AddNode(fullSheet1, 1, new TGeoTranslation(0, Ly1 / 2 + Ly2 + Ly1 + Ly2 / 2 + Ly1 / 2, zpos)); // Add Ly1 / 2 so the lower edge of the sheets start y=0
805
806 // -----------------
807 // Scintillator Fibers
808 // -----------------
809 // Lower set of fibers
810 TGeoVolume* fiber = gGeoManager->MakeTube("Fiber", scint, 0, fiber_radius, Lz / 2);
811 fiber->SetLineColor(kBlue);
812 mSensitive.push_back(fiber->GetName());
813 for (int i = 0; i < nholes; ++i) {
814 float x_fiber = start_x + i * hole_spacing + hole_diameter / 2;
815 pitchAssembly->AddNode(fiber, i, new TGeoTranslation(x_fiber, Ly1 / 2 - fiber_radius + Ly1 / 2, zpos));
816 }
817
818 // Upper set of fibers
819 for (int i = 0; i < nholes; ++i) {
820 float x_fiber = start_x2 + i * hole_spacing + hole_diameter / 2;
821 if (x_fiber > Lx / 2 - 0.05) {
822 break;
823 }
824 pitchAssembly->AddNode(fiber, i + nholes, new TGeoTranslation(x_fiber, Ly1 / 2 + Ly2 + Ly1 / 2 + Ly1 / 2 - fiber_radius + Ly1 / 2, zpos));
825 }
826
827 return pitchAssembly;
828}
829
830//_____________________________________________________________________________
832{
833 TGeoVolumeAssembly* volHCAL = new TGeoVolumeAssembly("HCAL");
834
835 // Dimensions
836 double Lx = 49.81; // cm
837 double Ly1 = 0.20; // cm (sheets with holes)
838 double Ly2 = 0.15; // cm (full sheets)
839 double Lz = 110.0; // cm
840
841 double fiber_radius = 0.05;
842
843 // HCal materials
844 int icomp = 0;
845 for (auto& comp : mGeoCompositions) {
846 Lz = comp->sizeZ();
847
848 if (comp->material() == "Scint") {
849 fiber_radius = comp->sizeX() / 2;
850 }
851 if (comp->material() == "CuHCAL" && icomp == 0) {
852 Lx = comp->sizeX();
853 Ly1 = comp->sizeY();
854 }
855 if (comp->material() == "CuHCAL" && icomp == 2) {
856 Ly2 = comp->sizeY();
857 }
858 icomp++;
859 }
860
861 double hole_diameter = fiber_radius * 2 + 0.01; // hole radius
862 double hole_spacing = mGeometry->getHCALPitchSize();
863 int nholes = (int)(Lx / hole_spacing); // Number of holes in one HCAL sheet
864
865 double beamPipeHole = mGeometry->getHCALBeamPipeHoleSize(); // cm The size of the beam pipe opening
866 int nBeamPipeHoles = (int)((Lx - beamPipeHole / 2) / hole_spacing); // Number of beam pipe holes
867
868 // Compute module height (two sheets with holes + two full sheets)
869 float pitch_height = Ly1 + Ly2 + Ly1 + Ly2;
870
871 int totalNumberOfPitches = mGeometry->getHCALTowersInY() * 2; // Number of pitches in the whole HCAL
872 int numberOfPitchesBeamPipe = (int)((beamPipeHole + 0.001) / pitch_height); // Number of pitches in the beam pipe region
873 int numberofPitchesOnYaxis = (totalNumberOfPitches - numberOfPitchesBeamPipe); // Number of pitches in the HCAL ouside the beam pipe region
874
875 TGeoVolumeAssembly* pitchAssembly = CreatePitchAssembly(Lx, Ly1, Ly2, Lz, hole_diameter, hole_spacing, nholes, fiber_radius, "Main");
876 pitchAssembly->SetVisibility(true);
877 TGeoVolumeAssembly* beamPipeAssembly = CreatePitchAssembly(Lx - beamPipeHole / 2, Ly1, Ly2, Lz, hole_diameter, hole_spacing, nBeamPipeHoles, fiber_radius, "BeamPipe");
878 beamPipeAssembly->SetVisibility(true);
879
880 TGeoVolumeAssembly* HalfHCAL = new TGeoVolumeAssembly("HalfHCAL");
881
882 for (int iPitch = 0; iPitch < numberofPitchesOnYaxis; iPitch++) {
883 float placement = iPitch * pitch_height - pitch_height * (totalNumberOfPitches) / 2.0;
884 if (placement < -beamPipeHole / 2.0) {
885 HalfHCAL->AddNode(pitchAssembly, iPitch, new TGeoTranslation(0, placement, 0.));
886 } else {
887 placement += beamPipeHole;
888 HalfHCAL->AddNode(pitchAssembly, iPitch, new TGeoTranslation(0, placement, 0.));
889 }
890 }
891
892 for (int iPitch = 0; iPitch < numberOfPitchesBeamPipe; iPitch++) {
893 float placement = iPitch * pitch_height - beamPipeHole / 2.0;
894 HalfHCAL->AddNode(beamPipeAssembly, iPitch, new TGeoTranslation(-beamPipeHole / 4, placement, 0.));
895 }
896
897 HalfHCAL->SetVisibility(true);
898 HalfHCAL->SetVisDaughters(true);
899
900 volHCAL->AddNode(HalfHCAL, 0, new TGeoTranslation(-Lx / 2, 0, 0.));
901 TGeoRotation* rotFlipZ = new TGeoRotation();
902 rotFlipZ->RotateY(180); // Flip around Y to reverse Z
903 TGeoCombiTrans* combHalf = new TGeoCombiTrans(Lx / 2, 0., 0., rotFlipZ);
904 volHCAL->AddNode(HalfHCAL, 1, combHalf);
905
906 gMC->Gspos("HCAL", 1, "FOCAL", 0, 0, mGeometry->getHCALCenterZ() - mGeometry->getFOCALSizeZ() / 2 + 0.01 + (mGeometry->getInsertFrontPadLayers() ? 2.0 : 0.0) - (mGeometry->getInsertHCalReadoutMaterial() ? 1.5 : 0.0), 0, "ONLY");
907}
908
909//_____________________________________________________________________________
911{
912 TGeoVolumeAssembly* volHCAL = new TGeoVolumeAssembly("HCAL");
913
915 Float_t pars[4]; // this is HMSC Assembly
916 pars[0] = mGeometry->getHCALTowerSize() / 2;
917 pars[1] = mGeometry->getHCALTowerSize() / 2;
918 pars[2] = mGeometry->getECALSizeZ() + mGeometry->getHCALSizeZ() / 2; // ECAL sizeZ is already added to the HCAL materials CenterZ, so it is also treated as offset
919 pars[3] = 0;
920
921 float offset = pars[2];
922
923 TGeoVolumeAssembly* volTower = new TGeoVolumeAssembly("Tower");
924
925 int iCu(0), iScint(0);
926 for (auto& icomp : mGeoCompositions) {
927
928 pars[0] = icomp->sizeX() / 2;
929 pars[1] = icomp->sizeY() / 2;
930 pars[2] = icomp->sizeZ() / 2;
931 pars[3] = 0;
932
933 // HCal materials
934
935 if (icomp->material() == "Pb") {
936 iCu++;
937 const TGeoMedium* medium = gGeoManager->GetMedium(getMediumID(ID_PB));
938 const TGeoBBox* HPadBox = new TGeoBBox("HPadBox", pars[0], pars[1], pars[2]);
939 TGeoVolume* HPad = new TGeoVolume("HPad", HPadBox, medium);
940 HPad->SetLineColor(kGray);
941 // mSensitiveHCAL.push_back(HPad->GetName());
942 mSensitive.push_back(HPad->GetName());
943 TGeoTranslation* trans = new TGeoTranslation(icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset);
944 volTower->AddNode(HPad, iCu, trans);
945 }
946 if (icomp->material() == "Scint") {
947 iScint++;
948 const TGeoMedium* medium = gGeoManager->GetMedium(getMediumID(ID_SC));
949 const TGeoBBox* HScintBox = new TGeoBBox("HScintBox", pars[0], pars[1], pars[2]);
950 TGeoVolume* HScint = new TGeoVolume("HScint", HScintBox, medium);
951 HScint->SetLineColor(kBlue);
952 // mSensitiveHCAL.push_back(HScint->GetName());
953 mSensitive.push_back(HScint->GetName());
954 TGeoTranslation* trans = new TGeoTranslation(icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset);
955 volTower->AddNode(HScint, iScint, trans);
956 }
957 if (icomp->material() == "CuHCAL") {
958 iCu++;
959 const TGeoMedium* medium = gGeoManager->GetMedium(getMediumID(ID_COPPER));
960 const TGeoBBox* HPadBox = new TGeoBBox("HPadBox", pars[0], pars[1], pars[2]);
961 TGeoVolume* HPad = new TGeoVolume("HPad", HPadBox, medium);
962 HPad->SetLineColor(kRed);
963 // mSensitiveHCAL.push_back(HPad->GetName());
964 mSensitive.push_back(HPad->GetName());
965 TGeoTranslation* trans = new TGeoTranslation(icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset);
966 volTower->AddNode(HPad, iCu, trans);
967 }
968 }
969 double TowerSize = mGeometry->getHCALTowerSize();
970
971 // Define the distance from the beam pipe in which towers will ommitted
972 double BeamPipeRadius = 3.6; // in cm
973 double TowerHalfDiag = TMath::Sqrt2() * 0.5 * TowerSize; // tower half diagonal
974 double MinRadius = BeamPipeRadius + TowerHalfDiag;
975
976 float SizeXHCAL = mGeometry->getHCALTowersInX() * TowerSize;
977 float SizeYHCAL = mGeometry->getHCALTowersInY() * TowerSize;
978
979 int nTowersX = mGeometry->getHCALTowersInX();
980 int nTowersY = mGeometry->getHCALTowersInY();
981
982 int Rows = 0;
983 int Columns = 0;
984 double RowPos = 0.;
985 int NumTowers = 1;
986
987 // Arranging towers
988 for (Rows = 0; Rows < nTowersY; Rows++) {
989 Columns = 0;
990 float ColumnPos = 0.;
991 RowPos = Rows * TowerSize;
992 for (Columns = 0; Columns < nTowersX; Columns++) {
993 ColumnPos = Columns * TowerSize;
994
995 TGeoTranslation* trans = new TGeoTranslation(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, RowPos - SizeYHCAL / 2 + TowerSize / 2, 0.);
996
997 // Remove the Towers that overlaps with the beam pipe
998 double RadialDistance = TMath::Power(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, 2) + TMath::Power(RowPos - SizeYHCAL / 2 + TowerSize / 2, 2);
999
1000 if (RadialDistance < MinRadius * MinRadius) {
1001 continue;
1002 }
1003
1004 // Adding the Tower to the HCAL
1005 volHCAL->AddNode(volTower, NumTowers, trans);
1006
1007 NumTowers++;
1008 }
1009 }
1010 LOG(info) << "Number of Towers is: " << (NumTowers - 1);
1011
1012 // Create an Aluminium plate at the back of HCal to simulate the electronics readout material
1013 // Hardcoded thickness of 1cm and placement at 2 cm behind HCAL
1014 TGeoBBox* alHcalBox = new TGeoBBox("AlHCalBox", SizeXHCAL / 2.0, SizeYHCAL / 2.0, 0.5 / 2.0);
1015 TGeoVolume* volumeAlHcalBox = new TGeoVolume("volAlHcalBox", alHcalBox, gGeoManager->GetMedium(getMediumID(ID_ALUMINIUM)));
1016 volumeAlHcalBox->SetLineColor(kOrange);
1017 if (mGeometry->getInsertHCalReadoutMaterial()) {
1018 TVirtualMC::GetMC()->Gspos("volAlHcalBox", 9999, "FOCAL", 0.0, 0.0, +1.0 * mGeometry->getFOCALSizeZ() / 2.0 + 1.0, 0, "ONLY");
1019 }
1020 TGeoBBox* alUnderBox = new TGeoBBox("AlUnderBox", SizeXHCAL / 2.0, 0.5, mGeometry->getFOCALSizeZ() / 2.0 + 1.5);
1021 TGeoVolume* volumeAlUnderBox = new TGeoVolume("volAlUnderBox", alUnderBox, gGeoManager->GetMedium(getMediumID(ID_ALUMINIUM)));
1022 volumeAlUnderBox->SetLineColor(kOrange);
1023 if (mGeometry->getInsertHCalReadoutMaterial()) {
1024 TVirtualMC::GetMC()->Gspos("volAlUnderBox", 9999, "FOCAL", 0.0, -1.0 * mGeometry->getFOCALSizeY() / 2 - 10.5, 0.0, 0, "ONLY");
1025 }
1026
1027 TGeoVolume* volFOCAL = gGeoManager->GetVolume("FOCAL");
1028 volFOCAL->AddNode(volHCAL, 1, new TGeoTranslation(0, 0, mGeometry->getHCALCenterZ() - mGeometry->getFOCALSizeZ() / 2 + 0.01 + (mGeometry->getInsertFrontPadLayers() ? 2.0 : 0.0) - (mGeometry->getInsertHCalReadoutMaterial() ? 1.5 : 0.0))); // 0.01 to avoid overlap with ECAL
1029}
1030
1031//_____________________________________________________________________________
1033{
1034 // using boost::algorithm::contains; // only when string operations
1035 Geometry* geom = getGeometry();
1036
1037 // Int_t *idtmed = fIdtmed->GetArray() - 3599; //599 -> 3599
1038
1042
1044 double pars[4]; // this is EMSC Assembly
1045 pars[0] = geom->getTowerSizeX() / 2. + geom->getTowerGapSizeX() / 2.;
1046 pars[1] = geom->getTowerSizeY() / 2. + geom->getTowerGapSizeY() / 2.;
1047 // pars[2] = mGeometry->GetFOCALSizeZ() / 2;
1048 pars[2] = geom->getECALSizeZ() / 2;
1049 pars[3] = 0;
1050 // this shifts all the pixel layers to the center near the beampipe
1051 double pixshift = geom->getTowerSizeX() - (geom->getGlobalPixelWaferSizeX() * geom->getNumberOfPIXsInX());
1052
1053 bool splitDet = mGeometry->getDetectorOpeningRight() > 0.0 || mGeometry->getDetectorOpeningLeft() > 0.0;
1054
1055 float offset = pars[2];
1056 // gMC->Gsvolu("EMSC1", "BOX", idtmed[3698], pars, 4);//Left towers (pixels shifted right)
1057 // gMC->Gsvolu("EMSC2", "BOX", idtmed[3698], pars, 4);//Right towers (pixels shifted left)
1058
1059 TVirtualMC::GetMC()->Gsvolu("EMSC1", "BOX", getMediumID(ID_AIR), pars, 4); // Left towers (pixels shifted right)
1060 TVirtualMC::GetMC()->Gsvolu("EMSC2", "BOX", getMediumID(ID_AIR), pars, 4); // Right towers (pixels shifted left)
1061 // mSensitiveECALPad.push_back("EMSC1");
1062 // mSensitiveECALPad.push_back("EMSC2");
1063 mSensitive.push_back("EMSC1");
1064 mSensitive.push_back("EMSC2");
1065
1066 // const Composition *icomp = new Composition(); //to be removed
1067 // for(int i = 0; i < 20; i++){ // old
1068 // icomp = geom->getComposition(i, 0); // obsolete
1069
1070 // loop over geometry composition elements
1071 for (auto& icomp : mGeoCompositions) {
1072
1073 pars[0] = icomp->sizeX() / 2.;
1074 pars[1] = icomp->sizeY() / 2.;
1075 pars[2] = icomp->sizeZ() / 2.;
1076 pars[3] = 0;
1077
1078 if (icomp->material() == "PureW") {
1079 // TVirtualMC::GetMC()->Gsvolu("EW1", "BOX", idtmed[3599], pars, 4);
1080 TVirtualMC::GetMC()->Gsvolu("EW1", "BOX", getMediumID(ID_TUNGSTEN), pars, 4);
1081 // mSensitiveECALPad.push_back("EW1");
1082 mSensitive.push_back("EW1");
1083 gGeoManager->GetVolume("EW1")->SetLineColor(kBlue);
1084 TVirtualMC::GetMC()->Gspos("EW1", icomp->id() + 1, "EMSC1",
1085 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1086 TVirtualMC::GetMC()->Gspos("EW1", icomp->id() + 1, "EMSC2",
1087 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1088 }
1089 if (icomp->material() == "Alloy") {
1090 // TVirtualMC::GetMC()->Gsvolu("EW1", "BOX", idtmed[3604], pars, 4);
1091 TVirtualMC::GetMC()->Gsvolu("EW1", "BOX", getMediumID(ID_ALLOY), pars, 4);
1092 // mSensitiveECALPad.push_back("EW1");
1093 mSensitive.push_back("EW1");
1094 TVirtualMC::GetMC()->Gspos("EW1", icomp->id() + 1, "EMSC1",
1095 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1096 TVirtualMC::GetMC()->Gspos("EW1", icomp->id() + 1, "EMSC2",
1097 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1098 }
1099
1100 if (icomp->material() == "G10") {
1101 // TVirtualMC::GetMC()->Gsvolu("G10RO1", "BOX", idtmed[3601], pars, 4);
1102 TVirtualMC::GetMC()->Gsvolu("G10RO1", "BOX", getMediumID(ID_G10), pars, 4);
1103 // mSensitiveECALPad.push_back("G10RO1");
1104 mSensitive.push_back("G10RO1");
1105 gGeoManager->GetVolume("G10RO1")->SetLineColor(kGreen);
1106 TVirtualMC::GetMC()->Gspos("G10RO1", icomp->id() + 1, "EMSC1",
1107 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1108 TVirtualMC::GetMC()->Gspos("G10RO1", icomp->id() + 1, "EMSC2",
1109 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1110 }
1111
1112 if (icomp->material() == "Cu") {
1113 // TVirtualMC::GetMC()->Gsvolu("EWCU", "BOX", idtmed[3602], pars, 4);
1114 TVirtualMC::GetMC()->Gsvolu("EWCU", "BOX", getMediumID(ID_COPPER), pars, 4);
1115 // mSensitiveECALPad.push_back("EWCU");
1116 mSensitive.push_back("EWCU");
1117 gGeoManager->GetVolume("EWCU")->SetLineColor(kViolet);
1118 TVirtualMC::GetMC()->Gspos("EWCU", icomp->id() + 1, "EMSC1",
1119 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1120 TVirtualMC::GetMC()->Gspos("EWCU", icomp->id() + 1, "EMSC2",
1121 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1122 }
1123
1124 if (icomp->material() == "Air") {
1125 // TVirtualMC::GetMC()->Gsvolu("EWAIR1", "BOX", idtmed[3698], pars, 4);
1126 TVirtualMC::GetMC()->Gsvolu("EWAIR1", "BOX", getMediumID(ID_AIR), pars, 4);
1127 // mSensitiveECALPad.push_back("EWAIR1");
1128 mSensitive.push_back("EWAIR1");
1129 gGeoManager->GetVolume("EWAIR1")->SetLineColor(kGray);
1130 TVirtualMC::GetMC()->Gspos("EWAIR1", icomp->id() + 1, "EMSC1",
1131 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1132 TVirtualMC::GetMC()->Gspos("EWAIR1", icomp->id() + 1, "EMSC2",
1133 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1134 }
1135
1136 if (icomp->material() == "Ceramic") {
1137 // TVirtualMC::GetMC()->Gsvolu("EWAIR1", "BOX", idtmed[3607], pars, 4);
1138 TVirtualMC::GetMC()->Gsvolu("EWAIR1", "BOX", getMediumID(ID_CERAMIC), pars, 4);
1139 // mSensitiveECALPad.push_back("EWAIR1");
1140 mSensitive.push_back("EWAIR1");
1141 TVirtualMC::GetMC()->Gspos("EWAIR1", icomp->id() + 1, "EMSC1",
1142 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1143 TVirtualMC::GetMC()->Gspos("EWAIR1", icomp->id() + 1, "EMSC2",
1144 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1145 }
1146
1147 if (icomp->material() == "SiPad") {
1148 // TVirtualMC::GetMC()->Gsvolu("EWSIPAD1", "BOX", idtmed[3600], pars, 4);
1149 TVirtualMC::GetMC()->Gsvolu("EWSIPAD1", "BOX", getMediumID(ID_SIPAD), pars, 4);
1150 // mSensitiveECALPad.push_back("EWSIPAD1");
1151 mSensitive.push_back("EWSIPAD1");
1152 gGeoManager->GetVolume("EWSIPAD1")->SetLineColor(kOrange - 7);
1153 int number = (icomp->id()) + (icomp->stack() << 12) + (icomp->layer() << 16);
1154 // cout<<" pad : "<< icomp->material()<<" "<<number<<" x: "<< pars[0] << " y: " << pars[1] <<" Z coord: " << icomp->centerZ()-offset <<endl;
1155 TVirtualMC::GetMC()->Gspos("EWSIPAD1", number + 1, "EMSC1",
1156 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1157 TVirtualMC::GetMC()->Gspos("EWSIPAD1", number + 1, "EMSC2",
1158 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1159 }
1160
1161 // Pixels (sensitive layer)
1162 if (icomp->material() == "SiPix") {
1163 // TVirtualMC::GetMC()->Gsvolu("EWSIPIX1", "BOX", idtmed[3600], pars, 4);
1164 TVirtualMC::GetMC()->Gsvolu("EWSIPIX1", "BOX", getMediumID(ID_SIPIX), pars, 4);
1165 // mSensitiveECALPix.push_back("EWSIPIX1");
1166 mSensitive.push_back("EWSIPIX1");
1167 gGeoManager->GetVolume("EWSIPIX1")->SetLineColor(kPink);
1168
1169 int number = (icomp->id()) + (icomp->stack() << 12) + (icomp->layer() << 16);
1170 TVirtualMC::GetMC()->Gspos("EWSIPIX1", number + 1, "EMSC1",
1171 icomp->centerX() - geom->getGlobalPixelOffsetX() + pixshift, icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1172 TVirtualMC::GetMC()->Gspos("EWSIPIX1", number + 1, "EMSC2",
1173 icomp->centerX() + geom->getGlobalPixelOffsetX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1174 }
1175
1176 // Passive silicon
1177 if (icomp->material() == "Si") {
1178 // TVirtualMC::GetMC()->Gsvolu("EWSI1", "BOX", idtmed[3610], pars, 4);
1179 TVirtualMC::GetMC()->Gsvolu("EWSI1", "BOX", getMediumID(ID_SIINSENS), pars, 4);
1180 // mSensitiveECALPix.push_back("EWSI1");
1181 mSensitive.push_back("EWSI1");
1182 gGeoManager->GetVolume("EWSI1")->SetLineColor(kPink);
1183 TVirtualMC::GetMC()->Gspos("EWSI1", icomp->id() + 1, "EMSC1",
1184 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1185 TVirtualMC::GetMC()->Gspos("EWSI1", icomp->id() + 1, "EMSC2",
1186 icomp->centerX(), icomp->centerY(), icomp->centerZ() - offset, 0, "ONLY");
1187 }
1188 } // end of loop over composition elements
1189
1190 // Add the coldplates to each of the left and right towers
1191 TGeoBBox* coldPlateBox = new TGeoBBox("ColdPlateBox", geom->getTowerSizeX() / 2.0, geom->getTowerGapSizeY() / 2.0, geom->getECALSizeZ() / 2.0);
1192 TGeoVolume* volumeColdPlate = nullptr;
1193
1194 if (geom->getTowerGapMaterial() == "Cu") { // Copper
1195 // if (contains(geom->getTowerGapMaterial(), "Cu")) { // Copper
1196 volumeColdPlate = new TGeoVolume("volColdPlate", coldPlateBox, gGeoManager->GetMedium(getMediumID(ID_COPPER)));
1197 } else if (geom->getTowerGapMaterial() == "Al") { // Aluminium
1198 // else if (contains(geom->getTowerGapMaterial(), "Al")) { // Aluminium
1199 volumeColdPlate = new TGeoVolume("volColdPlate", coldPlateBox, gGeoManager->GetMedium(getMediumID(ID_ALUMINIUM)));
1200 } else {
1201 volumeColdPlate = new TGeoVolume("volColdPlate", coldPlateBox, gGeoManager->GetMedium(getMediumID(ID_AIR)));
1202 }
1203 // mSensitiveECALPad.push_back(volumeColdPlate->GetName());
1204 mSensitive.push_back(volumeColdPlate->GetName());
1205 volumeColdPlate->SetLineColor(kOrange);
1206 TVirtualMC::GetMC()->Gspos("volColdPlate", 1, "EMSC1", 0.0, geom->getTowerSizeY() / 2.0 + geom->getTowerGapSizeY() / 2.0, 0.0, 0, "ONLY");
1207 TVirtualMC::GetMC()->Gspos("volColdPlate", 1, "EMSC2", 0.0, geom->getTowerSizeY() / 2.0 + geom->getTowerGapSizeY() / 2.0, 0.0, 0, "ONLY");
1208
1209 // Place the towers in the ECAL
1210 // --- Place the ECAL in FOCAL
1211 float fcal_pars[4];
1212 fcal_pars[0] = (geom->getFOCALSizeX() + 2. * geom->getMiddleTowerOffset() + mGeometry->getDetectorOpeningRight() + mGeometry->getDetectorOpeningLeft()) / 2.;
1213 fcal_pars[1] = geom->getFOCALSizeY() / 2.;
1214 fcal_pars[2] = geom->getECALSizeZ() / 2.;
1215 fcal_pars[3] = 0.;
1216
1217 // TVirtualMC::GetMC()->Gsvolu("ECAL", "BOX", idtmed[3698], fcal_pars, 4);
1218 TVirtualMC::GetMC()->Gsvolu("ECAL", "BOX", getMediumID(ID_AIR), fcal_pars, 4);
1219 // mSensitiveECALPad.push_back("ECAL");
1220 mSensitive.push_back("ECAL");
1221
1222 // Create SiPad box for the two sensitive layers to be placed in front of ECAL
1223 TGeoBBox* siPadBox = new TGeoBBox("SiPadBox", geom->getTowerSizeX() / 2. + geom->getTowerGapSizeX() / 2.,
1224 geom->getTowerSizeY() / 2. + geom->getTowerGapSizeY() / 2., 0.03 / 2.0);
1225 TGeoVolume* volumeSiPad = new TGeoVolume("volSiPad", siPadBox, gGeoManager->GetMedium(getMediumID(ID_SIPAD)));
1226 volumeSiPad->SetLineColor(kOrange + 7);
1227 // mSensitiveECALPad.push_back(volumeSiPad->GetName());
1228 if (geom->getInsertFrontPadLayers()) {
1229 mSensitive.push_back(volumeSiPad->GetName());
1230 }
1231
1232 double xp, yp, zp;
1233 int itowerx, itowery;
1234 // int number = i + j * geom->getNumberOfTowersInX();
1235 for (int number = 0; number < geom->getNumberOfTowersInX() * geom->getNumberOfTowersInY(); number++) {
1236 itowerx = number % geom->getNumberOfTowersInX();
1237 itowery = number / geom->getNumberOfTowersInX();
1238 // const auto towerCenter = geom->getGeoTowerCenter(number); //only ECAL part, second parameter = -1 by default
1239 // xp = std::get<0>towerCenter;
1240 // std::tie(xp, yp, zp) = geom->getGeoTowerCenter(number);
1241 auto [xp, yp, zp] = geom->getGeoTowerCenter(number); // only ECAL part, second parameter = -1 by default
1242
1243 if (itowerx == 0) {
1244 if (splitDet) {
1245 xp -= geom->getDetectorOpeningLeft();
1246 }
1247
1248 TVirtualMC::GetMC()->Gspos("EMSC1", number + 1, "ECAL", xp, yp, 0, 0, "ONLY");
1249 // Add the SiPad front volumes directly under the FOCAL volume
1250 if (geom->getInsertFrontPadLayers()) {
1251 TVirtualMC::GetMC()->Gspos("volSiPad", -1 * (number + 1), "FOCAL", xp, yp, -1.0 * geom->getFOCALSizeZ() / 2.0, 0, "ONLY");
1252 // mSensitiveECALPad.push_back("volSiPad");
1253 mSensitive.push_back("volSiPad");
1254 TVirtualMC::GetMC()->Gspos("volSiPad", -1 * (geom->getNumberOfTowersInX() * geom->getNumberOfTowersInY() + number + 1), "FOCAL", xp - 0.5, yp + 0.5, -1.0 * geom->getFOCALSizeZ() / 2.0 + 1.0, 0, "ONLY");
1255 // mSensitiveECALPad.push_back("volSiPad");
1256 mSensitive.push_back("volSiPad");
1257 }
1258 }
1259 if (itowerx == 1) {
1260 if (splitDet) {
1261 xp += geom->getDetectorOpeningRight();
1262 }
1263
1264 TVirtualMC::GetMC()->Gspos("EMSC2", number + 1, "ECAL", xp, yp, 0, 0, "ONLY");
1265 // Add the SiPad front volumes directly under the FOCAL volume
1266 if (geom->getInsertFrontPadLayers()) {
1267 TVirtualMC::GetMC()->Gspos("volSiPad", -1 * (number + 1), "FOCAL", xp, yp, -1.0 * geom->getFOCALSizeZ() / 2.0, 0, "ONLY");
1268 // mSensitiveECALPad.push_back("volSiPad");
1269 mSensitive.push_back("volSiPad");
1270 TVirtualMC::GetMC()->Gspos("volSiPad", -1 * (geom->getNumberOfTowersInX() * geom->getNumberOfTowersInY() + number + 1), "FOCAL", xp + 0.5, yp + 0.5, -1.0 * geom->getFOCALSizeZ() / 2.0 + 1.0, 0, "ONLY");
1271 // mSensitiveECALPad.push_back("volSiPad");
1272 mSensitive.push_back("volSiPad");
1273 }
1274 }
1275 } // end of loop over ECAL towers (TowersInX x TowersInY)
1276
1277 TVirtualMC::GetMC()->Gspos("ECAL", 1, "FOCAL", 0, 0, geom->getECALCenterZ() - geom->getFOCALSizeZ() / 2.0 + (geom->getInsertFrontPadLayers() ? 2.0 : 0.0) - (geom->getInsertHCalReadoutMaterial() ? 1.5 : 0.0), 0, "ONLY");
1278}
1279
1281{
1282 mCurrentPrimaryID = fMC->GetStack()->GetCurrentTrackNumber();
1283 LOG(debug) << "Starting primary " << mCurrentPrimaryID << " with energy " << fMC->GetStack()->GetCurrentTrack()->Energy();
1284}
1285
1287{
1288 LOG(debug) << "Finishing primary " << mCurrentPrimaryID;
1289 // Resetting primary and parent ID
1290 mCurrentPrimaryID = -1;
1291}
1292
1294{
1295 LOG(debug) << "We are in sensitive volume " << v->GetName() << ": " << TVirtualMC::GetMC()->CurrentVolPath();
1296
1297 double eloss = TVirtualMC::GetMC()->Edep() * 1e9; // energy in eV (GeV->eV)
1298 if (eloss < DBL_EPSILON) {
1299 return false; // only process hits which actually deposit some energy in the FOCAL
1300 }
1301
1302 // In case of new parent track create new track reference
1303 auto o2stack = static_cast<o2::data::Stack*>(TVirtualMC::GetMC()->GetStack());
1304 if (!mCurrentSuperparent->mHasTrackReference) {
1305 float x, y, z, px, py, pz, e;
1306 TVirtualMC::GetMC()->TrackPosition(x, y, z);
1307 TVirtualMC::GetMC()->TrackMomentum(px, py, pz, e);
1308 o2::TrackReference trackref(x, y, z, px, py, pz, TVirtualMC::GetMC()->TrackLength(), TVirtualMC::GetMC()->TrackTime(), mCurrentParentID, GetDetId());
1309 o2stack->addTrackReference(trackref);
1310 mCurrentSuperparent->mHasTrackReference = true;
1311 }
1312
1313 float posX, posY, posZ;
1314 TVirtualMC::GetMC()->TrackPosition(posX, posY, posZ);
1315
1316 auto [indetector, col, row, layer, segment] = mGeometry->getVirtualInfo(posX, posY, posZ);
1317
1318 if (!indetector) {
1319 // particle outside the detector
1320 return true;
1321 }
1322
1323 auto currenthit = FindHit(mCurrentParentID, col, row, layer);
1324 if (!currenthit) {
1325 // Condition for new hit:
1326 // - Processing different partent track (parent track must be produced outside FOCAL)
1327 // - Inside different cell
1328 // - First track of the event
1329 Double_t time = TVirtualMC::GetMC()->TrackTime() * 1e9; // time in ns
1330 LOG(debug3) << "Adding new hit for parent " << mCurrentParentID << " and cell Col: " << col << " Row: " << row << " segment: " << segment;
1331
1333 AddHit(mCurrentParentID, mCurrentPrimaryID, mCurrentSuperparent->mEnergy, row * col + col, o2::focal::Hit::Subsystem_t::EPADS, math_utils::Point3D<float>(posX, posY, posZ), time, eloss);
1334 o2stack->addHit(GetDetId());
1335 } else {
1336 LOG(debug3) << "Adding energy to the current hit";
1337 currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + eloss);
1338 }
1339 return true;
1340}
1341
1343{
1344 LOG(debug) << "We are in sensitive volume " << v->GetName() << ": " << fMC->CurrentVolPath();
1345
1346 double eloss = fMC->Edep() * 1e9; // energy in eV (GeV->eV)
1347 if (eloss < DBL_EPSILON) {
1348 return false; // only process hits which actually deposit some energy in the FOCAL
1349 }
1350
1351 // In case of new parent track create new track reference
1352 auto o2stack = static_cast<o2::data::Stack*>(fMC->GetStack());
1353 if (!mCurrentSuperparent->mHasTrackReference) {
1354 float x, y, z, px, py, pz, e;
1355 fMC->TrackPosition(x, y, z);
1356 fMC->TrackMomentum(px, py, pz, e);
1357 o2::TrackReference trackref(x, y, z, px, py, pz, fMC->TrackLength(), fMC->TrackTime(), mCurrentParentID, GetDetId());
1358 o2stack->addTrackReference(trackref);
1359 mCurrentSuperparent->mHasTrackReference = true;
1360 }
1361
1362 float posX, posY, posZ;
1363 fMC->TrackPosition(posX, posY, posZ);
1364
1365 auto [indetector, col, row, layer, segment] = mGeometry->getVirtualInfo(posX, posY, posZ);
1366
1367 if (!indetector) {
1368 // particle outside the detector
1369 return true;
1370 }
1371
1372 auto currenthit = FindHit(mCurrentParentID, col, row, layer);
1373 if (!currenthit) {
1374 // Condition for new hit:
1375 // - Processing different partent track (parent track must be produced outside FOCAL)
1376 // - Inside different cell
1377 // - First track of the event
1378 Double_t time = fMC->TrackTime() * 1e9; // time in ns
1379 LOG(debug3) << "Adding new hit for parent " << mCurrentParentID << " and cell Col: " << col << " Row: " << row << " segment: " << segment;
1380
1382 AddHit(mCurrentParentID, mCurrentPrimaryID, mCurrentSuperparent->mEnergy, row * col + col, o2::focal::Hit::Subsystem_t::HCAL, math_utils::Point3D<float>(posX, posY, posZ), time, eloss);
1383 o2stack->addHit(GetDetId());
1384 } else {
1385 LOG(debug3) << "Adding energy to the current hit";
1386 currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + eloss);
1387 }
1388 return true;
1389}
1390
1392{
1393 LOG(debug) << "We are in sensitive volume " << v->GetName() << ": " << fMC->CurrentVolPath();
1394
1395 double eloss = fMC->Edep() * 1e9; // energy in eV (GeV->eV)
1396 if (eloss < DBL_EPSILON) {
1397 return false; // only process hits which actually deposit some energy in the FOCAL
1398 }
1399
1400 // In case of new parent track create new track reference
1401 auto o2stack = static_cast<o2::data::Stack*>(fMC->GetStack());
1402 if (!mCurrentSuperparent->mHasTrackReference) {
1403 float x, y, z, px, py, pz, e;
1404 fMC->TrackPosition(x, y, z);
1405 fMC->TrackMomentum(px, py, pz, e);
1406 o2::TrackReference trackref(x, y, z, px, py, pz, fMC->TrackLength(), fMC->TrackTime(), mCurrentParentID, GetDetId());
1407 o2stack->addTrackReference(trackref);
1408 mCurrentSuperparent->mHasTrackReference = true;
1409 }
1410
1411 float posX, posY, posZ;
1412 fMC->TrackPosition(posX, posY, posZ);
1413
1414 auto [indetector, col, row, layer, segment] = mGeometry->getVirtualInfo(posX, posY, posZ);
1415
1416 if (!indetector) {
1417 // particle outside the detector
1418 return true;
1419 }
1420
1421 auto currenthit = FindHit(mCurrentParentID, col, row, layer);
1422 if (!currenthit) {
1423 // Condition for new hit:
1424 // - Processing different partent track (parent track must be produced outside FOCAL)
1425 // - Inside different cell
1426 // - First track of the event
1427 Double_t time = fMC->TrackTime() * 1e9; // time in ns
1428 LOG(debug3) << "Adding new hit for parent " << mCurrentParentID << " and cell Col: " << col << " Row: " << row << " segment: " << segment;
1429
1431 AddHit(mCurrentParentID, mCurrentPrimaryID, mCurrentSuperparent->mEnergy, row * col + col, o2::focal::Hit::Subsystem_t::EPIXELS, math_utils::Point3D<float>(posX, posY, posZ), time, eloss);
1432 o2stack->addHit(GetDetId());
1433 } else {
1434 LOG(debug3) << "Adding energy to the current hit";
1435 currenthit->SetEnergyLoss(currenthit->GetEnergyLoss() + eloss);
1436 }
1437 return true;
1438}
Definition of the Stack class.
int16_t time
Definition RawEventData.h:4
int32_t i
uint16_t pos
Definition RawData.h:3
uint32_t col
Definition RawData.h:4
std::ostringstream debug
int registerSensitiveVolumeAndGetVolID(std::string const &name)
Definition Detector.cxx:190
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
int getMediumID(int imed) const
Definition Detector.h:135
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
FOCAL detector simulation.
Definition Detector.h:55
bool ProcessHitsEPad(FairVolume *v=nullptr)
Processing hit creation in the ECAL Pad sensitive volume.
TGeoVolumeAssembly * CreatePitchAssembly(double Lx=498.1, double Ly1=2.0, double Ly2=1.5, double Lz=1100.0, double hole_diameter=1.1, double hole_spacing=4.0, int nholes=124, double fiber_radius=0.5, std::string suffix="")
Definition Detector.cxx:697
void addAlignableVolumesHCAL() const
Definition Detector.cxx:407
~Detector() override
Destructor.
Definition Detector.cxx:65
void InitializeO2Detector() override
Initializing detector.
Definition Detector.cxx:81
bool ProcessHitsHCAL(FairVolume *v=nullptr)
Processing hit creation in the HCAL sensitive volume.
virtual void CreateHCALSandwich()
Definition Detector.cxx:910
void Register() override
register container with hits
Definition Detector.cxx:197
virtual void CreateHCALSheets()
Definition Detector.cxx:831
Hit * AddHit(int trackID, int primary, double initialEnergy, int detID, o2::focal::Hit::Subsystem_t subsystem, const math_utils::Point3D< float > &pos, double time, double energyloss)
Add FOCAL hit.
Definition Detector.cxx:163
void FinishPrimary() override
Finish current primary.
virtual void CreateHCALSpaghetti()
Definition Detector.cxx:504
void ConstructGeometry() override
Definition Detector.cxx:428
void addAlignableVolumesECAL() const
Definition Detector.cxx:418
bool ProcessHitsEPix(FairVolume *v=nullptr)
Processing hit creation in the ECAL Pixel sensitive volume.
virtual void addAlignableVolumes() const override
declare alignable volumes of detector
Definition Detector.cxx:395
void CreateMaterials()
Creating detector materials for the FOCAL detector.
Definition Detector.cxx:216
Hit * FindHit(int parentID, int col, int row, int layer)
Try to find hit with same cell and parent track ID.
Definition Detector.cxx:176
Bool_t ProcessHits(FairVolume *v=nullptr) final
Processing hit creation in the FOCAL sensitive volume.
Definition Detector.cxx:108
void EndOfEvent() final
Steps to be carried out at the end of the event.
Definition Detector.cxx:195
void BeginPrimary() override
Begin primaray.
void CreateECALGeometry()
Generate ECAL geometry.
Detector()=default
Dummy constructor.
Parent * AddSuperparent(int trackID, int pdg, double energy)
Add new superparent to the container.
Definition Detector.cxx:188
void Reset() final
Clean point collection.
Definition Detector.cxx:202
Geometry * getGeometry(std::string name="")
Get the FOCAL geometry desciption.
Definition Detector.cxx:70
double getTowerGapSizeY() const
Definition Geometry.h:133
double getFOCALSizeY() const
Definition Geometry.cxx:860
float getDetectorOpeningRight() const
Definition Geometry.h:146
int getNumberOfTowersInX() const
Definition Geometry.h:129
bool getInsertFrontPadLayers() const
Definition Geometry.h:141
HCALDesgin getHCALDesign() const
Definition Geometry.h:169
bool getInsertHCalReadoutMaterial() const
Definition Geometry.h:142
double getTowerGapSizeX() const
Definition Geometry.h:132
double getHCALSizeZ() const
Definition Geometry.cxx:899
int getNumberOfTowersInY() const
Definition Geometry.h:130
float getDetectorOpeningLeft() const
Definition Geometry.h:147
double getECALSizeZ() const
Definition Geometry.cxx:878
int getHCALTowersInY() const
Definition Geometry.h:109
float getHCALPitchSize() const
Definition Geometry.h:143
std::tuple< bool, int, int, int, int > getVirtualInfo(double x, double y, double z) const
double getGlobalPixelWaferSizeX() const
Definition Geometry.h:135
double getFOCALSizeZ() const
Definition Geometry.cxx:866
float getHCALTowerSize() const
Definition Geometry.h:107
static Geometry * getInstance()
Definition Geometry.cxx:41
int getHCALTowersInX() const
Definition Geometry.h:108
double getHCALCenterZ() const
Definition Geometry.cxx:910
std::tuple< double, double, double > getGeoTowerCenter(int tower, int segment=-1) const
this gives global position of the center of tower
Definition Geometry.cxx:694
double getTowerSizeX() const
Definition Geometry.cxx:840
std::vector< const Composition * > getFOCALMicroModule(int layer) const
Definition Geometry.cxx:673
std::string_view getTowerGapMaterial() const
Definition Geometry.h:151
float getHCALBeamPipeHoleSize() const
Definition Geometry.h:144
double getFOCALZ0() const
Definition Geometry.h:128
float getMiddleTowerOffset() const
Definition Geometry.h:140
double getFOCALSizeX() const
Definition Geometry.cxx:854
double getTowerSizeY() const
Definition Geometry.cxx:847
int getNumberOfPIXsInX() const
Definition Geometry.h:105
double getECALCenterZ() const
Definition Geometry.cxx:890
double getGlobalPixelOffsetX() const
Definition Geometry.h:138
Common FOCAL hit class for the detector simulation.
Definition Hit.h:28
Subsystem_t
Subsystem index of the Hit.
Definition Hit.h:32
static ShmManager & Instance()
Definition ShmManager.h:61
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint segment
Definition glcorearb.h:4945
GLuint entry
Definition glcorearb.h:5735
const GLdouble * v
Definition glcorearb.h:832
GLuint const GLchar * name
Definition glcorearb.h:781
GLint y
Definition glcorearb.h:270
GLboolean * data
Definition glcorearb.h:298
GLintptr offset
Definition glcorearb.h:660
GLenum GLuint GLint GLint layer
Definition glcorearb.h:1310
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
D const SVectorGPU< T, D > & rhs
Definition SMatrixGPU.h:193
void freeSimVector(std::vector< T > *ptr)
std::vector< T > * createSimVector()
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Common utility functions.
Mapped information of a channel.
Definition Hit.h:41
Information about superparent (particle entering any FOCAL volume)
Definition Detector.h:33
double mEnergy
Total energy.
Definition Detector.h:35
bool mHasTrackReference
Flag indicating whether parent has a track reference.
Definition Detector.h:36
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
std::vector< int > row