Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • fmetzner/AnalysisTools
1 result
Select Git revision
Show changes
Commits on Source (2)
Showing
with 619 additions and 13 deletions
%% Cell type:code id:5eb05a44-f720-45b4-9362-be108931b133 tags:
``` python
import numpy as np
from numpy.linalg import eig
import uncertainties as unct
```
%% Cell type:code id:c5488e67-0c5a-4f24-9dba-e260c86928cd tags:
``` python
Cov = np.array([[1,0.25,0],[0.25,1,0],[0,0,1]])
```
%% Cell type:code id:d3a7c8c4-6719-4289-963d-0572f3bdabe5 tags:
``` python
print(Cov)
```
%% Output
[[1. 0.25 0. ]
[0.25 1. 0. ]
[0. 0. 1. ]]
%% Cell type:code id:d87f3791-3fdb-41f2-a54b-389ba6c9a991 tags:
``` python
w, v=eig(Cov)
```
%% Cell type:code id:cdd65172-bdbc-4c69-99ef-16069b06f82c tags:
``` python
print(w)
print(v)
```
%% Output
[1.25 0.75 1. ]
[[ 0.70710678 -0.70710678 0. ]
[ 0.70710678 0.70710678 0. ]
[ 0. 0. 1. ]]
%% Cell type:code id:a5cbe35a-b6c1-4d15-a8d8-ecd2d03430e2 tags:
``` python
for i in range(0,len(w)):
ev_wrong = np.sqrt(w[i]) * v[i]
print(ev_wrong)
for i in range(0,len(w)):
ev_correct = np.sqrt(w[i]) * v[:,i]
print(ev_correct)
```
%% Output
[ 0.79056942 -0.79056942 0. ]
[0.61237244 0.61237244 0. ]
[0. 0. 1.]
[0.79056942 0.79056942 0. ]
[-0.61237244 0.61237244 0. ]
[0. 0. 1.]
%% Cell type:code id:328c154e-9b9b-47d7-8484-286eae1dc53c tags:
``` python
for i in range(0,len(w)):
variation = np.zeros(len(w))
variation[i] = np.sqrt(w[i])
print(variation @ v)
for i in range(0,len(w)):
variation = np.zeros(len(w))
variation[i] = np.sqrt(w[i])
print(v @ variation)
```
%% Output
[ 0.79056942 -0.79056942 0. ]
[0.61237244 0.61237244 0. ]
[0. 0. 1.]
[0.79056942 0.79056942 0. ]
[-0.61237244 0.61237244 0. ]
[0. 0. 1.]
%% Cell type:code id:7cf38427-2395-4f40-b82e-f176a7e79590 tags:
``` python
(a,b,c) = unct.correlated_values([1,1,1], Cov)
s = a + b + c
print(s.std_dev())
```
%% Output
1.8708286933869707
<ipython-input-85-c53222b3ef90>:5: UserWarning: Obsolete: the std_dev attribute should not be called anymore: use .std_dev instead of .std_dev(). Code can be automatically updated with python -m uncertainties.1to2 -w ProgramDirectory.
print(s.std_dev())
%% Cell type:code id:1c1d71b6-8146-4c29-bd44-4cc846e3a007 tags:
``` python
values = np.array([1,1,1])
np.sum(values)
```
%% Output
3
%% Cell type:code id:7e3cde6f-e255-4c96-b626-e4b2283c28d7 tags:
``` python
err = 0
for i in range(0,len(w)):
var = values + np.sqrt(w[i]) * v[:,i]
df = np.sum(values) - np.sum(var)
err = np.sqrt(err**2 + df**2)
print(err)
```
%% Output
1.8708286933869702
%% Cell type:code id:785c6727-a9b7-40be-a09a-f2243f297162 tags:
``` python
err = 0
for i in range(0,len(w)):
var = values + np.sqrt(w[i]) * v[i]
df = np.sum(values) - np.sum(var)
err = np.sqrt(err**2 + df**2)
print(err)
```
%% Output
1.58113883008419
%% Cell type:code id:d98d47f5-cff0-4d89-984c-24a288935c7e tags:
``` python
```
......@@ -31,25 +31,25 @@ __all__ = [
_slow_pi_table_manager = TableManager.instance() # type: ignore
slow_pi_table_manager = _slow_pi_table_manager # type: TableManager
slow_pi_table_manager: TableManager = _slow_pi_table_manager
class SlowPionCorrection:
number_of_variations = 12 # type: int
number_of_variations: int = 12
def __init__(self) -> None:
self.weight_name = "slow_pi" # type: str
self.weight_name: str = "slow_pi"
self.table_name = "slowPionAllWeights" # type: str
self.base_path = os.path.dirname(os.path.abspath(__file__)) # type: PathType
self.table_name: str = "slow_pion_correction_table"
self.base_path: PathType = os.path.dirname(os.path.abspath(__file__))
self.columns = self.get_slow_pion_table_columns()
self.momentum_column = "momentum_lower" # type: str
self.charge_column = "charge" # type: str
self.exp_column = "exp_lower" # type: str
self.momentum_column: str = "momentum_lower"
self.charge_column: str = "charge"
self.exp_column: str = "exp_lower"
self.momentum_bins = self.get_momentum_bins() # type: List[float]
self.momentum_bins: List[float] = self.get_momentum_bins()
def get_slow_pion_correction(
self,
......@@ -155,8 +155,8 @@ class SlowPionCorrection:
return correction_matrix[mom_indices]
def _get_correction_matrix(self, sub_table_name: str, correction_col: str) -> np.ndarray:
matrix_name = f"{self.table_name}_{sub_table_name}_{correction_col}" # type: str
matrix = slow_pi_table_manager.get_correction_matrix(name=matrix_name) # type: np.ndarray
matrix_name: str = f"{self.table_name}_{sub_table_name}_{correction_col}"
matrix: np.ndarray = slow_pi_table_manager.get_correction_matrix(name=matrix_name)
if matrix is None:
table = self._get_slow_pi_table(sub_table_name=sub_table_name)
......@@ -166,7 +166,7 @@ class SlowPionCorrection:
return matrix
def _get_slow_pi_table(self, sub_table_name: str) -> pd.DataFrame:
table_path = os.path.join(self.base_path, f"{self.table_name}.csv") # type: PathType
table_path: PathType = os.path.join(self.base_path, f"{self.table_name}.csv")
table = slow_pi_table_manager.get_table(
name=self.table_name,
......@@ -191,7 +191,7 @@ class SlowPionCorrection:
# Ensure that momentum is unique in table:
assert not (table.groupby(self.momentum_column).size() != 1).any()
matrix = np.full(len(self.momentum_bins[:-1]), np.nan) # type: np.ndarray
matrix: np.ndarray = np.full(len(self.momentum_bins[:-1]), np.nan)
for mom_index, mom_bin in enumerate(self.momentum_bins[:-1]):
value_rows = table.loc[(table[self.momentum_column] == mom_bin)]
......
pmin, pmax, correction, err_uncorr, err_corr, err_syst
0.050, 0.075, 1.012, 0.204, 0.038, 0.001
0.075, 0.100, 1.302, 0.121, 0.053, 0.001
0.100, 0.125, 0.782, 0.090, 0.031, 0.001
0.125, 0.150, 0.880, 0.103, 0.033, 0.001
0.150, 0.175, 1.245, 0.093, 0.049, 0.001
0.175, 0.200, 1.116, 0.102, 0.039, 0.002
% correction factors for slow neutral pions, experiment range 7-27
% correction value = epsion_data / epsilon_MC
% pmin(GeV) pmax(GeV) correction err_uncorr err_corr err_syst
0.050 0.075 1.012 0.204 0.038 0.001
0.075 0.100 1.302 0.121 0.053 0.001
0.100 0.125 0.782 0.090 0.031 0.001
0.125 0.150 0.880 0.103 0.033 0.001
0.150 0.175 1.245 0.093 0.049 0.001
0.175 0.200 1.116 0.102 0.039 0.002
% 0 0.200 1.051 0.043 0.040 0.001
pmin, pmax, correction, err_uncorr, err_corr, err_syst
0.050, 0.075, 0.942, 0.086, 0.016, 0.001
0.075, 0.100, 1.067, 0.048, 0.019, 0.007
0.100, 0.125, 1.011, 0.039, 0.018, 0.001
0.125, 0.150, 0.969, 0.036, 0.017, 0.001
0.150, 0.175, 1.097, 0.033, 0.019, 0.000
0.175, 0.200, 0.986, 0.034, 0.017, 0.002
% correction factors for slow neutral pions, experiment range 31-65
% correction value = epsion_data / epsilon_MC
% pmin(GeV) pmax(GeV) correction err_uncorr err_corr err_syst
0.050 0.075 0.942 0.086 0.016 0.001
0.075 0.100 1.067 0.048 0.019 0.007
0.100 0.125 1.011 0.039 0.018 0.001
0.125 0.150 0.969 0.036 0.017 0.001
0.150 0.175 1.097 0.033 0.019 0.000
0.175 0.200 0.986 0.034 0.017 0.002
% 0 0.200 1.023 0.016 0.018 0.001
pmin, pmax, correction, err_uncorr, err_corr, err_syst
0.050, 0.075, 1.091, 0.182, 0.027, 0.013
0.075, 0.100, 0.937, 0.067, 0.024, 0.001
0.100, 0.125, 0.892, 0.053, 0.024, 0.001
0.125, 0.150, 1.058, 0.046, 0.028, 0.000
0.150, 0.175, 1.086, 0.042, 0.029, 0.000
0.175, 0.200, 1.028, 0.050, 0.027, 0.001
% correction factors for slow charged pions, experiment range 7-27
% correction value = epsion_data / epsilon_MC
% pmin(GeV) pmax(GeV) correction err_uncorr err_corr err_syst
0.050 0.075 1.091 0.182 0.027 0.013
0.075 0.100 0.937 0.067 0.024 0.001
0.100 0.125 0.892 0.053 0.024 0.001
0.125 0.150 1.058 0.046 0.028 0.000
0.150 0.175 1.086 0.042 0.029 0.000
0.175 0.200 1.028 0.050 0.027 0.001
% 0 0.200 1.020 0.022 0.027 0.000
pmin, pmax, correction, err_uncorr, err_corr, err_syst
0.050, 0.075, 0.832, 0.070, 0.009, 0.001
0.075, 0.100, 0.930, 0.027, 0.010, 0.001
0.100, 0.125, 0.967, 0.021, 0.010, 0.000
0.125, 0.150, 0.984, 0.017, 0.010, 0.000
0.150, 0.175, 1.009, 0.016, 0.011, 0.000
0.175, 0.200, 1.008, 0.018, 0.011, 0.000
% correction factors for slow charged pions, experiment range 31-65
% correction value = epsion_data / epsilon_MC
% pmin(GeV) pmax(GeV) correction err_uncorr err_corr err_syst
0.050 0.075 0.832 0.070 0.009 0.001
0.075 0.100 0.930 0.027 0.010 0.001
0.100 0.125 0.967 0.021 0.010 0.000
0.125 0.150 0.984 0.017 0.010 0.000
0.150 0.175 1.009 0.016 0.011 0.000
0.175 0.200 1.008 0.018 0.011 0.000
% 0 0.200 0.986 0.008 0.011 0.000
"""
Slow Pion correction table extension by calculating variations of the correction factors from the related uncertainties.
Calculates the variations of the slow pion correction table.
Sources:
General: https://belle.kek.jp/secured/wiki/doku.php?id=software:systematics
Data: https://belle.kek.jp/secured/belle_note/gn1176/
Belle Note: https://belle.kek.jp/secured/belle_note/gn1176/bn1176.pdf
Talk: https://kds.kek.jp/event/6029/contributions/138416/attachments/107716/127702/BAM-2010-11-15-WolfgangDungel.pdf
Felix Metzner 2023
"""
import os
import pandas as pd
from typing import List
from dataclasses import dataclass
from analysistools.utilities import PathType
__all__ = [
"SlowPionTablePreparer"
]
@dataclass(frozen=True)
class SlowPionTablesContainer:
pi_plus_old_svd: pd.DataFrame
pi_plus_new_svd: pd.DataFrame
pi_zero_old_svd: pd.DataFrame
pi_zero_new_svd: pd.DataFrame
@staticmethod
def load_tables() -> "SlowPionTablesContainer":
new_container = SlowPionTablesContainer(
pi_plus_old_svd=pi_plus_old_svd_table,
pi_plus_new_svd=pi_plus_new_svd_table,
pi_zero_old_svd=pi_zero_old_svd_table,
pi_zero_new_svd=pi_zero_new_svd_table,
)
return new_container
class SlowPionTablePreparer:
number_of_variations: int = 12
def __init__(self) -> None:
self.weight_name: str = "slow_pi"
self.table_name: str = "slowPionAllWeights"
self.base_path: PathType = os.path.dirname(os.path.abspath(__file__))
self.columns = self.get_slow_pion_table_columns()
self.momentum_column: str = "momentum_lower"
self.charge_column: str = "charge"
self.exp_column: str = "exp_lower"
self.momentum_bins: List[float] = self.get_momentum_bins()
def load_base_tables(self) -> SlowPionTablesContainer:
pass
def create_full_table_with_variations(self) -> pd.DataFrame:
pass
def get_varied_table(self) -> pd.DataFrame:
pass
def run_slow_pion_table_preparer() -> None:
pass
if __name__ == "__main__":
run_slow_pion_table_preparer()
%% Cell type:code id:c9827745fb474180 tags:
``` python
import numpy as np
import copy
```
%% Cell type:code id:6bc8fff55f63e701 tags:
``` python
def get_varied_central_values(
central: np.ndarray,
covariance_matrix: np.ndarray,
):
eigenvalues, eigenvectors = np.linalg.eig(covariance_matrix)
eigenvalue_variations = []
for index, eigenvalue in enumerate(eigenvalues):
variation = np.sqrt(eigenvalue)
up = central + variation * eigenvectors[:, index]
down = central - variation * eigenvectors[:, index]
eigenvalue_variations.append((up, down))
return eigenvalue_variations
```
%% Cell type:code id:283d796208b68fa8 tags:
``` python
```
%% Cell type:code id:36bdf3465840505b tags:
``` python
```
%% Cell type:code id:508f03ef694dafe4 tags:
``` python
def get_other_varied_central_values(
central: np.ndarray,
covariance_matrix: np.ndarray,
):
eigenvalues, eigenvectors_t = np.linalg.eig(covariance_matrix)
eigenvectors = eigenvectors_t.T
eigenvectors_inv = np.linalg.inv(eigenvectors)
fluctuations = np.sqrt(eigenvalues)
rho_tilde = np.matmul(central, eigenvectors_inv)
eigenvalue_variations = []
dim = fluctuations.shape[0]
for n in range(dim):
rho_tilde_temp = copy.deepcopy(rho_tilde)
rho_tilde_temp[n] += fluctuations[n] # err up
rho_prime = np.matmul(rho_tilde_temp, eigenvectors) # transfer to original basis
eigenvalue_variations.append(rho_prime)
rho_tilde_temp = copy.deepcopy(rho_tilde)
rho_tilde_temp[n-dim] -= fluctuations[n-dim] # err down
rho_prime = np.matmul(rho_tilde_temp, eigenvectors) # transfer to original basis
eigenvalue_variations.append(rho_prime)
return eigenvalue_variations
```
%% Cell type:code id:91c4dc7638d3e7b0 tags:
``` python
```
%% Cell type:code id:eab1ac53edcec01a tags:
``` python
central_vals = np.random.random(6)
cov_mat = np.random.random((6,6))
```
%% Cell type:code id:ef83100f015cabda tags:
``` python
central_vals
```
%% Output
array([0.62921316, 0.46960159, 0.66121828, 0.01703926, 0.3650362 ,
0.22322443])
%% Cell type:code id:2736cf0f1eed388a tags:
``` python
cov_mat[1, :]
```
%% Output
array([0.77512183, 0.58054596, 0.07501867, 0.73742079, 0.60000454,
0.11480988])
%% Cell type:code id:3fba6e7f499325d1 tags:
``` python
u_vec = np.array([0,1,0,0,0,0]).reshape(-1,1)
np.matmul(cov_mat, u_vec)
```
%% Output
array([[0.03680179],
[0.58054596],
[0.38189095],
[0.43605387],
[0.26525029],
[0.4395451 ]])
%% Cell type:code id:3b4e7bfbf9c92255 tags:
``` python
u_vec, np.transpose(u_vec)
```
%% Output
(array([[0],
[1],
[0],
[0],
[0],
[0]]),
array([[0, 1, 0, 0, 0, 0]]))
%% Cell type:code id:d290074692f79a04 tags:
``` python
```
%% Cell type:code id:125c62a61596964e tags:
``` python
np.matmul(u_vec, cov_mat)
```
%% Output
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/tmp/ipykernel_85113/682133678.py in <module>
----> 1 np.matmul(u_vec, cov_mat)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 6 is different from 1)
%% Cell type:code id:3fec546715825640 tags:
``` python
res_1 = get_varied_central_values(central_vals, cov_mat)
# res_1 = [res_1[j][i] for j in range(6) for i in (0, 1)]
res_1 = np.vstack(res_1)
res_1
```
%% Output
array([[ 1.22397216+0.j , 1.06804741+0.j ,
1.69066955+0.j , 0.67495151+0.j ,
0.87285429+0.j , 0.90788021+0.j ],
[ 0.03445415+0.j , -0.12884423+0.j ,
-0.368233 +0.j , -0.64087299+0.j ,
-0.1427819 +0.j , -0.46143135+0.j ],
[ 0.62921316-0.09775329j, 0.46960159-0.19403334j,
0.66121828+0.27223614j, 0.01703926+0.44250219j,
0.3650362 -0.02431763j, 0.22322443-0.4186711j ],
[ 0.62921316+0.09775329j, 0.46960159+0.19403334j,
0.66121828-0.27223614j, 0.01703926-0.44250219j,
0.3650362 +0.02431763j, 0.22322443+0.4186711j ],
[ 0.83978857-0.3582239j , 0.34905549-0.1772789j ,
0.6945451 -0.22426854j, -0.16904189+0.13850123j,
0.67935625+0.30825998j, -0.05564205+0.25028446j],
[ 0.41863775+0.3582239j , 0.59014768+0.1772789j ,
0.62789145+0.22426854j, 0.20312041-0.13850123j,
0.05071614-0.30825998j, 0.50209091-0.25028446j],
[ 0.83978857+0.3582239j , 0.34905549+0.1772789j ,
0.6945451 +0.22426854j, -0.16904189-0.13850123j,
0.67935625-0.30825998j, -0.05564205-0.25028446j],
[ 0.41863775-0.3582239j , 0.59014768-0.1772789j ,
0.62789145-0.22426854j, 0.20312041+0.13850123j,
0.05071614+0.30825998j, 0.50209091+0.25028446j],
[ 0.54011442+0.j , 0.86812894+0.j ,
0.69639245+0.j , -0.101313 +0.j ,
0.37307791+0.j , 0.1535801 +0.j ],
[ 0.7183119 +0.j , 0.07107423+0.j ,
0.6260441 +0.j , 0.13539151+0.j ,
0.35699448+0.j , 0.29286876+0.j ],
[ 0.79773146+0.j , -0.05256956+0.j ,
0.80039059+0.j , 0.07399786+0.j ,
0.25007745+0.j , 0.30629544+0.j ],
[ 0.46069486+0.j , 0.99177273+0.j ,
0.52204596+0.j , -0.03991934+0.j ,
0.47999494+0.j , 0.14015343+0.j ]])
%% Cell type:code id:44113454ce8d6f2f tags:
``` python
res_2 = get_other_varied_central_values(central_vals, cov_mat)
res_2 = np.vstack(res_2)
res_2
```
%% Output
array([[ 1.22397216+1.21393680e-16j, 1.06804741-2.61396639e-17j,
1.69066955-8.35812890e-18j, 0.67495151-1.02715743e-16j,
0.87285429+1.88942925e-16j, 0.90788021-1.21384824e-16j],
[ 0.03445415+1.21393680e-16j, -0.12884423-2.61396639e-17j,
-0.368233 -8.35812890e-18j, -0.64087299-1.02715743e-16j,
-0.1427819 +1.88942925e-16j, -0.46143135-1.21384824e-16j],
[ 0.62921316-9.77532894e-02j, 0.46960159-1.94033339e-01j,
0.66121828+2.72236135e-01j, 0.01703926+4.42502188e-01j,
0.3650362 -2.43176296e-02j, 0.22322443-4.18671096e-01j],
[ 0.62921316+9.77532894e-02j, 0.46960159+1.94033339e-01j,
0.66121828-2.72236135e-01j, 0.01703926-4.42502188e-01j,
0.3650362 +2.43176296e-02j, 0.22322443+4.18671096e-01j],
[ 0.83978857-3.58223901e-01j, 0.34905549-1.77278895e-01j,
0.6945451 -2.24268537e-01j, -0.16904189+1.38501228e-01j,
0.67935625+3.08259977e-01j, -0.05564205+2.50284457e-01j],
[ 0.41863775+3.58223901e-01j, 0.59014768+1.77278895e-01j,
0.62789145+2.24268537e-01j, 0.20312041-1.38501228e-01j,
0.05071614-3.08259977e-01j, 0.50209091-2.50284457e-01j],
[ 0.83978857+3.58223901e-01j, 0.34905549+1.77278895e-01j,
0.6945451 +2.24268537e-01j, -0.16904189-1.38501228e-01j,
0.67935625-3.08259977e-01j, -0.05564205-2.50284457e-01j],
[ 0.41863775-3.58223901e-01j, 0.59014768-1.77278895e-01j,
0.62789145-2.24268537e-01j, 0.20312041+1.38501228e-01j,
0.05071614+3.08259977e-01j, 0.50209091+2.50284457e-01j],
[ 0.54011442+1.21393680e-16j, 0.86812894-2.61396639e-17j,
0.69639245-8.35812890e-18j, -0.101313 -1.02715743e-16j,
0.37307791+1.88942925e-16j, 0.1535801 -1.21384824e-16j],
[ 0.7183119 +1.21393680e-16j, 0.07107423-2.61396639e-17j,
0.6260441 -8.35812890e-18j, 0.13539151-1.02715743e-16j,
0.35699448+1.88942925e-16j, 0.29286876-1.21384824e-16j],
[ 0.79773146+1.21393680e-16j, -0.05256956-2.61396639e-17j,
0.80039059-8.35812890e-18j, 0.07399786-1.02715743e-16j,
0.25007745+1.88942925e-16j, 0.30629544-1.21384824e-16j],
[ 0.46069486+1.21393680e-16j, 0.99177273-2.61396639e-17j,
0.52204596-8.35812890e-18j, -0.03991934-1.02715743e-16j,
0.47999494+1.88942925e-16j, 0.14015343-1.21384824e-16j]])
%% Cell type:code id:64809caedf704cbd tags:
``` python
np.diag(res_2)
```
%% Output
array([ 1.22397216+1.21393680e-16j, -0.12884423-2.61396639e-17j,
0.66121828+2.72236135e-01j, 0.01703926-4.42502188e-01j,
0.67935625+3.08259977e-01j, 0.50209091-2.50284457e-01j])
%% Cell type:code id:ba95733323f1bae3 tags:
``` python
np.diag(np.vstack(res_2))
```
%% Output
array([0.45187585-2.77555756e-17j, 0.355542 +2.63940903e-01j,
0.63257577-1.74676672e-01j, 0.93980196-5.55111512e-17j,
0.39819797-9.86662690e-03j, 0.67263833-1.55124705e-01j])
%% Cell type:code id:79a4df7e47b7c678 tags:
``` python
```
%% Output
array([ 0.45187585+0.j , -0.44006762+0.j ,
0.63257577+0.17467667j, 0.76618608+0.27279049j,
0.57583675+0.20229412j, 0.92065342-0.1288329j ])
%% Cell type:code id:eab954b03bbfab53 tags:
``` python
```