Ignore:
Timestamp:
07/29/02 09:10:23 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/WuerzburgSoft/Thomas/mphys
Files:
4 added
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/WuerzburgSoft/Thomas/mphys/Changelog

    r1448 r1449  
    11                                                                  -*-*- END -*-*-
     2 2002/06/29: Thomas Bretz
     3
     4   * MElectron.[h,cc]:
     5     - if the photon energy is to small return a straight line in
     6       Compton
     7     - changed lolimit of p_e integration -0.5
     8     - changed number of calculated points to 50 for fP, fQ
     9     - changed usage of TRandom to gRandom
     10     - changed some TVector and TMatrix objects to static
     11     - added new copy constructors
     12
     13   * MPairProduction.cc:
     14     - use the new copy constructor of MElectron
     15
     16   * MParticle.[h,cc]:
     17     - added new copy constructors
     18     - Added a TAxis object in Planck
     19     - replaced usage of TRandom by gRandom
     20     - replaced TVector, TMatrix by TVector3, TRotation
     21     - fixed a bug in SetNewPosition (fX was not incremented if the
     22       particle reached the observers plain)
     23     - incremented version-number
     24
     25   * MPhoton.[h,cc]:
     26     - added Fill member function
     27     - added new copy constructors
     28
     29   * MFit.[h,cc], MCascade.[h,cc]:
     30     - added
     31   
     32   * phys.C:
     33     - moved code to MCascade
     34
     35   * Makefile, PhysLinkDef.h:
     36     - added MFit, MCascade
     37
     38   * analp.C:
     39     - simplified
     40     - added fit for cutoff
     41
     42
     43
    244 2002/06/26: Thomas Bretz
    345

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

    r1448 r1449  
    129129    const Double_t E = k[0];
    130130
    131     Double_t omega  = epsilon*E/(E0*E0);
     131    Double_t flim;
     132    if (epsilon<1e-14)
     133    {
     134        const Double_t d = E/(E0*E0);
     135
     136        Double_t omega1 = 1e-13*d;
     137        Double_t omega2 = 1e-12*d;
     138
     139        const Double_t f1 = Flim(&omega1);
     140        const Double_t f2 = Flim(&omega2);
     141
     142        const Double_t m = log10(f2/f1);
     143        const Double_t t = pow(f2, 13)/pow(f1, 12);
     144
     145        flim = pow(epsilon, m) * t;
     146    }
     147    else
     148    {
     149        Double_t omega  = epsilon*E/(E0*E0);
     150        flim = Flim(&omega);
     151    }
    132152
    133153    const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1]
    134     return Flim(&omega)*n;
     154    return flim*n;
    135155}
    136156
    … …  
    250270    const Double_t E0 = 511e-6; //[GeV]
    251271
    252     const Double_t lolim = -log10(E)/7*4-13;
    253 
    254     TF1 fP("p_e", p_e, lolim, -10.6, 2);
    255     TF1 fQ("G",   G_q, 0, 1., 1);
    256 
     272    const Double_t lolim = -log10(E)/7*4-13.5;
     273
     274    static TF1 fP("p_e", p_e, lolim, -10.6, 2);
     275    static TF1 fQ("G",   G_q, 0, 1., 1);
     276
     277    fP.SetNpx(50);
     278    fQ.SetNpx(50);
     279
     280    fP.SetRange(lolim, -10.6);
    257281    fP.SetParameter(0, E);
    258282    fP.SetParameter(1, z);
    … …  
    288312MPhoton *MElectron::DoInvCompton(Double_t theta)
    289313{
    290     static TRandom rand(0);
    291 
    292314    const Double_t E0 = 511e-6; //[GeV]
    293315
    … …  
    304326    const Double_t f = fEnergy/e;
    305327
    306     Double_t t=-1;
     328    Double_t t=-10;
    307329    Double_t arg;
    308330    do
    309331    {
    310         if (t>0)
     332        if (t>-10)
    311333            cout << "~" << flush;
    312         t = rand.Uniform(TMath::Pi())+TMath::Pi()/2;
    313         Double_t er  = epsilon*(gamma-gammabeta*cos(t)); // photon energy rest frame
    314         arg = (f - E0/er - 1)/(sqrt(fEnergy*fEnergy-E0*E0)/e+1);
     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         */
    315367
    316368    } while (arg<-1 || arg>1);
    317369
    318370    const Double_t theta1s = acos(arg);
    319     const Double_t thetas = atan(sin(t)/(gamma*cos(t)-gammabeta));
     371    const Double_t thetas  = atan(sin(t)/(gamma*cos(t)-gammabeta));
    320372
    321373    const Double_t thetastar = thetas-theta1s;
    … …  
    325377    fEnergy -= e;
    326378
    327     const Double_t phi = rand.Uniform(TMath::Pi()*2);
    328 
     379    const Double_t phi = gRandom->Uniform(TMath::Pi()*2);
     380
     381    /*
    329382    MPhoton &p = *new MPhoton(e, fZ);
    330383    p = *this;
     384    */
     385    MPhoton &p = *new MPhoton(*this, e);
    331386    p.SetNewDirection(theta1, phi);
    332387
    … …  
    374429        return SetNewPosition();
    375430
    376     static TRandom rand(0);
    377     Double_t x = rand.Exp(GetInteractionLength());
     431    Double_t x = gRandom->Exp(GetInteractionLength());
    378432
    379433    // -----------------------------
    … …  
    381435    const Double_t E0 = 511e-6;            //[GeV]
    382436
    383     Double_t B_theta = rand.Uniform(TMath::Pi());
    384     Double_t B_phi   = rand.Uniform(TMath::Pi()*2);
    385 
    386     TMatrixD M(3,3);
     437    Double_t B_theta = gRandom->Uniform(TMath::Pi());
     438    Double_t B_phi   = gRandom->Uniform(TMath::Pi()*2);
     439
     440    static TMatrixD M(3,3);
    387441
    388442    M(0, 0) =  sin(B_theta)*cos(B_phi);
    … …  
    398452    M(2, 2) =  0;
    399453
    400     const Double_t beta = sqrt(1-E0/fEnergy*E0/fEnergy);
     454    const Double_t beta = sqrt(1.-E0/fEnergy*E0/fEnergy);
    401455
    402456    //
    … …  
    404458    // which the x-axis is identical with the B-Field
    405459    //
    406     TVectorD v(3);
     460    static TVectorD v(3);
    407461    v(0) = beta*sin(fTheta)*cos(fPsi);
    408462    v(1) = beta*sin(fTheta)*sin(fPsi);
    … …  
    418472    // -----------------------------
    419473
    420     TMatrixD N(3,3);
     474    static TMatrixD N(3,3);
    421475    N(0, 0) =  1;
    422476    N(1, 0) =  0;
    … …  
    436490                                    // v(2) = 0
    437491    // -----------------------------
    438     TVectorD p(3);
     492    static TVectorD p(3);
    439493
    440494    Double_t rho = 0;   
    … …  
    469523    }
    470524
    471     TVectorD s(3);
     525    static TVectorD s(3);
    472526    s(0) = fR*cos(fPhi);
    473527    s(1) = fR*sin(fPhi);
    … …  
    494548    // -----------------------------
    495549
    496     TVectorD w(3);
     550    static TVectorD w(3);
    497551    w(0) = beta_p/beta;
    498552    w(1) = beta_o/beta*cos(rho);

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

    r1430 r1449  
    1616        fZ      = z;
    1717    }
     18
     19    MElectron(MParticle &p, Bool_t type=kFALSE) : MParticle(p, type?MParticle::kEPositron:MParticle::kEElectron) { }
     20    MElectron(MParticle &p, Double_t e, Bool_t type=kFALSE) : MParticle(p, e, type?MParticle::kEPositron:MParticle::kEElectron) {}
    1821
    1922    void operator=(MParticle &p) { MParticle::operator=(p); }

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

    r1370 r1449  
    119119    // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush;
    120120
    121     MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ(), kTRUE);
    122     MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ(), kFALSE);
    123     p0 = *gamma; // Set Position and direction
    124     p1 = *gamma; // Set Position and direction
     121    MElectron &p0 = *new MElectron(*gamma, E*(1.-f), kTRUE);
     122    MElectron &p1 = *new MElectron(*gamma, E*(1.+f), kFALSE);
    125123
    126     static TRandom rand(0);
    127     Double_t rnd = rand.Uniform(TMath::Pi() * 2);
     124    Double_t rnd = gRandom->Uniform(TMath::Pi() * 2);
    128125    p0.SetNewDirection(atan(+tan1), rnd);
    129126    p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi());

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

    r1448 r1449  
    99#include <TRandom.h>
    1010#include <TMatrixD.h>
     11#include <TRotation.h>
    1112
    1213ClassImp(MParticle);
    … …  
    3940    const Double_t r  = x[0] /pc*ly*1e3;  // [m]
    4041
    41     const Double_t R = 1./(1-r*H0/c/2);   // [1]
     42    const Double_t R = 1./(1.-r*H0/c/2);   // [1]
    4243
    4344    return R*R - 1;
    … …  
    6768     const Double_t pc = 1./3.258;                    // [pc/ly]
    6869
    69      const Double_t R = c/H0 * 2 * (1 - 1./sqrt(z1)); // [m]
     70     const Double_t R = c/H0 * 2 * (1. - 1./sqrt(z1)); // [m]
    7071
    7172     return R * pc/ly/1e3;                            // [kpc]
    … …  
    8586    if (!isloaded)
    8687    {
    87         Double_t c   = 299792458;             // [m/s]
    88         Double_t e   = 1.602176462e-19;       // [C]
    89         Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
    90         Double_t hc  = h*c;                   // [GeVm]
     88        Double_t c     = 299792458;             // [m/s]
     89        Double_t e     = 1.602176462e-19;       // [C]
     90        Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
     91        Double_t hc    = h*c;                   // [GeVm]
    9192        Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
    9293
    … …  
    123124    }
    124125
     126    static TAxis &axez = *hist2->GetXaxis();
     127    static TAxis &axee = *hist2->GetYaxis();
     128
    125129    //
    126130    // y = (y1-y0)/(x1-x0) * (x-x0) + y0
    … …  
    129133    Double_t E = x[0];
    130134
    131     Int_t i = hist2->GetXaxis()->FindFixBin(z);
    132     Int_t j = hist2->GetYaxis()->FindFixBin(E);
    133 
    134     Double_t z1 = hist2->GetXaxis()->GetBinLowEdge(i+1);
    135     Double_t z0 = hist2->GetXaxis()->GetBinLowEdge(i);
    136 
    137     Double_t E1 = hist2->GetYaxis()->GetBinLowEdge(j+1);
    138     Double_t E0 = hist2->GetYaxis()->GetBinLowEdge(j);
     135    Int_t i = axez.FindFixBin(z);
     136    Int_t j = axee.FindFixBin(E);
     137
     138    Double_t z1 = axez.GetBinLowEdge(i+1);
     139    Double_t z0 = axez.GetBinLowEdge(i);
     140
     141    Double_t E1 = axee.GetBinLowEdge(j+1);
     142    Double_t E0 = axee.GetBinLowEdge(j);
    139143
    140144    Double_t n00 = hist2->GetBinContent(i,   j);
    … …  
    207211void MParticle::InitRandom()
    208212{
    209     static TRandom rnd(0);
    210     fPhi = rnd.Uniform(TMath::Pi()*2);
    211     fPsi = rnd.Uniform(TMath::Pi()*2);
     213    fPhi = gRandom->Uniform(TMath::Pi()*2);
     214    fPsi = gRandom->Uniform(TMath::Pi()*2);
    212215}
    213216
    214217void MParticle::SetNewDirection(Double_t theta, Double_t phi)
    215218{
    216     TMatrixD B(3, 3);
    217 
    218     B(0, 0) =  cos(fTheta)*cos(fPsi);
    219     B(1, 0) =  cos(fTheta)*sin(fPsi);
    220     B(2, 0) = -sin(fTheta);
    221 
    222     B(0, 1) = -sin(fPsi);
    223     B(1, 1) =  cos(fPsi);
     219    static TMatrixD B(3, 3);
     220
     221    const Double_t ct = cos(fTheta);
     222    const Double_t st = sin(fTheta);
     223    const Double_t cp = cos(fPsi);
     224    const Double_t sp = sin(fPsi);
     225
     226    /*
     227    TRotation B( ct*cp, ct*sp, -st,
     228                -sp,    cp,     0,
     229                 st*cp, st*sp,  ct);
     230                 */
     231
     232    // first row
     233    B(0, 0) =  ct*cp;
     234    B(1, 0) =  ct*sp;
     235    B(2, 0) = -st;
     236
     237    // second row
     238    B(0, 1) = -sp;
     239    B(1, 1) =  cp;
    224240    B(2, 1) =  0;
    225241
    226     B(0, 2) = sin(fTheta)*cos(fPsi);
    227     B(1, 2) = sin(fTheta)*sin(fPsi);
    228     B(2, 2) = cos(fTheta);
    229 
    230     // ------------------------------
    231 
    232     TVectorD r(3);
    233 
    234     r(0) = sin(theta)*cos(phi);
    235     r(1) = sin(theta)*sin(phi);
     242    // third row
     243    B(0, 2) = st*cp;
     244    B(1, 2) = st*sp;
     245    B(2, 2) = ct;
     246
     247    // ------------------------------
     248
     249    static TVectorD r(3);
     250
     251    const Double_t sint = sin(theta);
     252
     253    /*
     254     TVector3 r(sint*cos(phi), sint*sin(phi), cos(theta));
     255     */
     256
     257    r(0) = sint*cos(phi);
     258    r(1) = sint*sin(phi);
    236259    r(2) = cos(theta);
    237260
    … …  
    253276    Bool_t rc=kTRUE;
    254277
    255     TVectorD x(3);
    256 
    257     x(0) = sin(fTheta)*cos(fPsi);
    258     x(1) = sin(fTheta)*sin(fPsi);
    259     x(2) = cos(fTheta);
    260 
    261     // ------------------------------
    262 
    263     const Double_t R = RofZ(&fZ);
    264 
    265     const Double_t dx = R - dr*x(2);
    266 
     278    const Double_t st = sin(fTheta);
     279
     280    TVector3 d(st*cos(fPsi), st*sin(fPsi), cos(fTheta));
     281
     282    // ------------------------------
     283
     284    const Double_t R  = RofZ(&fZ);
     285    const Double_t dx = R - dr*d.z();
    267286    if (dx < 0)
    268287    {
    269         dr = R/x(2); // R>0 --> x(2)>0
     288        dr = R/d.z(); // R>0 --> x(2)>0
    270289        rc = kFALSE;
    271290    }
    272     else
    273         fX += dr*(1.-x(2));
    274 
    275     x  *= dr;
    276 
    277     // ------------------------------
    278 
    279     TVectorD r(3);
    280 
    281     r(0) = fR*cos(fPhi);
    282     r(1) = fR*sin(fPhi);
    283     r(2) = R;
    284 
    285     // ------------------------------
    286 
    287     r -= x;
     291
     292    fX += dr*(1.-d(2));
     293    d  *= dr;
     294
     295    // ------------------------------
     296
     297    TVector3 r(fR*cos(fPhi), fR*sin(fPhi), R);
     298
     299    r -= d;
     300
     301    // ------------------------------
    288302
    289303    if (fR!=0)
    290         fPhi = atan2(r(1), r(0));
    291 
    292     fR = sqrt(r(0)*r(0)+r(1)*r(1));
     304        fPhi = atan2(r.y(), r.x());
     305    fR = sqrt(r.x()*r.x()+r.y()*r.y());
    293306    fZ = ZofR(&r(2));
    294307
    … …  
    298311Bool_t MParticle::SetNewPosition()
    299312{
    300     static TRandom rand(0);
    301     Double_t r = rand.Exp(GetInteractionLength());
     313    Double_t r = gRandom->Exp(GetInteractionLength());
    302314    return SetNewPosition(r);
    303315}

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

    r1448 r1449  
    3737public:
    3838    MParticle(ParticleType_t t=kENone, const char *name=NULL, const char *title=NULL);
     39    MParticle(MParticle &p, ParticleType_t t=kENone) : fPType(t) { *this = p; }
     40    MParticle(MParticle &p, Double_t e, ParticleType_t t=kENone) : fPType(t) { *this = p; fEnergy = e; }
    3941
    4042    void InitRandom();
    … …  
    8991    Double_t GetPsi() const    { return fPsi; }
    9092
    91     ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters
     93    ClassDef(MParticle, 2) // Container which holds hostograms for the Hillas parameters
    9294};
    9395

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

    r1448 r1449  
    200200    DrawInteractionLength(fZ);
    201201}
     202
     203void MPhoton::Fill(TH1 &h, Double_t idx, Double_t w) const
     204{
     205    h.Fill(fEnergy, pow(fEnergy, idx)*w);
     206}

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

    r1428 r1449  
    55#include "MParticle.h"
    66#endif
     7
     8class TH1;
    79
    810class MPhoton : public MParticle
    … …  
    1719    }
    1820
     21        MPhoton(MParticle &p) : MParticle(p, MParticle::kEGamma) { }
     22        MPhoton(MParticle &p, Double_t e) : MParticle(p, e, MParticle::kEGamma) { }
     23
    1924    void operator=(MParticle &p) { MParticle::operator=(p); }
     25
     26    void Fill(TH1 &h, Double_t idx, Double_t w) const;
    2027
    2128    // ----------------------------------------------------------------

  • trunk/WuerzburgSoft/Thomas/mphys/Makefile

    r1364 r1449  
    2222#  connect the include files defined in the config.mk file
    2323#
    24 INCLUDES = -I. -I../mbase
     24INCLUDES = -I. -I../mbase -I../mhist
    2525
    2626#------------------------------------------------------------------------------
    … …  
    2929
    3030SRCFILES = MParticle.cc \
     31           MPairProduction.cc \
    3132           MPhoton.cc \
    32            MElectron.cc \
    33            MPairProduction.cc
     33           MFit.cc \
     34           MCascade.cc \
     35           MElectron.cc
    3436
    3537SRCS    = $(SRCFILES)

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

    r1364 r1449  
    88#pragma link C++ class MElectron+;
    99#pragma link C++ class MPhoton+;
     10#pragma link C++ class MCascade+;
     11#pragma link C++ class MFit+;
    1012#pragma link C++ class MPairProduction+;
    1113

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

    r1428 r1449  
    22
    33#include <float.h>
    4 
    5 Double_t GetMin(TH1 *h)
    6 {
    7     Double_t min = FLT_MAX;
    8     for (int i=1; i<=h->GetXaxis()->GetNbins(); i++)
    9     {
    10         if (h->GetBinContent(i)<min && h->GetBinContent(i)!=0)
    11             min = h->GetBinContent(i);
    12     }
    13     return min;
    14 }
    15 
    16 void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1)
     4#include <iomanip.h>
     5
     6#include <TH1.h>
     7#include <TH2.h>
     8#include <TLine.h>
     9#include <TLeaf.h>
     10#include <TChain.h>
     11#include <TGaxis.h>
     12#include <TStyle.h>
     13#include <TCanvas.h>
     14#include <TPolyLine.h>
     15#include <TPolyMarker.h>
     16
     17#include "MH.h"
     18#include "MBinning.h"
     19
     20#include "MFit.h"
     21#include "MPhoton.h"
     22
     23void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE)
    1724{
    1825    TString str("MPhoton.MParticle.");
     26    str += name;
    1927
    2028    MPhoton *p = new MPhoton;
    2129
    2230    chain->SetBranchAddress("MPhoton.", &p);
     31    chain->SetBranchStatus("*", 0);
     32    chain->SetBranchStatus(str, 1);
    2333
    2434    min =  FLT_MAX;
    2535    max = -FLT_MAX;
    2636
    27     Int_t n = chain->GetEntries();
     37    Int_t  numt = chain->GetTreeNumber();
     38    TLeaf *leaf = chain->GetLeaf(str);
     39
     40    if (!leaf && numt>=0)
     41        return;
     42
     43    Stat_t n = chain->GetEntries();
    2844    for (int i=0; i<n; i++)
    2945    {
    3046        chain->GetEntry(i);
    3147
    32         TLeaf *leaf = chain->GetLeaf(str+name);
    33         if (!leaf)
    34             return 0;
     48        if (numt != chain->GetTreeNumber())
     49        {
     50            if (!(leaf = chain->GetLeaf(str)))
     51                return;
     52
     53            numt = chain->GetTreeNumber();
     54        }
    3555
    3656        Double_t val = leaf->GetValue();
    … …  
    3959            continue;
    4060
    41         if (val < min)
     61        if (val < min && (zero || val!=0))
    4262            min = val;
    4363        if (val > max)
    … …  
    4868    max *= conv;
    4969
    50     chain.GetPlayer();
    51     TTreePlayer &play = *chain.GetPlayer();
    52 
    53     Int_t n;
    54     play.FindGoodLimits(nbins, n, min, max, kFALSE);
    55     nbins = n;
     70    Int_t newn=0;
     71    MH::FindGoodLimits(nbins, newn, min, max, kFALSE);
     72    nbins = newn;
    5673}
    5774
    5875void draw1h1426(Double_t E, Double_t plot=1)
    5976{
    60     plot=2;
    61     Double_t level = log10(E);
    62     level -= log10(5.73e-0*pow( 0.39e3, plot));
    63 
    64     cout << "Level: " << level << endl;
    65 
    66     TPolyMarker *m = new TPolyMarker(7);
    67     /*
    68      m->SetPoint(0,  0.28e3, 1.86e-0*pow( 0.28e3, plot)*level);
    69      m->SetPoint(1,  0.39e3, 5.73e-0*pow( 0.39e3, plot)*level);
    70      m->SetPoint(2,  0.80e3, 1.22e-0*pow( 0.80e3, plot)*level);
    71      m->SetPoint(3,  1.55e3, 5.10e-2*pow( 1.55e3, plot)*level);
    72      m->SetPoint(4,  2.82e3, 1.50e-2*pow( 2.82e3, plot)*level);
    73      m->SetPoint(5,  5.33e3, 7.70e-3*pow( 5.33e3, plot)*level);
    74      m->SetPoint(6, 10.00e3, 1.20e-4*pow(10.00e3, plot)*level);
    75      */
    76     m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot))+level);
    77     m->SetPoint(1, log10( 0.39e3), log10(5.73e-0*pow( 0.39e3, plot))+level);
    78     m->SetPoint(2, log10( 0.80e3), log10(1.22e-0*pow( 0.80e3, plot))+level);
    79     m->SetPoint(3, log10( 1.55e3), log10(5.10e-2*pow( 1.55e3, plot))+level);
    80     m->SetPoint(4, log10( 2.82e3), log10(1.50e-2*pow( 2.82e3, plot))+level);
    81     m->SetPoint(5, log10( 5.33e3), log10(7.70e-3*pow( 5.33e3, plot))+level);
    82     m->SetPoint(6, log10(10.00e3), log10(1.20e-4*pow(10.00e3, plot))+level);
     77    plot += 1;
     78
     79    const Double_t level = E / (5.73e-0*pow( 0.39e3, plot));
     80
     81    TPolyLine   *l = new TPolyLine(6);
     82    TPolyMarker *m = new TPolyMarker(6);
     83
     84    //  m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
     85    m->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
     86    m->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
     87    m->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
     88    m->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
     89    m->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
     90    m->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
    8391    m->SetMarkerStyle(kOpenStar);
    84 //    m->SetMarkerSize(1);
     92    m->SetMarkerColor(kMagenta);
     93    m->SetBit(kCanDelete);
    8594    m->Draw();
     95
     96    //  l->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
     97    l->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
     98    l->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
     99    l->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
     100    l->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
     101    l->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
     102    l->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
     103    l->SetLineColor(kMagenta);
     104    l->SetBit(kCanDelete);
     105    l->Draw();
    86106}
    87107
    88108void analp()
    89109{
     110    // -------------------------------------------------------------------
     111
    90112    TChain chain("Photons");
     113    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_01.root");
     114    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_02.root");
     115    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_03.root");
     116    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_04.root");
     117    chain.Add("cascade_0.1_18_1e2_1e5_B0_512_01.root");
     118    chain.Add("cascade_0.1_18_1e2_1e5_B0_512_03.root");
     119    chain.Add("cascade_0.1_18_1e2_1e5_B0_512_02.root");
    91120    /*
    92      chain.Add("delme3.root");
    93      chain.Add("delme3B.root");
    94      chain.Add("delme3C.root");
    95      chain.Add("delme3D.root");
    96      chain.Add("delme3E.root");
     121     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_3.root");
     122     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_4.root");
     123     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_5.root");
     124     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_6.root");
     125     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_7.root");
     126     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_8.root");
     127     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_9.root");
     128
     129     chain.Add("cascade_0.03_18_1e2_1e5_B0_512_01.root");
     130     chain.Add("cascade_0.03_18_1e2_1e5_B0_512_02.root");
    97131     */
    98     chain.Add("delme3H.root");
    99 
    100     Int_t nbins =  24;
    101     Double_t lo = 1e2;
    102     Double_t hi = 1e6;
    103 
    104     MBinning binsx;
    105     binsx.SetEdgesLog(nbins, lo, hi);
    106 
    107     // ------------------------------
    108 
    109     Double_t cutoff = 1e12;
    110 
    111     Double_t alpha = -1;
    112     Double_t plot  =  1;
    113 
    114     Int_t nbins = 16;
    115 
    116     // ------------------------------
    117 
    118     Int_t n = chain.GetEntries();
     132
     133    Int_t nbinse = 18;      // number of bins in your histogram
     134    Double_t lo  = 1e2/1e3; // lower limit of the energy histogram
     135    Double_t hi  = 1e5;     // upper limit of the energy histogram
     136
     137    Double_t minE = 1e2;    // minimum value produced in the simulation
     138
     139    Double_t cutoff = 1e12; // arbitrary energy cutoff (throw away the events)
     140
     141    Double_t alpha = -1;    // alpha of the energy spectrum (-1 is E^-2)
     142    Double_t plot  =  1;    // plot of the enrgy spectrum (1 is E^2)
     143
     144    Int_t nbins = 16;       // number of bins you want to have in psi/phi
     145
     146    // -------------------------------------------------------------------
     147
     148    MBinning binspolx;
     149    binspolx.SetEdges(16, -180, 180);
     150
     151    // -------------------------------------------------------------------
     152
     153    Stat_t n = chain.GetEntries();
    119154
    120155    cout << endl << "Found " << n << " entries." << endl;
    … …  
    123158        return;
    124159
    125     MBinning binspolx;
    126     MBinning binspola;
    127     MBinning binspolr;
    128     binspolx.SetEdges(16, -180, 180);
     160    // -----------------------------------
     161
     162    MBinning bins;
     163
     164    cout << "Determin min/max..." << endl;
    129165
    130166    Double_t min, max;
    131167    GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
    132     binspola.SetEdges(nbins, min, max*1.1);
    133168    cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;
    134169
     170    bins.SetEdges(nbins, min, max*1.1);
     171
     172    TH2D a, a2, a3;
     173    MH::SetBinning(&a,  &binspolx, &bins);
     174    MH::SetBinning(&a2, &binspolx, &bins);
     175    MH::SetBinning(&a3, &binspolx, &bins);
     176
     177    // -----------------------------------
     178
    135179    GetRange(&chain, "fR", nbins, min, max);
    136     binspolr.SetEdges(nbins, min, max*1.1);
    137180    cout << "fR,     " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;
    138181
    139     TH1D h, h2, h3;
    140     TH1D prim;
    141     MH::SetBinning(&h, &binsx);
    142     MH::SetBinning(&h2, &binsx);
    143     MH::SetBinning(&h3, &binsx);
     182    bins.SetEdges(nbins, min, max*1.1);
     183
     184    TH2D r, r2, r3;
     185    MH::SetBinning(&r,  &binspolx, &bins);
     186    MH::SetBinning(&r2, &binspolx, &bins);
     187    MH::SetBinning(&r3, &binspolx, &bins);
     188
     189    // -----------------------------------
     190
     191    MBinning binsx;
     192    binsx.SetEdgesLog(nbinse, lo, hi);
     193
     194    TH1D prim, h, h2, h3;
     195    MH::SetBinning(&h,    &binsx);
     196    MH::SetBinning(&h2,   &binsx);
     197    MH::SetBinning(&h3,   &binsx);
    144198    MH::SetBinning(&prim, &binsx);
    145 
    146     TH2D r, r2, r3;
    147     TH2D a, a2, a3;
    148     MH::SetBinning(&a, &binspolx, &binspola);
    149     MH::SetBinning(&r, &binspolx, &binspolr);
    150     MH::SetBinning(&a2, &binspolx, &binspola);
    151     MH::SetBinning(&r2, &binspolx, &binspolr);
    152     MH::SetBinning(&a3, &binspolx, &binspola);
    153     MH::SetBinning(&r3, &binspolx, &binspolr);
     199//    prim.Sumw2();
     200    h.Sumw2();
     201    h2.Sumw2();
     202    h3.Sumw2();
     203
     204    // -----------------------------------
    154205
    155206    MPhoton *p = new MPhoton;
    156207    chain.SetBranchAddress("MPhoton.", &p);
     208    chain.SetBranchStatus("*", 1);
    157209
    158210    chain.GetEntry(0);
    … …  
    165217    cout << "Z = " << z << endl;
    166218    cout << "R = " << R << "kpc" << endl;
     219
     220    // -----------------------------------
     221
     222    const Double_t skpc = 3600.*24.*365.*3.258*1e3;  // [s/kpc]
     223    GetRange(&chain, "fX", nbins, min, max, 1, kFALSE);
     224    min *= skpc;
     225    max *= skpc;
     226    if (min<1e-2)
     227        min=1e-2;
     228    bins.SetEdgesLog(nbins, min, max);
     229    cout << "dX,     " << nbins << ": " << min << " ... " << max << " [s]" << endl;
     230
     231    TH1D t;
     232    MH::SetBinning(&t, &bins);
     233    t.Sumw2();
     234
     235    // ----------------------------------------------------------------------
     236
     237    chain.SetBranchAddress("MPhoton.", &p);
     238    chain.SetBranchStatus("*", 1);
     239
     240    cout << "Filling: " << nbinse << " bins @ " << minE << " - " << hi << "... " << flush;
    167241
    168242    Double_t weight = 0;
    … …  
    193267            continue;
    194268
    195         h.Fill(Ep, pow(Ep,plot) * weight);
    196 
    197269        Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180;
    198270        Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;
    199 
    200         r.Fill(phi, isprim?0:p->GetR(),                 weight);
    201         a.Fill(psi, isprim?0:p->GetTheta()*kRad2Deg*60, weight);
    202271
    203272        if (isprim)
    204273        {
    205274            h3.Fill(Ep, pow(Ep,plot) * weight);
    206             r3.Fill(phi, 0,                 weight);
     275            r3.Fill(phi, 0, weight);
    207276            a3.Fill(psi, 0, weight);
    208277            isprim = kFALSE;
    … …  
    211280
    212281        h2.Fill(Ep, pow(Ep,plot) * weight);
    213         r2.Fill(phi, p->GetR(),                 weight);
     282        r2.Fill(phi, p->GetR(), weight);
    214283        a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
     284
     285        Double_t tm1 = sqrt(R*R+p->GetR()*p->GetR());
     286        Double_t tm2 = (R+p->GetX())/tm1 - 1;
     287        t.Fill(p->GetX()<=0 ? 0 : R*tm2*skpc, weight);
    215288    }
    216289
     290    // ----------------------------------------------------------------------
     291
    217292    delete p;
    218293
    219     TH1 *g=r.ProjectionX();
    220     cout << "Mean fPhi: " << g->GetMean() << " +- " << g->GetRMS() << endl;
    221     delete g;
    222     g=a.ProjectionX();
    223     cout << "Mean fPsi: " << g->GetMean() << " +- " << g->GetRMS() <<  endl;
    224     delete g;
     294    cout << "Done." << endl;
     295
     296    Double_t N = prim.GetBinContent(nbinse);
     297    Double_t Nerr = prim.GetBinError(nbinse);
     298    for (int i=1; i<=nbinse; i++)
     299    {
     300        Double_t E = prim.GetXaxis()->GetBinCenter(i);
     301        if (E>minE)
     302            continue;
     303
     304        Double_t E1 = prim.GetXaxis()->GetBinLowEdge(i+1);
     305        Double_t E0 = prim.GetXaxis()->GetBinLowEdge(i);
     306
     307        E = sqrt(E1*E0);
     308
     309        prim.SetBinContent(i, N);
     310        h3.SetBinContent(i, N);
     311        prim.SetBinError(i, Nerr);
     312        h3.SetBinError(i, Nerr);
     313    }
     314
     315    h.Add(&h2, &h3);
     316    r.Add(&r2, &r3);
     317    a.Add(&a2, &a3);
     318
     319    if (alpha==-plot)
     320    {
     321        cout << "Observed Energy: " << h.Integral() << "GeV" << endl;
     322        cout << "Emitted Energy:  " << prim.Integral() << "GeV" << endl;
     323        cout << "Ratio: " << h.Integral()/prim.Integral() << endl;
     324    }
     325
     326    cout << "Mean fPhi: " << r.GetMean() << " +- " << r.GetRMS() << endl;
     327    cout << "Mean fPsi: " << a.GetMean() << " +- " << a.GetRMS() <<  endl;
    225328
    226329    gStyle->SetOptStat(10);
     330    gStyle->SetPalette(1, 0);
    227331
    228332    TLine line;
    229333
     334    // ----------------------------------------------------------------------
     335
    230336//    delete gROOT->FindObject("Analysis Arrival");
    231     c = MH::MakeDefCanvas("Analysis Arrival", "");
     337    TCanvas *c = MH::MakeDefCanvas("Analysis Arrival", "");
    232338    c->Divide(2,2);
    233 
    234     gStyle->SetPalette(1, 0);
    235339
    236340    c->cd(1);
    … …  
    241345    r.SetXTitle("\\Phi [\\circ]");
    242346    r.SetYTitle("R [kpc]");
    243     TPad &pad = *gPad;
     347    TVirtualPad &pad = *gPad;
    244348    pad.Divide(2,2);
    245349    pad.cd(1);
    … …  
    248352    gPad->SetTheta(0);
    249353    gPad->SetPhi(90);
    250     g = r.DrawCopy("surf1pol");
     354    TH1 *g = r.DrawCopy("surf1pol");
    251355    pad.cd(2);
    252356    gPad->SetLogy();
    … …  
    275379    a.SetXTitle("\\Psi [\\circ]");
    276380    a.SetYTitle("\\Theta [']");
    277     TPad &pad = *gPad;
    278     pad.Divide(2,2);
    279     pad.cd(1);
     381    TVirtualPad &pad2 = *gPad;
     382    pad2.Divide(2,2);
     383    pad2.cd(1);
    280384    gPad->SetLogy();
    281385    gPad->SetLogz();
    … …  
    283387    gPad->SetPhi(90);
    284388    g = a.DrawCopy("surf1pol");
    285     pad.cd(2);
     389    pad2.cd(2);
    286390    gPad->SetLogy();
    287391    gPad->SetLogz();
    … …  
    289393    gPad->SetPhi(90);
    290394    g->Draw("surf1pol");
    291     pad.cd(3);
     395    pad2.cd(3);
    292396    gPad->SetLogy();
    293397    gPad->SetLogz();
    … …  
    295399    gPad->SetPhi(90);
    296400    g->Draw("surf1pol");
    297     pad.cd(4);
     401    pad2.cd(4);
    298402    gPad->SetLogy();
    299403    gPad->SetLogz();
    … …  
    307411    g->SetDirectory(NULL);
    308412    TH1 *g2=r2.ProjectionY();
    309     g->SetMinimum(GetMin(g2)/2);
     413    g->SetMinimum(MH::GetMinimumGT(*g2)/2);
    310414    g->SetBit(kCanDelete);
    311415    g->SetTitle(" Hit Observers Plain ");
    … …  
    329433    gPad->SetLogy();
    330434    g=a.ProjectionY();
    331     g->SetMinimum(GetMin(g)/2);
     435    g->SetMinimum(MH::GetMinimumGT(*g)/2);
    332436    g->SetDirectory(NULL);
    333437    g->SetBit(kCanDelete);
    … …  
    350454    g->Draw("same");
    351455
    352   //  delete gROOT->FindObject("Analysis Photons");
    353     TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
     456    // ----------------------------------------------------------------------
     457
     458  //  delete gROOT->FindObject("Analysis Photons");
     459    c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
    354460    c->Divide(1,2);
    355461
    … …  
    364470    h.GetXaxis()->SetLabelOffset(-0.015);
    365471    h.GetXaxis()->SetTitleOffset(1.1);
    366     h.SetMarkerStyle(kMultiply);
     472    h.SetMarkerStyle(kPlus);
    367473    prim.SetFillStyle(0);
    368474    prim.SetMarkerStyle(kPlus);
    369475    prim.SetMarkerColor(kRed);
    370476    prim.SetLineColor(kRed);
    371     TH1 &g1=*h.DrawCopy("P");
    372     g2=*prim.DrawCopy("Psame");
     477    TH1 &g1=*h.DrawCopy("P E4");
     478    g2=prim.DrawCopy("Psame E0");
    373479    g2->Draw("Csame");
    374480    g1.Draw("Csame");
    … …  
    378484    h2.SetMarkerColor(kBlue);
    379485    h2.SetLineColor(kBlue);
    380     TH1 &g3=*h2.DrawCopy("Psame");
     486    TH1 &g3=*h2.DrawCopy("Psame E2");
    381487    g3.Draw("Csame");
    382488    h3.SetFillStyle(0);
    … …  
    385491    h3.SetMarkerColor(kGreen);
    386492    h3.SetLineColor(kGreen);
    387     TH1 &g4=*h3.DrawCopy("Psame");
     493    TH1 &g4=*h3.DrawCopy("Psame E1");
    388494    g4.Draw("Csame");
    389495
    … …  
    394500    TH1D div;
    395501    MH::SetBinning(&div, &binsx);
     502    div.Sumw2();
    396503    div.Divide(&h, &prim);
    397     div.SetTitle("Spectrum / Primary Spectrum");
     504    div.SetTitle(" Spectrum / Primary Spectrum ");
    398505    div.GetXaxis()->SetLabelOffset(-0.015);
    399506    div.GetXaxis()->SetTitleOffset(1.1);
    400507    div.SetXTitle("E\\gamma [GeV]");
    401508    div.SetYTitle("Ratio");
     509    TH1 *gHist = div.DrawCopy("E4");
     510
     511    line.SetLineWidth(1);
     512    line.SetLineColor(kGreen);
     513    line.DrawLine(log10(lo), exp(-1), log10(hi), exp(-1));
    402514    line.SetLineColor(kBlue);
    403     line.SetLineWidth(1);
    404     div.DrawCopy();
    405515    line.DrawLine(log10(lo), 1, log10(hi), 1);
     516
     517    MFit fit("exp(-1)-[1]*log10(x/[0])");
     518    for (int i=0; i<gHist->GetEntries(); i++)
     519    {
     520        if (gHist->GetBinContent(i+1)<exp(-1))
     521        {
     522            fit.SetRange(gHist->GetBinCenter(i-2), gHist->GetBinCenter(i+2));
     523            cout << "Fitting from " << gHist->GetBinLowEdge(i-1) << " to ";
     524            cout << gHist->GetBinLowEdge(i+1) << endl;
     525            break;
     526        }
     527    }
     528    fit.SetParameter(0, "t", 750,  10, 1e5);
     529    fit.SetParameter(1, "m", 0.5, 0.1,  10);
     530    fit.FitLog(gHist);
     531    fit.Print();
     532    fit.DrawCopy("same");
     533
     534    cout << "Cutoff:  " << setprecision(2) << fit[0]/1e3 << "TeV +- ";
     535    cout << fit.GetParError(0) << endl;
     536
     537     // ----------------------------------------------------------------------
     538
     539    c = MH::MakeDefCanvas("Time Analysis", "", 580, 435);
     540    gPad->SetLogx();
     541    //    gPad->SetLogy();
     542    t.SetName("Times");
     543    t.SetTitle(" Arrival time distribution ");
     544    t.SetXTitle("\\Delta  t [s]");
     545    t.GetXaxis()->SetLabelOffset(-0.015);
     546    t.GetXaxis()->SetTitleOffset(1.1);
     547    t.DrawCopy("E4");
     548    line.DrawLine(log10(1),           0, log10(1),           t.GetMaximum()*1.05);
     549    line.DrawLine(log10(3600),        0, log10(3600),        t.GetMaximum()*1.05);
     550    line.DrawLine(log10(3600*24),     0, log10(3600*24),     t.GetMaximum()*1.05);
     551    line.DrawLine(log10(3600*24*365), 0, log10(3600*24*365), t.GetMaximum()*1.05);
    406552}

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

    r1448 r1449  
    1 #include <iostream.h>
     1#include "MCascade.h"
    22
    3 #include <fstream.h>
     3void phys()
     4{
     5    MCascade cascade;
    46
    5 #include <TStopwatch.h>
    6 #include <TF1.h>
    7 #include <TH1.h>
    8 #include <TH2.h>
    9 #include <TList.h>
    10 #include <TFile.h>
    11 #include <TTree.h>
    12 #include <TTimer.h>
    13 #include <TStyle.h>
    14 #include <TBranch.h>
    15 #include <TRandom.h>
    16 #include <TCanvas.h>
     7    cascade.SetSourceZ(0.1);             // Readshift of source
     8    cascade.SetB(0/*1e-6*/);             // [G] mean magnetic field
     9    cascade.SetBradius(50);              // [Mpc]
     10    cascade.SetRuntime(8*60);            // [min] maximum time to run the simulation
     11    cascade.SetEnergyBins(18, 1e2, 1e5); // [GeV]
     12    cascade.SetMaxInvCompton(512);       // [#] maximum number of inv. Compton (-1 means infinite)
     13    cascade.SetRatioInvCompton(0);       // [%] allowed Emis of electron (0 means disabled)
     14    cascade.SetSpectralIndex(-1);        // -1 means a E^2 plot
     15    cascade.SetDisplayIndex(1);          //  1 means a E^-2 spectrum
    1716
    18 #include "mbase/MParContainer.h"
    19 #include "mphys/MPhoton.h"
    20 #include "mphys/MElectron.h"
    21 #include "mphys/MPairProduction.h"
    22 #include "mhist/MBinning.h"
    23 #include "mhist/MH.h"
    24 
    25 // 2.96
    26 // 2.87
    27 // 2.73
    28 Double_t PrimSpect(Double_t *x, Double_t *k)
    29 {
    30     return pow(pow(10, x[0]), k[0]);
     17    //
     18    // Run the simulation: filename, eventdisply (on/off)
     19    //
     20    //cascade.Run("cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_256_1.root", kTRUE);
     21    cascade.Run("cascade_0.1_18_1e2_1e5_B0_512_02.root", kFALSE);
    3122}
    3223
    33 Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
    34 {
    35     Double_t Ep = pow(10, x[0]);
    36 
    37     Double_t res = MPhoton::Int2(&Ep, k);
    38     return res*1e55; //65/k[0];
    39 
    40     //return MPhoton::Planck(&Ep, &k[1]);
    41 }
    42 
    43 Double_t Sbar_sigmas(Double_t *x, Double_t *k)
    44 {
    45     Double_t sbar = pow(10, x[0]);
    46 
    47     Double_t s = 1./(sbar*4);
    48 
    49     Double_t sigma = MPhoton::Sigma_gg(&s);
    50 
    51     return sigma*sbar*1e28;
    52 }
    53 
    54 Double_t RandomTheta(Double_t Eg, Double_t Ep)
    55 {
    56     Double_t E0 = 511e-6; // [GeV]
    57 
    58     Double_t f = Eg/E0*Ep/E0;
    59 
    60     if (f<1)
    61         return 0;
    62 
    63     TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
    64 
    65     func.SetNpx(50);
    66 
    67     Double_t sbar  = pow(10, func.GetRandom());
    68     Double_t theta = acos(1.-sbar*2/f);
    69 
    70     return theta;
    71 }
    72 
    73 Double_t GetEnergy(Int_t n, Double_t lo, Double_t hi)
    74 {
    75     static int bin=0;
    76     static TRandom rnd(0);
    77 
    78     Double_t w = log10(hi/lo)/n;
    79 
    80     Double_t E = lo*pow(10, rnd.Uniform(w) + w*(n-bin-1));
    81 
    82     //cout << endl << w << " " << n << " " << w*bin << " " << E << endl;
    83 
    84     if (++bin==n)
    85         bin=0;
    86 
    87     return E;
    88 
    89     /*
    90     // TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
    91     // src.SetParameter(0, 0);
    92     // Double_t E = pow(10, src.GetRandom());
    93 
    94     //
    95     // 0: 11111111111111|22222222222222   2   2^0 + 1    1
    96     // 1: 11111111|22222222222|33333333   3   2^1 + 1    2
    97     // 2: 11111|22222|33333|44444|55555   5   2^2 + 1    7
    98     // 3:  |     | |   | |  |  |   |                    15
    99     //
    100 
    101     static int run=0;
    102     static int num=1;
    103 
    104     Double_t w = log10(hi/lo)/((1<<run) + 1);
    105 
    106     Double_t E = lo*pow(10, num*w);
    107 
    108     if (num == (1<<run))
    109     {
    110         run++;
    111         num=1;
    112     }
    113     else
    114         num++;
    115 
    116         return E;
    117         */
    118 }
    119 
    120 void DoIt()
    121 {
    122     // a -->  5
    123     // b --> 10
    124     // c -->  2
    125 
    126     // delme, delmeB starting integration at lolim
    127     // delme2, delme2B starting integration at -log10(Eg)-8
    128     // delme3,  lolim, 24, 1e2-1e6,   2x inv.Compton
    129     // delme3B, lolim, 24, 1e2-1e6,   4x inv.Compton
    130     // delme3C, lolim, 24, 1e2-1e6,   8x inv.Compton
    131     // delme3D, lolim, 24, 1e2-1e6,  16x inv.Compton
    132     // delme3E, lolim, 24, 1e2-1e6,  32x inv.Compton
    133 
    134     // delme3H, lolim, 24, 1e2-1e6, 256x inv.Compton 8*60
    135 
    136     Double_t startz = 0.03; //MParticle::ZofR(&R);
    137     Double_t R      = MParticle::RofZ(&startz); // [kpc]
    138 
    139     const char *filename = "cascade_0.03_18_1e2_1e5_B0_256_1.root";
    140     //const char *filename = "test.root";
    141 
    142     const Double_t B = 0;//1e-6*1e-4; // [T=1e4G] mean magnetic field
    143 
    144     Double_t runtime = 45*60;     // [s] maximum time to run the simulation
    145 
    146     Int_t nbins = 18;      // number of bins produced in energy spectrum
    147 
    148     Double_t lo = 1e2;     // lower limit of displayed spectrum
    149     Double_t hi = 1e5;     // upper limit of spectrum (cutoff)
    150 
    151     Int_t counter = 256;   // maximum number of inv. Compton (-1 means infinite)
    152 
    153     Double_t alpha = -1;   // -1 means a E^2 plot
    154     Double_t plot  =  1;   //  1 means a E^-2 spectrum
    155 
    156     Double_t bubbler = R-1e3*10; // [kpc]
    157     Double_t bubblez = MParticle::ZofR(&bubbler);
    158 
    159     // --------------------------------------------------------------------
    160 
    161     cout << "R = " << R << "kpc" << endl;
    162     cout << "Z = " << startz << endl;
    163 
    164     cout << "Setting up: Histograms... " << flush;
    165 
    166     MPairProduction pair;
    167 
    168     TH1D hist;
    169     TH1D histsrc;
    170 
    171     hist.SetMinimum(pow(lo, plot+alpha)/100);
    172     histsrc.SetMinimum(pow(lo, plot+alpha)/100);
    173 
    174     TList listg;
    175     TList liste;
    176 
    177     listg.SetOwner();
    178     liste.SetOwner();
    179 
    180     gStyle->SetOptStat(10);
    181 
    182     //
    183     // Don't change the order!!!
    184     //
    185     hist.SetFillStyle(0);
    186     hist.SetMarkerStyle(kPlus);
    187     histsrc.SetFillStyle(0);
    188     histsrc.SetMarkerStyle(kMultiply);
    189     histsrc.SetMarkerColor(kRed);
    190     histsrc.SetLineColor(kRed);
    191 
    192     MBinning bins;
    193     bins.SetEdgesLog(nbins, lo, hi);
    194 
    195     MH::SetBinning(&hist,    &bins);
    196     MH::SetBinning(&histsrc, &bins);
    197 
    198     TCanvas *c=MH::MakeDefCanvas();
    199 
    200     gPad->SetGrid();
    201     gPad->SetLogx();
    202     gPad->SetLogy();
    203 
    204     hist.SetName("Spectrum");
    205     hist.SetXTitle("E [GeV]");
    206     hist.SetYTitle(Form("E^{%.1f} Counts", plot));
    207     hist.GetXaxis()->SetLabelOffset(-0.015);
    208     hist.GetXaxis()->SetTitleOffset(1.1);
    209 
    210     hist.Draw("P");
    211     histsrc.Draw("Psame");
    212     histsrc.Draw("same");
    213     hist.Draw("same");
    214 
    215     // ------------------------------
    216 
    217     cout << "Output File... " << flush;
    218 
    219     TFile file(filename, "RECREATE", "Intergalactic cascade", 9);
    220     cout << "Trees... " << flush;
    221     TTree   *T1 = new TTree ("Photons",    "Photons from Cascade");
    222     TTree   *T2 = new TTree ("Electrons",  "Electrons in the Cascade");
    223     cout << "Branches... " << flush;
    224     MPhoton dummyp;
    225     void *ptr = &dummyp;
    226     TBranch *B1 = T1->Branch("MPhoton.",   "MPhoton",   &ptr);
    227     MPhoton dummye;
    228     ptr = &dummye;
    229     TBranch *B2 = T2->Branch("MElectron.", "MElectron", &ptr);
    230 
    231     // ------------------------------
    232 
    233     cout << "Timers... " << flush;
    234 
    235     TTimer timer("gSystem->ProcessEvents();", 250, kFALSE);
    236     timer.TurnOn();
    237 
    238     TStopwatch clock;
    239     clock.Start();
    240 
    241     cout << "Done. " << endl;
    242 
    243     Int_t n=0;
    244     Double_t starttime = TStopwatch::GetRealTime(); // s
    245     while (TStopwatch::GetRealTime()<starttime+runtime)
    246         for (int i=0; i<nbins; i++)
    247     {
    248         n++;
    249 
    250         Double_t E = GetEnergy(nbins, lo, hi);
    251         Double_t weight = pow(E, alpha);
    252         histsrc.Fill(E, pow(E, plot) * weight);
    253 
    254         MPhoton *gamma=new MPhoton(E, startz);
    255         gamma->SetIsPrimary();
    256         gamma->SetSrcR(R);
    257         gamma->InitRandom();
    258         listg.Add(gamma);
    259 
    260         B1->SetAddress(&gamma);
    261         T1->Fill();
    262 
    263         cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
    264 
    265         Double_t Esum=0;
    266         Double_t Emis=0;
    267         while (1)
    268         {
    269             if (listg.GetSize()!=0)
    270                 cout << " |P" << flush;
    271 
    272             TIter NextP(&listg);
    273             MPhoton *p = NULL;
    274             B1->SetAddress(&p);
    275             while ((p=(MPhoton*)NextP()))
    276             {
    277                 cout << ":" << flush;
    278                 Double_t Eg = p->GetEnergy();
    279 
    280                 Double_t E0 = 511e-6;
    281                 Double_t z  = p->GetZ();
    282                 Double_t lolim = E0*E0/Eg;
    283                 Double_t inf   = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
    284                 if (Eg<5e4)
    285                     inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
    286 
    287                 //TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
    288                 //Double_t E0 = 511e-6; // [GeV]
    289                 TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
    290                 phot.SetParameter(0, Eg);
    291                 while (1)
    292                 {
    293                     if (!p->SetNewPosition())
    294                     {
    295                         cout << "!" << flush;
    296 
    297                         hist.Fill(Eg, pow(Eg, plot)*weight);
    298 
    299                         p->SetIsPrimary(kFALSE);
    300                         T1->Fill();
    301 
    302                         Esum += Eg;
    303 
    304                         break;
    305                     }
    306 
    307                     //
    308                     // Sample phtoton from background and interaction angle
    309                     //
    310                     phot.SetParameter(1, p->GetZ());
    311                     Double_t pe = phot.GetRandom();
    312                     if (pe==0)
    313                     {
    314                         cout << "z" << flush;
    315                         continue;
    316                     }
    317 
    318                     Double_t Ep = pow(10, pe);
    319                     Double_t theta = RandomTheta(Eg, Ep);
    320                     if (theta==0)
    321                     {
    322                         cout << "t" << flush;
    323                         continue;
    324                     }
    325 
    326                     if (!pair.Process(p, Ep, theta, &liste))
    327                     {
    328                         cout << "0" << flush;
    329                         continue;
    330                     }
    331 
    332                     cout << "." << flush;
    333                     break;
    334                 }
    335                 delete listg.Remove(p);
    336             }
    337 
    338             if (liste.GetSize()==0 && listg.GetSize()==0)
    339                 break;
    340 
    341             if (liste.GetSize()!=0)
    342                 cout << " |E" << flush;
    343 
    344             TIter Next(&liste);
    345             MElectron *e = NULL;
    346             B2->SetAddress(&e);
    347             while ((e=(MElectron*)Next()))
    348             {
    349                 e->SetIsPrimary(kTRUE);
    350                 T2->Fill();
    351                 e->SetIsPrimary(kFALSE);
    352 
    353                 Double_t Ee = e->GetEnergy();
    354 
    355                 if (Ee<lo)
    356                     continue;
    357 
    358                 cout << ":" << flush;
    359                 int test = counter<0 ? -1 : 0;
    360                 while (test<0 ? true : (test++<counter))
    361                 {
    362                     if (!e->SetNewPositionB(e->GetZ()>bubblez ? B : 0))
    363                     {
    364                         cout << "!" << flush;
    365                         break;
    366                     }
    367 
    368                     // WRONG!
    369                     MPhoton *p = e->DoInvCompton(0);
    370 
    371                     T2->Fill();
    372 
    373                     cout << "." << flush;
    374                     listg.Add(p);
    375                     /*
    376                     if (fabs(p->GetTheta()*3437)<60 &&  // < 60min
    377                         p->GetEnergy()>lo)
    378                     {
    379                         cout << "." << flush;
    380                         listg.Add(p);
    381                     }
    382                     else
    383                     {
    384                         if (p->GetEnergy()<E*pow(10, alpha)/5 || p->GetEnergy()<=lo)
    385                             cout << "e" << flush;
    386                         else
    387                             cout << "t" << flush; // ignored
    388                         delete p;
    389                         Esum += p->GetEnergy();
    390                         Emis += p->GetEnergy();
    391                     }
    392                     */
    393                     if (fabs(e->GetTheta()*3437)>60)  // < 60min
    394                     {
    395                         cout << "T" << flush;
    396                         break;
    397                     }
    398 
    399                     if (e->GetEnergy()<Ee*1e-3) // <2e3
    400                     {
    401                         cout << "E" << flush;
    402                         break;
    403                     }
    404                 }
    405                 Esum += e->GetEnergy();
    406                 Emis += e->GetEnergy();
    407             }
    408             liste.Delete();
    409         }
    410         cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
    411 
    412         hist.SetTitle(Form("E^{%.1f}  z=%f  T=%d'%d\"  N=%d", alpha, startz,
    413                               (int)runtime/60, (int)runtime%60, n));
    414 
    415         timer.Stop();
    416         c->Update();
    417         timer.Start(250);
    418 
    419     }
    420     cout << endl;
    421 
    422     clock.Stop();
    423     clock.Print();
    424 
    425     timer.Stop();
    426 
    427     file.Write();
    428 
    429     cout << "Wrote: " << filename << endl;
    430     cout << "Created " << n << " gammas from source with E^" << alpha << endl;
    431     cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl;
    432 
    433     // ------------------------------
    434 
    435     gPad->Clear();
    436 
    437     hist.SetTitle(Form("E^{%.1f}  z=%f  T=%d'%d\"  N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
    438 
    439     TH1 &h1 = *hist.DrawCopy("P");
    440     TH1 &h2 = *histsrc.DrawCopy("Psame");
    441     h2.Draw("Csame");
    442     h1.Draw("Csame");
    443 }
    444 
    445 // -------------------------------------------------------------------
    446 void phys()
    447 {
    448     /*
    449     Double_t Eg = 10;
    450     Double_t E0 = 511e-6;
    451     Double_t z  = 0.03;
    452     Double_t lolim = E0*E0/Eg;
    453     Double_t inf   = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
    454     if (Eg<5e4)
    455         inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
    456     TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
    457     phot.SetParameter(0, Eg);
    458     phot.SetParameter(1, z);
    459 
    460     Double_t val[2] = {Eg,z};
    461     cout << phot.Integral(log10(lolim), log10(inf), val, 1e-2) << endl;
    462 
    463     cout << lolim << " " << inf << endl;
    464 
    465     phot.DrawCopy();
    466     return;
    467 */
    468 
    469     DoIt();
    470 
    471     /*
    472      Double_t E = 1e10;
    473      TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
    474      phot.SetParameter(0, E);
    475      phot.SetParameter(1, 5);
    476      phot.DrawCopy();
    477      return;
    478      */
    479 
    480     /*
    481      Double_t Eg = 1e5;
    482 
    483      Double_t E0    = 511e-6;                  // [GeV]
    484      Double_t lolim = E0*E0/Eg;
    485      Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
    486 
    487      // 1e5:   5e-13, 4e-12
    488      // 1e6:   5e-14, 2e-12
    489      // 1e8:   5e-16, 1e-12
    490      // 1e10:  5e-18, 1e-12
    491 
    492      // 1e5:   2.6e-12, 4e-12
    493      // 1e6:   2.6e-13, 2e-12
    494      // 1e8:   2.6e-15, 1e-12
    495      // 1e10:  1.6e-17, 1e-12
    496 
    497      TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
    498      f.SetParameter(0, Eg);
    499      f.SetParameter(1, 0);
    500 
    501      MH::MakeDefCanvas();
    502      gPad->SetLogx();
    503      f.DrawCopy();
    504     */
    505 }
    506 
Note: See TracChangeset for help on using the changeset viewer.