#Function that obtains any value in a graph def getVal(graph, val): leftPoint = -1 for i in range(1,graph.GetN()): #Here I check that if I give the exact value that it is already in the graph, it gives it and doesn't do the interpolation print(graph.GetPointX(i)) print(val) if int(graph.GetPointX(i))==int(val): print("llego aqui tb");return graph.GetPointY(i) elif graph.GetPointX(i) > val: break leftPoint = i # So now leftPoint is the point of the graph just before the value you want to interpolate print("llego aqui") leftY = graph.GetPointY(leftPoint) leftX = graph.GetPointX(leftPoint) rightY = graph.GetPointY(leftPoint+1) rightX= graph.GetPointX(leftPoint+1) #Careful: You have to treat the edge cases where val is smaller than the minimum and, worse, larger than the maximum #I would sinply ignore the first slope = (rightY-leftY)/(rightX-leftX) return leftY + slope * (val-leftX) #The previous function is no longer used #This function shifts the two arrays so I can put them in the TGraph in the good order to do later the interpolation def shift(arr1,arr2,arr3): aux1= array( 'f', [] ) aux2= array( 'f', [] ) aux3= array( 'f', [] ) for i in range(len(arr1)): a = min(arr1) index = arr1.index(a) aux1.append(a) aux2.append(arr2[index]) aux3.append(arr3[index]) arr1.pop(index) arr2.pop(index) arr3.pop(index) return aux1,aux2,aux3 #This function works out the value of track proportion from known values (well, actually it calls the getVal function and then obtain the expected values) def getValInterp(trackType,lifetime_val,mass_val): gr_lifetime = createGraph("lifetime","mass",trackType) gr_lifetime.Draw("ACP") gr_mass = createGraph("mass","lifetime",trackType) gr_mass.Draw("ACP") X = getValSpline(gr_lifetime,lifetime_val) Y = getValSpline(gr_mass,mass_val) AUX = getValSpline(gr_mass,2500) del gr_lifetime,gr_mass return Y*X/AUX #https://root.cern.ch/doc/master/classTSpline3.html #To obtain the spline function from a graph which has been drawn with the C option from ROOT import * def getValSpline(graph,val): spline = TSpline3("bla",graph) if spline.Eval(val)>1: returnVal = 1 elif spline.Eval(val)<0: returnVal = 0 else: returnVal = spline.Eval(val) #returnVal = 1 if spline.Eval(val) > 1 else spline.Eval(val) del spline return returnVal #Function that creates a graph for the quantity in common: lifetime or mass from Particle_data import particles_info, lhcb_info, labels from array import array from math import * import numpy as np def createGraph(prop,propEq,trackType,darkbosonType): #if propEq and prop in ["lifetime","mass"] or trackType in ["LL","DD","TT","rest"] == False: return "Invalid arguments"; return x_val = array( 'f', [] ) y_val = array( 'f', [] ) y_err = array( 'f', [] ) x_err = array( 'f', [] ) n = 0 darkboson = particles_info().scalarDark_boson if darkbosonType == "scalar" else particles_info().vectorDark_boson for key in list(darkboson.keys()): if int(darkboson[key][propEq]) in [100,2500]: x_val.append(float(darkboson[key][prop])) num = float(darkboson[key][trackType]) den = float(darkboson[key]["LL"])+float(darkboson[key]["DD"])+float(darkboson[key]["TT"])+float(darkboson[key]["rest"]) y_val.append(num/den) y_err.append(abs(num/den - max(TEfficiency.ClopperPearson(den,num,0.68,False),TEfficiency.ClopperPearson(den,num,0.68,True)))) n += 1 x_err.append(0.) x_val_shifted, y_val_shifted, y_err_shifted = shift(x_val,y_val,y_err) grph = TGraphErrors(n,x_val_shifted,y_val_shifted,x_err,y_err_shifted) return grph def create2Dgraph(lifetimes,masses,trackType,darkbosonType): darkboson = particles_info().scalarDark_boson_reconstr_hack #if darkbosonType == "scalar" else particles_info().vectorDark_boson numPoints = len(lifetimes)*len(masses) grph2d = TGraph2D(numPoints) n = 0 for i in range(len(masses)): for j in range(len(lifetimes)): decnum = masses[i]+j grph2d.SetPoint(n,masses[i],int(lifetimes[j]*1.e12),float(darkboson[str(decnum)][trackType])/(float(darkboson[str(decnum)]["LL"])+float(darkboson[str(decnum)]["DD"])+float(darkboson[str(decnum)]["TT"])+float(darkboson[str(decnum)]["rest"]))) n+=1 return grph2d if True: lifetime = [1.0000e-12,1.0000e-11,1.0000e-10,2.5000e-10,5.0000e-10,7.5000e-10,1.0000e-9,1.2500e-9,1.5000e-9,1.7500e-9,2.000e-9] mass = [500,1000,1500,2000,2500,3000,3500,4000,4500] C = TCanvas() C.SetRightMargin(0.20) C.SetBottomMargin(0.13) C.SetLeftMargin(0.15) gr2d = create2Dgraph(lifetime,mass,"TT","scalar") gr2d.SetTitle(';M (H_{0}) [MeV];#tau (H_{0}) [ps];"TT" vertex proportion ') gr2d.SetNpx(30) gr2d.SetNpy(30) gr2d.Draw("colz") h2dt = gr2d.GetHistogram() h2dt.GetYaxis().CenterTitle() h2dt.GetXaxis().CenterTitle() h2dt.GetZaxis().CenterTitle() h2dt.GetYaxis().SetTitleSize(0.045) h2dt.GetXaxis().SetTitleSize(0.045) h2dt.GetZaxis().SetTitleSize(0.045) h2dt.GetYaxis().SetLabelSize(0.04) h2dt.GetXaxis().SetLabelSize(0.04) h2dt.GetZaxis().SetLabelSize(0.04) h2dt.GetXaxis().SetLabelOffset(0.005) h2dt.GetYaxis().SetLabelOffset(0.005) h2dt.GetZaxis().SetLabelOffset(0.005) h2dt.GetXaxis().SetTitleOffset(1.3) h2dt.GetYaxis().SetTitleOffset(1.3) h2dt.GetZaxis().SetTitleOffset(1.3) h2dt.GetXaxis().SetLabelFont(132) h2dt.GetYaxis().SetLabelFont(132) h2dt.GetZaxis().SetLabelFont(132) h2dt.GetXaxis().SetTitleFont(132) h2dt.GetYaxis().SetTitleFont(132) h2dt.GetZaxis().SetTitleFont(132) h2dt.GetZaxis().SetRangeUser(0,0.7) l = TLegend(0.5,0.67,0.85,0.89) l.SetHeader("LHCb simulation","L") # l.AddEntry("l1","LHCb simulation","") # l.AddEntry(grphs1d[0],"H_{0} lifetime","l") l.SetFillStyle(0) l.SetBorderSize(0) l.SetTextFont(132) l.SetTextSize(0.05) l.Draw() # C.SetLogy() C.SaveAs("TT2D_hack_without_theta.png") C.SaveAs("TT2D_hack_without_theta.pdf") del C C = TCanvas() C.SetRightMargin(0.20) C.SetBottomMargin(0.13) C.SetLeftMargin(0.15) gr2dd = create2Dgraph(lifetime,mass,"DD","scalar") gr2dd.SetTitle(';M (H_{0}) [MeV];#tau (H_{0}) [ps];"DD" vertex proportion') gr2dd.SetNpx(30) gr2dd.SetNpy(30) gr2dd.Draw("colz") h2dd = gr2dd.GetHistogram() h2dd.GetYaxis().CenterTitle() h2dd.GetXaxis().CenterTitle() h2dd.GetZaxis().CenterTitle() h2dd.GetYaxis().SetTitleSize(0.045) h2dd.GetXaxis().SetTitleSize(0.045) h2dd.GetZaxis().SetTitleSize(0.045) h2dd.GetYaxis().SetLabelSize(0.04) h2dd.GetXaxis().SetLabelSize(0.04) h2dd.GetZaxis().SetLabelSize(0.04) h2dd.GetXaxis().SetLabelOffset(0.005) h2dd.GetYaxis().SetLabelOffset(0.005) h2dd.GetZaxis().SetLabelOffset(0.005) h2dd.GetXaxis().SetTitleOffset(1.3) h2dd.GetYaxis().SetTitleOffset(1.3) h2dd.GetZaxis().SetTitleOffset(1.3) h2dd.GetXaxis().SetLabelFont(132) h2dd.GetYaxis().SetLabelFont(132) h2dd.GetZaxis().SetLabelFont(132) h2dd.GetXaxis().SetTitleFont(132) h2dd.GetYaxis().SetTitleFont(132) h2dd.GetZaxis().SetTitleFont(132) h2dd.GetZaxis().SetRangeUser(0,0.7) C.SaveAs("DD2D_colz_rec.png") C.SaveAs("DD2D_colz_rec.pdf") del C C = TCanvas() C.SetRightMargin(0.20) C.SetBottomMargin(0.13) C.SetLeftMargin(0.15) gr2dl = create2Dgraph(lifetime,mass,"LL","scalar") gr2dl.SetTitle(';M (H_{0}) [MeV];#tau (H_{0}) [ps];"LL" vertex proportion') gr2dl.SetNpx(30) gr2dl.SetNpy(30) gr2dl.Draw("colz") h2dl = gr2dl.GetHistogram() h2dl.GetYaxis().CenterTitle() h2dl.GetXaxis().CenterTitle() h2dl.GetZaxis().CenterTitle() h2dl.GetYaxis().SetTitleSize(0.045) h2dl.GetXaxis().SetTitleSize(0.045) h2dl.GetZaxis().SetTitleSize(0.045) h2dl.GetYaxis().SetLabelSize(0.04) h2dl.GetXaxis().SetLabelSize(0.04) h2dl.GetZaxis().SetLabelSize(0.04) h2dl.GetXaxis().SetLabelOffset(0.005) h2dl.GetYaxis().SetLabelOffset(0.005) h2dl.GetZaxis().SetLabelOffset(0.005) h2dl.GetXaxis().SetTitleOffset(1.3) h2dl.GetYaxis().SetTitleOffset(1.3) h2dl.GetZaxis().SetTitleOffset(1.3) h2dl.GetXaxis().SetLabelFont(132) h2dl.GetYaxis().SetLabelFont(132) h2dl.GetZaxis().SetLabelFont(132) h2dl.GetXaxis().SetTitleFont(132) h2dl.GetYaxis().SetTitleFont(132) h2dl.GetZaxis().SetTitleFont(132) h2dl.GetZaxis().SetRangeUser(0,0.7) C.SaveAs("LL2D_colz_rec.png") C.SaveAs("LL2D_colz_rec.pdf") del C C = TCanvas() ht = h2dt.ProjectionY("",14,15) ht.SetMarkerColor(kRed) ht.SetLineStyle(1) ht.SetLineColor(kRed) ht.SetLineWidth(2) ht.GetYaxis().SetRangeUser(0,0.75) ht.SetTitle(";#tau (H_{0}) [ps];Vertex proportion") ht.SetStats(0) # # hs.Add(ht) hd = h2dd.ProjectionY("",14,15) hd.SetMarkerColor(kBlue) hd.SetLineStyle(1) hd.SetLineColor(kBlue) hd.SetLineWidth(2) hd.GetYaxis().SetRangeUser(0,0.75) # hs.Add(hd) hl = h2dl.ProjectionY("",14,15) hl.SetMarkerColor(kBlack) hl.SetLineStyle(1) hl.SetLineColor(kBlack) hl.SetLineWidth(2) hl.GetYaxis().SetRangeUser(0,0.75) ht.Draw("HIST") hd.Draw("SAME HIST") hl.Draw("SAME HIST") # hs.Add(hl) # hs.Draw() # l = TLegend(0.35,0.57,0.75,0.89) # l.SetHeader("LHCb simulation","L") # l.AddEntry(ht,"M(H_{0})=250 [MeV/c^{2}], vertex type: TT","LP") # l.AddEntry(hd,"M(H_{0})=250 [MeV/c^{2}], vertex type: DD","LP") # l.AddEntry(hl,"M(H_{0})=250 [MeV/c^{2}], vertex type: LL","LP") # # l.AddEntry(mg.GetListOfGraphs()[0],"H_{0} lifetime","l") # l.SetFillStyle(0) # l.SetBorderSize(0) # l.SetTextFont(132) # l.Draw() C.Modified() C.SaveAs("trial.png") # C = TCanvas() # gr2d = create2Dgraph(lifetime,mass,"LL","scalar") # gr2d.SetTitle(';M (H_{0}) [MeV];#tau (H_{0}) [ps];"LL" vertex proportion [%]') # gr2d.SetNpx(500) # gr2d.SetNpy(500) # gr2d.Draw("colz") # C.SetRightMargin(0.2) # C.SaveAs("LL2D_colz_rec.png") # C.SaveAs("LL2D_colz_rec.pdf") # del C,gr2d # tau = 50 # for i in range(20): # tau2 = tau + i*50 # m = 500 # for j in range(9): # m2 = m + j*500 # print("mass: ",m2," lifetime: ",tau2," LL prop: ",getValInterp("LL",tau2,m2)," DD prop: ",getValInterp("DD",tau2,m2)," TT prop: ",getValInterp("TT",tau2,m2)) # C = TCanvas() # C.SetLogx() # C.SetGrid() # mg = TMultiGraph(); # mg.SetTitle(';#tau (H) [ps] ;Vertex proportion [%]') # gr_scalar_LL = createGraph("lifetime","mass","LL","scalar") # gr_scalar_LL.SetTitle("LL, scalar dark boson") # gr_scalar_LL.SetMarkerStyle(kFullSquare) # gr_scalar_LL.SetMarkerSize(0.6) # gr_scalar_LL.SetMarkerColor(kRed) # gr_scalar_LL.SetLineColor(kRed) # gr_scalar_DD = createGraph("lifetime","mass","DD","scalar") # gr_scalar_DD.SetTitle("DD, scalar dark boson") # gr_scalar_DD.SetMarkerStyle(kFullSquare) # gr_scalar_DD.SetMarkerSize(0.6) # gr_scalar_DD.SetMarkerColor(kGreen) # gr_scalar_DD.SetLineColor(kGreen) # gr_scalar_TT = createGraph("lifetime","mass","TT","scalar") # gr_scalar_TT.SetTitle("TT, scalar dark boson") # gr_scalar_TT.SetMarkerStyle(kFullSquare) # gr_scalar_TT.SetMarkerSize(0.6) # gr_scalar_TT.SetMarkerColor(6) # gr_scalar_TT.SetLineColor(6) # gr_vector_LL = createGraph("lifetime","mass","LL","vector") # gr_vector_LL.SetTitle("LL, vector dark boson") # gr_vector_LL.SetMarkerStyle(kFullCircle) # gr_vector_LL.SetMarkerSize(0.6) # gr_vector_LL.SetMarkerColor(kBlue) # gr_vector_LL.SetLineColor(kBlue) # gr_vector_DD = createGraph("lifetime","mass","DD","vector") # gr_vector_DD.SetTitle("DD, vector dark boson") # gr_vector_DD.SetMarkerStyle(kFullCircle) # gr_vector_DD.SetMarkerSize(0.6) # gr_vector_DD.SetMarkerColor(kBlack) # gr_vector_DD.SetLineColor(kBlack) # gr_vector_TT = createGraph("lifetime","mass","TT","vector") # gr_vector_TT.SetTitle("TT, vector dark boson") # gr_vector_TT.SetMarkerStyle(kFullCircle) # gr_vector_TT.SetMarkerSize(0.6) # gr_vector_TT.SetMarkerColor(28) # gr_vector_TT.SetLineColor(28) # mg.Add(gr_scalar_LL) # mg.Add(gr_scalar_DD) # mg.Add(gr_scalar_TT) # mg.Add(gr_vector_LL) # mg.Add(gr_vector_DD) # mg.Add(gr_vector_TT) # mg.Draw("ALP") # # mg.SetTitle(';mg.Add(gr_scalar_LL)#tau (H) [ps] ;"LL" vertex proportion [%]') # # gr_scalar = createGraph("lifetime","mass","LL","scalar") # # gr_scalar.SetTitle("Scalar dark boson") # # gr_scalar.SetMarkerStyle(kFullSquare) # # gr_scalar.SetMarkerSize(0.6) # # gr_scalar.SetMarkerColor(kRed) # # gr_scalar.SetLineColor(kRed) # # gr_vector = createGraph("lifetime","mass","LL","vector") # # gr_vector.SetTitle("Vector dark boson") # # gr_vector.SetMarkerStyle(kFullCircle) # # gr_vector.SetMarkerSize(0.6) # # gr_vector.SetMarkerColor(kBlue) # # gr_vector.SetLineColor(kBlue) # # mg.Add(gr_scalar) # # mg.Add(gr_vector) # # l = TLegend(0.12,0.12,0.4,0.3,"","") # # l.SetBorderSize(0) # # l.AddEntry("l1","LHCb simulation","") # # l.AddEntry("l2","pp #sqrt(14) = 14 TeV","") # # # l.AddEntry(gr_scalar,"","PL") # # # l.AddEntry(gr_vector,"","PL") # # l.Draw() # C.Update() # l = C.BuildLegend(0.12,0.4,0.55,0.70,"LHCb simulation pp #sqrt{s} = 14 TeV M = 2500 MeV","PL") # l.SetBorderSize(0) # # l.AddEntry("l1","LHCb simulation","") # C.SaveAs("./lifetime_all.png") # C.SaveAs("./lifetime_all.pdf") # del C, gr_scalar_LL,gr_scalar_DD,gr_scalar_TT, gr_vector_LL,gr_vector_DD,gr_vector_TT