Changeset 1352 for trunk/WuerzburgSoft


Ignore:
Timestamp:
06/10/02 08:49:28 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

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

    r1349 r1352  
    3434#include <TF1.h>
    3535
     36#include "MPhoton.h"
     37
    3638ClassImp(MElectron);
    3739
    … …  
    238240
    239241
    240 Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k = NULL)
     242Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
    241243{
    242244    Double_t E = x[0];
    … …  
    255257    Double_t e = pow(10, fP.GetRandom());
    256258
     259    if (ep)
     260        *ep = e;
     261
    257262    Double_t omega = e*E/E0/E0;
    258263    Double_t Gamma = 4.*omega;
    … …  
    269274}
    270275
    271 Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z)
     276Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep)
    272277{
    273278    return EnergyLoss(&E, &z);
    274279}
    275280
    276 Double_t MElectron::GetEnergyLoss() const
    277 {
    278     return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ);
    279 }
     281Double_t MElectron::GetEnergyLoss(Double_t *ep) const
     282{
     283    return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep);
     284}
     285
     286#include <iostream.h>
     287
     288MPhoton *MElectron::DoInvCompton()
     289{
     290    Double_t E0 = 511e-6; //[GeV]
     291
     292    Double_t epsilon;
     293    Double_t e = GetEnergyLoss(&epsilon);
     294    /*
     295     Double_t gamma = fEnergy/E0;
     296     Double_t beta = sqrt(1-1/(gamma*gamma));
     297     Double_t bega = sqrt(gamma*gamma-1);
     298
     299     Double_t theta = TMath::Pi()/4;
     300
     301     // er: photon energy before interaction rest frame
     302     // e:  photon energy after interaction lab
     303     Double_t er = gamma*epsilon*(1-beta*cos(theta)); // photon energy rest frame
     304     Double_t r1 = fEnergy/e;
     305     Double_t r2 = E0/er;
     306     Double_t ctheta = (r1-r2-1)/(beta*r1-1);
     307     Double_t tg = sqrt(1-ctheta*ctheta)/gamma/(beta+ctheta);
     308    */
     309    fEnergy -= e;
     310
     311    MPhoton &p = *new MPhoton(e, fZ);
     312    p = *this;
     313
     314    /*
     315     cout << "(" << atan(tg)*180/TMath::Pi() << ")" << flush;
     316
     317     static TRandom rand(0);
     318     p->SetNewDirection(atan(tg), rand.Uniform(TMath::Pi()*2));
     319     */
     320
     321    // MISSING: Electron angle
     322    //
     323    // E1 = fEnergy  (after!)
     324    //
     325    // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0)
     326
     327    return &p;
     328}

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

    r1349 r1352  
    55#include "MParticle.h"
    66#endif
     7
     8class MPhoton;
    79
    810class MElectron : public MParticle
    … …  
    1416        fZ      = z;
    1517    }
     18
     19    void operator=(MParticle &p) { MParticle::operator=(p); }
    1620
    1721    // ----------------------------------------------------------------
    … …  
    3135    static Double_t p_e(Double_t *x, Double_t *k);
    3236    static Double_t G_q(Double_t *x, Double_t *k);
    33     static Double_t EnergyLoss(Double_t *E, Double_t *z=NULL);
    34     static Double_t GetEnergyLoss(Double_t E, Double_t z=0);
     37    static Double_t EnergyLoss(Double_t *E, Double_t *z=NULL, Double_t *ep=NULL);
     38    static Double_t GetEnergyLoss(Double_t E, Double_t z=0, Double_t *ep=NULL);
    3539
    36     Double_t GetEnergyLoss() const;
     40    Double_t GetEnergyLoss(Double_t *ep) const;
     41
     42    MPhoton *DoInvCompton();
    3743
    3844    ClassDef(MElectron, 1)

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

    r1349 r1352  
    3434
    3535#include "MParList.h"
    36 #include "MParticle.h"
     36#include "MPhoton.h"
    3737
    3838ClassImp(MGenIRPhoton);
    … …  
    104104{
    105105    const Double_t pi2 = TMath::Pi() * 2;
    106     MParticle *p = new MParticle(MParticle::kEGamma);
    107106
    108     p->SetAngle(fRand.Uniform(pi2));
    109     p->SetEnergy(fSrc->GetRandom()*1e-9);
     107    MPhoton *p = new MPhoton;
     108    /*
     109     MParticle *p = new MParticle(MParticle::kEGamma);
     110     p->SetAngle(fRand.Uniform(pi2));
     111     p->SetEnergy(fSrc->GetRandom()*1e-9);
     112     */
    110113
    111114    return p;

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

    r1349 r1352  
    7171Bool_t MGenPrimaryParticle::Process()
    7272{
    73     MParticle *p = new MParticle(MParticle::kEGamma);
    74     p->SetEnergy(fSrc->GetRandom());
     73//    MParticle *p = new MParticle(MParticle::kEGamma);
     74//    p->SetEnergy(fSrc->GetRandom());
    7575    return kTRUE;
    7676}

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

    r1349 r1352  
    3232#include <TF1.h>
    3333#include <TMath.h>
     34#include <TRandom.h>
    3435
    3536#include "TList.h"
    … …  
    7071
    7172// --------------------------------------------------------------------------
    72 Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, TList *list)
     73Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, const Double_t theta, TList *list)
    7374{
    7475    //
    … …  
    7879    const Double_t E0     = 511e-6;              // [GeV]
    7980
    80     const Double_t theta  = phot->GetAngle();    // [2pi]
     81    //const Double_t theta  = phot->GetAngle();    // [2pi]
    8182    const Double_t Ep     = phot->GetEnergy();   // [GeV]
    8283    const Double_t Eg     = gamma->GetEnergy();  // [GeV]
    … …  
    120121    const Double_t dE = E*Bcosd;
    121122
    122     MElectron *p[2];
    123     p[0] = new MElectron(E+dE, gamma->GetZ());
    124     p[1] = new MElectron(E-dE, gamma->GetZ());
     123    MElectron &p0 = *new MElectron(E+dE, gamma->GetZ());
     124    MElectron &p1 = *new MElectron(E-dE, gamma->GetZ());
     125    p0 = *gamma; // Set Position and direction
     126    p1 = *gamma; // Set Position and direction
    125127
    126128    const Double_t E1 = E0/(E+dE);
    … …  
    130132    const Double_t beta2 = sqrt(1.-E2*E2);
    131133
    132     p[0]->SetBeta(beta1);
    133     p[1]->SetBeta(beta2);
     134    p0.SetBeta(beta1);
     135    p1.SetBeta(beta2);
    134136
    135137    const Double_t Bscp = Bsind*cphi;
    … …  
    140142    const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));
    141143
    142     p[0]->SetAngle(atan(tan1));
    143     p[1]->SetAngle(atan(tan2));
     144    static TRandom rand(0);
     145    Double_t rnd = rand.Uniform(TMath::Pi() * 2);
     146    p0.SetNewDirection(atan(tan1), rnd);
     147    p1.SetNewDirection(atan(tan2), rnd+TMath::Pi());
    144148
    145     list->Add(p[0]);
    146     list->Add(p[1]);
     149    list->Add(&p0);
     150    list->Add(&p1);
    147151
    148152    return kTRUE;

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

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

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

    r1349 r1352  
    77#include "MParticle.h"
    88
     9#include <TRandom.h>
     10#include <TMatrixD.h>
     11
    912ClassImp(MParticle);
    1013
     14/***********************
     15 *
     16 * calc r from z:
     17 *
     18 *        R = 2*c/H0 *(1+z - sqrt(1+z))
     19 *
     20 *        z = -0.5 * (3 + R' +/- sqrt(R'+1))   R' = R*H0/c
     21 *
     22 *        H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
     23 *
     24 ************************
     25 */
     26
     27Double_t ZofR(Double_t *x, Double_t *k=NULL)
     28{
     29    const Double_t c  = 299792458;        // [m/s]
     30    const Double_t H0 = 50./3.0857e19;    // [km / Mpc s] --> [s^-1]
     31
     32    const Double_t ly = 3600.*24.*365.*c; // [m/ly]
     33    const Double_t pc = 1./3.258;         // [pc/ly]
     34    const Double_t r  = x[0] /pc*ly*1e3;  // [m]
     35
     36    const Double_t R = r*H0/c;            // [1]
     37
     38    return (R+1+sqrt(R*2+1))/2 - 1;
     39}
     40
     41Double_t RofZ(Double_t *x, Double_t *k=NULL)
     42{
     43    Double_t z1 = x[0] + 1;
     44
     45    const Double_t c  = 299792458;                 // [m/s]
     46    const Double_t H0 = 50./3.0857e19;             // [km / Mpc s] --> [s^-1]
     47
     48    const Double_t ly = 3600.*24.*365.*c;          // [m/ly]
     49    const Double_t pc = 1./3.258;                  // [pc/ly]
     50
     51    const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]
     52
     53    return  R * pc/ly/1e3;                   // [kpc]
     54}
     55
    1156MParticle::MParticle(ParticleType_t t, const char *name, const char *title)
    12     : fPType(t)
     57    : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0)
    1358{
    1459    //
    … …  
    2570}
    2671
     72void MParticle::SetNewDirection(Double_t theta, Double_t phi)
     73{
     74    TMatrixD B(3, 3);
     75
     76    B(0, 0) =  cos(fTheta)*cos(fPsi);
     77    B(1, 0) =  cos(fTheta)*sin(fPsi);
     78    B(2, 0) = -sin(fTheta);
     79
     80    B(0, 1) = -sin(fPsi);
     81    B(1, 1) =  cos(fPsi);
     82    B(2, 1) =  0;
     83
     84    B(0, 2) = sin(fTheta)*cos(fPsi);
     85    B(1, 2) = sin(fTheta)*sin(fPsi);
     86    B(2, 2) = cos(fTheta);
     87
     88    // ------------------------------
     89
     90    TVectorD r(3);
     91
     92    r(0) = sin(theta)*cos(phi);
     93    r(1) = sin(theta)*sin(phi);
     94    r(2) = cos(theta);
     95
     96    // ------------------------------
     97
     98    r *= B;
     99
     100    fTheta = sqrt(r(0)*r(0)+r(1)*r(1)); // Numerically bad: acos(r(2));
     101    fPsi   = atan2(r(1), r(0));
     102
     103    if (fTheta*2 > TMath::Pi())
     104        fTheta = fabs(fTheta-TMath::Pi());
     105}
     106
     107#include <iostream.h>
     108#include <iomanip.h>
     109
     110Bool_t MParticle::SetNewPosition(Double_t dr)
     111{
     112    Bool_t rc=kTRUE;
     113
     114    TVectorD x(3);
     115
     116    x(0) = sin(fTheta)*cos(fPsi);
     117    x(1) = sin(fTheta)*sin(fPsi);
     118    x(2) = cos(fTheta);
     119
     120    x *= dr;
     121
     122    // ------------------------------
     123
     124    const Double_t R = RofZ(&fZ);
     125
     126    if (x(2) > R*cos(fTheta))
     127    {
     128        x *= R/dr;
     129        rc = kFALSE;
     130    }
     131
     132    // ------------------------------
     133
     134    TVectorD r(3);
     135
     136    r(0) = fR*cos(fPhi);
     137    r(1) = fR*sin(fPhi);
     138    r(2) = R;
     139
     140    // ------------------------------
     141
     142    r -= x;
     143
     144    fR   = sqrt(r(0)*r(0)+r(1)*r(1));
     145    fPhi = atan2(r(1), r(0));
     146    fZ   = ZofR(&r(2));
     147
     148    return rc;
     149}
     150
     151Bool_t MParticle::SetNewPosition()
     152{
     153    static TRandom rand(0);
     154    Double_t r = rand.Exp(GetInteractionLength());
     155    return SetNewPosition(r);
     156}

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

    r1349 r1352  
    1717    Double_t fEnergy; // [GeV] Energy
    1818    Double_t fBeta;   // [1]   beta
    19     Double_t fAngle;  // [rad] Angle
     19
    2020    Double_t fZ;      // [1]   Red shift
     21    Double_t fR;      // [kpc] Radius from line of view
     22    Double_t fPhi;    // [rad] Phi@z from line of view
     23
     24    Double_t fTheta;  // [rad] Direction of momentum, angle || line of view
     25    Double_t fPsi;    // [rad] Direction of momentum, az. angle
    2126
    2227public:
    … …  
    2429     //    ~MParticle();
    2530
     31    void operator=(MParticle &p)
     32    {
     33        fZ = p.fZ;
     34        fR = p.fR;
     35        fPhi = p.fPhi;
     36        fTheta = p.fTheta;
     37        fPsi = p.fPsi;
     38    }
     39
    2640    //    void Fill(const MParContainer *inp);
    2741
    2842    //    void Draw(Option_t *opt=NULL);
     43    virtual Double_t GetInteractionLength() const = 0;
    2944
    3045    void SetEnergy(Double_t e) { fEnergy = e; }
    31     void SetAngle(Double_t a)  { fAngle = a; }
    3246    void SetBeta(Double_t b)   { fBeta = b; }
     47
    3348    void SetZ(Double_t z)      { fZ = z; }
     49    /*
     50     void SetPhi(Double_t p)    { fPhi = p; }
     51     void SetR(Double_t r)      { fR = r; }
     52     void SetTheta(Double_t t)  { fTheta = t; }
     53     void SetPsi(Double_t p)    { fPsi = p; }
     54     */
     55
     56    void SetNewDirection(Double_t theta, Double_t phi);
     57    Bool_t SetNewPosition(Double_t dr);
     58    Bool_t SetNewPosition();
    3459
    3560    Double_t GetEnergy() const { return fEnergy; }
    36     Double_t GetAngle() const  { return fAngle; }
     61
    3762    Double_t GetZ() const      { return fZ; }
     63    Double_t GetPhi() const    { return fPhi; }
     64    Double_t GetR() const      { return fR; }
     65
     66    Double_t GetTheta() const  { return fTheta; }
     67    Double_t GetPsi() const    { return fPsi; }
    3868
    3969    ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters

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

    r1349 r1352  
    1515    }
    1616
     17    void operator=(MParticle &p) { MParticle::operator=(p); }
     18
    1719    static Double_t Planck(Double_t *x, Double_t *k=NULL);
    1820    static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL);

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

    r1349 r1352  
    66#include <TF1.h>
    77#include <TH1.h>
     8#include <TH2.h>
    89#include <TList.h>
    910#include <TRandom.h>
    … …  
    689690void DoIt()
    690691{
    691     Double_t rcygnusx3 = 100; //30/3.258; // [~10kpc]
     692    Double_t rcygnusx3 = 5000; //100; //30/3.258; // [~10kpc]
    692693
    693694    Double_t startz = ZofR(&rcygnusx3);
    … …  
    697698
    698699    //Double_t nphot = 50;
    699     Double_t runtime = 18*60*60; // [s]
     700    Double_t runtime = 15*60; ///*18*60*/60; // [s]
    700701
    701702    Double_t lo = 2e4;
    … …  
    717718
    718719    gStyle->SetOptStat(10);
     720
     721    TH2D position;
     722    TH2D angle;
     723    TH1D histpos;
     724    TH1D hista;
     725
     726    MBinning binspolx;
     727    MBinning binspoly;
     728    binspolx.SetEdges(32, -180, 180);
     729    binspoly.SetEdgesLog(20, 1e-9, 1e-5);
     730    MH::SetBinning(&position, &binspolx, &binspoly);
     731    MH::SetBinning(&angle,    &binspolx, &binspoly);
     732    MH::SetBinning(&histpos,  &binspoly);
     733    MH::SetBinning(&hista,    &binspoly);
     734
     735    hista.SetTitle("Particle Angle");
     736    angle.SetTitle("Particle Angle");
     737    histpos.SetTitle("Particle Position");
     738    position.SetTitle("Particle Position");
    719739
    720740    TH1D hist;
    … …  
    735755    MH::SetBinning(&histsrc, &bins);
    736756
    737     TH1D hista;
    738     MBinning binsa;
    739     binsa.SetEdgesLog(16, 1e-12, 1e-8);
    740     MH::SetBinning(&hista, &binsa);
    741 
    742757    MH::MakeDefCanvas();
    743758
    … …  
    791806            while ((p=(MPhoton*)NextP()))
    792807            {
    793                 //
    794                 // Calculate way until interaction takes place
    795                 //
    796                 Double_t z = p->GetZ();
    797                 Double_t l = rand.Exp(p->GetInteractionLength());
    798                 Double_t r = RofZ(&z);
    799 
    800                 //
    801                 // If photon went along us...  (l==0: infinite)
    802                 //
    803                 if (r<l || l==0)
     808                if (!p->SetNewPosition())
    804809                {
    805                     cout << (l==0?">":"!") << flush;
     810                    cout << "!" << flush;
     811
    806812                    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
    807818                    fcas << p->GetEnergy() << endl;
     819
    808820                    listg.Remove(p);
    809821                    continue;
    810822                }
    811 
    812                 //
    813                 // Set new z value
    814                 //
    815                 r -= l;
    816                 z = ZofR(&r);
    817                 p->SetZ(z);
    818823
    819824                //
    … …  
    822827                MPhoton photon;
    823828
    824                 phot.SetParameter(0, z);
     829                phot.SetParameter(0, p->GetZ());
    825830                photon.SetEnergy(pow(10, phot.GetRandom()));
    826                 photon.SetAngle(rand.Uniform(TMath::Pi() * 2));
    827 
    828                 if (!pair.Process(p, &photon, &liste))
     831                Double_t theta = rand.Uniform(TMath::Pi() * 2);
     832
     833                if (!pair.Process(p, &photon, theta, &liste))
    829834                {
    830835                    cout << "0" << flush;;
    … …  
    832837                }
    833838
     839                listg.Remove(p);
    834840                cout << "." << flush;
    835 
    836                 listg.Remove(p);
    837841            }
    838842
    … …  
    847851            {
    848852                cout << ":" << flush;
    849 
    850                 hista.Fill(fabs(e->GetAngle()));
    851 
    852853                while(1)
    853854                {
    854                     Double_t E = e->GetEnergy();
    855                     Double_t z = e->GetZ();
    856                     Double_t l = rand.Exp(e->GetInteractionLength());
    857                     Double_t r = RofZ(&z) - l;
    858 
    859                     if (r<0)
     855                    if (!e->SetNewPosition())
    860856                    {
    861857                        cout << "!" << flush;
    … …  
    863859                    }
    864860
    865                     z = ZofR(&r);
    866                     e->SetZ(z);
    867 
    868                     Double_t ep = e->GetEnergyLoss();
    869 
    870                     if (E-ep<lo*5)
     861                    MPhoton *p = e->DoInvCompton();
     862                    if (e->GetEnergy()<lo*5)
    871863                    {
    872864                        cout << "0" << flush;
     865                        delete p;
    873866                        break;
    874867                    }
    875868
    876                     e->SetEnergy(E-ep);
    877 
    878                     if (ep<lo)
     869                    if (p->GetEnergy()<lo)
    879870                    {
    880871                        cout << "i" << flush; // ignored
     872                        delete p;
    881873                        continue;
    882874                    }
    883                     /*
    884                     Double_t eps = 1e-11; // photon vorher
    885 
    886                     Double_t E0 = 511e-6;
    887                     Double_t gamma = E/E0;
    888                     Double_t gamma2 = gamma*gamma;
    889                     Double_t beta = sqrt(1-1/gamma2);
    890 
    891                     Double_t theta = rand.Uniform(TMath::Pi()*2);
    892 
    893                     Double_t p3 = gamma2 * (1-cos(theta)) - ep/E0;
    894 
    895                     Double_t a  = 1- (1 + ep/(eps*p3));
    896 
    897                     cout << "/" << gamma << "," << p3 << "," << ep << "," << a;
    898                     cout << "(" << 90-180*(TMath::Pi()-acos(a))/TMath::Pi() <<")" << flush;
    899                     */
    900 
    901                     MPhoton *p = new MPhoton(ep, z);
     875
    902876                    listg.Add(p);
    903 
    904877                    cout << "." << flush;
    905878                }
    … …  
    912885        if (now!=timer || n<10)
    913886        {
    914             histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
     887            histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz,
     888                                  (int)runtime/60, (int)runtime%60, hist.GetEntries()));
    915889            gPad->Modified();
    916890            gPad->Update();
    … …  
    943917
    944918    MH::MakeDefCanvas();
    945     hista.SetXTitle("[rad]");
     919    position.SetXTitle("Pos: \\Phi [\\circ]");
     920    position.SetYTitle("Pos: R [kpc]");
     921    position.DrawCopy("surf1pol");
     922    gPad->SetLogy();
     923
     924    MH::MakeDefCanvas();
     925    angle.SetXTitle("Angle: \\Psi [\\circ]");
     926    angle.SetYTitle("Angle: \\Theta [\\circ]");
     927    angle.DrawCopy("surf1pol");
     928    gPad->SetLogy();
     929
     930    MH::MakeDefCanvas();
     931    histpos.SetXTitle("R [kpc]");
     932    histpos.SetYTitle("Counts");
     933    histpos.DrawCopy();
     934    gPad->SetLogx();
     935
     936    MH::MakeDefCanvas();
     937    hista.SetXTitle("\\Theta [\\circ]");
     938    hista.SetYTitle("Counts");
    946939    hista.DrawCopy();
    947940    gPad->SetLogx();
Note: See TracChangeset for help on using the changeset viewer.