-
Notifications
You must be signed in to change notification settings - Fork 496
Expand file tree
/
Copy patho2sim.C
More file actions
334 lines (299 loc) · 12.4 KB
/
o2sim.C
File metadata and controls
334 lines (299 loc) · 12.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
#include "build_geometry.C"
#if !defined(__CLING__) || defined(__ROOTCLING__)
#include <Generators/PrimaryGenerator.h>
#include <Generators/GeneratorFactory.h>
#include <Generators/Generator.h>
#include "SimulationDataFormat/O2DatabasePDG.h"
#include "SimulationDataFormat/MCEventHeader.h"
#include <SimConfig/SimConfig.h>
#include <SimConfig/SimParams.h>
#include <CommonUtils/ConfigurableParam.h>
#include <CommonUtils/RngHelper.h>
#include <TStopwatch.h>
#include <memory>
#include "DataFormatsParameters/GRPObject.h"
#include "DataFormatsParameters/GRPECSObject.h"
#include "DataFormatsParameters/GRPMagField.h"
#include "DataFormatsParameters/GRPLHCIFData.h"
#include "FairParRootFileIo.h"
#include "FairSystemInfo.h"
#include <SimSetup/SimSetup.h>
#include <Steer/O2RunSim.h>
#include <DetectorsBase/MaterialManager.h>
#include <CCDB/BasicCCDBManager.h>
#include <CommonUtils/NameConf.h>
#include "DetectorsBase/Aligner.h"
#include <FairRootFileSink.h>
#include <FairField.h>
#include <unistd.h>
#include <sstream>
#endif
#include "migrateSimFiles.C"
#include <boost/property_tree/ptree.hpp>
void check_notransport()
{
// Sometimes we just want to inspect
// the generator kinematics and prohibit any transport.
// We allow users to give the noGeant option. In this case
// we set the geometry cuts to almost zero. Other adjustments (disable physics procs) could
// be done on top.
// This is merely offered for user convenience as it can be done from outside as well.
auto& confref = o2::conf::SimConfig::Instance();
if (confref.isNoGeant()) {
LOG(info) << "Initializing without Geant transport by applying very tight geometry cuts";
o2::conf::ConfigurableParam::setValue("SimCutParams", "maxRTracking", 0.0000001); // 1 nanometer of tracking
o2::conf::ConfigurableParam::setValue("SimCutParams", "maxAbsZTracking", 0.0000001); // 1 nanometer of tracking
}
}
FairRunSim* o2sim_init(bool asservice, bool evalmat = false)
{
auto& confref = o2::conf::SimConfig::Instance();
// set the global information about the number of events to be generated
unsigned int nTotalEvents = confref.getNEvents();
o2::eventgen::Generator::setTotalNEvents(nTotalEvents);
// initialize CCDB service
auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
// fix the timestamp early
uint64_t timestamp = confref.getTimestamp();
uint64_t runEnd = timestamp + 3600000;
// see if we have a run number but not a timestamp
auto run_number = confref.getRunNumber();
if (run_number != -1) {
if (confref.getConfigData().mTimestampMode == o2::conf::TimeStampMode::kNow) {
// fix the time by talking to CCDB
auto [sor, eor] = ccdbmgr.getRunDuration(run_number);
LOG(info) << "Have run number. Fixing timestamp to " << sor;
timestamp = sor;
runEnd = eor;
}
}
ccdbmgr.setTimestamp(timestamp);
ccdbmgr.setURL(confref.getConfigData().mCCDBUrl);
// try to verify connection
if (!ccdbmgr.isHostReachable()) {
LOG(error) << "Could not setup CCDB connection";
} else {
LOG(info) << "Initialized CCDB Manager at URL: " << ccdbmgr.getURL();
LOG(info) << "Initialized CCDB Manager with timestamp : " << ccdbmgr.getTimestamp();
}
check_notransport();
// update the parameters from an INI/JSON file, if given (overrides code-based version)
o2::conf::ConfigurableParam::updateFromFile(confref.getConfigFile());
// update the parameters from stuff given at command line (overrides file-based version)
o2::conf::ConfigurableParam::updateFromString(confref.getKeyValueString());
// write the final configuration file
o2::conf::ConfigurableParam::writeINI(o2::base::NameConf::getMCConfigFileName(confref.getOutPrefix()));
// set seed
auto seed = o2::utils::RngHelper::setGRandomSeed(confref.getStartSeed());
LOG(info) << "RNG INITIAL SEED " << seed;
auto genconfig = confref.getGenerator();
FairRunSim* run = new o2::steer::O2RunSim(asservice, evalmat);
run->SetImportTGeoToVMC(false); // do not import TGeo to VMC since the latter is built together with TGeo
run->SetSimSetup([confref]() { o2::SimSetup::setup(confref.getMCEngine().c_str()); });
run->SetRunId(timestamp);
auto pid = getpid();
std::stringstream s;
s << confref.getOutPrefix();
if (asservice) {
s << "_" << pid;
}
s << ".root";
std::string outputfilename = s.str();
run->SetSink(new FairRootFileSink(outputfilename.c_str())); // Output file
run->SetName(confref.getMCEngine().c_str()); // Transport engine
run->SetIsMT(confref.getIsMT()); // MT mode
/** set event header **/
auto header = new o2::dataformats::MCEventHeader();
run->SetMCEventHeader(header);
// construct geometry / including magnetic field
auto flg = TGeoManager::LockDefaultUnits(false);
TGeoManager::SetDefaultUnits(TGeoManager::kRootUnits);
TGeoManager::LockDefaultUnits(flg);
build_geometry(run);
// setup generator
auto embedinto_filename = confref.getEmbedIntoFileName();
auto primGen = new o2::eventgen::PrimaryGenerator();
if (!embedinto_filename.empty()) {
primGen->embedInto(embedinto_filename);
}
if (!asservice) {
o2::eventgen::GeneratorFactory::setPrimaryGenerator(confref, primGen);
}
run->SetGenerator(primGen);
// Timer
TStopwatch timer;
timer.Start();
o2::detectors::DetID::mask_t detMask{};
o2::detectors::DetID::mask_t readoutDetMask{};
{
auto& modulelist = o2::conf::SimConfig::Instance().getActiveModules();
for (const auto& md : modulelist) {
int id = o2::detectors::DetID::nameToID(md.c_str());
if (id >= o2::detectors::DetID::First) {
detMask |= o2::detectors::DetID::getMask(id);
if (isReadout(md)) {
readoutDetMask |= o2::detectors::DetID::getMask(id);
}
}
}
if (readoutDetMask.none()) {
LOG(info) << "Hit creation disabled for all detectors";
}
// somewhat ugly, but this is the most straighforward way to make sure the detectors to align
// don't include detectors which are not activated
auto& aligner = o2::base::Aligner::Instance();
auto detMaskAlign = aligner.getDetectorsMask() & detMask;
aligner.setValue(fmt::format("{}.mDetectors", aligner.getName()), o2::detectors::DetID::getNames(detMaskAlign, ','));
}
// run init
run->Init();
// add ALICE particles to TDatabasePDG singleton
o2::O2DatabasePDG::addALICEParticles(TDatabasePDG::Instance());
long runStart = timestamp;
{
// store GRPobject
o2::parameters::GRPObject grp;
if (run_number != -1) {
grp.setRun(run_number);
} else {
grp.setRun(run->GetRunId());
}
uint64_t runStart = timestamp;
grp.setTimeStart(runStart);
grp.setTimeEnd(runEnd);
grp.setDetsReadOut(readoutDetMask);
// CTP is not a physical detector, just flag in the GRP if requested
if (isReadout("CTP")) {
grp.addDetReadOut(o2::detectors::DetID::CTP);
}
grp.print();
printf("VMC: %p\n", TVirtualMC::GetMC());
auto field = dynamic_cast<o2::field::MagneticField*>(run->GetField());
if (field) {
o2::units::Current_t currDip = field->getCurrentDipole();
o2::units::Current_t currL3 = field->getCurrentSolenoid();
grp.setL3Current(currL3);
grp.setDipoleCurrent(currDip);
grp.setFieldUniformity(field->IsUniform());
}
// save
std::string grpfilename = o2::base::NameConf::getGRPFileName(confref.getOutPrefix());
TFile grpF(grpfilename.c_str(), "recreate");
grpF.WriteObjectAny(&grp, grp.Class(), o2::base::NameConf::CCDBOBJECT.data());
}
// create GRPECS object
{
o2::parameters::GRPECSObject grp;
grp.setRun(run->GetRunId());
grp.setTimeStart(runStart);
grp.setTimeEnd(runEnd);
grp.setNHBFPerTF(128); // might be overridden later
grp.setDetsReadOut(readoutDetMask);
if (isReadout("CTP")) {
grp.addDetReadOut(o2::detectors::DetID::CTP);
}
grp.setIsMC(true);
grp.setRunType(o2::parameters::GRPECSObject::RunType::PHYSICS);
// grp.setDataPeriod("mc"); // decide what to put here
std::string grpfilename = o2::base::NameConf::getGRPECSFileName(confref.getOutPrefix());
TFile grpF(grpfilename.c_str(), "recreate");
grpF.WriteObjectAny(&grp, grp.Class(), o2::base::NameConf::CCDBOBJECT.data());
}
// create GRPMagField object
{
o2::parameters::GRPMagField grp;
auto field = dynamic_cast<o2::field::MagneticField*>(run->GetField());
if (!field) {
// this is not the ordinary Run3 MagneticField
// Let's see if it is another FairField implementation.
LOG(warn) << "No o2::field::MagneticField instance available; Not writing GRP - beware that propagation to other tasks may not work. Checking if it is at least a FairFied...";
if (!dynamic_cast<FairField*>(run->GetField())) {
LOGP(fatal, "Failed to get magnetic field from the FairRunSim");
} else {
LOG(warn) << " ... FairField found";
}
} else {
o2::units::Current_t currDip = field->getCurrentDipole();
o2::units::Current_t currL3 = field->getCurrentSolenoid();
grp.setL3Current(currL3);
grp.setDipoleCurrent(currDip);
grp.setFieldUniformity(field->IsUniform());
std::string grpfilename = o2::base::NameConf::getGRPMagFieldFileName(confref.getOutPrefix());
TFile grpF(grpfilename.c_str(), "recreate");
grpF.WriteObjectAny(&grp, grp.Class(), o2::base::NameConf::CCDBOBJECT.data());
}
}
// create GRPLHCIF object (just a placeholder, bunch filling will be set in digitization)
{
o2::parameters::GRPLHCIFData grp;
// eventually we need to set the beam info from the generator, at the moment put some plausible values
grp.setFillNumberWithTime(runStart, 0); // RS FIXME
grp.setInjectionSchemeWithTime(runStart, ""); // RS FIXME
grp.setBeamEnergyPerZWithTime(runStart, 6.8e3); // RS FIXME
grp.setAtomicNumberB1WithTime(runStart, 1.); // RS FIXME
grp.setAtomicNumberB2WithTime(runStart, 1.); // RS FIXME
grp.setCrossingAngleWithTime(runStart, 0.); // RS FIXME
grp.setBeamAZ();
std::string grpfilename = o2::base::NameConf::getGRPLHCIFFileName(confref.getOutPrefix());
TFile grpF(grpfilename.c_str(), "recreate");
grpF.WriteObjectAny(&grp, grp.Class(), o2::base::NameConf::CCDBOBJECT.data());
}
// print summary about cuts and processes used
auto& matmgr = o2::base::MaterialManager::Instance();
std::ofstream cutfile(o2::base::NameConf::getCutProcFileName(confref.getOutPrefix()));
matmgr.printCuts(cutfile);
matmgr.printProcesses(cutfile);
timer.Stop();
Double_t rtime = timer.RealTime();
Double_t ctime = timer.CpuTime();
// extract max memory usage for init
FairSystemInfo sysinfo;
LOG(info) << "Init: Real time " << rtime << " s, CPU time " << ctime << "s";
LOG(info) << "Init: Memory used " << sysinfo.GetMaxMemory() << " MB";
return run;
}
// only called from the normal o2sim
void o2sim_run(FairRunSim* run, bool asservice)
{
TStopwatch timer;
timer.Start();
auto& confref = o2::conf::SimConfig::Instance();
if (!asservice) {
run->Run(confref.getNEvents());
} else {
run->Run(1);
}
// Finish
timer.Stop();
Double_t rtime = timer.RealTime();
Double_t ctime = timer.CpuTime();
// extract max memory usage
FairSystemInfo sysinfo;
LOG(info) << "Macro finished succesfully.";
LOG(info) << "Real time " << rtime << " s, CPU time " << ctime << "s";
LOG(info) << "Memory used " << sysinfo.GetMaxMemory() << " MB";
// migrate to file format where hits sit in separate files
// (Note: The parallel version is doing this intrinsically;
// The serial version uses FairRootManager IO which handles a common file IO for all outputs)
if (!asservice) {
LOG(info) << "Migrating simulation output to separate hit file format";
migrateSimFiles(confref.getOutPrefix().c_str());
}
}
// asservice: in a parallel device-based context?
void o2sim(bool asservice = false, bool evalmat = false)
{
auto run = o2sim_init(asservice, evalmat);
o2sim_run(run, asservice);
delete run;
}