Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * 21 // * Parts of this code which have been devel 22 // * et al under contract (31-465) to the King 23 // * Science and Technology (KACST), the Natio 24 // * Mathematics and Physics (NCMP), Saudi Ara 25 // * 26 // * By using, copying, modifying or distri 27 // * any work based on the software) you ag 28 // * use in resulting scientific publicati 29 // * acceptance of all terms of the Geant4 Sof 30 // ******************************************* 31 // 32 /// \file hadronic/Hadr02/src/G4UrQMD1_3Model. 33 /// \brief Implementation of the G4UrQMD1_3Mod 34 // 35 // 36 ////////////////////////////////////////////// 37 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 38 // 39 // MODULE: G4UrQMD1_3Model.hh 40 // 41 // Version: 0.B 42 // Date: 25/01/12 43 // Authors: Kh. Abdel-Waged and Nuha Fe 44 // Revised by: V.V. Uzhinskii 45 // SPONSERED BY 46 // Customer: KAUST/NCMP 47 // Contract: 31-465 48 // 49 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 50 // 51 #ifdef G4_USE_URQMD 52 53 # include "G4UrQMD1_3Model.hh" 54 55 # include "G4UrQMD1_3Interface.hh" 56 //------------------------------- 57 # include "G4CollisionOutput.hh" 58 # include "G4DynamicParticle.hh" 59 # include "G4IonTable.hh" 60 # include "G4LorentzRotation.hh" 61 # include "G4Nucleus.hh" 62 # include "G4ParticleDefinition.hh" 63 # include "G4ParticleTable.hh" 64 # include "G4PhysicalConstants.hh" 65 # include "G4SystemOfUnits.hh" 66 # include "G4Track.hh" 67 # include "G4V3DNucleus.hh" 68 # include "globals.hh" 69 70 // AND-> 71 # include "G4Version.hh" 72 // AND<- 73 //----------------new_anti 74 # include "G4AntiAlpha.hh" 75 # include "G4AntiDeuteron.hh" 76 # include "G4AntiHe3.hh" 77 # include "G4AntiTriton.hh" 78 //--------------------------- 79 # include <fstream> 80 # include <string> 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 83 // 84 G4UrQMD1_3Model::G4UrQMD1_3Model(const G4Strin 85 : G4VIntraNuclearTransportModel(nam), verbos 86 { 87 if (verbose > 3) { 88 G4cout << " >>> G4UrQMD1_3Model default co 89 } 90 91 // 92 // Set the minimum and maximum range for the 93 94 // SetMinEnergy(100.0*MeV); 95 // SetMaxEnergy(200.0*GeV); 96 97 // 98 99 // 100 WelcomeMessage(); 101 // 102 CurrentEvent = 0; 103 // 104 105 InitialiseDataTables(); 106 107 // 108 } 109 110 //....oooOO0OOooo........oooOO0OOooo........oo 111 // Destructor 112 // 113 G4UrQMD1_3Model::~G4UrQMD1_3Model() {} 114 115 //....oooOO0OOooo........oooOO0OOooo........oo 116 117 G4ReactionProductVector* G4UrQMD1_3Model::Prop 118 { 119 return 0; 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oo 123 // 124 // ApplyYourself 125 // 126 // Member function to process an event, and ge 127 128 G4HadFinalState* G4UrQMD1_3Model::ApplyYoursel 129 130 { 131 // anti_new 132 // ------------------define anti_light_nuc 133 const G4ParticleDefinition* anti_deu = G4Ant 134 135 const G4ParticleDefinition* anti_he3 = G4Ant 136 137 const G4ParticleDefinition* anti_tri = G4Ant 138 139 const G4ParticleDefinition* anti_alp = G4Ant 140 //------------------------------------------ 141 // 142 // The secondaries will be returned in G4Had 143 // initialise this. The original track will 144 // secondaries followed. 145 // 146 theResult.Clear(); 147 theResult.SetStatusChange(stopAndKill); 148 149 G4DynamicParticle* cascadeParticle = 0; 150 // 151 // 152 // Get relevant information about the projec 153 // momentum, etc). 154 // 155 156 const G4ParticleDefinition* definitionP = th 157 const G4double AP = definitionP->GetBaryonNu 158 const G4double ZP = definitionP->GetPDGCharg 159 G4double AT = theTarget.GetN(); 160 G4double ZT = theTarget.GetZ(); 161 // ---------------------------------------- 162 G4int id = definitionP->GetPDGEncoding(); / 163 // ----------------------------------------- 164 G4int AP1 = G4lrint(AP); 165 G4int ZP1 = G4lrint(ZP); 166 G4int AT1 = G4lrint(AT); 167 G4int ZT1 = G4lrint(ZT); 168 // G4cout<<"------ap1--=="<<AP1<<"---zp1--- 169 // 170 // ***************************************** 171 // The following is the parameters necessary 172 // ----------------------------------------- 173 urqmdparams_.u_sptar = 0; //! 0= normal pro 174 urqmdparams_.u_spproj = 1; // projectile is 175 176 // new_anti 177 178 if (AP1 > 1 || definitionP == anti_deu || de 179 || definitionP == anti_alp) 180 { 181 urqmdparams_.u_ap = AP1; 182 urqmdparams_.u_zp = ZP1; 183 184 urqmdparams_.u_spproj = 0; 185 } 186 else if (id == 2212) { //! proton 187 urqmdparams_.u_ap = 1; 188 urqmdparams_.u_zp = 1; 189 } 190 else if (id == -2212) { //! anti-proton 191 urqmdparams_.u_ap = -1; 192 urqmdparams_.u_zp = -1; 193 } 194 else if (id == 2112) { //! neutron 195 urqmdparams_.u_ap = 1; 196 urqmdparams_.u_zp = -1; 197 } 198 else if (id == -2112) { //! anti-neutron 199 urqmdparams_.u_ap = -1; 200 urqmdparams_.u_zp = 1; 201 } 202 else if (id == 211) { //! pi+ 203 urqmdparams_.u_ap = 101; 204 urqmdparams_.u_zp = 2; 205 } 206 else if (id == -211) { //! pi- 207 urqmdparams_.u_ap = 101; 208 urqmdparams_.u_zp = -2; 209 } 210 else if (id == 321) { //! K+ 211 urqmdparams_.u_ap = 106; 212 urqmdparams_.u_zp = 1; 213 } 214 else if (id == -321) { //! K- 215 urqmdparams_.u_ap = -106; 216 urqmdparams_.u_zp = -1; 217 } 218 else if (id == 130 || id == 310) { // ! K0 219 urqmdparams_.u_ap = 106; 220 urqmdparams_.u_zp = -1; 221 } 222 else if (id == -130 || id == -310) { // ! K 223 urqmdparams_.u_ap = -106; 224 urqmdparams_.u_zp = 1; 225 } 226 else { 227 G4cout << " Sorry, No definition for parti 228 229 // AND-> 230 # if G4VERSION_NUMBER >= 950 231 // New signature (9.5) for G4Exception 232 // Using G4HadronicException 233 throw G4HadronicException(__FILE__, __LINE 234 # else 235 G4Exception(" "); 236 # endif 237 // AND<- 238 } // end if id 239 //------------------------------------------ 240 241 urqmdparams_.u_at = AT1; // Target identifi 242 urqmdparams_.u_zt = ZT1; 243 //------------------------------------------ 244 // identify Energy 245 // 246 G4ThreeVector Pbefore = theTrack.Get4Momentu 247 G4double T = theTrack.GetKineticEnergy(); 248 G4double E = theTrack.GetTotalEnergy(); 249 G4double TotalEbefore = E * AP1 + theTarget. 250 // -------------------------------------- 251 252 if (AP1 > 1) { 253 urqmdparams_.u_elab = T / (AP1 * GeV); // 254 255 E = E / AP1; // Units are GeV/nuc 256 } 257 else { 258 urqmdparams_.u_elab = T / GeV; // units a 259 260 TotalEbefore = E + theTarget.AtomicMass(AT 261 } 262 263 //------------------------------------------ 264 // identify impact parameter 265 urqmdparams_.u_imp = -(1.1 * std::pow(G4doub 266 // units are in fm for UrQMD; 267 //------------------------------------------ 268 ///////////////////////// initialise//////// 269 270 if (CurrentEvent == 0) { 271 G4cout << "\n creation of table, wait----- 272 273 G4cout << "\n" << G4endl; 274 275 G4int io = 0; 276 277 uinit_(&io); 278 279 G4cout << "\n end to create table " << G4 280 281 CurrentEvent = 1; 282 } 283 //////////////////////////////////////////// 284 285 // #ifdef debug_G4UrQMD1_3Model 286 287 G4cout << "UrQMDModel running-------------" 288 289 urqmd_(); 290 291 // #endif 292 293 // G4cout <<"Number of produced particles: 294 295 G4int n = sys_.npart; // no of produced par 296 if (n < 2) { 297 G4cout << "===============Warning========= 298 G4cout << "=============================== 299 300 G4cout << "Number of produced particles is 301 G4cout << "=============================== 302 303 // AND-> 304 # if G4VERSION_NUMBER >= 950 305 // New signature (9.5) for G4Exception 306 // Using G4HadronicException instead of ba 307 throw G4HadronicException(__FILE__, __LINE 308 # else 309 G4Exception(" "); // stop 310 # endif 311 // AND<- 312 } 313 else { 314 for (G4int i = 0; i < n; i++) { 315 G4int pid = pdgid_(&isys_.ityp[i], &isys 316 317 // Particle is a final state secondary a 318 // Determine what this secondary particl 319 // parameters. 320 // 321 322 G4ParticleDefinition* pd = G4ParticleTab 323 324 if (pd) { 325 G4double px = (coor_.px[i] + ffermi_.f 326 // units are in MeV/c for G4 327 G4double py = (coor_.py[i] + ffermi_.f 328 G4double pz = (coor_.pz[i] + ffermi_.f 329 330 G4double et = (coor_.p0[i]) * GeV; 331 332 // ------------------------------Us 333 334 G4LorentzVector lorenzvec = G4LorentzV 335 336 cascadeParticle = new G4DynamicParticl 337 338 theResult.AddSecondary(cascadeParticle 339 340 //==================================== 341 342 } // if 343 } // for 344 345 } // if warning 346 347 //========================================== 348 if (verbose >= 3) { 349 // 350 G4double TotalEafter = 0.0; 351 G4ThreeVector TotalPafter; 352 G4double charge = 0.0; 353 G4int baryon = 0; 354 G4int nSecondaries = theResult.GetNumberOf 355 356 for (G4int j = 0; j < nSecondaries; j++) { 357 TotalEafter += theResult.GetSecondary(j) 358 359 TotalPafter += theResult.GetSecondary(j) 360 361 G4ParticleDefinition* pd = theResult.Get 362 363 charge += pd->GetPDGCharge(); 364 baryon += pd->GetBaryonNumber(); 365 366 } // for secondaries 367 368 G4cout << "------------------------------- 369 << "------------------------------- 370 G4cout << "Total energy before collision 371 << " MeV" << G4endl; 372 G4cout << "Total energy after collision 373 << " MeV" << G4endl; 374 375 G4cout << "------------------------------- 376 377 G4cout << "Total momentum before collision 378 << " MeV/c" << G4endl; 379 G4cout << "Total momentum after collision 380 << " MeV/c" << G4endl; 381 G4cout << "------------------------------- 382 383 if (verbose >= 4) { 384 G4cout << "Total charge before collision 385 G4cout << "Total charge after collision 386 387 G4cout << "----------------------------- 388 389 G4cout << "Total baryon number before co 390 G4cout << "Total baryon number after col 391 G4cout << "----------------------------- 392 393 } // if verbose4 394 395 G4cout << "------------------------------- 396 << "------------------------------- 397 398 } // if verbose3 399 400 return &theResult; 401 } // G4hadfinal 402 403 //....oooOO0OOooo........oooOO0OOooo........oo 404 // 405 // WelcomeMessage 406 // 407 void G4UrQMD1_3Model::WelcomeMessage() const 408 { 409 G4cout << G4endl; 410 G4cout << " ******************************** 411 G4cout << " Interface to G4UrQMD_1.3 412 G4cout << " Version number : 00.00.0B 413 G4cout << " (Interface written by Kh. Abdel- 414 G4cout << G4endl; 415 G4cout << " ******************************** 416 G4cout << G4endl; 417 418 return; 419 } 420 421 //....oooOO0OOooo........oooOO0OOooo........oo 422 423 void G4UrQMD1_3Model::InitialiseDataTables() 424 { 425 // 426 // 427 // The next line is to make sure the block d 428 // executed. 429 // 430 431 g4urqmdblockdata_(); 432 433 //////////////////////////////////////////// 434 /////// Dynamic seed /////////////////////// 435 // G4int ranseed=-time_ (); 436 // Fixed seed ///////////////////////// 437 438 G4int ranseed = 1097569630; 439 440 G4cout << "\n seed: " << ranseed << G4endl; 441 442 sseed_(&ranseed); 443 444 loginit_(); 445 } 446 447 #endif // G4_USE_URQMD 448