Project
Loading...
Searching...
No Matches
MaterialManager.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
14
17#include "TVirtualMC.h"
18#include "TString.h" // for TString
19#include <TGeoMedium.h>
20#include <TGeoManager.h>
21#include <TList.h>
22#include <iostream>
23#include <utility>
24#include <fairlogger/Logger.h>
25#ifdef NDEBUG
26#undef NDEBUG
27#endif
28#include <cassert>
29#include <set>
30#include "rapidjson/document.h"
31#include "rapidjson/istreamwrapper.h"
32#include "rapidjson/ostreamwrapper.h"
33#include "rapidjson/prettywriter.h"
34#include <algorithm>
35#include <SimConfig/SimParams.h>
36
37using namespace o2::base;
38namespace rj = rapidjson;
39
40namespace
41{
43template <typename K, typename V>
44void writeSingleJSONParamBatch(std::unordered_map<K, const char*> const& idToName, std::map<K, V> const& valMap, V defaultValue, rapidjson::Value& parent, rapidjson::Document::AllocatorType& a)
45{
46 for (auto& itName : idToName) {
47 auto itVal = valMap.find(itName.first);
48 if (itVal != valMap.end()) {
49 parent.AddMember(rj::Value(itName.second, std::strlen(itName.second), a), rj::Value(itVal->second), a);
50 continue;
51 }
52 parent.AddMember(rj::Value(itName.second, std::strlen(itName.second), a), rj::Value(defaultValue), a);
53 }
54}
55
57static constexpr const char* jsonKeyID = "local_id";
58static constexpr const char* jsonKeyIDGlobal = "global_id";
59static constexpr const char* jsonKeyDefault = "default";
60static constexpr const char* jsonKeyCuts = "cuts";
61static constexpr const char* jsonKeyProcesses = "processes";
62static constexpr const char* jsonKeyEnableSpecialCuts = "enableSpecialCuts";
63static constexpr const char* jsonKeyEnableSpecialProcesses = "enableSpecialProcesses";
64} // namespace
65
66const std::unordered_map<EProc, const char*> MaterialManager::mProcessIDToName = {
67 {EProc::kPAIR, "PAIR"},
68 {EProc::kCOMP, "COMP"},
69 {EProc::kPHOT, "PHOT"},
70 {EProc::kPFIS, "PFIS"},
71 {EProc::kDRAY, "DRAY"},
72 {EProc::kANNI, "ANNI"},
73 {EProc::kBREM, "BREM"},
74 {EProc::kHADR, "HADR"},
75 {EProc::kMUNU, "MUNU"},
76 {EProc::kDCAY, "DCAY"},
77 {EProc::kLOSS, "LOSS"},
78 {EProc::kMULS, "MULS"},
79 {EProc::kCKOV, "CKOV"},
80 {EProc::kRAYL, "RAYL"},
81 {EProc::kLABS, "LABS"}};
82
83const std::unordered_map<ECut, const char*> MaterialManager::mCutIDToName = {
84 {ECut::kCUTGAM, "CUTGAM"},
85 {ECut::kCUTELE, "CUTELE"},
86 {ECut::kCUTNEU, "CUTNEU"},
87 {ECut::kCUTHAD, "CUTHAD"},
88 {ECut::kCUTMUO, "CUTMUO"},
89 {ECut::kBCUTE, "BCUTE"},
90 {ECut::kBCUTM, "BCUTM"},
91 {ECut::kDCUTE, "DCUTE"},
92 {ECut::kDCUTM, "DCUTM"},
93 {ECut::kPPCUTM, "PPCUTM"},
94 {ECut::kTOFMAX, "TOFMAX"}};
95
96// Constructing a map between module names and local material density values
97void MaterialManager::initDensityMap()
98{
100 if (globalDensityFactor < 0) {
101 LOG(fatal) << "Negative value "
102 << globalDensityFactor
103 << " found for global material density!\n";
104 }
105 std::string token;
106 std::istringstream input(
107 o2::conf::SimMaterialParams::Instance().localDensityFactor);
108 std::vector<std::string> inputModuleNames;
109 std::vector<std::string> inputDensityValues;
110 while (std::getline(input, token, ',')) {
111 std::size_t pos = token.find(':');
112 inputModuleNames.push_back(token.substr(0, pos));
113 inputDensityValues.push_back(token.substr(pos + 1));
114 }
115 for (std::size_t i = 0; i < inputModuleNames.size(); i++) {
116 if (std::stof(inputDensityValues[i]) < 0) {
117 LOG(fatal) << "Negative value " << std::stof(inputDensityValues[i])
118 << " found for material density in module "
119 << inputModuleNames[i] << "!\n";
120 }
121 mDensityMap[inputModuleNames[i]] = std::stof(inputDensityValues[i]);
122 }
123 mDensityMapInitialized = true;
124}
125
126float MaterialManager::getDensity(std::string const& modname, std::string const& matname)
127{
128 // This function returns the final density for a material of name matname inside module modname.
129 // The priority is
130 // - return density for a specific module + material if it exists in the lookup
131 // - return density for the module if it exists in the the lookup
132 // - return global density factor
133
134 auto debug = getenv("O2SIM_MATMGR_LOCALDENSITY_DEBUG");
135
136 if (!mDensityMapInitialized) {
137 initDensityMap();
138 }
139 // density on final material level
140 // (this works by a name lookup of pair "modname/matname")
141 std::string lookupstring = modname + "/" + matname;
142 auto iter = mDensityMap.find(lookupstring);
143 if (iter != mDensityMap.end()) {
144 if (debug) {
145 LOG(info) << "MatManager - " << modname << "/" << matname << " : applying density " << iter->second << " from material match";
146 }
147 return iter->second;
148 }
149 // density on module level
150 iter = mDensityMap.find(modname);
151 if (iter != mDensityMap.end()) {
152 if (debug) {
153 LOG(info) << "MatManager - " << modname << "/" << matname << " : applying density " << iter->second << " from module match";
154 }
155 return iter->second;
156 }
157 // global factor
159 if (debug && global != 1.0) {
160 LOG(info) << "MatManager - " << modname << "/" << matname << " : applying global density " << iter->second;
161 }
162 return global;
163}
164
165void MaterialManager::Material(const char* modname, Int_t imat, const char* name, Float_t a, Float_t z, Float_t dens,
166 Float_t radl, Float_t absl, Float_t* buf, Int_t nwbuf)
167{
168 TString uniquename = modname;
169 auto densityFactor = getDensity(modname, name);
170
171 uniquename.Append("_");
172 uniquename.Append(name);
173 if (TVirtualMC::GetMC()) {
174 // Check this!!!
175 int kmat = -1;
176 TVirtualMC::GetMC()->Material(kmat, uniquename.Data(), a, z,
177 dens * densityFactor, radl, absl, buf, nwbuf);
178 mMaterialMap[modname][imat] = kmat;
179 insertMaterialName(uniquename.Data(), kmat);
180 } else {
181 auto uid = gGeoManager->GetListOfMaterials()->GetSize();
182 auto mat = gGeoManager->Material(uniquename.Data(), a, z,
183 dens * densityFactor, uid, radl, absl);
184 mMaterialMap[modname][imat] = uid;
185 insertMaterialName(uniquename.Data(), uid);
186 }
187}
188
201void MaterialManager::Mixture(const char* modname, Int_t imat, const char* name, Float_t* a, Float_t* z, Float_t dens,
202 Int_t nlmat, Float_t* wmat)
203{
204 TString uniquename = modname;
205 auto densityFactor = getDensity(modname, name);
206 uniquename.Append("_");
207 uniquename.Append(name);
208
209 if (TVirtualMC::GetMC()) {
210 // Check this!!!
211 int kmat = -1;
212 TVirtualMC::GetMC()->Mixture(kmat, uniquename.Data(), a, z,
213 dens * densityFactor, nlmat, wmat);
214 mMaterialMap[modname][imat] = kmat;
215 insertMaterialName(uniquename.Data(), kmat);
216
217 } else {
218 auto uid = gGeoManager->GetListOfMaterials()->GetSize();
219 if (nlmat < 0) {
220 nlmat = -nlmat;
221 Double_t amol = 0;
222 Int_t i;
223 for (i = 0; i < nlmat; i++) {
224 amol += a[i] * wmat[i];
225 }
226 for (i = 0; i < nlmat; i++) {
227 wmat[i] *= a[i] / amol;
228 }
229 }
230 auto mix = gGeoManager->Mixture(uniquename.Data(), a, z,
231 dens * densityFactor, nlmat, wmat, uid);
232 mMaterialMap[modname][imat] = uid;
233 insertMaterialName(uniquename.Data(), uid);
234 }
235}
236
237void MaterialManager::Medium(const char* modname, Int_t numed, const char* name, Int_t nmat, Int_t isvol, Int_t ifield,
238 Float_t fieldm, Float_t tmaxfd, Float_t stemax, Float_t deemax, Float_t epsil,
239 Float_t stmin, Float_t* ubuf, Int_t nbuf)
240{
241 TString uniquename = modname;
242 uniquename.Append("_");
243 uniquename.Append(name);
244
245 if (TVirtualMC::GetMC()) {
246 // Check this!!!
247 int kmed = -1;
248 const int kmat = getMaterialID(modname, nmat);
249 TVirtualMC::GetMC()->Medium(kmed, uniquename.Data(), kmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil,
250 stmin, ubuf, nbuf);
251 mMediumMap[modname][numed] = kmed;
252 insertMediumName(uniquename.Data(), kmed);
253 insertTGeoMedium(modname, numed);
254 } else {
255 auto uid = gGeoManager->GetListOfMedia()->GetSize();
256 auto med = gGeoManager->Medium(uniquename.Data(), uid, getMaterialID(modname, nmat), isvol, ifield, fieldm, tmaxfd,
257 stemax, deemax, epsil, stmin);
258 mMediumMap[modname][numed] = uid;
259 insertMediumName(uniquename.Data(), uid);
260 insertTGeoMedium(modname, numed);
261 }
262}
263
264void MaterialManager::Processes(ESpecial special, int globalindex,
265 const std::initializer_list<std::pair<EProc, int>>& parIDValMap)
266{
267 for (auto& m : parIDValMap) {
268 Process(special, globalindex, m.first, m.second);
269 }
270}
271
272void MaterialManager::Cuts(ESpecial special, int globalindex,
273 const std::initializer_list<std::pair<ECut, Float_t>>& parIDValMap)
274{
275 for (auto& m : parIDValMap) {
276 Cut(special, globalindex, m.first, m.second);
277 }
278}
279
280void MaterialManager::Cut(ESpecial special, int globalindex, ECut cut, Float_t val)
281{
282 // this check is needed, in principal only for G3, otherwise SegFault
283 if (val < 0.) {
284 return;
285 }
286 // if low energy neutron transport is requested setting kCUTNEU will set to 0.005eV
287 if (mLowNeut && cut == ECut::kCUTNEU) {
288 LOG(info) << "Due to low energy neutrons, neutron cut value " << val << " discarded and reset to 5e-12";
289 val = 5.e-12;
290 }
291
292 auto it = mCutIDToName.find(cut);
293 if (it == mCutIDToName.end()) {
294 return;
295 }
296 if (special == ESpecial::kFALSE) {
297 auto ins = mDefaultCutMap.insert({cut, val});
298 if (ins.second) {
299 TVirtualMC::GetMC()->SetCut(it->second, val);
300 }
302 } else if (mApplySpecialCuts) {
303 auto ins = mMediumCutMap[globalindex].insert({cut, val});
304 if (ins.second) {
305 TVirtualMC::GetMC()->Gstpar(globalindex, it->second, val);
306 }
307 }
308}
309
310void MaterialManager::Process(ESpecial special, int globalindex, EProc process, int val)
311{
312 // this check is needed, in principal only for G3, otherwise SegFault
313 if (val < 0) {
314 return;
315 }
316 auto it = mProcessIDToName.find(process);
317 if (it == mProcessIDToName.end()) {
318 return;
319 }
320 if (special == ESpecial::kFALSE) {
321 auto ins = mDefaultProcessMap.insert({process, val});
322 if (ins.second) {
323 TVirtualMC::GetMC()->SetProcess(it->second, val);
324 }
326 } else if (mApplySpecialProcesses) {
327 auto ins = mMediumProcessMap[globalindex].insert({process, val});
328 if (ins.second) {
329 TVirtualMC::GetMC()->Gstpar(globalindex, it->second, val);
330 }
331 }
332}
333
335{
336 for (auto& p : mMaterialMap) {
337 auto name = p.first;
338 std::cout << "Materials for key " << name << "\n";
339 for (auto& e : p.second) {
340 std::cout << "internal id " << e.first << " to " << e.second << "\n";
341 }
342 }
343}
344
346{
347 for (auto& p : mMediumMap) {
348 auto name = p.first;
349 std::cout << "Tracking media for key " << name << "\n";
350 for (auto& e : p.second) {
351 auto medname = getMediumNameFromMediumID(e.second);
352 std::cout << medname << " ";
353 std::cout << "internal id " << e.first << " to " << e.second << "\n";
354 }
355 }
356}
357
359{
360 stream << "Summary of process settings per media.\n";
361 stream << "-- Default process settings:\n";
362 for (auto& p : mDefaultProcessMap) {
363 auto it = mProcessIDToName.find(p.first);
364 if (it != mProcessIDToName.end()) {
365 stream << "\t" << it->second << " = " << p.second << "\n";
366 }
367 }
368 if (mApplySpecialProcesses && mMediumProcessMap.size() > 0) {
369 stream << "-- Custom process settings for single media:\n";
370 for (auto& m : mMediumProcessMap) {
371 stream << "Global medium ID " << m.first << " (module = " << getModuleFromMediumID(m.first)
372 << ", medium name = " << getMediumNameFromMediumID(m.first) << "):\n";
373 for (auto& p : m.second) {
374 auto it = mProcessIDToName.find(p.first);
375 if (it != mProcessIDToName.end()) {
376 stream << "\t" << it->second << " = " << p.second << "\n";
377 }
378 }
379 }
380 }
381}
382
383void MaterialManager::printCuts(std::ostream& stream) const
384{
385 stream << "Summary of cut settings per media.\n";
386 stream << "-- Default cut settings:\n";
387 for (auto& c : mDefaultCutMap) {
388 auto it = mCutIDToName.find(c.first);
389 if (it != mCutIDToName.end()) {
390 stream << "\t" << it->second << " = " << c.second << "\n";
391 }
392 }
393 if (mApplySpecialCuts && mMediumCutMap.size() > 0) {
394 stream << "-- Custom cut settings for single media:\n";
395 for (auto& m : mMediumCutMap) {
396 stream << "Global medium ID " << m.first << " (module = " << getModuleFromMediumID(m.first)
397 << ", medium name = " << getMediumNameFromMediumID(m.first) << "):\n";
398 for (auto& c : m.second) {
399 auto it = mCutIDToName.find(c.first);
400 if (it != mCutIDToName.end()) {
401 stream << "\t" << it->second << " = " << c.second << "\n";
402 }
403 }
404 }
405 }
406}
407
408// insert material name
409void MaterialManager::insertMaterialName(const char* uniquename, int index)
410{
411 assert(mMaterialNameToGlobalIndexMap.find(uniquename) == mMaterialNameToGlobalIndexMap.end());
412 mMaterialNameToGlobalIndexMap[uniquename] = index;
413}
414
415// inserts the last
416void MaterialManager::insertTGeoMedium(std::string modname, int localindex)
417{
418 auto p = std::make_pair(modname, localindex);
419 assert(mTGeoMediumMap.find(p) == mTGeoMediumMap.end());
420 auto list = gGeoManager->GetListOfMedia();
421 mTGeoMediumMap[p] = (TGeoMedium*)list->At(list->GetEntries() - 1);
422
423 LOG(debug) << "mapping " << modname << " " << localindex << " to " << mTGeoMediumMap[p]->GetName();
424}
425
426void MaterialManager::insertMediumName(const char* uniquename, int index)
427{
428 assert(mMediumNameToGlobalIndexMap.find(uniquename) == mMediumNameToGlobalIndexMap.end());
429 mMediumNameToGlobalIndexMap[uniquename] = index;
430}
431
432// find TGeoMedium instance given a detector prefix and a local (vmc) medium index
433TGeoMedium* MaterialManager::getTGeoMedium(std::string const& modname, int localindex)
434{
435 auto p = std::make_pair(modname, localindex);
436 auto iter = mTGeoMediumMap.find(p);
437 if (iter == mTGeoMediumMap.end()) {
438 LOG(warning) << "No medium registered for " << modname << " index " << localindex << "\n";
439 return nullptr;
440 }
441 return iter->second;
442}
443
444// find TGeoMedium instance given the full medium name
445// mainly forwards directly to TGeoManager and raises a warning if medium is nullptr
446TGeoMedium* MaterialManager::getTGeoMedium(const char* mediumname)
447{
448 auto med = gGeoManager->GetMedium(mediumname);
449 assert(med != nullptr);
450 return med;
451}
452
454{
455 const std::string filenameIn = filename.empty() ? o2::MaterialManagerParam::Instance().inputFile : filename;
456 if (filenameIn.empty()) {
457 return;
458 }
459 std::ifstream is(filenameIn);
460 if (!is.is_open()) {
461 LOG(fatal) << "Cannot open MC cuts/processes file " << filenameIn;
462 return;
463 }
464 auto digestCutsFromJSON = [this](int globalindex, rj::Value& cuts) {
465 auto special = globalindex < 0 ? ESpecial::kFALSE : ESpecial::kTRUE;
466 for (auto& cut : cuts.GetObject()) {
467 auto name = cut.name.GetString();
468 bool found = false;
469 for (auto& cn : mCutIDToName) {
470 if (std::strcmp(name, cn.second) == 0) {
471 Cut(special, globalindex, cn.first, cut.value.GetFloat());
472 found = true;
473 }
474 }
475 if (!found) {
476 LOG(warn) << "Unknown cut parameter " << name;
477 }
478 }
479 };
480 auto digestProcessesFromJSON = [this](int globalindex, rj::Value& processes) {
481 auto special = globalindex < 0 ? ESpecial::kFALSE : ESpecial::kTRUE;
482 for (auto& proc : processes.GetObject()) {
483 auto name = proc.name.GetString();
484 for (auto& pn : mProcessIDToName) {
485 if (std::strcmp(name, pn.second) == 0) {
486 Process(special, globalindex, pn.first, proc.value.GetInt());
487 }
488 }
489 }
490 };
491
492 rj::IStreamWrapper isw(is);
493 rj::Document d;
494 d.ParseStream(isw);
495
496 if (special == ESpecial::kFALSE && d.HasMember(jsonKeyDefault)) {
497 // defaults
498 auto& defaultParams = d[jsonKeyDefault];
499 if (defaultParams.HasMember(jsonKeyCuts)) {
500 digestCutsFromJSON(-1, defaultParams[jsonKeyCuts]);
501 }
502 if (defaultParams.HasMember(jsonKeyProcesses)) {
503 digestProcessesFromJSON(-1, defaultParams[jsonKeyProcesses]);
504 }
505 } else if (special == ESpecial::kTRUE) {
506 // read whether to apply special cuts and processes at all
507 if (d.HasMember(jsonKeyEnableSpecialCuts)) {
508 enableSpecialCuts(d[jsonKeyEnableSpecialCuts].GetBool());
509 }
510 if (d.HasMember(jsonKeyEnableSpecialProcesses)) {
511 enableSpecialProcesses(d[jsonKeyEnableSpecialProcesses].GetBool());
512 }
513 // special
514 for (auto& m : d.GetObject()) {
515 if (m.name.GetString()[0] == '\0' || !m.value.IsArray()) {
516 // do not parse anything with empty key, these at the most meant to be comments
517 continue;
518 }
519 for (auto& batch : m.value.GetArray()) {
520 if (std::strcmp(m.name.GetString(), jsonKeyDefault) == 0) {
521 // don't do defaults here
522 continue;
523 }
524 // set via their global indices
525 auto index = getMediumID(m.name.GetString(), batch[jsonKeyID].GetInt());
526 if (index < 0) {
527 continue;
528 }
529 if (batch.HasMember(jsonKeyCuts)) {
530 digestCutsFromJSON(index, batch[jsonKeyCuts]);
531 }
532 if (batch.HasMember(jsonKeyProcesses)) {
533 digestProcessesFromJSON(index, batch[jsonKeyProcesses]);
534 }
535 }
536 }
537 }
538}
539
541{
542 const std::string filenameOut = filename.empty() ? o2::MaterialManagerParam::Instance().outputFile : filename;
543 if (filenameOut.empty()) {
544 return;
545 }
546
547 // write parameters as global AND module specific
548 std::ofstream os(filenameOut);
549 if (!os.is_open()) {
550 LOG(error) << "Cannot create file " << filenameOut;
551 return;
552 }
553
554 rj::Document d;
555 rj::Document::AllocatorType& a = d.GetAllocator();
556 d.SetObject();
557
558 // add each local medium with params per module
559 for (auto& itMed : mMediumMap) {
560 // prepare array for module
561 rj::Value toAdd(rj::kArrayType);
562 // extract each medium's local and global index
563 for (auto& locToGlob : itMed.second) {
564 auto globalindex = locToGlob.second;
565 auto itCut = mMediumCutMap.find(globalindex);
566 auto itProc = mMediumProcessMap.find(globalindex);
567 // prepare a batch summarising localID, globaldID, cuts and processes
568 rj::Value oLoc(rj::kObjectType);
569 // IDs
570 oLoc.AddMember(rj::Value(jsonKeyID, std::strlen(jsonKeyID), a), rj::Value(locToGlob.first), a);
571 oLoc.AddMember(rj::Value(jsonKeyIDGlobal, std::strlen(jsonKeyIDGlobal)), rj::Value(locToGlob.second), a);
572 // add medium and material name
573 auto mediumIt = mTGeoMediumMap.find({itMed.first, locToGlob.first});
574 const char* medName = mediumIt->second->GetName();
575 const char* matName = mediumIt->second->GetMaterial()->GetName();
576 // not using variables for key names cause they are only written for info but not read
577 oLoc.AddMember(rj::Value("medium_name", 11, a), rj::Value(medName, std::strlen(medName), a), a);
578 oLoc.AddMember(rj::Value("material_name", 13, a), rj::Value(matName, std::strlen(matName), a), a);
579 // prepare for cuts
580 if (itCut != mMediumCutMap.end()) {
581 rj::Value cutMap(rj::kObjectType);
582 writeSingleJSONParamBatch(mCutIDToName, itCut->second, -1.f, cutMap, a);
583 oLoc.AddMember(rj::Value(jsonKeyCuts, std::strlen(jsonKeyCuts), a), cutMap, a);
584 }
585 // prepare for processes
586 if (itProc != mMediumProcessMap.end()) {
587 rj::Value procMap(rj::kObjectType);
588 writeSingleJSONParamBatch(mProcessIDToName, itProc->second, -1, procMap, a);
589 oLoc.AddMember(rj::Value(jsonKeyProcesses, std::strlen(jsonKeyProcesses), a), procMap, a);
590 }
591 // append this medium to module array
592 toAdd.PushBack(oLoc, a);
593 }
594 // append the entire module array
595 d.AddMember(rj::Value(itMed.first.c_str(), itMed.first.size(), a), toAdd, a);
596 }
597 // also add default parameters
598 rj::Value cutMapDef(rj::kObjectType);
599 rj::Value procMapDef(rj::kObjectType);
600 writeSingleJSONParamBatch(mCutIDToName, mDefaultCutMap, -1.f, cutMapDef, a);
601 writeSingleJSONParamBatch(mProcessIDToName, mDefaultProcessMap, -1, procMapDef, a);
602 rj::Value defaultParams(rj::kObjectType);
603 defaultParams.AddMember(rj::Value(jsonKeyCuts, std::strlen(jsonKeyCuts), a), cutMapDef, a);
604 defaultParams.AddMember(rj::Value(jsonKeyProcesses, std::strlen(jsonKeyProcesses), a), procMapDef, a);
605 d.AddMember(rj::Value(jsonKeyDefault, std::strlen(jsonKeyDefault), a), defaultParams, a);
606
607 d.AddMember(rj::Value(jsonKeyEnableSpecialCuts, std::strlen(jsonKeyEnableSpecialCuts), a), rj::Value(mApplySpecialCuts), a);
608 d.AddMember(rj::Value(jsonKeyEnableSpecialProcesses, std::strlen(jsonKeyEnableSpecialProcesses), a), rj::Value(mApplySpecialProcesses), a);
609 // now write to file
610 rj::OStreamWrapper osw(os);
611 rj::PrettyWriter<rj::OStreamWrapper> writer(osw);
612 writer.SetIndent(' ', 2);
613 d.Accept(writer);
614}
615
616void MaterialManager::loadCutsAndProcessesFromFile(const char* modname, const char* filename)
617{
618 // Implementation of a method to set cuts and processes as done in AliRoot.
619 // The file is expected to contain columns of the form
620 // MODNAME LOCALMEDIUMID CUT0 ... CUT9 FLAG0 ... FLAG11
621 // where cuts and flags correspond to keys denoted in cutnames and procnames below.
622
623 const int NCUTS = 10; // number of cut columns expected in file
624 const int NFLAGS = 12; // number of process flag columns expected in file
625
627 using namespace o2::base;
628 ECut cutnames[NCUTS] = {ECut::kCUTGAM,
638
640 // missing STRA for the moment
641 EProc procnames[NFLAGS - 1] = {EProc::kANNI,
652
653 std::ifstream cutfile(filename);
654
655 if (!cutfile.is_open()) {
656 LOG(warn) << "File " << filename << " does not exist; Cannot apply cuts";
657 return;
658 }
659
660 // reading from file
661 float cut[NCUTS]; // to store cut values
662 int flag[NFLAGS]; // to store flags
663 int itmed, iret;
664 char line[256];
665 char detName[7];
666
667 while (cutfile.getline(line, 256)) {
668 // Initialise cuts and flags for this line
669 for (int i = 0; i < NCUTS; i++) {
670 cut[i] = -99;
671 }
672 for (int i = 0; i < NFLAGS; i++) {
673 flag[i] = -99;
674 }
675 if (strlen(line) == 0) {
676 continue;
677 }
678 // ignore comments marked by *
679 if (line[0] == '*') {
680 continue;
681 }
682 // Read the numbers
683 iret = sscanf(line, "%6s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d %d",
684 detName, &itmed, &cut[0], &cut[1], &cut[2], &cut[3], &cut[4], &cut[5], &cut[6], &cut[7], &cut[8],
685 &cut[9], &flag[0], &flag[1], &flag[2], &flag[3], &flag[4], &flag[5], &flag[6], &flag[7],
686 &flag[8], &flag[9], &flag[10], &flag[11]);
687 if (!iret) {
688 // nothing read
689 continue;
690 }
691
692 // apply cuts via material manager interface
693 for (int i = 0; i < NCUTS; ++i) {
694 if (cut[i] >= 0.) {
695 SpecialCut(detName, itmed, cutnames[i], cut[i]);
696 }
697 }
698
699 // apply process flags
700 for (int i = 0; i < NFLAGS - 1; ++i) {
701 if (flag[i] >= 0) {
702 SpecialProcess(detName, itmed, procnames[i], flag[i]);
703 }
704 }
705 } // end loop over lines
706}
707
711void MaterialManager::SpecialCuts(const char* modname, int localindex,
712 const std::initializer_list<std::pair<ECut, Float_t>>& parIDValMap)
713{
714 int globalindex = getMediumID(modname, localindex);
715 if (globalindex != -1) {
716 Cuts(ESpecial::kTRUE, globalindex, parIDValMap);
717 }
718}
719
720void MaterialManager::SpecialCut(const char* modname, int localindex, ECut parID, Float_t val)
721{
722 int globalindex = getMediumID(modname, localindex);
723 if (globalindex != -1) {
724 Cut(ESpecial::kTRUE, globalindex, parID, val);
725 } else {
726 LOG(warn) << "SpecialCut: NO GLOBALINDEX FOUND FOR " << modname << " " << localindex;
727 }
728}
729
731void MaterialManager::SpecialProcess(const char* modname, int localindex, EProc parID, int val)
732{
733 int globalindex = getMediumID(modname, localindex);
734 if (globalindex != -1) {
735 Process(ESpecial::kTRUE, globalindex, parID, val);
736 } else {
737 LOG(warn) << "SpecialProcess: NO GLOBALINDEX FOUND FOR " << modname << " " << localindex;
738 }
739}
740
741int MaterialManager::getMaterialID(const char* modname, int imat) const
742{
743 auto lookupiter = mMaterialMap.find(modname);
744 if (lookupiter == mMaterialMap.end()) {
745 return -1;
746 }
747 auto lookup = lookupiter->second;
748
749 auto iter = lookup.find(imat);
750 if (iter != lookup.end()) {
751 return iter->second;
752 }
753 return -1;
754}
755
756int MaterialManager::getMediumID(const char* modname, int imed) const
757{
758 auto lookupiter = mMediumMap.find(modname);
759 if (lookupiter == mMediumMap.end()) {
760 return -1;
761 }
762 auto lookup = lookupiter->second;
763
764 auto iter = lookup.find(imed);
765 if (iter != lookup.end()) {
766 return iter->second;
767 }
768 return -1;
769}
770
771// fill the medium index mapping into a standard vector
772// the vector gets sized properly and will be overridden
773void MaterialManager::getMediumIDMappingAsVector(const char* modname, std::vector<int>& mapping) const
774{
775 mapping.clear();
776
777 auto lookupiter = mMediumMap.find(modname);
778 if (lookupiter == mMediumMap.end()) {
779 return;
780 }
781 auto lookup = lookupiter->second;
782
783 // get the biggest mapped value (maps are sorted in keys)
784 auto maxkey = lookup.rbegin()->first;
785 // resize mapping and initialize with -1 by default
786 mapping.resize(maxkey + 1, -1);
787 // fill vector with entries from map
788 for (auto& p : lookup) {
789 mapping[p.first] = p.second;
790 }
791}
792
793void MaterialManager::getMediaWithSpecialProcess(EProc process, std::vector<int>& mediumProcessVector) const
794{
795 // clear
796 mediumProcessVector.clear();
797 // resize to maximum number of global IDs for which special processes are set. In case process is not
798 // implemented for a certain medium, value is -1
799 mediumProcessVector.resize(mMediumProcessMap.size(), -1);
800 // find media
801 for (auto& m : mMediumProcessMap) {
802 // loop over processes in medium
803 for (auto& p : m.second) {
804 // push medium ID if process is there
805 if (p.first == process) {
806 mediumProcessVector[m.first] = p.second;
807 break;
808 }
809 }
810 }
811}
812
813void MaterialManager::getMediaWithSpecialCut(ECut cut, std::vector<Float_t>& mediumCutVector) const
814{
815 // clear
816 mediumCutVector.clear();
817 // resize to maximum number of global IDs for which special cuts are set. In case cut is not implemented
818 // for a certain medium, value is -1.
819 mediumCutVector.resize(mMediumCutMap.size(), -1.);
820 // find media
821 for (auto& m : mMediumCutMap) {
822 // loop over cuts in medium
823 for (auto& c : m.second) {
824 // push medium ID if cut is there
825 if (c.first == cut) {
826 mediumCutVector[m.first] = c.second;
827 break;
828 }
829 }
830 }
831}
832
834void MaterialManager::getDefaultProcesses(std::vector<std::pair<EProc, int>>& processVector)
835{
836 processVector.clear();
837 for (auto& m : mDefaultProcessMap) {
838 processVector.emplace_back(m.first, m.second);
839 }
840}
842void MaterialManager::getDefaultCuts(std::vector<std::pair<ECut, Float_t>>& cutVector)
843{
844 cutVector.clear();
845 for (auto& m : mDefaultCutMap) {
846 cutVector.emplace_back(m.first, m.second);
847 }
848}
850void MaterialManager::getSpecialProcesses(int globalindex, std::vector<std::pair<EProc, int>>& processVector)
851{
852 processVector.clear();
853 if (mMediumProcessMap.find(globalindex) != mMediumProcessMap.end()) {
854 for (auto& m : mMediumProcessMap[globalindex]) {
855 processVector.emplace_back(m.first, m.second);
856 }
857 }
858}
860void MaterialManager::getSpecialProcesses(const char* modname, int localindex, std::vector<std::pair<EProc, int>>& processVector)
861{
862 int globalindex = getMediumID(modname, localindex);
863 if (globalindex != -1) {
864 getSpecialProcesses(globalindex, processVector);
865 }
866}
868void MaterialManager::getSpecialCuts(int globalindex, std::vector<std::pair<ECut, Float_t>>& cutVector)
869{
870 cutVector.clear();
871 if (mMediumCutMap.find(globalindex) != mMediumCutMap.end()) {
872 for (auto& m : mMediumCutMap[globalindex]) {
873 cutVector.emplace_back(m.first, m.second);
874 }
875 }
876}
878void MaterialManager::getSpecialCuts(const char* modname, int localindex, std::vector<std::pair<ECut, Float_t>>& cutVector)
879{
880 int globalindex = getMediumID(modname, localindex);
881 if (globalindex != -1) {
882 getSpecialCuts(globalindex, cutVector);
883 }
884}
885
886const char* MaterialManager::getModuleFromMediumID(int globalindex) const
887{
888 // loop over module names and corresponding local<->global mapping
889 for (auto& m : mMediumMap) {
890 for (auto& i : m.second) {
891 // is the global index there?
892 if (i.second == globalindex) {
893 // \note maybe unsafe in case mMediumMap is altered in the same scope where the returned C string is used
894 // since that points to memory of string it was derived from.
895 return m.first.c_str();
896 }
897 }
898 }
899 // module is UNKNOWN if global medium ID could not be found.
900 return "UNKNOWN";
901}
902
904const char* MaterialManager::getMediumNameFromMediumID(int globalindex) const
905{
906 // Get the name of the medium.
907 // TODO: avoid linear search
908 for (auto& n : mMediumNameToGlobalIndexMap) {
909 if (n.second == globalindex) {
910 // \note maybe unsafe in case mMediumMap is altered in the same scope where the returned C string is used since
911 // that points to memory of string it was derived from.
912 return n.first.c_str();
913 }
914 }
915 return "UNKNOWN";
916}
917
920void MaterialManager::printContainingMedia(std::string const& volumename)
921{
922 auto vol = gGeoManager->FindVolumeFast(volumename.c_str());
923 if (vol == nullptr) {
924 LOG(warn) << "No volume found; Cannot query medias";
925 }
926 std::set<TGeoMedium const*> media;
927
928 // a little well encapsulated helper to code the recursive visitor
929 // pure lambda cannot be recursive
930 std::function<void(TGeoVolume const* vol, std::set<TGeoMedium const*>& mediumset)> recursivevisitor;
931 recursivevisitor = [&recursivevisitor](TGeoVolume const* vol, std::set<TGeoMedium const*>& mediumset) {
932 // exclude assemblies
933 if (!vol->IsAssembly()) {
934 mediumset.insert(vol->GetMedium());
935 }
936 const int n = vol->GetNdaughters();
937 for (int i = 0; i < n; ++i) {
938 auto daughter = vol->GetNode(i)->GetVolume();
939 recursivevisitor(daughter, mediumset);
940 }
941 };
942 recursivevisitor(vol, media);
943
944 // simply print the media
945 for (auto m : media) {
946 std::cout << m->GetName() << "\n";
947 }
948}
949
uint8_t lookup(const char input) noexcept
int32_t i
ClassImp(IdPath)
uint16_t pos
Definition RawData.h:3
uint32_t c
Definition RawData.h:2
std::ostringstream debug
void printCuts(std::ostream &stream) const
Print all cuts for all media as well as defaults.
void getSpecialCuts(int globalindex, std::vector< std::pair< ECut, Float_t > > &cutVector)
Get special cuts for global medium ID.
void SpecialCuts(const char *modname, int localindex, const std::initializer_list< std::pair< ECut, Float_t > > &parIDValMap)
void getMediumIDMappingAsVector(const char *modname, std::vector< int > &mapping) const
void getDefaultProcesses(std::vector< std::pair< EProc, int > > &processVector)
Fill vector with default processes.
void printProcesses(std::ostream &stream) const
Print all processes for all media as well as defaults.
void enableSpecialCuts(bool val=true)
void writeCutsAndProcessesToJSON(std::string const &filename="")
void loadCutsAndProcessesFromFile(const char *modname, const char *filename)
load cuts and process flags from a data file (like AliRoot did)
void loadCutsAndProcessesFromJSON(ESpecial special=ESpecial::kFALSE, std::string const &filename="")
void getMediaWithSpecialProcess(EProc process, std::vector< int > &mediumProcessVector) const
get global medium IDs where special process is set along with process value
int getMediumID(const char *modname, int imed) const
void getSpecialProcesses(int globalindex, std::vector< std::pair< EProc, int > > &processVector)
Get special processes for global medium ID.
const char * getMediumNameFromMediumID(int globalindex) const
Get medium name from global medium ID.
void getMediaWithSpecialCut(ECut cut, std::vector< Float_t > &mediumCutVector) const
get global medium IDs where special cut is set along with cut value
static void printContainingMedia(std::string const &volumename)
void Mixture(const char *modname, Int_t imat, const char *name, Float_t *a, Float_t *z, Float_t dens, Int_t nlmat, Float_t *wmat)
void SpecialProcess(const char *modname, int localindex, EProc parID, int val)
Custom setting of process or cut given parameter name and value.
void Medium(const char *modname, 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)
void SpecialCut(const char *modname, int localindex, ECut parID, Float_t val)
Custom setting of process or cut given parameter name and value.
void Material(const char *modname, 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)
void enableSpecialProcesses(bool val=true)
Set flags whether to use special cuts and process settings.
const char * getModuleFromMediumID(int globalindex) const
Get module name which has medium of a certain global medium ID.
ESpecial
scoped enum to decide whether settings are done globally or for a certain medium
int getMaterialID(const char *modname, int imat) const
void getDefaultCuts(std::vector< std::pair< ECut, Float_t > > &cutVector)
Fill vector with default cuts.
TGeoMedium * getTGeoMedium(const std::string &modname, int localid)
GLdouble n
Definition glcorearb.h:1982
const GLfloat * m
Definition glcorearb.h:4066
GLuint index
Definition glcorearb.h:781
GLuint const GLchar * name
Definition glcorearb.h:781
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat * val
Definition glcorearb.h:1582
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLuint GLuint stream
Definition glcorearb.h:1806
GLenum GLuint GLenum GLsizei const GLchar * buf
Definition glcorearb.h:2514
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
ECut
cuts available
EProc
processes available
std::string filename()
Definition list.h:40
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"