Skip to content

Commit f2dbeb8

Browse files
committed
Intro of HitProcessingManager
1 parent 3c83db7 commit f2dbeb8

9 files changed

Lines changed: 494 additions & 4 deletions

File tree

DataFormats/simulation/include/SimulationDataFormat/MCInteractionRecord.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@ struct MCInteractionRecord {
2424
int orbit = 0; ///< LHC orbit
2525
int period = 0; ///< LHC period since beginning of run (if >0 -> time precision loss)
2626

27+
MCInteractionRecord() = default;
28+
2729
MCInteractionRecord(double tNS, int bcr = 0, int orb = 0, int per = 0) : timeNS(tNS), bc(bcr), orbit(orb), period(per)
2830
{
2931
}

Detectors/TPC/simulation/include/TPCSimulation/Point.h

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -90,6 +90,8 @@ class HitGroup : public o2::BaseHit {
9090
mHitsTVctr.emplace_back(time);
9191
mHitsEVctr.emplace_back(e);
9292
#endif
93+
mZAbsMax = std::max(std::abs(z), mZAbsMax);
94+
mZAbsMin = std::min(std::abs(z), mZAbsMin);
9395
}
9496

9597
size_t getSize() const {
@@ -135,8 +137,10 @@ class HitGroup : public o2::BaseHit {
135137
std::vector<float> mHitsXVctr;
136138
std::vector<float> mHitsYVctr;
137139
std::vector<float> mHitsZVctr;
138-
std::vector<float> mHitsTVctr;
139-
std::vector<short> mHitsEVctr;
140+
std::vector<float> mHitsTVctr;
141+
std::vector<short> mHitsEVctr;
142+
float mZAbsMin = 1E10; // minimal abs z position of all hits in this group
143+
float mZAbsMax = 0.; // maximal z position of all hits in this group
140144
#endif
141145
ClassDefNV(HitGroup, 1);
142146
};
Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
/*
2+
* filterHits.C
3+
*
4+
* Created on: Jan 30, 2018
5+
* Author: swenzel
6+
*/
7+
8+
#include <functional>
9+
#include <Steer/InteractionSampler.h>
10+
#include <TPCSimulation/Digitizer.h>
11+
#include <TPCSimulation/ToyDigitizerTask.h>
12+
13+
#include <cassert>
14+
15+
// an index to uniquely identify a single hit of TPC
16+
struct TPCHitGroupID {
17+
TPCHitGroupID() = default;
18+
TPCHitGroupID(int e, int gid) : entry{ e }, groupID{ gid } {}
19+
int entry = -1;
20+
int groupID = -1;
21+
};
22+
23+
TTree* tree;
24+
25+
void getHits(std::vector<std::vector<o2::TPC::HitGroup>*>& hitvectors, std::vector<TPCHitGroupID>& hitids,
26+
const std::vector<o2::MCInteractionRecord>& times, const char* branchname, float tmin /*NS*/,
27+
float tmax /*NS*/, std::function<float(float, float, float)>&& f)
28+
{
29+
// f is some function taking event time + z of hit and returns final "digit" time
30+
31+
// ingredients
32+
// *) tree
33+
34+
auto br = tree->GetBranch(branchname);
35+
if (!br) {
36+
std::cerr << "No branch found\n";
37+
return;
38+
}
39+
40+
auto nentries = br->GetEntries();
41+
hitvectors.resize(nentries, nullptr);
42+
// *) number entries/events in file
43+
44+
// do the filtering
45+
for (int entry = 0; entry < nentries; ++entry) {
46+
// std::cout << "ftimes " << f(times[entry].timeNS, 0, 250) << " " << f(times[entry].timeNS, 0, 0) << "\n";
47+
if (tmin > f(times[entry].timeNS, 0, 0)) {
48+
continue;
49+
}
50+
if (tmax < f(times[entry].timeNS, 0, 250)) {
51+
break;
52+
}
53+
// std::cout << "Keeping entry " << entry << " with time " << times[entry].timeNS << "\n";
54+
55+
// no filtering over hitgroups
56+
// std::vector<o2::TPC::HitGroup>* groups = nullptr;
57+
58+
// hitvectors.emplace_back(nullptr);
59+
60+
br->SetAddress(&hitvectors[entry]);
61+
br->GetEntry(entry);
62+
63+
int groupid = -1;
64+
auto groups = hitvectors[entry];
65+
for (auto& singlegroup : *groups) {
66+
groupid++;
67+
const auto zmax = singlegroup.mZAbsMax;
68+
const auto zmin = singlegroup.mZAbsMin;
69+
assert(zmin <= zmax);
70+
// auto tof = singlegroup.
71+
float tmaxtrack = f(times[entry].timeNS, 0., zmin);
72+
float tmintrack = f(times[entry].timeNS, 0., zmax);
73+
std::cout << tmintrack << " & " << tmaxtrack << "\n";
74+
assert(tmaxtrack > tmintrack);
75+
if (tmin > tmaxtrack || tmax < tmintrack) {
76+
std::cout << "DISCARDING " << groupid << " OF ENTRY " << entry << "\n";
77+
continue;
78+
}
79+
// need to record index of the group
80+
hitids.emplace_back(entry, groupid);
81+
}
82+
}
83+
}
84+
85+
// TPC hit selection lambda
86+
auto fTPC = [](float tNS, float tof, float z) {
87+
// returns time in NS
88+
return tNS + o2::TPC::Digitizer::getTime(z) * 1000 + tof;
89+
};
90+
91+
void loopHits(const std::vector<std::vector<o2::TPC::HitGroup>*>& hitvectors, const std::vector<TPCHitGroupID>& hitids,
92+
const std::vector<o2::MCInteractionRecord>& eventtimes)
93+
{
94+
float maxt = 0;
95+
float mint = 1E10;
96+
97+
for (auto& id : hitids) {
98+
auto entryvector = hitvectors[id.entry];
99+
auto& group = (*entryvector)[id.groupID];
100+
auto& MCrecord = eventtimes[id.entry];
101+
102+
for (size_t hitindex = 0; hitindex < group.getSize(); ++hitindex) {
103+
const auto& eh = group.getHit(hitindex);
104+
auto t = fTPC(MCrecord.timeNS, eh.GetTime(), eh.GetZ());
105+
maxt = std::max(t, maxt);
106+
mint = std::min(t, mint);
107+
// std::cout << "Hit from " << id.entry << " : " << eh.GetX() << " " << eh.GetY() << " " << eh.GetZ() << "\n";
108+
}
109+
}
110+
std::cout << " MINT " << mint << " MAXT " << maxt << "\n";
111+
}
112+
113+
void runTPCDigitization()
114+
{
115+
std::vector<std::vector<o2::TPC::HitGroup>*> hitvectors; // "TPCHitVector"
116+
std::vector<TPCHitGroupID> hitids; // "TPCHitIDs"
117+
std::vector<o2::MCInteractionRecord> times;
118+
119+
o2::TPC::ToyDigitizerTask task;
120+
121+
TFile file("o2sim.root");
122+
tree = (TTree*)file.Get("o2sim");
123+
124+
o2::steer::InteractionSampler sampler;
125+
times.resize(tree->GetEntries());
126+
sampler.generateCollisionTimes(times);
127+
128+
const auto TPCDRIFT = 100000;
129+
for (int driftinterval = 0;; ++driftinterval) {
130+
auto tmin = driftinterval * TPCDRIFT;
131+
auto tmax = (driftinterval + 1) * TPCDRIFT;
132+
133+
hitvectors.clear();
134+
hitids.clear();
135+
times.clear();
136+
137+
bool hasData = false;
138+
// loop over sectors
139+
for (int s = 0; s < 36; ++s) {
140+
std::stringstream sectornamestream;
141+
sectornamestream << "TPCHitsSector" << s;
142+
getHits(hitvectors, hitids, times, sectornamestream.str().c_str(), tmin, tmax, fTPC);
143+
144+
task.setData(&hitvectors, &times, s);
145+
task.Exec("");
146+
147+
hasData |= hitids.size();
148+
}
149+
150+
// condition to end:
151+
if (!hasData) {
152+
break;
153+
}
154+
}
155+
}
156+
157+
void filterHits()
158+
{
159+
// std::vector<std::vector<o2::TPC::HitGroup>*> hitvectors;
160+
// std::vector<TPCHitGroupID> hitids;
161+
// std::vector<o2::MCInteractionRecord> times;
162+
//
163+
// // 100us one TPC drift
164+
// // 100000ns one TPC drift
165+
// const auto TPCDRIFT = 100000;
166+
// auto tmin = 2.*TPCDRIFT;
167+
// auto tmax = 3.*TPCDRIFT;
168+
// getHits(hitvectors, hitids, times, "TPCHitsSector0", tmin, tmax, fTPC);
169+
//
170+
// std::cout << "obtained " << hitids.size() << " ids " << "\n";
171+
// std::cout << "spanning " << hitvectors.size() << " entries " << "\n";
172+
//
173+
// loopHits(hitvectors, hitids, times);
174+
Run();
175+
}
Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
/*
2+
* filterHits.C
3+
*
4+
* Created on: Jan 30, 2018
5+
* Author: swenzel
6+
*/
7+
#if !defined(__CLING__) || defined(__ROOTCLING__)
8+
#include <functional>
9+
#include <TPCSimulation/Digitizer.h>
10+
#include <TPCSimulation/DigitizerTask.h>
11+
#include <Steer/HitProcessingManager.h>
12+
#include "ITSMFTSimulation/Hit.h"
13+
#include <cassert>
14+
#endif
15+
16+
void getHits(const o2::steer::RunContext& context, std::vector<std::vector<o2::TPC::HitGroup>*>& hitvectors,
17+
std::vector<TPCHitGroupID>& hitids,
18+
19+
const char* branchname, float tmin /*NS*/, float tmax /*NS*/,
20+
std::function<float(float, float, float)>&& f)
21+
{
22+
// f is some function taking event time + z of hit and returns final "digit" time
23+
24+
// ingredients
25+
// *) tree
26+
27+
auto br = context.getBranch(branchname);
28+
if (!br) {
29+
std::cerr << "No branch found\n";
30+
return;
31+
}
32+
33+
auto& eventrecords = context.getEventRecords();
34+
35+
auto nentries = br->GetEntries();
36+
hitvectors.resize(nentries, nullptr);
37+
38+
// do the filtering
39+
for (int entry = 0; entry < nentries; ++entry) {
40+
if (tmin > f(eventrecords[entry].timeNS, 0, 0)) {
41+
continue;
42+
}
43+
if (tmax < f(eventrecords[entry].timeNS, 0, 250)) {
44+
break;
45+
}
46+
47+
br->SetAddress(&hitvectors[entry]);
48+
br->GetEntry(entry);
49+
50+
int groupid = -1;
51+
auto groups = hitvectors[entry];
52+
for (auto& singlegroup : *groups) {
53+
groupid++;
54+
const auto zmax = singlegroup.mZAbsMax;
55+
const auto zmin = singlegroup.mZAbsMin;
56+
assert(zmin <= zmax);
57+
// auto tof = singlegroup.
58+
float tmaxtrack = f(eventrecords[entry].timeNS, 0., zmin);
59+
float tmintrack = f(eventrecords[entry].timeNS, 0., zmax);
60+
std::cout << tmintrack << " & " << tmaxtrack << "\n";
61+
assert(tmaxtrack > tmintrack);
62+
if (tmin > tmaxtrack || tmax < tmintrack) {
63+
std::cout << "DISCARDING " << groupid << " OF ENTRY " << entry << "\n";
64+
continue;
65+
}
66+
// need to record index of the group
67+
hitids.emplace_back(entry, groupid);
68+
}
69+
}
70+
}
71+
72+
// TPC hit selection lambda
73+
auto fTPC = [](float tNS, float tof, float z) {
74+
// returns time in NS
75+
return tNS + o2::TPC::Digitizer::getTime(z) * 1000 + tof;
76+
};
77+
78+
void runTPCDigitization(const o2::steer::RunContext& context)
79+
{
80+
std::vector<std::vector<o2::TPC::HitGroup>*> hitvectorsleft; // "TPCHitVector"
81+
std::vector<std::vector<o2::TPC::HitGroup>*> hitvectorsright; // "TPCHitVector"
82+
std::vector<o2::TPC::TPCHitGroupID> hitidsleft; // "TPCHitIDs"
83+
std::vector<o2::TPC::TPCHitGroupID> hitidsright;
84+
85+
o2::TPC::DigitizerTask task;
86+
87+
const auto TPCDRIFT = 100000;
88+
for (int driftinterval = 0;; ++driftinterval) {
89+
auto tmin = driftinterval * TPCDRIFT;
90+
auto tmax = (driftinterval + 1) * TPCDRIFT;
91+
92+
hitvectorsleft.clear();
93+
hitidsleft.clear();
94+
hitvectorsright.clear();
95+
hitidsright.clear();
96+
97+
bool hasData = false;
98+
// loop over sectors
99+
for (int s = 0; s < 36; ++s) {
100+
task.setSector(s);
101+
std::stringstream sectornamestreamleft;
102+
sectornamestreamleft << "TPCHitsShiftedSector" << s;
103+
getHits(context, hitvectorsleft, hitidsleft, sectornamestreamleft.str().c_str(), tmin, tmax, fTPC);
104+
105+
std::stringstream sectornamestreamright;
106+
sectornamestreamright << "TPCHitsShiftedSector" << Shifted(s);
107+
getHits(context, hitvectorsright, hitidsright, sectornamestreamright.str().c_str(), tmin, tmax, fTPC);
108+
109+
task.setData(&hitvectorsleft, &hitvectorsright, &hitidsleft, &hitidsright, &context);
110+
task.Exec("");
111+
112+
hasData |= hitids.size();
113+
}
114+
115+
// condition to end:
116+
if (!hasData) {
117+
break;
118+
}
119+
task.FinishTask();
120+
}
121+
}
122+
123+
void runITSDigitization(const o2::steer::RunContext& context)
124+
{
125+
std::vector<o2::ITSMFT::Hit>* hitvector = nullptr;
126+
o2::ITS::DigitizerTask task(true);
127+
task.Init();
128+
auto br = context.getBranch("ITSHit");
129+
assert(br);
130+
br->SetAddress(&hitvector);
131+
for (int entry = 0; entry < context.getNEntries(); ++entry) {
132+
br->GetEntry(entry);
133+
task.setData(hitvector, &context);
134+
task.Exec("");
135+
}
136+
task.FinishTask();
137+
}
138+
139+
void runTPCDigitization_mgr()
140+
{
141+
auto& hitrunmgr = o2::steer::HitProcessingManager::instance();
142+
hitrunmgr.addInputFile("o2sim.root");
143+
hitrunmgr.registerRunFunction(runTPCDigitization);
144+
145+
// hitrunmgr.registerRunFunction(runITSDigitization);
146+
using Data_t = std::vector<o2::ITSMFT::Hit>;
147+
148+
o2::ITS::DigitizerTask task;
149+
TGeoManager::Import("O2geometry.root");
150+
task.Init();
151+
hitrunmgr.registerDefaultRunFunction(o2::steer::defaultRunFunction(task, "ITSHit"));
152+
153+
hitrunmgr.run();
154+
}

Detectors/TPC/simulation/src/DigitizerTask.cxx

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -72,8 +72,8 @@ InitStatus DigitizerTask::Init()
7272
}
7373

7474
/// Fetch the hits for the sector which is to be processed
75-
std::cout << "Processing sector " << mHitSector << " - loading HitSector " << mHitSector << " and "
76-
<< int(Sector::getRight(Sector(mHitSector))) << "\n";
75+
std::cout << "Processing sector " << mHitSector << " - loading HitSector "
76+
<< int(Sector::getLeft(Sector(mHitSector))) << " and " << mHitSector << "\n";
7777
std::stringstream sectornamestrleft, sectornamestrright;
7878
sectornamestrleft << "TPCHitsShiftedSector" << int(Sector::getLeft(Sector(mHitSector)));
7979
sectornamestrright << "TPCHitsShiftedSector" << mHitSector;

Steer/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ set(SRCS
1010
set(HEADERS
1111
include/${MODULE_NAME}/O2RunAna.h
1212
include/${MODULE_NAME}/InteractionSampler.h
13+
include/${MODULE_NAME}/HitProcessingManager.h
1314
)
1415

1516
set(LINKDEF src/SteerLinkDef.h)

0 commit comments

Comments
 (0)