/////////////////////////////////////////////////////////////////////// // // MParticle // /////////////////////////////////////////////////////////////////////// #include "MParticle.h" #include #include #include ClassImp(MParticle); /************************************************** * * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] * **************************************************/ Double_t MParticle::ZofR(Double_t *x, Double_t *k) { /* const Double_t c = 299792458; // [m/s] const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] const Double_t ly = 3600.*24.*365.*c; // [m/ly] const Double_t pc = 1./3.258; // [pc/ly] const Double_t r = x[0] /pc*ly*1e3; // [m] const Double_t R = r*H0/c; // [1] return (R+1+sqrt(R*2+1))/2 - 1; */ const Double_t c = 299792458; // [m/s] const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] const Double_t ly = 3600.*24.*365.*c; // [m/ly] const Double_t pc = 1./3.258; // [pc/ly] const Double_t r = x[0] /pc*ly*1e3; // [m] const Double_t R = 1./(1.-r*H0/c/2); // [1] return R*R - 1; } Double_t MParticle::RofZ(Double_t *x, Double_t *k) { /* Double_t z1 = x[0] + 1; const Double_t c = 299792458; // [m/s] const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] const Double_t ly = 3600.*24.*365.*c; // [m/ly] const Double_t pc = 1./3.258; // [pc/ly] const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m] return R * pc/ly/1e3; // [kpc] */ Double_t z1 = x[0] + 1; const Double_t c = 299792458; // [m/s] const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] const Double_t ly = 3600.*24.*365.*c; // [m/ly] const Double_t pc = 1./3.258; // [pc/ly] const Double_t R = c/H0 * 2 * (1. - 1./sqrt(z1)); // [m] return R * pc/ly/1e3; // [kpc] } #include #include #include "../mhist/MBinning.h" #include "../mhist/MH.h" TH2D *hist2; Double_t MParticle::Planck(Double_t *x, Double_t *k) { static Bool_t isloaded = kFALSE; if (!isloaded) { Double_t c = 299792458; // [m/s] Double_t e = 1.602176462e-19; // [C] Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] Double_t hc = h*c; // [GeVm] Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); ifstream fin("mphys/background.txt"); hist2 = new TH2D; MBinning binsz; MBinning binse; binsz.SetEdgesLog(100, 1e-6, 1); // --> 101 Edges / 100 bins binse.SetEdgesLog(100, 7e-15, 3e-8); // --> 101 Edges / 100 bins MH::SetBinning(hist2, &binsz, &binse); for (int y=0; y<101; y++) { Double_t val; fin >> val; for (int x=0; x<101; x++) { fin >> val; val += 9; Double_t z = binsz.GetEdges()[x]; Double_t f = z+1; Double_t newval = pow(10, val)/konst; hist2->SetBinContent(x, y, newval*f*f*f); } } isloaded = kTRUE; } static TAxis &axez = *hist2->GetXaxis(); static TAxis &axee = *hist2->GetYaxis(); // // y = (y1-y0)/(x1-x0) * (x-x0) + y0 // Double_t z = k ? k[0] : 0; Double_t E = x[0]; Int_t i = axez.FindFixBin(z); Int_t j = axee.FindFixBin(E); Double_t z1 = axez.GetBinLowEdge(i+1); Double_t z0 = axez.GetBinLowEdge(i); Double_t E1 = axee.GetBinLowEdge(j+1); Double_t E0 = axee.GetBinLowEdge(j); Double_t n00 = hist2->GetBinContent(i, j); Double_t n01 = hist2->GetBinContent(i+1, j); Double_t n10 = hist2->GetBinContent(i, j+1); Double_t n11 = hist2->GetBinContent(i+1, j+1); Double_t dz = (z-z0)/(z1-z0); Double_t dE = (E-E0)/(E1-E0); Double_t n0 = dz*(n01-n00)+n00; Double_t n1 = dz*(n11-n10)+n10; Double_t n = dE*(n1-n0)+n0; return n; /* // // TANJA2 // Double_t E1 = x[0]; Double_t E2 = x[0]/8; return (MParticle::Planck0(&E1, k)+MParticle::Planck0(&E2, k)/40e3)*0.7/0.4; */ /* // // TANJA // Double_t E1 = x[0]; Double_t E2 = x[0]/8; return Planck0(&E1, k)+Planck0(&E2, k)/5e3; */ } Double_t MParticle::Planck0(Double_t *x, Double_t *k) { // // Planck, per unit volume, per unit energy // // constants (see below) moved out of function // const Double_t E = x[0]; // [GeV] const Double_t z = k ? k[0] : 0; const Double_t T = 2.96*(z+1); // [K] const Double_t e = 1.602176462e-19; // [C] const Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] const Double_t EkT = E/kB/T; /* //Double_t c = 299792458; // [m/s] //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] //Double_t hc = h*c; // [GeVm] Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] */ return E*E / (exp(EkT)-1.); // [GeV^2] } MParticle::MParticle(ParticleType_t t, const char *name, const char *title) : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0), fX(0) { // // default constructor // } void MParticle::InitRandom() { fPhi = gRandom->Uniform(TMath::Pi()*2); fPsi = gRandom->Uniform(TMath::Pi()*2); } void MParticle::SetNewDirection(Double_t theta, Double_t phi) { static TMatrixD B(3, 3); const Double_t ct = cos(fTheta); const Double_t st = sin(fTheta); const Double_t cp = cos(fPsi); const Double_t sp = sin(fPsi); /* TRotation B( ct*cp, ct*sp, -st, -sp, cp, 0, st*cp, st*sp, ct); */ // first row B(0, 0) = ct*cp; B(1, 0) = ct*sp; B(2, 0) = -st; // second row B(0, 1) = -sp; B(1, 1) = cp; B(2, 1) = 0; // third row B(0, 2) = st*cp; B(1, 2) = st*sp; B(2, 2) = ct; // ------------------------------ static TVectorD r(3); const Double_t sint = sin(theta); /* TVector3 r(sint*cos(phi), sint*sin(phi), cos(theta)); */ r(0) = sint*cos(phi); r(1) = sint*sin(phi); r(2) = cos(theta); // ------------------------------ r *= B; fTheta = asin(sqrt(r(0)*r(0)+r(1)*r(1))); // Numerically bad: acos(r(2)); fPsi = atan2(r(1), r(0)); /* if (fTheta*2 > TMath::Pi()) fTheta = fabs(fTheta-TMath::Pi()); */ } Bool_t MParticle::SetNewPosition(Double_t dr) { Bool_t rc=kTRUE; const Double_t st = sin(fTheta); TVector3 d(st*cos(fPsi), st*sin(fPsi), cos(fTheta)); // ------------------------------ const Double_t R = RofZ(&fZ); const Double_t dx = R - dr*d.z(); if (dx < 0) { dr = R/d.z(); // R>0 --> x(2)>0 rc = kFALSE; } fX += dr*(1.-d(2)); d *= dr; // ------------------------------ TVector3 r(fR*cos(fPhi), fR*sin(fPhi), R); r -= d; // ------------------------------ if (fR!=0) fPhi = atan2(r.y(), r.x()); fR = sqrt(r.x()*r.x()+r.y()*r.y()); fZ = ZofR(&r(2)); return rc; } Bool_t MParticle::SetNewPosition() { Double_t r = gRandom->Exp(GetInteractionLength()); return SetNewPosition(r); }