détecteur à neutron scintillateur - monte carlo
import sys
import math
from random import *
from math import *
from ROOT import *
global nombre_muon
global nombre_plan
global longueur_z
global longueur_y
global longueur_x
global centre_plan
global normal_plan
global vitesse_lumiere
global vitesse_lumiere_vide
global position_source
global vitesse_photon
global speculaire
global run
global Lx
global Ly
global Lz
global sin_crit
global proba_muon
global i_scan_max
global delta_sample
global multi_hit
global indice_plan_bidon
global vitesse_milieu_bidon
global sin_crit_bidon
centre_plan={}
normal_plan={}
position_source=[]
vitesse_photon=[]
file_to_read='nebula.input'
print (file_to_read)
read_me=open(file_to_read,'r')
lines=read_me.readlines()
for line in lines:
cc=line.split('=')
if 'nombre_muon' in line:
nombre_muon=int(cc[1])
print ('nombre_muon',nombre_muon)
if 'nombre_plan' in line:
nombre_plan=int(cc[1])
print ('nombre_plan',nombre_plan)
if 'longueur_z' in line:
longueur_z=float(cc[1])
print ('longueur_z',longueur_z)
if 'longueur_x' in line:
longueur_x=float(cc[1])
print ('longueur_x',longueur_x)
if 'longueur_y' in line:
longueur_y=float(cc[1])
print ('longueur_y',longueur_y)
if 'indice_barre' in line:
indice_barre=float(cc[1])
if 'indice_pyra' in line:
indice_pyra=float(cc[1])
if 'speculaire' in line:
speculaire=float(cc[1])
if 'run' in line:
run=int(cc[1])
if 'hauteur_pyramide' in line:
hauteur_pyramide=float(cc[1])
if 'nombre_photon_par_cm' in line:
nombre_photon_par_cm=float(cc[1])
if 'nombre_photon_par_muon' in line:
nombre_photon_par_muon=float(cc[1])
if 'trigger' in line:
trigger=int(cc[1])
multi_hit=0
Lx=longueur_x/2.
Ly=longueur_y/2.
Lz=longueur_z/2.
##### longueur en cm, temps en ns, vitesse cm/ns
sin_crit=1./indice_barre
print (Lx,Ly,Lz)
print (sin_crit)
vitesse_lumiere_vide=30.
vitesse_lumiere=vitesse_lumiere_vide/indice_barre
#print 'celerite lumiere dans le milieu',vitesse_lumiere
###########################
def scalaire(k,vv):
prod_0=vv[0]*normal_plan[k][0]+vv[1]*normal_plan[k][1]+vv[2]*normal_plan[k][2]
return prod_0
########################################
def vitesse_gamma():
ctet=-1+random()*2
stet=sqrt(1-ctet*ctet)
phi=2*math.pi*random()
vz=vitesse_lumiere*ctet
vy=vitesse_lumiere*stet*cos(phi)
vx=vitesse_lumiere*stet*sin(phi)
return [vx,vy,vz]
########################################
def prod(ii,jj):
ax=ii[1]*jj[2]-ii[2]*jj[1]
ay=ii[2]*jj[0]-ii[0]*jj[2]
az=ii[0]*jj[1]-ii[1]*jj[0]
return[ax,ay,az]
######## ecrire les def ###############
def geometrie():
angle_pyra=68.
a=angle_pyra*3.1415/180.
G=4.0656
cos_a=cos(a)*G
sin_a=sin(a)*G
#print acos(.375)*180./3.1415
delta_z=hauteur_pyramide
delta_e=2.3 #2.3
############################
### points caracterisiriques:
A=[Lx,Ly,Lz]
B=[-Lx,Ly,Lz]
C=[Lx,-Ly,Lz]
D=[-Lx,-Ly,Lz]
E=[delta_e,delta_e,Lz+delta_z]
F=[delta_e,-delta_e,Lz+delta_z]
G=[-delta_e,delta_e,Lz+delta_z]
H=[-delta_e,-delta_e,Lz+delta_z]
###
Ap=[Lx,Ly,-Lz]
Bp=[-Lx,Ly,-Lz]
Cp=[Lx,-Ly,-Lz]
Dp=[-Lx,-Ly,-Lz]
Ep=[delta_e,delta_e,-Lz-delta_z]
Fp=[delta_e,-delta_e,-Lz-delta_z]
Gp=[-delta_e,delta_e,-Lz-delta_z]
Hp=[-delta_e,-delta_e,-Lz-delta_z]
###########################
######## plans de detection 0 et 1 ########################
centre_plan[0]=E
centre_plan[1]=Ep
##########################################################
### plans centraux ################################
centre_plan[2]=A ##### en haut x positif
centre_plan[3]=D ####### en bas y negatif
centre_plan[4]=B ########## y positif
centre_plan[5]=C ########## y negatif
#####################################################
#######################
centre_plan[6]=F
centre_plan[7]=G
centre_plan[8]=E
centre_plan[9]=H
####
centre_plan[10]=Fp
centre_plan[11]=Gp
centre_plan[12]=Ep
centre_plan[13]=Hp
############################
###########################
# vecteur normal pour chaque plan oriente vers linterieur
normal_plan[0]=[0,0,-1]
normal_plan[1]=[0,0,1]
normal_plan[2]=[-1,0,0]
normal_plan[3]=[1,0,0]
normal_plan[4]=[0,-1,0]
normal_plan[5]=[0,1,0]
###########################################
normal_plan[6]=prod_vec(F,E,F,C)
normal_plan[7]=prod_vec(G,H,G,B)
normal_plan[8]=prod_vec(G,B,G,E)
normal_plan[9]=prod_vec(F,C,F,H)
#############################################
normal_plan[10]=[normal_plan[6][0],normal_plan[6][1],-normal_plan[6][2]]
normal_plan[11]=[normal_plan[7][0],normal_plan[7][1],-normal_plan[7][2]]
normal_plan[12]=[normal_plan[8][0],normal_plan[8][1],-normal_plan[8][2]]
normal_plan[13]=[normal_plan[9][0],normal_plan[9][1],-normal_plan[9][2]]
###########################
indice_plan_bidon={}
vitesse_milieu_bidon={}
sin_crit_bidon={}
for j in range(nombre_plan):
print ('********************')
print (j,centre_plan[j][0],centre_plan[j][1],centre_plan[j][2])
print (j,normal_plan[j][0],normal_plan[j][1],normal_plan[j][2])
if j > 5:
indice_plan_bidon[j]=indice_pyra
else:
indice_plan_bidon[j]=indice_barre
vitesse_milieu_bidon[j]=vitesse_lumiere_vide/indice_plan_bidon[j]
sin_crit_bidon[j]=1./indice_plan_bidon[j]
return
######################################################
def distance(i,r,v):
method=1
if method == 0:
##### calcul distance minimu par rapport au plan ne depend du vecteur centre_plan si centre_plan est dans le plan considere
d_min=normal_plan[i][0]*(r[0]-centre_plan[i][0])+normal_plan[i][1]*(r[1]-centre_plan[i][1])+normal_plan[i][2]*(r[2]-centre_plan[i][2])
d_min=abs(d_min)
######### temps pour arriver d_min * projection de Vitesse le long de normal_plan
v_proj2=(normal_plan[i][0]*v[0])+(normal_plan[i][1]*v[1])+(normal_plan[i][2]*v[2])
v_proj=v_proj2
############ temps pour y arriver
t_inter=d_min/v_proj
########### point d interaction ###################
x_inter=r[0]+v[0]*t_inter
y_inter=r[1]+v[1]*t_inter
z_inter=r[2]+v[2]*t_inter
####################################################
if method == 1:
v_proj=(normal_plan[i][0]*v[0])+(normal_plan[i][1]*v[1])+(normal_plan[i][2]*v[2])
aux1=-normal_plan[i][0]*r[0]
aux2=-normal_plan[i][1]*r[1]
aux3=-normal_plan[i][2]*r[2]
aux4_centre=normal_plan[i][0]*centre_plan[i][0]+normal_plan[i][1]*centre_plan[i][1]+normal_plan[i][2]*centre_plan[i][2]
aux_v1=normal_plan[i][0]*v[0]
aux_v2=normal_plan[i][1]*v[1]
aux_v3=normal_plan[i][2]*v[2]
t_inter=(aux1+aux2+aux3+aux4_centre)/(aux_v1+aux_v2+aux_v3)
x_inter=r[0]+v[0]*t_inter
y_inter=r[1]+v[1]*t_inter
z_inter=r[2]+v[2]*t_inter
if method == 2:
v_proj=(normal_plan[i][0]*v[0])+(normal_plan[i][1]*v[1])+(normal_plan[i][2]*v[2])
aux=normal_plan[i][0]*(centre_plan[i][0]-r[0])
auy=normal_plan[i][1]*(centre_plan[i][1]-r[1])
auz=normal_plan[i][2]*(centre_plan[i][2]-r[2])
aux_v1=normal_plan[i][0]*v[0]
aux_v2=normal_plan[i][1]*v[1]
aux_v3=normal_plan[i][2]*v[2]
t_inter=(aux+auy+auz)/(aux_v1+aux_v2+aux_v3)
x_inter=r[0]+v[0]*t_inter
y_inter=r[1]+v[1]*t_inter
z_inter=r[2]+v[2]*t_inter
intersection=math.sqrt((x_inter-r[0])*(x_inter-r[0])+(y_inter-r[1])*(y_inter-r[1])+(z_inter-r[2])*(z_inter-r[2]))
#operation de symetrie cas speculaire
##################### traitement de la surface ##################
sin_inc=abs(v_proj/vitesse_lumiere)
if sin_inc < sin_crit:
speculaire=1.
transmission=1.
else:
speculaire=0.
transmission=1.
###################################################"
if speculaire == 1:
vx_inter= v[0]-2*scalaire(i,v)*normal_plan[i][0]
vy_inter= v[1]-2*scalaire(i,v)*normal_plan[i][1]
vz_inter= v[2]-2*scalaire(i,v)*normal_plan[i][2]
else:
diffu=0
while diffu==0:
ctet=-1+random()*2
stet=sqrt(1-ctet*ctet)
phi=2*math.pi*random()
#phi=math.pi*random()
vx=vitesse_lumiere*ctet
vy=vitesse_lumiere*stet*math.cos(phi)
vz=vitesse_lumiere*stet*math.sin(phi)
#print 'temps',t_inter
vitesse_photon=[vx,vy,vz]
if scalaire(i,vitesse_photon) > 0 :
diffu=1
vx_inter=vx
vy_inter=vy
vz_inter=vz
#print '~~~~~~~~~~~~~~~~~~~~~~'
#print 'ds distance','plan',i,'positions',x_inter,y_inter,z_inter,'temps',t_inter
#print '~~~~~~~~~~~~~~~~~~~~~~'
return [intersection,x_inter,y_inter,z_inter,vx_inter,vy_inter,vz_inter,transmission,t_inter,speculaire]
##########################################################
def muon(x_impact,y_impact,z_impact):
suite=True
while suite == True:
ran1=-1+2*random()
ran2=random()
if ran2 < ran1**2:
suite=False
# ran2 c'est cos(theta)
fac=sqrt(1-ran2**2)
delta_zmax=2.*Lx*fac
phi=2*math.pi*random()
position_muon=[x_impact,y_impact,z_impact]
vitesse_muon=[-vitesse_lumiere_vide*fac,vitesse_lumiere_vide*ran2*cos(phi),vitesse_lumiere_vide*ran2*sin(phi)]
#vitesse_muon=[-vitesse_lumiere_vide,0,0]
detect=False
while detect == False:
distance_limite=1.e5
for j in range(nombre_plan):
res=scalaire(j,vitesse_muon)
if res < 0:
sortie=distance(j,position_muon,vitesse_muon)
#print 'plan possible',j,'distance au plan',sortie[0]
if sortie[0] < distance_limite:
distance_limite=sortie[0]
plan_touche=j
x_plan=sortie[1]
y_plan=sortie[2]
z_plan= sortie[3]
vx_plan=sortie[4]
vy_plan=sortie[5]
vz_plan=sortie[6]
detect=True
t_muon=sortie[8]
#print 'point intersection muon',x_plan,y_plan,z_plan,'plan touche',j
#print 'distance',sortie[0]
#position muon c'est son point d'entree, egale a x_impact,etc... et la vitesse est celle echantillonee (aleatoirement), et le t_muon est le temps qu'il passe dans le detecteur
return [position_muon,vitesse_muon,t_muon,distance_limite]
######################################################################
def source(y_impact,z_impact):
x_new=position_courante[0]+vitesse_courante[0]*random()*t_inter_
y_new=position_courante[1]+vitesse_courante[1]*random()*t_inter_
z_new=position_courante[2]+vitesse_courante[2]*random()*t_inter_
#print x_new,y_new,z_new
ctet=-1+random()*2
stet=sqrt(1-ctet*ctet)
vz=vitesse_lumiere*ctet
vy=vitesse_lumiere*stet*cos(phi)
vx=vitesse_lumiere*stet*sin(phi)
vitesse_photon=[vx,vy,vz]
position_aleatoire=[x_new,y_new,z_new]
return [position_aleatoire,vitesse_photon,t_inter_,distance]
##############" un muon avec 90 pourcent des photons et des gammas avec 10 pourcent des photons
if random() < proba_muon:
### lumiere muon ########################
z_src=z_impact
y_src=y_impact
if abs(z_src) >= Lz:
x_src=(-1.+2.*random())*2.3
y_src=(-1.+2.*random())*2.3
else:
x_src=(-1+2*random())*Lx
#y_src=0
#z_src=((-delta_zmax/2./Lx)*x_src)+(delta_zmax/2.)
position_source=[x_src,y_src,z_src]
ctet=-1.+random()*2.
stet=sqrt(1-ctet*ctet)
phi=2*math.pi*random()
vz=vitesse_lumiere*ctet
vy=vitesse_lumiere*stet*cos(phi)
vx=vitesse_lumiere*stet*sin(phi)
vitesse_photon=[vx,vy,vz]
t_start=0.
return [position_source,vitesse_photon,t_start]
else:
#### lumiere gamma
x_src=(-1.+2.*random())*Lx
y_src=(-1.+2.*random())*Ly
z_src=(-1.+2.*random())*Lz
position_source=[x_src,y_src,z_src]
ctet=-1.+random()*2.
stet=sqrt(1-ctet*ctet)
phi=2*math.pi*random()
vz=vitesse_lumiere*ctet
vy=vitesse_lumiere*stet*cos(phi)
vx=vitesse_lumiere*stet*sin(phi)
vitesse_photon=[vx,vy,vz]
t_start=10.
return [position_source,vitesse_photon,t_start]
###########################
def prod_vec(I,J,K,L):
x=J[0]-I[0]
y=J[1]-I[1]
z=J[2]-I[2]
xp=L[0]-K[0]
yp=L[1]-K[1]
zp=L[2]-K[2]
aux=y*zp-yp*z
auy=z*xp-zp*x
auz=x*yp-xp*y
norm=sqrt(aux**2+auy**2+auz**2)
return [aux/norm,auy/norm,auz/norm]
##################################################
def analyse_temps(vecteur):
#########" recherche premier max liste a l 'endroit
yield_max=0
imax=-1000
for i in vecteur:
if vecteur[i] > yield_max:
yield_max=vecteur[i]
imax=i
t0=delta_sample*float(imax)
#### recherche minimum apes max
yield_max=0
imax=1000
t_delay=8
t_start=t0+t_delay
istart=int(t_start/delta_sample)
for i in range(istart,len(vecteur)):
if vecteur[i] > yield_max:
yield_max=vecteur[i]
imax=i
t1=delta_sample*float(imax)
######
return [t0,t1]
##########################################################################################################################
##########################################PROGRAMME PRINCIPAL#############################################################
##########################################################################################################################
if __name__=="__main__":
print ('debut du programme')
name_file='nebula_'+str(run)+'.root'
print (name_file)
anfile = TFile(name_file,"RECREATE","No comment")
ddd=TH1F('distance','distance',1500,0.,1500.)
x_y=TH2F('xy','xy',500,-10.,10.,500,-10.,10.)
x_z=TH2F('xz','xz',500,-120.,120.,500,-120.,120.)
y_z=TH2F('yz','yz',500,-120.,120.,500,-120.,120.)
e_t=TH2F('e_t','e_t',500,-30.,30.,100,0.,.02)
z_depart=TH1F('z0','z0',100,-100.,100.)
tof0=TH1F('temps t0','temps t0',1000,0.,50.)
tof1=TH1F('temps t1','temps t1',1000,0.,50.)
tof0_bis=TH1F('t0_2ns','t0_2ns',1500,0.,300.)
tof1_bis=TH1F('t1_2ns','t1_2ns',1500,0.,300.)
pos_x=TH1F('xini','xini',100,-1.5,1.5)
pos_y=TH1F('yini','yini',100,-1.5,1.5)
pos_z=TH1F('zini','zini',100,-1.5,1.5)
nhit0=TH1F('nhit0','nhit0',20,0,20)
nhit1=TH1F('nhit1','nhit1',20,0,20)
Q0=TH1F('Q0','Q0',200,0.,200.)
Q1=TH1F('Q1','Q1',200,0.,200.)
delta_z0=TH1F('delta_z0','delta_z0',400,-20.,20.)
delta_z1=TH1F('delta_z1','delta_z1',400,-20.,20.)
delta_t=TH1F('delta_t','delta_t',600,-30.,30.)
diff_time_weighted=TH1F('diff_time_weighted','diff_time_weighted',120,-15.,15.)
diff_time_unweighted=TH1F('diff_time_unweighted','diff_time_unweighted',120,-15.,15.)
det0_hit=TH1F('photon_0','photon_0',50,0.,50.)
det1_hit=TH1F('photon_1','photon_0',50,0.,50.)
Q0_t=TH2F('Q0_t','Q0_t',120,-15.,15.,100,0.,200.)
Q1_t=TH2F('Q1_t','Q1_t',120,-15.,15.,100,0.,200.)
Qtot_t=TH2F('Qtot_t','Qtot_t',120,-15.,15.,150,0.,300.)
bidim={}
trajet=TH3F('trajec','trajec',200,-100,100,120,-6.,6.,120,-6.,6.)
trajet_xy=TH2F('trajet_xy','trajet_xy',120,-30.,30.,120,-30.,30.)
trajet_yz=TH2F('trajet_yz','trajet_yz',120,-30.,30.,120,-30.,30.)
trajet_xz=TH2F('trajet_xz','trajet_xz',120,-30.,30.,120,-30.,30.)
for j in range(nombre_plan):
bidim[j]=TH2F(str(j),str(j),100,-10.,10.,100,-10.,10.)
hit_plan=TH1F('hit','hit',20,0.,20.)
print ('definition de la geometrie')
geometrie()
print ('fin de la geometrie')
#for j in range(nombre_plan):
#print j,centre_plan[j][0],centre_plan[j][1],centre_plan[j][2]
#for j in range(nombre_plan):
#print j,normal_plan[j][0],normal_plan[j][1],normal_plan[j][2]
out_ini={}
fac_correc=4. ### par comparaison avec donnees
tau_fluo=2.
tau_PM=4.
attenuation_lenght=380*fac_correc
distance_max= 20.*attenuation_lenght
detect_true=0
number_of_hit={}
for i in range(20):
number_of_hit[i]=0
n_detect0=0
n_detect1=0
time_hit0=[]
time_hit1=[]
weight_hit0=[]
weight_hit1=[]
energy_hit0=[]
energy_hit1=[]
proba_muon=1.
proba_gamma=1.-proba_muon
nombre_evt=int(nombre_muon*nombre_photon_par_muon)
nombre_photon_muon=nombre_evt*proba_muon
nombre_photon_gamma=nombre_evt*proba_gamma
quantum_muon=24./(nombre_photon_muon+.00001)
quantum_gamma=2./(nombre_photon_gamma+.00001)
z_geom=Lz+hauteur_pyramide
energy_muon_0=0
energy_gamma_0=0
energy_muon_1=0
energy_gamma_1=0
i_muon=0
muon_trigger=1
if muon_trigger == 1:
### trgger = 0 toute la barre
if trigger == 0:
z_fixe=(-1.+2.*random())*Lz
#### trigger = 1 TOP 90
if trigger == 1:
z_fixe=(-1.+2.*random())*1.
#### trigger = 2 close to det 0
if trigger == 2:
z_fixe=78+(-1.+2.*random())*5.
y_fixe=(-1.+2.*random())*Ly
x_fixe=Lx
else:
if random() < .5:
z_fixe=Lz+random()*hauteur_pyramide
else:
z_fixe=-Lz-random()*hauteur_pyramide
n0_mean=0
n1_mean=0
t0_mean=0
t1_mean=0
plan_touche=-100
w0_mean=0.
w1_mean=0.
weight=0
temps_parcours_final=0
# caracterise le muon
utile=muon(x_fixe,y_fixe,z_fixe)
weight=utile[3]*nombre_photon_par_cm/float(nombre_photon_par_muon)
##### liste des temps d arrivee pour chaque PM
##### echantillonnage en temps delta_sample= .5
delta_sample=.25
i_scan_max=250000
t_arrival_0={}
t_arrival_1={}
trigger_0=1000
trigger_1=1000
for i in range(i_scan_max):
t_arrival_0[i]=0
t_arrival_1[i]=0
for i in range(nombre_evt):
distance_parcourue=0.
transmission_totale=1
i_muon=i_muon+1
if i_muon <= nombre_photon_par_muon:
#out=source(y_fixe,z_fixe)
if plan_touche == 0:
if temps_parcours_final < trigger_0:
trigger_0=temps_parcours_final
i_arrival=int(temps_parcours_final/delta_sample)
t_arrival_0[i_arrival]=t_arrival_0[i_arrival]+weight
t0_mean=t0_mean+temps_parcours_final*weight
w0_mean=w0_mean+float(weight)
n0_mean=n0_mean+1
if plan_touche == 1:
if temps_parcours_final < trigger_1:
trigger_1=temps_parcours_final
i_arrival=int(temps_parcours_final/delta_sample)
t_arrival_1[i_arrival]=t_arrival_1[i_arrival]+weight
t1_mean=t1_mean+temps_parcours_final*weight
w1_mean=w1_mean+float(weight)
n1_mean=n1_mean+1
if i_muon == nombre_photon_par_muon:
i_muon=0
#### recherche du trigger == maximum de coups
###### recherche du premier maximum ##########################
pm_0=analyse_temps(t_arrival_0)
pm_1=analyse_temps(t_arrival_1)
print (pm_0,pm_1)
t0_0=pm_0[0]
t0_1=pm_0[1]
t1_0=pm_1[0]
t1_1=pm_1[1]
##### new event
for i in range(i_scan_max):
t_arrival_0[i]=0
t_arrival_1[i]=0
if muon_trigger == 1:
### trgger = 0 toute la barre
if trigger == 0:
z_fixe=(-1.+2.*random())*Lz
#### trigger = 1 TOP 90
if trigger == 1:
z_fixe=(-1.+2.*random())*1.
#### trigger = 2 close to det 0
if trigger == 2:
z_fixe=78+(-1.+2.*random())*5.
y_fixe=(-1.+2.*random())*Ly
x_fixe=Lx
utile=muon(x_fixe,y_fixe,z_fixe)
weight=float(utile[3])*float(nombre_photon_par_cm)/float(nombre_photon_par_muon)
else:
if random() < .5:
z_fixe=Lz+random()*hauteur_pyramide
else:
z_fixe=-Lz-random()*hauteur_pyramide
##### spectres finals
det0_hit.Fill(float(n0_mean),w0_mean)
det1_hit.Fill(float(n1_mean),w1_mean)
Q0.Fill(w0_mean)
Q1.Fill(w1_mean)
wtot=w0_mean+w1_mean
hit_tot=float(n0_mean)+random()+float(n1_mean)+random()
hit_left=float(n1_mean)+random()
#ddt=t0_0-t1_0
print (trigger_0,trigger_1)
ddt=trigger_0-trigger_1
temps_0=trigger_0+.15*(t0_0-trigger_0)
temps_1=trigger_1+.15*(t1_0-trigger_1)
ddt=t0_0-t1_0
#ddt=temps_0-temps_1
Q0_t.Fill(ddt,w0_mean)
Q1_t.Fill(ddt,w1_mean)
Qtot_t.Fill(ddt,wtot)
diff_time_unweighted.Fill(ddt)
if multi_hit != 0:
ddt=t0_0-t1_1
Q0_t.Fill(ddt,w0_mean)
Q1_t.Fill(ddt,w1_mean)
Qtot_t.Fill(ddt,wtot)
diff_time_unweighted.Fill(ddt)
ddt=t0_1-t1_0
Q0_t.Fill(ddt,w0_mean)
Q1_t.Fill(ddt,w1_mean)
Qtot_t.Fill(ddt,wtot)
diff_time_unweighted.Fill(ddt)
ddt=t0_1-t1_1
Q0_t.Fill(ddt,w0_mean)
Q1_t.Fill(ddt,w1_mean)
Qtot_t.Fill(ddt,wtot)
diff_time_unweighted.Fill(ddt)
t0_mean=0
t1_mean=0
w0_mean=0
w1_mean=0
n0_mean=0
n1_mean=0
plan_touche=-100
trigger_0=1000
trigger_1=1000
detect=0
n_hit=0
#### placer le photon a l aide de source()
#### utiliser utile pour fairele gamma
t_aleatoire=random()*utile[2]
x_ini=utile[0][0]+utile[1][0]*t_aleatoire
y_ini=utile[0][1]+utile[1][1]*t_aleatoire
z_ini=utile[0][2]+utile[1][2]*t_aleatoire
v_aux=utile[1][0]
v_auy=utile[1][1]
v_auz=utile[1][2]
utile_2=vitesse_gamma()
position_courante=[x_ini,y_ini,z_ini]
vitesse_courante=[utile_2[0],utile_2[1],utile_2[2]]
#######REMPLIR TRAJECTOIRE DANS ROOT#######
trajet.Fill(position_courante[2],position_courante[1],position_courante[0])
trajet_xy.Fill(v_aux,v_auy)
trajet_yz.Fill(v_auy,v_auz)
trajet_xz.Fill(v_aux,v_auz)
###########################################
t_creation=utile[2]
#print position_courante
#print vitesse_courante
while detect == 0:
distance_limite=1.e5
for j in range(nombre_plan):
res=scalaire(j,vitesse_courante)
#print 'scalaire',j,res
if res < 0:
sortie=distance(j,position_courante,vitesse_courante)
#print 'plan possible',j,'distance au plan',sortie[0]
if sortie[0] < distance_limite:
distance_limite=sortie[0]
plan_touche=j
x_plan=sortie[1]
y_plan=sortie[2]
z_plan= sortie[3]
vx_plan=sortie[4]
vy_plan=sortie[5]
vz_plan=sortie[6]
temps_parcours= distance_limite/vitesse_lumiere
transmission=sortie[7]
rebond_pm=sortie[9]
if j < 2:
transmission=1
#if plan_touche <= 1:
#print 'PM',plan_touche,rebond_pm
#print 'plan atteint',plan_touche,'x',x_plan,'y',y_plan,'z',z_plan,'vx',vx_plan,'vy_plan',vy_plan,'vz',vz_plan
number_of_hit[plan_touche]=number_of_hit[plan_touche]+1
kill_or_not=exp(-distance_limite/attenuation_lenght)
if random() > kill_or_not:
transmission=0
#sys.exit()
transmission_totale=transmission_totale*transmission
n_hit=n_hit+1
distance_parcourue=distance_parcourue + distance_limite
#transmit_pm=1. # .025
if plan_touche <= 1 and rebond_pm == 0:
detect = 1
if (plan_touche == 0 or plan_touche == 1 or distance_parcourue > distance_max or transmission == 0 and detect == 1):
#if plan_touche <= nombre_plan:
#print 'plan final',plan_touche
#detect=1
if (plan_touche == 0 or plan_touche == 1 and transmission_totale != 0):
#weight=exp(-distance_parcourue/attenuation_lenght)
#weight=weight*transmission_totale
#weight=1
z_depart.Fill(z_ini,weight)
detect_true=detect_true+weight
ddd.Fill(distance_parcourue)
x_y.Fill(x_ini,y_ini,weight)
x_z.Fill(x_ini,z_ini,weight)
y_z.Fill(y_ini,z_ini,weight)
#print plan_touche,x_plan,y_plan
bidim[plan_touche].Fill(x_plan,y_plan)
hit_plan.Fill(float(plan_touche))
temps_parcours_final=distance_parcourue/vitesse_lumiere
t_creation=-tau_fluo*log(random())-tau_PM*log(random())
temps_parcours_final=temps_parcours_final+t_creation
if (plan_touche==0 and transmission_totale != 0):
tof0.Fill(temps_parcours_final,weight)
tof0_bis.Fill(temps_parcours_final,weight)
nhit0.Fill(float(n_hit),weight)
time_hit0.append(temps_parcours_final)
weight_hit0.append(weight)
n_detect0=n_detect0+weight
if t_creation == 0:
energy_hit0.append(quantum_muon)
energy_muon_0=energy_muon_0+quantum_muon*weight
else:
energy_hit0.append(quantum_gamma)
energy_gamma_0=energy_gamma_0+quantum_gamma*weight
z_guess=vitesse_lumiere*temps_parcours_final
z_guess=z_guess-z_geom
delta_z0.Fill(z_guess,weight)
if (plan_touche==1 and transmission_totale != 0):
tof1.Fill(temps_parcours_final,weight)
tof1_bis.Fill(temps_parcours_final,weight)
nhit1.Fill(float(n_hit),weight)
time_hit1.append(temps_parcours_final)
weight_hit1.append(weight)
n_detect1=n_detect1+weight
if t_creation == 0:
energy_hit1.append(quantum_muon)
energy_muon_1=energy_muon_1+quantum_muon*weight
else:
energy_hit1.append(quantum_gamma)
energy_gamma_1=energy_gamma_1+quantum_gamma*weight
z_guess=vitesse_lumiere*temps_parcours_final
z_guess=z_guess-z_geom
delta_z1.Fill(z_guess,weight)
position_courante[0]=x_plan
position_courante[1]=y_plan
position_courante[2]=z_plan
vitesse_courante[0]=vx_plan
vitesse_courante[1]=vy_plan
vitesse_courante[2]=vz_plan
#print 'nouvelle',position_courante[0],position_courante[1],position_courante[2],vitesse_courante[0],vitesse_courante[1],vitesse_courante[2]
#print '#####################'
#print 'evt',i
print ('Nombre de photons propages',nombre_evt)
print ('nbe de photons detectes',detect_true,'efficacite totale',detect_true/float(nombre_evt))
print ('nbe de photons PM0',n_detect0,'efficacite',float(n_detect0)/float(nombre_evt))
print ('nbe de photons PM1',n_detect1,'efficacite',float(n_detect1)/float(nombre_evt))
print ('energie muon deposee dans plan 0',energy_muon_0 )
print ('energie muon deposee dans plan 1',energy_muon_1)
print ('energie gamma deposee dans plan 0',energy_gamma_0 )
print( 'energie gamma deposee dans plan 1',energy_gamma_1)
for i in range(nombre_plan):
print (i,number_of_hit[i])
print ('distrbution des temps')
#print time_hit0
#print time_hit1
print ('distrbution des poids')
#print weight_hit0
#print weight_hit1
i0max=len(time_hit0)
i1max=len(time_hit1)
imax=0
if i1max > imax :
i1max=imax
if i0max > imax :
i0max=imax
for i in range(i0max):
t_start=0
t0=time_hit0[i]+t_start
w0=weight_hit0[i]
e0=energy_hit0[i]
for j in range(i1max):
t_start=0
t1=time_hit1[j]+t_start
w1=weight_hit1[j]
delta=t0-t1
w_coinc=w0*w1
energie=e0+energy_hit1[i]
delta_t.Fill(delta,w_coinc)
e_t.Fill(delta,energie,w_coinc)
#write_calcul.close()
anfile.Write()
anfile.Close()
sys.exit(0)
#
Plot de l'animation