Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1364) @@ -1,3 +1,32 @@ -*-*- END -*-*- + 2002/06/14: Thomas Bretz + + * MElectron.[h,cc]: + - Moved Planck function to MParticle + - Added the DiSum function + - changed Li2 to use the DiSum function (speed reasons) + - changed the eps of all integrals to 1e-2 + - changed the p_e function to use the Compton function + + * MPhoton.[h,cc]: + - removed the planck function + - changed the eps of all integrals to 1e-2 + + * MParticle.[h,cc]: + - added the planck function + + * Makefile: + - removed unused source files + + * analp.C: + - did some small fixes + - changed the sclae of the radius- and phi-plots + + * phys.C: + - changed the integral eps to 1e-2 + + + + 2002/06/13: Thomas Bretz Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1364) @@ -42,62 +42,54 @@ ClassImp(MElectron); -Double_t MElectron::Planck(Double_t *x, Double_t *k=NULL) -{ - // - // Planck, per unit volume, per unit energy - // - // constants 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] -} - Double_t MElectron::Li(Double_t *x, Double_t *k) { const Double_t t = x[0]; - return log(1.-t)/t; } +Double_t DiSum(Double_t *x, Double_t *k=NULL) +{ + Double_t t = x[0]; + + const Double_t eps = fabs(t*1e-2); + + Double_t disum = t; + Double_t add = 0; + + 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 MElectron::Li2(Double_t *x, Double_t *k=NULL) { // // Dilog, Li2 - // - Double_t z = x[0]; + // ---------- + // + // Integral(0, 1) = konst; + // Double_t konst = 1./6*TMath::Pi()*TMath::Pi(); + // + // x[0]: z + // + const Double_t z = x[0]; + + if (fabs(z)<1) + return DiSum(x); TF1 IntLi("Li", Li, 0, z, 0); - const Double_t integ = IntLi.Integral(0, z); - - /* - if (fabs(z)<1) - { - Double_t disum = DiSum(&z); - cout << "Disum (" << z << ") " << disum << "=" << -integ << "\t" << disum-integ << endl; - return disum; - } - */ - - /* - Integral(0, 1) = konst; - Double_t konst = 1./6*TMath::Pi()*TMath::Pi(); - */ - + const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2); return -integ; } @@ -114,5 +106,5 @@ const Double_t d = w4*(w4 + 1); - Double_t s = -w*2*(1+1); // -2*omega*(1+beta) + Double_t s = -w*2*(1+1); // -2*omega*(1+beta) const Double_t li2 = Li2(&s); @@ -125,17 +117,15 @@ { const Double_t E0 = 511e-6; //[GeV] - const Double_t E02 = E0*E0; Double_t epsilon = x[0]; - Double_t E = k[0]; - // Double_t beta = k[1]; - Double_t z = k[2]; - - Double_t omega = epsilon*E/E02; - - const Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1] + Double_t z = k[1]; + + const Double_t E = k[0]; + + Double_t omega = epsilon*E/(E0*E0); + + const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1] return Flim(&omega)*n; } - Double_t MElectron::InteractionLength(Double_t *E, Double_t *k=NULL) @@ -170,8 +160,5 @@ const Double_t z = k ? k[0] : 0; - // Double_t beta = sqrt(1-E0/E*E0/E); - const Double_t beta = 1; //0.999999999999999999999999999; - - Double_t val[3] = { E[0], beta, z }; // E[GeV] + Double_t val[3] = { E[0], z }; // E[GeV] const Double_t from = 1e-17; @@ -184,9 +171,11 @@ ----------------------------------- */ - TF1 func("Compton", Compton, from, to, 3); // [0, inf] - - const Double_t integ = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf] + TF1 func("Compton", Compton, from, to, 2); // [0, inf] + + Double_t integ = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf] const Double_t aE = alpha/E[0]; // [1/GeV] + + const Double_t beta = 1; const Double_t konst = 2.*E02/hc/beta; // [1 / GeV m] @@ -196,5 +185,5 @@ const Double_t pc = 1./3.258; // [pc/ly] - return (1./ret)/ly*pc/1000; // [kpc] + return (1./ret)/ly*pc/1000; // [kpc] } @@ -211,21 +200,24 @@ // -------------------------------------------------------------------------- -Double_t MElectron::p_e(Double_t *x, Double_t *k) -{ - Double_t e = pow(10, x[0]); - Double_t z = k[1]; - - const Double_t E = k[0]; - - const Double_t E0 = 511e-6; //[GeV] - const Double_t E02 = E0*E0; - - Double_t omega = e*E/E02; - - const Double_t n = Planck(&e, &z); - - const Double_t F = Flim(&omega)/omega/omega; - - return n*F*1e26; +inline Double_t MElectron::p_e(Double_t *x, Double_t *k) +{ + Double_t e = pow(10, x[0]); + return Compton(&e, k); + /* + Double_t z = k[1]; + + const Double_t E = k[0]; + + const Double_t E0 = 511e-6; //[GeV] + const Double_t E02 = E0*E0; + + Double_t omega = e*E/E02; + + const Double_t n = Planck(&e, &z); + + const Double_t F = Flim(&omega)/omega/omega; + + return n*F*1e26; + */ } Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1364) @@ -21,5 +21,4 @@ // ---------------------------------------------------------------- - static Double_t Planck(Double_t *x, Double_t *k=NULL); static Double_t Li(Double_t *x, Double_t *k); static Double_t Li2(Double_t *x, Double_t *k=NULL); Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1364) @@ -71,13 +71,41 @@ 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 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] + return R * pc/ly/1e3; // [kpc] +} + +Double_t MParticle::Planck(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] + } Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1364) @@ -46,4 +46,10 @@ } + // ---------------------------------------------------------------- + + static Double_t Planck(Double_t *x, Double_t *k=NULL); + + // ---------------------------------------------------------------- + void SetIsPrimary(Bool_t is=kTRUE) { is ? SetBit(kEIsPrimary) : ResetBit(kEIsPrimary); } Bool_t IsPrimary() const { return TestBit(kEIsPrimary); } Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1364) @@ -38,32 +38,4 @@ ClassImp(MPhoton); - -Double_t MPhoton::Planck(Double_t *x, Double_t *k=NULL) -{ - // - // Planck, per unit volume, per unit energy - // - // constants moved out of function, see below - // - 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] -} Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL) @@ -125,6 +97,6 @@ TF1 f("int1", Int1, from, to, 2); - const Double_t int1 = f.Integral(from, to, val); // [m^2] - const Double_t planck = Planck(&Ep, &z); // [GeV^2] + const Double_t int1 = f.Integral(from, to, val, 1e-2); // [m^2] + const Double_t planck = MParticle::Planck(&Ep, &z); // [GeV^2] const Double_t res = planck * int1; @@ -157,5 +129,5 @@ TF1 f("int2", Int2, lolim, inf, 2); - Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2] + Double_t int2 = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2] if (int2==0) @@ -166,5 +138,5 @@ /* Planck constants: konst */ - Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + const Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); int2 *= konst; // [1 / m] Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1364) @@ -21,5 +21,4 @@ // ---------------------------------------------------------------- - static Double_t Planck(Double_t *x, Double_t *k=NULL); static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL); static Double_t Int1(Double_t *x, Double_t *k=NULL); Index: /trunk/WuerzburgSoft/Thomas/mphys/Makefile =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Makefile (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/Makefile (revision 1364) @@ -31,7 +31,5 @@ MPhoton.cc \ MElectron.cc \ - MGenIRPhoton.cc \ - MPairProduction.cc \ - MGenPrimaryParticle.cc + MPairProduction.cc SRCS = $(SRCFILES) Index: /trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 1364) @@ -8,6 +8,4 @@ #pragma link C++ class MElectron+; #pragma link C++ class MPhoton+; -#pragma link C++ class MGenPrimaryParticle+; -#pragma link C++ class MGenIRPhoton+; #pragma link C++ class MPairProduction+; Index: /trunk/WuerzburgSoft/Thomas/mphys/anale.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/anale.C (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/anale.C (revision 1364) @@ -4,7 +4,10 @@ { TChain chain("Electrons"); + //chain.Add("cascade.root"); + chain.Add("cascade_100kpc_01.root"); chain.Add("cascade_100kpc_02.root"); chain.Add("cascade_100kpc_03.root"); + MElectron *e = new MElectron; Index: /trunk/WuerzburgSoft/Thomas/mphys/analp.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/analp.C (revision 1364) @@ -4,7 +4,8 @@ { TChain chain("Photons"); - chain.Add("cascade_100kpc_01.root"); - chain.Add("cascade_100kpc_02.root"); - chain.Add("cascade_100kpc_03.root"); + chain.Add("cascade_500kpc_0*.root"); + chain.Add("cascade_500kpc_2*.root"); +// chain.Add("cascade_100kpc_0*.root"); +// chain.Add("cascade_100kpc_1*.root"); MPhoton *p = new MPhoton; @@ -15,6 +16,9 @@ cout << "Found " << n << " entries." << endl; + if (n==0) + return; + MBinning binsx; - binsx.SetEdgesLog(21, 1e4, 1e11); + binsx.SetEdgesLog(14, 1e4, 1e11); MBinning binspolx; @@ -22,6 +26,11 @@ MBinning binspoly2; binspolx.SetEdges(16, -180, 180); - binspoly1.SetEdges(20, 0, 5e-9); - binspoly2.SetEdges(20, 0, 2e-7); + + const Double_t ls = 299792458; // [m/s] + const Double_t ly = 3600.*24.*365.*ls; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + binspoly1.SetEdges(20, 0, 2e-6*3600); + binspoly2.SetEdges(20, 0, 1e-5*1e3/pc*365); TH1D h; @@ -59,6 +68,6 @@ h.Fill(Ep, Ep*Ep * weight); - r.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight); - a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight); + r.Fill(p->GetPhi()*kRad2Deg, p->GetR()*1e3/pc*365, weight); + a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg*3600, weight); } @@ -67,15 +76,63 @@ gStyle->SetOptStat(10); - TLine line; - line.SetLineColor(kBlack); - line.SetLineWidth(1); - - TCanvas *c = new TCanvas("c1", "name"); +// delete gROOT->FindObject("Analysis Arrival"); + c = MH::MakeDefCanvas("Analysis Arrival", ""); c->Divide(2,2); c->cd(1); - gPad->SetLogx(); gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(70); + gPad->SetPhi(90); + r.SetTitle(" Distance from Observer "); + r.GetXaxis()->SetLabelOffset(-0.015); + r.GetXaxis()->SetTitleOffset(1.1); + r.GetXaxis()->SetRangeUser(1e4, 1e9); + r.SetXTitle("\\Phi [\\circ]"); + r.SetYTitle("R [light days]"); + r.DrawCopy("surf1pol"); + c->cd(2); + gPad->SetLogy(); + gPad->SetLogz(); + gPad->SetTheta(70); + gPad->SetPhi(90); + a.SetTitle(" Hit Angle for Observer "); + a.GetXaxis()->SetLabelOffset(-0.015); + a.GetXaxis()->SetTitleOffset(1.1); + a.GetXaxis()->SetRangeUser(1e4, 1e9); + a.SetXTitle("\\Psi [\\circ]"); + a.SetYTitle("\\Theta [\"]"); + a.DrawCopy("surf1pol"); + + c->cd(3); + gPad->SetLogy(); + TH1 *g=r.ProjectionY("p1"); + g->SetBit(kCanDelete); + g->SetTitle(" Hit Observers Plain "); + g->GetXaxis()->SetLabelOffset(0); + g->GetXaxis()->SetTitleOffset(1.1); + g->GetYaxis()->SetTitleOffset(1.3); + g->SetXTitle("R [light days]"); + g->SetYTitle("Counts"); + g->Draw(); + + c->cd(4); + gPad->SetLogy(); + g=a.ProjectionY("p2"); + g->SetBit(kCanDelete); + g->SetTitle("Hit Angle"); + g->GetXaxis()->SetLabelOffset(0); + g->GetXaxis()->SetTitleOffset(1.1); + g->GetYaxis()->SetTitleOffset(1.3); + g->SetXTitle("\\Phi [\"]"); + g->SetYTitle("Counts"); + g->Draw(); + + // delete gROOT->FindObject("Analysis Photons"); + TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870); + c->Divide(1,2); + + c->cd(1); gPad->SetLogx(); gPad->SetLogy(); @@ -84,5 +141,5 @@ h.GetXaxis()->SetLabelOffset(-0.015); h.GetXaxis()->SetTitleOffset(1.1); - h.GetXaxis()->SetRangeUser(1e4, 1e9); +// h.GetXaxis()->SetRangeUser(1e4, 1e9); prim.SetMarkerStyle(kPlus); h.SetMarkerStyle(kMultiply); @@ -96,33 +153,4 @@ c->cd(2); - r.SetTitle(" Distance from Observer "); - r.GetXaxis()->SetLabelOffset(-0.015); - r.GetXaxis()->SetTitleOffset(1.1); - r.GetXaxis()->SetRangeUser(1e4, 1e9); - r.SetXTitle("\\Phi [\\circ]"); - r.SetYTitle("R [kpc]"); - r.DrawCopy("surf1pol"); - - c->cd(3); - a.SetTitle(" Hit Angle for Observer "); - a.GetXaxis()->SetLabelOffset(-0.015); - a.GetXaxis()->SetTitleOffset(1.1); - a.GetXaxis()->SetRangeUser(1e4, 1e9); - a.SetXTitle("\\Psi [\\circ]"); - a.SetYTitle("\\Theta [\\circ]"); - a.DrawCopy("surf1pol"); - - c->cd(4); - /* - TH1 *g=a.ProjectionY(); - g->SetBit(kCanDelete); - g->SetTitle("Hit Angle"); - g->GetXaxis()->SetLabelOffset(0); - g->GetXaxis()->SetTitleOffset(1.1); - g->GetYaxis()->SetTitleOffset(1.3); - g->SetXTitle("\\Phi [\\circ]"); - g->SetYTitle("Counts"); - g->Draw(); - */ gPad->SetLogx(); TH1D div; Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1363) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1364) @@ -35,6 +35,5 @@ Double_t res = MPhoton::Int2(&Ep, k); - return res*1e55; //65/k[0]; - // return MPhoton::Planck(&Ep, &k[1]); + return res*1e55; } @@ -125,5 +124,5 @@ MPairProduction pair; - Double_t runtime = 15*60; // [s] + Double_t runtime = 7*60*60; // [s] Double_t lo = 1e4; @@ -156,6 +155,6 @@ MBinning binspoly2; binspolx.SetEdges(16, -180, 180); - binspoly1.SetEdges(20, 0, 5e-9); - binspoly2.SetEdges(20, 0, 2e-7); + binspoly1.SetEdges(20, 0, 2e-6); + binspoly2.SetEdges(20, 0, 1e-5); MH::SetBinning(&angle, &binspolx, &binspoly1); MH::SetBinning(&position, &binspolx, &binspoly2); @@ -175,4 +174,6 @@ histsrc.SetFillStyle(0); histsrc.SetMarkerStyle(kMultiply); + histsrc.SetMarkerColor(kRed); + histsrc.SetLineColor(kRed); MBinning bins; @@ -272,5 +273,5 @@ phot.SetParameter(0, Eg); phot.SetParameter(1, p->GetZ()); - if (phot.Integral(-log10(Eg)-8, -10.5)==0) + if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0) { cout << "z" << flush;