Skip to content
Snippets Groups Projects

Repackage: split pidml package into pidplots and pidml

Merged Repackage: split pidml package into pidplots and pidml
repackage into main
Merged Connor Hainje requested to merge repackage into main
Compare
and
24 files
+ 2384
− 3803
Compare changes
  • Side-by-side
  • Inline
Files
24
notebooks/DStar_reconstruction_validation.ipynb 0 → 100644
+ 333
− 0
%% Cell type:code id: tags:
``` python
import root_pandas as rtp
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import math
from matplotlib.colors import LogNorm
import matplotlib.patches as patches
from matplotlib.ticker import NullFormatter
import seaborn as sns
from ROOT import TProfile2D
```
%% Cell type:code id: tags:
``` python
plt.rc('figure', figsize=[12.0, 8.0])
SMALL_SIZE = 20
MEDIUM_SIZE = 25
BIGGER_SIZE = 35
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
```
%% Cell type:code id: tags:
``` python
mc = rtp.read_root(['/group/belle2/dataprod/Systematics/production/MC14ri_1/Dst/Dst_mc14p1_MC14ri_1_merged_sWeights.root',
'/group/belle2/dataprod/Systematics/production/MC14ri_2/Dst/Dst_mc14p2_MC14ri_2_merged_sWeights.root']
)
d1 = rtp.read_root([
# '/group/belle2/dataprod/Systematics/production/proc12e7/Dst/Dst_p12e7_proc12e7_merged_sWeights.root',
# '/group/belle2/dataprod/Systematics/production/proc12e8/Dst/Dst_p12e8_proc12e8_merged_sWeights.root',
# '/group/belle2/dataprod/Systematics/production/proc12e10/Dst/Dst_p12e10_proc12e10_merged_sWeights.root',
'/group/belle2/dataprod/Systematics/production/proc12e12/Dst/Dst_p12e12_proc12e12_merged_sWeights.root'
])
```
%% Cell type:markdown id: tags:
## Truth selection
%% Cell type:code id: tags:
``` python
def TruthSelection(df, particleList):
head, *tail = particleList
goodData = df[df['{}_isSignal'.format(head)] == 1]
for particle in tail:
goodData = goodData[goodData['{}_isSignal'.format(particle)] == 1]
return goodData
```
%% Cell type:code id: tags:
``` python
mcDstar = TruthSelection(mc, ['DST_D0_K', 'DST_D0_pi', 'DST_pi'])
```
%% Cell type:markdown id: tags:
# Reconstructing the decay $D^{*} \rightarrow (D^{0}\rightarrow K \pi) \pi$
%% Cell type:markdown id: tags:
## Impact Parameter
%% Cell type:code id: tags:
``` python
goodData = d1[abs(d1["DST_D0_K_d0"])<2]
goodData = goodData[abs(goodData["DST_D0_K_z0"])<4]
goodData = goodData[abs(goodData["DST_D0_pi_d0"])<2]
goodData = goodData[abs(goodData["DST_D0_pi_z0"])<4]
goodData = goodData[abs(goodData["DST_pi_d0"])<2]
goodData = goodData[abs(goodData["DST_pi_z0"])<4]
#goodData2 = d2[abs(d2["K_d0"])<2]
#goodData2 = goodData2[abs(goodData2["K_z0"])<4]
#goodData2 = goodData2[abs(goodData2["pi_d0"])<2]
#goodData2 = goodData2[abs(goodData2["pi_z0"])<4]
#goodData2 = goodData2[abs(goodData2["pi_slow_d0"])<2]
#goodData2 = goodData2[abs(goodData2["pi_slow_z0"])<4]
goodMC = mc[abs(mc["DST_D0_K_d0"])<2]
goodMC = goodMC[abs(goodMC["DST_D0_K_z0"])<4]
goodMC = goodMC[abs(goodMC["DST_D0_pi_d0"])<2]
goodMC = goodMC[abs(goodMC["DST_D0_pi_z0"])<4]
goodMC = goodMC[abs(goodMC["DST_pi_d0"])<2]
goodMC = goodMC[abs(goodMC["DST_pi_z0"])<4]
```
%% Cell type:markdown id: tags:
## Track reconstruction quality
%% Cell type:code id: tags:
``` python
plt.hist(goodData["DST_D0_K_chiProb"], bins=np.linspace(0,1,100), density=1, alpha=0.4, label="K, proc12");
plt.hist(mcDstar["DST_D0_K_chiProb"], bins=np.linspace(0,1,100), density=1, alpha=0.4, label="K, MC14");
plt.xlabel("$\chi^2$ probability")
plt.legend();
plt.savefig("KaonTrackChi2.pdf")
```
%% Cell type:code id: tags:
``` python
plt.hist(goodData["DST_D0_K_nCDCHits"], bins=np.linspace(20,80,61)-0.5, density=1, alpha=0.4, label="K, proc12");
plt.hist(mcDstar["DST_D0_K_nCDCHits"], bins=np.linspace(20,80,61)-0.5, density=1, alpha=0.4, label="K, MC14");
plt.xlabel("Number of CDC hits")
plt.legend();
plt.savefig('Dstar_reco_nCDCHits.pdf');
```
%% Cell type:raw id: tags:
# plt.hist(goodData["DST_D0_K_chiProb"], bins=np.linspace(0, 1, 50))
# plt.hist(goodData["DST_D0_K_nCDCHits"], bins=np.linspace(20, 60, 41))
goodData = goodData[goodData['DST_D0_K_chiProb']>0.1]
goodData = goodData[goodData['DST_D0_pi_chiProb']>0.1]
#goodData = goodData[goodData['D0_chiProb']>0.2]
#goodData2 = goodData2[goodData2['K_chiProb']>0.1]
#goodData2 = goodData2[goodData2['pi_chiProb']>0.1]
#goodData2 = goodData2[goodData2['D0_chiProb']>0.2]
goodMC = goodMC[goodMC['DST_D0_K_chiProb']>0.1]
goodMC = goodMC[goodMC['DST_D0_pi_chiProb']>0.1]
#goodMC = goodMC[goodMC['D0_chiProb']>0.2]
%% Cell type:code id: tags:
``` python
goodData = goodData[goodData['DST_D0_K_nCDCHits']>40]
goodData = goodData[goodData['DST_D0_pi_nCDCHits']>40]
#goodData = goodData[goodData['DST_pi_nCDCHits']>40]
#goodData2 = goodData2[goodData2['DST_D0_K_nCDCHits']>40]
#goodData2 = goodData2[goodData2['DST_D0_pi_nCDCHits']>40]
#goodData2 = goodData2[goodData2['DST_pi_nCDCHits']>40]
goodMC = goodMC[goodMC['DST_D0_K_nCDCHits']>40]
goodMC = goodMC[goodMC['DST_D0_pi_nCDCHits']>40]
#goodMC = goodMC[goodMC['pi_slow_nCDCHits']>40]
```
%% Cell type:markdown id: tags:
## Exclude bad runs
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
# More cuts
%% Cell type:code id: tags:
``` python
goodData = goodData[goodData['DST_D0_p']>3.0]
#goodData2 = goodData2[goodData2['D0_p']>3.0]
goodMC = goodMC[goodMC['DST_D0_p']>3.0]
goodData = goodData[goodData['M']-goodData['DST_D0_M']<0.15]
#goodData2 = goodData2[goodData2['Dstar_M']-goodData2['D0_M']<0.15]
goodMC = goodMC[goodMC['M']-goodMC['DST_D0_M']<0.15]
plt.hist(goodData["DST_D0_M"], bins=np.linspace(1.80,1.92,100))
plt.xlabel("$D^0$ mass (GeV)")
plt.axvspan(1.845, 1.880, alpha=0.2, color='gray')
plt.savefig('Dstar_reco_D0Mass.pdf');
```
%% Cell type:code id: tags:
``` python
h1 = goodData["DST_D0_M"]
y1, x1 = np.histogram(h1, bins=np.linspace(1.80,1.92,100))
weights1 = np.ones_like(h1)/y1.max()
#h2 = goodData2["D0_M"]
#y2, x2 = np.histogram(h2, bins=np.linspace(1.80,1.92,100))
#weights2 = np.ones_like(h2)/y2.max()
h3 = goodMC["DST_D0_M"]
y3, x3 = np.histogram(h3, bins=np.linspace(1.80,1.92,100))
weights3 = np.ones_like(h3)/y3.max()
```
%% Cell type:code id: tags:
``` python
plt.hist(h1, alpha=0.4, bins=np.linspace(1.80,1.92,100), label="proc12", weights=weights1)
#plt.hist(h2, alpha=0.4, bins=np.linspace(1.80,1.92,100), label="", weights=weights2)
plt.hist(h3, alpha=0.4, bins=np.linspace(1.80,1.92,100), label="MC14", weights=weights3)
plt.xlabel("$D^0$ mass (GeV)")
plt.axvspan(1.845, 1.880, alpha=0.2, color='gray')
plt.legend()
plt.savefig('Dstar_reco_D0Mass_Comparison.pdf');
```
%% Cell type:code id: tags:
``` python
goodData_side = goodData[(goodData['DST_D0_M']>1.80) & (goodData['DST_D0_M']<1.92)]
goodData_side = goodData_side[np.logical_or(1.84487>goodData_side['DST_D0_M'], goodData_side['DST_D0_M']>1.880)]
h_side = goodData_side["DST_D0_M"]
y_side, x_side = np.histogram(h_side, bins=np.linspace(1.80,1.92,100))
goodData_center = goodData[(goodData['DST_D0_M']>1.845) & (goodData['DST_D0_M']<1.880)]
x_sidec = []
y_sidec = []
for i in range(len(x_side)-1):
if (y_side[i]>0): x_sidec.append((x_side[i+1]+x_side[i])/2.0)
if (y_side[i]>0): y_sidec.append(y_side[i])
ffit = np.polyfit(x_sidec, y_sidec, 2)
plt.hist(goodData['DST_D0_M'], alpha=0.4, bins=np.linspace(1.80,1.92,100))
out = plt.hist(goodData_center['DST_D0_M'], alpha=0.4, bins=np.linspace(1.80,1.92,100))
plt.xlabel("$D^0$ mass (GeV)")
plt.axvspan(1.845, 1.880, alpha=0.2, color='gray')
plt.plot(x_sidec,y_sidec,'o')
trendpoly = np.poly1d(ffit)
plt.plot(x_sidec,trendpoly(x_sidec))
plt.savefig('Dstar_reco_D0Mass_Fit.pdf');
from scipy import integrate as TG
aoc, err = TG.quad(trendpoly, 1.845, 1.880)
print ("BKG: ", aoc)
yields = sum(out[0][1:100]*np.diff(out[1][1:100]))
print ("BKG + SIG: ", yields)
print ("SIG/BKG: ", (yields-aoc)/aoc)
```
%% Cell type:code id: tags:
``` python
goodData = goodData[(goodData['DST_D0_M']>1.845) & (goodData['DST_D0_M']<1.880)]
#goodData2 = goodData2[(goodData2['D0_M']>1.845) & (goodData2['D0_M']<1.880)]
goodMC = goodMC[(goodMC['DST_D0_M']>1.845) & (goodMC['DST_D0_M']<1.880)]
print(len(goodData), "Events are left in goodData")
#print(len(goodData2), "Events are left in goodData2")
print(len(goodMC), "Events are left in goodMC")
```
%% Cell type:code id: tags:
``` python
plt.hist(goodData["M"], bins=np.linspace(1.98,2.04,100));
plt.xlabel("$D^*$ mass (GeV)")
plt.savefig('Dstar_reco_DstarMass.pdf');
```
%% Cell type:code id: tags:
``` python
h1 = goodData["M"]
y1, x1 = np.histogram(h1, bins=np.linspace(1.98,2.04,100))
weights1 = np.ones_like(h1)/y1.max()
#h2 = goodData2["Dstar_M"]
#y2, x2 = np.histogram(h2, bins=np.linspace(1.98,2.04,100))
#weights2 = np.ones_like(h2)/y2.max()
h3 = goodMC["M"]
y3, x3 = np.histogram(h3, bins=np.linspace(1.98,2.04,100))
weights3 = np.ones_like(h3)/y3.max()
```
%% Cell type:code id: tags:
``` python
plt.hist(h1, alpha=0.4, bins=np.linspace(1.98,2.04,100), label="proc12", weights=weights1)
#plt.hist(h2, alpha=0.4, bins=np.linspace(1.98,2.04,100), label="", weights=weights2)
plt.hist(h3, alpha=0.4, bins=np.linspace(1.98,2.04,100), label="MC14", weights=weights3)
plt.xlabel("$D^*$ mass (GeV)")
plt.legend()
plt.savefig('Dstar_reco_DstMass_Comparison.pdf');
```
%% Cell type:markdown id: tags:
## Diagnostic plots
%% Cell type:code id: tags:
``` python
fig, axs = plt.subplots(1, 2, figsize=(20, 10), sharey=False)
axs[0].hist(goodData['DST_D0_K_d0'], bins=np.linspace(-0.3,0.3,100), alpha=0.4, label='K')
axs[0].hist(goodData['DST_D0_pi_d0'], bins=np.linspace(-0.3,0.3,100), alpha=0.4, label='$\pi$')
axs[0].hist(goodData['DST_pi_d0'], bins=np.linspace(-0.3,0.3,100), alpha=0.4, label='$\pi_s$')
axs[0].set_xlabel("$d_0$ (cm)")
axs[0].legend()
axs[1].hist(goodData['DST_D0_K_z0'], bins=np.linspace(-0.4,0.4,100), alpha=0.4, label='K')
axs[1].hist(goodData['DST_D0_pi_z0'], bins=np.linspace(-0.4,0.4,100), alpha=0.4, label='$\pi$')
axs[1].hist(goodData['DST_pi_z0'], bins=np.linspace(-0.4,0.4,100), alpha=0.4, label='$\pi_s$')
axs[1].set_xlabel("$z_0$ (cm)")
axs[1].legend();
plt.savefig('Dstar_reco_IP.pdf');
```
%% Cell type:code id: tags:
``` python
plt.hist(goodData["DST_D0_K_pt"], bins=np.linspace(0,5,100), alpha=0.4, density=1, label="K");
plt.hist(goodData["DST_D0_pi_pt"], bins=np.linspace(0,5,100), alpha=0.4, density=1, label="$\pi$");
plt.hist(goodData["DST_D0_pt"], bins=np.linspace(0,5,100), alpha=0.4, density=1, label="$D^0$");
plt.xlabel("$p_T$ (GeV)")
plt.legend()
plt.savefig('Dstar_reco_pt_Comparison.pdf');
```
%% Cell type:code id: tags:
``` python
```

imprint | privacy | Conditions of Use