npdet_fields.cxx
Go to the documentation of this file.
1 #include "Math/Vector3D.h"
2 #include "Math/Vector4D.h"
3 #include "Math/VectorUtil.h"
4 #include "TCanvas.h"
5 #include "TLegend.h"
6 #include "TMath.h"
7 #include "TRandom3.h"
8 #include "TFile.h"
9 #include "TH1F.h"
10 #include "TH2F.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 #include "TTree.h"
14 #include "TF1.h"
15 
16 #include <algorithm>
17 #include <iterator>
18 #include <tuple>
19 #include <vector>
20 
21 #include "DD4hep/DD4hepUnits.h"
22 #include "DD4hep/Detector.h"
23 #include "DDG4/Geant4Data.h"
24 #include "DDRec/CellIDPositionConverter.h"
25 #include "DDRec/Surface.h"
26 #include "DDRec/SurfaceManager.h"
27 #include "TApplication.h"
28 #include "TGraph.h"
29 #include "TMultiGraph.h"
30 
31 #include "Math/DisplacementVector3D.h"
32 
33 #include <iostream>
34 #include <string>
35 #include "clipp.h"
36 
37 using namespace clipp;
38 using namespace ROOT::Math;
39 using namespace dd4hep;
40 
41 struct settings {
42  bool success = false;
43  std::string infile = "";
44  ROOT::Math::XYZVector start_position{0.0,0.0,0.0};
45  ROOT::Math::XYZVector direction{0.0,0.0,1.0};
46  double step_size = 0.5*cm;
47  bool rec = false;
48  bool utf16 = false;
49  bool search = false;
50  int search_level = 0;
51  std::vector<std::string> field_comps = {"Bz"};
52  std::string format = "json";
53  std::vector<std::string> axes;
54  int Nsteps = 200;
55  double x0 = 0.0;
56  double y0 = 0.0;
57  double z0 = 0.0;
58  double x1 = 0.0;
59  double y1 = 0.0;
60  double z1 = 0.0;
61  double x2 = 0.0;
62  double y2 = 0.0;
63  double z2 = 1.0;
64  bool by_step_size = false;
65  bool with_end_point = false;
66  bool verbose = false;
67 };
68 
69 settings cmdline_settings(int argc, char* argv[]) {
70  settings s;
71 
72  auto helpMode = "help mode:" % (
73  command("help")
74  );
75 
76  auto drawMode =
77  "draw mode:" %
78  (command("draw"), // values("component").set(s.field_comps),
79  option("--Nsteps") & number("Nsteps", s.step_size),
80  option("--step") & number("step", s.step_size),
81  option("--start") & (number("x", s.x0).if_missing([] { std::cout << "x missing!\n"; }),
82  number("y", s.y0).if_missing([] { std::cout << "y missing!\n"; }),
83  number("z", s.z0).if_missing([] { std::cout << "z missing!\n"; })),
84  (option("--end").set(s.with_end_point, true) & (number("x1", s.x1).if_missing([] { std::cout << "x1 missing!\n"; }),
85  number("y1", s.y1).if_missing([] { std::cout << "y1 missing!\n"; }),
86  number("z1", s.z1).if_missing([] { std::cout << "z1 missing!\n"; })) |
87  option("--direction") &
88  (number("x2", s.x2).if_missing([] { std::cout << "x2 missing!\n"; }),
89  number("y2", s.y2).if_missing([] { std::cout << "y2 missing!\n"; }),
90  number("z2", s.z2).if_missing([] { std::cout << "z2 missing!\n"; })))
91  // required("--vs","--Vs") & values("axes",s.axes ), % "axes"
92  //option("-l", "--level") & integers("level", s.search_level) % "search level"
93  );
94 
95  auto printMode = "integral mode (not implemented):" % (
96  command("integrate")
97  //option("--all") % "copy all"
98  //option("--replace") % "replace existing files",
99  //option("-f", "--force") % "don't ask for confirmation"
100  );
101 
102  auto firstOpt = "user interface options:" % (
103  //joinable(
104  option("-v", "--verbose").set(s.verbose) % "show detailed output",
105  option("-h", "--help") % "show help"
106  //option("-i", "--interactive") % "use interactive mode (not implemented)"
107  // )
108  );
109  auto lastOpt = value("file", s.infile).if_missing([] {
110  std::cout << "You need to provide an input xml filename as the last argument!\n";
111  }) % "compact description should probly have a \"fields\" sect4ion with more than one magnetic "
112  "field defined.";
113 
114  auto cli = (
115  firstOpt,
116  drawMode,// (| printMode | command("list")),
117  lastOpt
118  );
119  //auto cli = (
120  // command("help")
121  // command("search")
122  // );
123  if(!parse(argc, argv, cli)) {
124  s.success = false;
125  std::cout << make_man_page(cli, argv[0])
126  .prepend_section("DESCRIPTION", " Tool for quickly looking at magnetic fields.")
127  .append_section("EXAMPLES"," npdet_fields --start 0 0 -600 solid_sidis.xml\n");
128  return s;
129  }
130 
131  s.start_position.SetXYZ(s.x0*cm,s.y0*cm,s.z0*cm);
132  if(s.with_end_point){
133  s.direction = ROOT::Math::XYZVector(s.x1 * cm, s.y1 * cm, s.z1 * cm) - s.start_position;
134  s.step_size = s.direction.R()/s.Nsteps;
135  s.direction = s.direction.Unit();
136  } else {
137  s.direction.SetXYZ(s.x2 * cm, s.y2 * cm, s.z2 * cm);
138  //s.step_size = s.direction.R()/s.Nsteps;
139  s.direction = s.direction.Unit();
140  }
141  s.success = true;
142 
143  std::cout << "start [cm] : " << s.start_position/cm <<"\n";
144  std::cout << "dir [cm] : " << s.direction/cm <<"\n";
145  std::cout << "step [cm] : " << s.step_size/cm <<"\n";
146  std::cout << "nstep : " << s.Nsteps <<"\n";
147 
148  return s;
149 }
150 //______________________________________________________________________________
151 
152 TGraph* build_1D_field_graph(dd4hep::Detector& detector, const settings& s,
153  std::function<double(ROOT::Math::XYZVector)> B_comp);
154 
158 void mag_field(dd4hep::Detector& detector){
159 
160  double z0 = 0.0;
161 
162  double x_min = -100.0;
163  double x_max = 100.0;
164  double y_min = -100.0;
165  double y_max = 100.0;
166  double z_min = -100.0;
167  double z_max = 100.0;
168  int nx_bins = 200;
169  int ny_bins = 200;
170  int nz_bins = 200;
171  double dx = (x_max - x_min) / nx_bins;
172  double dy = (y_max - y_min) / ny_bins;
173  double dz = (z_max - z_min) / nz_bins;
174 
175  double y0 = 20.0 * cm;
176 
177  TH2F* hBx = new TH2F("hBx","; z [cm]; x [cm]; B_{x} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
178  TH2F* hBz = new TH2F("hBz","; z [cm]; x [cm]; B_{z} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
179  TH2F* hBr = new TH2F("hBr","; z [cm]; x [cm]; B_{r} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
180  TH2F* hBphi = new TH2F("hBphi","; z [cm]; x [cm]; B_{#phi} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
181 
182 
183  for(int ix = 0; ix< nx_bins; ix++) {
184  for(int iz = 0; iz< nz_bins; iz++) {
185 
186  double x = (x_min + ix*dx)*cm;
187  double z = (z_min + iz*dz)*cm;
188 
189  XYZVector pos{ x, y0, z };
190  XYZVector field = detector.field().magneticField(pos);
191 
192  hBz->Fill(z,x, field.z()/tesla);
193  hBx->Fill(z,x, field.x()/tesla);
194  hBr->Fill(z,x, field.rho()/tesla);
195  hBphi->Fill(z,x, field.phi()/tesla);
196  }
197 
198  }
199 
200  auto c = new TCanvas();
201  c->Divide(2,2);
202 
203  c->cd(1);
204  hBz->Draw("colz");
205 
206  c->cd(2);
207  hBx->Draw("colz");
208 
209  c->cd(3);
210  hBr->Draw("colz");
211 
212  c->cd(4);
213  hBphi->Draw("colz");
214 }
215 
216 int main (int argc, char *argv[]) {
217 
218  settings s = cmdline_settings(argc,argv);
219  if( !s.success ) {
220  return 1;
221  }
222 
223  // ------------------------
224  // TODO: CLI Checks
225 
226  using namespace dd4hep;
227 
228  // -------------------------
229  // Get the DD4hep instance
230  // Load the compact XML file
231  dd4hep::Detector& detector = dd4hep::Detector::getInstance();
232  detector.fromCompact(s.infile);
233 
234  int root_argc = 0;
235  char *root_argv[1];
236  argv[0] = "npdet_fields";
237 
238  TApplication theApp("tapp", &root_argc, root_argv);
239 
240  TMultiGraph* mg = new TMultiGraph();
241  for(const auto& comp : s.field_comps){
242  std::cout << comp << "\n";
243 
244  auto gr = build_1D_field_graph(detector, s, [](ROOT::Math::XYZVector B) { return B.z(); });
245  gr->SetLineColor(2);
246  gr->SetLineWidth(2);
247  gr->SetFillColor(0);
248  gr->SetTitle("B_{z}");
249  mg->Add(gr,"l");
250 
251  gr = build_1D_field_graph(detector, s,
252  [&](ROOT::Math::XYZVector B) { return B.Rho(); });
253  gr->SetLineColor(1);
254  gr->SetLineWidth(2);
255  gr->SetFillColor(0);
256  gr->SetTitle("B_{#perp}");
257  mg->Add(gr,"l");
258 
259  gr = build_1D_field_graph(detector, s, [&](ROOT::Math::XYZVector B) {
260  return (B - B.Unit() * (B.Dot(s.direction))).r();
261  });
262  gr->SetLineColor(4);
263  gr->SetLineWidth(2);
264  gr->SetFillColor(0);
265  gr->SetTitle("B_{perp}");
266  //mg->Add(gr,"l");
267 
268  }
269 
270  auto c = new TCanvas();
271  mg->Draw("a");
272  //mag_field( detector );
273  c->BuildLegend();
274 
275  std::cout << "input file : " << s.infile << "\n";
276  std::cout << " start : " << s.start_position << "\n";
277  std::cout << " direction : " << s.direction << "\n";
278  std::cout << " step : " << s.step_size << "\n";
279  std::cout << "\n";
280 
281 
282  theApp.Run();
283  return 0;
284 }
285 //______________________________________________________________________________
286 
287 TGraph* build_1D_field_graph(dd4hep::Detector& detector, const settings& s,
288  std::function<double(ROOT::Math::XYZVector)> B_comp) {
289 
290  int N_steps = s.Nsteps;
291  std::cout << "Steps : " << N_steps << "\n";
292 
293  auto dir = s.direction.Unit();
294 
295  auto gr = new TGraph();
296  for (int i = 0; i < N_steps; i++) {
297  ROOT::Math::XYZVector pos = s.start_position + double(i) * s.step_size * dir;
298  ROOT::Math::XYZVector field = detector.field().magneticField(pos);
299 
300  if (s.verbose) {
301  std::cout << "x : " << pos / mm << "\n";
302  std::cout << "B : " << field / tesla << "\n";
303  }
304  gr->SetPoint(i, pos.z(), B_comp(field) / tesla);
305  }
306 
307  return gr;
308 }
309 
main
int main(int argc, char *argv[])
Definition: npdet_fields.cxx:216
settings::step_size
double step_size
Definition: npdet_fields.cxx:46
clipp::value
parameter value(const doc_string &label, Targets &&... tgts)
makes required, blocking, repeatable value parameter; matches any non-empty string
Definition: clipp.h:2094
clipp
primary namespace
Definition: clipp.h:36
clipp::man_page::append_section
man_page & append_section(string title, string content)
Definition: clipp.h:6018
settings::z0
double z0
Definition: npdet_fields.cxx:57
settings::x0
double x0
Definition: npdet_fields.cxx:55
clipp::make_man_page
man_page make_man_page(const group &params, doc_string progname="", const doc_formatting &fmt=doc_formatting{})
generates man sections from command line parameters with sections "synopsis" and "options"
Definition: clipp.h:6086
settings::verbose
bool verbose
Definition: npdet_fields.cxx:66
settings::direction
ROOT::Math::XYZVector direction
Definition: npdet_fields.cxx:45
clipp::option
parameter option(String &&flag, Strings &&... flags)
makes optional, non-blocking exact match parameter
Definition: clipp.h:2078
vis_test.format
format
Definition: beamline/compact/scripts/vis_test.py:12
clipp::parse
parsing_result parse(arg_list args, const group &cli)
parses vector of arg strings and executes actions
Definition: clipp.h:4941
mag_field
void mag_field(dd4hep::Detector &detector)
Cell size example.
Definition: npdet_fields.cxx:158
clipp::detail::action_provider::if_missing
Derived & if_missing(simple_action a)
adds an action that will be called if a required parameter is missing
Definition: clipp.h:1252
settings::infile
std::string infile
Definition: dd_web_display.cxx:78
settings::y1
double y1
Definition: npdet_fields.cxx:59
settings::with_end_point
bool with_end_point
Definition: npdet_fields.cxx:65
cmdline_settings
settings cmdline_settings(int argc, char *argv[])
Definition: npdet_fields.cxx:69
settings::x2
double x2
Definition: npdet_fields.cxx:61
clipp::command
parameter command(String &&flag, Strings &&... flags)
makes required non-blocking exact match parameter
Definition: clipp.h:2048
settings::y2
double y2
Definition: npdet_fields.cxx:62
settings::y0
double y0
Definition: npdet_fields.cxx:56
settings::start_position
ROOT::Math::XYZVector start_position
Definition: npdet_fields.cxx:44
clipp::man_page::prepend_section
man_page & prepend_section(string title, string content)
Definition: clipp.h:6025
clipp::detail::action_provider::set
Derived & set(Target &t)
adds an action that will set the value of 't' from a 'const char*' arg
Definition: clipp.h:1215
build_1D_field_graph
TGraph * build_1D_field_graph(dd4hep::Detector &detector, const settings &s, std::function< double(ROOT::Math::XYZVector)> B_comp)
Definition: npdet_fields.cxx:287
settings::z1
double z1
Definition: npdet_fields.cxx:60
DDG4.Detector
Detector
Definition: DDG4.py:69
dd4hep
Namespace for the AIDA detector description toolkit.
Definition: Geant4Output2Podio.cxx:27
clipp::number
parameter number(const doc_string &label, Targets &&... tgts)
makes required, blocking value parameter; matches any string that represents a number
Definition: clipp.h:2286
settings::field_comps
std::vector< std::string > field_comps
Definition: npdet_fields.cxx:51
settings::axes
std::vector< std::string > axes
Definition: npdet_fields.cxx:53
settings::x1
double x1
Definition: npdet_fields.cxx:58
settings::z2
double z2
Definition: npdet_fields.cxx:63
Generated by   doxygen 1.8.18