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

Timestamp:

06/12/02 14:34:27 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

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

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -1 @@
-//#ifndef __CINT__
 #include <iostream.h>
+
 #include <fstream.h>
 
@@ -8 +8 @@
 #include <TH2.h>
 #include <TList.h>
+#include <TStyle.h>
 #include <TRandom.h>
-#include <TGraph.h>
 #include <TCanvas.h>
-#include <TStyle.h>
-//#endif
+
 #include "mbase/MParContainer.h"
 #include "mphys/MPhoton.h"
@@ -19 +18 @@
 #include "mhist/MBinning.h"
 #include "mhist/MH.h"
+
 // 2.96
 // 2.87
 // 2.73
-
-/***********************
- *
- * 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)
-{
-    Double_t c  = 299792458;        // [m/s]
-    Double_t H0 = 50./3.0857e19;    // [km / Mpc s] --> [s^-1]
-
-    Double_t ly = 3600.*24.*365.*c; // [m/ly]
-    Double_t pc = 1./3.258;         // [pc/ly]
-    Double_t r  = x[0] /pc*ly*1e3;  // [m]
-
-    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;
-
-    Double_t c  = 299792458;                 // [m/s]
-    Double_t H0 = 50./3.0857e19;             // [km / Mpc s] --> [s^-1]
-
-    Double_t ly = 3600.*24.*365.*c;          // [m/ly]
-    Double_t pc = 1./3.258;                  // [pc/ly]
-
-    Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]
-
-    return  R * pc/ly/1e3;                   // [kpc]
-}
-
-Double_t AngleDistrib(Double_t *x, Double_t *k)
-{
-
-    Double_t c = x[0]; // cos(alpha)
-    Double_t b = k[0]; // sqrt(1-1/s)
-
-    Double_t b2 = b*b;
-    Double_t b4 = b2*b2;
-
-    Double_t c2 = c*c;
-    Double_t c4 = c2*c2;
-
-    Double_t u = 1. - b4*c4 +2.*b2*(1.-b2)*(1.-c2);
-    Double_t d = 1.-b2*c2;
-
-    return u/(d*d);
-}
-
-
-// ================================================================
-Double_t DiSum(Double_t *x, Double_t *k=NULL)
-{
-    Double_t t = x[0];
-
-    Double_t disum = t;
-    Double_t add = 0;
-                   
-    Double_t eps = fabs(t*1e-1);
-
-    Int_t    n   = 2;
-    Double_t pow = t*t;   // t^2
-
-    do
-    {
-        add = pow/n/n;
-
-        pow *= t;       // pow = t^n
-        n++;
-
-        disum += add;
-
-    } while (fabs(add)>eps);
-
-    return disum;
-}
-
-Double_t F(Double_t *x, Double_t *k=NULL)
-{
-    Double_t o = x[0];
-    Double_t s = -2.*o;
-
-//    if (o<1e-10)
-//        return 2.125; //-3./8.;
-
-    return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8.
-}
-
-Double_t FlimLi(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam)  mit  b-->1   (Maple)
-{
-    Double_t w = x[0];
-
-    Double_t s   = -w*2*(1+1); // -2*omega*(1+beta)
-    Double_t li2 = MElectron::Li2(&s);
-
-    return fabs(li2);
-}
-
-#include <math.h>
-void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[],
-          Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[]))
-{
-    /*
-     * ---------------------------------------------------------
-     * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
-     * ---------------------------------------------------------
-     */
-    const Double_t a2  = 0.2;
-    const Double_t a3  = 0.3;
-    const Double_t a4  = 0.6;
-    const Double_t a5  = 1.0;
-    const Double_t a6  = 0.875;
-    const Double_t b21 = 0.2;
-    const Double_t b31 = 3.0/40.0;
-    const Double_t b32 = 9.0/40.0;
-    const Double_t b41 = 0.3;
-    const Double_t b42 = -0.9;
-    const Double_t b43 = 1.2;
-    const Double_t b51 = -11.0/54.0;
-    const Double_t b52 = 2.5;
-    const Double_t b53 = -70.0/27.0;
-    const Double_t b54 = 35.0/27.0;
-    const Double_t b61 = 1631.0/55296.0;
-    const Double_t b62 = 175.0/512.0;
-    const Double_t b63 = 575.0/13824.0;
-    const Double_t b64 = 44275.0/110592.0;
-    const Double_t b65 = 253.0/4096.0;
-    const Double_t c1  = 37.0/378.0;
-    const Double_t c3  = 250.0/621.0;
-    const Double_t c4  = 125.0/594.0;
-    const Double_t c6  = 512.0/1771.0;
-    const Double_t dc5 = -277.00/14336.0;
-
-    const Double_t dc1 = c1-2825.0/27648.0;
-    const Double_t dc3 = c3-18575.0/48384.0;
-    const Double_t dc4 = c4-13525.0/55296.0;
-    const Double_t dc6 = c6-0.25;
-
-    Double_t ak2[n];
-    Double_t ak3[n];
-    Double_t ak4[n];
-    Double_t ak5[n];
-    Double_t ak6[n];
-    Double_t ytemp[n];
-
-    for (int i=0; i<n; i++)
-        ytemp[i] = y[i]+b21*h*dydx[i];
-
-    (*derivs)(x+a2*h,ytemp,ak2);
-
-    for (int i=0; i<n; i++)
-        ytemp[i] = y[i]+h*(b31*dydx[i]+b32*ak2[i]);
-
-    (*derivs)(x+a3*h,ytemp,ak3);
-
-    for (int i=0; i<n; i++)
-        ytemp[i] = y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
-
-    (*derivs)(x+a4*h,ytemp,ak4);
-
-    for (int i=0; i<n; i++)
-        ytemp[i] = y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
-
-    (*derivs)(x+a5*h,ytemp,ak5);
-
-    for (int i=0; i<n; i++)
-        ytemp[i] = y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
-
-    (*derivs)(x+a6*h,ytemp,ak6);
-
-    for (int i=0; i<n; i++)
-        yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
-
-    for (int i=0; i<n; i++)
-        yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
-}
-
-Double_t FMIN(Double_t a, Double_t b)
-{
-    return a<b ? a : b;
-}
-
-Double_t FMAX(Double_t a, Double_t b)
-{
-    return a>b ? a : b;
-}
-
-void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps,
-            Double_t yscal[], Double_t *hdid, Double_t *hnext,
-            void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs
-{
-    /*
-     * ---------------------------------------------------------
-     * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
-     * ---------------------------------------------------------
-     */
-    const Double_t SAFETY = 0.9;
-    const Double_t PGROW  = -0.2;
-    const Double_t PSHRNK = -0.25;
-    const Double_t ERRCON = 1.89e-4;
-
-    Double_t yerr[n];
-    Double_t ytemp[n];
-
-    Double_t h = htry;
-    Double_t errmax;
-    while (1)
-    {
-        rkck(y, dydx, n, *x, h, ytemp, yerr, derivs);
-
-        errmax=0.0;
-
-        for (int i=0; i<n; i++)
-            errmax = FMAX(errmax, fabs(yerr[i]/yscal[i]) );
-
-        errmax /= eps;
-
-        if (errmax <= 1.0)
-            break;
-
-        Double_t htemp = SAFETY*h*pow(errmax,PSHRNK);
-
-        h = (h >= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h));
-
-        Double_t xnew= (*x) + h;
-
-        if (xnew != *x)
-            continue;
-
-        cout << "stepsize underflow in rkqs" << endl;
-        break;
-    }
-
-    *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h;
-
-    *x += (*hdid=h);
-
-    for (int i=0; i<n; i++)
-        y[i] = ytemp[i];
-}
-
-Double_t dEdt(Double_t E, Double_t t, Double_t z0=-1)
-{
-    /* ------------------------------------
-     * Lower limit as differential Equation
-     * ------------------------------------
-
-     Double_t T     = 2.96;                  // [K]
-     Double_t sigma = 6.653e-29;             // [m^2]
-     Double_t E0    = 511e-6;                // [GeV]
-     Double_t e     = 1.602176462e-19;       // [C]
-     Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
-     Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
-     Double_t hc    = h*c;                   // [GeVm]
-     Double_t pi    = TMath::Pi();           // [1]
-     Double_t khc   = pi*kB/hc;              // [1 / K m]
-     Double_t a     = 8./15 * pi * khc*khc*khc*khc * hc;         // [Gev / K^4 / m^3]
-
-     Double_t konst = 4./3. * sigma * a *T*T*T*T *c /E0 /E0;
-
-     return -konst *E*E;
-     */
-    Double_t pc = 1./3.258;                // [pc/ly]
-    Double_t c  = 299792458;               // [m/s]
-    Double_t ly = 3600.*24.*365.*c;        // [m/ly]
-    Double_t r  = t * c / ly * pc / 1000;  // [kpc]
-    Double_t R  = RofZ(&z0) - r;
-    Double_t z  = z0>=0 ? ZofR(&R) : 0;
-
-    Double_t e     = 1.602176462e-19;       // [C]
-    Double_t T     = 2.96*(z+1);            // [K]
-    Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
-    Double_t kBT   = kB*T;                  // [GeV]
-    Double_t alpha = 1./137;                // [|]
-    Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
-    Double_t E0    = 511e-6;                // [GeV]
-
-
-    Double_t k = TMath::Pi() * alpha * kBT;
-
-    Double_t ln = 4.*kBT/E0/E0;
-
-    return -1./3/h* k*k * (log(ln*E)-1.9805);
-}
-
-void dEdt(Double_t t, Double_t E[], Double_t dedt[])
-{
-    dedt[0] = dEdt(E[0], t);
-    //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl;
-}
-
-Int_t GetSeed()
-{
-    return (Int_t)fmod(TStopwatch::GetRealTime()*10, 65535);
-}
-
-void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="")
-{
-    Double_t t = 0;
-    Double_t E[1] = { E0 };
-    Double_t yscal[1] = { 1 };
-
-    Double_t dedt[1];
-
-    Double_t eps;// = 1e5;
-    Double_t step = 5e6;
-
-    Double_t hdid;
-    Double_t hnext;
-
-    cout << "Start: " << dedt[0] << endl;
-
-    Int_t n = 15000;
-    Double_t tres[n];
-    Double_t Eres[n];
-    int i;
-    for (i=0; i<n; i++)
-    {
-        tres[i] = t;
-
-        eps = E[0]/1e9; //9;
-
-        dedt[0] = dEdt(E[0], t);
-
-        SolvEq(E, dedt, 1, &t, step, eps, yscal, &hdid, &hnext, dEdt);
-
-        step = hnext;
-
-        Eres[i] = E[0];
-
-        if (i==0) cout << "Did: " << hdid << endl;
-
-        if (t>1.5e14 || E[0]<5e6)
-            break;
-//        cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl;
-    }
-
-    cout << i << endl;
-
-    TGraph grp(i<n?i:n, tres, Eres);
-    grp.Clone()->Draw(opt);
-
-}
-
-Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL)
-{
-    Double_t t  = x[0];
-    Double_t E  = k[0];
-    Double_t t0 = k[1];
-
-    Double_t c   = 299792458;               // [m/s]
-    Double_t ly  = 3600.*24.*365.*c;        // [m/ly]
-    Double_t pc  = 1./3.258;                // [pc/ly]
-
-    Double_t r   = t * c / ly * pc / 1000;  // [kpc]
-    Double_t R   = RofZ(&k[2]) - r;
-    Double_t z   = k[2]>=0 ? ZofR(&R) : 0;
-
-    Double_t T     = 2.96*(z+1);            // [K]
-    //    Double_t alpha = 1./137;                // [1]
-    Double_t sigma = 6.653e-29;             // [m^2]
-    Double_t E0    = 511e-6;                // [GeV]
-    Double_t e     = 1.602176462e-19;       // [C]
-    Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
-    Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
-    Double_t hc    = h*c;                   // [GeVm]
-    Double_t pi    = TMath::Pi();           // [1]
-    Double_t khc   = pi*kB/hc;              // [1 / K m]
-    Double_t a     = 8./15 * pi * khc*khc*khc*khc * hc;         // [Gev / K^4 / m^3]
-    Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s]
-
-    Double_t ret = 1./(konst*(t-t0) + 1./E);
-
-    return ret;
-}
-
-void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="")
-{
-    // 8.7
-
-    Double_t val[] = { E0, t0, z };
-    Double_t t = 1.5e14;
-    while (EnergyLossRateLoLim(&t, val)<1e4)
-        t -= 0.01e14;
-
-    TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2);
-    f1->SetParameter(0, E0);
-    f1->SetParameter(1, t0);
-
-    f1->Draw(opt);
-    f1->SetBit(kCanDelete);
-}
-
-//
-// (3) Energy loss rate of electrons and 'high energy particle'
-//
-Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL)
-{
-    Double_t ly = 3600.*24.*365.; // [s/ly]
-    Double_t pc = 1./3.258;       // [pc/ly]
-
-    TGraph *graph = new TGraph;
-    graph->SetPoint(0, 0, E);
-    graph->SetMaximum(E*3); // *MENU*
-
-    cout << "------ " << endl;
-
-    static TRandom rand;
-
-    Double_t x = 0;
-    for (int i=1; i<10; i++)
-    {
-        Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z));
-
-        if (z>=0)
-        {
-            Double_t r = RofZ(&z) - l;
-
-            cout << " " << i << ". R=" << RofZ(&z) << " l=" << l << " z=" << ZofR(&r) << endl;
-
-            z = r>=0 ? ZofR(&r) : 0;
-
-            if (z==0)
-                cout << "z<0" << endl;
-        }
-
-        x += l;
-
-        Double_t t = x/pc*ly*1000;
-        graph->SetPoint(i*2-1, t, E);
-
-        Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z);
-
-        E -= e1;
-
-        if (hist)
-            hist->Fill(e1);
-
-        cout << " Ep=" << e1 << flush;
-
-        graph->SetPoint(i*2, t, E);
-    }
-
-    graph->SetMinimum(E/3); // *MENU*
-    graph->Draw(opt);
-    graph->SetBit(kCanDelete);
-
-    //if (E<31500)
-    cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl;
-
-    return E;
-}
-
-void EnergyLossRate()
-{
-    if (gPad)
-        gPad->Clear();
-
-    Double_t E = 1.5e9; //  [GeV]
-    Double_t z = 0.03;
-
-    MBinning bins;
-    bins.SetEdgesLog(18, 0.1, 1e9);
-
-    TH1D hist;
-    hist.SetName("Phot");
-    hist.SetTitle("Photons from inverse Compton");
-    MH::SetBinning(&hist, &bins);
-
-    cout << "Working..." << flush;
-
-    for (int i=0; i<50; i++)
-    {
-        cout << i << "." << flush;
-        DrawDevelopment(E, z, i?"L":"AL", &hist);
-    }
-
-    //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen
-    DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen
-    DrawDevelopmentHiLim(E, z, "L");
-    gPad->SetLogy();
-
-    new TCanvas("Photons", "Photons created in inverse Compton");
-    hist.DrawCopy();
-    gPad->SetLogx();
-
-    cout << "...done." << endl;
-}
-
-void Scatter()
-{
-    Double_t E = 1.5e9;    // [GeV]
-
-    Int_t nbins = 100;
-    Double_t lo = E/1e6;
-    Double_t up = E*10;
-
-    Double_t binsize = log10(up/lo)/nbins;
-    Double_t *scale=new Double_t[nbins+1];
-    for (int i=0; i<nbins+1; i++)
-        scale[i] = pow(10, binsize*i) * lo;
-
-    MElectron epsilon;
-    epsilon.SetEnergy(E);
-
-/*
-    TH1F *hist = new TH1F("h", "h", nbins, scale);
-
-    hist->Draw();
-
-    TStopwatch clock;
-    clock.Start();
-
-    for (int i=0; i<1000; i++)
-    {
-        Double_t e1 = epsilon.GetRandom();
-
-        hist->Fill(e1);
-    }
-
-    clock.Stop();
-    clock.Print();
-
-    epsilon.DrawCopy("same");
-    gPad->SetLogx();
-    gPad->SetLogy();
-    */
-}
-
-void DrawILPhoton(Double_t z=0)
-{
-    if (!gPad)
-        new TCanvas("IL_Photon", "Mean Interaction Length Photon");
-    else
-        gPad->GetVirtCanvas()->cd(4);
-
-    TF1 f1("length", MPhoton::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 DrawILElectron(Double_t z=0)
-{
-    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 DrawRZ()
-{
-    new TCanvas("RZ", "r and z");
-
-    TF1 f1("ZofR", ZofR, 0, 4.5e5, 0);
-    TF1 f2("RofZ", RofZ, 0, 5,     0);
-
-    gPad->Divide(2,2);
-
-    gPad->GetVirtCanvas()->cd(1);
-    TH1 *h = f1.DrawCopy()->GetHistogram();
-
-    h->SetTitle("z(r)");
-    h->SetXTitle("r [kpc]");
-    h->SetYTitle("z");
-
-    gPad->Modified();
-    gPad->Update();
-
-    gPad->GetVirtCanvas()->cd(2);
-    h = f2.DrawCopy()->GetHistogram();
-
-    h->SetTitle("r(z)");
-    h->SetXTitle("z");
-    h->SetYTitle("r [kpc]");
-
-    gPad->Modified();
-    gPad->Update();
-}
-
 Double_t PrimSpect(Double_t *x, Double_t *k)
 {
@@ -649 +27 @@
 }
 
-/*
-Double_t Planck(Double_t *x, Double_t *k=NULL)
-{
-    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;
-
-    return E*E*E / (exp(EkT)-1.); // [GeV^2]
-}
-*/
-
-Double_t Planck(Double_t *x, Double_t *k=NULL)
-{
-    //
-    // Planck, per unit volume, per unit energy
-    //
-    // constants moved out of function, see below
-    //
-
-    //
-    // This is to get a correct random value in a TF1!
-    //
-    Double_t E   = pow(10, 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;
-
-    return E*E / (exp(EkT)-1.); // [GeV^2]
+Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
+{
+    Double_t Ep  = pow(10, x[0]);
+    Double_t res = MPhoton::Int2(&Ep, k);
+
+    return res*1e55; //65/k[0];
+    // return MPhoton::Planck(&Ep, &k[1]);
+}
+
+Double_t Sbar_sigmas(Double_t *x, Double_t *k)
+{
+    Double_t sbar = pow(10, x[0]);
+
+    Double_t s = 1./(sbar*4);
+
+    Double_t sigma = MPhoton::Sigma_gg(&s);
+
+    return sigma*sbar*1e28;
+}
+
+Double_t RandomTheta(Double_t Eg, Double_t Ep)
+{
+    Double_t E0 = 511e-6; // [GeV]
+
+    Double_t f = Eg/E0*Ep/E0;
+
+    if (f<1)
+        return 0;
+
+    TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
+
+    Double_t sbar  = pow(10, func.GetRandom());
+    Double_t theta = acos(1.-sbar*2/f);
+
+    return theta;
 }
 
 void DoIt()
 {
-    Double_t rcygnusx3 = 5000; //100; //30/3.258; // [~10kpc]
-
-    Double_t startz = ZofR(&rcygnusx3);
+    Double_t R = 100; // [kpc]
+    Double_t startz = MParticle::ZofR(&R);
+
+    cout << "R = " << R << "kpc" << endl;
+    cout << "Z = " << startz << endl;
 
     static TRandom rand(0);
     MPairProduction pair;
 
-    //Double_t nphot = 50;
-    Double_t runtime = 15*60; ///*18*60*/60; // [s]
-
-    Double_t lo = 2e5;
-    Double_t hi = 2e10;
+    Double_t runtime = 30*60; // [s]
+
+    Double_t lo = 1e4; 
+    Double_t hi = 1e10;
     Double_t alpha = -2;
 
-    Double_t nbins = 4*log10(hi/lo);
-
-    TF1 phot("PhotonSpectrum", Planck, -15, -10.5, 1);
+    Double_t nbins = 24; //4*log10(hi/lo);
+
     TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
-
-    src.SetParameter(0, alpha);
+    src.SetParameter(0, 0);
+
+    TH1D hist;
+    TH1D histsrc;
+
+    hist.SetMinimum(lo*lo * pow(lo, alpha)/10);
+    histsrc.SetMinimum(lo*lo * pow(lo, alpha)/10);
 
     TList listg;
@@ -726 +107 @@
     MBinning binspolx;
     MBinning binspoly;
-    binspolx.SetEdges(32, -180, 180);
-    binspoly.SetEdgesLog(20, 1e-9, 1e-5);
+    binspolx.SetEdges(16, -180, 180);
+    binspoly.SetEdgesLog(20, 1e-13, 1e-4);
     MH::SetBinning(&position, &binspolx, &binspoly);
     MH::SetBinning(&angle,    &binspolx, &binspoly);
@@ -737 +118 @@
     histpos.SetTitle("Particle Position");
     position.SetTitle("Particle Position");
-
-    TH1D hist;
-    TH1D histsrc;
 
     //
@@ -761 +139 @@
     gPad->SetLogy();
 
+    histsrc.SetName("Spectrum");
+    histsrc.SetXTitle("E [GeV]");
+    histsrc.SetYTitle("E^{2}\\dotCounts");
+    histsrc.GetXaxis()->SetLabelOffset(-0.015);
+    histsrc.GetXaxis()->SetTitleOffset(1.1);
+
     histsrc.Draw("P");
     hist.Draw("Psame");
     histsrc.Draw("Csame");
     hist.Draw("Csame");
-
-    histsrc.SetName("Spectrum");
-    histsrc.SetXTitle("E [GeV]");
-    histsrc.SetYTitle("E*Counts");
-    histsrc.GetXaxis()->SetLabelOffset(-0.015);
-    histsrc.GetXaxis()->SetTitleOffset(1.1);
 
     ofstream fsrc("source.dat");
@@ -789 +167 @@
 
         Double_t E = pow(10, src.GetRandom());
-
-        histsrc.Fill(E, E);
+        Double_t weight = pow(E, alpha);
+
+        histsrc.Fill(E, E*E * weight);
+
         fsrc << E << endl;
 
@@ -800 +180 @@
         while (1)
         {
-            cout << " |P:" << flush;
+            if (listg.GetSize()!=0)
+                cout << " |P:" << flush;
 
             TIter NextP(&listg);
@@ -806 +187 @@
             while ((p=(MPhoton*)NextP()))
             {
+                Double_t Eg = p->GetEnergy();
+
                 if (!p->SetNewPosition())
                 {
                     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;
+                    hist.Fill(Eg, Eg*Eg*weight);
+                    position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight);
+                    angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
+                    histpos.Fill(p->GetR(), weight);
+                    hista.Fill(p->GetTheta()*kRad2Deg, weight);
+
+                    fcas << Eg << endl;
                     delete listg.Remove(p);
                     continue;
@@ -824 +207 @@
                 // Sample phtoton from background and interaction angle
                 //
-                MPhoton photon;
-
-                phot.SetParameter(0, p->GetZ());
-                photon.SetEnergy(pow(10, phot.GetRandom()));
-                Double_t theta = rand.Uniform(TMath::Pi() * 2);
-
-                if (!pair.Process(p, &photon, theta, &liste))
-                {
-                    cout << "0" << flush;;
+                TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
+                phot.SetParameter(0, Eg);
+                phot.SetParameter(1, p->GetZ());
+                if (phot.Integral(-log10(Eg)-8, -10.5)==0)
+                {
+                    // ???????????????????
+                    hist.Fill(Eg, Eg*Eg*weight);
+                    angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
+                    hista.Fill(p->GetTheta()*kRad2Deg, weight);
+                    delete listg.Remove(p);
+                    cout << "z" << flush;
+                    continue;
+                }
+
+                Double_t Ep = pow(10, phot.GetRandom());
+
+                Double_t theta = RandomTheta(Eg, Ep);
+                if (theta==0)
+                {
+                    cout << "t" << flush;
+                    continue;
+                }
+
+                if (!pair.Process(p, Ep, theta, &liste))
+                {
+                    cout << "0" << flush;
                     continue;
                 }
@@ -843 +243 @@
                 break;
 
-            cout << " |E" << flush;
+            if (liste.GetSize()!=0)
+                cout << " |E" << flush;
 
             TIter Next(&liste);
@@ -849 +250 @@
             while ((e=(MElectron*)Next()))
             {
+                if (e->GetEnergy()<lo)
+                    continue;
+
                 cout << ":" << flush;
                 while(1)
@@ -858 +262 @@
                     }
 
-                    MPhoton *p = e->DoInvCompton();
-                    if (e->GetEnergy()<lo*5)
+                    // WRONG!
+                    Double_t theta = rand.Uniform(TMath::Pi()*2);
+                    MPhoton *p = e->DoInvCompton(theta);
+
+                    if (fabs(p->GetTheta()*3437)<1 &&  // < 1min
+                        p->GetEnergy()>lo)
                     {
-                        cout << "0" << flush;
-                        delete p;
-                        break;
+                        cout << "." << flush;
+                        listg.Add(p);
                     }
-
-                    if (p->GetEnergy()<lo)
+                    else
                     {
                         cout << "i" << flush; // ignored
                         delete p;
+                    }
+
+                    if (fabs(e->GetTheta()*3437)<1 &&  // < 1min
+                        e->GetEnergy()>lo*2)
                         continue;
-                    }
-
-                    listg.Add(p);
-                    cout << "." << flush;
+
+                    cout << "0" << flush;
+                    break;
                 }
             }
@@ -940 +349 @@
     hista.DrawCopy();
     gPad->SetLogx();
-
 }
 
@@ -946 +354 @@
 void phys()
 {
+    DoIt();
+
     /*
-    Double_t Eg = 1e5;
-
-    Double_t E0    = 511e-6;                  // [GeV]
-    Double_t lolim = E0*E0/Eg;
-    Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
-
-    // 1e5:   5e-13, 4e-12
-    // 1e6:   5e-14, 2e-12
-    // 1e8:   5e-16, 1e-12
-    // 1e10:  5e-18, 1e-12
-
-    // 1e5:   2.6e-12, 4e-12
-    // 1e6:   2.6e-13, 2e-12
-    // 1e8:   2.6e-15, 1e-12
-    // 1e10:  1.6e-17, 1e-12
-
-    TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
-    f.SetParameter(0, Eg);
-    f.SetParameter(1, 0);
-
-    MH::MakeDefCanvas();
-    gPad->SetLogx();
-    f.DrawCopy();
-    return;
+     Double_t E = 1e10;
+     TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
+     phot.SetParameter(0, E);
+     phot.SetParameter(1, 5);
+     phot.DrawCopy();
+     return;
+     */
+
+    /*
+     Double_t Eg = 1e5;
+
+     Double_t E0    = 511e-6;                  // [GeV]
+     Double_t lolim = E0*E0/Eg;
+     Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
+
+     // 1e5:   5e-13, 4e-12
+     // 1e6:   5e-14, 2e-12
+     // 1e8:   5e-16, 1e-12
+     // 1e10:  5e-18, 1e-12
+
+     // 1e5:   2.6e-12, 4e-12
+     // 1e6:   2.6e-13, 2e-12
+     // 1e8:   2.6e-15, 1e-12
+     // 1e10:  1.6e-17, 1e-12
+
+     TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
+     f.SetParameter(0, Eg);
+     f.SetParameter(1, 0);
+
+     MH::MakeDefCanvas();
+     gPad->SetLogx();
+     f.DrawCopy();
     */
-//    MH::MakeDefCanvas();
-//    DrawILPhoton();
-
-//    return ;
-
-    DoIt();
-return;
-    EnergyLossRate();
-
-    MH::MakeDefCanvas();
-    DrawILElectron();
-
-    DrawRZ();
-
-    return;
-
-    /*
-    cout << "Starting..." << endl;
-
-    MEvtLoop   magic;
-    MParList   plist;
-    MTaskList  tlist;
-    plist.AddToList(&tlist);
-
-    MHInput  *input  = new MHInput;
-    MHOutput *output = new MHOutput;
-
-    plist.AddToList(input);
-    plist.AddToList(output);
-
-    MGenerateInput geninp;
-    tlist.AddToList(&geninp);
-
-    MFillH finput("MSinglePairInput", input);
-    tlist.AddToList(&finput);
-
-    MSinglePairProduction prod;
-    tlist.AddToList(&prod);
-
-    MFillH foutput("MSinglePairOutput", output);
-    tlist.AddToList(&foutput);
-
-    magic.SetParList(&plist);
-    if (!magic.Eventloop(10000))
-        return;
-
-    input->Draw();
-    output->Draw();
-
-    return;
-*/
-
-    /*
-    MEvtLoop   magic;
-    MParList   plist;
-    MTaskList  tlist;
-    plist.AddToList(&tlist);
-
-    MHInput  *input  = new MHInput;
-    MHOutput *output = new MHOutput;
-
-    plist.AddToList(input);
-    plist.AddToList(output);
-
-    MGenerateInput geninp;
-    tlist.AddToList(&geninp);
-
-    MFillH finput("MSinglePairInput", input);
-    tlist.AddToList(&finput);
-
-    MSinglePairProduction prod;
-    tlist.AddToList(&prod);
-
-    MFillH foutput("MSinglePairOutput", output);
-    tlist.AddToList(&foutput);
-
-    magic.SetParList(&plist);
-    if (!magic.Eventloop(10000))
-        return;
-
-    input->Draw();
-    output->Draw();
-
-    return;
-   */
-}
-
+}
+