Project
Loading...
Searching...
No Matches
DigitizationContext.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
16#include <TChain.h>
17#include <TFile.h>
18#include <iostream>
19#include <numeric> // for iota
20#include <MathUtils/Cartesian.h>
22#include <filesystem>
23
24using namespace o2::steer;
25
26void DigitizationContext::printCollisionSummary(bool withQED, int truncateOutputTo) const
27{
28 std::cout << "Summary of DigitizationContext --\n";
29 std::cout << "Maximal parts per collision " << mMaxPartNumber << "\n";
30 std::cout << "Collision parts taken from simulations specified by prefix:\n";
31 for (int p = 0; p < mSimPrefixes.size(); ++p) {
32 std::cout << "Part " << p << " : " << mSimPrefixes[p] << "\n";
33 }
34 std::cout << "QED information included " << isQEDProvided() << "\n";
35 if (withQED) {
36 std::cout << "Number of Collisions " << mEventRecords.size() << "\n";
37 std::cout << "Number of QED events " << mEventRecordsWithQED.size() - mEventRecords.size() << "\n";
38 // loop over combined stuff
39 for (int i = 0; i < mEventRecordsWithQED.size(); ++i) {
40 if (truncateOutputTo >= 0 && i > truncateOutputTo) {
41 std::cout << "--- Output truncated to " << truncateOutputTo << " ---\n";
42 break;
43 }
44 std::cout << "Record " << i << " TIME " << mEventRecordsWithQED[i];
45 for (auto& e : mEventPartsWithQED[i]) {
46 std::cout << " (" << e.sourceID << " , " << e.entryID << ")";
47 }
48 std::cout << "\n";
49 }
50 } else {
51 std::cout << "Number of Collisions " << mEventRecords.size() << "\n";
52 if (mEventPartsWithQED.size() > 0) {
53 auto num_qed_events = mEventPartsWithQED.size() - mEventRecords.size();
54 if (num_qed_events > 0) {
55 std::cout << "Number of QED events (but not shown) " << num_qed_events << "\n";
56 // find first and last QED collision so that we can give the range in orbits where these
57 // things are included
58 auto firstQEDcoll_iter = std::find_if(mEventPartsWithQED.begin(), mEventPartsWithQED.end(),
59 [](const std::vector<EventPart>& vec) {
60 return std::find_if(vec.begin(), vec.end(), [](EventPart const& p) { return p.sourceID == 99; }) != vec.end();
61 });
62
63 auto lastColl_iter = std::find_if(mEventPartsWithQED.rbegin(), mEventPartsWithQED.rend(),
64 [](const std::vector<EventPart>& vec) {
65 return std::find_if(vec.begin(), vec.end(), [](EventPart const& p) { return p.sourceID == 99; }) != vec.end();
66 });
67
68 auto firstindex = std::distance(mEventPartsWithQED.begin(), firstQEDcoll_iter);
69 auto lastindex = std::distance(mEventPartsWithQED.begin(), lastColl_iter.base()) - 1;
70 std::cout << "QED from: " << mEventRecordsWithQED[firstindex] << " ---> " << mEventRecordsWithQED[lastindex] << "\n";
71 }
72 }
73 for (int i = 0; i < mEventRecords.size(); ++i) {
74 if (truncateOutputTo >= 0 && i > truncateOutputTo) {
75 std::cout << "--- Output truncated to " << truncateOutputTo << " ---\n";
76 break;
77 }
78 std::cout << "Collision " << i << " TIME " << mEventRecords[i];
79 for (auto& e : mEventParts[i]) {
80 std::cout << " (" << e.sourceID << " , " << e.entryID << ")";
81 }
82 if (i < mInteractionVertices.size()) {
83 std::cout << " sampled vertex : " << mInteractionVertices[i];
84 }
85 std::cout << "\n";
86 }
87 }
88}
89
90void DigitizationContext::setSimPrefixes(std::vector<std::string> const& prefixes)
91{
92 mSimPrefixes = prefixes;
93}
94
95bool DigitizationContext::initSimChains(o2::detectors::DetID detid, std::vector<TChain*>& simchains) const
96{
97 if (!(simchains.size() == 0)) {
98 // nothing to do ... already setup
99 return false;
100 }
101
102 // check that all files are present, otherwise quit
103 for (int source = 0; source < mSimPrefixes.size(); ++source) {
104 if (!std::filesystem::exists(o2::base::DetectorNameConf::getHitsFileName(detid, mSimPrefixes[source].data()))) {
105 LOG(info) << "Not hit file present for " << detid.getName() << " (exiting SimChain setup)";
106 return false;
107 }
108 }
109
110 simchains.emplace_back(new TChain("o2sim"));
111 // add the main (background) file
112 simchains.back()->AddFile(o2::base::DetectorNameConf::getHitsFileName(detid, mSimPrefixes[0].data()).c_str());
113
114 for (int source = 1; source < mSimPrefixes.size(); ++source) {
115 simchains.emplace_back(new TChain("o2sim"));
116 // add signal files
117 simchains.back()->AddFile(o2::base::DetectorNameConf::getHitsFileName(detid, mSimPrefixes[source].data()).c_str());
118 }
119
120 // QED part
121 if (mEventRecordsWithQED.size() > 0) {
122 if (mSimPrefixes.size() >= QEDSOURCEID) {
123 LOG(fatal) << "Too many signal chains; crashes with QED source ID";
124 }
125
126 // it might be better to use an unordered_map for the simchains but this requires interface changes
127 simchains.resize(QEDSOURCEID + 1, nullptr);
128 simchains[QEDSOURCEID] = new TChain("o2sim");
129 simchains[QEDSOURCEID]->AddFile(o2::base::DetectorNameConf::getHitsFileName(detid, mQEDSimPrefix).c_str());
130 }
131
132 return true;
133}
134
138bool DigitizationContext::initSimKinematicsChains(std::vector<TChain*>& simkinematicschains) const
139{
140 if (!(simkinematicschains.size() == 0)) {
141 // nothing to do ... already setup
142 return false;
143 }
144
145 simkinematicschains.emplace_back(new TChain("o2sim"));
146 // add the main (background) file
147 simkinematicschains.back()->AddFile(o2::base::NameConf::getMCKinematicsFileName(mSimPrefixes[0].data()).c_str());
148
149 for (int source = 1; source < mSimPrefixes.size(); ++source) {
150 simkinematicschains.emplace_back(new TChain("o2sim"));
151 // add signal files
152 simkinematicschains.back()->AddFile(o2::base::NameConf::getMCKinematicsFileName(mSimPrefixes[source].data()).c_str());
153 }
154
155 // we add QED, if used in the digitization context
156 if (mEventRecordsWithQED.size() > 0) {
157 if (mSimPrefixes.size() >= QEDSOURCEID) {
158 LOG(fatal) << "Too many signal chains; crashes with QED source ID";
159 }
160
161 // it might be better to use an unordered_map for the simchains but this requires interface changes
162 simkinematicschains.resize(QEDSOURCEID + 1, nullptr);
163 simkinematicschains[QEDSOURCEID] = new TChain("o2sim");
164 simkinematicschains[QEDSOURCEID]->AddFile(o2::base::DetectorNameConf::getMCKinematicsFileName(mQEDSimPrefix).c_str());
165 }
166
167 return true;
168}
169
171{
172 if (mMaxPartNumber == 1) {
173 return true;
174 }
175
176 auto checkVertexPair = [](math_utils::Point3D<double> const& p1, math_utils::Point3D<double> const& p2) -> bool {
177 return (p2 - p1).Mag2() < 1E-6;
178 };
179
180 std::vector<TChain*> kinematicschain;
181 std::vector<TBranch*> headerbranches;
182 std::vector<o2::dataformats::MCEventHeader*> headers;
183 std::vector<math_utils::Point3D<double>> vertices;
184 initSimKinematicsChains(kinematicschain);
185 bool consistent = true;
186 if (kinematicschain.size() > 0) {
187 headerbranches.resize(kinematicschain.size(), nullptr);
188 headers.resize(kinematicschain.size(), nullptr);
189 // loop over all collisions in this context
190 int collisionID = 0;
191 for (auto& collision : getEventParts()) {
192 collisionID++;
193 vertices.clear();
194 for (auto& part : collision) {
195 const auto source = part.sourceID;
196 const auto entry = part.entryID;
197 auto chain = kinematicschain[source];
198 if (!headerbranches[source]) {
199 headerbranches[source] = chain->GetBranch("MCEventHeader.");
200 headerbranches[source]->SetAddress(&headers[source]);
201 }
202 // get the MCEventHeader to read out the vertex
203 headerbranches[source]->GetEntry(entry);
204 auto header = headers[source];
205 vertices.emplace_back(header->GetX(), header->GetY(), header->GetZ());
206 }
207 // analyse vertex matching
208 bool thiscollision = true;
209 const auto& p1 = vertices[0];
210 for (int j = 1; j < vertices.size(); ++j) {
211 const auto& p2 = vertices[j];
212 bool thischeck = checkVertexPair(p1, p2);
213 thiscollision &= thischeck;
214 }
215 if (verbose && !thiscollision) {
216 std::stringstream text;
217 text << "Found inconsistent vertices for digit collision ";
218 text << collisionID << " : ";
219 for (auto& p : vertices) {
220 text << p << " ";
221 }
222 LOG(error) << text.str();
223 }
224 consistent &= thiscollision;
225 }
226 }
227 return consistent;
228}
229
231{
232 if (!mGRP) {
233 // we take the GRP from the background file
234 // maybe we should add a check that all GRPs are consistent ..
235 mGRP = o2::parameters::GRPObject::loadFrom(mSimPrefixes[0]);
236 }
237 return *mGRP;
238}
239
240void DigitizationContext::saveToFile(std::string_view filename) const
241{
242 // checks if the path content of filename exists ... otherwise it is created before creating the ROOT file
243 auto ensure_path_exists = [](std::string_view filename) {
244 try {
245 // Extract the directory path from the filename
246 std::filesystem::path file_path(filename);
247 std::filesystem::path dir_path = file_path.parent_path();
248
249 // Check if the directory path is empty (which means filename was just a name without path)
250 if (dir_path.empty()) {
251 // nothing to do
252 return true;
253 }
254
255 // Create directories if they do not exist
256 if (!std::filesystem::exists(dir_path)) {
257 if (std::filesystem::create_directories(dir_path)) {
258 // std::cout << "Directories created successfully: " << dir_path.string() << std::endl;
259 return true;
260 } else {
261 std::cerr << "Failed to create directories: " << dir_path.string() << std::endl;
262 return false;
263 }
264 }
265 return true;
266 } catch (const std::filesystem::filesystem_error& ex) {
267 std::cerr << "Filesystem error: " << ex.what() << std::endl;
268 return false;
269 } catch (const std::exception& ex) {
270 std::cerr << "General error: " << ex.what() << std::endl;
271 return false;
272 }
273 };
274
275 if (!ensure_path_exists(filename)) {
276 LOG(error) << "Filename contains path component which could not be created";
277 return;
278 }
279
280 TFile file(filename.data(), "RECREATE");
281 if (file.IsOpen()) {
282 auto cl = TClass::GetClass(typeid(*this));
283 file.WriteObjectAny(this, cl, "DigitizationContext");
284 file.Close();
285 } else {
286 LOG(error) << "Could not write to file " << filename.data();
287 }
288}
289
291{
292 std::string tmpFile;
293 if (filename == "") {
295 filename = tmpFile;
296 }
297 DigitizationContext* incontext = nullptr;
298 TFile file(filename.data(), "OPEN");
299 file.GetObject("DigitizationContext", incontext);
300 return incontext;
301}
302
303void DigitizationContext::fillQED(std::string_view QEDprefix, int max_events, double qedrate)
304{
305 if (mEventRecords.size() <= 1) {
306 // nothing to do
307 return;
308 }
309
310 o2::steer::InteractionSampler qedInteractionSampler;
311 qedInteractionSampler.setBunchFilling(mBCFilling);
312
313 // get first and last "hadronic" interaction records and let
314 // QED events range from the first bunch crossing to the last bunch crossing
315 // in this range
316 auto first = mEventRecords.front();
317 auto last = mEventRecords.back();
318 first.bc = 0;
320 LOG(info) << "QED RATE " << qedrate;
321 qedInteractionSampler.setInteractionRate(qedrate);
322 qedInteractionSampler.setFirstIR(first);
323 qedInteractionSampler.init();
324 qedInteractionSampler.print();
325 std::vector<o2::InteractionTimeRecord> qedinteractionrecords;
327 LOG(info) << "GENERATING COL TIMES";
328 t = qedInteractionSampler.generateCollisionTime();
329 while ((t = qedInteractionSampler.generateCollisionTime()) < last) {
330 qedinteractionrecords.push_back(t);
331 }
332 LOG(info) << "DONE GENERATING COL TIMES";
333 fillQED(QEDprefix, qedinteractionrecords, max_events, false);
334}
335
336void DigitizationContext::fillQED(std::string_view QEDprefix, std::vector<o2::InteractionTimeRecord> const& irecord, int max_events, bool fromKinematics)
337{
338 mQEDSimPrefix = QEDprefix;
339
340 std::vector<std::vector<o2::steer::EventPart>> qedEventParts;
341
342 int numberQEDevents = max_events; // if this is -1 there will be no limitation
343 if (fromKinematics) {
344 // we need to fill the QED parts (using a simple round robin scheme)
345 auto qedKinematicsName = o2::base::NameConf::getMCKinematicsFileName(mQEDSimPrefix);
346 // find out how many events are stored
347 TFile f(qedKinematicsName.c_str(), "OPEN");
348 auto t = (TTree*)f.Get("o2sim");
349 if (!t) {
350 LOG(error) << "No QED kinematics found";
351 throw std::runtime_error("No QED kinematics found");
352 }
353 numberQEDevents = t->GetEntries();
354 }
355
356 int eventID = 0;
357 for (auto& tmp : irecord) {
358 std::vector<EventPart> qedpart;
359 qedpart.emplace_back(QEDSOURCEID, eventID++);
360 qedEventParts.push_back(qedpart);
361 if (eventID == numberQEDevents) {
362 eventID = 0;
363 }
364 }
365
366 // we need to do the interleaved event records for detectors consuming both
367 // normal and QED events
368 // --> merge both; sort first according to times and sort second one according to same order
369 auto combinedrecords = mEventRecords;
370 combinedrecords.insert(combinedrecords.end(), irecord.begin(), irecord.end());
371 auto combinedparts = mEventParts;
372 combinedparts.insert(combinedparts.end(), qedEventParts.begin(), qedEventParts.end());
373
374 // get sorted index vector based on event records
375 std::vector<size_t> idx(combinedrecords.size());
376 std::iota(idx.begin(), idx.end(), 0);
377
378 std::stable_sort(idx.begin(), idx.end(),
379 [&combinedrecords](size_t i1, size_t i2) { return combinedrecords[i1] < combinedrecords[i2]; });
380
381 mEventRecordsWithQED.clear();
382 mEventPartsWithQED.clear();
383 for (int i = 0; i < idx.size(); ++i) {
384 mEventRecordsWithQED.push_back(combinedrecords[idx[i]]);
385 mEventPartsWithQED.push_back(combinedparts[idx[i]]);
386 }
387}
388
389namespace
390{
391// a common helper for timeframe structure
392std::vector<std::pair<int, int>> getTimeFrameBoundaries(std::vector<o2::InteractionTimeRecord> const& irecords, long startOrbit, long orbitsPerTF)
393{
394 std::vector<std::pair<int, int>> result;
395
396 // the goal is to determine timeframe boundaries inside the interaction record vectors
397 // determine if we can do anything
398 if (irecords.size() == 0) {
399 // nothing to do
400 return result;
401 }
402
403 if (irecords.back().orbit < startOrbit) {
404 LOG(error) << "start orbit larger than last collision entry";
405 return result;
406 }
407
408 // skip to the first index falling within our constrained
409 int left = 0;
410 while (left < irecords.size() && irecords[left].orbit < startOrbit) {
411 left++;
412 }
413
414 // now we can start (2 pointer approach)
415 auto right = left;
416 int timeframe_count = 1;
417 while (right < irecords.size()) {
418 if (irecords[right].orbit >= startOrbit + timeframe_count * orbitsPerTF) {
419 // we finished one timeframe
420 result.emplace_back(std::pair<int, int>(left, right - 1));
421 timeframe_count++;
422 left = right;
423 }
424 right++;
425 }
426 // finished last timeframe
427 result.emplace_back(std::pair<int, int>(left, right - 1));
428 return result;
429}
430
431// a common helper for timeframe structure - includes indices for orbits-early (orbits from last timeframe still affecting current one)
432std::vector<std::tuple<int, int, int>> getTimeFrameBoundaries(std::vector<o2::InteractionTimeRecord> const& irecords,
433 long startOrbit,
434 long orbitsPerTF,
435 float orbitsEarly)
436{
437 // we could actually use the other method first ... then do another pass to fix the early-index ... or impact index
438 auto true_indices = getTimeFrameBoundaries(irecords, startOrbit, orbitsPerTF);
439
440 std::vector<std::tuple<int, int, int>> indices_with_early{};
441 for (int ti = 0; ti < true_indices.size(); ++ti) {
442 // for each timeframe we copy the true indices
443 auto& tf_range = true_indices[ti];
444
445 // init new index without fixing the early index yet
446 indices_with_early.push_back(std::make_tuple(tf_range.first, tf_range.second, -1));
447
448 // from the second timeframe on we can determine the index in the previous timeframe
449 // which matches our criterion
450 if (orbitsEarly > 0. && ti > 0) {
451 auto& prev_tf_range = true_indices[ti - 1];
452 // in this range search the smallest index which precedes
453 // timeframe ti by not more than "orbitsEarly" orbits
454 // (could probably use binary search, in case optimization becomes necessary)
455 int earlyOrbitIndex = -1; // init to start of this timeframe ... there may not be early orbits
456
457 // this is the orbit of the ti-th timeframe start
458 auto orbit_timeframe_start = startOrbit + ti * orbitsPerTF;
459
460 auto orbit_timeframe_early_fractional = 1. * orbit_timeframe_start - orbitsEarly;
461 auto orbit_timeframe_early_integral = static_cast<long>(std::floor(orbit_timeframe_early_fractional));
462
463 auto bc_early = (uint32_t)((orbit_timeframe_early_fractional - orbit_timeframe_early_integral) * o2::constants::lhc::LHCMaxBunches);
464
465 // this is the interaction record of the ti-th timeframe start
466 o2::InteractionRecord timeframe_start_record(0, orbit_timeframe_start);
467 // this is the interaction record in some previous timeframe after which interactions could still
468 // influence the ti-th timeframe according to orbitsEarly
469 o2::InteractionRecord timeframe_early_record(bc_early, orbit_timeframe_early_integral);
470
471 auto differenceInBCNS_max = timeframe_start_record.differenceInBCNS(timeframe_early_record);
472
473 for (int j = prev_tf_range.second; j >= prev_tf_range.first; --j) {
474 // determine difference in timing in NS; compare that with the limit given by orbitsEarly
475 auto timediff_NS = timeframe_start_record.differenceInBCNS(irecords[j]);
476 if (timediff_NS < differenceInBCNS_max) {
477 earlyOrbitIndex = j;
478 } else {
479 break;
480 }
481 }
482 std::get<2>(indices_with_early.back()) = earlyOrbitIndex;
483 }
484 }
485 return indices_with_early;
486}
487
488} // namespace
489
490void DigitizationContext::applyMaxCollisionFilter(std::vector<std::tuple<int, int, int>>& timeframeindices, long startOrbit, long orbitsPerTF, int maxColl, double orbitsEarly)
491{
492 // the idea is to go through each timeframe and throw away collisions beyond a certain count
493 // then the indices should be condensed
494
495 std::vector<std::vector<o2::steer::EventPart>> newparts;
496 std::vector<o2::InteractionTimeRecord> newrecords;
497
498 std::unordered_map<int, int> currMaxId; // the max id encountered for a source
499 std::unordered_map<int, std::unordered_map<int, int>> reIndexMap; // for each source, a map of old to new index for the event parts
500
501 if (maxColl == -1) {
502 maxColl = mEventRecords.size();
503 }
504
505 // the actual first actual timeframe
506 int first_timeframe = orbitsEarly > 0. ? 1 : 0;
507
508 // mapping of old to new indices
509 std::unordered_map<size_t, size_t> indices_old_to_new;
510
511 // now we can go through the structure timeframe by timeframe
512 for (int tf_id = first_timeframe; tf_id < timeframeindices.size(); ++tf_id) {
513 auto& tf_indices = timeframeindices[tf_id];
514
515 auto firstindex = std::get<0>(tf_indices); // .first;
516 auto lastindex = std::get<1>(tf_indices); // .second;
517 auto previndex = std::get<2>(tf_indices);
518
519 LOG(info) << "timeframe indices " << previndex << " : " << firstindex << " : " << lastindex;
520
521 int collCount = 0; // counting collisions within timeframe
522 // copy to new structure
523 for (int index = previndex >= 0 ? previndex : firstindex; index <= lastindex; ++index) {
524 if (collCount >= maxColl) {
525 continue;
526 }
527
528 // look if this index was already done?
529 // avoid duplicate entries in transformed records
530 if (indices_old_to_new.find(index) != indices_old_to_new.end()) {
531 continue;
532 }
533
534 // we put these events under a certain condition
535 bool keep = index < firstindex || collCount < maxColl;
536
537 if (!keep) {
538 continue;
539 }
540
541 if (index >= firstindex) {
542 collCount++;
543 }
544
545 // we must also make sure that we don't duplicate the records
546 // moreover some records are merely put as precoll of tf2 ---> so they shouldn't be part of tf1 in the final
547 // extraction, ouch !
548 // maybe we should combine the filter and individual tf extraction in one step !!
549 indices_old_to_new[index] = newrecords.size();
550 newrecords.push_back(mEventRecords[index]);
551 newparts.push_back(mEventParts[index]);
552
553 // reindex the event parts to achieve compactification (and initial linear increase)
554 for (auto& part : newparts.back()) {
555 auto source = part.sourceID;
556 auto entry = part.entryID;
557 // have we remapped this entry before?
558 if (reIndexMap.find(source) != reIndexMap.end() && reIndexMap[source].find(entry) != reIndexMap[source].end()) {
559 part.entryID = reIndexMap[source][entry];
560 } else {
561 // assign the next free index
562 if (currMaxId.find(source) == currMaxId.end()) {
563 currMaxId[source] = 0;
564 }
565 part.entryID = currMaxId[source];
566 // cache this assignment
567 reIndexMap[source][entry] = currMaxId[source];
568 currMaxId[source] += 1;
569 }
570 }
571 } // ends one timeframe
572
573 // correct the timeframe indices
574 if (indices_old_to_new.find(firstindex) != indices_old_to_new.end()) {
575 std::get<0>(tf_indices) = indices_old_to_new[firstindex]; // start
576 }
577 if (indices_old_to_new.find(lastindex) != indices_old_to_new.end()) {
578 std::get<1>(tf_indices) = indices_old_to_new[lastindex]; // end;
579 } else {
580 std::get<1>(tf_indices) = newrecords.size(); // end;
581 }
582 if (indices_old_to_new.find(previndex) != indices_old_to_new.end()) {
583 std::get<2>(tf_indices) = indices_old_to_new[previndex]; // previous or "early" index
584 }
585 }
586 // reassignment
587 mEventRecords = newrecords;
588 mEventParts = newparts;
589}
590
591std::vector<std::tuple<int, int, int>> DigitizationContext::calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly) const
592{
593 auto timeframeindices = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF, orbitsEarly);
594 LOG(info) << "Fixed " << timeframeindices.size() << " timeframes ";
595 for (auto p : timeframeindices) {
596 LOG(info) << std::get<0>(p) << " " << std::get<1>(p) << " " << std::get<2>(p);
597 }
598
599 return timeframeindices;
600}
601
602std::unordered_map<int, int> DigitizationContext::getCollisionIndicesForSource(int source) const
603{
604 // go through all collisions and pick those that have the give source
605 // then keep only the first collision index
606 std::unordered_map<int, int> result;
607 const auto& parts = getEventParts(false);
608 for (int collindex = 0; collindex < parts.size(); ++collindex) {
609 for (auto& eventpart : parts[collindex]) {
610 if (eventpart.sourceID == source) {
611 result[eventpart.entryID] = collindex;
612 }
613 }
614 }
615 return result;
616}
617
618int DigitizationContext::findSimPrefix(std::string const& prefix) const
619{
620 auto iter = std::find(mSimPrefixes.begin(), mSimPrefixes.end(), prefix);
621 if (iter != mSimPrefixes.end()) {
622 return std::distance(mSimPrefixes.begin(), iter);
623 }
624 return -1;
625}
626
627namespace
628{
629struct pair_hash {
630 template <class T1, class T2>
631 std::size_t operator()(const std::pair<T1, T2>& pair) const
632 {
633 return std::hash<T1>()(pair.first) ^ std::hash<T2>()(pair.second);
634 }
635};
636} // namespace
637
639{
640 // mapping of source x event --> index into mInteractionVertices
641 std::unordered_map<std::pair<int, int>, int, pair_hash> vertex_cache;
642
643 mInteractionVertices.clear();
644 int collcount = 0;
645
646 std::unordered_set<int> vset; // used to detect vertex incompatibilities
647 for (auto& coll : mEventParts) {
648 collcount++;
649 vset.clear();
650
651 // first detect if any of these entries already has an associated vertex
652 for (auto& part : coll) {
653 auto source = part.sourceID;
654 auto event = part.entryID;
655 auto cached_iter = vertex_cache.find(std::pair<int, int>(source, event));
656
657 if (cached_iter != vertex_cache.end()) {
658 vset.emplace(cached_iter->second);
659 }
660 }
661
662 // make sure that there is no conflict
663 if (vset.size() > 1) {
664 LOG(fatal) << "Impossible conflict during interaction vertex sampling";
665 }
666
667 int cacheindex = -1;
668 if (vset.size() == 1) {
669 // all of the parts need to be assigned the same existing vertex
670 cacheindex = *(vset.begin());
671 mInteractionVertices.push_back(mInteractionVertices[cacheindex]);
672 } else {
673 // we need to sample a new point
674 mInteractionVertices.emplace_back(meanv.sample());
675 cacheindex = mInteractionVertices.size() - 1;
676 }
677 // update the cache
678 for (auto& part : coll) {
679 auto source = part.sourceID;
680 auto event = part.entryID;
681 vertex_cache[std::pair<int, int>(source, event)] = cacheindex;
682 }
683 }
684}
685
686DigitizationContext DigitizationContext::extractSingleTimeframe(int timeframeid, std::vector<std::tuple<int, int, int>> const& timeframeindices, std::vector<int> const& sources_to_offset)
687{
688 DigitizationContext r; // make a return object
689 if (timeframeindices.size() == 0) {
690 LOG(error) << "Timeframe index structure empty; Returning empty object.";
691 return r;
692 }
693 r.mSimPrefixes = mSimPrefixes;
694 r.mMuBC = mMuBC;
695 r.mBCFilling = mBCFilling;
696 try {
697 auto tf_ranges = timeframeindices.at(timeframeid);
698
699 auto startindex = std::get<0>(tf_ranges);
700 auto endindex = std::get<1>(tf_ranges);
701 auto earlyindex = std::get<2>(tf_ranges);
702
703 if (earlyindex >= 0) {
704 startindex = earlyindex;
705 }
706 std::copy(mEventRecords.begin() + startindex, mEventRecords.begin() + endindex, std::back_inserter(r.mEventRecords));
707 std::copy(mEventParts.begin() + startindex, mEventParts.begin() + endindex, std::back_inserter(r.mEventParts));
708 if (mInteractionVertices.size() >= endindex) {
709 std::copy(mInteractionVertices.begin() + startindex, mInteractionVertices.begin() + endindex, std::back_inserter(r.mInteractionVertices));
710 }
711
712 // let's assume we want to fix the ids for source = source_id
713 // Then we find the first index that has this source_id and take the corresponding number
714 // as offset. Thereafter we subtract this offset from all known event parts.
715 auto perform_offsetting = [&r](int source_id) {
716 auto indices_for_source = r.getCollisionIndicesForSource(source_id);
717 int minvalue = std::numeric_limits<int>::max();
718 for (auto& p : indices_for_source) {
719 if (p.first < minvalue) {
720 minvalue = p.first;
721 }
722 }
723 // now fix them
724 for (auto& p : indices_for_source) {
725 auto index_into_mEventParts = p.second;
726 for (auto& part : r.mEventParts[index_into_mEventParts]) {
727 if (part.sourceID == source_id) {
728 part.entryID -= minvalue;
729 }
730 }
731 }
732 };
733 for (auto source_id : sources_to_offset) {
734 perform_offsetting(source_id);
735 }
736
737 } catch (std::exception) {
738 LOG(warn) << "No such timeframe id in collision context. Returing empty object";
739 }
740 // fix number of collisions
741 r.setNCollisions(r.mEventRecords.size());
742 return r;
743}
Definition of the Names Generator class.
uint64_t orbit
Definition RawEventData.h:6
int32_t i
GPUChain * chain
constexpr int p2()
constexpr int p1()
constexpr to accelerate the coordinates changing
uint32_t j
Definition RawData.h:0
static std::string getHitsFileName(DId d, const std::string_view prefix=STANDARDSIMPREFIX)
static std::string getMCKinematicsFileName(const std::string_view prefix=STANDARDSIMPREFIX)
Definition NameConf.h:46
static std::string getCollisionContextFileName(const std::string_view prefix="")
Definition NameConf.cxx:52
math_utils::Point3D< float > sample() const
sample a vertex from the MeanVertex parameters
Static class with identifiers, bitmasks and names for ALICE detectors.
Definition DetID.h:58
static constexpr const char * getName(ID id)
names of defined detectors
Definition DetID.h:145
static GRPObject * loadFrom(const std::string &grpFileName="")
DigitizationContext extractSingleTimeframe(int timeframeid, std::vector< std::tuple< int, int, int > > const &timeframeindices, std::vector< int > const &sources_to_offset)
bool checkVertexCompatibility(bool verbose=false) const
Check collision parts for vertex consistency.
void fillQED(std::string_view QEDprefix, int max_events, double qedrate)
add QED contributions to context, giving prefix; maximal event number and qed interaction rate
bool initSimChains(o2::detectors::DetID detid, std::vector< TChain * > &simchains) const
std::unordered_map< int, int > getCollisionIndicesForSource(int source) const
void printCollisionSummary(bool withQED=false, int truncateOutputTo=-1) const
int findSimPrefix(std::string const &prefix) const
void applyMaxCollisionFilter(std::vector< std::tuple< int, int, int > > &timeframeindices, long startOrbit, long orbitsPerTF, int maxColl, double orbitsEarly=0.)
void setSimPrefixes(std::vector< std::string > const &p)
bool initSimKinematicsChains(std::vector< TChain * > &simkinematicschains) const
o2::parameters::GRPObject const & getGRP() const
returns the GRP object associated to this context
void sampleInteractionVertices(o2::dataformats::MeanVertexObject const &v)
std::vector< std::vector< o2::steer::EventPart > > & getEventParts(bool withQED=false)
void saveToFile(std::string_view filename) const
std::vector< std::tuple< int, int, int > > calcTimeframeIndices(long startOrbit, long orbitsPerTF, double orbitsEarly=0.) const
get timeframe structure --> index markers where timeframe starts/ends/is_influenced_by
static DigitizationContext * loadFromFile(std::string_view filename="")
void setFirstIR(const o2::InteractionRecord &ir)
const o2::InteractionTimeRecord & generateCollisionTime()
void setBunchFilling(const BunchFilling &bc)
void setInteractionRate(float rateHz)
struct _cl_event * event
Definition glcorearb.h:2982
GLuint64EXT * result
Definition glcorearb.h:5662
GLuint entry
Definition glcorearb.h:5735
GLuint index
Definition glcorearb.h:781
GLdouble GLdouble right
Definition glcorearb.h:4077
GLdouble f
Definition glcorearb.h:310
GLsizei GLsizei GLchar * source
Definition glcorearb.h:798
GLint left
Definition glcorearb.h:1979
GLboolean * data
Definition glcorearb.h:298
GLboolean r
Definition glcorearb.h:1233
constexpr int LHCMaxBunches
std::string filename()
std::vector< o2::ctf::BufferType > vec
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"