Project

General

Profile

Bug #1951 » check_geometry.py

Rogers, Chris, 22 January 2018 13:12

 
1
import Configuration
2
import maus_cpp.globals
3
import maus_cpp.material
4
import xboa.common
5
import ROOT
6
import time
7

    
8
VERBOSE = False
9

    
10
def initialise_maus():
11
    configuration = Configuration.Configuration().\
12
                                          getConfigJSON(command_line_args=True)
13
    maus_cpp.globals.birth(configuration)
14

    
15
MATERIAL_LIST = []
16
def material_to_colour(material):
17
    global MATERIAL_LIST
18
    if material[0:3] == "G4_":
19
        material = material[3:]
20
    if material not in MATERIAL_LIST:
21
        MATERIAL_LIST.append(material)
22
    if material in ("Galactic", "AIR", "He"):
23
        return None
24
    if material in ("Fe"): # "kill volumes"
25
        return 1
26
    if material in ("MYLAR", "POLYSTYRENE", "RenCast6400", "NYLON-6-6", "POLYCARBONATE", "POLYVINYL_TOLUENE", "POLYURETHANE", "G10", "TUFNOL"):
27
        return 8
28
    if material in ("Zn", "Cu", "W", "Al", "ALUMINUM", "TUNGSTEN", "BRASS", "STEEL", "IRON"):
29
        return 2
30
    if material in ("lH2", "MICE_LITHIUM_HYDRIDE", "LITHIUM_HYDRIDE"):
31
        return 4
32
    print "UNRECOGNISED MATERIAL", material
33
    return 1
34

    
35
def get_materials(radius, z_start, z_end, z_step):
36
    x = radius
37
    material = None
38
    material_start = []
39
    n_steps = int((z_end-z_start)/z_step)
40
    for i in range(n_steps):
41
        z = z_step*i+z_start
42
        maus_cpp.material.set_position(x, 0., z)
43
        material_data = maus_cpp.material.get_material_data()
44
        if VERBOSE == True:
45
            print x, z, material
46
        new_material = material_data['name']
47
        if new_material != material:
48
            material = new_material
49
            material_start.append({"x":x, "z":z, "material":material})
50
    return material_start
51

    
52
ROOT_GRAPHS = []
53
def plot_materials(r_start, r_end, r_step, z_start, z_end, z_step):
54
    global ROOT_GRAPHS
55
    canvas = xboa.common.make_root_canvas("materials")
56
    canvas.SetWindowSize(1900, 1000)
57
    n_steps = int((r_end-r_start)/r_step)
58
    hist = ROOT.TH2D("materials", ";z [mm]; x [mm]", 1000, z_start, z_end, 1000, r_start, r_end)
59
    hist.SetStats(False)
60
    hist.Draw()
61
    ROOT_GRAPHS.append(hist)
62
    for i in range(n_steps):
63
        r = r_step*i+r_start
64
        materials = get_materials(r, z_start,z_end, z_step)
65
        print "At radius", r, "found", len(materials), "materials using", len(ROOT_GRAPHS), "root objects"
66
        for i, material in enumerate(materials):
67
            colour = material_to_colour(material["material"])
68
            if colour == None:
69
                continue
70
            z_min = material["z"]
71
            radius = material["x"]
72
            if i+1 >= len(materials):
73
                z_max = z_end+1
74
            else:
75
                z_max = materials[i+1]["z"]
76
            if i == 0:
77
                z_min -= 1
78
            graph = ROOT.TGraph(2)
79
            graph.SetPoint(0, z_min, radius)
80
            graph.SetPoint(1, z_max, radius)
81
            graph.SetLineColor(colour)
82
            graph.SetMarkerColor(colour)
83
            graph.SetMarkerStyle(6)
84
            graph.SetLineWidth(2)
85
            graph.Draw("plsame")
86
            ROOT_GRAPHS.append(graph)
87
            if i % 10 == 0:
88
                canvas.Update()
89

    
90
    canvas.Update()
91
    for format in "png", "eps", "root":
92
        canvas.Print("plots/materials."+format)
93

    
94
def main():
95
    initialise_maus()
96
    old_time = time.time()
97
    plot_materials(0.0, 150.5, 1., 20200., 20500., 0.1)
98
    print "Plotting took", time.time() - old_time, "seconds"
99
    print "Found the following materials", MATERIAL_LIST
100

    
101
if __name__ == "__main__":
102
    main()
103
    raw_input()
(3-3/5)