Changeset 1356 for trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (modified) (3 diffs)

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);
+}