Index: trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1356) @@ -1,3 +1,22 @@ -*-*- END -*-*- + 2002/06/12: Thomas Bretz + + * MElectron.[h,cc]: + - added a primitive theta dependancy to DoInvCompton + - added DrawInteractionLength + + * MParticle.[h,cc]: + - added RofZ and ZofR as static functions + + * MPhoton.[h,cc]: + - added DrawInteractionLength + - fixed the bug causing the 'strange factor': sigma_gg needs + the sqrt of s. + + * energyloss.C: + - added + + + 2002/06/10: Thomas Bretz Index: trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1356) @@ -33,4 +33,8 @@ #include +#include +#include +#include +#include #include "MPhoton.h" @@ -286,36 +290,51 @@ #include -MPhoton *MElectron::DoInvCompton() -{ +MPhoton *MElectron::DoInvCompton(Double_t theta) +{ + static TRandom rand(0); + Double_t E0 = 511e-6; //[GeV] Double_t epsilon; Double_t e = GetEnergyLoss(&epsilon); + + // er: photon energy before interaction, rest frame + // e: photon energy after interaction, lab + + Double_t gamma = fEnergy/E0; + Double_t beta = sqrt(1.-1./(gamma*gamma)); + //Double_t gammabeta = sqrt(gamma*gamma-1); + + Double_t f = fEnergy/e; + + Double_t t; + Double_t arg; + do + { + t = rand.Uniform(TMath::Pi()*2); + Double_t er = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame + arg = (f - E0/er - 1)/(f*beta+1); + cout << "~" << flush; + + } while (arg<-1 || arg>1); + + Double_t theta1s = acos(arg); + Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta))); + + Double_t thetastar = thetas-theta1s; + + Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta))); + /* - Double_t gamma = fEnergy/E0; - Double_t beta = sqrt(1-1/(gamma*gamma)); - Double_t bega = sqrt(gamma*gamma-1); - - Double_t theta = TMath::Pi()/4; - - // er: photon energy before interaction rest frame - // e: photon energy after interaction lab - Double_t er = gamma*epsilon*(1-beta*cos(theta)); // photon energy rest frame - Double_t r1 = fEnergy/e; - Double_t r2 = E0/er; - Double_t ctheta = (r1-r2-1)/(beta*r1-1); - Double_t tg = sqrt(1-ctheta*ctheta)/gamma/(beta+ctheta); - */ + cout << "(" << theta1/TMath::Pi()*180 << ","; + cout << theta1s/TMath::Pi()*180<< ","; + cout << arg << ")" << flush; + */ + fEnergy -= e; MPhoton &p = *new MPhoton(e, fZ); p = *this; - - /* - cout << "(" << atan(tg)*180/TMath::Pi() << ")" << flush; - - static TRandom rand(0); - p->SetNewDirection(atan(tg), rand.Uniform(TMath::Pi()*2)); - */ + p.SetNewDirection(theta1, rand.Uniform(TMath::Pi()*2)); // MISSING: Electron angle @@ -327,2 +346,32 @@ return &p; } + +void MElectron::DrawInteractionLength(Double_t z) +{ + 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 MElectron::DrawInteractionLength() const +{ + DrawInteractionLength(fZ); +} Index: trunk/WuerzburgSoft/Thomas/mphys/MElectron.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1356) @@ -40,5 +40,10 @@ Double_t GetEnergyLoss(Double_t *ep) const; - MPhoton *DoInvCompton(); + MPhoton *DoInvCompton(Double_t theta); + + // ---------------------------------------------------------------- + + static void DrawInteractionLength(Double_t z); + void DrawInteractionLength() const; ClassDef(MElectron, 1) Index: trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1356) @@ -25,5 +25,5 @@ */ -Double_t ZofR(Double_t *x, Double_t *k=NULL) +Double_t MParticle::ZofR(Double_t *x, Double_t *k) { const Double_t c = 299792458; // [m/s] @@ -39,5 +39,5 @@ } -Double_t RofZ(Double_t *x, Double_t *k=NULL) +Double_t MParticle::RofZ(Double_t *x, Double_t *k) { Double_t z1 = x[0] + 1; Index: trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1356) @@ -16,5 +16,5 @@ protected: Double_t fEnergy; // [GeV] Energy - Double_t fBeta; // [1] beta +// Double_t fBeta; // [1] beta Double_t fZ; // [1] Red shift @@ -38,4 +38,7 @@ } + 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); @@ -44,5 +47,5 @@ void SetEnergy(Double_t e) { fEnergy = e; } - void SetBeta(Double_t b) { fBeta = b; } + // void SetBeta(Double_t b) { fBeta = b; } void SetZ(Double_t z) { fZ = z; } Index: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1356) @@ -33,4 +33,7 @@ #include +#include +#include +#include ClassImp(MPhoton); @@ -43,12 +46,12 @@ // constants moved out of function, see below // - 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; + 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; /* @@ -66,27 +69,22 @@ Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL) { - Double_t s = x[0]; // omega: CM mass - - Double_t E0 = 511e-6; // [GeV] - Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2 - - Double_t m = E0/s; - - Double_t m2 = m*m; - Double_t beta = sqrt(1.-m2); - Double_t beta2 = 1.-m2; - - Double_t p1 = r0*r0*TMath::Pi()/2; - - // ----- Extreme Relativistic ----- + const Double_t m2 = x[0]; // m2: (E0/sqrt(s))^2 + + const Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2 + + const Double_t beta2 = 1.-m2; + const Double_t beta = sqrt(beta2); + + const Double_t p1 = r0*r0*TMath::Pi()/2; + + // ----- Extreme Relativistic ------- // return p1*2 * m*m*m* (log(2./m)-1); - // -------------------------------- - - Double_t p2 = m2; - Double_t p3 = 3.-beta2*beta2; - Double_t p4 = log((1.+beta)/(1.-beta)); - Double_t p5 = beta*2*(1.+m2); - - Double_t sigma = p1*p2*(p3*p4-p5); // [m^2] + // ---------------------------------- + + const Double_t p2 = 3.-beta2*beta2; + const Double_t p3 = log((1.+beta)/(1.-beta)); + const Double_t p4 = beta*2*(1.+m2); + + const Double_t sigma = p1*m2*(p2*p3-p4); // [m^2] return sigma; @@ -95,17 +93,16 @@ Double_t MPhoton::Int1(Double_t *x, Double_t *k=NULL) { - Double_t costheta = x[0]; - - Double_t Eg = k[0]; - Double_t Ep = k[1]; - - Double_t E0 = 511e-6; // [GeV] - - Double_t s = Eg*Ep/E0*(1.-costheta)*2; - - if (s1) // Why is this necessary??? return 0; - Double_t sigma = Sigma_gg(&s); // [m^2] + const Double_t sigma = Sigma_gg(&s); // [m^2] return sigma/2 * (1.-costheta); // [m^2] @@ -114,23 +111,22 @@ Double_t MPhoton::Int2(Double_t *x, Double_t *k) { - Double_t E0 = 511e-6; // [GeV] + const Double_t E0 = 511e-6; // [GeV] Double_t Ep = x[0]; - Double_t Eg = k[0]; - Double_t z = k[1]; + const Double_t Eg = k[0]; + Double_t val[2] = { Eg, Ep }; - Double_t from = -1.0; - Double_t to = 1.-E0*E0/2./Eg/Ep; // Originally Was: 1. + const Double_t from = -1.0; + const Double_t to = 1.-E0*E0/(2.*Eg*Ep); // Originally Was: 1. TF1 f("int1", Int1, from, to, 2); - Double_t int1 = f.Integral(from, to, val); // [m^2] - Double_t planck = Planck(&Ep, &z); // [GeV^2] - - Double_t res = planck * int1; - - res *= Eg/E0*1e-9; // FIXME!!!!!!!!!! WHICH FACTOR IS THIS???? + + const Double_t int1 = f.Integral(from, to, val); // [m^2] + const Double_t planck = Planck(&Ep, &z); // [GeV^2] + + const Double_t res = planck * int1; return res; // [GeV^2 m^2] @@ -156,6 +152,6 @@ Double_t val[2] = { Eg, z }; - Double_t lolim = E0*E0/Eg; - Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg); + Double_t lolim = E0*E0 > 1e-8 ? E0*E0/Eg : 1e-8/Eg; + Double_t inf = 3e-11; TF1 f("int2", Int2, lolim, inf, 2); @@ -192,2 +188,32 @@ } +void MPhoton::DrawInteractionLength(Double_t z) +{ + if (!gPad) + new TCanvas("ILPhoton", "Mean Interaction Length Photon"); + else + gPad->GetVirtCanvas()->cd(4); + + TF1 f1("length", 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 MPhoton::DrawInteractionLength() const +{ + DrawInteractionLength(fZ); +} Index: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1355) +++ trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1356) @@ -17,4 +17,6 @@ void operator=(MParticle &p) { MParticle::operator=(p); } + // ---------------------------------------------------------------- + static Double_t Planck(Double_t *x, Double_t *k=NULL); static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL); @@ -26,4 +28,9 @@ Double_t GetInteractionLength() const; + // ---------------------------------------------------------------- + + static void DrawInteractionLength(Double_t z); + void DrawInteractionLength() const; + ClassDef(MPhoton, 1) }; Index: trunk/WuerzburgSoft/Thomas/mphys/energyloss.C =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/energyloss.C (revision 1356) +++ trunk/WuerzburgSoft/Thomas/mphys/energyloss.C (revision 1356) @@ -0,0 +1,489 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "mbase/MParContainer.h" +#include "mphys/MPhoton.h" +#include "mphys/MElectron.h" +#include "mphys/MPairProduction.h" +#include "mhist/MBinning.h" +#include "mhist/MH.h" + +// 2.96 +// 2.87 +// 2.73 + +// ================================================================ +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. +} + +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 ? MParticle::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; +} + +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 = MParticle::RofZ(&k[2]) - r; + Double_t z = k[2]>=0 ? MParticle::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 = MParticle::RofZ(&z) - l; + + cout << " " << i << ". R=" << MParticle::RofZ(&z) << " l=" << l << " z=" << MParticle::ZofR(&r) << endl; + + z = r>=0 ? MParticle::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 DrawRZ() +{ + new TCanvas("RZ", "r and z"); + + TF1 f1("ZofR", MParticle::ZofR, 0, 4.5e5, 0); + TF1 f2("RofZ", MParticle::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(); +} + +// ------------------------------------------------------------------- + +void energyloss() +{ + EnergyLossRate(); + + MH::MakeDefCanvas(); + + DrawRZ(); + + return; + +}