Changeset 1356 for trunk/WuerzburgSoft/Thomas/mphys

Timestamp:

06/12/02 09:35:01 (23 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (1 diff)
• MElectron.cc (modified) (3 diffs)
• MElectron.h (modified) (1 diff)
• MParticle.cc (modified) (2 diffs)
• MParticle.h (modified) (3 diffs)
• MPhoton.cc (modified) (7 diffs)
• MPhoton.h (modified) (2 diffs)
• energyloss.C (added)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+ 2002/06/12: Thomas Bretz
+
+   * MElectron.[h,cc]:
+     - added a primitive theta dependancy to DoInvCompton
+     - added DrawInteractionLength
+     
+   * MParticle.[h,cc]:
+     - added RofZ and ZofR as static functions
+
+   * MPhoton.[h,cc]:
+     - added DrawInteractionLength
+     - fixed the bug causing the 'strange factor': sigma_gg needs
+       the sqrt of s.
+
+   * energyloss.C:
+     - added
+
+
+
  2002/06/10: Thomas Bretz
 

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -33 +33 @@
 
 #include <TF1.h>
+#include <TH1.h>
+#include <TPad.h>
+#include <TCanvas.h>
+#include <TRandom.h>
 
 #include "MPhoton.h"
@@ -286 +290 @@
 #include <iostream.h>
 
-MPhoton *MElectron::DoInvCompton()
-{
+MPhoton *MElectron::DoInvCompton(Double_t theta)
+{
+    static TRandom rand(0);
+
     Double_t E0 = 511e-6; //[GeV]
 
     Double_t epsilon;
     Double_t e = GetEnergyLoss(&epsilon);
+
+    // er: photon energy before interaction, rest frame
+    // e:  photon energy after interaction, lab
+
+    Double_t gamma     = fEnergy/E0;
+    Double_t beta      = sqrt(1.-1./(gamma*gamma));
+    //Double_t gammabeta = sqrt(gamma*gamma-1);
+
+    Double_t f = fEnergy/e;
+
+    Double_t t;
+    Double_t arg;
+    do
+    {
+        t = rand.Uniform(TMath::Pi()*2);
+        Double_t er  = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame
+        arg = (f - E0/er - 1)/(f*beta+1);
+        cout << "~" << flush;
+
+    } while (arg<-1 || arg>1);
+
+    Double_t theta1s = acos(arg);
+    Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));
+
+    Double_t thetastar = thetas-theta1s;
+
+    Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
+
     /*
-     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);
-    */
+     cout << "(" << theta1/TMath::Pi()*180 << ",";
+     cout << theta1s/TMath::Pi()*180<< ",";
+     cout << arg << ")" << flush;
+     */
+
     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));
-     */
+    p.SetNewDirection(theta1, rand.Uniform(TMath::Pi()*2));
 
     // MISSING: Electron angle
@@ -327 +346 @@
     return &p;
 }
+
+void MElectron::DrawInteractionLength(Double_t z)
+{
+    if (!gPad)
+        new TCanvas("IL_Electron", "Mean Interaction Length Electron");
+    else
+        gPad->GetVirtCanvas()->cd(3);
+
+    TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0);
+    f2.SetParameter(0, z);
+
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+    f2.SetLineWidth(1);
+
+    TH1 &h=*f2.DrawCopy()->GetHistogram();
+
+    h.SetTitle("Mean Interaction Length (Electron)");
+    h.SetXTitle("E [GeV]");
+    h.SetYTitle("x [kpc]");
+
+    gPad->Modified();
+    gPad->Update();
+}
+
+void MElectron::DrawInteractionLength() const
+{
+    DrawInteractionLength(fZ);
+}

trunk/WuerzburgSoft/Thomas/mphys/MElectron.h

@@ -40 +40 @@
     Double_t GetEnergyLoss(Double_t *ep) const;
 
