multifrag-test/reallyBasicTest.py at master · ffavela/multifrag-test · 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
146
from multifragStuff import *

print(colored(easyStr,"blue"))
#print all energies in MeV
beamE=67.2
m6He=5606.55669604
m4He=3728.401315862896
m3He=2809.41351708

m3H=2809.43210768
m2H=1876.12392312
m1H=938.783071355
m1n=939.565417991

#Types of dictionaries
tP="particle"
tS="set" #particle set
tD="detector"
tI="initial"


##################################################
##########Particle dict ternary 6He case tritons##########
##################################################
initDict={"type":tI,"name":"a+t","massP":m4He,
          "massT":m3H,"ELab":beamE}
helium6Dict={"type":tP,"name":"6He","mass":m6He,"exE":14.2}
# helium6Dict={"type":tP,"name":"6He","mass":m6He,"exE":15.5}
#The actual state seems to be 15.8 for iii, with a contribution in
#18.67 and 9.4 in the spectra fig4(b).

#Same fig the peaks at 13.73,28.11,9.5 and 30.8 seem to come from the
#18.6 *6He excitation level.

#Notice the overlap of peaks below the 10MeV part.

#Was not able to reproduce any *4He decay fragmentation.

protonDict={"type":tP,"name":"p","mass":m1H}
tritium1Dict={"type":tP,"name":"t1","mass":m3H,"exE":0.0}
tritium2Dict={"type":tP,"name":"t2","mass":m3H,"exE":0.0}

#Defining the detectors
# d1Dict={"type":tD,"name":"d1","angles":[radians(20),radians(180)]}
d2Dict={"type":tD,"name":"d2","angles":[radians(10),radians(0)]}

#Completing the dictionaries
protonDict["dictList"]=[{},{}]
# tritium1Dict["dictList"]=[d1Dict,{}]
# tritium2Dict["dictList"]=[d2Dict,{}]

# helium6Dict["dictList"]=[tritium1Dict,tritium2Dict]
helium6Dict["dictList"]=[d2Dict,{}]

initDict["dictList"]=[protonDict,helium6Dict]

# ######################################################
# ##########Particle end dict ternary 6He case tritons######
# ##################################################

makeTreeCompletion(initDict)
# makeInitialTreeCompletion(initDict)
printTree(initDict)

print("")
print(colored("The energy print out function","yellow"))
printLastNodes(initDict)

genSimpVCMD=getGeneralSimplifiedVCMD(initDict)
print(genSimpVCMD)
print(colored(easyStr,"red"))

###########################################
##Some extra tests#########################
###########################################

# mIn=m4He+m3H

# mOut=m6He+m1H


# print("mIn,mOut = %0.4f,%0.4f" % (mIn,mOut))
# print("beamE = %0.4f" % beamE)
# ECMListIn=getAllEcms(m4He,m3H,beamE)
# print(ECMListIn)
# ECMSysIn=ECMListIn[3]
# ECMAvailIn=ECMListIn[2]
# print("ECMSysIn,ECMAvailIn = %0.4f,%0.4f" %(ECMSysIn,ECMAvailIn))

# ETotIn=ECMSysIn+ECMAvailIn
# print("ETotIn = %f" %(ETotIn))

# ECMSysOut=1.0*(mIn/mOut)*ECMSysIn
# Q=getQVal(mIn,0,mOut,0)
# print("Q=%0.4f" %(Q))
# Ex=18.6
# # Ex=0.0
# Q+=Ex
# ECMAvailOut=ECMSysIn*(1.0-1.0*mIn/mOut)+ECMAvailIn+Q
# print("ECMSysOut,ECMAvailOut = %0.4f,%0.4f" %(ECMSysOut,ECMAvailOut))
# ETotOut=ECMSysOut+ECMAvailOut-Q
# print("ETotOut = %f" %(ETotOut))

# print("Getting the respective values for the ECM of the 6He and proton")
# ECM6He,ECM1H=getEcmsFromECM2(m6He,m3H,ECMAvailOut)
# print("ECM6He,ECM1H = %0.4f,%0.4f" %(ECM6He,ECM1H))

# eTempSum=ECM6He+ECM1H
# print("eTempSum = %0.4f" %(eTempSum))


############################################
####The plotting part#######################
############################################

fig = plt.figure()
ax = fig.gca(projection='3d')


fig_size = plt.rcParams["figure.figsize"]
# print("The figsize is = ",fig_size)
# Set figure width to 9 and height to 9, a square
fig_size[0] = 8
fig_size[1] = 4
plt.rcParams["figure.figsize"] = fig_size
# print("The new figsize is = ",fig_size)

plt.xlim(-15, 15)
plt.ylim(-15, 15)

plotAllLines(initDict,ax)
ax.legend()

modifyAx4Arrows(ax,genSimpVCMD)

u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

ax.set_zlim3d(-15, 15)

plt.show()

############################################
####The plotting part finish################
############################################