Changeset 1352 for trunk/WuerzburgSoft/Thomas/mphys/phys.C

Timestamp:

06/10/02 08:49:28 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/phys.C (modified) (12 diffs)

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -6 +6 @@
 #include <TF1.h>
 #include <TH1.h>
+#include <TH2.h>
 #include <TList.h>
 #include <TRandom.h>
@@ -689 +690 @@
 void DoIt()
 {
-    Double_t rcygnusx3 = 100; //30/3.258; // [~10kpc]
+    Double_t rcygnusx3 = 5000; //100; //30/3.258; // [~10kpc]
 
     Double_t startz = ZofR(&rcygnusx3);
@@ -697 +698 @@
 
     //Double_t nphot = 50;
-    Double_t runtime = 18*60*60; // [s]
+    Double_t runtime = 15*60; ///*18*60*/60; // [s]
 
     Double_t lo = 2e4;
@@ -717 +718 @@
 
     gStyle->SetOptStat(10);
+
+    TH2D position;
+    TH2D angle;
+    TH1D histpos;
+    TH1D hista;
+
+    MBinning binspolx;
+    MBinning binspoly;
+    binspolx.SetEdges(32, -180, 180);
+    binspoly.SetEdgesLog(20, 1e-9, 1e-5);
+    MH::SetBinning(&position, &binspolx, &binspoly);
+    MH::SetBinning(&angle,    &binspolx, &binspoly);
+    MH::SetBinning(&histpos,  &binspoly);
+    MH::SetBinning(&hista,    &binspoly);
+
+    hista.SetTitle("Particle Angle");
+    angle.SetTitle("Particle Angle");
+    histpos.SetTitle("Particle Position");
+    position.SetTitle("Particle Position");
 
     TH1D hist;
@@ -735 +755 @@
     MH::SetBinning(&histsrc, &bins);
 
-    TH1D hista;
-    MBinning binsa;
-    binsa.SetEdgesLog(16, 1e-12, 1e-8);
-    MH::SetBinning(&hista, &binsa);
-
     MH::MakeDefCanvas();
 
@@ -791 +806 @@
             while ((p=(MPhoton*)NextP()))
             {
-                //
-                // Calculate way until interaction takes place
-                //
-                Double_t z = p->GetZ();
-                Double_t l = rand.Exp(p->GetInteractionLength());
-                Double_t r = RofZ(&z);
-
-                //
-                // If photon went along us...  (l==0: infinite)
-                //
-                if (r<l || l==0)
+                if (!p->SetNewPosition())
                 {
-                    cout << (l==0?">":"!") << flush;
+                    cout << "!" << flush;
+
                     hist.Fill(p->GetEnergy(), p->GetEnergy());
+                    position.Fill(p->GetPhi()*kRad2Deg, p->GetR());
+                    angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg);
+                    histpos.Fill(p->GetR());
+                    hista.Fill(p->GetTheta()*kRad2Deg);
+
                     fcas << p->GetEnergy() << endl;
+
                     listg.Remove(p);
                     continue;
                 }
-
-                //
-                // Set new z value
-                //
-                r -= l;
-                z = ZofR(&r);
-                p->SetZ(z);
 
                 //
@@ -822 +827 @@
                 MPhoton photon;
 
-                phot.SetParameter(0, z);
+                phot.SetParameter(0, p->GetZ());
                 photon.SetEnergy(pow(10, phot.GetRandom()));
-                photon.SetAngle(rand.Uniform(TMath::Pi() * 2));
-
-                if (!pair.Process(p, &photon, &liste))
+                Double_t theta = rand.Uniform(TMath::Pi() * 2);
+
+                if (!pair.Process(p, &photon, theta, &liste))
                 {
                     cout << "0" << flush;;
@@ -832 +837 @@
                 }
 
+                listg.Remove(p);
                 cout << "." << flush;
-
-                listg.Remove(p);
             }
 
@@ -847 +851 @@
             {
                 cout << ":" << flush;
-
-                hista.Fill(fabs(e->GetAngle()));
-
                 while(1)
                 {
-                    Double_t E = e->GetEnergy();
-                    Double_t z = e->GetZ();
-                    Double_t l = rand.Exp(e->GetInteractionLength());
-                    Double_t r = RofZ(&z) - l;
-
-                    if (r<0)
+                    if (!e->SetNewPosition())
                     {
                         cout << "!" << flush;
@@ -863 +859 @@
                     }
 
-                    z = ZofR(&r);
-                    e->SetZ(z);
-
-                    Double_t ep = e->GetEnergyLoss();
-
-                    if (E-ep<lo*5)
+                    MPhoton *p = e->DoInvCompton();
+                    if (e->GetEnergy()<lo*5)
                     {
                         cout << "0" << flush;
+                        delete p;
                         break;
                     }
 
-                    e->SetEnergy(E-ep);
-
-                    if (ep<lo)
+                    if (p->GetEnergy()<lo)
                     {
                         cout << "i" << flush; // ignored
+                        delete p;
                         continue;
                     }
-                    /*
-                    Double_t eps = 1e-11; // photon vorher
-
-                    Double_t E0 = 511e-6;
-                    Double_t gamma = E/E0;
-                    Double_t gamma2 = gamma*gamma;
-                    Double_t beta = sqrt(1-1/gamma2);
-
-                    Double_t theta = rand.Uniform(TMath::Pi()*2);
-
-                    Double_t p3 = gamma2 * (1-cos(theta)) - ep/E0;
-
-                    Double_t a  = 1- (1 + ep/(eps*p3));
-
-                    cout << "/" << gamma << "," << p3 << "," << ep << "," << a;
-                    cout << "(" << 90-180*(TMath::Pi()-acos(a))/TMath::Pi() <<")" << flush;
-                    */
-
-                    MPhoton *p = new MPhoton(ep, z);
+
                     listg.Add(p);
-
                     cout << "." << flush;
                 }
@@ -912 +885 @@
         if (now!=timer || n<10)
         {
-            histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
+            histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz,
+                                  (int)runtime/60, (int)runtime%60, hist.GetEntries()));
             gPad->Modified();
             gPad->Update();
@@ -943 +917 @@
 
     MH::MakeDefCanvas();
-    hista.SetXTitle("[rad]");
+    position.SetXTitle("Pos: \\Phi [\\circ]");
+    position.SetYTitle("Pos: R [kpc]");
+    position.DrawCopy("surf1pol");
+    gPad->SetLogy();
+
+    MH::MakeDefCanvas();
+    angle.SetXTitle("Angle: \\Psi [\\circ]");
+    angle.SetYTitle("Angle: \\Theta [\\circ]");
+    angle.DrawCopy("surf1pol");
+    gPad->SetLogy();
+
+    MH::MakeDefCanvas();
+    histpos.SetXTitle("R [kpc]");
+    histpos.SetYTitle("Counts");
+    histpos.DrawCopy();
+    gPad->SetLogx();
+
+    MH::MakeDefCanvas();
+    hista.SetXTitle("\\Theta [\\circ]");
+    hista.SetYTitle("Counts");
     hista.DrawCopy();
     gPad->SetLogx();