Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TTree.h"
#include "TF1.h"
#include "ROOT/TDataFrame.hxx"
#include <vector>
#include <tuple>
#include <algorithm>
#include <iterator>
// DD4hep
// -----
// In .rootlogon.C
// gSystem->Load("libDDDetectors");
// gSystem->Load("libDDG4IO");
// gInterpreter->AddIncludePath("/opt/software/local/include");
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DD4hep/DD4hepUnits.h"
/** Cell size example.
*
*/
void mag_field(double z0 = 0.0 ){
using namespace ROOT::Math;
using namespace dd4hep;
// -------------------------
// Get the DD4hep instance
// Load the compact XML file
dd4hep::Detector& detector = dd4hep::Detector::getInstance();
detector.fromCompact("beamline_test.xml");
double x_min = -100.0;
double x_max = 100.0;
double y_min = -100.0;
double y_max = 100.0;
double z_min = -100.0;
double z_max = 100.0;
int nx_bins = 200;
int ny_bins = 200;
int nz_bins = 200;
double dx = (x_max -x_min)/nx_bins;
double dy = (y_max -y_min)/ny_bins;
double dz = (z_max -z_min)/nz_bins;
double y0 = 20.0*cm;
TH2F* hBx = new TH2F("hBx","; z [cm]; x [cm]; B_{x} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
TH2F* hBz = new TH2F("hBz","; z [cm]; x [cm]; B_{z} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
TH2F* hBr = new TH2F("hBr","; z [cm]; x [cm]; B_{r} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
TH2F* hBphi = new TH2F("hBphi","; z [cm]; x [cm]; B_{#phi} [T]",nz_bins,z_min,z_max,nx_bins,x_min,x_max);
for(int ix = 0; ix< nx_bins; ix++) {
for(int iz = 0; iz< nz_bins; iz++) {
double x = (x_min + ix*dx)*cm;
double z = (z_min + iz*dz)*cm;
XYZVector pos{ x, y0, z };
XYZVector field = detector.field().magneticField(pos);
hBz->Fill(z,x, field.z()/tesla);
hBx->Fill(z,x, field.x()/tesla);
hBr->Fill(z,x, field.rho()/tesla);
hBphi->Fill(z,x, field.phi()/tesla);
}
//for(int iy = 0; iy< ny_bins; iy++) {
// double x = (x_min + ix*dx)*cm;
// double y = (y_min + iy*dy)*cm;
// XYyVector pos{ x, y, z0 };
// XYyVector field = detector.field().magneticField(pos);
// hBz->Fill(z,x, field.z()/tesla);
// hBx->Fill(z,x, field.x()/tesla);
// hBr->Fill(z,x, field.rho()/tesla);
//}
}
auto c = new TCanvas();
c->Divide(2,2);
c->cd(1);
hBz->Draw("colz");
c->cd(2);
hBx->Draw("colz");
c->cd(3);
hBr->Draw("colz");
c->cd(4);
hBphi->Draw("colz");
}