python_simulation/AngularDistribution.py at master · sgiani/python_simulation · 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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#Angular distribution module
"""
.. module:: GenerateEvent
   :synopsis: Definition of Angular Distribution functions for the photons
.. moduleauthor:: Sebastiana Giani <sebastiana.giani@epfl.ch>
"""

import math
from ROOT import TF1
import unittest


def function_1 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -0.7367+5.0550*x[0]-0.0554*x[0]*x[0]


def function_2 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -0.2846+13.3962*x[0]-0.0838*x[0]*x[0]


def function_3 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -1.9196+25.6074*x[0]-0.2702*x[0]*x[0]


def function_4 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -1.2325+34.1110*x[0]-0.3271*x[0]*x[0]

def function_5 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -6.7858+44.746*x[0]-0.4284*x[0]*x[0]


def function_6 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<= 27:
        return -8.5570+56.2115*x[0]-0.5866*x[0]*x[0]
    else:
        return 7451.11-335.88*x[0]+3.6112*x[0]*x[0]


def function_7 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -6.0329+62.0612*x[0]-0.4463*x[0]*x[0]
    else:
        return 6296.19-257.260*x[0]+2.5809*x[0]*x[0]


def function_8 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=27:
        return -2.3125+71.1939*x[0]-0.4842*x[0]*x[0]
    else:
        return 4405.57-141.793*x[0]+1.1392*x[0]*x[0]


def function_9 (x):

    if x[0]<=0.25:
        return 0.0
    if x[0]<=27:
        return -18.2144+82.3692*x[0]-0.9137*x[0]*x[0]
    else:
        return 2221.03-40.1972*x[0]+0.1705*x[0]*x[0]


def function_10 (x):

    if x[0]<=0.2:
        return 0.0
    if x[0]<=16:
        return -0.4247+18.1516*x[0]+0.4020*x[0]*x[0]
    else:
        return 1473.86*math.exp(-x[0]/11.7009)


def choose_angle_from_distribution(radius):

    """

     Given the position of the photon rispect to the center of the fiber(radius), this function return us the emission angle
     of the photons.Depending by the radius, a random number is returned from the distributions obtained by a Geant Simulation"""
    #print radius

    angle=0
    if radius>=0 and radius<=12.5:
        f1=TF1("f1",function_1,0,27)
        angle=f1.GetRandom()
    elif radius>12.5 and radius<=25:
        f2=TF1("f2",function_2,0,27)
        angle=f2.GetRandom()
    elif radius>25 and radius<=37.5:
        f1=TF1("f1",function_3,0,27)
        angle=f1.GetRandom()
    elif radius>37.5 and radius<=50:
        f1=TF1("f1",function_4,0,27)
        angle=f1.GetRandom()
    elif radius>50 and radius<=62.5:
        f1=TF1("f1",function_5,0,27)
        angle=f1.GetRandom()
    elif radius>62.5 and radius<=75:
        f1=TF1("f1",function_6,0,37)
        angle=f1.GetRandom()
    elif radius>75 and radius<=87.5:
        f1=TF1("f1",function_7,0,44)
        angle=f1.GetRandom()
    elif radius>87.5 and radius<=100:
        f1=TF1("f1",function_8,0,60)
        angle=f1.GetRandom()
    elif radius>100 and radius<=112.5:
        f1=TF1("f1",function_9,0,86)
        angle=f1.GetRandom()
    elif radius>=112.5 and radius<=125:
        f1=TF1("f1",function_10,0,70)
        angle=f1.GetRandom()

    #print angle
    return math.radians(angle)