Changeset 1363 for trunk/WuerzburgSoft/Thomas/mphys
- Timestamp:
- 06/13/02 16:36:19 (23 years ago)
- Location:
- trunk/WuerzburgSoft/Thomas/mphys
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/WuerzburgSoft/Thomas/mphys/Changelog
r1360 r1363 4 4 * MParticle.h: 5 5 - Implemented Set/IsPrimary 6 - removed pure abstract method (root can't read this) 7 - added default constructor 8 9 * MElectron.cc: 10 - changed many variables to const 11 - changed the output 12 13 * MPairProduction.cc: 14 - changed many variables to const 15 16 * MPhoton.[h,cc]: 17 - changed the upper and lower limit of the integration in Int2 18 - added upper and lower limit of draw to DrawInteractionLength 6 19 7 20 * energyloss.C: … … 11 24 - changed energy randomizer 12 25 - added root file output 26 - removed log scale from point spread plots 27 28 * analp.C: 29 - added macro to analize the photons from the root file 30 31 * anale.C: 32 - added macro to analize the electrons from the root file 13 33 14 34 -
trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc
r1358 r1363 49 49 // constants moved out of function 50 50 // 51 Double_t E = x[0]; // [GeV]52 Double_t z = k ? k[0] : 0;53 54 Double_t T = 2.96*(z+1); // [K]55 Double_t e = 1.602176462e-19; // [C]56 Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K]57 58 Double_t EkT = E/kB/T;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 59 60 60 /* … … 71 71 Double_t MElectron::Li(Double_t *x, Double_t *k) 72 72 { 73 Double_t t = x[0];73 const Double_t t = x[0]; 74 74 75 75 return log(1.-t)/t; … … 84 84 85 85 TF1 IntLi("Li", Li, 0, z, 0); 86 Double_t integ = IntLi.Integral(0, z);86 const Double_t integ = IntLi.Integral(0, z); 87 87 88 88 /* … … 105 105 Double_t MElectron::Flim(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam) mit b-->1 (Maple) 106 106 { 107 Double_t w = x[0];108 109 Double_t w4 = w*4;110 Double_t wsqr = w*w;111 112 Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);113 Double_t u2 = -w4*(wsqr*2 + w*9 + 2);114 Double_t d = w4*(w4 + 1);107 const Double_t w = x[0]; 108 109 const Double_t w4 = w*4; 110 const Double_t wsqr = w*w; 111 112 const Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1); 113 const Double_t u2 = -w4*(wsqr*2 + w*9 + 2); 114 const Double_t d = w4*(w4 + 1); 115 115 116 116 Double_t s = -w*2*(1+1); // -2*omega*(1+beta) 117 Double_t li2 = Li2(&s);118 119 Double_t res = (u1+u2)/d + li2;117 const Double_t li2 = Li2(&s); 118 119 const Double_t res = (u1+u2)/d + li2; 120 120 121 121 return res; //<1e-10? 0 : res; … … 124 124 Double_t MElectron::Compton(Double_t *x, Double_t *k) 125 125 { 126 Double_t E0 = 511e-6; //[GeV]127 Double_t E02 = E0*E0;126 const Double_t E0 = 511e-6; //[GeV] 127 const Double_t E02 = E0*E0; 128 128 129 129 Double_t epsilon = x[0]; … … 134 134 Double_t omega = epsilon*E/E02; 135 135 136 Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]136 const Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1] 137 137 return Flim(&omega)*n; 138 138 } … … 160 160 // F(o-)~F(0) = 14/8 = 1.75 161 161 162 Double_t E0 = 511e-6; // [GeV]163 Double_t E02 = E0*E0; // [GeV^2]164 Double_t c = 299792458; // [m/s]165 Double_t e = 1.602176462e-19; // [C]166 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]167 Double_t hc = h*c; // [GeVm]168 Double_t alpha = 1./137.; // [1]169 170 Double_t z = k ? k[0] : 0;162 const Double_t E0 = 511e-6; // [GeV] 163 const Double_t E02 = E0*E0; // [GeV^2] 164 const Double_t c = 299792458; // [m/s] 165 const Double_t e = 1.602176462e-19; // [C] 166 const Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] 167 const Double_t hc = h*c; // [GeVm] 168 const Double_t alpha = 1./137.; // [1] 169 170 const Double_t z = k ? k[0] : 0; 171 171 172 172 // Double_t beta = sqrt(1-E0/E*E0/E); 173 Double_t beta = 1; //0.999999999999999999999999999;173 const Double_t beta = 1; //0.999999999999999999999999999; 174 174 175 175 Double_t val[3] = { E[0], beta, z }; // E[GeV] 176 176 177 Double_t from = 1e-17;178 Double_t to = 1e-11;177 const Double_t from = 1e-17; 178 const Double_t to = 1e-11; 179 179 180 180 /* -------------- old ---------------- … … 186 186 TF1 func("Compton", Compton, from, to, 3); // [0, inf] 187 187 188 Double_t integ = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]189 190 Double_t aE = alpha/E[0]; // [1/GeV]191 192 Double_t konst = 2.*E02/hc/beta; // [1 / GeV m]193 Double_t ret = konst * (aE*aE) * integ; // [1 / m]194 195 Double_t ly = 3600.*24.*365.*c; // [m/ly]196 Double_t pc = 1./3.258; // [pc/ly]188 const Double_t integ = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf] 189 190 const Double_t aE = alpha/E[0]; // [1/GeV] 191 192 const Double_t konst = 2.*E02/hc/beta; // [1 / GeV m] 193 const Double_t ret = konst * (aE*aE) * integ; // [1 / m] 194 195 const Double_t ly = 3600.*24.*365.*c; // [m/ly] 196 const Double_t pc = 1./3.258; // [pc/ly] 197 197 198 198 return (1./ret)/ly*pc/1000; // [kpc] … … 214 214 { 215 215 Double_t e = pow(10, x[0]); 216 Double_t E = k[0];217 216 Double_t z = k[1]; 218 217 219 Double_t E0 = 511e-6; //[GeV] 220 Double_t E02 = E0*E0; 218 const Double_t E = k[0]; 219 220 const Double_t E0 = 511e-6; //[GeV] 221 const Double_t E02 = E0*E0; 221 222 222 223 Double_t omega = e*E/E02; 223 224 224 Double_t n = Planck(&e, &z);225 226 Double_t F = Flim(&omega)/omega/omega;225 const Double_t n = Planck(&e, &z); 226 227 const Double_t F = Flim(&omega)/omega/omega; 227 228 228 229 return n*F*1e26; … … 231 232 Double_t MElectron::G_q(Double_t *x, Double_t *k) 232 233 { 233 Double_t q = x[0];234 Double_t Gamma = k[0];235 236 Double_t Gq = Gamma*q;237 238 Double_t s1 = 2.*q*log(q);239 Double_t s2 = (1.+2.*q);240 Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;234 const Double_t q = x[0]; 235 const Double_t Gamma = k[0]; 236 237 const Double_t Gq = Gamma*q; 238 239 const Double_t s1 = 2.*q*log(q); 240 const Double_t s2 = (1.+2.*q); 241 const Double_t s3 = (Gq*Gq)/(1.+Gq)/2.; 241 242 242 243 return s1+(s2+s3)*(1.-q); … … 246 247 Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep) 247 248 { 248 Double_t E = x[0];249 Double_t z = k ? k[0] : 0;250 251 Double_t E0 = 511e-6; //[GeV]252 253 Double_t lolim = -log10(E)/7*4-13;249 const Double_t E = x[0]; 250 const Double_t z = k ? k[0] : 0; 251 252 const Double_t E0 = 511e-6; //[GeV] 253 254 const Double_t lolim = -log10(E)/7*4-13; 254 255 255 256 TF1 fP("p_e", p_e, lolim, -10.8, 2); … … 259 260 fP.SetParameter(1, z); 260 261 261 Double_t e = pow(10, fP.GetRandom());262 const Double_t e = pow(10, fP.GetRandom()); 262 263 263 264 if (ep) 264 265 *ep = e; 265 266 266 Double_t omega = e*E/E0/E0;267 Double_t Gamma = 4.*omega;267 const Double_t omega = e*E/E0/E0; 268 const Double_t Gamma = 4.*omega; 268 269 269 270 // --old-- fQ.SetRange(1e-6, 1./(1.+ 2.*Gamma)); 270 271 fQ.SetParameter(0, Gamma); 271 272 272 Double_t q = fQ.GetRandom();273 Double_t Gq = Gamma*q;274 275 Double_t e1 = Gq*E/(1.+Gq);273 const Double_t q = fQ.GetRandom(); 274 const Double_t Gq = Gamma*q; 275 276 const Double_t e1 = Gq*E/(1.+Gq); 276 277 277 278 return e1; … … 292 293 static TRandom rand(0); 293 294 294 Double_t E0 = 511e-6; //[GeV]295 const Double_t E0 = 511e-6; //[GeV] 295 296 296 297 Double_t epsilon; 297 Double_t e = GetEnergyLoss(&epsilon);298 const Double_t e = GetEnergyLoss(&epsilon); 298 299 299 300 // er: photon energy before interaction, rest frame 300 301 // e: photon energy after interaction, lab 301 302 302 Double_t gamma = fEnergy/E0; 303 Double_t beta = sqrt(1.-1./(gamma*gamma)); 304 //Double_t gammabeta = sqrt(gamma*gamma-1); 305 306 Double_t f = fEnergy/e; 307 308 Double_t t; 303 const Double_t gamma = fEnergy/E0; 304 const Double_t beta = sqrt(1.-1./(gamma*gamma)); 305 306 const Double_t f = fEnergy/e; 307 308 Double_t t=-1; 309 309 Double_t arg; 310 310 do 311 311 { 312 t = rand.Uniform(TMath::Pi()*2); 312 if (t>0) 313 cout << "~" << flush; 314 t = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4; 313 315 Double_t er = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame 314 316 arg = (f - E0/er - 1)/(f*beta+1); 315 cout << "~" << flush;316 317 317 318 } while (arg<-1 || arg>1); 318 319 319 Double_t theta1s = acos(arg); 320 Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta))); 321 322 Double_t thetastar = thetas-theta1s; 323 324 Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta))); 325 326 /* 327 cout << "(" << theta1/TMath::Pi()*180 << ","; 328 cout << theta1s/TMath::Pi()*180<< ","; 329 cout << arg << ")" << flush; 330 */ 320 const Double_t theta1s = acos(arg); 321 const Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta))); 322 323 const Double_t thetastar = thetas-theta1s; 324 325 const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta))); 331 326 332 327 fEnergy -= e; … … 335 330 p = *this; 336 331 p.SetNewDirection(theta1, rand.Uniform(TMath::Pi()*2)); 337 338 // MISSING: Electron angle339 //340 // E1 = fEnergy (after!)341 //342 // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0)343 332 344 333 return &p; -
trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc
r1358 r1363 79 79 // theta: interaction angle [rad] 80 80 // 81 82 83 81 const Double_t E0 = 511e-6; // [GeV] 84 82 const Double_t Eg = gamma->GetEnergy(); // [GeV] … … 98 96 const Double_t GammaH = (Eg+Ep)/sqrt(s); 99 97 100 Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));98 const Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah))); 101 99 102 100 fAngle->SetParameter(0, sqrt(sqrbetae)); … … 104 102 const Double_t alpha = psi-acos(fAngle->GetRandom()); 105 103 106 Double_t salpha = sin(alpha);107 Double_t calpha = cos(alpha);104 const Double_t salpha = sin(alpha); 105 const Double_t calpha = cos(alpha); 108 106 109 107 const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi) 110 108 111 Double_t bb = sqrt(sqrbetah/sqrbetae);109 const Double_t bb = sqrt(sqrbetah/sqrbetae); 112 110 113 Double_t s1 = calpha/GammaH;114 Double_t s2 = tphi*s1 - salpha - bb;111 const Double_t s1 = calpha/GammaH; 112 const Double_t s2 = tphi*s1 - salpha - bb; 115 113 116 Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;117 Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;114 const Double_t tan1 = ((salpha+bb)*tphi+s1)/s2; 115 const Double_t tan2 = ((salpha-bb)*tphi+s1)/s2; 118 116 119 117 const Double_t E = (Eg+Ep)/2;; -
trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc
r1356 r1363 152 152 Double_t val[2] = { Eg, z }; 153 153 154 Double_t lolim = E0*E0 > 1e-8 ? E0*E0/Eg : 1e-8/Eg;155 Double_t inf = 3e-11;154 Double_t lolim = E0*E0/Eg; 155 Double_t inf = Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1); 156 156 157 157 TF1 f("int2", Int2, lolim, inf, 2); … … 188 188 } 189 189 190 void MPhoton::DrawInteractionLength(Double_t z )190 void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to) 191 191 { 192 192 if (!gPad) … … 195 195 gPad->GetVirtCanvas()->cd(4); 196 196 197 TF1 f1("length", InteractionLength, 1e4, 1e11, 1);197 TF1 f1("length", InteractionLength, from, to, 1); 198 198 f1.SetParameter(0, z); 199 199 -
trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h
r1360 r1363 32 32 // ---------------------------------------------------------------- 33 33 34 static void DrawInteractionLength(Double_t z );34 static void DrawInteractionLength(Double_t z, Double_t from=1e4, Double_t to=1e11); 35 35 void DrawInteractionLength() const; 36 36 -
trunk/WuerzburgSoft/Thomas/mphys/anale.C
r1362 r1363 4 4 { 5 5 TChain chain("Electrons"); 6 chain.Add("cascade01.root"); 7 chain.Add("cascade02.root"); 6 chain.Add("cascade_100kpc_01.root"); 7 chain.Add("cascade_100kpc_02.root"); 8 chain.Add("cascade_100kpc_03.root"); 8 9 9 10 MElectron *e = new MElectron; -
trunk/WuerzburgSoft/Thomas/mphys/analp.C
r1362 r1363 4 4 { 5 5 TChain chain("Photons"); 6 chain.Add("cascade01.root"); 7 chain.Add("cascade02.root"); 6 chain.Add("cascade_100kpc_01.root"); 7 chain.Add("cascade_100kpc_02.root"); 8 chain.Add("cascade_100kpc_03.root"); 8 9 9 10 MPhoton *p = new MPhoton; -
trunk/WuerzburgSoft/Thomas/mphys/energyloss.C
r1360 r1363 477 477 // ------------------------------------------------------------------- 478 478 479 Double_t func(Double_t *x, Double_t *k) 480 { 481 return MPhoton::Int2(x, k)*1e68; 482 } 483 479 484 void energyloss() 480 485 { 486 /* Double_t E0 = 511e-6; // [GeV] 487 488 Double_t Eg = 1e4; //3.6e4; 489 Double_t z = 5; 490 491 Double_t val[2] = { Eg, z }; 492 493 Double_t lolim = E0*E0/Eg; 494 Double_t inf = Eg<1e6 ? 3e-11*z : 3e-12*z; 495 496 cout << Eg << " " << z << " " << lolim << " " << inf << endl; 497 498 TF1 f("int2", func, lolim, inf, 2); 499 500 Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2] 501 502 f.SetParameter(0, Eg); 503 f.SetParameter(1, z); 504 505 cout << int2 << endl; 506 507 new TCanvas("ILPhoton", "Mean Interaction Length Photon"); 508 509 gPad->SetLogx(); 510 // gPad->SetLogy(); 511 gPad->SetGrid(); 512 513 f.SetLineWidth(1); 514 f.DrawCopy(); 515 516 gPad->Modified(); 517 gPad->Update(); 518 519 // return; 520 */ MPhoton p; 521 p.DrawInteractionLength(); 522 523 return; 481 524 // EnergyLossRate(); 482 525 -
trunk/WuerzburgSoft/Thomas/mphys/phys.C
r1361 r1363 116 116 void DoIt() 117 117 { 118 // Double_t R = 1798; // [kpc]119 Double_t startz = 0.003; //MParticle::ZofR(&R);120 121 //cout << "R = " << R << "kpc" << endl;118 Double_t R = 500; // [kpc] 119 Double_t startz = MParticle::ZofR(&R); 120 121 cout << "R = " << R << "kpc" << endl; 122 122 cout << "Z = " << startz << endl; 123 123 … … 131 131 Double_t alpha = -2; 132 132 133 Int_t nbins = 18; //4*log10(hi/lo);133 Int_t nbins = 21; //4*log10(hi/lo); 134 134 135 135 TH1D hist; … … 153 153 154 154 MBinning binspolx; 155 MBinning binspoly; 155 MBinning binspoly1; 156 MBinning binspoly2; 156 157 binspolx.SetEdges(16, -180, 180); 157 binspoly.SetEdgesLog(20, 1e-10, 1e-3); 158 MH::SetBinning(&position, &binspolx, &binspoly); 159 MH::SetBinning(&angle, &binspolx, &binspoly); 160 MH::SetBinning(&histpos, &binspoly); 161 MH::SetBinning(&hista, &binspoly); 158 binspoly1.SetEdges(20, 0, 5e-9); 159 binspoly2.SetEdges(20, 0, 2e-7); 160 MH::SetBinning(&angle, &binspolx, &binspoly1); 161 MH::SetBinning(&position, &binspolx, &binspoly2); 162 MH::SetBinning(&hista, &binspoly1); 163 MH::SetBinning(&histpos, &binspoly2); 162 164 163 165 hista.SetTitle("Particle Angle"); … … 257 259 hista.Fill(p->GetTheta()*kRad2Deg, weight); 258 260 261 p->SetIsPrimary(kFALSE); 259 262 T->Fill(); 260 263 … … 393 396 position.SetYTitle("Pos: R [kpc]"); 394 397 position.DrawCopy("surf1pol"); 395 gPad->SetLogy();398 //gPad->SetLogy(); 396 399 397 400 MH::MakeDefCanvas(); … … 399 402 angle.SetYTitle("Angle: \\Theta [\\circ]"); 400 403 angle.DrawCopy("surf1pol"); 401 gPad->SetLogy();404 //gPad->SetLogy(); 402 405 403 406 MH::MakeDefCanvas();
Note:
See TracChangeset
for help on using the changeset viewer.