- Timestamp:
- 06/20/02 16:24:20 (22 years ago)
- Location:
- trunk/WuerzburgSoft/Thomas/mphys
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/WuerzburgSoft/Thomas/mphys/Changelog
r1372 r1373 8 8 * MParticle.cc: 9 9 - made TRandom static (InitRandom) 10 11 * mphys.C: 12 - added support for Emis and Esum 10 13 11 14 -
trunk/WuerzburgSoft/Thomas/mphys/phys.C
r1370 r1373 118 118 void DoIt() 119 119 { 120 Double_t R = 500; //MParticle::RofZ(&startz); // [kpc]120 Double_t R = 100; //MParticle::RofZ(&startz); // [kpc] 121 121 Double_t startz = MParticle::ZofR(&R); 122 122 123 const char *filename = " cascade_500kpc_21d_B1e-18_03.root";124 125 const Double_t B = 1e-18;123 const char *filename = "data/cascade_delme.root"; 124 125 const Double_t B = 0; 126 126 127 127 cout << "R = " << R << "kpc" << endl; … … 131 131 MPairProduction pair; 132 132 133 Double_t runtime = 1; // [s]133 Double_t runtime = 60; // [s] 134 134 135 135 Double_t lo = 1e4; … … 246 246 cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush; 247 247 248 Double_t Esum=0; 249 Double_t Emis=0; 250 248 251 while (1) 249 252 { … … 270 273 p->SetIsPrimary(kFALSE); 271 274 T1->Fill(); 275 276 Esum += Eg; 272 277 273 278 delete listg.Remove(p); … … 335 340 { 336 341 cout << "!" << flush; 342 Esum += e->GetEnergy(); 343 Emis += e->GetEnergy(); 337 344 break; 338 345 } 339 346 340 347 // WRONG! 341 Double_t theta = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;348 Double_t theta = 0; //rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4; 342 349 MPhoton *p = e->DoInvCompton(theta); 343 350 … … 357 364 cout << "t" << flush; // ignored 358 365 delete p; 366 Esum += p->GetEnergy(); 367 Emis += p->GetEnergy(); 359 368 } 360 369 361 if (fabs(e->GetTheta()*3437)<60 && // < 60min 362 e->GetEnergy()>5*lo) 363 continue; 364 365 cout << "0" << flush; 366 break; 370 if (fabs(e->GetTheta()*3437)>60) // < 60min 371 { 372 cout << "T" << flush; 373 Esum += e->GetEnergy(); 374 Emis += e->GetEnergy(); 375 break; 376 } 377 378 if (e->GetEnergy()<lo)//*5) 379 { 380 cout << "E" << flush; 381 Esum += e->GetEnergy(); 382 Emis += e->GetEnergy(); 383 break; 384 } 367 385 } 368 386 } 369 387 liste.Delete(); 370 388 } 371 cout << endl;389 cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl; 372 390 373 391 Int_t now = (int)(TStopwatch::GetRealTime()- starttime)/5;
Note:
See TracChangeset
for help on using the changeset viewer.