Ignore:
Timestamp:
08/19/02 08:58:09 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/WuerzburgSoft/Thomas/mphys/Changelog

    r1476 r1507  
    11                                                                  -*-*- END -*-*-
     2 2002/08/19: Thomas Bretz
     3
     4   * MCascade.[h,cc]:
     5     - removed argument of DoInvCompton
     6     - delete GRandom only if set
     7     - set gRandom to 0 after deletion
     8
     9   * MElectron.[h,cc]:
     10     - added Sigma_ge
     11     - removed argument from DoInvCompton
     12     - Implemented Omega_sigmae and RandomThetaE
     13     - Changed old wrong dertermination of interaction angle to
     14       correct simulation.
     15     - Reengenieered some formulas to be more correct in terms of numerics
     16
     17
     18
    219 2002/08/02: Thomas Bretz
    320
  • trunk/WuerzburgSoft/Thomas/mphys/MCascade.cc

    r1476 r1507  
    5050}
    5151
    52 Double_t RandomTheta(Double_t Eg, Double_t Ep)
     52Double_t RandomThetaG(Double_t Eg, Double_t Ep)
    5353{
    5454    Double_t E0 = 511e-6; // [GeV]
     
    5959        return 0;
    6060
    61     static TF1 func("RndTheta", Sbar_sigmas, 0, 0, 0);
     61    static TF1 func("RndThetaG", Sbar_sigmas, 0, 0, 0);
    6262
    6363    func.SetRange(0, log10(f));
     
    139139        }
    140140
    141         // WRONG!
    142         MPhoton *p = e.DoInvCompton(0);
     141        MPhoton *p = e.DoInvCompton();
    143142
    144143        fBranchElectrons->GetTree()->Fill();
     
    200199
    201200        Double_t Ep = pow(10, pe);
    202         Double_t theta = RandomTheta(Eg, Ep);
     201        Double_t theta = RandomThetaG(Eg, Ep);
    203202        if (theta==0)
    204203        {
     
    316315MCascade::MCascade()
    317316{
     317    if (gRandom)
     318        delete gRandom;
     319
    318320    TRandom r(0);
    319     delete gRandom;
    320321    gRandom = new TRandom3(r.GetSeed());
    321322
     
    332333{
    333334    delete gRandom;
     335    gRandom = 0;
    334336}
    335337
     
    373375
    374376    // ------------------------------
     377
     378    cout << endl;
    375379
    376380    cout << "R = " << fSrcR << "kpc" << endl;
  • trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

    r1449 r1507  
    4747ClassImp(MElectron);
    4848
     49Double_t MElectron::Sigma_ge(Double_t *x, Double_t *k)
     50{
     51    const Double_t E0 = 511e-6;                 // [GeV]
     52    const Double_t c  = 299792458;              // [m/s]
     53    const Double_t e  = 1.602176462e-19;        // [C]
     54    const Double_t h  = 1e-9/e*6.62606876e-34;  // [GeVs]
     55    const Double_t a  = 1./137;                 // [1]
     56
     57    const Double_t o = x[0];
     58
     59    const Double_t o1  = o+1;
     60    const Double_t o21 = o*2+1;
     61
     62    const Double_t s1 = TMath::Pi()*2;
     63    const Double_t s2 = a*h*c/E0;             // [m]
     64    const Double_t s3 = o1/(o*o*o);
     65    const Double_t s4 = o*2*o1/o21;
     66    const Double_t s5 = log(o21);
     67    const Double_t s6 = s5/o/2;
     68    const Double_t s7 = (o*3+1)/(o21*o21);
     69
     70    return s2*s2/s1*(s3*(s4-s5)+s6-s7); // [m^2]
     71}
     72
    4973Double_t MElectron::Li(Double_t *x, Double_t *k)
    5074{
     
    310334}
    311335
    312 MPhoton *MElectron::DoInvCompton(Double_t theta)
     336
     337Double_t Omega_sigmae(Double_t *x, Double_t *k)
     338{
     339    Double_t sbar = pow(10,x[0]);
     340
     341    Double_t omega = (sbar-1)/2;
     342
     343    Double_t sigma = MElectron::Sigma_ge(&omega);
     344
     345    return (sbar-1)*sigma*1e28;
     346}
     347
     348Double_t RandomThetaE(Double_t Ee, Double_t Ep)
     349{
     350    Double_t E0 = 511e-6; // [GeV]
     351
     352    Double_t f = 2*Ee/E0*Ep/E0;
     353
     354    static TF1 func("RndThetaE", Omega_sigmae, 0, 0, 0);
     355
     356    Double_t beta = sqrt(1-E0/Ee*E0/Ee);
     357
     358    //func.SetRange(0, log10(1+f*(1+beta)));
     359    func.SetRange(log10(1+f*(1-beta)), log10(1+f*(1+beta)));
     360    func.SetNpx(50);
     361
     362    Double_t sbar = pow(10, func.GetRandom());
     363
     364    Double_t bcost = 1 - (sbar-1)/f;
     365    return bcost;
     366
     367    Double_t theta = acos(bcost/beta);
     368
     369    return theta;
     370}
     371
     372MPhoton *MElectron::DoInvCompton()
    313373{
    314374    const Double_t E0 = 511e-6; //[GeV]
     
    317377    const Double_t e = GetEnergyLoss(&epsilon);
    318378
    319     // er: photon energy before interaction, rest frame
    320     // e:  photon energy after interaction, lab
    321 
    322     const Double_t gamma = fEnergy/E0;
    323     // const Double_t beta  = sqrt(1.-1./(gamma*gamma));
     379    // epsilon: photon energy befor interaction, lab
     380    // e:       photon energy after interaction, lab
     381
     382    const Double_t gamma     = fEnergy/E0;
    324383    const Double_t gammabeta = sqrt(gamma*gamma-1);
    325 
    326     const Double_t f = fEnergy/e;
    327 
    328     Double_t t=-10;
    329     Double_t arg;
    330     do
    331     {
    332         if (t>-10)
    333             cout << "~" << flush;
    334 
    335         //
    336         // Create an angle uniformly distributed in the range of possible
    337         // interaction angles due to the known energies.
    338         //
    339         // The distibution is a theta-function which is incorrect, but
    340         // should be correct in a first order...
    341         //
    342         t = gRandom->Uniform(TMath::Pi())+TMath::Pi()/2;
    343         Double_t er = epsilon*(gamma-gammabeta*cos(t)); // photon energy rest frame
    344         Double_t n  = sqrt(fEnergy*fEnergy-E0*E0)/e+1;
    345         arg = (f - E0/er - 1)/n;
    346 
    347         /*  ------------------ some hints ------------------------------
    348          -1 < arg < 1
    349          -n < f - r/er - 1 < n
    350          1-f-n < r/er < 1-f+n
    351          r/(1-f+n) < er < r/(1-f-n)
    352          r/(1-f+n) < gamma-gammabeta*cos(t) < r/(1-f-n)
    353          r/(1-f+n)-gamma < -gammabeta*cos(t) < r/(1-f-n)-gamma
    354          gamma-r/(1-f-n) < gammabeta*cos(t) < gamma-r/(1-f+n)
    355          (gamma-r/(1-f-n))/gammabeta < cos(t) < (gamma-r/(1-f+n))/gammabeta
    356          acos((gamma-r/(1-f+n))/gammabeta) < t < acos((gamma-r/(1-f-n))/gammabeta)
    357 
    358 
    359          (gamma+E0/(n*arg+1-f)/epsilon)/gammabeta = cos(t);
    360          1 = (epsilon/E0*cos(t)*gammabeta-gamma)*(n*arg+1-f);
    361 
    362          arg=1:   (gamma+E0/epsilon*(1-f+n))/gammabeta = cos(t);
    363          arg=-1 : (gamma+E0/epsilon*(1-f-n))/gammabeta = cos(t);
    364 
    365          (E0/(1-f-n)/epsilon+gamma)/gammabeta < cos(t) < (E0/(1-f+n)/epsilon+gamma)/gammabeta
    366          */
    367 
    368     } while (arg<-1 || arg>1);
     384    const Double_t beta      = gammabeta/gamma;
     385
     386    const Double_t bcost = RandomThetaE(fEnergy, epsilon);
     387    const Double_t t     = acos(bcost/beta);
     388
     389    const Double_t f = epsilon/fEnergy;
     390    const Double_t r = gamma*(1-bcost);
     391
     392    Double_t arg = (1 - 1/(gamma*r) - f)/(beta-f);
     393
     394    if (arg<-1 || arg>1)
     395        cout << "<" << (int)(t*180/TMath::Pi()) << "°>" << flush;
     396
     397    //
     398    // Due to numerical uncertanties arg can be something like:
     399    // 1+1e-15 which we cannot allow
     400    //
     401    if (arg<-1)
     402        arg = -1;
     403    if (arg>1)
     404        arg = 1;
    369405
    370406    const Double_t theta1s = acos(arg);
    371     const Double_t thetas  = atan(sin(t)/(gamma*cos(t)-gammabeta));
     407    const Double_t thetas  = atan(-sin(t)/(beta*r)/*(1-bcost)/gammabeta*/);
    372408
    373409    const Double_t thetastar = thetas-theta1s;
    374410
    375     const Double_t theta1 = atan(sin(thetastar)/(gamma*cos(thetastar)+gammabeta));
     411    //    const Double_t theta1 = atan(sin(thetastar)/(gamma*cos(thetastar)+gammabeta));
     412    const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
    376413
    377414    fEnergy -= e;
     
    379416    const Double_t phi = gRandom->Uniform(TMath::Pi()*2);
    380417
    381     /*
    382     MPhoton &p = *new MPhoton(e, fZ);
    383     p = *this;
    384     */
    385418    MPhoton &p = *new MPhoton(*this, e);
    386419    p.SetNewDirection(theta1, phi);
    387420
    388     const Double_t beta2  = sqrt(1.-E0/fEnergy*E0/fEnergy);
    389     const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2);
     421    /*
     422     const Double_t beta2  = sqrt(1.-E0/fEnergy*E0/fEnergy);
     423     const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2);
     424     */
     425
     426    const Double_t div    = sqrt(fEnergy*fEnergy-E0*E0);
     427    const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/div);
    390428
    391429    SetNewDirection(theta2, phi);
  • trunk/WuerzburgSoft/Thomas/mphys/MElectron.h

    r1449 r1507  
    2424    // ----------------------------------------------------------------
    2525
     26    static Double_t Sigma_ge(Double_t *x, Double_t *k=NULL);
     27
     28    // ----------------------------------------------------------------
     29
    2630    static Double_t DiSum(Double_t *x, Double_t *k=NULL);
    2731    static Double_t Li(Double_t *x, Double_t *k);
     
    4549    // ----------------------------------------------------------------
    4650
    47     MPhoton *DoInvCompton(Double_t theta);
     51    MPhoton *DoInvCompton();
    4852    Bool_t SetNewPositionB(Double_t B);
    4953
  • trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h

    r1449 r1507  
    1919    }
    2020
    21         MPhoton(MParticle &p) : MParticle(p, MParticle::kEGamma) { }
    22         MPhoton(MParticle &p, Double_t e) : MParticle(p, e, MParticle::kEGamma) { }
     21    MPhoton(MParticle &p) : MParticle(p, MParticle::kEGamma) { }
     22    MPhoton(MParticle &p, Double_t e) : MParticle(p, e, MParticle::kEGamma) { }
    2323
    2424    void operator=(MParticle &p) { MParticle::operator=(p); }
Note: See TracChangeset for help on using the changeset viewer.