Ignore:
Timestamp:
06/12/02 14:34:27 (23 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/WuerzburgSoft/Thomas/mphys/Changelog

    r1356 r1357  
    88   * MParticle.[h,cc]:
    99     - added RofZ and ZofR as static functions
     10     - changed from MParContainer to TObject
     11     - removed unnecessary or commented inline functions
    1012
    1113   * MPhoton.[h,cc]:
    … …  
    1618   * energyloss.C:
    1719     - added
     20
     21   * MPairProduction.[h,cc]:
     22     - simplified formula
     23     - changed arguments to Energy of Photon and interaction angle
     24
     25   * phys.C:
     26     - added correct photon energy for pair production
     27     - added correct angle for pair production
     28     - changed output to E^2*Counts
     29     - changed gamma production from E^-2 to Uniform by weighting
     30     - removed all unecessary functions and code (s.energyloss.C)
    1831
    1932

  • trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

    r1352 r1357  
    7171
    7272// --------------------------------------------------------------------------
    73 Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, const Double_t theta, TList *list)
     73Bool_t MPairProduction::Process(MParticle *gamma, const Double_t Ep, const Double_t theta, TList *list)
    7474{
    7575    //
    7676    //  gamma: primary particle from source
    7777    //  phot:  infrared photon from background. (angle = interaction angle)
    78     //
     78    //  Ep: Energy photon [GeV]
     79    //  theta: interaction angle [rad]
     80    //
     81
     82
    7983    const Double_t E0     = 511e-6;              // [GeV]
    80 
    81     //const Double_t theta  = phot->GetAngle();    // [2pi]
    82     const Double_t Ep     = phot->GetEnergy();   // [GeV]
    8384    const Double_t Eg     = gamma->GetEnergy();  // [GeV]
    8485
    8586    const Double_t ctheta = cos(theta);
    8687
    87     const Double_t s      = (Ep*Eg/(2*E0*E0))*(1-ctheta); //[1]
    88     if (s<1)
     88    const Double_t s      = Ep*Eg*2*(1-ctheta); //[1]
     89    const Double_t S      = s/(E0*E0*4); //[1]
     90    if (S<1)
    8991        return kFALSE;
    9092
     93    const Double_t stheta = sin(theta);
     94
     95    const Double_t sqrbetah = s/((Eg+Ep)*(Eg+Ep)) + 1;
     96    const Double_t sqrbetae = 1.-1./S;
     97
     98    const Double_t GammaH = (Eg+Ep)/sqrt(s);
     99
     100    Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));
     101
     102    fAngle->SetParameter(0, sqrt(sqrbetae));
     103
     104    const Double_t alpha = psi-acos(fAngle->GetRandom());
     105
     106    Double_t salpha = sin(alpha);
     107    Double_t calpha = cos(alpha);
     108
     109    const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi)
     110
     111    Double_t bb = sqrt(sqrbetah/sqrbetae);
     112
     113    Double_t s1 = calpha/GammaH;
     114    Double_t s2 = tphi*s1 - salpha - bb;
     115
     116    Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
     117    Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;
     118
     119    const Double_t E = (Eg+Ep)/2;;
     120    const Double_t f = sqrt(sqrbetah*sqrbetae)*salpha;
     121
     122    cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush;
     123
     124    MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ());
     125    MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ());
     126    p0 = *gamma; // Set Position and direction
     127    p1 = *gamma; // Set Position and direction
     128
     129    static TRandom rand(0);
     130    Double_t rnd = rand.Uniform(TMath::Pi() * 2);
     131    p0.SetNewDirection(atan(+tan1), rnd);
     132    p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi());
     133
     134    list->Add(&p0);
     135    list->Add(&p1);
     136
     137    /*
    91138    const Double_t Epg    = Ep/Eg; // 1+Epg ~ 1
    92139    const Double_t Egp    = Eg/Ep; // 1+Egp ~ Egp
    … …  
    120167    const Double_t E = sqrt(s)*E0/gamma1;
    121168    const Double_t dE = E*Bcosd;
     169
     170    const Double_t E1 = E0/(E+dE);
     171    const Double_t E2 = E0/(E-dE);
     172
     173    const Double_t beta1 = sqrt(1.-E1*E1);
     174    const Double_t beta2 = sqrt(1.-E2*E2);
     175
     176    const Double_t Bscp = Bsind*cphi;
     177    const Double_t spg  = sphi/gamma1;
     178    const Double_t cpg  = cphi/gamma1;
     179
     180    const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp));
     181    const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));
    122182
    123183    MElectron &p0 = *new MElectron(E+dE, gamma->GetZ());
    … …  
    126186    p1 = *gamma; // Set Position and direction
    127187
    128     const Double_t E1 = E0/(E+dE);
    129     const Double_t E2 = E0/(E-dE);
    130 
    131     const Double_t beta1 = sqrt(1.-E1*E1);
    132     const Double_t beta2 = sqrt(1.-E2*E2);
    133 
    134188    p0.SetBeta(beta1);
    135189    p1.SetBeta(beta2);
    136 
    137     const Double_t Bscp = Bsind*cphi;
    138     const Double_t spg  = sphi/gamma1;
    139     const Double_t cpg  = cphi/gamma1;
    140 
    141     const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp));
    142     const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));
    143190
    144191    static TRandom rand(0);
    … …  
    149196    list->Add(&p0);
    150197    list->Add(&p1);
    151 
     198    */
    152199    return kTRUE;
    153200}

  • trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h

    r1352 r1357  
    1919    virtual ~MPairProduction();
    2020
    21     Bool_t Process(MParticle *g, MParticle *p, const Double_t theta, TList *l);
     21    Bool_t Process(MParticle *g, const Double_t Ep, const Double_t theta, TList *l);
    2222
    2323    ClassDef(MPairProduction, 0) //

  • trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

    r1356 r1357  
    2424 ************************
    2525 */
     26#include <iostream.h>
     27#include <iomanip.h>
    2628
    2729Double_t MParticle::ZofR(Double_t *x, Double_t *k)
    … …  
    6567    //   set the name and title of this object
    6668    //
    67    
    68     fName  = name  ? name  : "MParticle";
    69     fTitle = title ? title : "Container for a particle";
     69
     70    /*
     71     fName  = name  ? name  : "MParticle";
     72     fTitle = title ? title : "Container for a particle";
     73     */
    7074}
    7175
    … …  
    104108        fTheta = fabs(fTheta-TMath::Pi());
    105109}
    106 
    107 #include <iostream.h>
    108 #include <iomanip.h>
    109110
    110111Bool_t MParticle::SetNewPosition(Double_t dr)

  • trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

    r1356 r1357  
    22#define MARS_MParticle
    33
    4 #ifndef MARS_MParContainer
    5 #include "MParContainer.h"
     4#ifndef ROOT_TObject
     5#include <TObject.h>
    66#endif
    77
    8 class MParticle : public MParContainer
     8/*
     9 #ifndef MARS_MParContainer
     10 #include "MParContainer.h"
     11 #endif
     12 */
     13
     14class MParticle : public TObject
    915{
    1016public:
    … …  
    1622protected:
    1723    Double_t fEnergy; // [GeV] Energy
    18 //    Double_t fBeta;   // [1]   beta
    1924
    2025    Double_t fZ;      // [1]   Red shift
    … …  
    2732public:
    2833    MParticle(ParticleType_t t, const char *name=NULL, const char *title=NULL);
    29      //    ~MParticle();
     34
     35    static Double_t ZofR(Double_t *x, Double_t *k=NULL);
     36    static Double_t RofZ(Double_t *x, Double_t *k=NULL);
    3037
    3138    void operator=(MParticle &p)
    … …  
    3845    }
    3946
    40     static Double_t MParticle::ZofR(Double_t *x, Double_t *k=NULL);
    41     static Double_t MParticle::RofZ(Double_t *x, Double_t *k=NULL);
    42 
    43     //    void Fill(const MParContainer *inp);
    44 
    45     //    void Draw(Option_t *opt=NULL);
    4647    virtual Double_t GetInteractionLength() const = 0;
    4748
    4849    void SetEnergy(Double_t e) { fEnergy = e; }
    49   //  void SetBeta(Double_t b)   { fBeta = b; }
    50 
    5150    void SetZ(Double_t z)      { fZ = z; }
    52     /*
    53      void SetPhi(Double_t p)    { fPhi = p; }
    54      void SetR(Double_t r)      { fR = r; }
    55      void SetTheta(Double_t t)  { fTheta = t; }
    56      void SetPsi(Double_t p)    { fPsi = p; }
    57      */
    5851
    5952    void SetNewDirection(Double_t theta, Double_t phi);

  • trunk/WuerzburgSoft/Thomas/mphys/phys.C

    r1355 r1357  
    1 //#ifndef __CINT__
    21#include <iostream.h>
     2
    33#include <fstream.h>
    44
    … …  
    88#include <TH2.h>
    99#include <TList.h>
     10#include <TStyle.h>
    1011#include <TRandom.h>
    11 #include <TGraph.h>
    1212#include <TCanvas.h>
    13 #include <TStyle.h>
    14 //#endif
     13
    1514#include "mbase/MParContainer.h"
    1615#include "mphys/MPhoton.h"
    … …  
    1918#include "mhist/MBinning.h"
    2019#include "mhist/MH.h"
     20
    2121// 2.96
    2222// 2.87
    2323// 2.73
    24 
    25 /***********************
    26  *
    27  * calc r from z:
    28  *
    29  *        R = 2*c/H0 *(1+z - sqrt(1+z))
    30  *
    31  *        z = -0.5 * (3 + R' +/- sqrt(R'+1))   R' = R*H0/c
    32  *
    33  *        H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
    34  *
    35  ************************
    36  */
    37 
    38 Double_t ZofR(Double_t *x, Double_t *k=NULL)
    39 {
    40     Double_t c  = 299792458;        // [m/s]
    41     Double_t H0 = 50./3.0857e19;    // [km / Mpc s] --> [s^-1]
    42 
    43     Double_t ly = 3600.*24.*365.*c; // [m/ly]
    44     Double_t pc = 1./3.258;         // [pc/ly]
    45     Double_t r  = x[0] /pc*ly*1e3;  // [m]
    46 
    47     Double_t R = r*H0/c;            // [1]
    48 
    49     return (R+1+sqrt(R*2+1))/2 - 1;
    50 }
    51 
    52 Double_t RofZ(Double_t *x, Double_t *k=NULL)
    53 {
    54     Double_t z1 = x[0] + 1;
    55 
    56     Double_t c  = 299792458;                 // [m/s]
    57     Double_t H0 = 50./3.0857e19;             // [km / Mpc s] --> [s^-1]
    58 
    59     Double_t ly = 3600.*24.*365.*c;          // [m/ly]
    60     Double_t pc = 1./3.258;                  // [pc/ly]
    61 
    62     Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]
    63 
    64     return  R * pc/ly/1e3;                   // [kpc]
    65 }
    66 
    67 Double_t AngleDistrib(Double_t *x, Double_t *k)
    68 {
    69 
    70     Double_t c = x[0]; // cos(alpha)
    71     Double_t b = k[0]; // sqrt(1-1/s)
    72 
    73     Double_t b2 = b*b;
    74     Double_t b4 = b2*b2;
    75 
    76     Double_t c2 = c*c;
    77     Double_t c4 = c2*c2;
    78 
    79     Double_t u = 1. - b4*c4 +2.*b2*(1.-b2)*(1.-c2);
    80     Double_t d = 1.-b2*c2;
    81 
    82     return u/(d*d);
    83 }
    84 
    85 
    86 // ================================================================
    87 Double_t DiSum(Double_t *x, Double_t *k=NULL)
    88 {
    89     Double_t t = x[0];
    90 
    91     Double_t disum = t;
    92     Double_t add = 0;
    93                    
    94     Double_t eps = fabs(t*1e-1);
    95 
    96     Int_t    n   = 2;
    97     Double_t pow = t*t;   // t^2
    98 
    99     do
    100     {
    101         add = pow/n/n;
    102 
    103         pow *= t;       // pow = t^n
    104         n++;
    105 
    106         disum += add;
    107 
    108     } while (fabs(add)>eps);
    109 
    110     return disum;
    111 }
    112 
    113 Double_t F(Double_t *x, Double_t *k=NULL)
    114 {
    115     Double_t o = x[0];
    116     Double_t s = -2.*o;
    117 
    118 //    if (o<1e-10)
    119 //        return 2.125; //-3./8.;
    120 
    121     return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8.
    122 }
    123 
    124 Double_t FlimLi(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam)  mit  b-->1   (Maple)
    125 {
    126     Double_t w = x[0];
    127 
    128     Double_t s   = -w*2*(1+1); // -2*omega*(1+beta)
    129     Double_t li2 = MElectron::Li2(&s);
    130 
    131     return fabs(li2);
    132 }
    133 
    134 #include <math.h>
    135 void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[],
    136           Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[]))
    137 {
    138     /*
    139      * ---------------------------------------------------------
    140      * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
    141      * ---------------------------------------------------------
    142      */
    143     const Double_t a2  = 0.2;
    144     const Double_t a3  = 0.3;
    145     const Double_t a4  = 0.6;
    146     const Double_t a5  = 1.0;
    147     const Double_t a6  = 0.875;
    148     const Double_t b21 = 0.2;
    149     const Double_t b31 = 3.0/40.0;
    150     const Double_t b32 = 9.0/40.0;
    151     const Double_t b41 = 0.3;
    152     const Double_t b42 = -0.9;
    153     const Double_t b43 = 1.2;
    154     const Double_t b51 = -11.0/54.0;
    155     const Double_t b52 = 2.5;
    156     const Double_t b53 = -70.0/27.0;
    157     const Double_t b54 = 35.0/27.0;
    158     const Double_t b61 = 1631.0/55296.0;
    159     const Double_t b62 = 175.0/512.0;
    160     const Double_t b63 = 575.0/13824.0;
    161     const Double_t b64 = 44275.0/110592.0;
    162     const Double_t b65 = 253.0/4096.0;
    163     const Double_t c1  = 37.0/378.0;
    164     const Double_t c3  = 250.0/621.0;
    165     const Double_t c4  = 125.0/594.0;
    166     const Double_t c6  = 512.0/1771.0;
    167     const Double_t dc5 = -277.00/14336.0;
    168 
    169     const Double_t dc1 = c1-2825.0/27648.0;
    170     const Double_t dc3 = c3-18575.0/48384.0;
    171     const Double_t dc4 = c4-13525.0/55296.0;
    172     const Double_t dc6 = c6-0.25;
    173 
    174     Double_t ak2[n];
    175     Double_t ak3[n];
    176     Double_t ak4[n];
    177     Double_t ak5[n];
    178     Double_t ak6[n];
    179     Double_t ytemp[n];
    180 
    181     for (int i=0; i<n; i++)
    182         ytemp[i] = y[i]+b21*h*dydx[i];
    183 
    184     (*derivs)(x+a2*h,ytemp,ak2);
    185 
    186     for (int i=0; i<n; i++)
    187         ytemp[i] = y[i]+h*(b31*dydx[i]+b32*ak2[i]);
    188 
    189     (*derivs)(x+a3*h,ytemp,ak3);
    190 
    191     for (int i=0; i<n; i++)
    192         ytemp[i] = y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
    193 
    194     (*derivs)(x+a4*h,ytemp,ak4);
    195 
    196     for (int i=0; i<n; i++)
    197         ytemp[i] = y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
    198 
    199     (*derivs)(x+a5*h,ytemp,ak5);
    200 
    201     for (int i=0; i<n; i++)
    202         ytemp[i] = y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
    203 
    204     (*derivs)(x+a6*h,ytemp,ak6);
    205 
    206     for (int i=0; i<n; i++)
    207         yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
    208 
    209     for (int i=0; i<n; i++)
    210         yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
    211 }
    212 
    213 Double_t FMIN(Double_t a, Double_t b)
    214 {
    215     return a<b ? a : b;
    216 }
    217 
    218 Double_t FMAX(Double_t a, Double_t b)
    219 {
    220     return a>b ? a : b;
    221 }
    222 
    223 void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps,
    224             Double_t yscal[], Double_t *hdid, Double_t *hnext,
    225             void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs
    226 {
    227     /*
    228      * ---------------------------------------------------------
    229      * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
    230      * ---------------------------------------------------------
    231      */
    232     const Double_t SAFETY = 0.9;
    233     const Double_t PGROW  = -0.2;
    234     const Double_t PSHRNK = -0.25;
    235     const Double_t ERRCON = 1.89e-4;
    236 
    237     Double_t yerr[n];
    238     Double_t ytemp[n];
    239 
    240     Double_t h = htry;
    241     Double_t errmax;
    242     while (1)
    243     {
    244         rkck(y, dydx, n, *x, h, ytemp, yerr, derivs);
    245 
    246         errmax=0.0;
    247 
    248         for (int i=0; i<n; i++)
    249             errmax = FMAX(errmax, fabs(yerr[i]/yscal[i]) );
    250 
    251         errmax /= eps;
    252 
    253         if (errmax <= 1.0)
    254             break;
    255 
    256         Double_t htemp = SAFETY*h*pow(errmax,PSHRNK);
    257 
    258         h = (h >= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h));
    259 
    260         Double_t xnew= (*x) + h;
    261 
    262         if (xnew != *x)
    263             continue;
    264 
    265         cout << "stepsize underflow in rkqs" << endl;
    266         break;
    267     }
    268 
    269     *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h;
    270 
    271     *x += (*hdid=h);
    272 
    273     for (int i=0; i<n; i++)
    274         y[i] = ytemp[i];
    275 }
    276 
    277 Double_t dEdt(Double_t E, Double_t t, Double_t z0=-1)
    278 {
    279     /* ------------------------------------
    280      * Lower limit as differential Equation
    281      * ------------------------------------
    282 
    283      Double_t T     = 2.96;                  // [K]
    284      Double_t sigma = 6.653e-29;             // [m^2]
    285      Double_t E0    = 511e-6;                // [GeV]
    286      Double_t e     = 1.602176462e-19;       // [C]
    287      Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
    288      Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
    289      Double_t hc    = h*c;                   // [GeVm]
    290      Double_t pi    = TMath::Pi();           // [1]
    291      Double_t khc   = pi*kB/hc;              // [1 / K m]
    292      Double_t a     = 8./15 * pi * khc*khc*khc*khc * hc;         // [Gev / K^4 / m^3]
    293 
    294      Double_t konst = 4./3. * sigma * a *T*T*T*T *c /E0 /E0;
    295 
    296      return -konst *E*E;
    297      */
    298     Double_t pc = 1./3.258;                // [pc/ly]
    299     Double_t c  = 299792458;               // [m/s]
    300     Double_t ly = 3600.*24.*365.*c;        // [m/ly]
    301     Double_t r  = t * c / ly * pc / 1000;  // [kpc]
    302     Double_t R  = RofZ(&z0) - r;
    303     Double_t z  = z0>=0 ? ZofR(&R) : 0;
    304 
    305     Double_t e     = 1.602176462e-19;       // [C]
    306     Double_t T     = 2.96*(z+1);            // [K]
    307     Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
    308     Double_t kBT   = kB*T;                  // [GeV]
    309     Double_t alpha = 1./137;                // [|]
    310     Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
    311     Double_t E0    = 511e-6;                // [GeV]
    312 
    313 
    314     Double_t k = TMath::Pi() * alpha * kBT;
    315 
    316     Double_t ln = 4.*kBT/E0/E0;
    317 
    318     return -1./3/h* k*k * (log(ln*E)-1.9805);
    319 }
    320 
    321 void dEdt(Double_t t, Double_t E[], Double_t dedt[])
    322 {
    323     dedt[0] = dEdt(E[0], t);
    324     //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl;
    325 }
    326 
    327 Int_t GetSeed()
    328 {
    329     return (Int_t)fmod(TStopwatch::GetRealTime()*10, 65535);
    330 }
    331 
    332 void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="")
    333 {
    334     Double_t t = 0;
    335     Double_t E[1] = { E0 };
    336     Double_t yscal[1] = { 1 };
    337 
    338     Double_t dedt[1];
    339 
    340     Double_t eps;// = 1e5;
    341     Double_t step = 5e6;
    342 
    343     Double_t hdid;
    344     Double_t hnext;
    345 
    346     cout << "Start: " << dedt[0] << endl;
    347 
    348     Int_t n = 15000;
    349     Double_t tres[n];
    350     Double_t Eres[n];
    351     int i;
    352     for (i=0; i<n; i++)
    353     {
    354         tres[i] = t;
    355 
    356         eps = E[0]/1e9; //9;
    357 
    358         dedt[0] = dEdt(E[0], t);
    359 
    360         SolvEq(E, dedt, 1, &t, step, eps, yscal, &hdid, &hnext, dEdt);
    361 
    362         step = hnext;
    363 
    364         Eres[i] = E[0];
    365 
    366         if (i==0) cout << "Did: " << hdid << endl;
    367 
    368         if (t>1.5e14 || E[0]<5e6)
    369             break;
    370 //        cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl;
    371     }
    372 
    373     cout << i << endl;
    374 
    375     TGraph grp(i<n?i:n, tres, Eres);
    376     grp.Clone()->Draw(opt);
    377 
    378 }
    379 
    380 Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL)
    381 {
    382     Double_t t  = x[0];
    383     Double_t E  = k[0];
    384     Double_t t0 = k[1];
    385 
    386     Double_t c   = 299792458;               // [m/s]
    387     Double_t ly  = 3600.*24.*365.*c;        // [m/ly]
    388     Double_t pc  = 1./3.258;                // [pc/ly]
    389 
    390     Double_t r   = t * c / ly * pc / 1000;  // [kpc]
    391     Double_t R   = RofZ(&k[2]) - r;
    392     Double_t z   = k[2]>=0 ? ZofR(&R) : 0;
    393 
    394     Double_t T     = 2.96*(z+1);            // [K]
    395     //    Double_t alpha = 1./137;                // [1]
    396     Double_t sigma = 6.653e-29;             // [m^2]
    397     Double_t E0    = 511e-6;                // [GeV]
    398     Double_t e     = 1.602176462e-19;       // [C]
    399     Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
    400     Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
    401     Double_t hc    = h*c;                   // [GeVm]
    402     Double_t pi    = TMath::Pi();           // [1]
    403     Double_t khc   = pi*kB/hc;              // [1 / K m]
    404     Double_t a     = 8./15 * pi * khc*khc*khc*khc * hc;         // [Gev / K^4 / m^3]
    405     Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s]
    406 
    407     Double_t ret = 1./(konst*(t-t0) + 1./E);
    408 
    409     return ret;
    410 }
    411 
    412 void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="")
    413 {
    414     // 8.7
    415 
    416     Double_t val[] = { E0, t0, z };
    417     Double_t t = 1.5e14;
    418     while (EnergyLossRateLoLim(&t, val)<1e4)
    419         t -= 0.01e14;
    420 
    421     TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2);
    422     f1->SetParameter(0, E0);
    423     f1->SetParameter(1, t0);
    424 
    425     f1->Draw(opt);
    426     f1->SetBit(kCanDelete);
    427 }
    428 
    429 //
    430 // (3) Energy loss rate of electrons and 'high energy particle'
    431 //
    432 Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL)
    433 {
    434     Double_t ly = 3600.*24.*365.; // [s/ly]
    435     Double_t pc = 1./3.258;       // [pc/ly]
    436 
    437     TGraph *graph = new TGraph;
    438     graph->SetPoint(0, 0, E);
    439     graph->SetMaximum(E*3); // *MENU*
    440 
    441     cout << "------ " << endl;
    442 
    443     static TRandom rand;
    444 
    445     Double_t x = 0;
    446     for (int i=1; i<10; i++)
    447     {
    448         Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z));
    449 
    450         if (z>=0)
    451         {
    452             Double_t r = RofZ(&z) - l;
    453 
    454             cout << " " << i << ". R=" << RofZ(&z) << " l=" << l << " z=" << ZofR(&r) << endl;
    455 
    456             z = r>=0 ? ZofR(&r) : 0;
    457 
    458             if (z==0)
    459                 cout << "z<0" << endl;
    460         }
    461 
    462         x += l;
    463 
    464         Double_t t = x/pc*ly*1000;
    465         graph->SetPoint(i*2-1, t, E);
    466 
    467         Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z);
    468 
    469         E -= e1;
    470 
    471         if (hist)
    472             hist->Fill(e1);
    473 
    474         cout << " Ep=" << e1 << flush;
    475 
    476         graph->SetPoint(i*2, t, E);
    477     }
    478 
    479     graph->SetMinimum(E/3); // *MENU*
    480     graph->Draw(opt);
    481     graph->SetBit(kCanDelete);
    482 
    483     //if (E<31500)
    484     cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl;
    485 
    486     return E;
    487 }
    488 
    489 void EnergyLossRate()
    490 {
    491     if (gPad)
    492         gPad->Clear();
    493 
    494     Double_t E = 1.5e9; //  [GeV]
    495     Double_t z = 0.03;
    496 
    497     MBinning bins;
    498     bins.SetEdgesLog(18, 0.1, 1e9);
    499 
    500     TH1D hist;
    501     hist.SetName("Phot");
    502     hist.SetTitle("Photons from inverse Compton");
    503     MH::SetBinning(&hist, &bins);
    504 
    505     cout << "Working..." << flush;
    506 
    507     for (int i=0; i<50; i++)
    508     {
    509         cout << i << "." << flush;
    510         DrawDevelopment(E, z, i?"L":"AL", &hist);
    511     }
    512 
    513     //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen
    514     DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen
    515     DrawDevelopmentHiLim(E, z, "L");
    516     gPad->SetLogy();
    517 
    518     new TCanvas("Photons", "Photons created in inverse Compton");
    519     hist.DrawCopy();
    520     gPad->SetLogx();
    521 
    522     cout << "...done." << endl;
    523 }
    524 
    525 void Scatter()
    526 {
    527     Double_t E = 1.5e9;    // [GeV]
    528 
    529     Int_t nbins = 100;
    530     Double_t lo = E/1e6;
    531     Double_t up = E*10;
    532 
    533     Double_t binsize = log10(up/lo)/nbins;
    534     Double_t *scale=new Double_t[nbins+1];
    535     for (int i=0; i<nbins+1; i++)
    536         scale[i] = pow(10, binsize*i) * lo;
    537 
    538     MElectron epsilon;
    539     epsilon.SetEnergy(E);
    540 
    541 /*
    542     TH1F *hist = new TH1F("h", "h", nbins, scale);
    543 
    544     hist->Draw();
    545 
    546     TStopwatch clock;
    547     clock.Start();
    548 
    549     for (int i=0; i<1000; i++)
    550     {
    551         Double_t e1 = epsilon.GetRandom();
    552 
    553         hist->Fill(e1);
    554     }
    555 
    556     clock.Stop();
    557     clock.Print();
    558 
    559     epsilon.DrawCopy("same");
    560     gPad->SetLogx();
    561     gPad->SetLogy();
    562     */
    563 }
    564 
    565 void DrawILPhoton(Double_t z=0)
    566 {
    567     if (!gPad)
    568         new TCanvas("IL_Photon", "Mean Interaction Length Photon");
    569     else
    570         gPad->GetVirtCanvas()->cd(4);
    571 
    572     TF1 f1("length", MPhoton::InteractionLength, 1e4, 1e11, 1);
    573     f1.SetParameter(0, z);
    574 
    575     gPad->SetLogx();
    576     gPad->SetLogy();
    577     gPad->SetGrid();
    578     f1.SetMaximum(1e5);
    579     f1.SetLineWidth(1);
    580 
    581     TH1 *h=f1.DrawCopy()->GetHistogram();
    582 
    583     h->SetTitle("Mean Interaction Length (Photon)");
    584     h->SetXTitle("E [GeV]");
    585     h->SetYTitle("x [kpc]");
    586 
    587     gPad->Modified();
    588     gPad->Update();
    589 }
    590 
    591 void DrawILElectron(Double_t z=0)
    592 {
    593     if (!gPad)
    594         new TCanvas("IL_Electron", "Mean Interaction Length Electron");
    595     else
    596         gPad->GetVirtCanvas()->cd(3);
    597 
    598     TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0);
    599     f2.SetParameter(0, z);
    600 
    601     gPad->SetLogx();
    602     gPad->SetLogy();
    603     gPad->SetGrid();
    604     f2.SetLineWidth(1);
    605 
    606     TH1 *h=f2.DrawCopy()->GetHistogram();
    607 
    608     h->SetTitle("Mean Interaction Length (Electron)");
    609     h->SetXTitle("E [GeV]");
    610     h->SetYTitle("x [kpc]");
    611 
    612     gPad->Modified();
    613     gPad->Update();
    614 }
    615 
    616 void DrawRZ()
    617 {
    618     new TCanvas("RZ", "r and z");
    619 
    620     TF1 f1("ZofR", ZofR, 0, 4.5e5, 0);
    621     TF1 f2("RofZ", RofZ, 0, 5,     0);
    622 
    623     gPad->Divide(2,2);
    624 
    625     gPad->GetVirtCanvas()->cd(1);
    626     TH1 *h = f1.DrawCopy()->GetHistogram();
    627 
    628     h->SetTitle("z(r)");
    629     h->SetXTitle("r [kpc]");
    630     h->SetYTitle("z");
    631 
    632     gPad->Modified();
    633     gPad->Update();
    634 
    635     gPad->GetVirtCanvas()->cd(2);
    636     h = f2.DrawCopy()->GetHistogram();
    637 
    638     h->SetTitle("r(z)");
    639     h->SetXTitle("z");
    640     h->SetYTitle("r [kpc]");
    641 
    642     gPad->Modified();
    643     gPad->Update();
    644 }
    645 
    64624Double_t PrimSpect(Double_t *x, Double_t *k)
    64725{
    … …  
    64927}
    65028
    651 /*
    652 Double_t Planck(Double_t *x, Double_t *k=NULL)
    653 {
    654     Double_t E   = x[0];                    // [GeV]
    655     Double_t z   = k ? k[0] : 0;
    656 
    657     Double_t T   = 2.96*(z+1);              // [K]
    658     Double_t e   = 1.602176462e-19;         // [C]
    659     Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
    660 
    661     Double_t EkT = E/kB/T;
    662 
    663     return E*E*E / (exp(EkT)-1.); // [GeV^2]
    664 }
    665 */
    666 
    667 Double_t Planck(Double_t *x, Double_t *k=NULL)
    668 {
    669     //
    670     // Planck, per unit volume, per unit energy
    671     //
    672     // constants moved out of function, see below
    673     //
    674 
    675     //
    676     // This is to get a correct random value in a TF1!
    677     //
    678     Double_t E   = pow(10, x[0]);           // [GeV]
    679     Double_t z   = k ? k[0] : 0;
    680 
    681     Double_t T   = 2.96*(z+1);              // [K]
    682     Double_t e   = 1.602176462e-19;         // [C]
    683     Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
    684 
    685     Double_t EkT = E/kB/T;
    686 
    687     return E*E / (exp(EkT)-1.); // [GeV^2]
     29Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
     30{
     31    Double_t Ep  = pow(10, x[0]);
     32    Double_t res = MPhoton::Int2(&Ep, k);
     33
     34    return res*1e55; //65/k[0];
     35    // return MPhoton::Planck(&Ep, &k[1]);
     36}
     37
     38Double_t Sbar_sigmas(Double_t *x, Double_t *k)
     39{
     40    Double_t sbar = pow(10, x[0]);
     41
     42    Double_t s = 1./(sbar*4);
     43
     44    Double_t sigma = MPhoton::Sigma_gg(&s);
     45
     46    return sigma*sbar*1e28;
     47}
     48
     49Double_t RandomTheta(Double_t Eg, Double_t Ep)
     50{
     51    Double_t E0 = 511e-6; // [GeV]
     52
     53    Double_t f = Eg/E0*Ep/E0;
     54
     55    if (f<1)
     56        return 0;
     57
     58    TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
     59
     60    Double_t sbar  = pow(10, func.GetRandom());
     61    Double_t theta = acos(1.-sbar*2/f);
     62
     63    return theta;
    68864}
    68965
    69066void DoIt()
    69167{
    692     Double_t rcygnusx3 = 5000; //100; //30/3.258; // [~10kpc]
    693 
    694     Double_t startz = ZofR(&rcygnusx3);
     68    Double_t R = 100; // [kpc]
     69    Double_t startz = MParticle::ZofR(&R);
     70
     71    cout << "R = " << R << "kpc" << endl;
     72    cout << "Z = " << startz << endl;
    69573
    69674    static TRandom rand(0);
    69775    MPairProduction pair;
    69876
    699     //Double_t nphot = 50;
    700     Double_t runtime = 15*60; ///*18*60*/60; // [s]
    701 
    702     Double_t lo = 2e5;
    703     Double_t hi = 2e10;
     77    Double_t runtime = 30*60; // [s]
     78
     79    Double_t lo = 1e4;
     80    Double_t hi = 1e10;
    70481    Double_t alpha = -2;
    70582
    706     Double_t nbins = 4*log10(hi/lo);
    707 
    708     TF1 phot("PhotonSpectrum", Planck, -15, -10.5, 1);
     83    Double_t nbins = 24; //4*log10(hi/lo);
     84
    70985    TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
    710 
    711     src.SetParameter(0, alpha);
     86    src.SetParameter(0, 0);
     87
     88    TH1D hist;
     89    TH1D histsrc;
     90
     91    hist.SetMinimum(lo*lo * pow(lo, alpha)/10);
     92    histsrc.SetMinimum(lo*lo * pow(lo, alpha)/10);
    71293
    71394    TList listg;
    … …  
    726107    MBinning binspolx;
    727108    MBinning binspoly;
    728     binspolx.SetEdges(32, -180, 180);
    729     binspoly.SetEdgesLog(20, 1e-9, 1e-5);
     109    binspolx.SetEdges(16, -180, 180);
     110    binspoly.SetEdgesLog(20, 1e-13, 1e-4);
    730111    MH::SetBinning(&position, &binspolx, &binspoly);
    731112    MH::SetBinning(&angle,    &binspolx, &binspoly);
    … …  
    737118    histpos.SetTitle("Particle Position");
    738119    position.SetTitle("Particle Position");
    739 
    740     TH1D hist;
    741     TH1D histsrc;
    742120
    743121    //
    … …  
    761139    gPad->SetLogy();
    762140
     141    histsrc.SetName("Spectrum");
     142    histsrc.SetXTitle("E [GeV]");
     143    histsrc.SetYTitle("E^{2}\\dotCounts");
     144    histsrc.GetXaxis()->SetLabelOffset(-0.015);
     145    histsrc.GetXaxis()->SetTitleOffset(1.1);
     146
    763147    histsrc.Draw("P");
    764148    hist.Draw("Psame");
    765149    histsrc.Draw("Csame");
    766150    hist.Draw("Csame");
    767 
    768     histsrc.SetName("Spectrum");
    769     histsrc.SetXTitle("E [GeV]");
    770     histsrc.SetYTitle("E*Counts");
    771     histsrc.GetXaxis()->SetLabelOffset(-0.015);
    772     histsrc.GetXaxis()->SetTitleOffset(1.1);
    773151
    774152    ofstream fsrc("source.dat");
    … …  
    789167
    790168        Double_t E = pow(10, src.GetRandom());
    791 
    792         histsrc.Fill(E, E);
     169        Double_t weight = pow(E, alpha);
     170
     171        histsrc.Fill(E, E*E * weight);
     172
    793173        fsrc << E << endl;
    794174
    … …  
    800180        while (1)
    801181        {
    802             cout << " |P:" << flush;
     182            if (listg.GetSize()!=0)
     183                cout << " |P:" << flush;
    803184
    804185            TIter NextP(&listg);
    … …  
    806187            while ((p=(MPhoton*)NextP()))
    807188            {
     189                Double_t Eg = p->GetEnergy();
     190
    808191                if (!p->SetNewPosition())
    809192                {
    810193                    cout << "!" << flush;
    811194
    812                     hist.Fill(p->GetEnergy(), p->GetEnergy());
    813                     position.Fill(p->GetPhi()*kRad2Deg, p->GetR());
    814                     angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg);
    815                     histpos.Fill(p->GetR());
    816                     hista.Fill(p->GetTheta()*kRad2Deg);
    817 
    818                     fcas << p->GetEnergy() << endl;
     195                    hist.Fill(Eg, Eg*Eg*weight);
     196                    position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight);
     197                    angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
     198                    histpos.Fill(p->GetR(), weight);
     199                    hista.Fill(p->GetTheta()*kRad2Deg, weight);
     200
     201                    fcas << Eg << endl;
    819202                    delete listg.Remove(p);
    820203                    continue;
    … …  
    824207                // Sample phtoton from background and interaction angle
    825208                //
    826                 MPhoton photon;
    827 
    828                 phot.SetParameter(0, p->GetZ());
    829                 photon.SetEnergy(pow(10, phot.GetRandom()));
    830                 Double_t theta = rand.Uniform(TMath::Pi() * 2);
    831 
    832                 if (!pair.Process(p, &photon, theta, &liste))
    833                 {
    834                     cout << "0" << flush;;
     209                TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
     210                phot.SetParameter(0, Eg);
     211                phot.SetParameter(1, p->GetZ());
     212                if (phot.Integral(-log10(Eg)-8, -10.5)==0)
     213                {
     214                    // ???????????????????
     215                    hist.Fill(Eg, Eg*Eg*weight);
     216                    angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
     217                    hista.Fill(p->GetTheta()*kRad2Deg, weight);
     218                    delete listg.Remove(p);
     219                    cout << "z" << flush;
     220                    continue;
     221                }
     222
     223                Double_t Ep = pow(10, phot.GetRandom());
     224
     225                Double_t theta = RandomTheta(Eg, Ep);
     226                if (theta==0)
     227                {
     228                    cout << "t" << flush;
     229                    continue;
     230                }
     231
     232                if (!pair.Process(p, Ep, theta, &liste))
     233                {
     234                    cout << "0" << flush;
    835235                    continue;
    836236                }
    … …  
    843243                break;
    844244
    845             cout << " |E" << flush;
     245            if (liste.GetSize()!=0)
     246                cout << " |E" << flush;
    846247
    847248            TIter Next(&liste);
    … …  
    849250            while ((e=(MElectron*)Next()))
    850251            {
     252                if (e->GetEnergy()<lo)
     253                    continue;
     254
    851255                cout << ":" << flush;
    852256                while(1)
    … …  
    858262                    }
    859263
    860                     MPhoton *p = e->DoInvCompton();
    861                     if (e->GetEnergy()<lo*5)
     264                    // WRONG!
     265                    Double_t theta = rand.Uniform(TMath::Pi()*2);
     266                    MPhoton *p = e->DoInvCompton(theta);
     267
     268                    if (fabs(p->GetTheta()*3437)<1 &&  // < 1min
     269                        p->GetEnergy()>lo)
    862270                    {
    863                         cout << "0" << flush;
    864                         delete p;
    865                         break;
     271                        cout << "." << flush;
     272                        listg.Add(p);
    866273                    }
    867 
    868                     if (p->GetEnergy()<lo)
     274                    else
    869275                    {
    870276                        cout << "i" << flush; // ignored
    871277                        delete p;
     278                    }
     279
     280                    if (fabs(e->GetTheta()*3437)<1 &&  // < 1min
     281                        e->GetEnergy()>lo*2)
    872282                        continue;
    873                     }
    874 
    875                     listg.Add(p);
    876                     cout << "." << flush;
     283
     284                    cout << "0" << flush;
     285                    break;
    877286                }
    878287            }
    … …  
    940349    hista.DrawCopy();
    941350    gPad->SetLogx();
    942 
    943351}
    944352
    … …  
    946354void phys()
    947355{
     356    DoIt();
     357
    948358    /*
    949     Double_t Eg = 1e5;
    950 
    951     Double_t E0    = 511e-6;                  // [GeV]
    952     Double_t lolim = E0*E0/Eg;
    953     Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
    954 
    955     // 1e5:   5e-13, 4e-12
    956     // 1e6:   5e-14, 2e-12
    957     // 1e8:   5e-16, 1e-12
    958     // 1e10:  5e-18, 1e-12
    959 
    960     // 1e5:   2.6e-12, 4e-12
    961     // 1e6:   2.6e-13, 2e-12
    962     // 1e8:   2.6e-15, 1e-12
    963     // 1e10:  1.6e-17, 1e-12
    964 
    965     TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
    966     f.SetParameter(0, Eg);
    967     f.SetParameter(1, 0);
    968 
    969     MH::MakeDefCanvas();
    970     gPad->SetLogx();
    971     f.DrawCopy();
    972     return;
     359     Double_t E = 1e10;
     360     TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
     361     phot.SetParameter(0, E);
     362     phot.SetParameter(1, 5);
     363     phot.DrawCopy();
     364     return;
     365     */
     366
     367    /*
     368     Double_t Eg = 1e5;
     369
     370     Double_t E0    = 511e-6;                  // [GeV]
     371     Double_t lolim = E0*E0/Eg;
     372     Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
     373
     374     // 1e5:   5e-13, 4e-12
     375     // 1e6:   5e-14, 2e-12
     376     // 1e8:   5e-16, 1e-12
     377     // 1e10:  5e-18, 1e-12
     378
     379     // 1e5:   2.6e-12, 4e-12
     380     // 1e6:   2.6e-13, 2e-12
     381     // 1e8:   2.6e-15, 1e-12
     382     // 1e10:  1.6e-17, 1e-12
     383
     384     TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
     385     f.SetParameter(0, Eg);
     386     f.SetParameter(1, 0);
     387
     388     MH::MakeDefCanvas();
     389     gPad->SetLogx();
     390     f.DrawCopy();
    973391    */
    974 //    MH::MakeDefCanvas();
    975 //    DrawILPhoton();
    976 
    977 //    return ;
    978 
    979     DoIt();
    980 return;
    981     EnergyLossRate();
    982 
    983     MH::MakeDefCanvas();
    984     DrawILElectron();
    985 
    986     DrawRZ();
    987 
    988     return;
    989 
    990     /*
    991     cout << "Starting..." << endl;
    992 
    993     MEvtLoop   magic;
    994     MParList   plist;
    995     MTaskList  tlist;
    996     plist.AddToList(&tlist);
    997 
    998     MHInput  *input  = new MHInput;
    999     MHOutput *output = new MHOutput;
    1000 
    1001     plist.AddToList(input);
    1002     plist.AddToList(output);
    1003 
    1004     MGenerateInput geninp;
    1005     tlist.AddToList(&geninp);
    1006 
    1007     MFillH finput("MSinglePairInput", input);
    1008     tlist.AddToList(&finput);
    1009 
    1010     MSinglePairProduction prod;
    1011     tlist.AddToList(&prod);
    1012 
    1013     MFillH foutput("MSinglePairOutput", output);
    1014     tlist.AddToList(&foutput);
    1015 
    1016     magic.SetParList(&plist);
    1017     if (!magic.Eventloop(10000))
    1018         return;
    1019 
    1020     input->Draw();
    1021     output->Draw();
    1022 
    1023     return;
    1024 */
    1025 
    1026     /*
    1027     MEvtLoop   magic;
    1028     MParList   plist;
    1029     MTaskList  tlist;
    1030     plist.AddToList(&tlist);
    1031 
    1032     MHInput  *input  = new MHInput;
    1033     MHOutput *output = new MHOutput;
    1034 
    1035     plist.AddToList(input);
    1036     plist.AddToList(output);
    1037 
    1038     MGenerateInput geninp;
    1039     tlist.AddToList(&geninp);
    1040 
    1041     MFillH finput("MSinglePairInput", input);
    1042     tlist.AddToList(&finput);
    1043 
    1044     MSinglePairProduction prod;
    1045     tlist.AddToList(&prod);
    1046 
    1047     MFillH foutput("MSinglePairOutput", output);
    1048     tlist.AddToList(&foutput);
    1049 
    1050     magic.SetParList(&plist);
    1051     if (!magic.Eventloop(10000))
    1052         return;
    1053 
    1054     input->Draw();
    1055     output->Draw();
    1056 
    1057     return;
    1058    */
    1059 }
    1060 
     392}
     393
Note: See TracChangeset for help on using the changeset viewer.