Added to repository
This commit is contained in:
paraiso
364
geant4/LEMuSR/MEYER/meyer.for
Normal file
364
geant4/LEMuSR/MEYER/meyer.for
Normal file
@ -0,0 +1,364 @@
|
||||
|
||||
c-------------------------------------------------------------------------------
|
||||
c Konstanten und Variable fuer Berechnung der Winkelaufstreuung in Triggerfolie
|
||||
c mittels Meyer-Formel (L.Meyer, phys.stat.sol. (b) 44, 253 (1971)):
|
||||
|
||||
real g1, g2 ! Tabellierte Funktionen der Referenz
|
||||
real effRedThick ! effektive reduzierte Dicke ('tau' der Referenz)
|
||||
|
||||
|
||||
c - Parameter:
|
||||
|
||||
real Z1, Z2 ! die atomaren Nummern von Projektil und Target
|
||||
real a0 ! Bohrscher Radius in cm
|
||||
real screeningPar ! Screeningparameter 'a' in cm fuer Teilchen der
|
||||
! Kernladungszahl Z1=1 in Kohlenstoff (Z2 = 6)
|
||||
! bei Streichung von Z1 (vgl. Referenz, S. 268)
|
||||
|
||||
real r0Meyer ! r0(C) berechnet aus dem screeningParameter 'a'
|
||||
! und dem ebenfalls bei Meyer angegebenem
|
||||
! Verhaeltnis a/r0=0.26 (vgl. Referenz, S. 263 oben)
|
||||
real eSquare ! elektrische Ladung zum Quadrat in keV*cm
|
||||
real HWHM2sigma ! Umrechnungsfaktor von (halber!) Halbwertsbreite
|
||||
! nach Sigma der Gaussfunktion
|
||||
|
||||
real Na ! die Avogadrokonstante
|
||||
real mMolC ! molare Masse von C in ug
|
||||
real Pi ! die Kreiszahl
|
||||
|
||||
parameter (Z1 = 1, Z2 = 6, a0 = 5.29E-9, ScreeningPar = 2.5764E-9)
|
||||
parameter (r0Meyer = 9.909E-9, eSquare = 1.44E-10, HWHM2sigma = 1./1.17741)
|
||||
parameter (Na = 6.022e23, mMolC = 12.011e6, Pi = 3.141592654)
|
||||
|
||||
|
||||
c - Bei der Berechnung von Sigma auftretende Vorfaktoren.
|
||||
c (Meyer_faktor 1 wird benoetigt fuer Berechnung der reduzierten Dicke aus der
|
||||
c 'ug/cm2'-Angabe der Foliendicke. Meyer_faktor2 und Meyer_faktor3 werden
|
||||
c direkt fuer die Berechnung von sigma aus den beiden tabellierten Funktionen
|
||||
c g1 und g2 verwendet):
|
||||
|
||||
real Meyer_Faktor1, Meyer_Faktor2, Meyer_Faktor3
|
||||
|
||||
parameter (Meyer_faktor1 = Pi*screeningPar*screeningPar * Na/mMolC)
|
||||
! Na/mMolC = 1/m(C-Atom)
|
||||
parameter (Meyer_faktor2 = (2*Z1*Z2 * eSquare)/ScreeningPar * 180./Pi
|
||||
+ * HWHM2sigma)
|
||||
parameter (Meyer_faktor3 = (screeningPar/r0Meyer) * (screeningPar/r0Meyer))
|
||||
|
||||
|
||||
c-------------------------------------------------------------------------------
|
||||
c Kommentar zur Berechnung der Winkelaufstreuung nach Meyer:
|
||||
c
|
||||
c Als Bedingung fuer die Gueltigkeit der Rechnung wird verlangt, dass
|
||||
c
|
||||
c (1) die Anzahl n der Stoesse >> 20*(a/r0)^(4/3) sein muss. Fuer Protonen auf
|
||||
c Graphit ist laut Referenz a/r0 gleich 0.26 (mit Dichte von 3.5 g/ccm habe
|
||||
c ich einen Wert von 0.29 abgeschaetzt). Fuer Myonen hat man den selben
|
||||
c Wert zu nehmen. Damit ergibt sich die Forderung, dass n >> 3.5 sein muss.
|
||||
c
|
||||
c (2) unabhaengig von (1) n >> 5 sein muss, was (1) also mit einschliesst.
|
||||
c
|
||||
c Mit n = Pi*r0*r0*Teilchen/Flaeche ergibt sich fuer eine Foliendicke von
|
||||
c 3 ug/cm^2 als Abschaetzung fuer n ein Wert von 37. (r0 ueber r0 = 0.5 N^(1/3)
|
||||
c und 3.5 g/ccm zu 8.9e-9 cm abgeschaetzt). D.h., dass die Bedingungen in
|
||||
c unserem Fall gut erfuellt sind.
|
||||
c In dem Paper wird eine Formel fuer Halbwertsbreiten angegeben. Ich habe nicht
|
||||
c kontrolliert, in wie weit die Form der Verteilung tatsaechlich einer Gauss-
|
||||
c verteilung entspricht. Zumindest im Bereich der Vorwaertsstreuung sollte
|
||||
c die in diesem Programm verwendete Gaussverteilung aber eine sehr gute
|
||||
c Naeherung abgeben. Abweichungen bei groesseren Winkeln koennten jedoch u. U.
|
||||
c die absolute Streuintensitaet in Vorwaertsrichtung verfaelschen.
|
||||
|
||||
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
||||
c HIER GEHT DER PROGRAMMTEXT RICHTIG LOS
|
||||
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
c===============================================================================
|
||||
|
||||
options /extend_source
|
||||
|
||||
subroutine Get_F_Function_Meyer(tau,Ekin)
|
||||
c =========================================
|
||||
|
||||
implicit none
|
||||
|
||||
real tau
|
||||
real Ekin
|
||||
|
||||
real thetaSchlange,thetaSchlangeMax
|
||||
real theta,thetaMax,thetaStep
|
||||
real f1,f2,F
|
||||
|
||||
|
||||
c------------------------------------
|
||||
c - Parameter:
|
||||
|
||||
real Z1, Z2 ! die atomaren Nummern von Projektil und Target
|
||||
c real a0 ! Bohrscher Radius in cm
|
||||
real screeningPar ! Screeningparameter 'a' in cm fuer Teilchen der
|
||||
! Kernladungszahl Z1=1 in Kohlenstoff (Z2 = 6)
|
||||
! bei Streichung von Z1 (vgl. Referenz, S. 268)
|
||||
|
||||
real r0Meyer ! r0(C) berechnet aus dem screeningParameter 'a'
|
||||
! und dem ebenfalls bei Meyer angegebenem
|
||||
! Verhaeltnis a/r0=0.26 (vgl. Referenz, S. 263 oben)
|
||||
real eSquare ! elektrische Ladung zum Quadrat in keV*cm
|
||||
|
||||
real Pi ! die Kreiszahl
|
||||
|
||||
c parameter (a0 = 5.29E-9)
|
||||
parameter (Z1 = 1, Z2 = 6, ScreeningPar = 2.5764E-9)
|
||||
parameter (r0Meyer = 9.909E-9, eSquare = 1.44E-10)
|
||||
parameter (Pi = 3.141592654)
|
||||
|
||||
real Meyer_Faktor3
|
||||
real Meyer_Faktor4
|
||||
real zzz ! 'Hilfsparameter'
|
||||
real Meyer_Faktor5
|
||||
|
||||
parameter (Meyer_faktor3 = (screeningPar/r0Meyer) * (screeningPar/r0Meyer))
|
||||
parameter (Meyer_faktor4 = screeningPar / (2.*Z1*Z2*eSquare) * Pi/180.)
|
||||
parameter (zzz = screeningPar / (2.*Z1*Z2*eSquare))
|
||||
parameter (Meyer_faktor5 = zzz*zzz / (8*Pi*Pi))
|
||||
|
||||
c------------------------------------
|
||||
|
||||
integer nBin,nBinMax
|
||||
parameter (nBinMax=201)
|
||||
real value(0:nBinMax) /0.,nBinMax*0./
|
||||
real area(nBinMax) / nBinMax*0./
|
||||
real integ(0:nBinMax) /0.,nBinMax*0./
|
||||
common /MeyerTable/ value,area,integ,thetaStep,nBin
|
||||
|
||||
integer i
|
||||
real rhelp
|
||||
|
||||
integer HB_memsize
|
||||
parameter(HB_memsize=500000)
|
||||
real memory(HB_memsize)
|
||||
COMMON /PAWC/ memory
|
||||
|
||||
|
||||
c nur noch fuer Testzwecke:
|
||||
|
||||
real fValues(203)
|
||||
real fValuesFolded(203)
|
||||
|
||||
integer idh
|
||||
parameter (idh = 50)
|
||||
|
||||
INCLUDE 'mutrack$sourcedirectory:COM_DIRS.INC'
|
||||
character filename*20 ! Name der Ausgabe-Dateien
|
||||
COMMON /filename/ filename
|
||||
|
||||
c-------------------------------------------------------------------------------
|
||||
|
||||
c Festlegen des maximalen Theta-Wertes sowie der Schrittweite:
|
||||
|
||||
if (tau.LT.0.2) then
|
||||
write(*,*) 'Subroutine ''Get_F_Function_Meyer'':'
|
||||
write(*,*) 'Effektive Dicke ist kleiner als 0.2 => kann ich nicht ... => STOP'
|
||||
call exit
|
||||
elseif (tau.LE.2.) then
|
||||
! => Tabelle A
|
||||
thetaSchlangeMax = 4.0
|
||||
elseif (tau.LE.8.) then
|
||||
! => Tabelle B
|
||||
thetaSchlangeMax = 7.0
|
||||
elseif (tau.LE.20.) then
|
||||
! => Tabelle C
|
||||
thetaSchlangeMax = 20.0
|
||||
else
|
||||
write(*,*) 'Subroutine ''Get_F_Function_Meyer'':'
|
||||
write(*,*) 'Effektive Dicke ist groesser als 20 => kann ich nicht ... => STOP'
|
||||
call exit
|
||||
endif
|
||||
|
||||
thetaMax = thetaSchlangeMax / Meyer_Faktor4 / Ekin
|
||||
if (thetaMax.GT.50) then
|
||||
thetaStep = .5
|
||||
elseif (thetaMax.GT.25) then
|
||||
thetaStep = .25
|
||||
elseif (thetaMax.GT.12.5) then
|
||||
thetaStep = .125
|
||||
else
|
||||
thetaStep = .0625
|
||||
endif
|
||||
|
||||
|
||||
c Tabelle der F-Werte erstellen:
|
||||
|
||||
nBin = 0
|
||||
do theta = thetaStep, thetaMax, thetaStep
|
||||
|
||||
! Berechne aus theta das 'reduzierte' thetaSchlange (dabei gleich
|
||||
! noch von degree bei theta in Radiant bei thetaSchlange umrechnen):
|
||||
|
||||
thetaSchlange = Meyer_faktor4 * Ekin * theta
|
||||
|
||||
! Auslesen der Tabellenwerte fuer die f-Funktionen:
|
||||
|
||||
call F_Functions_Meyer(tau,thetaSchlange,f1,f2)
|
||||
if (thetaSchlange.EQ.-1) then
|
||||
! wir sind jenseits von thetaSchlangeMax
|
||||
goto 10
|
||||
endif
|
||||
|
||||
! Berechnen der Streuintensitaet:
|
||||
F = Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2)
|
||||
|
||||
nBin = nBin + 1
|
||||
if (nBin.GT.nBinMax) then
|
||||
write(*,*) 'nBin > nBinMax => EXIT'
|
||||
call exit
|
||||
endif
|
||||
value(nBin) = sind(theta)*F
|
||||
|
||||
fValues(nBin+1) = F ! fuer Testzwecke
|
||||
fValuesFolded(nBin+1) = sind(theta)*F ! fuer Testzwecke
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
c Berechnen der Flaecheninhalte der einzelnen Kanaele sowie der Integrale:
|
||||
|
||||
10 do i = 1, nBin
|
||||
area(i) = (value(i)+value(i-1))/2. * thetaStep
|
||||
integ(i) = integ(i-1) + area(i)
|
||||
enddo
|
||||
|
||||
|
||||
c Normiere totale Flaeche auf 1:
|
||||
|
||||
rHelp = integ(nBin)
|
||||
do i = 1, nBin
|
||||
value(i) = value(i) / rHelp
|
||||
area(i) = area(i) / rHelp
|
||||
integ(i) = integ(i) / rHelp
|
||||
enddo
|
||||
|
||||
|
||||
c vorerst noch: gib Tabelle in Datei und Histogrammfile aus:
|
||||
|
||||
! Berechne die Werte fuer theta=0:
|
||||
|
||||
call F_Functions_Meyer(tau,0.,f1,f2)
|
||||
F = Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2)
|
||||
fValues(1) = F
|
||||
fValuesFolded(1) = 0.
|
||||
|
||||
! Gib die Werte in das Tabellenfile aus:
|
||||
|
||||
c theta = 0.
|
||||
c open (10,file=outDir//':'//filename//'.TAB',status='NEW')
|
||||
c do i = 1, nBin+1
|
||||
c write(10,*) theta, fValues(i), fValuesFolded(i)
|
||||
c theta = theta + thetaStep
|
||||
c enddo
|
||||
c close (10)
|
||||
|
||||
|
||||
! Buchen und Fuellen der Histogramme:
|
||||
|
||||
call HBOOK1(idh,'F',nBin+1,-0.5*thetaStep,(real(nBin)+0.5)*thetaStep,0.)
|
||||
call HPAK(idh,fValues)
|
||||
call HRPUT(idh,outDir//':'//filename//'.RZ','N')
|
||||
call HDELET(idh)
|
||||
|
||||
call HBOOK1(idh+1,'F*sin([q])',nBin+1,-0.5*thetaStep,(real(nBin)+0.5)*thetaStep,0.)
|
||||
call HPAK(idh+1,fValuesFolded)
|
||||
call HRPUT(idh+1,outDir//':'//filename//'.RZ','U')
|
||||
call HDELET(idh+1)
|
||||
|
||||
|
||||
END
|
||||
|
||||
|
||||
c===============================================================================
|
||||
|
||||
options /extend_source
|
||||
|
||||
subroutine throwMeyerAngle (theta)
|
||||
c ==================================
|
||||
|
||||
implicit none
|
||||
|
||||
real lowerbound,y1,y2,f,root,radiant,fraction
|
||||
integer bin,nBin
|
||||
integer nBinMax
|
||||
parameter (nBinMax=201)
|
||||
|
||||
real theta,thetaStep
|
||||
real value(0:nBinMax) /0.,nBinMax*0./
|
||||
real area(nBinMax) / nBinMax*0./
|
||||
real integ(0:nBinMax) /0.,nBinMax*0./
|
||||
common /MeyerTable/ value,area,integ,thetaStep,nBin
|
||||
|
||||
real rhelp
|
||||
|
||||
real random
|
||||
integer seed
|
||||
common /seed/ seed
|
||||
|
||||
|
||||
c bin: Nummer des Bins, innerhalb dessen das Integral den Wert von
|
||||
c random erreicht oder ueberschreitet:
|
||||
|
||||
random = ran(seed)
|
||||
|
||||
bin = 1
|
||||
do while (random.GT.integ(bin))
|
||||
bin = bin + 1
|
||||
if (bin.GT.nBin) then
|
||||
write(*,*) 'error 1'
|
||||
call exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
fraction = (random-integ(bin-1)) / (integ(bin)-integ(bin-1))
|
||||
y1 = value(bin-1)
|
||||
y2 = value(bin)
|
||||
f = thetaStep / (y2-y1)
|
||||
rHelp = y1*f
|
||||
|
||||
radiant = rHelp*rHelp + fraction*thetaStep*(y1+y2)*f
|
||||
root = SQRT(radiant)
|
||||
lowerBound = real(bin-1)*thetaStep
|
||||
if (f.GT.0) then
|
||||
theta = lowerBound - rHelp + root
|
||||
else
|
||||
theta = lowerBound - rHelp - root
|
||||
endif
|
||||
|
||||
|
||||
END
|
||||
|
||||
|
||||
c===============================================================================
|
||||
|
||||
options /extend_source
|
||||
|
||||
subroutine F_Functions_Meyer(tau,thetaSchlange,f1,f2)
|
||||
c =====================================================
|
||||
|
||||
implicit none
|
||||
|
||||
c Diese Routine gibt in Abhaengigkeit von 'thetaSchlange' und 'tau'
|
||||
c Funktionswerte fuer f1 und f2 zurueck. f1 und f2 entsprechen dabei den
|
||||
c bei Meyer angegebenen Funktion gleichen Namens. Die in dieser Routine
|
||||
c verwendeten Tabellen sind eben dieser Referenz entnommen:
|
||||
c L.Meyer, phys.stat.sol. (b) 44, 253 (1971)
|
||||
|
||||
real tau,thetaSchlange
|
||||
real f1, f2, f1_(2), f2_(2)
|
||||
|
||||
integer column_,column,row
|
||||
|
||||
integer iColumn
|
||||
real weightCol, weightRow
|
||||
|
||||
c-------------------------------------------------------------------------------
|
||||
Reference in New Issue
Block a user