DD4hep  1.30.0
Detector Description Toolkit for High Energy Physics
DrawBField.cpp
Go to the documentation of this file.
1 //==========================================================================
2 // AIDA Detector description implementation
3 //--------------------------------------------------------------------------
4 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
5 // All rights reserved.
6 //
7 // For the licensing terms see $DD4hepINSTALL/LICENSE.
8 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
9 //
10 // Author : M.Frank
11 //
12 //==========================================================================
13 
14 // Framework includes
15 #include <DD4hep/Printout.h>
16 #include <DD4hep/FieldTypes.h>
17 #include <DD4hep/MatrixHelpers.h>
18 #include <DD4hep/DetFactoryHelper.h>
19 #include <cmath>
20 
21 using namespace dd4hep;
22 
23 #include <TH2.h>
24 #include <TPad.h>
25 #include <TRint.h>
26 #include <TAxis.h>
27 #include <TStyle.h>
28 #include <TArrow.h>
29 #include <TCanvas.h>
30 #include <THistPainter.h>
31 
32 #include "Hoption.h"
33 #include "Hparam.h"
34 
35 extern Hoption_t Hoption;
36 extern Hparam_t Hparam;
37 
38 
39 namespace {
40  void help_bfield(int argc, char** argv) {
41  std::cout <<
42  "Usage: DD4hep_DrawBField -arg [-arg] \n"
43  " The plugin implicitly assumes that the geometry is already loaded. \n"
44  " -area <string> Define area for which the B-field will be drawn\n"
45  " -nx: <number> Number of bins in X direction \n"
46  " -ny: <number> Number of bins in Y direction \n"
47  " -dx: <range> Definition of the range in X \n"
48  " -dy: <range> Definition of the range in Y \n"
49  " -dz: <range> Definition of the range in Z \n"
50  " Range definition as tuple: min,max \n"
51  " min: lower boundary \n"
52  " max: upper boundary \n"
53  "\tArguments given: " << arguments(argc,argv) << std::endl << std::flush;
54  ::exit(EINVAL);
55  }
56 
57  struct range_t {
58  double rmin=0e0, rmax=0e0;
59  };
60  struct field_t {
61  Position position;
62  Direction bfield;
63  };
64 
65  range_t get_range(const std::string& s, int argc, char** argv) {
66  int ival=0;
67  range_t range;
68  std::string val;
69  for(std::size_t i=0; i<s.length()+1; ++i) {
70  char c = s.c_str()[i];
71  if ( c == ',' || i == s.length() ) {
72  switch(ival) {
73  case 0:
74  range.rmin = _toDouble(val);
75  break;
76  case 1:
77  range.rmax = _toDouble(val);
78  break;
79  default:
80  except("","+++ Too many value for range descriptor provided: %s", s.c_str());
81  break;
82  }
83  val = "";
84  ++ival;
85  continue;
86  }
87  val += c;
88  }
89  if ( ival != 2 ) {
90  help_bfield(argc, argv);
91  }
92  return range;
93  }
94 
95  class MyHistPainter : public THistPainter {
96  public:
97  using THistPainter::THistPainter;
98  TH2* paintArrows(TH2* histo, const std::vector<field_t>& points) {
99  fYaxis->SetTitle("Y coordinate [cm]");
100  fXaxis->SetTitle("X coordinate [cm]");
101  for(const auto& p : points) {
102  double strength2 = p.bfield.mag2();
103  if ( strength2 > 1e-8 ) {
104  double strength = sqrt(p.bfield.X()*p.bfield.X()+p.bfield.Y()*p.bfield.Y());
105  int ix = fXaxis->FindBin(p.position.X());
106  int iy = fYaxis->FindBin(p.position.Y());
107  double xlo = fXaxis->GetBinLowEdge(ix);
108  double xhi = fXaxis->GetBinLowEdge(ix+1);
109  double ylo = fYaxis->GetBinLowEdge(iy);
110  double yhi = fYaxis->GetBinLowEdge(iy+1);
111  auto norm_field = p.bfield * ((xhi-xlo) / strength);
112  double x1 = xlo + (xhi-xlo)/2.0;
113  double x2 = x1 + norm_field.X();
114  double y1 = ylo + (yhi-ylo)/2.0;
115  double y2 = y1 + norm_field.Y();
116  auto* arrow = new TArrow(x1, y1, x2, y2, 0.015, "|>");
117  arrow->SetAngle(25);
118  arrow->SetFillColor(kRed);
119  arrow->SetLineColor(kRed);
120  arrow->Draw();
121  histo->Fill(p.position.X(), p.position.X(), strength);
122  ::fprintf(stdout, "Arrow %+15.8e %+15.8e %+15.8e %+15.8e %+15.8e %+15.8e\n",
123  x1, y1, x2, y2, x2-x1, y2-y1);
124  //delete arrow;
125  }
126  }
127  PaintPalette();
128  return histo;
129  }
130  };
131 
133 
145  static long draw_bfield(Detector& description, int argc, char** argv) {
146  std::vector<std::string> srange_x, srange_y, srange_z;
147  std::size_t nbin_x = 100, nbin_y = 100, nbin_z = 1;
148  std::string output;
149  double z_value = 0e0;
150  bool draw = false;
151 
152  printout(WARNING,"DrawBField","This is a TEST. It does not function properly!");
153 
154  for(int i = 0; i < argc && argv[i]; ++i) {
155  if ( 0 == ::strncmp("-rx",argv[i],4) )
156  srange_x.push_back(argv[++i]);
157  else if ( 0 == ::strncmp("-ry",argv[i],3) )
158  srange_y.push_back(argv[++i]);
159  else if ( 0 == ::strncmp("-rz",argv[i],3) )
160  srange_z.push_back(argv[++i]);
161  else if ( 0 == ::strncmp("-nx",argv[i],3) )
162  nbin_x = _toULong(argv[++i]);
163  else if ( 0 == ::strncmp("-ny",argv[i],3) )
164  nbin_y = _toULong(argv[++i]);
165  else if ( 0 == ::strncmp("-nz",argv[i],3) )
166  nbin_z = _toULong(argv[++i]);
167  else if ( 0 == ::strncmp("-z",argv[i],2) )
168  z_value = _toDouble(argv[++i]);
169  else if ( 0 == ::strncmp("-draw",argv[i],4) )
170  draw = true;
171  else if ( 0 == ::strncmp("-output",argv[i],4) )
172  output = argv[++i];
173  else if ( 0 == ::strncmp("-help",argv[i],2) )
174  help_bfield(argc, argv);
175  }
176  if ( srange_x.empty() || srange_y.empty() ) {
177  help_bfield(argc, argv);
178  }
179  std::vector<range_t> range_x, range_y, range_z;
180  for(const auto& r : srange_x) {
181  range_t rg = get_range(r, argc, argv);
182  range_x.push_back(rg);
183  }
184  for(const auto& r : srange_y) {
185  range_t rg = get_range(r, argc, argv);
186  range_y.push_back(rg);
187  }
188  for(const auto& r : srange_z) {
189  range_t rg = get_range(r, argc, argv);
190  range_z.push_back(rg);
191  }
192 
193  double ma = std::numeric_limits<double>::max();
194  range_t envelope_x { ma, -ma };
195  range_t envelope_y { ma, -ma };
196  range_t envelope_z { ma, -ma };
197  for( const auto& range : range_x ) {
198  envelope_x.rmin = std::min(range.rmin, envelope_x.rmin);
199  envelope_x.rmax = std::max(range.rmax, envelope_x.rmax);
200  }
201  for( const auto& range : range_y ) {
202  envelope_y.rmin = std::min(range.rmin, envelope_y.rmin);
203  envelope_y.rmax = std::max(range.rmax, envelope_y.rmax);
204  }
205  for( const auto& range : range_z ) {
206  envelope_z.rmin = std::min(range.rmin, envelope_z.rmin);
207  envelope_z.rmax = std::max(range.rmax, envelope_z.rmax);
208  }
209 
210  double dx = (envelope_x.rmax - envelope_x.rmin) / double(nbin_x);
211  double dy = (envelope_y.rmax - envelope_y.rmin) / double(nbin_y);
212  double dz = nbin_z == 1 ? 0e0 : (envelope_z.rmax - envelope_z.rmin) / double(nbin_z);
213  printf("Range(x) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
214  nbin_x, envelope_x.rmin/cm, envelope_x.rmax/cm, dx/cm);
215  printf("Range(y) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
216  nbin_y, envelope_y.rmin/cm, envelope_y.rmax/cm, dy/cm);
217  if ( nbin_z > 1 ) printf("Range(z) min:%4ld bins %+15.8e cm max:%+15.8e cm dx:%+15.8e cm\n",
218  nbin_z, envelope_z.rmin/cm, envelope_z.rmax/cm, dz/cm);
219 
220  FILE* out_file = ::fopen(output.empty() ? "/dev/null" : output.c_str(), "w");
221  ::fprintf(out_file,"#######################################################################################################\n");
222  ::fprintf(out_file," x[cm] y[cm] z[cm] Bx[Tesla] By[Tesla] Bz[Tesla] \n");
223  std::vector<field_t> field_values;
224  for( std::size_t i = 0; i < nbin_x; ++i ) {
225  float x = envelope_x.rmin + double(i)*dx + dx/2e0;
226  for( std::size_t j = 0; j < nbin_y; ++j ) {
227  float y = envelope_y.rmin + double(j)*dy + dy/2e0;
228  for( std::size_t k = 0; k < nbin_z; ++k ) {
229  float z = nbin_z == 1 ? z_value : envelope_z.rmin + double(k)*dz + dz/2e0;
230  field_t value;
231  value.position = { x, y, z };
232  value.bfield = { 0e0, 0e0, 0e0 };
233  value.bfield = description.field().magneticField(value.position);
234  ::fprintf(out_file, " %+15.8e %+15.8e %+15.8e %+15.8e %+15.8e %+15.8e\n",
235  value.position.X()/cm, value.position.Y()/cm, value.position.Z()/cm,
236  value.bfield.X()/dd4hep::tesla, value.bfield.Y()/dd4hep::tesla, value.bfield.Z()/dd4hep::tesla);
237  field_values.emplace_back(value);
238  }
239  }
240  }
241  ::fclose(out_file);
242  if ( draw ) {
243  if ( 0 == gApplication ) {
244  std::pair<int, char**> a(argc,argv);
245  gApplication = new TRint("DD4hepRint", &a.first, a.second);
246  printout(INFO,"DD4hepRint","++ Created ROOT interpreter instance for DD4hepUI.");
247  }
248  auto* histo = new TH2F("Bfield", "B-Field strength in Tesla",
249  nbin_x, envelope_x.rmin, envelope_x.rmax,
250  nbin_y, envelope_y.rmin, envelope_y.rmax);
251  MyHistPainter paint;
252  paint.SetHistogram(histo);
253  TCanvas* c1 = new TCanvas("B-Field");
254  c1->SetWindowSize(1000,1000);
255  //paint.Paint();
256  histo->SetStats(kFALSE);
257  //histo->Draw();
258  paint.paintArrows(histo, field_values);
259  histo->Draw("COLZ SAME");
260  gPad->SetGridx();
261  gPad->SetGridy();
262  if ( !gApplication->IsRunning() ) {
263  gApplication->Run();
264  }
265  }
266  return 1;
267  }
268 }
269 DECLARE_APPLY(DD4hep_DrawBField,draw_bfield)
MatrixHelpers.h
Hoption
Hoption_t Hoption
DECLARE_APPLY
#define DECLARE_APPLY(name, func)
Definition: Factories.h:281
dd4hep::Direction
Position Direction
Definition: Fields.h:28
dd4hep::_toULong
unsigned long _toULong(const std::string &value)
String conversions: string to long integer value.
Definition: Handle.cpp:102
Hparam
Hparam_t Hparam
dd4hep::Detector::field
virtual OverlayedField field() const =0
Return handle to the combined electromagentic field description.
dd4hep::Position
ROOT::Math::XYZVector Position
Definition: Objects.h:81
DetFactoryHelper.h
dd4hep::_toDouble
double _toDouble(const std::string &value)
String conversions: string to double value.
Definition: Handle.cpp:116
dd4hep
Namespace for the AIDA detector description toolkit.
Definition: AlignmentsCalib.h:28
dd4hep::Detector
The main interface to the dd4hep detector description package.
Definition: Detector.h:90
dd4hep::OverlayedField::magneticField
void magneticField(const Position &pos, double *field) const
Returns the 3 magnetic field components (x, y, z).
Definition: Fields.cpp:140
FieldTypes.h
Printout.h
Generated on Thu Nov 7 2024 13:31:56 for DD4hep by   doxygen 1.8.18