Belle II Software  release-05-01-25
ClusterEfficiency.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 """
5 <header>
6  <contact>ritter</contact>
7  <description>
8  Plot the efficiency to find a cluster for a given truehit for all layers of
9  the VXD.
10  </description>
11 </header>
12 """
13 import sys
14 import math
15 from basf2 import *
16 from ROOT import PyConfig
17 PyConfig.IgnoreCommandLineOptions = 1
18 from svd import *
19 from pxd import *
20 import ROOT
21 from ROOT import Belle2
22 set_log_level(LogLevel.ERROR)
23 # set some random seed
24 set_random_seed(10346)
25 
26 # momenta to generate the plots for
27 momenta = [3.0]
28 # theta parameters
29 theta_params = [90, 0.1]
30 
31 
32 class ClusterEfficiency(Module):
33  """
34  Plot Efficiency to find a U and V cluster for each given truehit.
35  """
36  def initialize(self):
37  """
38  Create ROOT TProfiles for all layers and momenta.
39  """
40 
41  # let's create a root file to store all profiles
42 
43  self.rfile = ROOT.TFile("ClusterEfficiency.root", "RECREATE")
44  self.rfile.cd()
45 
46  self.profiles = []
47 
48  self.eff = {}
49 
50  def profile(name, title, text, contact):
51  """small helper function to create a phi profile and set the names
52  and descriptions"""
53  prof = ROOT.TProfile(name, title, 60, -180, 180)
54  prof.GetListOfFunctions().Add(ROOT.TNamed("Description", text))
55  prof.GetListOfFunctions().Add(ROOT.TNamed("Contact", contact))
56  prof.GetListOfFunctions().Add(ROOT.TNamed("Check", "Should be close to 1 everywhere"))
57  # make a list of all profiles
58  self.profiles.append(prof)
59  return prof
60 
61  # name, title and decription templates
62  prof_name = "ClusterEfficiency_layer{layer}_{p:.1f}GeV"
63  prof_title = "Cluster Efficiency in #phi, layer={layer}, "\
64  + "p={p:.1f} GeV;#phi in degrees;efficiency"
65  prof_text = "Efficiency to find a U and V cluster for any given truehit "\
66  + "in layer {layer} when simulation muons with p={p:.1f} GeV uniformly "\
67  + "in phi and theta={theta[0]}+-{theta[1]} degree. phi is the angle "\
68  + "of the generated particle, not the truehit position."
69  # contact person for PXD (layer<3 = True) or SVD (layer<3 = False)
70  prof_contact = {
71  True: "Benjamin Schwenker <Benjamin.Schwenker@phys.uni-goettingen.de>",
72  False: "Andrzej Bozek <bozek@belle2.ifj.edu.pl>",
73  }
74 
75  # create all profiles
76  for layer in range(1, 7):
77  self.eff[layer] = {}
78  for p in momenta:
79  name = prof_name.format(p=p, layer=layer)
80  title = prof_title.format(p=p, layer=layer)
81  text = prof_text.format(p=p, layer=layer, theta=theta_params)
82  self.eff[layer][p] = profile(name, title, text, prof_contact[layer < 3])
83 
84  def terminate(self):
85  """
86  Format all profiles and write the ROOT file.
87  """
88  # determine min value in all profiles to set all profiles to a common
89  # range. We always want to show at least the range from 0.9
90  minVal = 0.9
91  for prof in self.profiles:
92  minBin = prof.GetMinimumBin()
93  # minimum value is bin content - error and a 5% margin
94  minVal = min(minVal, (prof.GetBinContent(minBin) - prof.GetBinError(minBin)) * 0.95)
95 
96  # no let's set the range of all profiles to be equal
97  for prof in self.profiles:
98  prof.SetMinimum(max(0, minVal))
99  prof.SetMaximum(1.02)
100  prof.SetStats(False)
101 
102  self.rfile.Write()
103  self.rfile.Close()
104 
105  def fill_truehits(self, phi, p, truehits):
106  """
107  Loop over all truehits for a track with the given angle phi and momentum
108  p and update the efficiency profiles.
109  """
110  # iterate over all truehits and look for the clusters
111  for i, truehit in enumerate(truehits):
112  # we ignore all truehits which have an negative weight
113  if truehits.weight(i) < 0:
114  continue
115 
116  # get layer number and a list of all clusters related to the truehit
117  layer = truehit.getSensorID().getLayerNumber()
118  if isinstance(truehit, Belle2.SVDTrueHit):
119  is_pxd = False
120  clusters = truehit.getRelationsFrom("SVDClusters")
121  else:
122  is_pxd = True
123  clusters = truehit.getRelationsFrom("PXDClusters")
124 
125  # now check if we find a cluster
126  has_cluster = {False: 0, True: 0}
127  for j, cls in enumerate(clusters):
128  # we ignore all clusters where less then 100 electrons come from
129  # our truehit
130  if clusters.weight(j) < 100:
131  continue
132 
133  if is_pxd:
134  # for pxt we always have both directions so set efficiency
135  # for this TrueHit to 1
136  has_cluster[True] = 1
137  has_cluster[False] = 1
138  else:
139  # for SVD we remember that we set the efficiency for the
140  # direction of the cluster
141  has_cluster[cls.isUCluster()] = 1
142 
143  # and we fill the profile
144  self.eff[layer][p].Fill(phi, has_cluster[True] & has_cluster[False])
145 
146  def event(self):
147  """
148  Update the efficienies by iterating over all primary particles
149  """
150  mcparticles = Belle2.PyStoreArray("MCParticles")
151  for mcp in mcparticles:
152  # we only look at primary particles
153  if not mcp.hasStatus(1):
154  continue
155 
156  # let's get the momentum and determine which of the ones we wanted
157  # to generate it is
158  p = mcp.getMomentum()
159  p_gen = None
160  for i in momenta:
161  if abs(p.Mag() - i) < 0.05:
162  p_gen = i
163  break
164 
165  # meh, something strange with the momentum, ignore this one
166  if p_gen is None:
167  B2WARNING("Strange particle momentum: %f, expected one of %s" %
168  (p.Mag(), ", ".join(str() for p in momenta)))
169  continue
170 
171  # and check all truehits
172  pxdtruehits = mcp.getRelationsTo("PXDTrueHits")
173  self.fill_truehits(math.degrees(p.Phi()), p_gen, pxdtruehits)
174  svdtruehits = mcp.getRelationsTo("SVDTrueHits")
175  self.fill_truehits(math.degrees(p.Phi()), p_gen, svdtruehits)
176 
177 # Now let's create a path to simulate our events. We need a bit of statistics but
178 # that's not too bad since we only simulate single muons
179 main = create_path()
180 main.add_module("EventInfoSetter", evtNumList=[10000])
181 main.add_module("Gearbox")
182 # we only need the vxd for this
183 main.add_module("Geometry", components=['MagneticFieldConstant4LimitedRSVD',
184  'BeamPipe', 'PXD', 'SVD'])
185 particlegun = main.add_module("ParticleGun")
186 particlegun.param({
187  "nTracks": 1,
188  "pdgCodes": [13, -13],
189  # generate discrete momenta with equal weights
190  "momentumGeneration": 'discrete',
191  "momentumParams": momenta + [1]*len(momenta),
192  "thetaGeneration": 'normal',
193  "thetaParams": theta_params,
194 })
195 main.add_module("FullSim")
196 add_pxd_simulation(main)
197 add_svd_simulation(main)
198 add_pxd_reconstruction(main)
199 add_svd_reconstruction(main)
200 
201 clusterefficiency = ClusterEfficiency()
202 main.add_module(clusterefficiency)
203 main.add_module("Progress")
204 
205 process(main)
206 print(statistics)
ClusterEfficiency.ClusterEfficiency.rfile
rfile
Output file to store all plots.
Definition: ClusterEfficiency.py:43
ClusterEfficiency.ClusterEfficiency.fill_truehits
def fill_truehits(self, phi, p, truehits)
Definition: ClusterEfficiency.py:105
ClusterEfficiency.ClusterEfficiency.initialize
def initialize(self)
Definition: ClusterEfficiency.py:36
Belle2::SVDTrueHit
Class SVDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: SVDTrueHit.h:35
ClusterEfficiency.ClusterEfficiency.eff
eff
layer/momentum hierachy of all profiles
Definition: ClusterEfficiency.py:48
ClusterEfficiency.ClusterEfficiency.event
def event(self)
Definition: ClusterEfficiency.py:146
ClusterEfficiency.ClusterEfficiency
Definition: ClusterEfficiency.py:32
ClusterEfficiency.ClusterEfficiency.profiles
profiles
flat list of all profiles for easy access
Definition: ClusterEfficiency.py:46
Belle2::PyStoreArray
a (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:58
ClusterEfficiency.ClusterEfficiency.terminate
def terminate(self)
Definition: ClusterEfficiency.py:84