Belle II Software light-2406-ragdoll
CreateFieldMapModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <geometry/modules/CreateFieldMapModule.h>
10#include <framework/geometry/BFieldManager.h>
11#include <framework/utilities/Utils.h>
12#include <framework/gearbox/Unit.h>
13
14#include <Math/Vector3D.h>
15#include <Math/VectorUtil.h>
16#include <TFile.h>
17#include <TH2D.h>
18#include <TTree.h>
19
20#include <boost/algorithm/string.hpp>
21#include <boost/format.hpp>
22
23using namespace Belle2;
24
25//-----------------------------------------------------------------
26// Register the Module
27//-----------------------------------------------------------------
28REG_MODULE(CreateFieldMap);
29
30//-----------------------------------------------------------------
31// Implementation
32//-----------------------------------------------------------------
34{
35 // Set module properties
36 setDescription("Create Field maps of the Belle II magnetic field used in the simulation");
37
38 // Parameter definitions
39 addParam("filename", m_filename, "ROOT filename for the map", std::string("FieldMap.root"));
40 addParam("type", m_type, "type of the fieldmap (xy, zx, zy, zr)", std::string("zx"));
41 addParam("nU", m_nU, "number of steps along first coordinate", 1800);
42 addParam("minU", m_minU, "minimum value for the first coordinate", -400.);
43 addParam("maxU", m_maxU, "maximum value for the first coordinate", 500.);
44 addParam("nV", m_nV, "number of steps along second coordinate", 1600);
45 addParam("minV", m_minV, "minimum value for the second coordinate", -400.);
46 addParam("maxV", m_maxV, "maximum value for the second coordinate", 400.);
47 addParam("phi", m_phi, "phi angle for the scan in radians", 0.);
48 addParam("wOffset", m_wOffset, "value of third coordinate", 0.);
49 addParam("nPhi", m_nPhi, "number of phi steps for zr averaging", m_nPhi);
50 addParam("saveAllPoints", m_createTree, "save all sampled points in a TTree, "
51 "WARNING: output can be huge", false);
52}
53
55{
56 if (m_filename.empty()) {
57 B2WARNING("CreateFieldMap: No Filename given, just sampling for fun");
58 }
59 boost::to_lower(m_type);
60 boost::trim(m_type);
61 if (m_type != "xy" && m_type != "zx" && m_type != "zy" && m_type != "zr") {
62 B2ERROR("CreateFieldMap type '" << m_type << "' not valid, use one of 'xy', 'zx', 'zy' or 'zr'");
63 }
64 if (m_nU < 0 || m_nV < 0) {
65 B2ERROR("Number of steps has to be positive");
66 }
67 if (m_maxU == m_minU) {
68 B2ERROR("Range for first coordinate is zero");
69 } else if (m_maxU < m_minU) {
70 std::swap(m_maxU, m_minU);
71 }
72 if (m_maxV == m_minV) {
73 B2ERROR("Range for second coordinate is zero");
74 } else if (m_maxV < m_minV) {
75 std::swap(m_maxV, m_minV);
76 }
77 if (m_type == "zr") {
78 if (m_minV < 0) {
79 B2ERROR("R values (minV, maxV) must be positive");
80 }
81 if (m_nPhi <= 0) {
82 B2ERROR("Number of steps in phi must be positive");
83 }
84 }
85}
86
88{
89 const bool save{!m_filename.empty()};
90 TFile* outfile{nullptr};
91 if (save) {
92 //Create histograms
93 outfile = new TFile(m_filename.c_str(), "RECREATE");
94 outfile->cd();
95 }
96
97 //Determine type of scan
98 EFieldTypes type = c_XY;
99 if (m_type == "zx")
100 type = c_ZX;
101 else if (m_type == "zy")
102 type = c_ZY;
103 else if (m_type == "zr")
104 type = c_ZR;
105
106 //some values needed for output
107 int lastPercent(-1);
108 uint64_t nSteps = m_nU * m_nV;
109 uint64_t curStep{0};
110 double startTime{0};
111 auto showProgress = [&]() {
112 if (curStep == 0) startTime = Utils::getClock();
113 const int64_t donePercent = 100 * ++curStep / nSteps;
114 if (donePercent > lastPercent) {
115 const double totalTime = Utils::getClock() - startTime;
116 const double perStep = totalTime / curStep;
117 if (donePercent == 100) {
118 B2INFO(boost::format("BField %s Scan: %d samples, %.3f us per sample, total: %.2f seconds")
119 % m_type % curStep
120 % (perStep / Unit::us) % (totalTime / Unit::s));
121 } else {
122 B2INFO(boost::format("BField %s Scan: %3d%%, %.3f us per sample")
123 % m_type % donePercent
124 % (perStep / Unit::us));
125 }
126 lastPercent = donePercent;
127 }
128 };
129
130 struct { float x{0}, y{0}, z{0}, bx{0}, by{0}, bz{0}; } field_point;
131 TTree* all_values{nullptr};
132 if (save && m_createTree) {
133 all_values = new TTree("bfield_values", "All B field values");
134 all_values->Branch("x", &field_point.x, "x/F");
135 all_values->Branch("y", &field_point.y, "y/F");
136 all_values->Branch("z", &field_point.z, "z/F");
137 all_values->Branch("bx", &field_point.bx, "bx/F");
138 all_values->Branch("by", &field_point.by, "by/F");
139 all_values->Branch("bz", &field_point.bz, "bz/F");
140 }
141 auto fillTree = [&](const ROOT::Math::XYZVector & p, const ROOT::Math::XYZVector & b) {
142 if (!all_values) return;
143 field_point.x = p.X();
144 field_point.y = p.Y();
145 field_point.z = p.Z();
146 field_point.bx = b.X();
147 field_point.by = b.Y();
148 field_point.bz = b.Z();
149 all_values->Fill();
150 };
151
152 if (type == c_ZR) {
153 nSteps *= m_nPhi;
154 TH2D* h_b = new TH2D("B", "$B$ average;$z$/cm;$r$/cm", m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
155 TH2D* h_br = new TH2D("Br", "$B_r$ average;$z$/cm;$r$/cm", m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
156 TH2D* h_bz = new TH2D("Bz", "$B_z$ average;$z$/cm;$r$/cm", m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
157 for (int iU = 0; iU < m_nU; ++iU) {
158 for (int iV = 0; iV < m_nV; ++iV) {
159 //find value for first and second coordinate
160 const double u = h_b->GetXaxis()->GetBinCenter(iU + 1);
161 const double v = h_b->GetYaxis()->GetBinCenter(iV + 1);
162 //Determine global coordinates
163 for (int iPhi = 0; iPhi < m_nPhi; ++iPhi) {
164 ROOT::Math::XYZVector pos(v, 0, u);
165 pos = ROOT::Math::VectorUtil::RotateZ(pos, 2 * M_PI * iPhi / m_nPhi);
166 //Obtain magnetic field
167 ROOT::Math::XYZVector bfield = BFieldManager::getFieldInTesla(pos);
168 //And fill histograms
169 if (save) {
170 h_br->Fill(u, v, bfield.Rho());
171 h_bz->Fill(u, v, bfield.Z());
172 h_b->Fill(u, v, bfield.R());
173 fillTree(pos, bfield);
174 }
175 showProgress();
176 }
177 }
178 }
179 //Write histograms and close file.
180 if (save) {
181 for (TH2D* h : {h_br, h_bz, h_b}) {
182 h->Scale(1. / m_nPhi);
183 h->Write();
184 }
185 }
186 } else {
187 const std::string nu = m_type.substr(0, 1);
188 const std::string nv = m_type.substr(1, 1);
189 TH2D* h_b = new TH2D("B", ("$B$;$" + nu + "$/cm;$" + nv + "$/cm").c_str(), m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
190 TH2D* h_bx = new TH2D("Bx", ("$B_x$;$" + nu + "$/cm;$" + nv + "$/cm").c_str(), m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
191 TH2D* h_by = new TH2D("By", ("$B_y$;$" + nu + "$/cm;$" + nv + "$/cm").c_str(), m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
192 TH2D* h_bz = new TH2D("Bz", ("$B_z$;$" + nu + "$/cm;$" + nv + "$/cm").c_str(), m_nU, m_minU, m_maxU, m_nV, m_minV, m_maxV);
193
194 //Loop over all bins
195 for (int iU = 0; iU < m_nU; ++iU) {
196 for (int iV = 0; iV < m_nV; ++iV) {
197 ROOT::Math::XYZVector pos(0, 0, 0);
198 //find value for first and second coordinate
199 const double u = h_bx->GetXaxis()->GetBinCenter(iU + 1);
200 const double v = h_bx->GetYaxis()->GetBinCenter(iV + 1);
201 //Determine global coordinates
202 switch (type) {
203 case c_XY:
204 pos.SetXYZ(u, v, m_wOffset);
205 break;
206 case c_ZX:
207 pos.SetXYZ(v, m_wOffset, u);
208 break;
209 case c_ZY:
210 pos.SetXYZ(m_wOffset, v, u);
211 break;
212 default:
213 break;
214 }
215 pos = ROOT::Math::VectorUtil::RotateZ(pos, m_phi);
216 //Obtain magnetic field
217 ROOT::Math::XYZVector bfield = BFieldManager::getFieldInTesla(pos);
218 //And fill histograms
219 if (save) {
220 h_bx->Fill(u, v, bfield.X());
221 h_by->Fill(u, v, bfield.Y());
222 h_bz->Fill(u, v, bfield.Z());
223 h_b->Fill(u, v, bfield.R());
224 fillTree(pos, bfield);
225 }
226 showProgress();
227 }
228 }
229 //Write histograms.
230 if (save) {
231 for (TH2D* h : {h_bx, h_by, h_bz, h_b}) {
232 h->Write();
233 }
234 }
235 }
236 if (save) {
237 if (all_values) all_values->Write();
238 outfile->Close();
239 delete outfile;
240 }
241 //histograms seem to be deleted when file is closed
242}
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:61
int m_nPhi
number of steps in phi
virtual void initialize() override
Check input parameters.
double m_minV
start value for the first coordinate
double m_wOffset
offset on the third coordinate
double m_minU
start value for the first coordinate
double m_phi
phi rotation when sampling magnetic field (value of pi will convert ZX to ZY scan)
int m_nV
number of steps along the second coordinate
int m_nU
number of steps along the first coordinate
std::string m_type
type of the fieldmap (zx, zy, zr)
double m_maxV
end value for the first coordinate
virtual void beginRun() override
Create the fieldmap.
std::string m_filename
output filename for the fieldmap
EFieldTypes
Types of Fieldmap to be created.
@ c_ZR
scan along Z and R, averaging over phi
double m_maxU
end value for the first coordinate
bool m_createTree
if true create a TTree with all sampled points
CreateFieldMapModule()
Constructor: Sets the description, the properties and the parameters of the module.
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
static const double us
[microsecond]
Definition: Unit.h:97
static const double s
[second]
Definition: Unit.h:95
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
double getClock()
Return current value of the real-time clock.
Definition: Utils.cc:66
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24