Zufallszahlen

1. Erzeugen von Zufallszahlen Ein gängiges Verfahren zur Erzeugung von Pseudo-Zufallszahlen ist die lineare Kongruenzmethode:
   
   $${\cal{I}}_{n+1} = (a \cdot {\cal{I}}_{n} + c ) \; mod \; m$$
   
   
   Grundsätzlich sind solche Folgen periodisch, die Periodenlänge ist im besten Falle $m$. Entscheidend für die Güte des Generators ist die Wahl der Parameter $a$ und $b$, sowie des `Seed'-Werts ${\cal{I}}_{0}$.

   Im folgenden sollen die Eigenschaften von Generatoren mit verschiedenen Parametersets untersucht werden.

   Dazu benötigt man zunächst einen Zufallsgenerator, den man einfach als C++ Klasse implementieren kann (file MyRnd.C).

   Die Implementation des Zufallsgenerators als C++ Klasse hat wesentliche Vorteile gegenüber einer konventionellen Implementation als Funktion: alle Instanzen von MyRnd sind völlig unabhängig, d.h. man kann beliebig viele Zufallsgeneratoren mit unterschiedlichen Seeds und Parametern verwenden.

   Einfache Tests von Zufallsgeneratoren sind die Überprüfung der Gleichverteilung in 1, 2 oder 3 Dimensionen.

   Beispielprogramm tRnd.C:

   Aufgabe a:
   Testen Sie verschiedene Parameter für den Zufallsgenerator:
   (a) $a=137,\; c = 187,\; m = 256$
   (b) $a=193,\; c = 73,\; m = 256$
   (c) $a=65539,\; c = 0,\; m = 2147483648$
   

   Als seed sollte man eine Zahl $> 0$ verwenden.

   (a) und (b) sind `kurze' Generatoren. Sie nutzen die Länge $m$ voll aus, zeigen aber deutliche Unterschiede in ihrer Qualität in 2 D.
   (c) ist der Generator randu von IBM, der in den 60er und 70er Jahren weitverbreitet war. Er hat zwar eine hohe Periode ($2^{29}$), aber $a = 2^{16}+3$ ist lausig gewählt, schon in 3 D zeigen sich deutliche Abweichungen von der Gleichverteilung. Viele wissenschaftliche Arbeiten haben darunter gelitten ...

   Zum Vergleich können Sie auch den ROOT Standard Generator mit gRandom->Rndm() betrachten.

2. Monte Carlo Berechung der Kreiszahl $\pi$
   Berechnen sie die Kreiszahl pi mit Hilfe der Monte Carlo Integration wie in der Vorlesung gesprochen.
   • Verwenden sie hierzu die Tropfenmethode. Erzeugen sie zwei unabhnigige gleichförming verteilte Zufallszahlen mit:
   • TRandom ran
   • ran.SetSeed("1'')
   • x = ran.Uniform()
   • Testen sie ob $r^2=x^2+y^2$ innerhalb oder außhalb der Kreisfläche liegt
   • Berechnen sie $\pi$ für 100, 1000 und 10000 Tropfen und den zugehörigen Fehler.
   • Füllen sie die erzeugten Zufallszahlen in ein Histogramm histo = TH2F('test' ,'Berechnung von #pi', 100,-0.1,1.1, 100,-0.1,1.1 )

   • Lösungsbeispiel: pi.C, pi.py
   

3. Simulation des radioaktiven Zerfalls
   Mit gleichverteilten Zufallszahlen lässt sich der radioaktive Zerfall leicht simulieren, die Wahrscheinlichkeit daß ein Kern zerfällt ist
   
   $$p = \lambda \Delta t$$
   
   
   In tdecay.C wird ein Zerfallsexperiment simuliert, wobei man die Zerfallskonstante $\lambda$, das Zeitintervall $dt$, die Gesamtzeit $t$ und die Zahl der Nukleonen $n0$ zu Beginn vorgeben kann.

   Testen Sie zunächst das Programm für verschiedene Parameter. Dann ändern Sie das Programm so ab, daß Sie viele solcher Experimente durchführen und zusätzlich die Zahl der zerfallenen Teilchen für einige Zeitintervalle histogrammieren.

   
   // simulate radioactive decay
    void tdecay( Float_t lambda = 0.5, Float_t dt = 0.01, Float_t t = 20,
         Int_t n0 = 1000 )
    {
      TH1F * hdec = new TH1F("hdec", "radioactive decay", t/dt, 0, t);    
      // loop over time
      for ( Float_t tnow = 0.; tnow < t; tnow += dt ) {
        // loop over remaining particles
        Int_t ndec = 0;
        for ( int ie = 0; ie < n0; ie++ ) {
          if ( gRandom->Rndm() < dt*lambda ) {
            ndec ++;
            hdec->Fill(tnow);
        }}
        n0 -= ndec; // subtract decayed particles
      }
    }
   
   

4. Zusatzaufgaben:
   Wie sieht die analytische Transformation aus, um Zufallszahlen gemäss $f(x) = 2 \cdot x$ bzw. $e^{-x}$ zu erzeugen ? Probieren Sie es in ROOT.

   Erzeugen Sie eine beliebige Zufallsverteilung nach dem Hit&Miss Verfahren, z.B. für die Breit-Wigner Funktion (BreitWig.C).

   Hinweis: Funktionen als Macros
   Neben der direkten Art Funktionen zu definieren, wie in
   TF1 *f1 = new TF1("f1","sin(x)", 0, 10 );
   kann man auch eine Funtkion als Adresse angeben, z.B.:

   
   // file BreitWig.C
    Double_t BreitWig( Double_t *x, Double_t *par) 
    {
      // Breit- Wigner function
      Double_t mw = par[0], gw = par[1], mw2, gw2, eb2;
      mw2 = mw*mw;
      gw2 = gw*gw;
      eb2 = x[0]*x[0];
      return( gw2*mw2 / ( pow( eb2 - mw2, 2 ) + mw2 * gw2 ) );
    }
    root [0] .L BreitWig.C 
    root [3] TF1 *f1 = new TF1("f1",BreitWig,50,100,2)
    root [4] f1->SetParameters( 80., 2.);    
    root [5] f1->Draw();    
   
   

   Als 5. Argument von TF1(..) muss die Zahl der Parameter stehen. Die Parameter können dann per f1->SetParameters() gesetzt werden.
   Zufallszahlen die gemäss einer Funktion verteilt sind können in ROOT mit der Methode TF1::GetRandom() erzeugt werden, also z.B.:
   root [5] f1->GetRandom()