- Timestamp:
- 06/14/02 09:03:11 (22 years ago)
- Location:
- trunk/WuerzburgSoft/Thomas/mphys
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/WuerzburgSoft/Thomas/mphys/Changelog
r1363 r1364 1 1 -*-*- END -*-*- 2 2002/06/14: Thomas Bretz 3 4 * MElectron.[h,cc]: 5 - Moved Planck function to MParticle 6 - Added the DiSum function 7 - changed Li2 to use the DiSum function (speed reasons) 8 - changed the eps of all integrals to 1e-2 9 - changed the p_e function to use the Compton function 10 11 * MPhoton.[h,cc]: 12 - removed the planck function 13 - changed the eps of all integrals to 1e-2 14 15 * MParticle.[h,cc]: 16 - added the planck function 17 18 * Makefile: 19 - removed unused source files 20 21 * analp.C: 22 - did some small fixes 23 - changed the sclae of the radius- and phi-plots 24 25 * phys.C: 26 - changed the integral eps to 1e-2 27 28 29 30 2 31 2002/06/13: Thomas Bretz 3 32 -
trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc
r1363 r1364 42 42 ClassImp(MElectron); 43 43 44 Double_t MElectron::Planck(Double_t *x, Double_t *k=NULL)45 {46 //47 // Planck, per unit volume, per unit energy48 //49 // constants moved out of function50 //51 const Double_t E = x[0]; // [GeV]52 const Double_t z = k ? k[0] : 0;53 54 const Double_t T = 2.96*(z+1); // [K]55 const Double_t e = 1.602176462e-19; // [C]56 const Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]57 58 const Double_t EkT = E/kB/T;59 60 /*61 //Double_t c = 299792458; // [m/s]62 //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]63 //Double_t hc = h*c; // [GeVm]64 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);65 return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]66 */67 68 return E*E / (exp(EkT)-1.); // [GeV^2]69 }70 71 44 Double_t MElectron::Li(Double_t *x, Double_t *k) 72 45 { 73 46 const Double_t t = x[0]; 74 75 47 return log(1.-t)/t; 76 48 } 77 49 50 Double_t DiSum(Double_t *x, Double_t *k=NULL) 51 { 52 Double_t t = x[0]; 53 54 const Double_t eps = fabs(t*1e-2); 55 56 Double_t disum = t; 57 Double_t add = 0; 58 59 Int_t n = 2; 60 Double_t pow = t*t; // t^2 61 62 do 63 { 64 add = pow/n/n; 65 66 pow *= t; // pow = t^n 67 n++; 68 69 disum += add; 70 71 } while (fabs(add)>eps); 72 73 return disum; 74 } 75 78 76 Double_t MElectron::Li2(Double_t *x, Double_t *k=NULL) 79 77 { 80 78 // 81 79 // Dilog, Li2 82 // 83 Double_t z = x[0]; 80 // ---------- 81 // 82 // Integral(0, 1) = konst; 83 // Double_t konst = 1./6*TMath::Pi()*TMath::Pi(); 84 // 85 // x[0]: z 86 // 87 const Double_t z = x[0]; 88 89 if (fabs(z)<1) 90 return DiSum(x); 84 91 85 92 TF1 IntLi("Li", Li, 0, z, 0); 86 const Double_t integ = IntLi.Integral(0, z); 87 88 /* 89 if (fabs(z)<1) 90 { 91 Double_t disum = DiSum(&z); 92 cout << "Disum (" << z << ") " << disum << "=" << -integ << "\t" << disum-integ << endl; 93 return disum; 94 } 95 */ 96 97 /* 98 Integral(0, 1) = konst; 99 Double_t konst = 1./6*TMath::Pi()*TMath::Pi(); 100 */ 101 93 const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2); 102 94 return -integ; 103 95 } … … 114 106 const Double_t d = w4*(w4 + 1); 115 107 116 Double_t s 108 Double_t s = -w*2*(1+1); // -2*omega*(1+beta) 117 109 const Double_t li2 = Li2(&s); 118 110 … … 125 117 { 126 118 const Double_t E0 = 511e-6; //[GeV] 127 const Double_t E02 = E0*E0;128 119 129 120 Double_t epsilon = x[0]; 130 Double_t E = k[0];131 // Double_t beta = k[1]; 132 Double_t z = k[2];133 134 Double_t omega = epsilon*E/ E02;135 136 const Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]121 Double_t z = k[1]; 122 123 const Double_t E = k[0]; 124 125 Double_t omega = epsilon*E/(E0*E0); 126 127 const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1] 137 128 return Flim(&omega)*n; 138 129 } 139 140 130 141 131 Double_t MElectron::InteractionLength(Double_t *E, Double_t *k=NULL) … … 170 160 const Double_t z = k ? k[0] : 0; 171 161 172 // Double_t beta = sqrt(1-E0/E*E0/E); 173 const Double_t beta = 1; //0.999999999999999999999999999; 174 175 Double_t val[3] = { E[0], beta, z }; // E[GeV] 162 Double_t val[3] = { E[0], z }; // E[GeV] 176 163 177 164 const Double_t from = 1e-17; … … 184 171 ----------------------------------- 185 172 */ 186 TF1 func("Compton", Compton, from, to, 3);// [0, inf]187 188 const Double_t integ = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]173 TF1 func("Compton", Compton, from, to, 2); // [0, inf] 174 175 Double_t integ = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf] 189 176 190 177 const Double_t aE = alpha/E[0]; // [1/GeV] 178 179 const Double_t beta = 1; 191 180 192 181 const Double_t konst = 2.*E02/hc/beta; // [1 / GeV m] … … 196 185 const Double_t pc = 1./3.258; // [pc/ly] 197 186 198 return (1./ret)/ly*pc/1000; // [kpc]187 return (1./ret)/ly*pc/1000; // [kpc] 199 188 } 200 189 … … 211 200 // -------------------------------------------------------------------------- 212 201 213 Double_t MElectron::p_e(Double_t *x, Double_t *k) 214 { 215 Double_t e = pow(10, x[0]); 216 Double_t z = k[1]; 217 218 const Double_t E = k[0]; 219 220 const Double_t E0 = 511e-6; //[GeV] 221 const Double_t E02 = E0*E0; 222 223 Double_t omega = e*E/E02; 224 225 const Double_t n = Planck(&e, &z); 226 227 const Double_t F = Flim(&omega)/omega/omega; 228 229 return n*F*1e26; 202 inline Double_t MElectron::p_e(Double_t *x, Double_t *k) 203 { 204 Double_t e = pow(10, x[0]); 205 return Compton(&e, k); 206 /* 207 Double_t z = k[1]; 208 209 const Double_t E = k[0]; 210 211 const Double_t E0 = 511e-6; //[GeV] 212 const Double_t E02 = E0*E0; 213 214 Double_t omega = e*E/E02; 215 216 const Double_t n = Planck(&e, &z); 217 218 const Double_t F = Flim(&omega)/omega/omega; 219 220 return n*F*1e26; 221 */ 230 222 } 231 223 -
trunk/WuerzburgSoft/Thomas/mphys/MElectron.h
r1356 r1364 21 21 // ---------------------------------------------------------------- 22 22 23 static Double_t Planck(Double_t *x, Double_t *k=NULL);24 23 static Double_t Li(Double_t *x, Double_t *k); 25 24 static Double_t Li2(Double_t *x, Double_t *k=NULL); -
trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc
r1359 r1364 71 71 Double_t z1 = x[0] + 1; 72 72 73 const Double_t c = 299792458; // [m/s]74 const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]75 76 const Double_t ly = 3600.*24.*365.*c; // [m/ly]77 const Double_t pc = 1./3.258; // [pc/ly]73 const Double_t c = 299792458; // [m/s] 74 const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] 75 76 const Double_t ly = 3600.*24.*365.*c; // [m/ly] 77 const Double_t pc = 1./3.258; // [pc/ly] 78 78 79 79 const Double_t R = c/H0 * 2 * (1 - 1./sqrt(z1)); // [m] 80 80 81 return R * pc/ly/1e3; // [kpc] 81 return R * pc/ly/1e3; // [kpc] 82 } 83 84 Double_t MParticle::Planck(Double_t *x, Double_t *k) 85 { 86 // 87 // Planck, per unit volume, per unit energy 88 // 89 // constants (see below) moved out of function 90 // 91 const Double_t E = x[0]; // [GeV] 92 const Double_t z = k ? k[0] : 0; 93 94 const Double_t T = 2.96*(z+1); // [K] 95 const Double_t e = 1.602176462e-19; // [C] 96 const Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] 97 98 const Double_t EkT = E/kB/T; 99 100 /* 101 //Double_t c = 299792458; // [m/s] 102 //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] 103 //Double_t hc = h*c; // [GeVm] 104 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); 105 return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] 106 */ 107 108 return E*E / (exp(EkT)-1.); // [GeV^2] 109 82 110 } 83 111 -
trunk/WuerzburgSoft/Thomas/mphys/MParticle.h
r1362 r1364 46 46 } 47 47 48 // ---------------------------------------------------------------- 49 50 static Double_t Planck(Double_t *x, Double_t *k=NULL); 51 52 // ---------------------------------------------------------------- 53 48 54 void SetIsPrimary(Bool_t is=kTRUE) { is ? SetBit(kEIsPrimary) : ResetBit(kEIsPrimary); } 49 55 Bool_t IsPrimary() const { return TestBit(kEIsPrimary); } -
trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc
r1363 r1364 38 38 39 39 ClassImp(MPhoton); 40 41 Double_t MPhoton::Planck(Double_t *x, Double_t *k=NULL)42 {43 //44 // Planck, per unit volume, per unit energy45 //46 // constants moved out of function, see below47 //48 const Double_t E = x[0]; // [GeV]49 const Double_t z = k ? k[0] : 0;50 51 const Double_t T = 2.96*(z+1); // [K]52 const Double_t e = 1.602176462e-19; // [C]53 const Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]54 55 const Double_t EkT = E/kB/T;56 57 /*58 //Double_t c = 299792458; // [m/s]59 //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]60 //Double_t hc = h*c; // [GeVm]61 62 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);63 return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]64 */65 66 return E*E / (exp(EkT)-1.); // [GeV^2]67 }68 40 69 41 Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL) … … 125 97 TF1 f("int1", Int1, from, to, 2); 126 98 127 const Double_t int1 = f.Integral(from, to, val ); // [m^2]128 const Double_t planck = Planck(&Ep, &z);// [GeV^2]99 const Double_t int1 = f.Integral(from, to, val, 1e-2); // [m^2] 100 const Double_t planck = MParticle::Planck(&Ep, &z); // [GeV^2] 129 101 130 102 const Double_t res = planck * int1; … … 157 129 TF1 f("int2", Int2, lolim, inf, 2); 158 130 159 Double_t int2 = f.Integral(lolim, inf, val ); //[GeV^3 m^2]131 Double_t int2 = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2] 160 132 161 133 if (int2==0) … … 166 138 167 139 /* Planck constants: konst */ 168 Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);140 const Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); 169 141 int2 *= konst; // [1 / m] 170 142 -
trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h
r1363 r1364 21 21 // ---------------------------------------------------------------- 22 22 23 static Double_t Planck(Double_t *x, Double_t *k=NULL);24 23 static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL); 25 24 static Double_t Int1(Double_t *x, Double_t *k=NULL); -
trunk/WuerzburgSoft/Thomas/mphys/Makefile
r1349 r1364 31 31 MPhoton.cc \ 32 32 MElectron.cc \ 33 MGenIRPhoton.cc \ 34 MPairProduction.cc \ 35 MGenPrimaryParticle.cc 33 MPairProduction.cc 36 34 37 35 SRCS = $(SRCFILES) -
trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h
r1349 r1364 8 8 #pragma link C++ class MElectron+; 9 9 #pragma link C++ class MPhoton+; 10 #pragma link C++ class MGenPrimaryParticle+;11 #pragma link C++ class MGenIRPhoton+;12 10 #pragma link C++ class MPairProduction+; 13 11 -
trunk/WuerzburgSoft/Thomas/mphys/anale.C
r1363 r1364 4 4 { 5 5 TChain chain("Electrons"); 6 //chain.Add("cascade.root"); 7 6 8 chain.Add("cascade_100kpc_01.root"); 7 9 chain.Add("cascade_100kpc_02.root"); 8 10 chain.Add("cascade_100kpc_03.root"); 11 9 12 10 13 MElectron *e = new MElectron; -
trunk/WuerzburgSoft/Thomas/mphys/analp.C
r1363 r1364 4 4 { 5 5 TChain chain("Photons"); 6 chain.Add("cascade_100kpc_01.root"); 7 chain.Add("cascade_100kpc_02.root"); 8 chain.Add("cascade_100kpc_03.root"); 6 chain.Add("cascade_500kpc_0*.root"); 7 chain.Add("cascade_500kpc_2*.root"); 8 // chain.Add("cascade_100kpc_0*.root"); 9 // chain.Add("cascade_100kpc_1*.root"); 9 10 10 11 MPhoton *p = new MPhoton; … … 15 16 cout << "Found " << n << " entries." << endl; 16 17 18 if (n==0) 19 return; 20 17 21 MBinning binsx; 18 binsx.SetEdgesLog( 21, 1e4, 1e11);22 binsx.SetEdgesLog(14, 1e4, 1e11); 19 23 20 24 MBinning binspolx; … … 22 26 MBinning binspoly2; 23 27 binspolx.SetEdges(16, -180, 180); 24 binspoly1.SetEdges(20, 0, 5e-9); 25 binspoly2.SetEdges(20, 0, 2e-7); 28 29 const Double_t ls = 299792458; // [m/s] 30 const Double_t ly = 3600.*24.*365.*ls; // [m/ly] 31 const Double_t pc = 1./3.258; // [pc/ly] 32 33 binspoly1.SetEdges(20, 0, 2e-6*3600); 34 binspoly2.SetEdges(20, 0, 1e-5*1e3/pc*365); 26 35 27 36 TH1D h; … … 59 68 60 69 h.Fill(Ep, Ep*Ep * weight); 61 r.Fill(p->GetPhi()*kRad2Deg, p->GetR() , weight);62 a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg , weight);70 r.Fill(p->GetPhi()*kRad2Deg, p->GetR()*1e3/pc*365, weight); 71 a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg*3600, weight); 63 72 } 64 73 … … 67 76 gStyle->SetOptStat(10); 68 77 69 TLine line; 70 line.SetLineColor(kBlack); 71 line.SetLineWidth(1); 72 73 TCanvas *c = new TCanvas("c1", "name"); 78 // delete gROOT->FindObject("Analysis Arrival"); 79 c = MH::MakeDefCanvas("Analysis Arrival", ""); 74 80 c->Divide(2,2); 75 81 76 82 c->cd(1); 77 gPad->SetLogx();78 83 gPad->SetLogy(); 84 gPad->SetLogz(); 85 gPad->SetTheta(70); 86 gPad->SetPhi(90); 87 r.SetTitle(" Distance from Observer "); 88 r.GetXaxis()->SetLabelOffset(-0.015); 89 r.GetXaxis()->SetTitleOffset(1.1); 90 r.GetXaxis()->SetRangeUser(1e4, 1e9); 91 r.SetXTitle("\\Phi [\\circ]"); 92 r.SetYTitle("R [light days]"); 93 r.DrawCopy("surf1pol"); 79 94 95 c->cd(2); 96 gPad->SetLogy(); 97 gPad->SetLogz(); 98 gPad->SetTheta(70); 99 gPad->SetPhi(90); 100 a.SetTitle(" Hit Angle for Observer "); 101 a.GetXaxis()->SetLabelOffset(-0.015); 102 a.GetXaxis()->SetTitleOffset(1.1); 103 a.GetXaxis()->SetRangeUser(1e4, 1e9); 104 a.SetXTitle("\\Psi [\\circ]"); 105 a.SetYTitle("\\Theta [\"]"); 106 a.DrawCopy("surf1pol"); 107 108 c->cd(3); 109 gPad->SetLogy(); 110 TH1 *g=r.ProjectionY("p1"); 111 g->SetBit(kCanDelete); 112 g->SetTitle(" Hit Observers Plain "); 113 g->GetXaxis()->SetLabelOffset(0); 114 g->GetXaxis()->SetTitleOffset(1.1); 115 g->GetYaxis()->SetTitleOffset(1.3); 116 g->SetXTitle("R [light days]"); 117 g->SetYTitle("Counts"); 118 g->Draw(); 119 120 c->cd(4); 121 gPad->SetLogy(); 122 g=a.ProjectionY("p2"); 123 g->SetBit(kCanDelete); 124 g->SetTitle("Hit Angle"); 125 g->GetXaxis()->SetLabelOffset(0); 126 g->GetXaxis()->SetTitleOffset(1.1); 127 g->GetYaxis()->SetTitleOffset(1.3); 128 g->SetXTitle("\\Phi [\"]"); 129 g->SetYTitle("Counts"); 130 g->Draw(); 131 132 // delete gROOT->FindObject("Analysis Photons"); 133 TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870); 134 c->Divide(1,2); 135 136 c->cd(1); 80 137 gPad->SetLogx(); 81 138 gPad->SetLogy(); … … 84 141 h.GetXaxis()->SetLabelOffset(-0.015); 85 142 h.GetXaxis()->SetTitleOffset(1.1); 86 h.GetXaxis()->SetRangeUser(1e4, 1e9);143 // h.GetXaxis()->SetRangeUser(1e4, 1e9); 87 144 prim.SetMarkerStyle(kPlus); 88 145 h.SetMarkerStyle(kMultiply); … … 96 153 97 154 c->cd(2); 98 r.SetTitle(" Distance from Observer ");99 r.GetXaxis()->SetLabelOffset(-0.015);100 r.GetXaxis()->SetTitleOffset(1.1);101 r.GetXaxis()->SetRangeUser(1e4, 1e9);102 r.SetXTitle("\\Phi [\\circ]");103 r.SetYTitle("R [kpc]");104 r.DrawCopy("surf1pol");105 106 c->cd(3);107 a.SetTitle(" Hit Angle for Observer ");108 a.GetXaxis()->SetLabelOffset(-0.015);109 a.GetXaxis()->SetTitleOffset(1.1);110 a.GetXaxis()->SetRangeUser(1e4, 1e9);111 a.SetXTitle("\\Psi [\\circ]");112 a.SetYTitle("\\Theta [\\circ]");113 a.DrawCopy("surf1pol");114 115 c->cd(4);116 /*117 TH1 *g=a.ProjectionY();118 g->SetBit(kCanDelete);119 g->SetTitle("Hit Angle");120 g->GetXaxis()->SetLabelOffset(0);121 g->GetXaxis()->SetTitleOffset(1.1);122 g->GetYaxis()->SetTitleOffset(1.3);123 g->SetXTitle("\\Phi [\\circ]");124 g->SetYTitle("Counts");125 g->Draw();126 */127 155 gPad->SetLogx(); 128 156 TH1D div; -
trunk/WuerzburgSoft/Thomas/mphys/phys.C
r1363 r1364 35 35 Double_t res = MPhoton::Int2(&Ep, k); 36 36 37 return res*1e55; //65/k[0]; 38 // return MPhoton::Planck(&Ep, &k[1]); 37 return res*1e55; 39 38 } 40 39 … … 125 124 MPairProduction pair; 126 125 127 Double_t runtime = 15*60; // [s]126 Double_t runtime = 7*60*60; // [s] 128 127 129 128 Double_t lo = 1e4; … … 156 155 MBinning binspoly2; 157 156 binspolx.SetEdges(16, -180, 180); 158 binspoly1.SetEdges(20, 0, 5e-9);159 binspoly2.SetEdges(20, 0, 2e-7);157 binspoly1.SetEdges(20, 0, 2e-6); 158 binspoly2.SetEdges(20, 0, 1e-5); 160 159 MH::SetBinning(&angle, &binspolx, &binspoly1); 161 160 MH::SetBinning(&position, &binspolx, &binspoly2); … … 175 174 histsrc.SetFillStyle(0); 176 175 histsrc.SetMarkerStyle(kMultiply); 176 histsrc.SetMarkerColor(kRed); 177 histsrc.SetLineColor(kRed); 177 178 178 179 MBinning bins; … … 272 273 phot.SetParameter(0, Eg); 273 274 phot.SetParameter(1, p->GetZ()); 274 if (phot.Integral(-log10(Eg)-8, -10.5 )==0)275 if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0) 275 276 { 276 277 cout << "z" << flush;
Note:
See TracChangeset
for help on using the changeset viewer.