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-------------------------------------------------------------------------------