Newer
Older
Master_thesis / data / CLs / finished / f1d1 / -+ / D-True / 2288030 / raremodel-nb.py
@Sascha Liechti Sascha Liechti on 8 Oct 2019 74 KB ...
  1. #!/usr/bin/env python
  2. # coding: utf-8
  3.  
  4. # # Import
  5.  
  6. # In[ ]:
  7.  
  8.  
  9. import os
  10.  
  11. # os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
  12. import random
  13. import numpy as np
  14. from pdg_const1 import pdg
  15. import matplotlib
  16. import matplotlib.pyplot as plt
  17. import pickle as pkl
  18. import sys
  19. import time
  20. from helperfunctions import display_time, prepare_plot
  21. import cmath as c
  22. import scipy.integrate as integrate
  23. from scipy.optimize import fminbound
  24. from array import array as arr
  25. import collections
  26. from itertools import compress
  27. import tensorflow as tf
  28. import zfit
  29. from zfit import ztf
  30. # from IPython.display import clear_output
  31. import os
  32. import tensorflow_probability as tfp
  33. tfd = tfp.distributions
  34.  
  35.  
  36. # In[ ]:
  37.  
  38.  
  39. # chunksize = 10000
  40. # zfit.run.chunking.active = True
  41. # zfit.run.chunking.max_n_points = chunksize
  42.  
  43.  
  44. # # Build model and graphs
  45. # ## Create graphs
  46.  
  47. # In[ ]:
  48.  
  49.  
  50.  
  51.  
  52.  
  53. # In[ ]:
  54.  
  55.  
  56. def formfactor(q2, subscript, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2): #returns real value
  57. #check if subscript is viable
  58.  
  59. if subscript != "0" and subscript != "+" and subscript != "T":
  60. raise ValueError('Wrong subscript entered, choose either 0, + or T')
  61.  
  62. #get constants
  63.  
  64. mK = ztf.constant(pdg['Ks_M'])
  65. mbstar0 = ztf.constant(pdg["mbstar0"])
  66. mbstar = ztf.constant(pdg["mbstar"])
  67.  
  68.  
  69. mmu = ztf.constant(pdg['muon_M'])
  70. mb = ztf.constant(pdg['bquark_M'])
  71. ms = ztf.constant(pdg['squark_M'])
  72. mB = ztf.constant(pdg['Bplus_M'])
  73.  
  74. #N comes from derivation in paper
  75.  
  76. N = 3
  77.  
  78. #some helperfunctions
  79.  
  80. tpos = (mB - mK)**2
  81. tzero = (mB + mK)*(ztf.sqrt(mB)-ztf.sqrt(mK))**2
  82.  
  83. z_oben = ztf.sqrt(tpos - q2) - ztf.sqrt(tpos - tzero)
  84. z_unten = ztf.sqrt(tpos - q2) + ztf.sqrt(tpos - tzero)
  85. z = tf.divide(z_oben, z_unten)
  86.  
  87. #calculate f0
  88.  
  89. if subscript == "0":
  90. prefactor = 1/(1 - q2/(mbstar0**2))
  91. _sum = 0
  92. b0 = [b0_0, b0_1, b0_2]
  93.  
  94. for i in range(N):
  95. _sum += b0[i]*(tf.pow(z,i))
  96.  
  97. return ztf.to_complex(prefactor * _sum)
  98.  
  99. #calculate f+ or fT
  100.  
  101. else:
  102. prefactor = 1/(1 - q2/(mbstar**2))
  103. _sum = 0
  104.  
  105. if subscript == "T":
  106. bT = [bT_0, bT_1, bT_2]
  107. for i in range(N):
  108. _sum += bT[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))
  109. else:
  110. bplus = [bplus_0, bplus_1, bplus_2]
  111. for i in range(N):
  112. _sum += bplus[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))
  113.  
  114. return ztf.to_complex(prefactor * _sum)
  115.  
  116. def resonance(q, _mass, width, phase, scale):
  117.  
  118. q2 = tf.pow(q, 2)
  119.  
  120. mmu = ztf.constant(pdg['muon_M'])
  121.  
  122. p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2))
  123.  
  124. p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2)
  125.  
  126. gamma_j = tf.divide(p, q) * _mass * width / p0
  127.  
  128. #Calculate the resonance
  129.  
  130. _top = tf.complex(_mass * width, ztf.constant(0.0))
  131.  
  132. _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)
  133.  
  134. com = _top/_bottom
  135.  
  136. #Rotate by the phase
  137.  
  138. r = ztf.to_complex(scale*tf.abs(com))
  139.  
  140. _phase = tf.angle(com)
  141.  
  142. _phase += phase
  143.  
  144. com = r * tf.exp(tf.complex(ztf.constant(0.0), _phase))
  145.  
  146. return com
  147.  
  148.  
  149. def axiv_nonres(q, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):
  150.  
  151. GF = ztf.constant(pdg['GF'])
  152. alpha_ew = ztf.constant(pdg['alpha_ew'])
  153. Vtb = ztf.constant(pdg['Vtb'])
  154. Vts = ztf.constant(pdg['Vts'])
  155. C10eff = ztf.constant(pdg['C10eff'])
  156.  
  157. mmu = ztf.constant(pdg['muon_M'])
  158. mb = ztf.constant(pdg['bquark_M'])
  159. ms = ztf.constant(pdg['squark_M'])
  160. mK = ztf.constant(pdg['Ks_M'])
  161. mB = ztf.constant(pdg['Bplus_M'])
  162.  
  163. q2 = tf.pow(q, 2)
  164.  
  165. #Some helperfunctions
  166.  
  167. beta = 1. - 4. * mmu**2. / q2
  168.  
  169. kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)
  170.  
  171. #prefactor in front of whole bracket
  172.  
  173. prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)
  174.  
  175. #left term in bracket
  176.  
  177. bracket_left = 2./3. * tf.pow(kabs,2) * tf.pow(beta,2) * tf.pow(tf.abs(ztf.to_complex(C10eff)*formfactor(q2, "+", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)
  178.  
  179. #middle term in bracket
  180.  
  181. _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)
  182.  
  183. _under = q2 * mB**2.
  184.  
  185. bracket_middle = _top/_under *tf.pow(tf.abs(ztf.to_complex(C10eff) * formfactor(q2, "0", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)), 2)
  186. #Note sqrt(q2) comes from derivation as we use q2 and plot q
  187.  
  188. return prefactor1 * (bracket_left + bracket_middle) * 2 * q
  189.  
  190. def vec(q, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2):
  191. q2 = tf.pow(q, 2)
  192.  
  193. GF = ztf.constant(pdg['GF'])
  194. alpha_ew = ztf.constant(pdg['alpha_ew'])
  195. Vtb = ztf.constant(pdg['Vtb'])
  196. Vts = ztf.constant(pdg['Vts'])
  197. C7eff = ztf.constant(pdg['C7eff'])
  198.  
  199. mmu = ztf.constant(pdg['muon_M'])
  200. mb = ztf.constant(pdg['bquark_M'])
  201. ms = ztf.constant(pdg['squark_M'])
  202. mK = ztf.constant(pdg['Ks_M'])
  203. mB = ztf.constant(pdg['Bplus_M'])
  204.  
  205. #Some helperfunctions
  206.  
  207. beta = 1. - 4. * mmu**2. / q2
  208.  
  209. kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
  210. #prefactor in front of whole bracket
  211.  
  212. prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)
  213.  
  214. #right term in bracket
  215.  
  216. prefactor2 = tf.pow(kabs,2) * (1. - 1./3. * beta)
  217.  
  218. abs_bracket = tf.pow(tf.abs(c9eff(q, funcs) * formfactor(q2, "+", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2) + ztf.to_complex(2.0 * C7eff * (mb + ms)/(mB + mK)) * formfactor(q2, "T", b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)),2)
  219.  
  220. bracket_right = prefactor2 * abs_bracket
  221.  
  222. #Note sqrt(q2) comes from derivation as we use q2 and plot q
  223.  
  224. return prefactor1 * bracket_right * 2 * q
  225.  
  226. def c9eff(q, funcs):
  227.  
  228. C9eff_nr = ztf.to_complex(ztf.constant(pdg['C9eff']))
  229.  
  230. c9 = C9eff_nr + funcs
  231.  
  232. return c9
  233.  
  234.  
  235. # In[ ]:
  236.  
  237.  
  238. def G(y):
  239. def inner_rect_bracket(q):
  240. return tf.log(ztf.to_complex((1+tf.sqrt(q))/(1-tf.sqrt(q)))-tf.complex(ztf.constant(0), -1*ztf.constant(np.pi)))
  241. def inner_right(q):
  242. return ztf.to_complex(2 * tf.atan(1/tf.sqrt(tf.math.real(-q))))
  243. big_bracket = tf.where(tf.math.real(y) > ztf.constant(0.0), inner_rect_bracket(y), inner_right(y))
  244. return ztf.to_complex(tf.sqrt(tf.abs(y))) * big_bracket
  245.  
  246. def h_S(m, q):
  247. return ztf.to_complex(2) - G(ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2)))
  248.  
  249. def h_P(m, q):
  250. return ztf.to_complex(2/3) + (ztf.to_complex(1) - ztf.to_complex(4*tf.pow(m, 2)) / ztf.to_complex(tf.pow(q, 2))) * h_S(m,q)
  251.  
  252. def two_p_ccbar(mD, m_D_bar, m_D_star, q):
  253. #Load constants
  254. nu_D_bar = ztf.to_complex(pdg["nu_D_bar"])
  255. nu_D = ztf.to_complex(pdg["nu_D"])
  256. nu_D_star = ztf.to_complex(pdg["nu_D_star"])
  257. phase_D_bar = ztf.to_complex(pdg["phase_D_bar"])
  258. phase_D = ztf.to_complex(pdg["phase_D"])
  259. phase_D_star = ztf.to_complex(pdg["phase_D_star"])
  260. #Calculation
  261. left_part = nu_D_bar * tf.exp(tf.complex(ztf.constant(0.0), phase_D_bar)) * h_S(m_D_bar, q)
  262. right_part_D = nu_D * tf.exp(tf.complex(ztf.constant(0.0), phase_D)) * h_P(m_D, q)
  263. right_part_D_star = nu_D_star * tf.exp(tf.complex(ztf.constant(0.0), phase_D_star)) * h_P(m_D_star, q)
  264.  
  265. return left_part + right_part_D + right_part_D_star
  266.  
  267.  
  268. # ## Build pdf
  269.  
  270. # In[ ]:
  271.  
  272.  
  273. class total_pdf_cut(zfit.pdf.ZPDF):
  274. _N_OBS = 1 # dimension, can be omitted
  275. _PARAMS = ['b0_0', 'b0_1', 'b0_2',
  276. 'bplus_0', 'bplus_1', 'bplus_2',
  277. 'bT_0', 'bT_1', 'bT_2',
  278. 'rho_mass', 'rho_scale', 'rho_phase', 'rho_width',
  279. 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',
  280. 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',
  281. 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',
  282. 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',
  283. 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',
  284. 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',
  285. 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',
  286. 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width',
  287. 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',
  288. 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',
  289. 'tau_mass', 'C_tt']
  290. # the name of the parameters
  291.  
  292. def _unnormalized_pdf(self, x):
  293. x = x.unstack_x()
  294. b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]
  295. bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]
  296. bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]
  297. def rho_res(q):
  298. return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],
  299. phase = self.params['rho_phase'], width = self.params['rho_width'])
  300. def omega_res(q):
  301. return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],
  302. phase = self.params['omega_phase'], width = self.params['omega_width'])
  303. def phi_res(q):
  304. return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],
  305. phase = self.params['phi_phase'], width = self.params['phi_width'])
  306.  
  307. def jpsi_res(q):
  308. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'],
  309. scale = self.params['jpsi_scale'],
  310. phase = self.params['jpsi_phase'],
  311. width = self.params['jpsi_width'])
  312. def psi2s_res(q):
  313. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'],
  314. scale = self.params['psi2s_scale'],
  315. phase = self.params['psi2s_phase'],
  316. width = self.params['psi2s_width'])
  317. def p3770_res(q):
  318. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'],
  319. scale = self.params['p3770_scale'],
  320. phase = self.params['p3770_phase'],
  321. width = self.params['p3770_width'])
  322. def p4040_res(q):
  323. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'],
  324. scale = self.params['p4040_scale'],
  325. phase = self.params['p4040_phase'],
  326. width = self.params['p4040_width'])
  327. def p4160_res(q):
  328. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'],
  329. scale = self.params['p4160_scale'],
  330. phase = self.params['p4160_phase'],
  331. width = self.params['p4160_width'])
  332. def p4415_res(q):
  333. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'],
  334. scale = self.params['p4415_scale'],
  335. phase = self.params['p4415_phase'],
  336. width = self.params['p4415_width'])
  337. def P2_D(q):
  338. Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))
  339. DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*(ztf.to_complex(h_P(self.params['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_mass'], q)))
  340. return Dbar_contrib + DDstar_contrib
  341. def ttau_cusp(q):
  342. return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))
  343.  
  344. funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(x) + P2_D(x) + ttau_cusp(x)
  345.  
  346. vec_f = vec(x, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)
  347.  
  348. axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)
  349.  
  350. tot = vec_f + axiv_nr
  351. #Cut out jpsi and psi2s
  352. tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)
  353. tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)
  354. return tot
  355. class total_pdf_full(zfit.pdf.ZPDF):
  356. _N_OBS = 1 # dimension, can be omitted
  357. _PARAMS = ['b0_0', 'b0_1', 'b0_2',
  358. 'bplus_0', 'bplus_1', 'bplus_2',
  359. 'bT_0', 'bT_1', 'bT_2',
  360. 'rho_mass', 'rho_scale', 'rho_phase', 'rho_width',
  361. 'jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',
  362. 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',
  363. 'p3770_mass', 'p3770_scale', 'p3770_phase', 'p3770_width',
  364. 'p4040_mass', 'p4040_scale', 'p4040_phase', 'p4040_width',
  365. 'p4160_mass', 'p4160_scale', 'p4160_phase', 'p4160_width',
  366. 'p4415_mass', 'p4415_scale', 'p4415_phase', 'p4415_width',
  367. 'omega_mass', 'omega_scale', 'omega_phase', 'omega_width',
  368. 'phi_mass', 'phi_scale', 'phi_phase', 'phi_width',
  369. 'Dbar_mass', 'Dbar_scale', 'Dbar_phase',
  370. 'Dstar_mass', 'DDstar_scale', 'DDstar_phase', 'D_mass',
  371. 'tau_mass', 'C_tt']
  372. # the name of the parameters
  373.  
  374. def _unnormalized_pdf(self, x):
  375. x = x.unstack_x()
  376. b0 = [self.params['b0_0'], self.params['b0_1'], self.params['b0_2']]
  377. bplus = [self.params['bplus_0'], self.params['bplus_1'], self.params['bplus_2']]
  378. bT = [self.params['bT_0'], self.params['bT_1'], self.params['bT_2']]
  379. def rho_res(q):
  380. return resonance(q, _mass = self.params['rho_mass'], scale = self.params['rho_scale'],
  381. phase = self.params['rho_phase'], width = self.params['rho_width'])
  382. def omega_res(q):
  383. return resonance(q, _mass = self.params['omega_mass'], scale = self.params['omega_scale'],
  384. phase = self.params['omega_phase'], width = self.params['omega_width'])
  385. def phi_res(q):
  386. return resonance(q, _mass = self.params['phi_mass'], scale = self.params['phi_scale'],
  387. phase = self.params['phi_phase'], width = self.params['phi_width'])
  388.  
  389. def jpsi_res(q):
  390. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['jpsi_mass'], 2)) * resonance(q, _mass = self.params['jpsi_mass'],
  391. scale = self.params['jpsi_scale'],
  392. phase = self.params['jpsi_phase'],
  393. width = self.params['jpsi_width'])
  394. def psi2s_res(q):
  395. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['psi2s_mass'], 2)) * resonance(q, _mass = self.params['psi2s_mass'],
  396. scale = self.params['psi2s_scale'],
  397. phase = self.params['psi2s_phase'],
  398. width = self.params['psi2s_width'])
  399. def p3770_res(q):
  400. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p3770_mass'], 2)) * resonance(q, _mass = self.params['p3770_mass'],
  401. scale = self.params['p3770_scale'],
  402. phase = self.params['p3770_phase'],
  403. width = self.params['p3770_width'])
  404. def p4040_res(q):
  405. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4040_mass'], 2)) * resonance(q, _mass = self.params['p4040_mass'],
  406. scale = self.params['p4040_scale'],
  407. phase = self.params['p4040_phase'],
  408. width = self.params['p4040_width'])
  409. def p4160_res(q):
  410. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4160_mass'], 2)) * resonance(q, _mass = self.params['p4160_mass'],
  411. scale = self.params['p4160_scale'],
  412. phase = self.params['p4160_phase'],
  413. width = self.params['p4160_width'])
  414. def p4415_res(q):
  415. return ztf.to_complex(tf.pow(q, 2) / tf.pow(self.params['p4415_mass'], 2)) * resonance(q, _mass = self.params['p4415_mass'],
  416. scale = self.params['p4415_scale'],
  417. phase = self.params['p4415_phase'],
  418. width = self.params['p4415_width'])
  419. def P2_D(q):
  420. Dbar_contrib = ztf.to_complex(self.params['Dbar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['Dbar_phase']))*ztf.to_complex(h_S(self.params['Dbar_mass'], q))
  421. DDstar_contrib = ztf.to_complex(self.params['DDstar_scale'])*tf.exp(tf.complex(ztf.constant(0.0), self.params['DDstar_phase']))*(ztf.to_complex(h_P(self.params['Dstar_mass'], q)) + ztf.to_complex(h_P(self.params['D_mass'], q)))
  422. return Dbar_contrib + DDstar_contrib
  423. def ttau_cusp(q):
  424. return ztf.to_complex(self.params['C_tt'])*(ztf.to_complex((h_S(self.params['tau_mass'], q))) - ztf.to_complex(h_P(self.params['tau_mass'], q)))
  425.  
  426. funcs = rho_res(x) + omega_res(x) + phi_res(x) + jpsi_res(x) + psi2s_res(x) + p3770_res(x) + p4040_res(x)+ p4160_res(x) + p4415_res(x) + P2_D(x) + ttau_cusp(x)
  427.  
  428. vec_f = vec(x, funcs, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)
  429.  
  430. axiv_nr = axiv_nonres(x, b0_0, b0_1, b0_2, bplus_0, bplus_1, bplus_2, bT_0, bT_1, bT_2)
  431.  
  432. tot = vec_f + axiv_nr
  433. #Cut out jpsi and psi2s
  434. # tot = tf.where(tf.math.logical_or(x < ztf.constant(jpsi_mass-60.), x > ztf.constant(jpsi_mass+70.)), tot, 0.0*tot)
  435. # tot = tf.where(tf.math.logical_or(x < ztf.constant(psi2s_mass-50.), x > ztf.constant(psi2s_mass+50.)), tot, 0.0*tot)
  436. return tot
  437.  
  438.  
  439. # ## Setup parameters
  440.  
  441. # In[ ]:
  442.  
  443.  
  444. # formfactors
  445.  
  446. b0_0 = zfit.Parameter("b0_0", ztf.constant(0.292), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)
  447. b0_1 = zfit.Parameter("b0_1", ztf.constant(0.281), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)
  448. b0_2 = zfit.Parameter("b0_2", ztf.constant(0.150), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)
  449.  
  450. bplus_0 = zfit.Parameter("bplus_0", ztf.constant(0.466), lower_limit = -2.0, upper_limit= 2.0)
  451. bplus_1 = zfit.Parameter("bplus_1", ztf.constant(-0.885), lower_limit = -2.0, upper_limit= 2.0)
  452. bplus_2 = zfit.Parameter("bplus_2", ztf.constant(-0.213), lower_limit = -2.0, upper_limit= 2.0)
  453.  
  454. bT_0 = zfit.Parameter("bT_0", ztf.constant(0.460), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)
  455. bT_1 = zfit.Parameter("bT_1", ztf.constant(-1.089), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)
  456. bT_2 = zfit.Parameter("bT_2", ztf.constant(-1.114), floating = False) #, lower_limit = -2.0, upper_limit= 2.0)
  457.  
  458.  
  459. #rho
  460.  
  461. rho_mass, rho_width, rho_phase, rho_scale = pdg["rho"]
  462.  
  463. rho_m = zfit.Parameter("rho_m", ztf.constant(rho_mass), floating = False) #lower_limit = rho_mass - rho_width, upper_limit = rho_mass + rho_width)
  464. rho_w = zfit.Parameter("rho_w", ztf.constant(rho_width), floating = False)
  465. rho_p = zfit.Parameter("rho_p", ztf.constant(rho_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  466. rho_s = zfit.Parameter("rho_s", ztf.constant(rho_scale), lower_limit=rho_scale-np.sqrt(rho_scale), upper_limit=rho_scale+np.sqrt(rho_scale))
  467.  
  468. #omega
  469.  
  470. omega_mass, omega_width, omega_phase, omega_scale = pdg["omega"]
  471.  
  472. omega_m = zfit.Parameter("omega_m", ztf.constant(omega_mass), floating = False)
  473. omega_w = zfit.Parameter("omega_w", ztf.constant(omega_width), floating = False)
  474. omega_p = zfit.Parameter("omega_p", ztf.constant(omega_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  475. omega_s = zfit.Parameter("omega_s", ztf.constant(omega_scale), lower_limit=omega_scale-np.sqrt(omega_scale), upper_limit=omega_scale+np.sqrt(omega_scale))
  476.  
  477.  
  478. #phi
  479.  
  480. phi_mass, phi_width, phi_phase, phi_scale = pdg["phi"]
  481.  
  482. phi_m = zfit.Parameter("phi_m", ztf.constant(phi_mass), floating = False)
  483. phi_w = zfit.Parameter("phi_w", ztf.constant(phi_width), floating = False)
  484. phi_p = zfit.Parameter("phi_p", ztf.constant(phi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  485. phi_s = zfit.Parameter("phi_s", ztf.constant(phi_scale), lower_limit=phi_scale-np.sqrt(phi_scale), upper_limit=phi_scale+np.sqrt(phi_scale))
  486.  
  487. #jpsi
  488.  
  489. jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg["jpsi"]
  490.  
  491. jpsi_m = zfit.Parameter("jpsi_m", ztf.constant(jpsi_mass), floating = False)
  492. jpsi_w = zfit.Parameter("jpsi_w", ztf.constant(jpsi_width), floating = False)
  493. jpsi_p = zfit.Parameter("jpsi_p", ztf.constant(jpsi_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  494. jpsi_s = zfit.Parameter("jpsi_s", ztf.constant(jpsi_scale), floating = False) #, lower_limit=jpsi_scale-np.sqrt(jpsi_scale), upper_limit=jpsi_scale+np.sqrt(jpsi_scale))
  495.  
  496. #psi2s
  497.  
  498. psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
  499.  
  500. psi2s_m = zfit.Parameter("psi2s_m", ztf.constant(psi2s_mass), floating = False)
  501. psi2s_w = zfit.Parameter("psi2s_w", ztf.constant(psi2s_width), floating = False)
  502. psi2s_p = zfit.Parameter("psi2s_p", ztf.constant(psi2s_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  503. psi2s_s = zfit.Parameter("psi2s_s", ztf.constant(psi2s_scale), floating = False) #, lower_limit=psi2s_scale-np.sqrt(psi2s_scale), upper_limit=psi2s_scale+np.sqrt(psi2s_scale))
  504.  
  505. #psi(3770)
  506.  
  507. p3770_mass, p3770_width, p3770_phase, p3770_scale = pdg["p3770"]
  508.  
  509. p3770_m = zfit.Parameter("p3770_m", ztf.constant(p3770_mass), floating = False)
  510. p3770_w = zfit.Parameter("p3770_w", ztf.constant(p3770_width), floating = False)
  511. p3770_p = zfit.Parameter("p3770_p", ztf.constant(p3770_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  512. p3770_s = zfit.Parameter("p3770_s", ztf.constant(p3770_scale), lower_limit=p3770_scale-np.sqrt(p3770_scale), upper_limit=p3770_scale+np.sqrt(p3770_scale))
  513.  
  514. #psi(4040)
  515.  
  516. p4040_mass, p4040_width, p4040_phase, p4040_scale = pdg["p4040"]
  517.  
  518. p4040_m = zfit.Parameter("p4040_m", ztf.constant(p4040_mass), floating = False)
  519. p4040_w = zfit.Parameter("p4040_w", ztf.constant(p4040_width), floating = False)
  520. p4040_p = zfit.Parameter("p4040_p", ztf.constant(p4040_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  521. p4040_s = zfit.Parameter("p4040_s", ztf.constant(p4040_scale), lower_limit=p4040_scale-np.sqrt(p4040_scale), upper_limit=p4040_scale+np.sqrt(p4040_scale))
  522.  
  523. #psi(4160)
  524.  
  525. p4160_mass, p4160_width, p4160_phase, p4160_scale = pdg["p4160"]
  526.  
  527. p4160_m = zfit.Parameter("p4160_m", ztf.constant(p4160_mass), floating = False)
  528. p4160_w = zfit.Parameter("p4160_w", ztf.constant(p4160_width), floating = False)
  529. p4160_p = zfit.Parameter("p4160_p", ztf.constant(p4160_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  530. p4160_s = zfit.Parameter("p4160_s", ztf.constant(p4160_scale), lower_limit=p4160_scale-np.sqrt(p4160_scale), upper_limit=p4160_scale+np.sqrt(p4160_scale))
  531.  
  532. #psi(4415)
  533.  
  534. p4415_mass, p4415_width, p4415_phase, p4415_scale = pdg["p4415"]
  535.  
  536. p4415_m = zfit.Parameter("p4415_m", ztf.constant(p4415_mass), floating = False)
  537. p4415_w = zfit.Parameter("p4415_w", ztf.constant(p4415_width), floating = False)
  538. p4415_p = zfit.Parameter("p4415_p", ztf.constant(p4415_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)
  539. p4415_s = zfit.Parameter("p4415_s", ztf.constant(p4415_scale), lower_limit=p4415_scale-np.sqrt(p4415_scale), upper_limit=p4415_scale+np.sqrt(p4415_scale))
  540.  
  541.  
  542. # ## Dynamic generation of 2 particle contribution
  543.  
  544. # In[ ]:
  545.  
  546.  
  547. m_c = 1300
  548.  
  549. Dbar_phase = 0.0
  550. DDstar_phase = 0.0
  551. Dstar_mass = pdg['Dst_M']
  552. Dbar_mass = pdg['D0_M']
  553. D_mass = pdg['D0_M']
  554.  
  555. Dbar_s = zfit.Parameter("Dbar_s", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)
  556. Dbar_m = zfit.Parameter("Dbar_m", ztf.constant(Dbar_mass), floating = False)
  557. Dbar_p = zfit.Parameter("Dbar_p", ztf.constant(Dbar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)
  558. DDstar_s = zfit.Parameter("DDstar_s", ztf.constant(0.0), lower_limit=-0.3, upper_limit=0.3)#, floating = False)
  559. Dstar_m = zfit.Parameter("Dstar_m", ztf.constant(Dstar_mass), floating = False)
  560. D_m = zfit.Parameter("D_m", ztf.constant(D_mass), floating = False)
  561. DDstar_p = zfit.Parameter("DDstar_p", ztf.constant(DDstar_phase), lower_limit=-2*np.pi, upper_limit=2*np.pi)#, floating = False)
  562.  
  563.  
  564. # ## Tau parameters
  565.  
  566. # In[ ]:
  567.  
  568.  
  569. tau_m = zfit.Parameter("tau_m", ztf.constant(pdg['tau_M']), floating = False)
  570. Ctt = zfit.Parameter("Ctt", ztf.constant(0.0), lower_limit=-2.5, upper_limit=2.5)
  571.  
  572.  
  573. # ## Load data
  574.  
  575. # In[ ]:
  576.  
  577.  
  578. x_min = 2*pdg['muon_M']
  579. x_max = (pdg["Bplus_M"]-pdg["Ks_M"]-0.1)
  580.  
  581. # # Full spectrum
  582.  
  583. obs_toy = zfit.Space('q', limits = (x_min, x_max))
  584.  
  585. # Jpsi and Psi2s cut out
  586.  
  587. obs1 = zfit.Space('q', limits = (x_min, jpsi_mass - 60.))
  588. obs2 = zfit.Space('q', limits = (jpsi_mass + 70., psi2s_mass - 50.))
  589. obs3 = zfit.Space('q', limits = (psi2s_mass + 50., x_max))
  590.  
  591. obs_fit = obs1 + obs2 + obs3
  592.  
  593. # with open(r"./data/slim_points/slim_points_toy_0_range({0}-{1}).pkl".format(int(x_min), int(x_max)), "rb") as input_file:
  594. # part_set = pkl.load(input_file)
  595.  
  596. # x_part = part_set['x_part']
  597.  
  598. # x_part = x_part.astype('float64')
  599.  
  600. # data = zfit.data.Data.from_numpy(array=x_part, obs=obs)
  601.  
  602.  
  603. # ## Setup pdf
  604.  
  605. # In[ ]:
  606.  
  607.  
  608. total_f = total_pdf_cut(obs=obs_toy, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,
  609. psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,
  610. p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,
  611. p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,
  612. p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,
  613. p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,
  614. rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,
  615. omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,
  616. phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w,
  617. Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,
  618. Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,
  619. tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,
  620. bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,
  621. bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)
  622.  
  623. total_f_fit = total_pdf_full(obs=obs_fit, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w,
  624. psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w,
  625. p3770_mass = p3770_m, p3770_scale = p3770_s, p3770_phase = p3770_p, p3770_width = p3770_w,
  626. p4040_mass = p4040_m, p4040_scale = p4040_s, p4040_phase = p4040_p, p4040_width = p4040_w,
  627. p4160_mass = p4160_m, p4160_scale = p4160_s, p4160_phase = p4160_p, p4160_width = p4160_w,
  628. p4415_mass = p4415_m, p4415_scale = p4415_s, p4415_phase = p4415_p, p4415_width = p4415_w,
  629. rho_mass = rho_m, rho_scale = rho_s, rho_phase = rho_p, rho_width = rho_w,
  630. omega_mass = omega_m, omega_scale = omega_s, omega_phase = omega_p, omega_width = omega_w,
  631. phi_mass = phi_m, phi_scale = phi_s, phi_phase = phi_p, phi_width = phi_w,
  632. Dstar_mass = Dstar_m, DDstar_scale = DDstar_s, DDstar_phase = DDstar_p, D_mass = D_m,
  633. Dbar_mass = Dbar_m, Dbar_scale = Dbar_s, Dbar_phase = Dbar_p,
  634. tau_mass = tau_m, C_tt = Ctt, b0_0 = b0_0, b0_1 = b0_1, b0_2 = b0_2,
  635. bplus_0 = bplus_0, bplus_1 = bplus_1, bplus_2 = bplus_2,
  636. bT_0 = bT_0, bT_1 = bT_1, bT_2 = bT_2)
  637. # print(total_pdf.obs)
  638.  
  639. # print(calcs_test)
  640.  
  641. # for param in total_f.get_dependents():
  642. # print(zfit.run(param))
  643.  
  644.  
  645. # In[ ]:
  646.  
  647.  
  648. # total_f_fit.normalization(obs_fit)
  649.  
  650.  
  651. # ## Test if graphs actually work and compute values
  652.  
  653. # In[ ]:
  654.  
  655.  
  656. # def total_test_tf(xq):
  657.  
  658. # def jpsi_res(q):
  659. # return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)
  660.  
  661. # def psi2s_res(q):
  662. # return resonance(q, psi2s_m, psi2s_s, psi2s_p, psi2s_w)
  663.  
  664. # def cusp(q):
  665. # return bifur_gauss(q, cusp_m, sig_L, sig_R, cusp_s)
  666.  
  667. # funcs = jpsi_res(xq) + psi2s_res(xq) + cusp(xq)
  668.  
  669. # vec_f = vec(xq, funcs)
  670.  
  671. # axiv_nr = axiv_nonres(xq)
  672.  
  673. # tot = vec_f + axiv_nr
  674. # return tot
  675.  
  676. # def jpsi_res(q):
  677. # return resonance(q, jpsi_m, jpsi_s, jpsi_p, jpsi_w)
  678.  
  679. # calcs = zfit.run(total_test_tf(x_part))
  680.  
  681.  
  682.  
  683. test_q = np.linspace(x_min, x_max, int(2e6))
  684.  
  685. probs = total_f_fit.pdf(test_q, norm_range=False)
  686.  
  687. calcs_test = zfit.run(probs)
  688.  
  689. Ctt.set_value(0.5)
  690.  
  691. calcs_test1 = zfit.run(probs)
  692.  
  693. Ctt.set_value(0.0)
  694.  
  695. Dbar_s.set_value(0.3)
  696.  
  697. DDstar_s.set_value(0.3)
  698.  
  699. calcs_test2 = zfit.run(probs)
  700. # res_y = zfit.run(jpsi_res(test_q))
  701. # b0 = [b0_0, b0_1, b0_2]
  702. # bplus = [bplus_0, bplus_1, bplus_2]
  703. # bT = [bT_0, bT_1, bT_2]
  704. # f0_y = zfit.run(tf.math.real(formfactor(test_q,"0", b0, bplus, bT)))
  705. # fplus_y = zfit.run(tf.math.real(formfactor(test_q,"+", b0, bplus, bT)))
  706. # fT_y = zfit.run(tf.math.real(formfactor(test_q,"T", b0, bplus, bT)))
  707.  
  708.  
  709. # In[ ]:
  710.  
  711.  
  712. plt.clf()
  713. # plt.plot(x_part, calcs, '.')
  714. plt.plot(test_q, calcs_test, label = 'pdf (Ctt = 0.0)')
  715. plt.plot(test_q, calcs_test1, label = 'pdf (Ctt = 0.5)')
  716. plt.plot(test_q, calcs_test2, label = 'pdf (D-contribs = 0.3)')
  717. # plt.plot(test_q, f0_y, label = '0')
  718. # plt.plot(test_q, fT_y, label = 'T')
  719. # plt.plot(test_q, fplus_y, label = '+')
  720. # plt.plot(test_q, res_y, label = 'res')
  721. plt.legend()
  722. plt.ylim(0.0, 1.5e-6)
  723. # plt.yscale('log')
  724. # plt.xlim(770, 785)
  725. plt.savefig('test.png')
  726. # print(jpsi_width)
  727.  
  728.  
  729. # In[ ]:
  730.  
  731.  
  732.  
  733.  
  734. # probs = mixture.prob(test_q)
  735. # probs_np = zfit.run(probs)
  736. # probs_np *= np.max(calcs_test) / np.max(probs_np)
  737. # plt.figure()
  738. # plt.semilogy(test_q, probs_np,label="importance sampling")
  739. # plt.semilogy(test_q, calcs_test, label = 'pdf')
  740.  
  741.  
  742. # In[ ]:
  743.  
  744.  
  745. # 0.213/(0.00133+0.213+0.015)
  746.  
  747.  
  748. # ## Adjust scaling of different parts
  749.  
  750. # In[ ]:
  751.  
  752.  
  753. total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)
  754. # inte = total_f.integrate(limits = (950., 1050.), norm_range=False)
  755. # inte_fl = zfit.run(inte)
  756. # print(inte_fl/4500)
  757. # print(pdg["jpsi_BR"]/pdg["NR_BR"], inte_fl*pdg["psi2s_auc"]/pdg["NR_auc"])
  758.  
  759.  
  760. # In[ ]:
  761.  
  762.  
  763. # # print("jpsi:", inte_fl)
  764. # # print("Increase am by factor:", np.sqrt(pdg["jpsi_BR"]/pdg["NR_BR"]*pdg["NR_auc"]/inte_fl))
  765. # # print("New amp:", pdg["jpsi"][3]*np.sqrt(pdg["jpsi_BR"]/pdg["NR_BR"]*pdg["NR_auc"]/inte_fl))
  766.  
  767. # # print("psi2s:", inte_fl)
  768. # # print("Increase am by factor:", np.sqrt(pdg["psi2s_BR"]/pdg["NR_BR"]*pdg["NR_auc"]/inte_fl))
  769. # # print("New amp:", pdg["psi2s"][3]*np.sqrt(pdg["psi2s_BR"]/pdg["NR_BR"]*pdg["NR_auc"]/inte_fl))
  770.  
  771. # name = "phi"
  772.  
  773. # print(name+":", inte_fl)
  774. # print("Increase am by factor:", np.sqrt(pdg[name+"_BR"]/pdg["NR_BR"]*pdg["NR_auc"]/inte_fl))
  775. # print("New amp:", pdg[name][0]*np.sqrt(pdg[name+"_BR"]/pdg["NR_BR"]*pdg["NR_auc"]/inte_fl))
  776.  
  777.  
  778. # print(x_min)
  779. # print(x_max)
  780. # # total_f.update_integration_options(draws_per_dim=2000000, mc_sampler=None)
  781. # total_f.update_integration_options(mc_sampler=lambda dim, num_results,
  782. # dtype: tf.random_uniform(maxval=1., shape=(num_results, dim), dtype=dtype),
  783. # draws_per_dim=1000000)
  784. # # _ = []
  785.  
  786. # # for i in range(10):
  787.  
  788. # # inte = total_f.integrate(limits = (x_min, x_max))
  789. # # inte_fl = zfit.run(inte)
  790. # # print(inte_fl)
  791. # # _.append(inte_fl)
  792.  
  793. # # print("mean:", np.mean(_))
  794.  
  795. # _ = time.time()
  796.  
  797. # inte = total_f.integrate(limits = (x_min, x_max))
  798. # inte_fl = zfit.run(inte)
  799. # print(inte_fl)
  800. # print("Time taken: {}".format(display_time(int(time.time() - _))))
  801.  
  802. # print(pdg['NR_BR']/pdg['NR_auc']*inte_fl)
  803. # print(0.25**2*4.2/1000)
  804.  
  805.  
  806. # # Sampling
  807. # ## Mixture distribution for sampling
  808.  
  809. # In[ ]:
  810.  
  811.  
  812.  
  813. # print(list_of_borders[:9])
  814. # print(list_of_borders[-9:])
  815.  
  816.  
  817. class UniformSampleAndWeights(zfit.util.execution.SessionHolderMixin):
  818. def __call__(self, limits, dtype, n_to_produce):
  819. # n_to_produce = tf.cast(n_to_produce, dtype=tf.int32)
  820. low, high = limits.limit1d
  821. low = tf.cast(low, dtype=dtype)
  822. high = tf.cast(high, dtype=dtype)
  823. # uniform = tfd.Uniform(low=low, high=high)
  824. # uniformjpsi = tfd.Uniform(low=tf.constant(3080, dtype=dtype), high=tf.constant(3112, dtype=dtype))
  825. # uniformpsi2s = tfd.Uniform(low=tf.constant(3670, dtype=dtype), high=tf.constant(3702, dtype=dtype))
  826.  
  827. # list_of_borders = []
  828. # _p = []
  829. # splits = 10
  830.  
  831. # _ = np.linspace(x_min, x_max, splits)
  832.  
  833. # for i in range(splits):
  834. # list_of_borders.append(tf.constant(_[i], dtype=dtype))
  835. # _p.append(tf.constant(1/splits, dtype=dtype))
  836. # mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=_p[:(splits-1)]),
  837. # components_distribution=tfd.Uniform(low=list_of_borders[:(splits-1)],
  838. # high=list_of_borders[-(splits-1):]))
  839. mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.05, dtype=dtype),
  840. tf.constant(0.93, dtype=dtype),
  841. tf.constant(0.05, dtype=dtype),
  842. tf.constant(0.065, dtype=dtype),
  843. tf.constant(0.04, dtype=dtype),
  844. tf.constant(0.05, dtype=dtype)]),
  845. components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype),
  846. tf.constant(3090, dtype=dtype),
  847. tf.constant(3681, dtype=dtype),
  848. tf.constant(3070, dtype=dtype),
  849. tf.constant(1000, dtype=dtype),
  850. tf.constant(3660, dtype=dtype)],
  851. high=[tf.constant(x_max, dtype=dtype),
  852. tf.constant(3102, dtype=dtype),
  853. tf.constant(3691, dtype=dtype),
  854. tf.constant(3110, dtype=dtype),
  855. tf.constant(1040, dtype=dtype),
  856. tf.constant(3710, dtype=dtype)]))
  857. # dtype = tf.float64
  858. # mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.04, dtype=dtype),
  859. # tf.constant(0.90, dtype=dtype),
  860. # tf.constant(0.02, dtype=dtype),
  861. # tf.constant(0.07, dtype=dtype),
  862. # tf.constant(0.02, dtype=dtype)]),
  863. # components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype),
  864. # tf.constant(3089, dtype=dtype),
  865. # tf.constant(3103, dtype=dtype),
  866. # tf.constant(3681, dtype=dtype),
  867. # tf.constant(3691, dtype=dtype)],
  868. # high=[tf.constant(3089, dtype=dtype),
  869. # tf.constant(3103, dtype=dtype),
  870. # tf.constant(3681, dtype=dtype),
  871. # tf.constant(3691, dtype=dtype),
  872. # tf.constant(x_max, dtype=dtype)]))
  873. # mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))
  874. # sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)
  875. sample = mixture.sample((n_to_produce, 1))
  876. # sample = tf.random.uniform((n_to_produce, 1), dtype=dtype)
  877. weights = mixture.prob(sample)[:,0]
  878. # weights = tf.broadcast_to(tf.constant(1., dtype=dtype), shape=(n_to_produce,))
  879. # sample = tf.expand_dims(sample, axis=-1)
  880. # print(sample, weights)
  881. # weights = tf.ones(shape=(n_to_produce,), dtype=dtype)
  882. weights_max = None
  883. thresholds = tf.random_uniform(shape=(n_to_produce,), dtype=dtype)
  884. return sample, thresholds, weights, weights_max, n_to_produce
  885.  
  886.  
  887. # In[ ]:
  888.  
  889.  
  890. # total_f._sample_and_weights = UniformSampleAndWeights
  891.  
  892.  
  893. # In[ ]:
  894.  
  895.  
  896. # 0.00133/(0.00133+0.213+0.015)*(x_max-3750)/(x_max-x_min)
  897.  
  898.  
  899. # In[ ]:
  900.  
  901.  
  902. # zfit.settings.set_verbosity(10)
  903.  
  904.  
  905. # In[ ]:
  906.  
  907.  
  908. # # zfit.run.numeric_checks = False
  909.  
  910. # nr_of_toys = 1
  911. # nevents = int(pdg["number_of_decays"])
  912. # nevents = pdg["number_of_decays"]
  913. # event_stack = 1000000
  914. # # zfit.settings.set_verbosity(10)
  915. # calls = int(nevents/event_stack + 1)
  916.  
  917. # total_samp = []
  918.  
  919. # start = time.time()
  920.  
  921. # sampler = total_f.create_sampler(n=event_stack)
  922.  
  923. # for toy in range(nr_of_toys):
  924. # dirName = 'data/zfit_toys/toy_{0}'.format(toy)
  925. # if not os.path.exists(dirName):
  926. # os.mkdir(dirName)
  927. # print("Directory " , dirName , " Created ")
  928.  
  929. # for call in range(calls):
  930.  
  931. # sampler.resample(n=event_stack)
  932. # s = sampler.unstack_x()
  933. # sam = zfit.run(s)
  934. # # clear_output(wait=True)
  935.  
  936. # c = call + 1
  937. # print("{0}/{1} of Toy {2}/{3}".format(c, calls, toy+1, nr_of_toys))
  938. # print("Time taken: {}".format(display_time(int(time.time() - start))))
  939. # print("Projected time left: {}".format(display_time(int((time.time() - start)/(c+calls*(toy))*((nr_of_toys-toy)*calls-c)))))
  940.  
  941. # with open("data/zfit_toys/toy_{0}/{1}.pkl".format(toy, call), "wb") as f:
  942. # pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)
  943.  
  944.  
  945. # In[ ]:
  946.  
  947.  
  948. # with open(r"data/zfit_toys/toy_0/0.pkl", "rb") as input_file:
  949. # sam = pkl.load(input_file)
  950. # print(sam[:10])
  951.  
  952. # with open(r"data/zfit_toys/toy_0/1.pkl", "rb") as input_file:
  953. # sam2 = pkl.load(input_file)
  954. # print(sam2[:10])
  955.  
  956. # print(np.sum(sam-sam2))
  957.  
  958.  
  959. # In[ ]:
  960.  
  961.  
  962. # print("Time to generate full toy: {} s".format(int(time.time()-start)))
  963.  
  964. # total_samp = []
  965.  
  966. # for call in range(calls):
  967. # with open(r"data/zfit_toys/toy_0/{}.pkl".format(call), "rb") as input_file:
  968. # sam = pkl.load(input_file)
  969. # total_samp = np.append(total_samp, sam)
  970.  
  971. # total_samp = total_samp.astype('float64')
  972.  
  973. # data2 = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)
  974.  
  975. # data3 = zfit.data.Data.from_numpy(array=total_samp, obs=obs)
  976.  
  977. # print(total_samp[:nevents].shape)
  978.  
  979.  
  980. # In[ ]:
  981.  
  982.  
  983. # plt.clf()
  984.  
  985. # bins = int((x_max-x_min)/7)
  986.  
  987. # # calcs = zfit.run(total_test_tf(samp))
  988. # print(total_samp[:nevents].shape)
  989.  
  990. # plt.hist(total_samp[:nevents], bins = bins, range = (x_min,x_max), label = 'data')
  991. # # plt.plot(test_q, calcs_test*nevents , label = 'pdf')
  992.  
  993. # # plt.plot(sam, calcs, '.')
  994. # # plt.plot(test_q, calcs_test)
  995. # # plt.yscale('log')
  996. # plt.ylim(0, 200)
  997. # # plt.xlim(3080, 3110)
  998.  
  999. # plt.legend()
  1000.  
  1001. # plt.savefig('test2.png')
  1002.  
  1003.  
  1004. # In[ ]:
  1005.  
  1006.  
  1007. # sampler = total_f.create_sampler(n=nevents)
  1008. # nll = zfit.loss.UnbinnedNLL(model=total_f, data=sampler, fit_range = (x_min, x_max))
  1009.  
  1010. # # for param in pdf.get_dependents():
  1011. # # param.set_value(initial_value)
  1012.  
  1013. # sampler.resample(n=nevents)
  1014.  
  1015. # # Randomise initial values
  1016. # # for param in pdf.get_dependents():
  1017. # # param.set_value(random value here)
  1018.  
  1019. # # Minimise the NLL
  1020. # minimizer = zfit.minimize.MinuitMinimizer(verbosity = 10)
  1021. # minimum = minimizer.minimize(nll)
  1022.  
  1023.  
  1024. # In[ ]:
  1025.  
  1026.  
  1027. # jpsi_width
  1028.  
  1029.  
  1030. # In[ ]:
  1031.  
  1032.  
  1033. # plt.hist(sample, weights=1 / prob(sample))
  1034.  
  1035.  
  1036. # # Fitting
  1037.  
  1038. # In[ ]:
  1039.  
  1040.  
  1041. # start = time.time()
  1042.  
  1043. # for param in total_f.get_dependents():
  1044. # param.randomize()
  1045. # # for param in total_f.get_dependents():
  1046. # # print(zfit.run(param))
  1047. # nll = zfit.loss.UnbinnedNLL(model=total_f, data=data2, fit_range = (x_min, x_max))
  1048.  
  1049. # minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)
  1050. # # minimizer._use_tfgrad = False
  1051. # result = minimizer.minimize(nll)
  1052.  
  1053. # # param_errors = result.error()
  1054.  
  1055. # # for var, errors in param_errors.items():
  1056. # # print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower']))
  1057.  
  1058. # print("Function minimum:", result.fmin)
  1059. # # print("Results:", result.params)
  1060. # print("Hesse errors:", result.hesse())
  1061.  
  1062.  
  1063. # In[ ]:
  1064.  
  1065.  
  1066. # print("Time taken for fitting: {}".format(display_time(int(time.time()-start))))
  1067.  
  1068. # # probs = total_f.pdf(test_q)
  1069.  
  1070. # calcs_test = zfit.run(probs)
  1071. # res_y = zfit.run(jpsi_res(test_q))
  1072.  
  1073.  
  1074. # In[ ]:
  1075.  
  1076.  
  1077. # plt.clf()
  1078. # # plt.plot(x_part, calcs, '.')
  1079. # plt.plot(test_q, calcs_test, label = 'pdf')
  1080. # # plt.plot(test_q, res_y, label = 'res')
  1081. # plt.legend()
  1082. # plt.ylim(0.0, 10e-6)
  1083. # # plt.yscale('log')
  1084. # # plt.xlim(3080, 3110)
  1085. # plt.savefig('test3.png')
  1086. # # print(jpsi_width)
  1087.  
  1088.  
  1089. # In[ ]:
  1090.  
  1091.  
  1092. # _tot = 4.37e-7+6.02e-5+4.97e-6
  1093. # _probs = []
  1094. # _probs.append(6.02e-5/_tot)
  1095. # _probs.append(4.97e-6/_tot)
  1096. # _probs.append(4.37e-7/_tot)
  1097. # print(_probs)
  1098.  
  1099.  
  1100. # In[ ]:
  1101.  
  1102.  
  1103. # dtype = 'float64'
  1104. # # mixture = tfd.Uniform(tf.constant(x_min, dtype=dtype), tf.constant(x_max, dtype=dtype))
  1105. # mixture = tfd.MixtureSameFamily(mixture_distribution=tfd.Categorical(probs=[tf.constant(0.007, dtype=dtype),
  1106. # tf.constant(0.917, dtype=dtype),
  1107. # tf.constant(0.076, dtype=dtype)]),
  1108. # components_distribution=tfd.Uniform(low=[tf.constant(x_min, dtype=dtype),
  1109. # tf.constant(3080, dtype=dtype),
  1110. # tf.constant(3670, dtype=dtype)],
  1111. # high=[tf.constant(x_max, dtype=dtype),
  1112. # tf.constant(3112, dtype=dtype),
  1113. # tf.constant(3702, dtype=dtype)]))
  1114. # # for i in range(10):
  1115. # # print(zfit.run(mixture.prob(mixture.sample((10, 1)))))
  1116.  
  1117.  
  1118. # In[ ]:
  1119.  
  1120.  
  1121. # print((zfit.run(jpsi_p)%(2*np.pi))/np.pi)
  1122. # print((zfit.run(psi2s_p)%(2*np.pi))/np.pi)
  1123.  
  1124.  
  1125. # In[ ]:
  1126.  
  1127.  
  1128. # def jpsi_res(q):
  1129. # return resonance(q, _mass = jpsi_mass, scale = jpsi_scale,
  1130. # phase = jpsi_phase, width = jpsi_width)
  1131.  
  1132. # def psi2s_res(q):
  1133. # return resonance(q, _mass = psi2s_mass, scale = psi2s_scale,
  1134. # phase = psi2s_phase, width = psi2s_width)
  1135. # def p3770_res(q):
  1136. # return resonance(q, _mass = p3770_mass, scale = p3770_scale,
  1137. # phase = p3770_phase, width = p3770_width)
  1138. # def p4040_res(q):
  1139. # return resonance(q, _mass = p4040_mass, scale = p4040_scale,
  1140. # phase = p4040_phase, width = p4040_width)
  1141. # def p4160_res(q):
  1142. # return resonance(q, _mass = p4160_mass, scale = p4160_scale,
  1143. # phase = p4160_phase, width = p4160_width)
  1144. # def p4415_res(q):
  1145. # return resonance(q, _mass = p4415_mass, scale = p4415_scale,
  1146. # phase = p4415_phase, width = p4415_width)
  1147.  
  1148.  
  1149. # In[ ]:
  1150.  
  1151.  
  1152. # 0.15**2*4.2/1000
  1153. # result.hesse()
  1154.  
  1155.  
  1156. # ## Constraints
  1157.  
  1158. # In[ ]:
  1159.  
  1160.  
  1161. # 1. Constraint - Real part of sum of Psi contrib and D contribs
  1162.  
  1163. sum_list = []
  1164.  
  1165. sum_list.append(ztf.to_complex(jpsi_s) * tf.exp(tf.complex(ztf.constant(0.0), jpsi_p)) * ztf.to_complex(jpsi_w / (tf.pow(jpsi_m,3))))
  1166. sum_list.append(ztf.to_complex(psi2s_s) * tf.exp(tf.complex(ztf.constant(0.0), psi2s_p)) * ztf.to_complex(psi2s_w / (tf.pow(psi2s_m,3))))
  1167. sum_list.append(ztf.to_complex(p3770_s) * tf.exp(tf.complex(ztf.constant(0.0), p3770_p)) * ztf.to_complex(p3770_w / (tf.pow(p3770_m,3))))
  1168. sum_list.append(ztf.to_complex(p4040_s) * tf.exp(tf.complex(ztf.constant(0.0), p4040_p)) * ztf.to_complex(p4040_w / (tf.pow(p4040_m,3))))
  1169. sum_list.append(ztf.to_complex(p4160_s) * tf.exp(tf.complex(ztf.constant(0.0), p4160_p)) * ztf.to_complex(p4160_w / (tf.pow(p4160_m,3))))
  1170. sum_list.append(ztf.to_complex(p4415_s) * tf.exp(tf.complex(ztf.constant(0.0), p4415_p)) * ztf.to_complex(p4415_w / (tf.pow(p4415_m,3))))
  1171. sum_list.append(ztf.to_complex(DDstar_s) * tf.exp(tf.complex(ztf.constant(0.0), DDstar_p)) * (ztf.to_complex(1.0 / (10.0*tf.pow(Dstar_m,2)) + 1.0 / (10.0*tf.pow(D_m,2)))))
  1172. sum_list.append(ztf.to_complex(Dbar_s) * tf.exp(tf.complex(ztf.constant(0.0), Dbar_p)) * ztf.to_complex(1.0 / (6.0*tf.pow(Dbar_m,2))))
  1173.  
  1174. sum_ru_1 = ztf.to_complex(ztf.constant(0.0))
  1175.  
  1176. for part in sum_list:
  1177. sum_ru_1 += part
  1178.  
  1179. sum_1 = tf.math.real(sum_ru_1)
  1180. # constraint1 = zfit.constraint.GaussianConstraint(params = sum_1, mu = ztf.constant(1.7*10**-8),
  1181. # sigma = ztf.constant(2.2*10**-8))
  1182.  
  1183. constraint1 = tf.pow((sum_1-ztf.constant(1.7*10**-8))/ztf.constant(2.2*10**-8),2)/ztf.constant(2.)
  1184.  
  1185. # 2. Constraint - Abs. of sum of Psi contribs and D contribs
  1186.  
  1187. sum_2 = tf.abs(sum_ru_1)
  1188. constraint2 = tf.cond(tf.greater_equal(sum_2, 5.0e-8), lambda: 100000., lambda: 0.)
  1189.  
  1190. # 3. Constraint - Maximum eta of D contribs
  1191.  
  1192. constraint3_0 = tf.cond(tf.greater_equal(tf.abs(Dbar_s), 0.2), lambda: 100000., lambda: 0.)
  1193.  
  1194. constraint3_1 = tf.cond(tf.greater_equal(tf.abs(DDstar_s), 0.2), lambda: 100000., lambda: 0.)
  1195.  
  1196. # 4. Constraint - Formfactor multivariant gaussian covariance fplus
  1197.  
  1198. Cov_matrix = [[ztf.constant( 1.), ztf.constant( 0.45), ztf.constant( 0.19), ztf.constant(0.857), ztf.constant(0.598), ztf.constant(0.531), ztf.constant(0.752), ztf.constant(0.229), ztf.constant(0,117)],
  1199. [ztf.constant( 0.45), ztf.constant( 1.), ztf.constant(0.677), ztf.constant(0.708), ztf.constant(0.958), ztf.constant(0.927), ztf.constant(0.227), ztf.constant(0.443), ztf.constant(0.287)],
  1200. [ztf.constant( 0.19), ztf.constant(0.677), ztf.constant( 1.), ztf.constant(0.595), ztf.constant(0.770), ztf.constant(0.819),ztf.constant(-0.023), ztf.constant( 0.07), ztf.constant(0.196)],
  1201. [ztf.constant(0.857), ztf.constant(0.708), ztf.constant(0.595), ztf.constant( 1.), ztf.constant( 0.83), ztf.constant(0.766), ztf.constant(0.582), ztf.constant(0.237), ztf.constant(0.192)],
  1202. [ztf.constant(0.598), ztf.constant(0.958), ztf.constant(0.770), ztf.constant( 0.83), ztf.constant( 1.), ztf.constant(0.973), ztf.constant(0.324), ztf.constant(0.372), ztf.constant(0.272)],
  1203. [ztf.constant(0.531), ztf.constant(0.927), ztf.constant(0.819), ztf.constant(0.766), ztf.constant(0.973), ztf.constant( 1.), ztf.constant(0.268), ztf.constant(0.332), ztf.constant(0.269)],
  1204. [ztf.constant(0.752), ztf.constant(0.227),ztf.constant(-0.023), ztf.constant(0.582), ztf.constant(0.324), ztf.constant(0.268), ztf.constant( 1.), ztf.constant( 0.59), ztf.constant(0.515)],
  1205. [ztf.constant(0.229), ztf.constant(0.443), ztf.constant( 0.07), ztf.constant(0.237), ztf.constant(0.372), ztf.constant(0.332), ztf.constant( 0.59), ztf.constant( 1.), ztf.constant(0.897)],
  1206. [ztf.constant(0.117), ztf.constant(0.287), ztf.constant(0.196), ztf.constant(0.192), ztf.constant(0.272), ztf.constant(0.269), ztf.constant(0.515), ztf.constant(0.897), ztf.constant( 1.)]]
  1207.  
  1208. def triGauss(val1,val2,val3,m = Cov_matrix):
  1209.  
  1210. mean1 = ztf.constant(0.466)
  1211. mean2 = ztf.constant(-0.885)
  1212. mean3 = ztf.constant(-0.213)
  1213. sigma1 = ztf.constant(0.014)
  1214. sigma2 = ztf.constant(0.128)
  1215. sigma3 = ztf.constant(0.548)
  1216. x1 = (val1-mean1)/sigma1
  1217. x2 = (val2-mean2)/sigma2
  1218. x3 = (val3-mean3)/sigma3
  1219. rho12 = m[0][1]
  1220. rho13 = m[0][2]
  1221. rho23 = m[1][2]
  1222. w = x1*x1*(rho23*rho23-1) + x2*x2*(rho13*rho13-1)+x3*x3*(rho12*rho12-1)+2*(x1*x2*(rho12-rho13*rho23)+x1*x3*(rho13-rho12*rho23)+x2*x3*(rho23-rho12*rho13))
  1223. d = 2*(rho12*rho12+rho13*rho13+rho23*rho23-2*rho12*rho13*rho23-1)
  1224. fcn = -w/d
  1225. chisq = -2*fcn
  1226. return chisq
  1227.  
  1228. constraint4 = triGauss(bplus_0, bplus_1, bplus_2)
  1229.  
  1230. # mean1 = ztf.constant(0.466)
  1231. # mean2 = ztf.constant(-0.885)
  1232. # mean3 = ztf.constant(-0.213)
  1233. # sigma1 = ztf.constant(0.014)
  1234. # sigma2 = ztf.constant(0.128)
  1235. # sigma3 = ztf.constant(0.548)
  1236. # constraint4_0 = tf.pow((bplus_0-mean1)/sigma1,2)/ztf.constant(2.)
  1237. # constraint4_1 = tf.pow((bplus_1-mean2)/sigma2,2)/ztf.constant(2.)
  1238. # constraint4_2 = tf.pow((bplus_2-mean3)/sigma3,2)/ztf.constant(2.)
  1239.  
  1240. # 5. Constraint - Abs. of sum of light contribs
  1241.  
  1242. sum_list_5 = []
  1243.  
  1244. sum_list_5.append(rho_s*rho_w/rho_m)
  1245. sum_list_5.append(omega_s*omega_w/omega_m)
  1246. sum_list_5.append(phi_s*phi_w/phi_m)
  1247.  
  1248.  
  1249. sum_ru_5 = ztf.constant(0.0)
  1250.  
  1251. for part in sum_list_5:
  1252. sum_ru_5 += part
  1253.  
  1254. constraint5 = tf.cond(tf.greater_equal(tf.abs(sum_ru_5), ztf.constant(0.02)), lambda: 100000., lambda: 0.)
  1255.  
  1256. # 6. Constraint on phases of Jpsi and Psi2s for cut out fit
  1257.  
  1258.  
  1259. # constraint6_0 = zfit.constraint.GaussianConstraint(params = jpsi_p, mu = ztf.constant(pdg["jpsi_phase_unc"]),
  1260. # sigma = ztf.constant(jpsi_phase))
  1261. # constraint6_1 = zfit.constraint.GaussianConstraint(params = psi2s_p, mu = ztf.constant(pdg["psi2s_phase_unc"]),
  1262. # sigma = ztf.constant(psi2s_phase))
  1263.  
  1264. constraint6_0 = tf.pow((jpsi_p-ztf.constant(jpsi_phase))/ztf.constant(pdg["jpsi_phase_unc"]),2)/ztf.constant(2.)
  1265. constraint6_1 = tf.pow((psi2s_p-ztf.constant(psi2s_phase))/ztf.constant(pdg["psi2s_phase_unc"]),2)/ztf.constant(2.)
  1266.  
  1267. # 7. Constraint on Ctt with higher limits
  1268.  
  1269. constraint7 = tf.cond(tf.greater_equal(Ctt*Ctt, 0.25), lambda: 100000., lambda: 0.)
  1270.  
  1271. constraint7dtype = tf.float64
  1272.  
  1273. # zfit.run(constraint6_0)
  1274.  
  1275. # ztf.convert_to_tensor(constraint6_0)
  1276.  
  1277. #List of all constraints
  1278.  
  1279. constraints = [constraint1, constraint2, constraint3_0, constraint3_1,# constraint4, #constraint4_0, constraint4_1, constraint4_2,
  1280. constraint6_0, constraint6_1]#, constraint7]
  1281.  
  1282.  
  1283. # ## Reset params
  1284.  
  1285. # In[ ]:
  1286.  
  1287.  
  1288. param_values_dic = {
  1289. 'jpsi_m': jpsi_mass,
  1290. 'jpsi_s': jpsi_scale,
  1291. 'jpsi_p': jpsi_phase,
  1292. 'jpsi_w': jpsi_width,
  1293. 'psi2s_m': psi2s_mass,
  1294. 'psi2s_s': psi2s_scale,
  1295. 'psi2s_p': psi2s_phase,
  1296. 'psi2s_w': psi2s_width,
  1297. 'p3770_m': p3770_mass,
  1298. 'p3770_s': p3770_scale,
  1299. 'p3770_p': p3770_phase,
  1300. 'p3770_w': p3770_width,
  1301. 'p4040_m': p4040_mass,
  1302. 'p4040_s': p4040_scale,
  1303. 'p4040_p': p4040_phase,
  1304. 'p4040_w': p4040_width,
  1305. 'p4160_m': p4160_mass,
  1306. 'p4160_s': p4160_scale,
  1307. 'p4160_p': p4160_phase,
  1308. 'p4160_w': p4160_width,
  1309. 'p4415_m': p4415_mass,
  1310. 'p4415_s': p4415_scale,
  1311. 'p4415_p': p4415_phase,
  1312. 'p4415_w': p4415_width,
  1313. 'rho_m': rho_mass,
  1314. 'rho_s': rho_scale,
  1315. 'rho_p': rho_phase,
  1316. 'rho_w': rho_width,
  1317. 'omega_m': omega_mass,
  1318. 'omega_s': omega_scale,
  1319. 'omega_p': omega_phase,
  1320. 'omega_w': omega_width,
  1321. 'phi_m': phi_mass,
  1322. 'phi_s': phi_scale,
  1323. 'phi_p': phi_phase,
  1324. 'phi_w': phi_width,
  1325. 'Dstar_m': Dstar_mass,
  1326. 'DDstar_s': 0.0,
  1327. 'DDstar_p': 0.0,
  1328. 'D_m': D_mass,
  1329. 'Dbar_m': Dbar_mass,
  1330. 'Dbar_s': 0.0,
  1331. 'Dbar_p': 0.0,
  1332. 'tau_m': pdg['tau_M'],
  1333. 'Ctt': 0.0,
  1334. 'b0_0': 0.292,
  1335. 'b0_1': 0.281,
  1336. 'b0_2': 0.150,
  1337. 'bplus_0': 0.466,
  1338. 'bplus_1': -0.885,
  1339. 'bplus_2': -0.213,
  1340. 'bT_0': 0.460,
  1341. 'bT_1': -1.089,
  1342. 'bT_2': -1.114}
  1343.  
  1344.  
  1345. def reset_param_values(variation = 0.05):
  1346. for param in total_f_fit.get_dependents():
  1347. if param.floating:
  1348. param.set_value(param_values_dic[param.name] + random.uniform(-1, 1) * param_values_dic[param.name]* variation)
  1349. # print(param.name)
  1350. # for param in totalf.get_dependents():
  1351. # param.set_value()
  1352.  
  1353.  
  1354. # # Analysis
  1355.  
  1356. # In[ ]:
  1357.  
  1358.  
  1359. # # zfit.run.numeric_checks = False
  1360.  
  1361. # fitting_range = 'cut'
  1362. # total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7
  1363. # cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR
  1364.  
  1365. # Ctt_list = []
  1366. # Ctt_error_list = []
  1367.  
  1368. # nr_of_toys = 1
  1369. # if fitting_range == 'cut':
  1370. # nevents = int(pdg["number_of_decays"]*cut_BR)
  1371. # else:
  1372. # nevents = int(pdg["number_of_decays"])
  1373. # # nevents = pdg["number_of_decays"]
  1374. # event_stack = 1000000
  1375. # # nevents *= 41
  1376. # # zfit.settings.set_verbosity(10)
  1377. # calls = int(nevents/event_stack + 1)
  1378.  
  1379. # total_samp = []
  1380.  
  1381. # start = time.time()
  1382.  
  1383. # sampler = total_f.create_sampler(n=event_stack)
  1384.  
  1385. # for toy in range(nr_of_toys):
  1386. # ### Generate data
  1387. # # clear_output(wait=True)
  1388. # print("Toy {}: Generating data...".format(toy))
  1389. # dirName = 'data/zfit_toys/toy_{0}'.format(toy)
  1390. # if not os.path.exists(dirName):
  1391. # os.mkdir(dirName)
  1392. # print("Directory " , dirName , " Created ")
  1393. # reset_param_values()
  1394. # if fitting_range == 'cut':
  1395. # sampler.resample(n=nevents)
  1396. # s = sampler.unstack_x()
  1397. # sam = zfit.run(s)
  1398. # calls = 0
  1399. # c = 1
  1400. # else:
  1401. # for call in range(calls):
  1402.  
  1403. # sampler.resample(n=event_stack)
  1404. # s = sampler.unstack_x()
  1405. # sam = zfit.run(s)
  1406.  
  1407. # c = call + 1
  1408.  
  1409. # with open("data/zfit_toys/toy_{0}/{1}.pkl".format(toy, call), "wb") as f:
  1410. # pkl.dump(sam, f, pkl.HIGHEST_PROTOCOL)
  1411. # print("Toy {}: Data generation finished".format(toy))
  1412. # ### Load data
  1413. # print("Toy {}: Loading data...".format(toy))
  1414. # if fitting_range == 'cut':
  1415. # total_samp = sam
  1416. # else:
  1417. # for call in range(calls):
  1418. # with open(r"data/zfit_toys/toy_0/{}.pkl".format(call), "rb") as input_file:
  1419. # sam = pkl.load(input_file)
  1420. # total_samp = np.append(total_samp, sam)
  1421.  
  1422. # total_samp = total_samp.astype('float64')
  1423. # if fitting_range == 'full':
  1424.  
  1425. # data = zfit.data.Data.from_numpy(array=total_samp[:int(nevents)], obs=obs)
  1426. # print("Toy {}: Loading data finished".format(toy))
  1427.  
  1428. # ### Fit data
  1429.  
  1430. # print("Toy {}: Fitting pdf...".format(toy))
  1431.  
  1432. # for param in total_f.get_dependents():
  1433. # param.randomize()
  1434.  
  1435. # nll = zfit.loss.UnbinnedNLL(model=total_f, data=data, fit_range = (x_min, x_max), constraints = constraints)
  1436.  
  1437. # minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)
  1438. # # minimizer._use_tfgrad = False
  1439. # result = minimizer.minimize(nll)
  1440.  
  1441. # print("Toy {}: Fitting finished".format(toy))
  1442.  
  1443. # print("Function minimum:", result.fmin)
  1444. # print("Hesse errors:", result.hesse())
  1445.  
  1446. # params = result.params
  1447. # Ctt_list.append(params[Ctt]['value'])
  1448. # Ctt_error_list.append(params[Ctt]['minuit_hesse']['error'])
  1449.  
  1450. # #plotting the result
  1451.  
  1452. # plotdirName = 'data/plots'.format(toy)
  1453.  
  1454. # if not os.path.exists(plotdirName):
  1455. # os.mkdir(plotdirName)
  1456. # # print("Directory " , dirName , " Created ")
  1457. # probs = total_f.pdf(test_q, norm_range=False)
  1458. # calcs_test = zfit.run(probs)
  1459. # plt.clf()
  1460. # plt.plot(test_q, calcs_test, label = 'pdf')
  1461. # plt.legend()
  1462. # plt.ylim(0.0, 6e-6)
  1463. # plt.savefig(plotdirName + '/toy_fit_full_range{}.png'.format(toy))
  1464.  
  1465. # print("Toy {0}/{1}".format(toy+1, nr_of_toys))
  1466. # print("Time taken: {}".format(display_time(int(time.time() - start))))
  1467. # print("Projected time left: {}".format(display_time(int((time.time() - start)/(c+calls*(toy))*((nr_of_toys-toy)*calls-c)))))
  1468. # if fitting_range == 'cut':
  1469. # _1 = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 60.)))
  1470. # tot_sam_1 = total_samp[_1]
  1471. # _2 = np.where((total_samp >= (jpsi_mass + 70.)) & (total_samp <= (psi2s_mass - 50.)))
  1472. # tot_sam_2 = total_samp[_2]
  1473.  
  1474. # _3 = np.where((total_samp >= (psi2s_mass + 50.)) & (total_samp <= x_max))
  1475. # tot_sam_3 = total_samp[_3]
  1476.  
  1477. # tot_sam = np.append(tot_sam_1, tot_sam_2)
  1478. # tot_sam = np.append(tot_sam, tot_sam_3)
  1479. # data = zfit.data.Data.from_numpy(array=tot_sam[:int(nevents)], obs=obs_fit)
  1480. # print("Toy {}: Loading data finished".format(toy))
  1481. # ### Fit data
  1482.  
  1483. # print("Toy {}: Fitting pdf...".format(toy))
  1484.  
  1485. # for param in total_f_fit.get_dependents():
  1486. # param.randomize()
  1487.  
  1488. # nll = zfit.loss.UnbinnedNLL(model=total_f_fit, data=data, constraints = constraints)
  1489.  
  1490. # minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)
  1491. # # minimizer._use_tfgrad = False
  1492. # result = minimizer.minimize(nll)
  1493.  
  1494. # print("Function minimum:", result.fmin)
  1495. # print("Hesse errors:", result.hesse())
  1496.  
  1497. # params = result.params
  1498. # if result.converged:
  1499. # Ctt_list.append(params[Ctt]['value'])
  1500. # Ctt_error_list.append(params[Ctt]['minuit_hesse']['error'])
  1501.  
  1502. # #plotting the result
  1503.  
  1504. # plotdirName = 'data/plots'.format(toy)
  1505.  
  1506. # if not os.path.exists(plotdirName):
  1507. # os.mkdir(plotdirName)
  1508. # # print("Directory " , dirName , " Created ")
  1509. # plt.clf()
  1510. # plt.hist(tot_sam, bins = int((x_max-x_min)/7.), label = 'toy data')
  1511. # plt.savefig(plotdirName + '/toy_histo_cut_region{}.png'.format(toy))
  1512.  
  1513. # probs = total_f_fit.pdf(test_q, norm_range=False)
  1514. # calcs_test = zfit.run(probs)
  1515. # plt.clf()
  1516. # plt.plot(test_q, calcs_test, label = 'pdf')
  1517. # plt.axvline(x=jpsi_mass-60.,color='red', linewidth=0.7, linestyle = 'dotted')
  1518. # plt.axvline(x=jpsi_mass+70.,color='red', linewidth=0.7, linestyle = 'dotted')
  1519. # plt.axvline(x=psi2s_mass-50.,color='red', linewidth=0.7, linestyle = 'dotted')
  1520. # plt.axvline(x=psi2s_mass+50.,color='red', linewidth=0.7, linestyle = 'dotted')
  1521. # plt.legend()
  1522. # plt.ylim(0.0, 1.5e-6)
  1523. # plt.savefig(plotdirName + '/toy_fit_cut_region{}.png'.format(toy))
  1524. # print("Toy {0}/{1}".format(toy+1, nr_of_toys))
  1525. # print("Time taken: {}".format(display_time(int(time.time() - start))))
  1526. # print("Projected time left: {}".format(display_time(int((time.time() - start)/(toy+1))*((nr_of_toys-toy-1)))))
  1527.  
  1528.  
  1529. # In[ ]:
  1530.  
  1531.  
  1532. # with open("data/results/Ctt_list.pkl", "wb") as f:
  1533. # pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)
  1534. # with open("data/results/Ctt_error_list.pkl", "wb") as f:
  1535. # pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)
  1536.  
  1537.  
  1538. # In[ ]:
  1539.  
  1540.  
  1541. # print('{0}/{1} fits converged'.format(len(Ctt_list), nr_of_toys))
  1542. # print('Mean Ctt value = {}'.format(np.mean(Ctt_list)))
  1543. # print('Mean Ctt error = {}'.format(np.mean(Ctt_error_list)))
  1544. # print('95 Sensitivy = {}'.format(((2*np.mean(Ctt_error_list))**2)*4.2/1000))
  1545.  
  1546.  
  1547. # In[ ]:
  1548.  
  1549.  
  1550. # plt.hist(tot_sam, bins = int((x_max-x_min)/7.))
  1551.  
  1552. # plt.show()
  1553.  
  1554. # # _ = np.where((total_samp >= x_min) & (total_samp <= (jpsi_mass - 50.)))
  1555.  
  1556. # tot_sam.shape
  1557.  
  1558.  
  1559. # In[ ]:
  1560.  
  1561.  
  1562. # Ctt.floating = False
  1563.  
  1564.  
  1565. # In[ ]:
  1566.  
  1567.  
  1568. # zfit.run(nll.value())
  1569.  
  1570.  
  1571. # In[ ]:
  1572.  
  1573.  
  1574. # result.fmin
  1575.  
  1576.  
  1577. # In[ ]:
  1578.  
  1579.  
  1580. # BR_steps = np.linspace(0.0, 1e-3, 11)
  1581. pull_dic = {}
  1582.  
  1583. mi = 1e-4
  1584. ma = 3e-3
  1585. ste = 20
  1586.  
  1587. for param in total_f_fit.get_dependents():
  1588. if param.floating:
  1589. pull_dic[param.name] = []
  1590. for step in range(2*ste):
  1591. pull_dic[param.name].append([])
  1592.  
  1593.  
  1594. def save_pulls(step):
  1595. for param in total_f_fit.get_dependents():
  1596. if param.floating:
  1597. pull_dic[param.name][step].append((params[param]['value'] - param_values_dic[param.name])/params[param]['minuit_hesse']['error'])
  1598.  
  1599.  
  1600.  
  1601. # for key in pull_dic.keys():
  1602. # print(np.shape(pull_dic[key]))
  1603. # save_pulls(New_step=True)
  1604. # params[Ctt]['value']
  1605.  
  1606.  
  1607. # In[ ]:
  1608.  
  1609.  
  1610. # for param in total_f_fit.get_dependents():
  1611. # if param.floating:
  1612. # print(param.name)
  1613.  
  1614. # print(params[Ctt])
  1615.  
  1616.  
  1617. # # CLS Code
  1618.  
  1619. # In[ ]:
  1620.  
  1621.  
  1622. # zfit.run.numeric_checks = False
  1623.  
  1624. load = False
  1625.  
  1626. bo = True
  1627.  
  1628. D_contribs = True
  1629.  
  1630. if not D_contribs:
  1631. Dbar_s.floating = False
  1632. Dbar_p.floating = False
  1633. DDstar_s.floating = False
  1634. DDstar_p.floating = False
  1635.  
  1636. bo_set = 1
  1637.  
  1638. fitting_range = 'cut'
  1639. total_BR = 1.7e-10 + 4.9e-10 + 2.5e-9 + 6.02e-5 + 4.97e-6 + 1.38e-9 + 4.2e-10 + 2.6e-9 + 6.1e-10 + 4.37e-7
  1640. cut_BR = 1.0 - (6.02e-5 + 4.97e-6)/total_BR
  1641.  
  1642. Ctt_list = []
  1643. Ctt_error_list = []
  1644.  
  1645. nr_of_toys = 1
  1646. nevents = int(pdg["number_of_decays"]*cut_BR)
  1647. # nevents = pdg["number_of_decays"]
  1648. event_stack = 1000000
  1649. # nevents *= 41
  1650. # zfit.settings.set_verbosity(10)
  1651.  
  1652. # mi = 1e-4
  1653. # ma = 3e-3
  1654. # ste = 13
  1655.  
  1656. BR_steps = np.linspace(mi, ma, ste)
  1657.  
  1658. BR_steps[0] = 0.0
  1659.  
  1660. print(BR_steps)
  1661.  
  1662. Ctt_steps = np.sqrt(BR_steps/4.2*1000)
  1663.  
  1664. print(Ctt_steps)
  1665.  
  1666. # total_samp = []
  1667.  
  1668. start = time.time()
  1669.  
  1670. Nll_list = []
  1671.  
  1672. sampler = total_f.create_sampler(n=nevents, fixed_params = False)
  1673. sampler.set_data_range(obs_fit)
  1674.  
  1675. __ = -1
  1676.  
  1677. #-----------------------------------------------------
  1678.  
  1679. if not load:
  1680. for Ctt_step in Ctt_steps:
  1681. __ += 1
  1682. for i in range(2):
  1683. Ctt_list.append([])
  1684. Ctt_error_list.append([])
  1685. Nll_list.append([])
  1686.  
  1687. for param in total_f_fit.get_dependents():
  1688. if param.floating:
  1689. pull_dic[param.name].append([])
  1690. for toy in range(nr_of_toys):
  1691. newset = True
  1692. while newset:
  1693. for floaty in [True, False]:
  1694. Ctt.floating = floaty
  1695. for bo_step in range(bo_set):
  1696.  
  1697. print('Step: {0}/{1}'.format(int(__), ste))
  1698. print('Current Ctt: {0}'.format(Ctt_step))
  1699. print('Ctt floating: {0}'.format(floaty))
  1700. reset_param_values(variation = 0.0)
  1701.  
  1702. if floaty:
  1703. print('Toy {0}/{1} - Fit {2}/{3}'.format(toy, nr_of_toys, bo_step, bo_set))
  1704. Ctt.set_value(Ctt_step)
  1705.  
  1706. else:
  1707. Ctt.set_value(0.0)
  1708. print('Toy {0}/{1} - Fit {2}/{3}'.format(toy, nr_of_toys, bo_step, bo_set))
  1709.  
  1710. if newset:
  1711. sampler.resample(n=nevents)
  1712. data = sampler
  1713. newset = False
  1714.  
  1715. ### Fit data
  1716. if floaty:
  1717. #plt.clf()
  1718. #plt.title('Ctt value: {:.2f}'.format(Ctt_step))
  1719. #plt.hist(zfit.run(data), bins = int((x_max-x_min)/7), range = (x_min, x_max))
  1720. #plt.savefig('data/CLs/plots/set_histo{}.png'.format(__))
  1721. _step = 2*__
  1722. else:
  1723. _step = 2*__+1
  1724.  
  1725. nll = zfit.loss.UnbinnedNLL(model=total_f_fit, data=data, constraints = constraints)
  1726.  
  1727. minimizer = zfit.minimize.MinuitMinimizer(verbosity = 5)
  1728. # minimizer._use_tfgrad = False
  1729. result = minimizer.minimize(nll)
  1730.  
  1731. print("Function minimum:", result.fmin)
  1732. print("Hesse errors:", result.hesse())
  1733.  
  1734. params = result.params
  1735.  
  1736. if result.converged:
  1737. save_pulls(step = _step)
  1738.  
  1739. if floaty:
  1740. Nll_list[-2].append(result.fmin)
  1741. Ctt_list[-2].append(params[Ctt]['value'])
  1742. Ctt_error_list[-2].append(params[Ctt]['minuit_hesse']['error'])
  1743.  
  1744. else:
  1745. Nll_list[-1].append(result.fmin)
  1746. Ctt_list[-1].append(0.0)
  1747. Ctt_error_list[-1].append(0.0)
  1748.  
  1749. else:
  1750. for _ in [1,2]:
  1751. del Nll_list[-_][toy*bo_set:]
  1752. # print(np.shape(Nll_list[-_]))
  1753. del Ctt_list[-_][toy*bo_set:]
  1754. del Ctt_error_list[-_][toy*bo_set:]
  1755. for param in total_f_fit.get_dependents():
  1756. if param.floating:
  1757. del pull_dic[param.name][_step+1-_][toy*bo_set:]
  1758. newset = True
  1759. break
  1760. if not result.converged:
  1761. break
  1762. print()
  1763. print('Time taken: {}'.format(display_time(int(time.time()-start))))
  1764. print('Estimated time left: {}'.format(display_time(int((time.time()-start)/(__+(toy+1)/nr_of_toys)*(ste-__-(nr_of_toys-toy-1)/nr_of_toys)))))
  1765.  
  1766.  
  1767. # In[ ]:
  1768.  
  1769.  
  1770. if load:
  1771. phase_combi = '--'
  1772. if D_contribs:
  1773. D_dir = 'D-True'
  1774. else:
  1775. D_dir = 'D-False'
  1776.  
  1777. _dir = 'data/CLs/finished/f1d1/{}/{}/'.format(phase_combi, D_dir)
  1778. jobs = os.listdir(_dir)
  1779. First = True
  1780. # print(jobs)
  1781. for job in jobs:
  1782. dirName = _dir + str(job) + '/data/CLs'
  1783. if not os.path.exists("{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys)):
  1784. # print(job)
  1785. continue
  1786. with open(r"{}/variab.pkl".format(dirName), "rb") as input_file:
  1787. variab = pkl.load(input_file)
  1788. # print(variab)
  1789. ### sanity check:
  1790. if variab['mi'] != mi or variab['ma'] != ma or variab['ste'] != ste or bo_set != bo_set:
  1791. print('Fitting parameters of data dont equal the ones given -- Job {} skipped!'.format(job))
  1792. with open(r"{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "rb") as input_file:
  1793. _Nll_list = pkl.load(input_file)
  1794. with open(r"{}/{}-{}_{}s{}b{}t--Ctt_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "rb") as input_file:
  1795. _Ctt_list = pkl.load(input_file)
  1796. with open(r"{}/{}-{}_{}s{}b{}t--Ctt_error_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "rb") as input_file:
  1797. _Ctt_error_list = pkl.load(input_file)
  1798. with open(r"{}/{}-{}_{}s{}b{}t--pull_dic.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "rb") as input_file:
  1799. _pull_dic = pkl.load(input_file)
  1800. with open(r"{}/{}-{}_{}s--CLs_list.pkl".format(dirName, mi,ma,ste), "rb") as input_file:
  1801. _CLs_list = pkl.load(input_file)
  1802. if First:
  1803. Nll_list = _Nll_list
  1804. Ctt_list = _Ctt_list
  1805. Ctt_error_list = _Ctt_error_list
  1806. pull_dic = _pull_dic
  1807. # print(_pull_dic)
  1808. CLs_list = _CLs_list
  1809. First = False
  1810. else:
  1811. for step in range(2*ste):
  1812. # print(Nll_list[step], step)
  1813. Nll_list[step].extend(_Nll_list[step])
  1814. Ctt_list[step].extend(_Ctt_list[step])
  1815. Ctt_error_list[step].extend(_Ctt_error_list[step])
  1816. for key in pull_dic.keys():
  1817. # print(key, np.shape(pull_dic[key]))
  1818. pull_dic[key][step].extend(_pull_dic[key][step])
  1819. for step in range(ste):
  1820. CLs_list[step].extend(_CLs_list[step])
  1821.  
  1822. # print('----------------------')
  1823.  
  1824.  
  1825. # In[ ]:
  1826.  
  1827.  
  1828. dirName = 'data/CLs'
  1829.  
  1830. # if bo and not load:
  1831. # for s in range(2*ste):
  1832. # Nll_list[s] = [np.min(Nll_list[s])]
  1833.  
  1834.  
  1835. if not load:
  1836. if not os.path.exists(dirName):
  1837. os.mkdir(dirName)
  1838. print("Directory " , dirName , " Created ")
  1839.  
  1840. with open("{}/{}-{}_{}s{}b{}t--CLs_Nll_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "wb") as f:
  1841. pkl.dump(Nll_list, f, pkl.HIGHEST_PROTOCOL)
  1842. with open("{}/{}-{}_{}s{}b{}t--Ctt_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "wb") as f:
  1843. pkl.dump(Ctt_list, f, pkl.HIGHEST_PROTOCOL)
  1844. with open("{}/{}-{}_{}s{}b{}t--Ctt_error_list.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "wb") as f:
  1845. pkl.dump(Ctt_error_list, f, pkl.HIGHEST_PROTOCOL)
  1846. with open("{}/{}-{}_{}s{}b{}t--pull_dic.pkl".format(dirName, mi,ma,ste,bo_set,nr_of_toys), "wb") as f:
  1847. pkl.dump(pull_dic, f, pkl.HIGHEST_PROTOCOL)
  1848. variab = {'mi': mi,
  1849. 'ma': ma,
  1850. 'ste': ste,
  1851. 'bo_set': bo_set,
  1852. 'nr_of_toys': nr_of_toys}
  1853. with open("{}/variab.pkl".format(dirName), "wb") as f:
  1854. pkl.dump(variab, f, pkl.HIGHEST_PROTOCOL)
  1855. CLs_values = []
  1856. toy_size = bo_set
  1857.  
  1858. print(np.shape(Nll_list))
  1859. print(Nll_list[0:1])
  1860. for step in range(ste):
  1861. CLs_values.append([])
  1862. for toy in range(nr_of_toys):
  1863. float_min = np.min(Nll_list[2*step][toy*bo_set:(toy+1)*bo_set])
  1864. fix_min = np.min(Nll_list[2*step+1][toy*bo_set:(toy+1)*bo_set])
  1865. CLs_values[step].append(float_min-fix_min)
  1866. print(np.shape(CLs_values))
  1867. with open("{}/{}-{}_{}s--CLs_list.pkl".format(dirName, mi,ma,ste), "wb") as f:
  1868. pkl.dump(CLs_values, f, pkl.HIGHEST_PROTOCOL)
  1869.  
  1870.  
  1871. # In[ ]:
  1872.  
  1873.  
  1874. # print(variab['mi'] != mi)
  1875.  
  1876.  
  1877. # ## Plot
  1878.  
  1879. # In[ ]:
  1880.  
  1881.  
  1882. l = []
  1883.  
  1884.  
  1885. if load:
  1886. CLs_values = -1*np.array(CLs_list)
  1887.  
  1888. if not os.path.exists('data/CLs/plots'):
  1889. os.mkdir('data/CLs/plots')
  1890. print("Directory " , 'data/CLs/plots' , " Created ")
  1891.  
  1892. print(np.shape(CLs_values))
  1893. for step in range(1,ste):
  1894. plt.clf()
  1895. plt.title('Ctt value: {:.2f}'.format(Ctt_steps[step]))
  1896. plt.hist(CLs_values[0], bins = 150, range = (-25, 50), label = 'Ctt fixed to 0')
  1897. plt.hist(CLs_values[step], bins = 150, range = (-25, 50), label = 'Ctt floating')
  1898. plt.axvline(x=np.mean(CLs_values[0]),color='red', linewidth=1.0, linestyle = 'dotted')
  1899. plt.legend()
  1900. plt.savefig('data/CLs/plots/CLs-BR({:.1E}).png'.format(BR_steps[step]))
  1901. l.append(len(np.where(np.array(CLs_values[step]) < np.mean(CLs_values[0]))[0]))
  1902. for step in range(2*ste):
  1903. if step%2 == 0:
  1904. floaty = True
  1905. else:
  1906. floaty = False
  1907. for key in pull_dic.keys():
  1908. if not os.path.exists('data/CLs/plots/{}'.format(key)):
  1909. os.mkdir('data/CLs/plots/{}'.format(key))
  1910. plt.clf()
  1911. plt.title('Pull {} - Ctt value {:.2f} - floating {}'.format(key, Ctt_steps[int(step/2)], floaty))
  1912. plt.hist(pull_dic[key][step], bins = 50, range = (-5,5))
  1913. plt.xlabel('Pull')
  1914. plt.savefig('data/CLs/plots/{}/{:.2f}Ctt{}s{}f.png'.format(key, Ctt_steps[int(step/2)], step, floaty))
  1915.  
  1916.  
  1917. # In[ ]:
  1918.  
  1919.  
  1920. for s in range(len(l)):
  1921. print('BR: {:.4f}'.format(BR_steps[s+1]))
  1922. print(2*l[s]/len(CLs_values[s]))
  1923. print()
  1924.  
  1925.  
  1926. # In[ ]:
  1927.  
  1928.  
  1929. # for step in range(2*ste):
  1930. # for key in pull_dic.keys():
  1931. # print(pull_dic[key][step])
  1932.  
  1933.  
  1934. # In[ ]:
  1935.  
  1936.  
  1937. # for param in total_f_fit.get_dependents():
  1938. # if param.floating:
  1939. # print(params[param]['value'])
  1940.  
  1941.  
  1942. # In[ ]:
  1943.  
  1944.  
  1945. print(display_time(int(time.time()-start)))
  1946.  
  1947.  
  1948. # In[ ]:
  1949.  
  1950.  
  1951. # variab['mi'] =! mi
  1952.  
  1953.  
  1954. # In[ ]:
  1955.  
  1956.  
  1957.  
  1958.  
  1959.  
  1960. # In[ ]:
  1961.  
  1962.  
  1963.  
  1964.  
  1965.  
  1966. # In[ ]:
  1967.  
  1968.  
  1969.  
  1970.  
  1971.  
  1972. # In[ ]:
  1973.  
  1974.  
  1975.  
  1976.  
  1977.  
  1978. # In[ ]:
  1979.  
  1980.  
  1981.  
  1982.  
  1983.  
  1984. # In[ ]:
  1985.  
  1986.  
  1987.  
  1988.  
  1989.  
  1990. # In[ ]:
  1991.  
  1992.  
  1993.  
  1994.  
  1995.  
  1996. # In[ ]:
  1997.  
  1998.  
  1999.  
  2000.  
  2001.  
  2002. # In[ ]:
  2003.  
  2004.  
  2005.  
  2006.  
  2007.  
  2008. # In[ ]:
  2009.  
  2010.  
  2011.  
  2012.