Skip to content

Commit c64a091

Browse files
authored
O2sim: Custom stepping + configurable geometry tracking cuts (#1712)
O2sim: Custom stepping + configurable geometry tracking cuts This commit is introducing the possibility to have tracking cuts based on geometry and other stepping status. Custom filtering is possible inside our overriden stepping function. Introducing a macro to compare information on hits, potentially usefull to judge if a tracking cut is not changing behaviour.
1 parent d8767ca commit c64a091

7 files changed

Lines changed: 293 additions & 5 deletions

File tree

Common/SimConfig/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,12 @@ set(SRCS
66
src/SimConfig.cxx
77
src/ConfigurableParam.cxx
88
src/ConfigurableParamHelper.cxx
9+
src/SimCutParams.cxx
910
)
1011

1112
set(HEADERS
1213
include/${MODULE_NAME}/SimConfig.h
14+
include/${MODULE_NAME}/SimCutParams.h
1315
include/${MODULE_NAME}/ConfigurableParam.h
1416
include/${MODULE_NAME}/ConfigurableParamHelper.h
1517
)
Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
#ifndef O2_SIMCONFIG_SIMCUTPARAM_H_
12+
#define O2_SIMCONFIG_SIMCUTPARAM_H_
13+
14+
#include "SimConfig/ConfigurableParam.h"
15+
#include "SimConfig/ConfigurableParamHelper.h"
16+
17+
namespace o2
18+
{
19+
namespace conf
20+
{
21+
// parameters to influence/set the tracking cuts in simulation
22+
// (mostly used in O2MCApplication stepping)
23+
struct SimCutParams : public o2::conf::ConfigurableParamHelper<SimCutParams> {
24+
bool stepFiltering = true; // if we activate the step filtering in O2BaseMCApplication
25+
26+
double maxRTracking = 1E20; // max R tracking cut in cm (in the VMC sense) -- applied in addition to cutting in the stepping function
27+
double maxAbsZTracking = 1E20; // max |Z| tracking cut in cm (in the VMC sense) -- applied in addition to cutting in the stepping function
28+
29+
O2ParamDef(SimCutParams, "SimCutParams");
30+
};
31+
} // namespace conf
32+
} // namespace o2
33+
34+
#endif /* O2_SIMCONFIG_SIMCUTPARAM_H_ */

Common/SimConfig/src/SimConfigLinkDef.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,4 +17,6 @@
1717
#pragma link C++ class o2::conf::SimConfig+;
1818
#pragma link C++ class o2::conf::SimConfigData+;
1919

20+
#pragma link C++ class o2::conf::SimCutParams + ;
21+
#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::conf::SimCutParams > +;
2022
#endif
Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
/*
2+
* SimCutParams.cxx
3+
*
4+
* Created on: Feb 18, 2019
5+
* Author: sandro
6+
*/
7+
8+
// Copyright CERN and copyright holders of ALICE O2. This software is
9+
// distributed under the terms of the GNU General Public License v3 (GPL
10+
// Version 3), copied verbatim in the file "COPYING".
11+
//
12+
// See http://alice-o2.web.cern.ch/license for full licensing information.
13+
//
14+
// In applying this license CERN does not waive the privileges and immunities
15+
// granted to it by virtue of its status as an Intergovernmental Organization
16+
// or submit itself to any jurisdiction.
17+
18+
#include "SimConfig/SimCutParams.h"
19+
O2ParamImpl(o2::conf::SimCutParams);

Steer/include/Steer/O2MCApplicationBase.h

Lines changed: 16 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,10 @@
1818
#ifndef STEER_INCLUDE_STEER_O2MCAPPLICATIONBASE_H_
1919
#define STEER_INCLUDE_STEER_O2MCAPPLICATIONBASE_H_
2020

21-
#include "FairMCApplication.h"
21+
#include <FairMCApplication.h>
2222
#include "Rtypes.h" // for Int_t, Bool_t, Double_t, etc
23+
#include <TVirtualMC.h>
24+
#include "SimConfig/SimCutParams.h"
2325

2426
namespace o2
2527
{
@@ -32,14 +34,23 @@ namespace steer
3234
class O2MCApplicationBase : public FairMCApplication
3335
{
3436
public:
35-
using FairMCApplication::FairMCApplication;
37+
O2MCApplicationBase() : FairMCApplication(), mCutParams(o2::conf::SimCutParams::Instance()) {}
38+
O2MCApplicationBase(const char* name, const char* title, TObjArray* ModList, const char* MatName) : FairMCApplication(name, title, ModList, MatName), mCutParams(o2::conf::SimCutParams::Instance())
39+
{
40+
}
41+
3642
~O2MCApplicationBase() override = default;
3743

44+
void Stepping() override;
45+
3846
// specific implementation of our hard geometry limits
39-
double TrackingRmax() const override { return 1E20; }
40-
double TrackingZmax() const override { return 1E20; }
47+
double TrackingRmax() const override { return mCutParams.maxRTracking; }
48+
double TrackingZmax() const override { return mCutParams.maxAbsZTracking; }
49+
50+
protected:
51+
o2::conf::SimCutParams const& mCutParams; // reference to parameter system
4152

42-
ClassDefOverride(O2MCApplicationBase, 1) //Interface to MonteCarlo application
53+
ClassDefOverride(O2MCApplicationBase, 1)
4354
};
4455

4556
} // end namespace steer

Steer/src/O2MCApplication.cxx

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,14 @@ void TypedVectorAttach(const char* name, FairMQChannel& channel, FairMQParts& pa
4343
}
4444
}
4545

