FairShip
Loading...
Searching...
No Matches
checkMagFields.py
Go to the documentation of this file.
1# SPDX-License-Identifier: LGPL-3.0-or-later
2# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
3
4import sys
5
6import ROOT
7import rootUtils as ut
8import shipunit as u
9
10
11def run():
12 fGeo = ROOT.gGeoManager
13 run = sys.modules["__main__"].run
14 if hasattr(sys.modules["__main__"], "h"):
15 h = sys.modules["__main__"].h
16 else:
17 h = {}
18 grid = 120, 100, 1500
19 xmin, ymin, zmin = -4 * u.m, -5 * u.m, -100 * u.m
20 xmax, ymax, zmax = 4 * u.m, 5 * u.m, 50 * u.m
21 dx, dy, dz = (xmax - xmin) / grid[0], (ymax - ymin) / grid[1], (zmax - zmin) / grid[2]
22 ut.bookHist(
23 h, "Bx-", "Bx- component ;z [cm];x [cm];y [cm]", grid[2], zmin, zmax, grid[0], xmin, xmax, grid[1], ymin, ymax
24 )
25 ut.bookHist(
26 h, "By-", "By- component ;z [cm];x [cm];y [cm]", grid[2], zmin, zmax, grid[0], xmin, xmax, grid[1], ymin, ymax
27 )
28 ut.bookHist(
29 h, "Bz-", "Bz- component ;z [cm];x [cm];y [cm]", grid[2], zmin, zmax, grid[0], xmin, xmax, grid[1], ymin, ymax
30 )
31 ut.bookHist(
32 h, "Bx+", "Bx+ component ;z [cm];x [cm];y [cm]", grid[2], zmin, zmax, grid[0], xmin, xmax, grid[1], ymin, ymax
33 )
34 ut.bookHist(
35 h, "By+", "By+ component ;z [cm];x [cm];y [cm]", grid[2], zmin, zmax, grid[0], xmin, xmax, grid[1], ymin, ymax
36 )
37 ut.bookHist(
38 h, "Bz+", "Bz+ component ;z [cm];x [cm];y [cm]", grid[2], zmin, zmax, grid[0], xmin, xmax, grid[1], ymin, ymax
39 )
40 h["Bx-"].SetMarkerColor(ROOT.kGreen - 3)
41 h["Bx+"].SetMarkerColor(ROOT.kGreen + 3)
42 h["By-"].SetMarkerColor(ROOT.kBlue - 3)
43 h["By+"].SetMarkerColor(ROOT.kBlue + 3)
44 h["Bz-"].SetMarkerColor(ROOT.kCyan - 2)
45 h["Bz+"].SetMarkerColor(ROOT.kCyan + 2)
46 for ix in range(grid[0]):
47 for iy in range(grid[1]):
48 for iz in range(grid[2]):
49 x, y, z = xmin + ix * dx, ymin + iy * dy, zmin + iz * dz
50 n = fGeo.FindNode(x, y, z)
51 f = n.GetVolume().GetField()
52 if f:
53 if f.GetFieldValue()[0] < 0:
54 h["Bx-"].Fill(z, x, y, -f.GetFieldValue()[0] / u.tesla)
55 if f.GetFieldValue()[0] > 0:
56 h["Bx+"].Fill(z, x, y, f.GetFieldValue()[0] / u.tesla)
57 if f.GetFieldValue()[1] < 0:
58 h["By-"].Fill(z, x, y, -f.GetFieldValue()[1] / u.tesla)
59 if f.GetFieldValue()[1] > 0:
60 h["By+"].Fill(z, x, y, f.GetFieldValue()[1] / u.tesla)
61 if f.GetFieldValue()[2] < 0:
62 h["Bz-"].Fill(z, x, y, -f.GetFieldValue()[2] / u.tesla)
63 if f.GetFieldValue()[2] > 0:
64 h["Bz+"].Fill(z, x, y, f.GetFieldValue()[2] / u.tesla)
65 f = run.GetField()
66 if f.GetBx(x, y, z) < 0:
67 h["Bx-"].Fill(z, x, y, -f.GetBx(x, y, z) / u.tesla)
68 if f.GetBx(x, y, z) > 0:
69 h["Bx+"].Fill(z, x, y, f.GetBx(x, y, z) / u.tesla)
70 if f.GetBy(x, y, z) < 0:
71 h["By-"].Fill(z, x, y, -f.GetBy(x, y, z) / u.tesla)
72 if f.GetBy(x, y, z) > 0:
73 h["By+"].Fill(z, x, y, f.GetBy(x, y, z) / u.tesla)
74 for x in h:
75 hi = h[x]
76 if hi.ClassName() == "TH3F":
77 h[x + "_xz"] = h[x].Project3D("xy")
78 h[x + "_xz"].SetTitle(hi.GetTitle() + " top view")
79 h[x + "_yz"] = h[x].Project3D("xz")
80 h[x + "_yz"].SetTitle(hi.GetTitle() + " side view")
81 for x in h:
82 h[x].SetStats(0)
83 h[x].SetMarkerSize(3)
84 txt = {"y": [" Up", " Down"], "x": [" Right", " Left"]}
85 for pol in ["y", "x"]:
86 for p in ["_xz", "_yz"]:
87 cn = "c" + pol + p
88 ut.bookCanvas(h, key=cn, title="field check", nx=1600, ny=1200, cx=1, cy=1)
89 h[cn].cd(1)
90 h["B" + pol + "+" + p].Draw()
91 h["B" + pol + "-" + p].Draw("same")
92 h["B" + pol + "L" + p] = ROOT.TLegend(0.79, 0.72, 0.91, 0.87)
93 h["B" + pol + "L" + p].AddEntry(h["B" + pol + "+"], "B" + pol + txt[pol][0], "PM")
94 h["B" + pol + "L" + p].AddEntry(h["B" + pol + "-"], "B" + pol + txt[pol][1], "PM")
95 h["B" + pol + "L" + p].Draw()
96 h[cn].Update()
97 h[cn].Print("FieldB" + pol + "Proj" + p + ".png")