Changeset 1352 for trunk/WuerzburgSoft/Thomas/mphys

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (added)
• MElectron.cc (modified) (4 diffs)
• MElectron.h (modified) (3 diffs)
• MGenIRPhoton.cc (modified) (2 diffs)
• MGenPrimaryParticle.cc (modified) (1 diff)
• MPairProduction.cc (modified) (6 diffs)
• MPairProduction.h (modified) (1 diff)
• MParticle.cc (modified) (2 diffs)
• MParticle.h (modified) (2 diffs)
• MPhoton.h (modified) (1 diff)
• phys.C (modified) (12 diffs)

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -34 +34 @@
 #include <TF1.h>
 
+#include "MPhoton.h"
+
 ClassImp(MElectron);
 
@@ -238 +240 @@
 
 
-Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k = NULL)
+Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
 {
     Double_t E = x[0];
@@ -255 +257 @@
     Double_t e = pow(10, fP.GetRandom());
 
+    if (ep)
+        *ep = e;
+
     Double_t omega = e*E/E0/E0;
     Double_t Gamma = 4.*omega;
@@ -269 +274 @@
 }
 
-Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z)
+Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep)
 {
     return EnergyLoss(&E, &z);
 }
 
-Double_t MElectron::GetEnergyLoss() const
-{
-    return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ);
-}
+Double_t MElectron::GetEnergyLoss(Double_t *ep) const
+{
+    return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep);
+}
+
+#include <iostream.h>
+
+MPhoton *MElectron::DoInvCompton()
+{
+    Double_t E0 = 511e-6; //[GeV]
+
+    Double_t epsilon;
+    Double_t e = GetEnergyLoss(&epsilon);
+    /*
+     Double_t gamma = fEnergy/E0;
+     Double_t beta = sqrt(1-1/(gamma*gamma));
+     Double_t bega = sqrt(gamma*gamma-1);
+
+     Double_t theta = TMath::Pi()/4;
+
+     // er: photon energy before interaction rest frame
+     // e:  photon energy after interaction lab
+     Double_t er = gamma*epsilon*(1-beta*cos(theta)); // photon energy rest frame
+     Double_t r1 = fEnergy/e;
+     Double_t r2 = E0/er;
+     Double_t ctheta = (r1-r2-1)/(beta*r1-1);
+     Double_t tg = sqrt(1-ctheta*ctheta)/gamma/(beta+ctheta);
+    */
+    fEnergy -= e;
+
+    MPhoton &p = *new MPhoton(e, fZ);
+    p = *this;
+
+    /*
+     cout << "(" << atan(tg)*180/TMath::Pi() << ")" << flush;
+
+     static TRandom rand(0);
+     p->SetNewDirection(atan(tg), rand.Uniform(TMath::Pi()*2));
+     */
+
+    // MISSING: Electron angle
+    //
+    // E1 = fEnergy  (after!)
+    //
+    // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0)
+
+    return &p;
+}

trunk/WuerzburgSoft/Thomas/mphys/MElectron.h

@@ -5 +5 @@
 #include "MParticle.h"
 #endif
+
+class MPhoton;
 
 class MElectron : public MParticle
@@ -14 +16 @@
         fZ      = z;
     }
+
+    void operator=(MParticle &p) { MParticle::operator=(p); }
 
     // ----------------------------------------------------------------
@@ -31 +35 @@
     static Double_t p_e(Double_t *x, Double_t *k);
     static Double_t G_q(Double_t *x, Double_t *k);
-    static Double_t EnergyLoss(Double_t *E, Double_t *z=NULL);
-    static Double_t GetEnergyLoss(Double_t E, Double_t z=0);
+    static Double_t EnergyLoss(Double_t *E, Double_t *z=NULL, Double_t *ep=NULL);
+    static Double_t GetEnergyLoss(Double_t E, Double_t z=0, Double_t *ep=NULL);
 
-    Double_t GetEnergyLoss() const;
+    Double_t GetEnergyLoss(Double_t *ep) const;
+
+    MPhoton *DoInvCompton();
 
     ClassDef(MElectron, 1)

trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc

@@ -34 +34 @@
 
 #include "MParList.h"
-#include "MParticle.h"
+#include "MPhoton.h"
 
 ClassImp(MGenIRPhoton);