46+
void O2MCApplicationBase::Stepping()
47+
{
48+
// dispatch first to stepping function in FairRoot
49+
FairMCApplication::Stepping();
50+
51+
// now the right place to put our specific step/filtering criteria:
52+
}
53+
4654
void O2MCApplication::initLate()
4755
{
4856
o2::utils::ShmManager::Instance().occupySegment();

macro/compareHits.C

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
#if !defined(__CLING__) || defined(__ROOTCLING__)
2+
#include "TFile.h"
3+
#include "TTree.h"
4+
#include "ITSMFTSimulation/Hit.h"
5+
#include "TOFSimulation/Detector.h"
6+
#include "EMCALBase/Hit.h"
7+
#include "TRDSimulation/Detector.h"
8+
#endif
9+
10+
template <typename Hit, typename Accumulator>
11+
Accumulator analyse(TTree* tr, const char* brname)
12+
{
13+
Accumulator prop;
14+
auto br = tr->GetBranch(brname);
15+
if (!br) {
16+
return prop;
17+
}
18+
auto entries = br->GetEntries();
19+
std::vector<Hit>* hitvector = nullptr;
20+
br->SetAddress(&hitvector);
21+
22+
for (int i = 0; i < entries; ++i) {
23+
br->GetEntry(i);
24+
for (auto& hit : *hitvector) {
25+
prop.addHit(hit);
26+
}
27+
}
28+
prop.normalize();
29+
return prop;
30+
};
31+
32+
template <typename T>
33+
struct HitStats {
34+
int NHits = 0;
35+
double XAvg = 0.; // avg 1st moment
36+
double YAvg = 0.;
37+
double ZAvg = 0.;
38+
double X2Avg = 0.; // avg 2nd moment
39+
double Y2Avg = 0.;
40+
double Z2Avg = 0.;
41+
double EAvg = 0.; // average total energy
42+
double E2Avg = 0.;
43+
double TAvg = 0.; // average T
44+
double T2Avg = 0.; // average T^2
45+
46+
void print() const
47+
{
48+
std::cout << NHits << " "
49+
<< XAvg << " "
50+
<< YAvg << " "
51+
<< ZAvg << " "
52+
<< X2Avg << " "
53+
<< Y2Avg << " "
54+
<< Z2Avg << " "
55+
<< EAvg << " "
56+
<< E2Avg << " "
57+
<< TAvg << " "
58+
<< T2Avg << "\n";
59+
}
60+
61+
// adds a hit to the statistics
62+
void addHit(T const& hit)
63+
{
64+
NHits++;
65+
auto x = hit.GetX();
66+
XAvg += x;
67+
X2Avg += x * x;
68+
auto y = hit.GetY();
69+
YAvg += y;
70+
Y2Avg += y * y;
71+
auto z = hit.GetZ();
72+
ZAvg += z;
73+
Z2Avg += z * z;
74+
auto e = hit.GetEnergyLoss();
75+
EAvg += e;
76+
E2Avg += e * e;
77+
auto t = hit.GetTime();
78+
TAvg += t;
79+
T2Avg += t * t;
80+
}
81+
82+
void normalize()
83+
{
84+
XAvg /= NHits;
85+
YAvg /= NHits;
86+
ZAvg /= NHits;
87+
X2Avg /= NHits;
88+
Y2Avg /= NHits;
89+
Z2Avg /= NHits;
90+
EAvg /= NHits;
91+
E2Avg /= NHits;
92+
TAvg /= NHits;
93+
T2Avg /= NHits;
94+
}
95+
};
96+
97+
struct ITSHitStats {
98+
int NHits = 0;
99+
double XAvg = 0.; // avg 1st moment
100+
double YAvg = 0.;
101+
double ZAvg = 0.;
102+
double X2Avg = 0.; // avg 2nd moment
103+
double Y2Avg = 0.;
104+
double Z2Avg = 0.;
105+
double EAvg = 0.; // average total energy
106+
double E2Avg = 0.;
107+
double TAvg = 0.; // average T
108+
double T2Avg = 0.; // average T^2
109+
110+
void print() const
111+
{
112+
std::cout << NHits << " "
113+
<< XAvg << " "
114+
<< YAvg << " "
115+
<< ZAvg << " "
116+
<< X2Avg << " "
117+
<< Y2Avg << " "
118+
<< Z2Avg << " "
119+
<< EAvg << " "
120+
<< E2Avg << " "
121+
<< TAvg << " "
122+
<< T2Avg << "\n";
123+
}
124+
125+
// adds a hit to the statistics
126+
void addHit(o2::ITSMFT::Hit const& hit)
127+
{
128+
NHits++;
129+
auto x = hit.GetStartX();
130+
XAvg += x;
131+
X2Avg += x * x;
132+
auto y = hit.GetStartY();
133+
YAvg += y;
134+
Y2Avg += y * y;
135+
auto z = hit.GetStartZ();
136+
ZAvg += z;
137+
Z2Avg += z * z;
138+
auto e = hit.GetTotalEnergy();
139+
EAvg += e;
140+
E2Avg += e * e;
141+
auto t = hit.GetTime();
142+
TAvg += t;
143+
T2Avg += t * t;
144+
}
145+
146+
void normalize()
147+
{
148+
XAvg /= NHits;
149+
YAvg /= NHits;
150+
ZAvg /= NHits;
151+
X2Avg /= NHits;
152+
Y2Avg /= NHits;
153+
Z2Avg /= NHits;
154+
EAvg /= NHits;
155+
E2Avg /= NHits;
156+
TAvg /= NHits;
157+
T2Avg /= NHits;
158+
}
159+
}; // end struct
160+
161+
// do comparison for ITS
162+
void checkITS(TTree* reftree, TTree* testtree)
163+
{
164+
auto refresult = analyse<o2::ITSMFT::Hit, ITSHitStats>(reftree, "ITSHit");
165+
refresult.print();
166+
auto testresult = analyse<o2::ITSMFT::Hit, ITSHitStats>(testtree, "ITSHit");
167+
testresult.print();
168+
}
169+
170+
// do comparison for TOF
171+
void checkTOF(TTree* reftree, TTree* testtree)
172+
{
173+
auto refresult = analyse<o2::tof::HitType, HitStats<o2::tof::HitType>>(reftree, "TOFHit");
174+
refresult.print();
175+
auto testresult = analyse<o2::tof::HitType, HitStats<o2::tof::HitType>>(testtree, "TOFHit");
176+
testresult.print();
177+
}
178+
179+
// do comparison for EMC
180+
void checkEMC(TTree* reftree, TTree* testtree)
181+
{
182+
auto refresult = analyse<o2::EMCAL::Hit, HitStats<o2::EMCAL::Hit>>(reftree, "EMCHit");
183+
refresult.print();
184+
auto testresult = analyse<o2::EMCAL::Hit, HitStats<o2::EMCAL::Hit>>(testtree, "EMCHit");
185+
testresult.print();
186+
}
187+
188+
// do comparison for TRD
189+
void checkTRD(TTree* reftree, TTree* testtree)
190+
{
191+
auto refresult = analyse<o2::trd::HitType, HitStats<o2::trd::HitType>>(reftree, "TRDHit");
192+
refresult.print();
193+
auto testresult = analyse<o2::trd::HitType, HitStats<o2::trd::HitType>>(testtree, "TRDHit");
194+
testresult.print();
195+
}
196+
197+
// Simple macro to compare properties of simulated hits
198+
// of two runs (for example reference run and test run).
199+
// Used for instance to validate MC optimizations or to study
200+
// the effect on the physics at the lowest level.
201+
void compareHits(const char* reffilename = "o2sim.root", const char* testfilename = "o2sim.root")
202+
{
203+
TFile rf(reffilename, "OPEN");
204+
TFile tf(testfilename, "OPEN");
205+
auto reftree = (TTree*)rf.Get("o2sim");
206+
auto testtree = (TTree*)tf.Get("o2sim");
207+
208+
checkITS(reftree, testtree);
209+
checkTOF(reftree, testtree);
210+
checkEMC(reftree, testtree);
211+
checkTRD(reftree, testtree);
212+
}

0 commit comments

Comments
 (0)