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

Legend:

Unmodified
Added
Removed

  • trunk/WuerzburgSoft/Thomas/mphys/Changelog

    r1363 r1364  
    11                                                                  -*-*- END -*-*-
     2 2002/06/14: Thomas Bretz
     3
     4   * MElectron.[h,cc]:
     5     - Moved Planck function to MParticle
     6     - Added the DiSum function
     7     - changed Li2 to use the DiSum function (speed reasons)
     8     - changed the eps of all integrals to 1e-2
     9     - changed the p_e function to use the Compton function
     10
     11   * MPhoton.[h,cc]:
     12     - removed the planck function
     13     - changed the eps of all integrals to 1e-2
     14
     15   * MParticle.[h,cc]:
     16     - added the planck function
     17
     18   * Makefile:
     19     - removed unused source files
     20
     21   * analp.C:
     22     - did some small fixes
     23     - changed the sclae of the radius- and phi-plots
     24
     25   * phys.C:
     26     - changed the integral eps to 1e-2
     27     
     28
     29
     30
    231 2002/06/13: Thomas Bretz
    332

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

    r1363 r1364  
    4242ClassImp(MElectron);
    4343
    44 Double_t MElectron::Planck(Double_t *x, Double_t *k=NULL)
    45 {
    46     //
    47     // Planck, per unit volume, per unit energy
    48     //
    49     // constants moved out of function
    50     //
    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;
    59 
    60     /*
    61      //Double_t c   = 299792458;             // [m/s]
    62      //Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
    63      //Double_t hc  = h*c;                   // [GeVm]
    64      Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
    65      return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]
    66      */
    67 
    68     return E*E / (exp(EkT)-1.); // [GeV^2]
    69 }
    70 
    7144Double_t MElectron::Li(Double_t *x, Double_t *k)
    7245{
    7346    const Double_t t = x[0];
    74 
    7547    return log(1.-t)/t;
    7648}
    7749
     50Double_t DiSum(Double_t *x, Double_t *k=NULL)
     51{
     52    Double_t t = x[0];
     53
     54    const Double_t eps = fabs(t*1e-2);
     55
     56    Double_t disum = t;
     57    Double_t add = 0;
     58
     59    Int_t    n   = 2;
     60    Double_t pow = t*t;   // t^2
     61
     62    do
     63    {
     64        add = pow/n/n;
     65
     66        pow *= t;       // pow = t^n
     67        n++;
     68
     69        disum += add;
     70
     71    } while (fabs(add)>eps);
     72
     73    return disum;
     74}
     75
    7876Double_t MElectron::Li2(Double_t *x, Double_t *k=NULL)
    7977{
    8078    //
    8179    // Dilog, Li2
    82     //
    83     Double_t z = x[0];
     80    // ----------
     81    //
     82    //   Integral(0, 1) = konst;
     83    //   Double_t konst = 1./6*TMath::Pi()*TMath::Pi();
     84    //
     85    //   x[0]: z
     86    //
     87    const Double_t z = x[0];
     88
     89    if (fabs(z)<1)
     90        return DiSum(x);
    8491
    8592    TF1 IntLi("Li", Li, 0, z, 0);
    86     const Double_t integ = IntLi.Integral(0, z);
    87 
    88     /*
    89      if (fabs(z)<1)
    90      {
    91      Double_t disum = DiSum(&z);
    92      cout << "Disum (" << z << ") " << disum << "=" << -integ << "\t" << disum-integ << endl;
    93      return disum;
    94      }
    95     */
    96 
    97     /*
    98       Integral(0, 1) = konst;
    99       Double_t konst = 1./6*TMath::Pi()*TMath::Pi();
    100       */
    101 
     93    const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2);
    10294    return -integ;
    10395}
    … …  
    114106    const Double_t d  = w4*(w4 + 1);
    115107
    116     Double_t s   = -w*2*(1+1); // -2*omega*(1+beta)
     108    Double_t s = -w*2*(1+1); // -2*omega*(1+beta)
    117109    const Double_t li2 = Li2(&s);
    118110
    … …  
    125117{
    126118    const Double_t E0  = 511e-6; //[GeV]
    127     const Double_t E02 = E0*E0;
    128119
    129120    Double_t epsilon = x[0];
    130     Double_t E       = k[0];
    131     // Double_t beta    = k[1];
    132     Double_t z       = k[2];
    133 
    134     Double_t omega  = epsilon*E/E02;
    135 
    136     const Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]
     121    Double_t z       = k[1];
     122
     123    const Double_t E = k[0];
     124
     125    Double_t omega  = epsilon*E/(E0*E0);
     126
     127    const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1]
    137128    return Flim(&omega)*n;
    138129}
    139 
    140130
    141131Double_t MElectron::InteractionLength(Double_t *E, Double_t *k=NULL)
    … …  
    170160    const Double_t z      = k ? k[0] : 0;
    171161
    172     // Double_t beta = sqrt(1-E0/E*E0/E);
    173     const Double_t beta   = 1; //0.999999999999999999999999999;
    174 
    175     Double_t val[3] = { E[0], beta, z };             // E[GeV]
     162    Double_t val[3] = { E[0], z };                         // E[GeV]
    176163
    177164    const Double_t from   = 1e-17;
    … …  
    184171       -----------------------------------
    185172    */
    186     TF1 func("Compton", Compton, from, to, 3);      // [0, inf]
    187 
    188     const Double_t integ  = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]
     173    TF1 func("Compton", Compton, from, to, 2);            // [0, inf]
     174
     175    Double_t integ  = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf]
    189176
    190177    const Double_t aE     = alpha/E[0];                   // [1/GeV]
     178
     179    const Double_t beta   = 1;
    191180
    192181    const Double_t konst  = 2.*E02/hc/beta;               // [1 / GeV m]
    … …  
    196185    const Double_t pc     = 1./3.258;                     // [pc/ly]
    197186
    198     return (1./ret)/ly*pc/1000;                     // [kpc]
     187    return (1./ret)/ly*pc/1000;                           // [kpc]
    199188}
    200189
    … …  
    211200// --------------------------------------------------------------------------
    212201
    213 Double_t MElectron::p_e(Double_t *x, Double_t *k)
    214 {
    215     Double_t e  = pow(10, x[0]);
    216     Double_t z  = k[1];
    217 
    218     const Double_t E  = k[0];
    219 
    220     const Double_t E0  = 511e-6; //[GeV]
    221     const Double_t E02 = E0*E0;
    222 
    223     Double_t omega = e*E/E02;
    224 
    225     const Double_t n = Planck(&e, &z);
    226 
    227     const Double_t F = Flim(&omega)/omega/omega;
    228 
    229     return n*F*1e26;
     202inline Double_t MElectron::p_e(Double_t *x, Double_t *k)
     203{
     204    Double_t e = pow(10, x[0]);
     205    return Compton(&e, k);
     206    /*
     207     Double_t z  = k[1];
     208
     209     const Double_t E  = k[0];
     210
     211     const Double_t E0  = 511e-6; //[GeV]
     212     const Double_t E02 = E0*E0;
     213
     214     Double_t omega = e*E/E02;
     215
     216     const Double_t n = Planck(&e, &z);
     217
     218     const Double_t F = Flim(&omega)/omega/omega;
     219
     220     return n*F*1e26;
     221     */
    230222}
    231223

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

    r1356 r1364  
    2121    // ----------------------------------------------------------------
    2222
    23     static Double_t Planck(Double_t *x, Double_t *k=NULL);
    2423    static Double_t Li(Double_t *x, Double_t *k);
    2524    static Double_t Li2(Double_t *x, Double_t *k=NULL);

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

    r1359 r1364  
    7171     Double_t z1 = x[0] + 1;
    7272
    73      const Double_t c  = 299792458;                 // [m/s]
    74      const Double_t H0 = 50./3.0857e19;             // [km / Mpc s] --> [s^-1]
    75 
    76      const Double_t ly = 3600.*24.*365.*c;          // [m/ly]
    77      const Double_t pc = 1./3.258;                  // [pc/ly]
     73     const Double_t c  = 299792458;                   // [m/s]
     74     const Double_t H0 = 50./3.0857e19;               // [km / Mpc s] --> [s^-1]
     75
     76     const Double_t ly = 3600.*24.*365.*c;            // [m/ly]
     77     const Double_t pc = 1./3.258;                    // [pc/ly]
    7878
    7979     const Double_t R = c/H0 * 2 * (1 - 1./sqrt(z1)); // [m]
    8080
    81      return  R * pc/ly/1e3;                   // [kpc]
     81     return  R * pc/ly/1e3;                           // [kpc]
     82}
     83
     84Double_t MParticle::Planck(Double_t *x, Double_t *k)
     85{
     86    //
     87    // Planck, per unit volume, per unit energy
     88    //
     89    // constants (see below) moved out of function
     90    //
     91    const Double_t E   = x[0];                    // [GeV]
     92    const Double_t z   = k ? k[0] : 0;
     93
     94    const Double_t T   = 2.96*(z+1);              // [K]
     95    const Double_t e   = 1.602176462e-19;         // [C]
     96    const Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
     97
     98    const Double_t EkT = E/kB/T;
     99
     100    /*
     101     //Double_t c   = 299792458;             // [m/s]
     102     //Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
     103     //Double_t hc  = h*c;                   // [GeVm]
     104     Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
     105     return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]
     106     */
     107
     108    return E*E / (exp(EkT)-1.); // [GeV^2]
     109
    82110}
    83111

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

    r1362 r1364  
    4646    }
    4747
     48    // ----------------------------------------------------------------
     49
     50    static Double_t Planck(Double_t *x, Double_t *k=NULL);
     51
     52    // ----------------------------------------------------------------
     53
    4854    void SetIsPrimary(Bool_t is=kTRUE) { is ? SetBit(kEIsPrimary) : ResetBit(kEIsPrimary); }
    4955    Bool_t IsPrimary() const { return TestBit(kEIsPrimary); }

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

    r1363 r1364  
    3838
    3939ClassImp(MPhoton);
    40 
    41 Double_t MPhoton::Planck(Double_t *x, Double_t *k=NULL)
    42 {
    43     //
    44     // Planck, per unit volume, per unit energy
    45     //
    46     // constants moved out of function, see below
    47     //
    48     const Double_t E   = x[0];                    // [GeV]
    49     const Double_t z   = k ? k[0] : 0;
    50 
    51     const Double_t T   = 2.96*(z+1);              // [K]
    52     const Double_t e   = 1.602176462e-19;         // [C]
    53     const Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
    54 
    55     const Double_t EkT = E/kB/T;
    56 
    57     /*
    58      //Double_t c   = 299792458;             // [m/s]
    59      //Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
    60      //Double_t hc  = h*c;                   // [GeVm]
    61 
    62      Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
    63      return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]
    64      */
    65 
    66     return E*E / (exp(EkT)-1.); // [GeV^2]
    67 }
    6840
    6941Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL)
    … …  
    12597    TF1 f("int1", Int1, from, to, 2);
    12698
    127     const Double_t int1   = f.Integral(from, to, val);  // [m^2]
    128     const Double_t planck = Planck(&Ep, &z);            // [GeV^2]
     99    const Double_t int1   = f.Integral(from, to, val, 1e-2);  // [m^2]
     100    const Double_t planck = MParticle::Planck(&Ep, &z); // [GeV^2]
    129101
    130102    const Double_t res = planck * int1;
    … …  
    157129    TF1 f("int2", Int2, lolim, inf, 2);
    158130
    159     Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2]
     131    Double_t int2 = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2]
    160132
    161133    if (int2==0)
    … …  
    166138
    167139    /* Planck constants: konst */
    168     Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
     140    const Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
    169141    int2 *= konst;           // [1 / m]
    170142

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

    r1363 r1364  
    2121    // ----------------------------------------------------------------
    2222
    23     static Double_t Planck(Double_t *x, Double_t *k=NULL);
    2423    static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL);
    2524    static Double_t Int1(Double_t *x, Double_t *k=NULL);

  • trunk/WuerzburgSoft/Thomas/mphys/Makefile

    r1349 r1364  
    3131           MPhoton.cc \
    3232           MElectron.cc \
    33            MGenIRPhoton.cc \
    34            MPairProduction.cc \
    35            MGenPrimaryParticle.cc
     33           MPairProduction.cc
    3634
    3735SRCS    = $(SRCFILES)

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

    r1349 r1364  
    88#pragma link C++ class MElectron+;
    99#pragma link C++ class MPhoton+;
    10 #pragma link C++ class MGenPrimaryParticle+;
    11 #pragma link C++ class MGenIRPhoton+;
    1210#pragma link C++ class MPairProduction+;
    1311

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

    r1363 r1364  
    44{
    55    TChain chain("Electrons");
     6    //chain.Add("cascade.root");
     7
    68    chain.Add("cascade_100kpc_01.root");
    79    chain.Add("cascade_100kpc_02.root");
    810    chain.Add("cascade_100kpc_03.root");
     11
    912
    1013    MElectron *e = new MElectron;

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

    r1363 r1364  
    44{
    55    TChain chain("Photons");
    6     chain.Add("cascade_100kpc_01.root");
    7     chain.Add("cascade_100kpc_02.root");
    8     chain.Add("cascade_100kpc_03.root");
     6    chain.Add("cascade_500kpc_0*.root");
     7    chain.Add("cascade_500kpc_2*.root");
     8//    chain.Add("cascade_100kpc_0*.root");
     9//    chain.Add("cascade_100kpc_1*.root");
    910
    1011    MPhoton *p = new MPhoton;
    … …  
    1516    cout << "Found " << n << " entries." << endl;
    1617
     18    if (n==0)
     19        return;
     20
    1721    MBinning binsx;
    18     binsx.SetEdgesLog(21, 1e4, 1e11);
     22    binsx.SetEdgesLog(14, 1e4, 1e11);
    1923
    2024    MBinning binspolx;
    … …  
    2226    MBinning binspoly2;
    2327    binspolx.SetEdges(16, -180, 180);
    24     binspoly1.SetEdges(20, 0, 5e-9);
    25     binspoly2.SetEdges(20, 0, 2e-7);
     28
     29    const Double_t ls = 299792458;        // [m/s]
     30    const Double_t ly = 3600.*24.*365.*ls; // [m/ly]
     31    const Double_t pc = 1./3.258;         // [pc/ly]
     32
     33    binspoly1.SetEdges(20, 0, 2e-6*3600);
     34    binspoly2.SetEdges(20, 0, 1e-5*1e3/pc*365);
    2635
    2736    TH1D h;
    … …  
    5968
    6069        h.Fill(Ep, Ep*Ep * weight);
    61         r.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight);
    62         a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
     70        r.Fill(p->GetPhi()*kRad2Deg, p->GetR()*1e3/pc*365, weight);
     71        a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg*3600, weight);
    6372    }
    6473
    … …  
    6776    gStyle->SetOptStat(10);
    6877
    69     TLine line;
    70     line.SetLineColor(kBlack);
    71     line.SetLineWidth(1);
    72 
    73     TCanvas *c = new TCanvas("c1", "name");
     78//    delete gROOT->FindObject("Analysis Arrival");
     79    c = MH::MakeDefCanvas("Analysis Arrival", "");
    7480    c->Divide(2,2);
    7581
    7682    c->cd(1);
    77     gPad->SetLogx();
    7883    gPad->SetLogy();
     84    gPad->SetLogz();
     85    gPad->SetTheta(70);
     86    gPad->SetPhi(90);
     87    r.SetTitle(" Distance from Observer ");
     88    r.GetXaxis()->SetLabelOffset(-0.015);
     89    r.GetXaxis()->SetTitleOffset(1.1);
     90    r.GetXaxis()->SetRangeUser(1e4, 1e9);
     91    r.SetXTitle("\\Phi [\\circ]");
     92    r.SetYTitle("R [light days]");
     93    r.DrawCopy("surf1pol");
    7994
     95    c->cd(2);
     96    gPad->SetLogy();
     97    gPad->SetLogz();
     98    gPad->SetTheta(70);
     99    gPad->SetPhi(90);
     100    a.SetTitle(" Hit Angle for Observer ");
     101    a.GetXaxis()->SetLabelOffset(-0.015);
     102    a.GetXaxis()->SetTitleOffset(1.1);
     103    a.GetXaxis()->SetRangeUser(1e4, 1e9);
     104    a.SetXTitle("\\Psi [\\circ]");
     105    a.SetYTitle("\\Theta [\"]");
     106    a.DrawCopy("surf1pol");
     107
     108    c->cd(3);
     109    gPad->SetLogy();
     110    TH1 *g=r.ProjectionY("p1");
     111    g->SetBit(kCanDelete);
     112    g->SetTitle(" Hit Observers Plain ");
     113    g->GetXaxis()->SetLabelOffset(0);
     114    g->GetXaxis()->SetTitleOffset(1.1);
     115    g->GetYaxis()->SetTitleOffset(1.3);
     116    g->SetXTitle("R [light days]");
     117    g->SetYTitle("Counts");
     118    g->Draw();
     119
     120    c->cd(4);
     121    gPad->SetLogy();
     122    g=a.ProjectionY("p2");
     123    g->SetBit(kCanDelete);
     124    g->SetTitle("Hit Angle");
     125    g->GetXaxis()->SetLabelOffset(0);
     126    g->GetXaxis()->SetTitleOffset(1.1);
     127    g->GetYaxis()->SetTitleOffset(1.3);
     128    g->SetXTitle("\\Phi [\"]");
     129    g->SetYTitle("Counts");
     130    g->Draw();
     131
     132  //  delete gROOT->FindObject("Analysis Photons");
     133    TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
     134    c->Divide(1,2);
     135
     136    c->cd(1);
    80137    gPad->SetLogx();
    81138    gPad->SetLogy();
    … …  
    84141    h.GetXaxis()->SetLabelOffset(-0.015);
    85142    h.GetXaxis()->SetTitleOffset(1.1);
    86     h.GetXaxis()->SetRangeUser(1e4, 1e9);
     143//    h.GetXaxis()->SetRangeUser(1e4, 1e9);
    87144    prim.SetMarkerStyle(kPlus);
    88145    h.SetMarkerStyle(kMultiply);
    … …  
    96153 
    97154    c->cd(2);
    98     r.SetTitle(" Distance from Observer ");
    99     r.GetXaxis()->SetLabelOffset(-0.015);
    100     r.GetXaxis()->SetTitleOffset(1.1);
    101     r.GetXaxis()->SetRangeUser(1e4, 1e9);
    102     r.SetXTitle("\\Phi [\\circ]");
    103     r.SetYTitle("R [kpc]");
    104     r.DrawCopy("surf1pol");
    105 
    106     c->cd(3);
    107     a.SetTitle(" Hit Angle for Observer ");
    108     a.GetXaxis()->SetLabelOffset(-0.015);
    109     a.GetXaxis()->SetTitleOffset(1.1);
    110     a.GetXaxis()->SetRangeUser(1e4, 1e9);
    111     a.SetXTitle("\\Psi [\\circ]");
    112     a.SetYTitle("\\Theta [\\circ]");
    113     a.DrawCopy("surf1pol");
    114 
    115     c->cd(4);
    116     /*
    117      TH1 *g=a.ProjectionY();
    118      g->SetBit(kCanDelete);
    119      g->SetTitle("Hit Angle");
    120      g->GetXaxis()->SetLabelOffset(0);
    121      g->GetXaxis()->SetTitleOffset(1.1);
    122      g->GetYaxis()->SetTitleOffset(1.3);
    123      g->SetXTitle("\\Phi [\\circ]");
    124      g->SetYTitle("Counts");
    125      g->Draw();
    126      */
    127155    gPad->SetLogx();
    128156    TH1D div;

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

    r1363 r1364  
    3535    Double_t res = MPhoton::Int2(&Ep, k);
    3636
    37     return res*1e55; //65/k[0];
    38     // return MPhoton::Planck(&Ep, &k[1]);
     37    return res*1e55;
    3938}
    4039
    … …  
    125124    MPairProduction pair;
    126125
    127     Double_t runtime = 15*60; // [s]
     126    Double_t runtime = 7*60*60; // [s]
    128127
    129128    Double_t lo = 1e4;
    … …  
    156155    MBinning binspoly2;
    157156    binspolx.SetEdges(16, -180, 180);
    158     binspoly1.SetEdges(20, 0, 5e-9);
    159     binspoly2.SetEdges(20, 0, 2e-7);
     157    binspoly1.SetEdges(20, 0, 2e-6);
     158    binspoly2.SetEdges(20, 0, 1e-5);
    160159    MH::SetBinning(&angle,    &binspolx, &binspoly1);
    161160    MH::SetBinning(&position, &binspolx, &binspoly2);
    … …  
    175174    histsrc.SetFillStyle(0);
    176175    histsrc.SetMarkerStyle(kMultiply);
     176    histsrc.SetMarkerColor(kRed);
     177    histsrc.SetLineColor(kRed);
    177178
    178179    MBinning bins;
    … …  
    272273                phot.SetParameter(0, Eg);
    273274                phot.SetParameter(1, p->GetZ());
    274                 if (phot.Integral(-log10(Eg)-8, -10.5)==0)
     275                if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0)
    275276                {
    276277                    cout << "z" << flush;
Note: See TracChangeset for help on using the changeset viewer.