@@ -104 +104 @@
 {
     const Double_t pi2 = TMath::Pi() * 2;
-    MParticle *p = new MParticle(MParticle::kEGamma);
 
-    p->SetAngle(fRand.Uniform(pi2));
-    p->SetEnergy(fSrc->GetRandom()*1e-9);
+    MPhoton *p = new MPhoton;
+    /*
+     MParticle *p = new MParticle(MParticle::kEGamma);
+     p->SetAngle(fRand.Uniform(pi2));
+     p->SetEnergy(fSrc->GetRandom()*1e-9);
+     */
 
     return p;

trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc

@@ -71 +71 @@
 Bool_t MGenPrimaryParticle::Process()
 {
-    MParticle *p = new MParticle(MParticle::kEGamma);
-    p->SetEnergy(fSrc->GetRandom());
+//    MParticle *p = new MParticle(MParticle::kEGamma);
+//    p->SetEnergy(fSrc->GetRandom());
     return kTRUE;
 } 

trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

@@ -32 +32 @@
 #include <TF1.h>
 #include <TMath.h>
+#include <TRandom.h>
 
 #include "TList.h"
@@ -70 +71 @@
 
 // --------------------------------------------------------------------------
-Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, TList *list)
+Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, const Double_t theta, TList *list)
 {
     //
@@ -78 +79 @@
     const Double_t E0     = 511e-6;              // [GeV]
 
-    const Double_t theta  = phot->GetAngle();    // [2pi]
+    //const Double_t theta  = phot->GetAngle();    // [2pi]
     const Double_t Ep     = phot->GetEnergy();   // [GeV]
     const Double_t Eg     = gamma->GetEnergy();  // [GeV]
@@ -120 +121 @@
     const Double_t dE = E*Bcosd;
 
-    MElectron *p[2];
-    p[0] = new MElectron(E+dE, gamma->GetZ());
-    p[1] = new MElectron(E-dE, gamma->GetZ());
+    MElectron &p0 = *new MElectron(E+dE, gamma->GetZ());
+    MElectron &p1 = *new MElectron(E-dE, gamma->GetZ());
+    p0 = *gamma; // Set Position and direction
+    p1 = *gamma; // Set Position and direction
 
     const Double_t E1 = E0/(E+dE);
@@ -130 +132 @@
     const Double_t beta2 = sqrt(1.-E2*E2);
 
-    p[0]->SetBeta(beta1);
-    p[1]->SetBeta(beta2);
+    p0.SetBeta(beta1);
+    p1.SetBeta(beta2);
 
     const Double_t Bscp = Bsind*cphi;
@@ -140 +142 @@
     const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));
 
-    p[0]->SetAngle(atan(tan1));
-    p[1]->SetAngle(atan(tan2));
+    static TRandom rand(0);
+    Double_t rnd = rand.Uniform(TMath::Pi() * 2);
+    p0.SetNewDirection(atan(tan1), rnd);
+    p1.SetNewDirection(atan(tan2), rnd+TMath::Pi());
 
-    list->Add(p[0]);
-    list->Add(p[1]);
+    list->Add(&p0);
+    list->Add(&p1);
 
     return kTRUE;

trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h

@@ -19 +19 @@
     virtual ~MPairProduction();
 
-    Bool_t Process(MParticle *g, MParticle *p, TList *l);
+    Bool_t Process(MParticle *g, MParticle *p, const Double_t theta, TList *l);
 
     ClassDef(MPairProduction, 0) //

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -7 +7 @@
 #include "MParticle.h"
 
+#include <TRandom.h>
+#include <TMatrixD.h>
+
 ClassImp(MParticle);
 
+/***********************
+ *
+ * calc r from z:
+ *
+ *        R = 2*c/H0 *(1+z - sqrt(1+z))
+ *
+ *        z = -0.5 * (3 + R' +/- sqrt(R'+1))   R' = R*H0/c
+ *
+ *        H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
+ *
+ ************************
+ */
+
+Double_t ZofR(Double_t *x, Double_t *k=NULL)
+{
+    const Double_t c  = 299792458;        // [m/s]
+    const Double_t H0 = 50./3.0857e19;    // [km / Mpc s] --> [s^-1]
+
+    const Double_t ly = 3600.*24.*365.*c; // [m/ly]
+    const Double_t pc = 1./3.258;         // [pc/ly]
+    const Double_t r  = x[0] /pc*ly*1e3;  // [m]
+
+    const Double_t R = r*H0/c;            // [1]
+
+    return (R+1+sqrt(R*2+1))/2 - 1;
+}
+
+Double_t RofZ(Double_t *x, Double_t *k=NULL)
+{
+    Double_t z1 = x[0] + 1;
+
+    const Double_t c  = 299792458;                 // [m/s]
+    const Double_t H0 = 50./3.0857e19;             // [km / Mpc s] --> [s^-1]
+
+    const Double_t ly = 3600.*24.*365.*c;          // [m/ly]
+    const Double_t pc = 1./3.258;                  // [pc/ly]
+
+    const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]
+
+    return  R * pc/ly/1e3;                   // [kpc]
+}
+
 MParticle::MParticle(ParticleType_t t, const char *name, const char *title)
-    : fPType(t)
+    : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0)
 {
     //
@@ -25 +70 @@
 }
 