-    MPhoton *DoInvCompton();
+    MPhoton *DoInvCompton(Double_t theta);
+
+    // ----------------------------------------------------------------
+
+    static void DrawInteractionLength(Double_t z);
+    void DrawInteractionLength() const;
 
     ClassDef(MElectron, 1)

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -25 +25 @@
  */
 
-Double_t ZofR(Double_t *x, Double_t *k=NULL)
+Double_t MParticle::ZofR(Double_t *x, Double_t *k)
 {
     const Double_t c  = 299792458;        // [m/s]
@@ -39 +39 @@
 }
 
-Double_t RofZ(Double_t *x, Double_t *k=NULL)
+Double_t MParticle::RofZ(Double_t *x, Double_t *k)
 {
     Double_t z1 = x[0] + 1;

trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

@@ -16 +16 @@
 protected:
     Double_t fEnergy; // [GeV] Energy
-    Double_t fBeta;   // [1]   beta
+//    Double_t fBeta;   // [1]   beta
 
     Double_t fZ;      // [1]   Red shift
@@ -38 +38 @@
     }
 
+    static Double_t MParticle::ZofR(Double_t *x, Double_t *k=NULL);
+    static Double_t MParticle::RofZ(Double_t *x, Double_t *k=NULL);
+
     //    void Fill(const MParContainer *inp);
 
@@ -44 +47 @@
 
     void SetEnergy(Double_t e) { fEnergy = e; }
-    void SetBeta(Double_t b)   { fBeta = b; }
+  //  void SetBeta(Double_t b)   { fBeta = b; }
 
     void SetZ(Double_t z)      { fZ = z; }

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc

@@ -33 +33 @@
 
 #include <TF1.h>
+#include <TH1.h>
+#include <TPad.h>
+#include <TCanvas.h>
 
 ClassImp(MPhoton);
@@ -43 +46 @@
     // constants moved out of function, see below
     //
