Belle II Software  release-05-01-25
MaterialRay.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 from basf2 import logging, LogLevel, create_path, process
5 import ROOT as root
6 import numpy as np
7 import matplotlib as mpl
8 mpl.use("Agg")
9 from matplotlib import pyplot as pl
10 from root_numpy import hist2array
11 
12 # Don't show all the info messages
13 logging.log_level = LogLevel.ERROR
14 
15 # create path
16 main = create_path()
17 # We need to process one event so we need the EventInfoSetter
18 main.add_module("EventInfoSetter", evtNumList=[1])
19 # We need the geometry parameters
20 main.add_module("Gearbox")
21 # as well as the geometry. We assign regions to each creator which allows to
22 # split the material budget by region instead of material.
23 geometry = main.add_module("Geometry", logLevel=LogLevel.INFO, assignRegions=True)
24 # MaterialScan uses the Geant4 setup which is created by the FullSim module so
25 # we need this as well
26 main.add_module('FullSim')
27 # And finally the MaterialScan module
28 materialscan = main.add_module("MaterialScan", logLevel=LogLevel.INFO)
29 materialscan.param({
30  # Filename for output File
31  'Filename': 'MaterialRay.root',
32  'spherical': False,
33  'planar': False,
34 
35  'ray': True,
36  'ray.theta': 140,
37  'ray.phi': 0,
38  # Alternatively one can set the direction of the ray directly
39  # 'ray.direction': [1, 0, 0],
40  # Max depth for the scan in cm (x axis limit on the histogram). 0 = no limit
41  'ray.maxDepth': 0,
42  # Bin width (roughly) for the output histogram in cm
43  'ray.sampleDepth': 0.1,
44  # Opening angle for the rays: If <=0 only one ray is shot. Otherwise all
45  # rays will be distributed randomly in a cone with this opening angle in
46  # degrees, distributed uniformly in the area (flat in cos(theta)
47  'ray.opening': 0.05,
48  # Number of rays to shoot if opening>0
49  'ray.count': 1000,
50 })
51 
52 # Do the MaterialScan, can take some time depending on the number of steps
53 process(main)
54 
55 # Done, start plotting
56 rmaterial_file = root.TFile("MaterialRay.root")
57 
58 
59 def plot_hist(region, **argk):
60  """Get the histogram for a given region from the file and plot it. Return the histogram data and edges"""
61  h = rmaterial_file.Get("Ray/%s_x0" % region)
62  if not h:
63  return None
64  data, edges = hist2array(h, return_edges=True)
65  # Now we have to play around a bit: pl.plot(edges[:-1], data, drawstyle="steps-post") does work but
66  # the last bin does not get a top line so we need to add a 0 to the
67  # histogram to go back down to where we were
68  data = np.append(data, 0)
69  # now plot
70  pl.plot(edges[0], data, drawstyle="steps-post", **argk)
71  return data, edges[0]
72 
73 
74 # Get the all materials plot:
75 all_data, mat_edges = plot_hist("All_Regions", label="All regions", c="k")
76 
77 # Now plot the different regions separately
78 for region in ["BeamPipe", "PXD", "SVD", "CDC", "ARICH", "TOP", "ECL", "EKLM", "BKLM", "STR", "COIL"]:
79  plot_hist(region, label=region)
80 
81 # Log scale anyone?
82 # pl.yscale("log")
83 pl.ylim(ymin=0)
84 pl.xlim(xmin=mat_edges[0], xmax=mat_edges[-1]+1)
85 pl.xlabel("Flight length [cm]")
86 pl.ylabel("Radiation length [$X_0 / %.3g$ cm]" % (mat_edges[1] - mat_edges[0]))
87 pl.legend(loc="best")
88 pl.tight_layout()
89 pl.savefig("MaterialRay.pdf")