+void MParticle::SetNewDirection(Double_t theta, Double_t phi)
+{
+    TMatrixD B(3, 3);
+
+    B(0, 0) =  cos(fTheta)*cos(fPsi);
+    B(1, 0) =  cos(fTheta)*sin(fPsi);
+    B(2, 0) = -sin(fTheta);
+
+    B(0, 1) = -sin(fPsi);
+    B(1, 1) =  cos(fPsi);
+    B(2, 1) =  0;
+
+    B(0, 2) = sin(fTheta)*cos(fPsi);
+    B(1, 2) = sin(fTheta)*sin(fPsi);
+    B(2, 2) = cos(fTheta);
+
+    // ------------------------------
+
+    TVectorD r(3);
+
+    r(0) = sin(theta)*cos(phi);
+    r(1) = sin(theta)*sin(phi);
+    r(2) = cos(theta);
+
+    // ------------------------------
+
+    r *= B;
+
+    fTheta = sqrt(r(0)*r(0)+r(1)*r(1)); // Numerically bad: acos(r(2));
+    fPsi   = atan2(r(1), r(0));
+
+    if (fTheta*2 > TMath::Pi())
+        fTheta = fabs(fTheta-TMath::Pi());
+}
+
+#include <iostream.h>
+#include <iomanip.h>
+
+Bool_t MParticle::SetNewPosition(Double_t dr)
+{
+    Bool_t rc=kTRUE;
+
+    TVectorD x(3);
+
+    x(0) = sin(fTheta)*cos(fPsi);
+    x(1) = sin(fTheta)*sin(fPsi);
+    x(2) = cos(fTheta);
+
+    x *= dr;
+
+    // ------------------------------
+
+    const Double_t R = RofZ(&fZ);
+
+    if (x(2) > R*cos(fTheta))
+    {
+        x *= R/dr;
+        rc = kFALSE;
+    }
+
+    // ------------------------------
+
+    TVectorD r(3);
+
+    r(0) = fR*cos(fPhi);
+    r(1) = fR*sin(fPhi);
+    r(2) = R;
+
+    // ------------------------------
+
+    r -= x;
+
+    fR   = sqrt(r(0)*r(0)+r(1)*r(1));
+    fPhi = atan2(r(1), r(0));
+    fZ   = ZofR(&r(2));
+
+    return rc;
+}
+
+Bool_t MParticle::SetNewPosition()
+{
+    static TRandom rand(0);
+    Double_t r = rand.Exp(GetInteractionLength());
+    return SetNewPosition(r);
+}

trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

@@ -17 +17 @@
     Double_t fEnergy; // [GeV] Energy
     Double_t fBeta;   // [1]   beta
-    Double_t fAngle;  // [rad] Angle
+
     Double_t fZ;      // [1]   Red shift
+    Double_t fR;      // [kpc] Radius from line of view
+    Double_t fPhi;    // [rad] Phi@z from line of view
+
+    Double_t fTheta;  // [rad] Direction of momentum, angle || line of view
+    Double_t fPsi;    // [rad] Direction of momentum, az. angle
 
 public:
@@ -24 +29 @@
      //    ~MParticle();
 
+    void operator=(MParticle &p)
+    {
+        fZ = p.fZ;
+        fR = p.fR;
+        fPhi = p.fPhi;
+        fTheta = p.fTheta;
+        fPsi = p.fPsi;
+    }
+
     //    void Fill(const MParContainer *inp);
 
     //    void Draw(Option_t *opt=NULL);
+    virtual Double_t GetInteractionLength() const = 0;
 
     void SetEnergy(Double_t e) { fEnergy = e; }
-    void SetAngle(Double_t a)  { fAngle = a; }
     void SetBeta(Double_t b)   { fBeta = b; }
+
     void SetZ(Double_t z)      { fZ = z; }
+    /*
+     void SetPhi(Double_t p)    { fPhi = p; }
+     void SetR(Double_t r)      { fR = r; }
+     void SetTheta(Double_t t)  { fTheta = t; }
+     void SetPsi(Double_t p)    { fPsi = p; }
+     */
+
+    void SetNewDirection(Double_t theta, Double_t phi);
+    Bool_t SetNewPosition(Double_t dr);
+    Bool_t SetNewPosition();
 
     Double_t GetEnergy() const { return fEnergy; }
-    Double_t GetAngle() const  { return fAngle; }
+
     Double_t GetZ() const      { return fZ; }
+    Double_t GetPhi() const    { return fPhi; }
+    Double_t GetR() const      { return fR; }
+
+    Double_t GetTheta() const  { return fTheta; }
+    Double_t GetPsi() const    { return fPsi; }
 
     ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h

@@ -15 +15 @@
     }
 
+    void operator=(MParticle &p) { MParticle::operator=(p); }
+
     static Double_t Planck(Double_t *x, Double_t *k=NULL);
     static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL);

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();