-    Double_t E   = x[0];                    // [GeV]
-    Double_t z   = k ? k[0] : 0;
-
-    Double_t T   = 2.96*(z+1);              // [K]
-    Double_t e   = 1.602176462e-19;         // [C]
-    Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
-
-    Double_t EkT = E/kB/T;
+    const Double_t E   = x[0];                    // [GeV]
+    const Double_t z   = k ? k[0] : 0;
+
+    const Double_t T   = 2.96*(z+1);              // [K]
+    const Double_t e   = 1.602176462e-19;         // [C]
+    const Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
+
+    const Double_t EkT = E/kB/T;
 
     /*
@@ -66 +69 @@
 Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL)
 {
-    Double_t s = x[0]; // omega: CM mass
-
-    Double_t E0 = 511e-6;         // [GeV]
-    Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2
-
-    Double_t m  = E0/s;
-
-    Double_t m2 = m*m;
-    Double_t beta  = sqrt(1.-m2);
-    Double_t beta2 = 1.-m2;
-
-    Double_t p1 = r0*r0*TMath::Pi()/2;
-
-    // ----- Extreme Relativistic -----
+    const Double_t m2 = x[0]; // m2: (E0/sqrt(s))^2
+
+    const Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2
+
+    const Double_t beta2 = 1.-m2;
+    const Double_t beta  = sqrt(beta2);
+
+    const Double_t p1 = r0*r0*TMath::Pi()/2;
+
+    // ----- Extreme Relativistic -------
     // return p1*2 * m*m*m* (log(2./m)-1);
-    // --------------------------------
-
-    Double_t p2 = m2;
-    Double_t p3 = 3.-beta2*beta2;
-    Double_t p4 = log((1.+beta)/(1.-beta));
-    Double_t p5 = beta*2*(1.+m2);
-
-    Double_t sigma = p1*p2*(p3*p4-p5); // [m^2]
+    // ----------------------------------
+
+    const Double_t p2 = 3.-beta2*beta2;
+    const Double_t p3 = log((1.+beta)/(1.-beta));
+    const Double_t p4 = beta*2*(1.+m2);
+
+    const Double_t sigma = p1*m2*(p2*p3-p4); // [m^2]
 
     return sigma;
@@ -95 +93 @@
 Double_t MPhoton::Int1(Double_t *x, Double_t *k=NULL)
 {
-    Double_t costheta = x[0];
-
-    Double_t Eg = k[0];
-    Double_t Ep = k[1];
-
-    Double_t E0 = 511e-6; // [GeV]
-
-    Double_t s = Eg*Ep/E0*(1.-costheta)*2;
-
-    if (s<E0)   // Why is this necessary???
+    const Double_t costheta = x[0];
+
+    const Double_t Eg = k[0];
+    const Double_t Ep = k[1];
+
+    const Double_t E0 = 511e-6; // [GeV]
+
+    Double_t s = E0/Eg*E0/Ep/(1.-costheta)/2;
+    if (s>1)   // Why is this necessary???
         return 0;
 
-    Double_t sigma = Sigma_gg(&s);  // [m^2]
+    const Double_t sigma = Sigma_gg(&s);  // [m^2]
 
     return sigma/2 * (1.-costheta); // [m^2]
@@ -114 +111 @@
 Double_t MPhoton::Int2(Double_t *x, Double_t *k)
 {
-    Double_t E0 = 511e-6; // [GeV]
+    const Double_t E0 = 511e-6; // [GeV]
 
     Double_t Ep = x[0];
-    Double_t Eg = k[0];
-
     Double_t z  = k[1];
 
+    const Double_t Eg = k[0];
+
     Double_t val[2] = { Eg, Ep };
 
-    Double_t from = -1.0;
-    Double_t to   = 1.-E0*E0/2./Eg/Ep; // Originally Was: 1.
+    const Double_t from = -1.0;
+    const Double_t to   = 1.-E0*E0/(2.*Eg*Ep); // Originally Was: 1.
 
     TF1 f("int1", Int1, from, to, 2);
-    Double_t int1   = f.Integral(from, to, val);  // [m^2]
-    Double_t planck = Planck(&Ep, &z);            // [GeV^2]
-
-    Double_t res = planck * int1;
-
-    res *= Eg/E0*1e-9; // FIXME!!!!!!!!!! WHICH FACTOR IS THIS????
+
+    const Double_t int1   = f.Integral(from, to, val);  // [m^2]
+    const Double_t planck = Planck(&Ep, &z);            // [GeV^2]
+
+    const Double_t res = planck * int1;
 
     return res; // [GeV^2 m^2]
@@ -156 +152 @@
     Double_t val[2] = { Eg, z };
 
-    Double_t lolim = E0*E0/Eg;
-    Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
+    Double_t lolim = E0*E0 > 1e-8 ? E0*E0/Eg : 1e-8/Eg;
+    Double_t inf   = 3e-11; 
 
     TF1 f("int2", Int2, lolim, inf, 2);
@@ -192 +188 @@
 }
 
+void MPhoton::DrawInteractionLength(Double_t z)
+{
+    if (!gPad)
+        new TCanvas("ILPhoton", "Mean Interaction Length Photon");
+    else
+        gPad->GetVirtCanvas()->cd(4);
+
+    TF1 f1("length", InteractionLength, 1e4, 1e11, 1);
+    f1.SetParameter(0, z);
+
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+    f1.SetMaximum(1e5);
+    f1.SetLineWidth(1);
+
+    TH1 &h=*f1.DrawCopy()->GetHistogram();
+
+    h.SetTitle("Mean Interaction Length (Photon)");
+    h.SetXTitle("E [GeV]");
+    h.SetYTitle("x [kpc]");
+
+    gPad->Modified();
+    gPad->Update();
+}
+
+void MPhoton::DrawInteractionLength() const
+{
+    DrawInteractionLength(fZ);
+}

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h

@@ -17 +17 @@
     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);
@@ -26 +28 @@
     Double_t GetInteractionLength() const;
 
+    // ----------------------------------------------------------------
+
+    static void DrawInteractionLength(Double_t z);
+    void DrawInteractionLength() const;
+
     ClassDef(MPhoton, 1)
 };