In [1]:
from FutureColliderTools import *
from FutureColliderDataLoader import *
import numpy as np
import ROOT
import sys
sys.path.append('/home/ismith/ANA/FourVector')
from FourVector import FourVector
from ThreeVector import ThreeVector
ROOT.enableJSVis()
ROOT.gStyle.SetOptStat(0)
In [2]:
DataTuple = {}
DataTuple["13512010"] = ("root://eoslhcb.cern.ch//eos/lhcb/wg/semileptonic/Bs2KmunuAna/Tuples/Bs2KMuNu_MCTUPLES_RAW_08Feb17/DTT_13512010_Bs_Kmunu_DecProdCut_Up_Py8_MC12.root", "Bs2KMuNuMCTuple/MCDecayTree" )
DataTuple["13512010_combi"] = ("root://eoslhcb.cern.ch//eos/lhcb/wg/semileptonic/Bs2KmunuAna/Tuples/Bs2KMuNu_MCTUPLES_RAW_08Feb17/DTT_13512010_Bs_Kmunu_DecProdCut_Up_Py8_MC12.root", "Bs2KMuNuMCTuple/MCDecayTree" )
DataTuple["13774000_combi"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsMuNu/MCDecayTree")
DataTuple["11876001"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV11876001.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["12513010"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV12513010.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["15574010"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV15574010.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["13114001"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13114001.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["13114006"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13114006.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["13144032"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13144032.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["13146004"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13146004.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["12145052"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV12145052.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["22114001"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV22114001.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["23573002"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV23573002.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["13774011"] = ("/afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13774011.root", "MCDecayTreeTuple/MCDecayTree")
DataTuple["13774000_Ds"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsMuNu/MCDecayTree")
DataTuple["13774000_Dsstar"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsStarMuNu_g/MCDecayTree")
DataTuple["13774000_DsstarTau"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsStarTauNu_g/MCDecayTree")
DataTuple["13774000_DsTau"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsTauNu/MCDecayTree")
DataTuple["13774000_Dsstar0"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsStar0MuNu_pi0/MCDecayTree")
DataTuple["13774000_Dsstar1"] = ("root://eoslhcb.cern.ch//eos/lhcb/user/i/ismith/MCDecayTree/DTT_MC11_Bs2DsMuNu_13774000_Cocktail_SIGNAL.root", "TupleBs2DsStar12460MuNu_pi0/MCDecayTree")
In [3]:
sigma_PV_LHCb = np.array( [0.01, 0.01, 0.07])
sigma_SV_LHCb = np.array( [0.02, 0.02, 0.31])
In [4]:
import pickle
f_Histogram = ROOT.TFile.Open("Source_Histograms_KMu_{0}_LHCb.root".format(0.5), "RECREATE")
for EvType, FileTree in DataTuple.iteritems():
combinatorial = "combi" in EvType
PickleName = "/home/ismith/tmp/Bs_K_" + EvType + ".pickle"
InputData = None
try:
InputData,K_PE,Mu_PE,B_PE,B_Origin,B_End = LoadData_KMuNu(EvType, FileTree)
except:
continue
DataType = InputData[K_PE[0]].dtype
B_SV = InputData[B_Origin].copy().view(( DataType, 3 ))
B_PV = InputData[B_End].copy().view(( DataType, 3 ))
K = FourVector(InputData[K_PE])
Mu = FourVector(InputData[Mu_PE])
B = FourVector(InputData[B_PE])
if combinatorial:
B_SV = B_SV[1:-1]
B_PV = B_PV[0:-2]
K = K [1:-1]
Mu = Mu[0:-2]
B = B[0:-2]
print len( K), len( Mu )
B_SV_New = SmearVertex( B_SV, sigma_SV_LHCb * 0.5)
B_PV_New = SmearVertex( B_PV, sigma_PV_LHCb * 0.5)
B_Direction = ThreeVector(B_SV_New - B_PV_New)
Y = K + Mu
Cuts = ()
Cuts += ( K.P() > 10000, )
Cuts += ( K.Pt() > 500, )
Cuts += ( Mu.P() > 6000., )
Cuts += ( Mu.Pt() > 1500., )
Cuts += ( -np.cos( B_Direction.Angle( Y.Vect() ) ) > 0.997, )
PassCuts = np.all( Cuts, axis = 0)
"""
if EvType == "13512010":
print Cuts
print B_Direction.Angle(Y.Vect() )
print B_Direction.Dot( Y.Vect() )
"""
print EvType, 1.0* np.sum(PassCuts) / K.P().shape[0]
h_MCORR = ROOT.TH1F("MCORR_" + EvType, "", 100, 2500, 6000)
h_MM2 = ROOT.TH1F("MissingMass2_" + EvType, "", 100, -10e6, 20e6)
h_Q21 = ROOT.TH1F("QSQ_SOL1_" + EvType, "", 100, 0, 25e6)
h_Q22 = ROOT.TH1F("QSQ_SOL2_" + EvType, "", 100, 0, 25e6)
h_Q20 = ROOT.TH1F("QSQ_SOLTrue_" + EvType, "", 100, 0, 25e6)
QSQ = (B - K).M2()
MCORR = GetCorrectedMass(Y[PassCuts], B_Direction[PassCuts])
MM2 = GetMissingMass2(K[PassCuts], Mu[PassCuts], B_Direction[PassCuts])
Q21, Q22 = GetQ2(Y[PassCuts], Mu[PassCuts], B_Direction[PassCuts])
nev = MM2.shape[0]
#for MC, MM in zip(MCORR, MM2):
h_MCORR.FillN( nev, MCORR, np.ones(nev))
h_MM2 .FillN( nev, MM2, np.ones(nev))
h_Q21 .FillN( nev, Q21,np.ones(nev))
h_Q22 .FillN( nev, Q22,np.ones(nev))
h_Q20 .FillN( nev, QSQ,np.ones(nev))
h_MCORR.GetXaxis().SetTitle("m_{corr}(B_{s})")
h_MM2.GetXaxis().SetTitle("Missing Mass^{2}")
h_Q21.GetXaxis().SetTitle("q^{2}~Solution~1")
h_Q22.GetXaxis().SetTitle("q^{2}~Solution~2")
f_Histogram.cd()
SaveHistogram( f_Histogram, h_MCORR )
SaveHistogram( f_Histogram, h_MM2 )
SaveHistogram( f_Histogram, h_Q21 )
SaveHistogram( f_Histogram, h_Q22 )
SaveHistogram( f_Histogram, h_Q20 )
f_Histogram.Close()
In [5]:
f_Histogram = ROOT.TFile.Open("Source_Histograms_DsMu_{0}_LHCb.root".format(0.5), "RECREATE")
for EvType, FileTree in DataTuple.iteritems():
PickleName = "/home/ismith/tmp/Bs_Ds_" + EvType + ".pickle"
InputData = None
try:
InputData,K1_PE,K2_PE,Pi_PE,Mu_PE,B_PE,B_Origin,B_End = LoadData_DsMuNu(EvType, FileTree)
except:
continue
DataType = InputData[K_PE[0]].dtype
K1 = FourVector(InputData[K1_PE])
K2 = FourVector(InputData[K2_PE])
Pi = FourVector(InputData[Pi_PE])
Mu = FourVector(InputData[Mu_PE])
B = FourVector(InputData[B_PE])
Ds = K1 + K2 + Pi
Y = Ds + Mu
QSQ = ( B - Ds ).M2()
B_SV = InputData[B_Origin].copy().view(( DataType, 3 ))
B_PV = InputData[B_End].copy().view(( DataType, 3 ))
B_SV_New = SmearVertex( B_SV, sigma_SV_LHCb * 0.5)
B_PV_New = SmearVertex( B_PV, sigma_PV_LHCb * 0.5)
B_Direction = ThreeVector(B_SV_New - B_PV_New)
Cuts = ()
Cuts += ( K1.P() > 10000., )
Cuts += ( K1.Pt() > 500, )
Cuts += ( K2.P() > 10000, )
Cuts += ( K2.Pt() > 500, )
Cuts += ( Pi.P() > 2000, )
Cuts += ( Pi.Pt() > 250, )
Cuts += ( Mu.P() > 6000., )
Cuts += ( Mu.Pt() > 1500., )
Cuts += ( np.abs( Ds.M() - 1968 ) < 80, )
Cuts += ( Y.M() > 2200, )
Cuts += ( Y.M() < 8000, )
Cuts += ( -np.cos( B_Direction.Angle( Y.Vect() ) ) > 0.999, )
PassCuts = np.all( Cuts, axis = 0)
print EvType, 1.0*np.sum(PassCuts)/ K1.P().shape[0]
#print Cuts
#print PassCuts
h_MCORR = ROOT.TH1F("MCORR_" + EvType, "", 100, 2500, 6000)
h_MM2 = ROOT.TH1F("MissingMass2_" + EvType, "", 100, -10e6, 20e6)
h_Q21 = ROOT.TH1F("QSQ_SOL1_" + EvType, "", 100, 0, 12e6)
h_Q22 = ROOT.TH1F("QSQ_SOL2_" + EvType, "", 100, 0, 12e6)
h_Q20 = ROOT.TH1F("QSQ_SOLTrue_" + EvType, "", 100, 0, 12e6)
Q20 = (B - Ds).M2()
MCORR = GetCorrectedMass(Y[PassCuts], B_Direction[PassCuts])
MM2 = GetMissingMass2(Ds[PassCuts], Mu[PassCuts], B_Direction[PassCuts])
Q21, Q22 = GetQ2(Y[PassCuts], Mu[PassCuts], B_Direction[PassCuts])
nev = MM2.shape[0]
#for MC, MM in zip(MCORR, MM2):
h_MCORR.FillN( nev, MCORR, np.ones(nev))
h_MM2 .FillN( nev, MM2, np.ones(nev))
h_Q21 .FillN( nev, Q21,np.ones(nev))
h_Q22 .FillN( nev, Q22,np.ones(nev))
h_Q20 .FillN( nev, QSQ,np.ones(nev))
h_MCORR.GetXaxis().SetTitle("m_{corr}(B_{s})")
h_MM2.GetXaxis().SetTitle("Missing Mass^{2}")
h_Q21.GetXaxis().SetTitle("q^{2}~Solution~1")
h_Q22.GetXaxis().SetTitle("q^{2}~Solution~2")
f_Histogram.cd()
SaveHistogram( f_Histogram, h_MCORR )
SaveHistogram( f_Histogram, h_MM2 )
SaveHistogram( f_Histogram, h_Q21 )
SaveHistogram( f_Histogram, h_Q22 )
SaveHistogram( f_Histogram, h_Q20 )
f_Histogram.Close()
In [6]:
Signal_Yields = {}
Signal_Yields["13512010"] = 17000. # Signal
Signal_Yields["13774000_Dsstar"] = 40000.
Signal_Yields["13774000_Ds"] = 20000.
Signal_Yields["13774000_DsTau"] = 3000
Signal_Yields["13774000_DsstarTau"] = 5000
Signal_Yields["13774000_Dsstar1"] = 1000
Signal_Yields["13774000_Dsstar0"] = 1000
Signal_Yields["13114006"] = 1000
Signal_Yields["12513010"] = 5000
Signal_Yields["23573002"] = 1000
Signal_Yields["12145052"] = 5000
Signal_Yields["15574010"] = 2000
Signal_Yields["11876001"] = 1000
Control_Yields = {}
Control_Yields["13774000_Dsstar"] = 4.42e5
Control_Yields["13774000_Ds"] = 1.97e5
Control_Yields["13774000_DsTau"] = 5e3
Control_Yields["13774000_DsstarTau"] = 8e3
Control_Yields["13774000_Dsstar1"] = 2e4
Control_Yields["13774000_Dsstar0"] = 5e4
Control_Yields["13146004"] = 1000
Control_Yields["11876001"] = 1500 #Double D
In [7]:
DataFile = ROOT.TFile.Open("Data_Histograms_{0}_LHCb.root".format(0.5), "RECREATE")
SignalFile = ROOT.TFile.Open("Source_Histograms_KMu_{0}_LHCb.root".format(0.5), "READ")
ControlFile = ROOT.TFile.Open("Source_Histograms_DsMu_{0}_LHCb.root".format(0.5), "READ")
h_MCORR = ROOT.TH1F("MCORR_Data_KMuNu", "", 100, 2500, 6000)
for hname, Yield in Signal_Yields.iteritems():
print "MCORR_" + hname
temphist = SignalFile.Get("MCORR_" + hname)
#print temphist
for ev in xrange(int(Yield)):
MCORR = temphist.GetRandom()
h_MCORR.Fill(MCORR)
DataFile.cd()
h_MCORR.Write()
h_MCORR = ROOT.TH1F("MCORR_Data_DsMuNu", "", 100, 2500, 6000)
for hname, Yield in Control_Yields.iteritems():
print "MCORR_" + hname
temphist = ControlFile.Get("MCORR_" + hname)
for ev in xrange(int(Yield)):
MCORR = temphist.GetRandom()
h_MCORR.Fill(MCORR)
DataFile.cd()
h_MCORR.Write()
DataFile.Close()
In [ ]: