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)


Welcome to JupyROOT 6.08/02

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()


14966 14966
13774011 0.196445275959
1418250 1418250
13774000_Dsstar 0.219750396616
617460 617460
13774000_Ds 0.209266997052
577349 577349
12513010 0.150834244105
617458 617458
13774000_combi 0.0874197111382
67603 67603
13774000_Dsstar1 0.174474505569
18298 18298
13146004 0.221608919008
21721 21721
12145052 0.140417107868
16681 16681
15574010 0.143336730412
2288 2288
11876001 0.0380244755245
1537542 1537542
13512010 0.391815638207
1537540 1537540
13512010_combi 0.166602494894
41776 41776
13774000_DsTau 0.0942885867484
19570 19570
13114006 0.259427695452
27181 27181
23573002 0.0202347227843
79367 79367
13774000_DsstarTau 0.0972696460746
205840 205840
13774000_Dsstar0 0.179309172173
21005 21005
13114001 0.257224470364
FutureColliderTools.py:158: RuntimeWarning: invalid value encountered in sqrt
  px_rest = np.sqrt( E_nurest * E_nurest  - Ppmu.PY() * Ppmu.PY() );
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13144032.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV22114001.root does not exist

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()


13774011 0.107643993051
13774000_Dsstar 0.114962101181
13774000_Ds 0.116830240016
13774000_combi 0.116746024034
13774000_Dsstar1 0.0905729035694
13146004 0.00120231719314
11876001 0.0262237762238
13774000_DsTau 0.0539544235925
13774000_DsstarTau 0.0526919248554
13774000_Dsstar0 0.0955499417023
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13144032.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV12513010.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV12145052.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV15574010.root does not exist
Error in <TUnixSystem::GetHostByName>: getaddrinfo failed for 'eoslhcb.cern.ch': 'Temporary failure in name resolution'
Error in <TNetXNGFile::Open>: [FATAL] Invalid address
Error in <TUnixSystem::GetHostByName>: getaddrinfo failed for 'eoslhcb.cern.ch': 'Temporary failure in name resolution'
Error in <TNetXNGFile::Open>: [FATAL] Invalid address
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13114006.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV23573002.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV22114001.root does not exist
Error in <TFile::TFile>: file /afs/cern.ch/work/i/ismith/NTuples/VubFromBplus/Backgrounds/GenLevel/DV13114001.root does not exist

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()


MCORR_13774000_Dsstar
MCORR_15574010
MCORR_12513010
MCORR_11876001
MCORR_12145052
MCORR_13774000_Ds
MCORR_23573002
MCORR_13512010
MCORR_13774000_DsTau
MCORR_13114006
MCORR_13774000_DsstarTau
MCORR_13774000_Dsstar0
MCORR_13774000_Dsstar1
MCORR_13774000_Dsstar
MCORR_13146004
MCORR_13774000_Ds
MCORR_11876001
MCORR_13774000_DsTau
MCORR_13774000_DsstarTau
MCORR_13774000_Dsstar0
MCORR_13774000_Dsstar1

In [ ]: