SSI2016_project/plot.py at master · lime86/SSI2016_project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
import sys
import os
import glob
import array
import argparse
import collections

from ROOT import *
#import TMath
from math import *


canvas = TCanvas("name","I am a canvas",1000,500)
canvas.Divide(2,1)
mg = []
legend = []

#energy GeV
E = [5, 11, 200]

#mm
d_dut = [20, 150]

dxbyx0_tel = 0.001
dxbyx0_dut = 0.1


#GeV
def theta(E,dx):
	_theta = (13.6e-3 / E) * sqrt(dx) * (1+ 0.038 * log(dx))
	return _theta


def myfunction(x, param):
	_E=param[0]
	_dut=param[1]
	f = 1./(3.5e-3)
	#g = 1./(5e-3)
	g = 0
	u = 1./(x*theta(_E,dxbyx0_tel))
	v = 1./(_dut*theta(_E,dxbyx0_dut))

#	function = sqrt( (2*v**2 + 3*f**2)/(12*f**2*v**2))
#	function = sqrt( \
#			( f**6 + 2.*(u**4) * (v**2) + (f**4) * ( 8.*(u**2) + 2.*u*v + 3.*(v**2) ) + (f**2) * (u**2) * ( 3.*(u**2) + 6.*u*v + 17.*(v**2) ) )/ \
#			((2.*(g**2) * (u**4) * (v**2) + (f**6) * ( g**2 + 6.*(v**2) ) + (f**4) * ( 44.*(u**2) * (v**2) + (g**2) * (8.*(u**2) + 2.*u*v +3.*(v**2) ))) \
#			+ (f**2) * (u**2) * ( 12.*(u**2) * (v**2) + (g**2) * ( 3.*(u**2) + 6.*u*v + 17.*(v**2))))

	num = ( f**6 + 2.*(u**4) * (v**2) + (f**4) * ( 8.*(u**2) + 2.*u*v + 3.*(v**2) ) + (f**2) * (u**2) * ( 3.*(u**2) + 6.*u*v + 17.*(v**2) ) )
	denom = ((f**6) * ( g**2 + 6.*(v**2) ) + (f**4) * ( 44.*(u**2) * (v**2) + (g**2) * (8.*(u**2) + 2.*u*v +3.*(v**2) ))) + (f**2) * (u**2) * ( 12.*(u**2) * (v**2) + (g**2) * ( 3.*(u**2) + 6.*u*v + 17.*(v**2)))
	function = sqrt(num/denom)
	#print num, denom
	return function

n=0
flist = []
thetalist = []

for j in d_dut:
	_mg = TMultiGraph()
	_legend = TLegend(0.74,0.78,0.96,0.96)
	color = 1
	for i in E:
		thetalist.append(theta(i,dxbyx0_tel))
		n += 1
		param = []
		param.append(i)
		param.append(j)
		#print param
		graph = TGraph(200)
		for k in range(5,500):
			graph.SetPoint(k, k, myfunction(k,param)*1000)

		graph.SetMarkerColor(color)
		graph.SetMarkerStyle(7)
		#graph.SetMarkerColor(color)
		label = str(i) + " GeV"
		_legend.AddEntry(graph,label,"p")
		_mg.Add(graph)
		color += 1


	mg.append(_mg)
	legend.append(_legend)
	n = d_dut.index(j)

	canvas.cd(n+1)
	mg[n].Draw("AP")
	legend[n].Draw("SAME")

	graphtitle = "d_{dut} = " + str(j) + "mm"
	xaxis = "d (mm)"
	yaxis = "#sigma_{dut} (#mum)"
	mg[n].SetTitle(graphtitle)
	mg[n].GetXaxis().SetTitle(xaxis)
	mg[n].GetYaxis().SetTitle(yaxis)
	mg[n].SetMinimum(1)
	mg[n].SetMaximum(gPad.GetUymax()*1.1)

	canvas.cd(n+1).Update()

	#canvas.cd(n+1).SetLogy(1)


raw_input()