Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1356) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1357) @@ -8,4 +8,6 @@ * MParticle.[h,cc]: - added RofZ and ZofR as static functions + - changed from MParContainer to TObject + - removed unnecessary or commented inline functions * MPhoton.[h,cc]: @@ -16,4 +18,15 @@ * energyloss.C: - added + + * MPairProduction.[h,cc]: + - simplified formula + - changed arguments to Energy of Photon and interaction angle + + * phys.C: + - added correct photon energy for pair production + - added correct angle for pair production + - changed output to E^2*Counts + - changed gamma production from E^-2 to Uniform by weighting + - removed all unecessary functions and code (s.energyloss.C) Index: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1356) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1357) @@ -71,22 +71,69 @@ // -------------------------------------------------------------------------- -Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, const Double_t theta, TList *list) +Bool_t MPairProduction::Process(MParticle *gamma, const Double_t Ep, const Double_t theta, TList *list) { // // gamma: primary particle from source // phot: infrared photon from background. (angle = interaction angle) - // + // Ep: Energy photon [GeV] + // theta: interaction angle [rad] + // + + const Double_t E0 = 511e-6; // [GeV] - - //const Double_t theta = phot->GetAngle(); // [2pi] - const Double_t Ep = phot->GetEnergy(); // [GeV] const Double_t Eg = gamma->GetEnergy(); // [GeV] const Double_t ctheta = cos(theta); - const Double_t s = (Ep*Eg/(2*E0*E0))*(1-ctheta); //[1] - if (s<1) + const Double_t s = Ep*Eg*2*(1-ctheta); //[1] + const Double_t S = s/(E0*E0*4); //[1] + if (S<1) return kFALSE; + const Double_t stheta = sin(theta); + + const Double_t sqrbetah = s/((Eg+Ep)*(Eg+Ep)) + 1; + const Double_t sqrbetae = 1.-1./S; + + const Double_t GammaH = (Eg+Ep)/sqrt(s); + + Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah))); + + fAngle->SetParameter(0, sqrt(sqrbetae)); + + const Double_t alpha = psi-acos(fAngle->GetRandom()); + + Double_t salpha = sin(alpha); + Double_t calpha = cos(alpha); + + const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi) + + Double_t bb = sqrt(sqrbetah/sqrbetae); + + Double_t s1 = calpha/GammaH; + Double_t s2 = tphi*s1 - salpha - bb; + + Double_t tan1 = ((salpha+bb)*tphi+s1)/s2; + Double_t tan2 = ((salpha-bb)*tphi+s1)/s2; + + const Double_t E = (Eg+Ep)/2;; + const Double_t f = sqrt(sqrbetah*sqrbetae)*salpha; + + cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush; + + MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ()); + MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ()); + p0 = *gamma; // Set Position and direction + p1 = *gamma; // Set Position and direction + + static TRandom rand(0); + Double_t rnd = rand.Uniform(TMath::Pi() * 2); + p0.SetNewDirection(atan(+tan1), rnd); + p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi()); + + list->Add(&p0); + list->Add(&p1); + + /* const Double_t Epg = Ep/Eg; // 1+Epg ~ 1 const Double_t Egp = Eg/Ep; // 1+Egp ~ Egp @@ -120,4 +167,17 @@ const Double_t E = sqrt(s)*E0/gamma1; const Double_t dE = E*Bcosd; + + const Double_t E1 = E0/(E+dE); + const Double_t E2 = E0/(E-dE); + + const Double_t beta1 = sqrt(1.-E1*E1); + const Double_t beta2 = sqrt(1.-E2*E2); + + const Double_t Bscp = Bsind*cphi; + const Double_t spg = sphi/gamma1; + const Double_t cpg = cphi/gamma1; + + const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp)); + const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp)); MElectron &p0 = *new MElectron(E+dE, gamma->GetZ()); @@ -126,19 +186,6 @@ p1 = *gamma; // Set Position and direction - const Double_t E1 = E0/(E+dE); - const Double_t E2 = E0/(E-dE); - - const Double_t beta1 = sqrt(1.-E1*E1); - const Double_t beta2 = sqrt(1.-E2*E2); - p0.SetBeta(beta1); p1.SetBeta(beta2); - - const Double_t Bscp = Bsind*cphi; - const Double_t spg = sphi/gamma1; - const Double_t cpg = cphi/gamma1; - - const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp)); - const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp)); static TRandom rand(0); @@ -149,5 +196,5 @@ list->Add(&p0); list->Add(&p1); - + */ return kTRUE; } Index: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 1356) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 1357) @@ -19,5 +19,5 @@ virtual ~MPairProduction(); - Bool_t Process(MParticle *g, MParticle *p, const Double_t theta, TList *l); + Bool_t Process(MParticle *g, const Double_t Ep, const Double_t theta, TList *l); ClassDef(MPairProduction, 0) // Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1356) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1357) @@ -24,4 +24,6 @@ ************************ */ +#include +#include Double_t MParticle::ZofR(Double_t *x, Double_t *k) @@ -65,7 +67,9 @@ // set the name and title of this object // - - fName = name ? name : "MParticle"; - fTitle = title ? title : "Container for a particle"; + + /* + fName = name ? name : "MParticle"; + fTitle = title ? title : "Container for a particle"; + */ } @@ -104,7 +108,4 @@ fTheta = fabs(fTheta-TMath::Pi()); } - -#include -#include Bool_t MParticle::SetNewPosition(Double_t dr) Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1356) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1357) @@ -2,9 +2,15 @@ #define MARS_MParticle -#ifndef MARS_MParContainer -#include "MParContainer.h" +#ifndef ROOT_TObject +#include #endif -class MParticle : public MParContainer +/* + #ifndef MARS_MParContainer + #include "MParContainer.h" + #endif + */ + +class MParticle : public TObject { public: @@ -16,5 +22,4 @@ protected: Double_t fEnergy; // [GeV] Energy -// Double_t fBeta; // [1] beta Double_t fZ; // [1] Red shift @@ -27,5 +32,7 @@ public: MParticle(ParticleType_t t, const char *name=NULL, const char *title=NULL); - // ~MParticle(); + + static Double_t ZofR(Double_t *x, Double_t *k=NULL); + static Double_t RofZ(Double_t *x, Double_t *k=NULL); void operator=(MParticle &p) @@ -38,22 +45,8 @@ } - static Double_t MParticle::ZofR(Double_t *x, Double_t *k=NULL); - static Double_t MParticle::RofZ(Double_t *x, Double_t *k=NULL); - - // void Fill(const MParContainer *inp); - - // void Draw(Option_t *opt=NULL); virtual Double_t GetInteractionLength() const = 0; void SetEnergy(Double_t e) { fEnergy = e; } - // void SetBeta(Double_t b) { fBeta = b; } - void SetZ(Double_t z) { fZ = z; } - /* - void SetPhi(Double_t p) { fPhi = p; } - void SetR(Double_t r) { fR = r; } - void SetTheta(Double_t t) { fTheta = t; } - void SetPsi(Double_t p) { fPsi = p; } - */ void SetNewDirection(Double_t theta, Double_t phi); Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1356) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1357) @@ -1,4 +1,4 @@ -//#ifndef __CINT__ #include + #include @@ -8,9 +8,8 @@ #include #include +#include #include -#include #include -#include -//#endif + #include "mbase/MParContainer.h" #include "mphys/MPhoton.h" @@ -19,629 +18,8 @@ #include "mhist/MBinning.h" #include "mhist/MH.h" + // 2.96 // 2.87 // 2.73 - -/*********************** - * - * calc r from z: - * - * R = 2*c/H0 *(1+z - sqrt(1+z)) - * - * z = -0.5 * (3 + R' +/- sqrt(R'+1)) R' = R*H0/c - * - * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] - * - ************************ - */ - -Double_t ZofR(Double_t *x, Double_t *k=NULL) -{ - Double_t c = 299792458; // [m/s] - Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] - - Double_t ly = 3600.*24.*365.*c; // [m/ly] - Double_t pc = 1./3.258; // [pc/ly] - Double_t r = x[0] /pc*ly*1e3; // [m] - - Double_t R = r*H0/c; // [1] - - return (R+1+sqrt(R*2+1))/2 - 1; -} - -Double_t RofZ(Double_t *x, Double_t *k=NULL) -{ - Double_t z1 = x[0] + 1; - - Double_t c = 299792458; // [m/s] - Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] - - Double_t ly = 3600.*24.*365.*c; // [m/ly] - Double_t pc = 1./3.258; // [pc/ly] - - Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m] - - return R * pc/ly/1e3; // [kpc] -} - -Double_t AngleDistrib(Double_t *x, Double_t *k) -{ - - Double_t c = x[0]; // cos(alpha) - Double_t b = k[0]; // sqrt(1-1/s) - - Double_t b2 = b*b; - Double_t b4 = b2*b2; - - Double_t c2 = c*c; - Double_t c4 = c2*c2; - - Double_t u = 1. - b4*c4 +2.*b2*(1.-b2)*(1.-c2); - Double_t d = 1.-b2*c2; - - return u/(d*d); -} - - -// ================================================================ -Double_t DiSum(Double_t *x, Double_t *k=NULL) -{ - Double_t t = x[0]; - - Double_t disum = t; - Double_t add = 0; - - Double_t eps = fabs(t*1e-1); - - Int_t n = 2; - Double_t pow = t*t; // t^2 - - do - { - add = pow/n/n; - - pow *= t; // pow = t^n - n++; - - disum += add; - - } while (fabs(add)>eps); - - return disum; -} - -Double_t F(Double_t *x, Double_t *k=NULL) -{ - Double_t o = x[0]; - Double_t s = -2.*o; - -// if (o<1e-10) -// return 2.125; //-3./8.; - - return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8. -} - -Double_t FlimLi(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam) mit b-->1 (Maple) -{ - Double_t w = x[0]; - - Double_t s = -w*2*(1+1); // -2*omega*(1+beta) - Double_t li2 = MElectron::Li2(&s); - - return fabs(li2); -} - -#include -void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[], - Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[])) -{ - /* - * --------------------------------------------------------- - * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method - * --------------------------------------------------------- - */ - const Double_t a2 = 0.2; - const Double_t a3 = 0.3; - const Double_t a4 = 0.6; - const Double_t a5 = 1.0; - const Double_t a6 = 0.875; - const Double_t b21 = 0.2; - const Double_t b31 = 3.0/40.0; - const Double_t b32 = 9.0/40.0; - const Double_t b41 = 0.3; - const Double_t b42 = -0.9; - const Double_t b43 = 1.2; - const Double_t b51 = -11.0/54.0; - const Double_t b52 = 2.5; - const Double_t b53 = -70.0/27.0; - const Double_t b54 = 35.0/27.0; - const Double_t b61 = 1631.0/55296.0; - const Double_t b62 = 175.0/512.0; - const Double_t b63 = 575.0/13824.0; - const Double_t b64 = 44275.0/110592.0; - const Double_t b65 = 253.0/4096.0; - const Double_t c1 = 37.0/378.0; - const Double_t c3 = 250.0/621.0; - const Double_t c4 = 125.0/594.0; - const Double_t c6 = 512.0/1771.0; - const Double_t dc5 = -277.00/14336.0; - - const Double_t dc1 = c1-2825.0/27648.0; - const Double_t dc3 = c3-18575.0/48384.0; - const Double_t dc4 = c4-13525.0/55296.0; - const Double_t dc6 = c6-0.25; - - Double_t ak2[n]; - Double_t ak3[n]; - Double_t ak4[n]; - Double_t ak5[n]; - Double_t ak6[n]; - Double_t ytemp[n]; - - for (int i=0; ib ? a : b; -} - -void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps, - Double_t yscal[], Double_t *hdid, Double_t *hnext, - void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs -{ - /* - * --------------------------------------------------------- - * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method - * --------------------------------------------------------- - */ - const Double_t SAFETY = 0.9; - const Double_t PGROW = -0.2; - const Double_t PSHRNK = -0.25; - const Double_t ERRCON = 1.89e-4; - - Double_t yerr[n]; - Double_t ytemp[n]; - - Double_t h = htry; - Double_t errmax; - while (1) - { - rkck(y, dydx, n, *x, h, ytemp, yerr, derivs); - - errmax=0.0; - - for (int i=0; i= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h)); - - Double_t xnew= (*x) + h; - - if (xnew != *x) - continue; - - cout << "stepsize underflow in rkqs" << endl; - break; - } - - *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h; - - *x += (*hdid=h); - - for (int i=0; i=0 ? ZofR(&R) : 0; - - Double_t e = 1.602176462e-19; // [C] - Double_t T = 2.96*(z+1); // [K] - Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] - Double_t kBT = kB*T; // [GeV] - Double_t alpha = 1./137; // [|] - Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] - Double_t E0 = 511e-6; // [GeV] - - - Double_t k = TMath::Pi() * alpha * kBT; - - Double_t ln = 4.*kBT/E0/E0; - - return -1./3/h* k*k * (log(ln*E)-1.9805); -} - -void dEdt(Double_t t, Double_t E[], Double_t dedt[]) -{ - dedt[0] = dEdt(E[0], t); - //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl; -} - -Int_t GetSeed() -{ - return (Int_t)fmod(TStopwatch::GetRealTime()*10, 65535); -} - -void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="") -{ - Double_t t = 0; - Double_t E[1] = { E0 }; - Double_t yscal[1] = { 1 }; - - Double_t dedt[1]; - - Double_t eps;// = 1e5; - Double_t step = 5e6; - - Double_t hdid; - Double_t hnext; - - cout << "Start: " << dedt[0] << endl; - - Int_t n = 15000; - Double_t tres[n]; - Double_t Eres[n]; - int i; - for (i=0; i1.5e14 || E[0]<5e6) - break; -// cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl; - } - - cout << i << endl; - - TGraph grp(iDraw(opt); - -} - -Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL) -{ - Double_t t = x[0]; - Double_t E = k[0]; - Double_t t0 = k[1]; - - Double_t c = 299792458; // [m/s] - Double_t ly = 3600.*24.*365.*c; // [m/ly] - Double_t pc = 1./3.258; // [pc/ly] - - Double_t r = t * c / ly * pc / 1000; // [kpc] - Double_t R = RofZ(&k[2]) - r; - Double_t z = k[2]>=0 ? ZofR(&R) : 0; - - Double_t T = 2.96*(z+1); // [K] - // Double_t alpha = 1./137; // [1] - Double_t sigma = 6.653e-29; // [m^2] - Double_t E0 = 511e-6; // [GeV] - Double_t e = 1.602176462e-19; // [C] - Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] - Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] - Double_t hc = h*c; // [GeVm] - Double_t pi = TMath::Pi(); // [1] - Double_t khc = pi*kB/hc; // [1 / K m] - Double_t a = 8./15 * pi * khc*khc*khc*khc * hc; // [Gev / K^4 / m^3] - Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s] - - Double_t ret = 1./(konst*(t-t0) + 1./E); - - return ret; -} - -void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="") -{ - // 8.7 - - Double_t val[] = { E0, t0, z }; - Double_t t = 1.5e14; - while (EnergyLossRateLoLim(&t, val)<1e4) - t -= 0.01e14; - - TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2); - f1->SetParameter(0, E0); - f1->SetParameter(1, t0); - - f1->Draw(opt); - f1->SetBit(kCanDelete); -} - -// -// (3) Energy loss rate of electrons and 'high energy particle' -// -Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL) -{ - Double_t ly = 3600.*24.*365.; // [s/ly] - Double_t pc = 1./3.258; // [pc/ly] - - TGraph *graph = new TGraph; - graph->SetPoint(0, 0, E); - graph->SetMaximum(E*3); // *MENU* - - cout << "------ " << endl; - - static TRandom rand; - - Double_t x = 0; - for (int i=1; i<10; i++) - { - Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z)); - - if (z>=0) - { - Double_t r = RofZ(&z) - l; - - cout << " " << i << ". R=" << RofZ(&z) << " l=" << l << " z=" << ZofR(&r) << endl; - - z = r>=0 ? ZofR(&r) : 0; - - if (z==0) - cout << "z<0" << endl; - } - - x += l; - - Double_t t = x/pc*ly*1000; - graph->SetPoint(i*2-1, t, E); - - Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z); - - E -= e1; - - if (hist) - hist->Fill(e1); - - cout << " Ep=" << e1 << flush; - - graph->SetPoint(i*2, t, E); - } - - graph->SetMinimum(E/3); // *MENU* - graph->Draw(opt); - graph->SetBit(kCanDelete); - - //if (E<31500) - cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl; - - return E; -} - -void EnergyLossRate() -{ - if (gPad) - gPad->Clear(); - - Double_t E = 1.5e9; // [GeV] - Double_t z = 0.03; - - MBinning bins; - bins.SetEdgesLog(18, 0.1, 1e9); - - TH1D hist; - hist.SetName("Phot"); - hist.SetTitle("Photons from inverse Compton"); - MH::SetBinning(&hist, &bins); - - cout << "Working..." << flush; - - for (int i=0; i<50; i++) - { - cout << i << "." << flush; - DrawDevelopment(E, z, i?"L":"AL", &hist); - } - - //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen - DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen - DrawDevelopmentHiLim(E, z, "L"); - gPad->SetLogy(); - - new TCanvas("Photons", "Photons created in inverse Compton"); - hist.DrawCopy(); - gPad->SetLogx(); - - cout << "...done." << endl; -} - -void Scatter() -{ - Double_t E = 1.5e9; // [GeV] - - Int_t nbins = 100; - Double_t lo = E/1e6; - Double_t up = E*10; - - Double_t binsize = log10(up/lo)/nbins; - Double_t *scale=new Double_t[nbins+1]; - for (int i=0; iDraw(); - - TStopwatch clock; - clock.Start(); - - for (int i=0; i<1000; i++) - { - Double_t e1 = epsilon.GetRandom(); - - hist->Fill(e1); - } - - clock.Stop(); - clock.Print(); - - epsilon.DrawCopy("same"); - gPad->SetLogx(); - gPad->SetLogy(); - */ -} - -void DrawILPhoton(Double_t z=0) -{ - if (!gPad) - new TCanvas("IL_Photon", "Mean Interaction Length Photon"); - else - gPad->GetVirtCanvas()->cd(4); - - TF1 f1("length", MPhoton::InteractionLength, 1e4, 1e11, 1); - f1.SetParameter(0, z); - - gPad->SetLogx(); - gPad->SetLogy(); - gPad->SetGrid(); - f1.SetMaximum(1e5); - f1.SetLineWidth(1); - - TH1 *h=f1.DrawCopy()->GetHistogram(); - - h->SetTitle("Mean Interaction Length (Photon)"); - h->SetXTitle("E [GeV]"); - h->SetYTitle("x [kpc]"); - - gPad->Modified(); - gPad->Update(); -} - -void DrawILElectron(Double_t z=0) -{ - if (!gPad) - new TCanvas("IL_Electron", "Mean Interaction Length Electron"); - else - gPad->GetVirtCanvas()->cd(3); - - TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0); - f2.SetParameter(0, z); - - gPad->SetLogx(); - gPad->SetLogy(); - gPad->SetGrid(); - f2.SetLineWidth(1); - - TH1 *h=f2.DrawCopy()->GetHistogram(); - - h->SetTitle("Mean Interaction Length (Electron)"); - h->SetXTitle("E [GeV]"); - h->SetYTitle("x [kpc]"); - - gPad->Modified(); - gPad->Update(); -} - -void DrawRZ() -{ - new TCanvas("RZ", "r and z"); - - TF1 f1("ZofR", ZofR, 0, 4.5e5, 0); - TF1 f2("RofZ", RofZ, 0, 5, 0); - - gPad->Divide(2,2); - - gPad->GetVirtCanvas()->cd(1); - TH1 *h = f1.DrawCopy()->GetHistogram(); - - h->SetTitle("z(r)"); - h->SetXTitle("r [kpc]"); - h->SetYTitle("z"); - - gPad->Modified(); - gPad->Update(); - - gPad->GetVirtCanvas()->cd(2); - h = f2.DrawCopy()->GetHistogram(); - - h->SetTitle("r(z)"); - h->SetXTitle("z"); - h->SetYTitle("r [kpc]"); - - gPad->Modified(); - gPad->Update(); -} - Double_t PrimSpect(Double_t *x, Double_t *k) { @@ -649,65 +27,68 @@ } -/* -Double_t Planck(Double_t *x, Double_t *k=NULL) -{ - Double_t E = x[0]; // [GeV] - Double_t z = k ? k[0] : 0; - - Double_t T = 2.96*(z+1); // [K] - Double_t e = 1.602176462e-19; // [C] - Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] - - Double_t EkT = E/kB/T; - - return E*E*E / (exp(EkT)-1.); // [GeV^2] -} -*/ - -Double_t Planck(Double_t *x, Double_t *k=NULL) -{ - // - // Planck, per unit volume, per unit energy - // - // constants moved out of function, see below - // - - // - // This is to get a correct random value in a TF1! - // - Double_t E = pow(10, x[0]); // [GeV] - Double_t z = k ? k[0] : 0; - - Double_t T = 2.96*(z+1); // [K] - Double_t e = 1.602176462e-19; // [C] - Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] - - Double_t EkT = E/kB/T; - - return E*E / (exp(EkT)-1.); // [GeV^2] +Double_t PhotonSpect(Double_t *x, Double_t *k=NULL) +{ + Double_t Ep = pow(10, x[0]); + Double_t res = MPhoton::Int2(&Ep, k); + + return res*1e55; //65/k[0]; + // return MPhoton::Planck(&Ep, &k[1]); +} + +Double_t Sbar_sigmas(Double_t *x, Double_t *k) +{ + Double_t sbar = pow(10, x[0]); + + Double_t s = 1./(sbar*4); + + Double_t sigma = MPhoton::Sigma_gg(&s); + + return sigma*sbar*1e28; +} + +Double_t RandomTheta(Double_t Eg, Double_t Ep) +{ + Double_t E0 = 511e-6; // [GeV] + + Double_t f = Eg/E0*Ep/E0; + + if (f<1) + return 0; + + TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0); + + Double_t sbar = pow(10, func.GetRandom()); + Double_t theta = acos(1.-sbar*2/f); + + return theta; } void DoIt() { - Double_t rcygnusx3 = 5000; //100; //30/3.258; // [~10kpc] - - Double_t startz = ZofR(&rcygnusx3); + Double_t R = 100; // [kpc] + Double_t startz = MParticle::ZofR(&R); + + cout << "R = " << R << "kpc" << endl; + cout << "Z = " << startz << endl; static TRandom rand(0); MPairProduction pair; - //Double_t nphot = 50; - Double_t runtime = 15*60; ///*18*60*/60; // [s] - - Double_t lo = 2e5; - Double_t hi = 2e10; + Double_t runtime = 30*60; // [s] + + Double_t lo = 1e4; + Double_t hi = 1e10; Double_t alpha = -2; - Double_t nbins = 4*log10(hi/lo); - - TF1 phot("PhotonSpectrum", Planck, -15, -10.5, 1); + Double_t nbins = 24; //4*log10(hi/lo); + TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV - - src.SetParameter(0, alpha); + src.SetParameter(0, 0); + + TH1D hist; + TH1D histsrc; + + hist.SetMinimum(lo*lo * pow(lo, alpha)/10); + histsrc.SetMinimum(lo*lo * pow(lo, alpha)/10); TList listg; @@ -726,6 +107,6 @@ MBinning binspolx; MBinning binspoly; - binspolx.SetEdges(32, -180, 180); - binspoly.SetEdgesLog(20, 1e-9, 1e-5); + binspolx.SetEdges(16, -180, 180); + binspoly.SetEdgesLog(20, 1e-13, 1e-4); MH::SetBinning(&position, &binspolx, &binspoly); MH::SetBinning(&angle, &binspolx, &binspoly); @@ -737,7 +118,4 @@ histpos.SetTitle("Particle Position"); position.SetTitle("Particle Position"); - - TH1D hist; - TH1D histsrc; // @@ -761,14 +139,14 @@ gPad->SetLogy(); + histsrc.SetName("Spectrum"); + histsrc.SetXTitle("E [GeV]"); + histsrc.SetYTitle("E^{2}\\dotCounts"); + histsrc.GetXaxis()->SetLabelOffset(-0.015); + histsrc.GetXaxis()->SetTitleOffset(1.1); + histsrc.Draw("P"); hist.Draw("Psame"); histsrc.Draw("Csame"); hist.Draw("Csame"); - - histsrc.SetName("Spectrum"); - histsrc.SetXTitle("E [GeV]"); - histsrc.SetYTitle("E*Counts"); - histsrc.GetXaxis()->SetLabelOffset(-0.015); - histsrc.GetXaxis()->SetTitleOffset(1.1); ofstream fsrc("source.dat"); @@ -789,6 +167,8 @@ Double_t E = pow(10, src.GetRandom()); - - histsrc.Fill(E, E); + Double_t weight = pow(E, alpha); + + histsrc.Fill(E, E*E * weight); + fsrc << E << endl; @@ -800,5 +180,6 @@ while (1) { - cout << " |P:" << flush; + if (listg.GetSize()!=0) + cout << " |P:" << flush; TIter NextP(&listg); @@ -806,15 +187,17 @@ while ((p=(MPhoton*)NextP())) { + Double_t Eg = p->GetEnergy(); + if (!p->SetNewPosition()) { cout << "!" << flush; - hist.Fill(p->GetEnergy(), p->GetEnergy()); - position.Fill(p->GetPhi()*kRad2Deg, p->GetR()); - angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg); - histpos.Fill(p->GetR()); - hista.Fill(p->GetTheta()*kRad2Deg); - - fcas << p->GetEnergy() << endl; + hist.Fill(Eg, Eg*Eg*weight); + position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight); + angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight); + histpos.Fill(p->GetR(), weight); + hista.Fill(p->GetTheta()*kRad2Deg, weight); + + fcas << Eg << endl; delete listg.Remove(p); continue; @@ -824,13 +207,30 @@ // Sample phtoton from background and interaction angle // - MPhoton photon; - - phot.SetParameter(0, p->GetZ()); - photon.SetEnergy(pow(10, phot.GetRandom())); - Double_t theta = rand.Uniform(TMath::Pi() * 2); - - if (!pair.Process(p, &photon, theta, &liste)) - { - cout << "0" << flush;; + TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2); + phot.SetParameter(0, Eg); + phot.SetParameter(1, p->GetZ()); + if (phot.Integral(-log10(Eg)-8, -10.5)==0) + { + // ??????????????????? + hist.Fill(Eg, Eg*Eg*weight); + angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight); + hista.Fill(p->GetTheta()*kRad2Deg, weight); + delete listg.Remove(p); + cout << "z" << flush; + continue; + } + + Double_t Ep = pow(10, phot.GetRandom()); + + Double_t theta = RandomTheta(Eg, Ep); + if (theta==0) + { + cout << "t" << flush; + continue; + } + + if (!pair.Process(p, Ep, theta, &liste)) + { + cout << "0" << flush; continue; } @@ -843,5 +243,6 @@ break; - cout << " |E" << flush; + if (liste.GetSize()!=0) + cout << " |E" << flush; TIter Next(&liste); @@ -849,4 +250,7 @@ while ((e=(MElectron*)Next())) { + if (e->GetEnergy()DoInvCompton(); - if (e->GetEnergy()DoInvCompton(theta); + + if (fabs(p->GetTheta()*3437)<1 && // < 1min + p->GetEnergy()>lo) { - cout << "0" << flush; - delete p; - break; + cout << "." << flush; + listg.Add(p); } - - if (p->GetEnergy()GetTheta()*3437)<1 && // < 1min + e->GetEnergy()>lo*2) continue; - } - - listg.Add(p); - cout << "." << flush; + + cout << "0" << flush; + break; } } @@ -940,5 +349,4 @@ hista.DrawCopy(); gPad->SetLogx(); - } @@ -946,115 +354,40 @@ void phys() { + DoIt(); + /* - Double_t Eg = 1e5; - - Double_t E0 = 511e-6; // [GeV] - Double_t lolim = E0*E0/Eg; - Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg); - - // 1e5: 5e-13, 4e-12 - // 1e6: 5e-14, 2e-12 - // 1e8: 5e-16, 1e-12 - // 1e10: 5e-18, 1e-12 - - // 1e5: 2.6e-12, 4e-12 - // 1e6: 2.6e-13, 2e-12 - // 1e8: 2.6e-15, 1e-12 - // 1e10: 1.6e-17, 1e-12 - - TF1 f("int2", MPhoton::Int2, lolim, inf, 2); - f.SetParameter(0, Eg); - f.SetParameter(1, 0); - - MH::MakeDefCanvas(); - gPad->SetLogx(); - f.DrawCopy(); - return; + Double_t E = 1e10; + TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2); + phot.SetParameter(0, E); + phot.SetParameter(1, 5); + phot.DrawCopy(); + return; + */ + + /* + Double_t Eg = 1e5; + + Double_t E0 = 511e-6; // [GeV] + Double_t lolim = E0*E0/Eg; + Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg); + + // 1e5: 5e-13, 4e-12 + // 1e6: 5e-14, 2e-12 + // 1e8: 5e-16, 1e-12 + // 1e10: 5e-18, 1e-12 + + // 1e5: 2.6e-12, 4e-12 + // 1e6: 2.6e-13, 2e-12 + // 1e8: 2.6e-15, 1e-12 + // 1e10: 1.6e-17, 1e-12 + + TF1 f("int2", MPhoton::Int2, lolim, inf, 2); + f.SetParameter(0, Eg); + f.SetParameter(1, 0); + + MH::MakeDefCanvas(); + gPad->SetLogx(); + f.DrawCopy(); */ -// MH::MakeDefCanvas(); -// DrawILPhoton(); - -// return ; - - DoIt(); -return; - EnergyLossRate(); - - MH::MakeDefCanvas(); - DrawILElectron(); - - DrawRZ(); - - return; - - /* - cout << "Starting..." << endl; - - MEvtLoop magic; - MParList plist; - MTaskList tlist; - plist.AddToList(&tlist); - - MHInput *input = new MHInput; - MHOutput *output = new MHOutput; - - plist.AddToList(input); - plist.AddToList(output); - - MGenerateInput geninp; - tlist.AddToList(&geninp); - - MFillH finput("MSinglePairInput", input); - tlist.AddToList(&finput); - - MSinglePairProduction prod; - tlist.AddToList(&prod); - - MFillH foutput("MSinglePairOutput", output); - tlist.AddToList(&foutput); - - magic.SetParList(&plist); - if (!magic.Eventloop(10000)) - return; - - input->Draw(); - output->Draw(); - - return; -*/ - - /* - MEvtLoop magic; - MParList plist; - MTaskList tlist; - plist.AddToList(&tlist); - - MHInput *input = new MHInput; - MHOutput *output = new MHOutput; - - plist.AddToList(input); - plist.AddToList(output); - - MGenerateInput geninp; - tlist.AddToList(&geninp); - - MFillH finput("MSinglePairInput", input); - tlist.AddToList(&finput); - - MSinglePairProduction prod; - tlist.AddToList(&prod); - - MFillH foutput("MSinglePairOutput", output); - tlist.AddToList(&foutput); - - magic.SetParList(&plist); - if (!magic.Eventloop(10000)) - return; - - input->Draw(); - output->Draw(); - - return; - */ -} - +} +