Ignore:
Timestamp:
06/13/02 16:36:19 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/WuerzburgSoft/Thomas/mphys/Changelog

    r1360 r1363  
    44   * MParticle.h:
    55     - Implemented Set/IsPrimary
     6     - removed pure abstract method (root can't read this)
     7     - added default constructor
     8
     9   * MElectron.cc:
     10     - changed many variables to const
     11     - changed the output
     12
     13   * MPairProduction.cc:
     14     - changed many variables to const
     15
     16   * MPhoton.[h,cc]:
     17     - changed the upper and lower limit of the integration in Int2
     18     - added upper and lower limit of draw to DrawInteractionLength
    619
    720   * energyloss.C:
    … …  
    1124     - changed energy randomizer
    1225     - added root file output
     26     - removed log scale from point spread plots
     27
     28   * analp.C:
     29     - added macro to analize the photons from the root file
     30
     31   * anale.C:
     32     - added macro to analize the electrons from the root file
    1333
    1434

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

    r1358 r1363  
    4949    // constants moved out of function
    5050    //
    51     Double_t E   = x[0];                    // [GeV]
    52     Double_t z   = k ? k[0] : 0;
    53 
    54     Double_t T   = 2.96*(z+1);              // [K]
    55     Double_t e   = 1.602176462e-19;         // [C]
    56     Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
    57 
    58     Double_t EkT = E/kB/T;
     51    const Double_t E   = x[0];                    // [GeV]
     52    const Double_t z   = k ? k[0] : 0;
     53
     54    const Double_t T   = 2.96*(z+1);              // [K]
     55    const Double_t e   = 1.602176462e-19;         // [C]
     56    const Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
     57
     58    const Double_t EkT = E/kB/T;
    5959
    6060    /*
    … …  
    7171Double_t MElectron::Li(Double_t *x, Double_t *k)
    7272{
    73     Double_t t = x[0];
     73    const Double_t t = x[0];
    7474
    7575    return log(1.-t)/t;
    … …  
    8484
    8585    TF1 IntLi("Li", Li, 0, z, 0);
    86     Double_t integ = IntLi.Integral(0, z);
     86    const Double_t integ = IntLi.Integral(0, z);
    8787
    8888    /*
    … …  
    105105Double_t MElectron::Flim(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam)  mit  b-->1   (Maple)
    106106{
    107     Double_t w = x[0];
    108 
    109     Double_t w4   = w*4;
    110     Double_t wsqr = w*w;
    111 
    112     Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
    113     Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
    114     Double_t d  = w4*(w4 + 1);
     107    const Double_t w = x[0];
     108
     109    const Double_t w4   = w*4;
     110    const Double_t wsqr = w*w;
     111
     112    const Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
     113    const Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
     114    const Double_t d  = w4*(w4 + 1);
    115115
    116116    Double_t s   = -w*2*(1+1); // -2*omega*(1+beta)
    117     Double_t li2 = Li2(&s);
    118 
    119     Double_t res = (u1+u2)/d + li2;
     117    const Double_t li2 = Li2(&s);
     118
     119    const Double_t res = (u1+u2)/d + li2;
    120120
    121121    return res; //<1e-10? 0 : res;
    … …  
    124124Double_t MElectron::Compton(Double_t *x, Double_t *k)
    125125{
    126     Double_t E0  = 511e-6; //[GeV]
    127     Double_t E02 = E0*E0;
     126    const Double_t E0  = 511e-6; //[GeV]
     127    const Double_t E02 = E0*E0;
    128128
    129129    Double_t epsilon = x[0];
    … …  
    134134    Double_t omega  = epsilon*E/E02;
    135135
    136     Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]
     136    const Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]
    137137    return Flim(&omega)*n;
    138138}
    … …  
    160160    // F(o-)~F(0)  = 14/8 = 1.75
    161161
    162     Double_t E0     = 511e-6;                        // [GeV]
    163     Double_t E02    = E0*E0;                         // [GeV^2]
    164     Double_t c      = 299792458;                     // [m/s]
    165     Double_t e      = 1.602176462e-19;               // [C]
    166     Double_t h      = 1e-9/e*6.62606876e-34;         // [GeVs]
    167     Double_t hc     = h*c;                           // [GeVm]
    168     Double_t alpha  = 1./137.;                       // [1]
    169 
    170     Double_t z      = k ? k[0] : 0;
     162    const Double_t E0     = 511e-6;                        // [GeV]
     163    const Double_t E02    = E0*E0;                         // [GeV^2]
     164    const Double_t c      = 299792458;                     // [m/s]
     165    const Double_t e      = 1.602176462e-19;               // [C]
     166    const Double_t h      = 1e-9/e*6.62606876e-34;         // [GeVs]
     167    const Double_t hc     = h*c;                           // [GeVm]
     168    const Double_t alpha  = 1./137.;                       // [1]
     169
     170    const Double_t z      = k ? k[0] : 0;
    171171
    172172    // Double_t beta = sqrt(1-E0/E*E0/E);
    173     Double_t beta   = 1; //0.999999999999999999999999999;
     173    const Double_t beta   = 1; //0.999999999999999999999999999;
    174174
    175175    Double_t val[3] = { E[0], beta, z };             // E[GeV]
    176176
    177     Double_t from   = 1e-17;
    178     Double_t to     = 1e-11;
     177    const Double_t from   = 1e-17;
     178    const Double_t to     = 1e-11;
    179179
    180180    /* -------------- old ----------------
    … …  
    186186    TF1 func("Compton", Compton, from, to, 3);      // [0, inf]
    187187
    188     Double_t integ  = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]
    189 
    190     Double_t aE     = alpha/E[0];                   // [1/GeV]
    191 
    192     Double_t konst  = 2.*E02/hc/beta;               // [1 / GeV m]
    193     Double_t ret    = konst * (aE*aE) * integ;      // [1 / m]
    194 
    195     Double_t ly     = 3600.*24.*365.*c;             // [m/ly]
    196     Double_t pc     = 1./3.258;                     // [pc/ly]
     188    const Double_t integ  = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]
     189
     190    const Double_t aE     = alpha/E[0];                   // [1/GeV]
     191
     192    const Double_t konst  = 2.*E02/hc/beta;               // [1 / GeV m]
     193    const Double_t ret    = konst * (aE*aE) * integ;      // [1 / m]
     194
     195    const Double_t ly     = 3600.*24.*365.*c;             // [m/ly]
     196    const Double_t pc     = 1./3.258;                     // [pc/ly]
    197197
    198198    return (1./ret)/ly*pc/1000;                     // [kpc]
    … …  
    214214{
    215215    Double_t e  = pow(10, x[0]);
    216     Double_t E  = k[0];
    217216    Double_t z  = k[1];
    218217
    219     Double_t E0  = 511e-6; //[GeV]
    220     Double_t E02 = E0*E0;
     218    const Double_t E  = k[0];
     219
     220    const Double_t E0  = 511e-6; //[GeV]
     221    const Double_t E02 = E0*E0;
    221222
    222223    Double_t omega = e*E/E02;
    223224
    224     Double_t n = Planck(&e, &z);
    225 
    226     Double_t F = Flim(&omega)/omega/omega;
     225    const Double_t n = Planck(&e, &z);
     226
     227    const Double_t F = Flim(&omega)/omega/omega;
    227228
    228229    return n*F*1e26;
    … …  
    231232Double_t MElectron::G_q(Double_t *x, Double_t *k)
    232233{
    233     Double_t q     = x[0];
    234     Double_t Gamma = k[0];
    235 
    236     Double_t Gq = Gamma*q;
    237 
    238     Double_t s1 = 2.*q*log(q);
    239     Double_t s2 = (1.+2.*q);
    240     Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;
     234    const Double_t q     = x[0];
     235    const Double_t Gamma = k[0];
     236
     237    const Double_t Gq = Gamma*q;
     238
     239    const Double_t s1 = 2.*q*log(q);
     240    const Double_t s2 = (1.+2.*q);
     241    const Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;
    241242
    242243    return s1+(s2+s3)*(1.-q);
    … …  
    246247Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
    247248{
    248     Double_t E = x[0];
    249     Double_t z = k ? k[0] : 0;
    250 
    251     Double_t E0 = 511e-6; //[GeV]
    252 
    253     Double_t lolim = -log10(E)/7*4-13;
     249    const Double_t E = x[0];
     250    const Double_t z = k ? k[0] : 0;
     251
     252    const Double_t E0 = 511e-6; //[GeV]
     253
     254    const Double_t lolim = -log10(E)/7*4-13;
    254255
    255256    TF1 fP("p_e", p_e, lolim, -10.8, 2);
    … …  
    259260    fP.SetParameter(1, z);
    260261
    261     Double_t e = pow(10, fP.GetRandom());
     262    const Double_t e = pow(10, fP.GetRandom());
    262263
    263264    if (ep)
    264265        *ep = e;
    265266
    266     Double_t omega = e*E/E0/E0;
    267     Double_t Gamma = 4.*omega;
     267    const Double_t omega = e*E/E0/E0;
     268    const Double_t Gamma = 4.*omega;
    268269
    269270    // --old-- fQ.SetRange(1e-6, 1./(1.+ 2.*Gamma));
    270271    fQ.SetParameter(0, Gamma);
    271272
    272     Double_t q  = fQ.GetRandom();
    273     Double_t Gq = Gamma*q;
    274 
    275     Double_t e1 = Gq*E/(1.+Gq);
     273    const Double_t q  = fQ.GetRandom();
     274    const Double_t Gq = Gamma*q;
     275
     276    const Double_t e1 = Gq*E/(1.+Gq);
    276277
    277278    return e1;
    … …  
    292293    static TRandom rand(0);
    293294
    294     Double_t E0 = 511e-6; //[GeV]
     295    const Double_t E0 = 511e-6; //[GeV]
    295296
    296297    Double_t epsilon;
    297     Double_t e = GetEnergyLoss(&epsilon);
     298    const Double_t e = GetEnergyLoss(&epsilon);
    298299
    299300    // er: photon energy before interaction, rest frame
    300301    // e:  photon energy after interaction, lab
    301302
    302     Double_t gamma     = fEnergy/E0;
    303     Double_t beta      = sqrt(1.-1./(gamma*gamma));
    304     //Double_t gammabeta = sqrt(gamma*gamma-1);
    305 
    306     Double_t f = fEnergy/e;
    307 
    308     Double_t t;
     303    const Double_t gamma     = fEnergy/E0;
     304    const Double_t beta      = sqrt(1.-1./(gamma*gamma));
     305
     306    const Double_t f = fEnergy/e;
     307
     308    Double_t t=-1;
    309309    Double_t arg;
    310310    do
    311311    {
    312         t = rand.Uniform(TMath::Pi()*2);
     312        if (t>0)
     313            cout << "~" << flush;
     314        t = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
    313315        Double_t er  = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame
    314316        arg = (f - E0/er - 1)/(f*beta+1);
    315         cout << "~" << flush;
    316317
    317318    } while (arg<-1 || arg>1);
    318319
    319     Double_t theta1s = acos(arg);
    320     Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));
    321 
    322     Double_t thetastar = thetas-theta1s;
    323 
    324     Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
    325 
    326     /*
    327      cout << "(" << theta1/TMath::Pi()*180 << ",";
    328      cout << theta1s/TMath::Pi()*180<< ",";
    329      cout << arg << ")" << flush;
    330      */
     320    const Double_t theta1s = acos(arg);
     321    const Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));
     322
     323    const Double_t thetastar = thetas-theta1s;
     324
     325    const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
    331326
    332327    fEnergy -= e;
    … …  
    335330    p = *this;
    336331    p.SetNewDirection(theta1, rand.Uniform(TMath::Pi()*2));
    337 
    338     // MISSING: Electron angle
    339     //
    340     // E1 = fEnergy  (after!)
    341     //
    342     // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0)
    343332
    344333    return &p;

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

    r1358 r1363  
    7979    //  theta: interaction angle [rad]
    8080    //
    81 
    82 
    8381    const Double_t E0     = 511e-6;              // [GeV]
    8482    const Double_t Eg     = gamma->GetEnergy();  // [GeV]
    … …  
    9896    const Double_t GammaH = (Eg+Ep)/sqrt(s);
    9997
    100     Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));
     98    const Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));
    10199
    102100    fAngle->SetParameter(0, sqrt(sqrbetae));
    … …  
    104102    const Double_t alpha = psi-acos(fAngle->GetRandom());
    105103
    106     Double_t salpha = sin(alpha);
    107     Double_t calpha = cos(alpha);
     104    const Double_t salpha = sin(alpha);
     105    const Double_t calpha = cos(alpha);
    108106
    109107    const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi)
    110108
    111     Double_t bb = sqrt(sqrbetah/sqrbetae);
     109    const Double_t bb = sqrt(sqrbetah/sqrbetae);
    112110
    113     Double_t s1 = calpha/GammaH;
    114     Double_t s2 = tphi*s1 - salpha - bb;
     111    const Double_t s1 = calpha/GammaH;
     112    const Double_t s2 = tphi*s1 - salpha - bb;
    115113
    116     Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
    117     Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;
     114    const Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
     115    const Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;
    118116
    119117    const Double_t E = (Eg+Ep)/2;;

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

    r1356 r1363  
    152152    Double_t val[2] = { Eg, z };
    153153
    154     Double_t lolim = E0*E0 > 1e-8 ? E0*E0/Eg : 1e-8/Eg;
    155     Double_t inf   = 3e-11;
     154    Double_t lolim = E0*E0/Eg;
     155    Double_t inf   = Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1);
    156156
    157157    TF1 f("int2", Int2, lolim, inf, 2);
    … …  
    188188}
    189189
    190 void MPhoton::DrawInteractionLength(Double_t z)
     190void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to)
    191191{
    192192    if (!gPad)
    … …  
    195195        gPad->GetVirtCanvas()->cd(4);
    196196
    197     TF1 f1("length", InteractionLength, 1e4, 1e11, 1);
     197    TF1 f1("length", InteractionLength, from, to, 1);
    198198    f1.SetParameter(0, z);
    199199

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

    r1360 r1363  
    3232    // ----------------------------------------------------------------
    3333
    34     static void DrawInteractionLength(Double_t z);
     34    static void DrawInteractionLength(Double_t z, Double_t from=1e4, Double_t to=1e11);
    3535    void DrawInteractionLength() const;
    3636

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

    r1362 r1363  
    44{
    55    TChain chain("Electrons");
    6     chain.Add("cascade01.root");
    7     chain.Add("cascade02.root");
     6    chain.Add("cascade_100kpc_01.root");
     7    chain.Add("cascade_100kpc_02.root");
     8    chain.Add("cascade_100kpc_03.root");
    89
    910    MElectron *e = new MElectron;

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

    r1362 r1363  
    44{
    55    TChain chain("Photons");
    6     chain.Add("cascade01.root");
    7     chain.Add("cascade02.root");
     6    chain.Add("cascade_100kpc_01.root");
     7    chain.Add("cascade_100kpc_02.root");
     8    chain.Add("cascade_100kpc_03.root");
    89
    910    MPhoton *p = new MPhoton;

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

    r1360 r1363  
    477477// -------------------------------------------------------------------
    478478
     479Double_t func(Double_t *x, Double_t *k)
     480{
     481    return MPhoton::Int2(x, k)*1e68;
     482}
     483
    479484void energyloss()
    480485{
     486/*    Double_t E0    = 511e-6;                // [GeV]
     487
     488    Double_t Eg = 1e4; //3.6e4;
     489    Double_t z  = 5;
     490
     491    Double_t val[2] = { Eg, z };
     492
     493    Double_t lolim = E0*E0/Eg;
     494    Double_t inf   = Eg<1e6 ? 3e-11*z : 3e-12*z;
     495
     496    cout << Eg << " " << z << " " << lolim << " " << inf << endl;
     497
     498    TF1 f("int2", func, lolim, inf, 2);
     499
     500    Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2]
     501
     502    f.SetParameter(0, Eg);
     503    f.SetParameter(1, z);
     504
     505    cout << int2 << endl;
     506
     507    new TCanvas("ILPhoton", "Mean Interaction Length Photon");
     508
     509    gPad->SetLogx();
     510//    gPad->SetLogy();
     511    gPad->SetGrid();
     512
     513    f.SetLineWidth(1);
     514    f.DrawCopy();
     515
     516    gPad->Modified();
     517    gPad->Update();
     518
     519    //    return;
     520  */  MPhoton p;
     521    p.DrawInteractionLength();
     522
     523    return;
    481524//    EnergyLossRate();
    482525

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

    r1361 r1363  
    116116void DoIt()
    117117{
    118     //    Double_t R = 1798; // [kpc]
    119     Double_t startz = 0.003; //MParticle::ZofR(&R);
    120 
    121     //    cout << "R = " << R << "kpc" << endl;
     118    Double_t R = 500; // [kpc]
     119    Double_t startz = MParticle::ZofR(&R);
     120
     121    cout << "R = " << R << "kpc" << endl;
    122122    cout << "Z = " << startz << endl;
    123123
    … …  
    131131    Double_t alpha = -2;
    132132
    133     Int_t nbins = 18; //4*log10(hi/lo);
     133    Int_t nbins = 21; //4*log10(hi/lo);
    134134
    135135    TH1D hist;
    … …  
    153153
    154154    MBinning binspolx;
    155     MBinning binspoly;
     155    MBinning binspoly1;
     156    MBinning binspoly2;
    156157    binspolx.SetEdges(16, -180, 180);
    157     binspoly.SetEdgesLog(20, 1e-10, 1e-3);
    158     MH::SetBinning(&position, &binspolx, &binspoly);
    159     MH::SetBinning(&angle,    &binspolx, &binspoly);
    160     MH::SetBinning(&histpos,  &binspoly);
    161     MH::SetBinning(&hista,    &binspoly);
     158    binspoly1.SetEdges(20, 0, 5e-9);
     159    binspoly2.SetEdges(20, 0, 2e-7);
     160    MH::SetBinning(&angle,    &binspolx, &binspoly1);
     161    MH::SetBinning(&position, &binspolx, &binspoly2);
     162    MH::SetBinning(&hista,    &binspoly1);
     163    MH::SetBinning(&histpos,  &binspoly2);
    162164
    163165    hista.SetTitle("Particle Angle");
    … …  
    257259                    hista.Fill(p->GetTheta()*kRad2Deg, weight);
    258260
     261                    p->SetIsPrimary(kFALSE);
    259262                    T->Fill();
    260263
    … …  
    393396    position.SetYTitle("Pos: R [kpc]");
    394397    position.DrawCopy("surf1pol");
    395     gPad->SetLogy();
     398    //gPad->SetLogy();
    396399
    397400    MH::MakeDefCanvas();
    … …  
    399402    angle.SetYTitle("Angle: \\Theta [\\circ]");
    400403    angle.DrawCopy("surf1pol");
    401     gPad->SetLogy();
     404    //gPad->SetLogy();
    402405
    403406    MH::MakeDefCanvas();
Note: See TracChangeset for help on using the changeset viewer.