- Timestamp:
- Dec 22, 2010, 3:52:27 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/incl/src/G4InclCascadeInterface.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4InclCascadeInterface.cc,v 1.1 2 2010/10/26 02:47:59 kaitanie Exp $26 // $Id: G4InclCascadeInterface.cc,v 1.15 2010/11/17 20:19:09 kaitanie Exp $ 27 27 // Translation of INCL4.2/ABLA V3 28 28 // Pekka Kaitaniemi, HIP (translation) … … 34 34 35 35 #include "G4InclCascadeInterface.hh" 36 #include "G4FermiBreakUp.hh" 36 37 #include "math.h" 37 38 #include "G4GenericIon.hh" 38 39 #include "CLHEP/Random/Random.h" 39 40 40 G4InclCascadeInterface::G4InclCascadeInterface() 41 G4InclCascadeInterface::G4InclCascadeInterface(const G4String& nam) 42 :G4VIntraNuclearTransportModel(nam) 41 43 { 42 44 hazard = new G4Hazard(); … … 45 47 46 48 varntp = new G4VarNtp(); 47 inclInput= 0;49 calincl = 0; 48 50 ws = new G4Ws(); 49 51 mat = new G4Mat(); 50 incl = new G4Incl(hazard, inclInput, ws, mat, varntp); 52 incl = new G4Incl(hazard, calincl, ws, mat, varntp); 53 54 theExcitationHandler = new G4ExcitationHandler; 55 thePrecoModel = new G4PreCompoundModel(theExcitationHandler); 56 57 if(!getenv("G4INCLABLANOFERMIBREAKUP")) { // Use Fermi Break-up by default if it is NOT explicitly disabled 58 incl->setUseFermiBreakUp(true); 59 } 51 60 52 61 verboseLevel = 0; … … 55 64 G4InclCascadeInterface::~G4InclCascadeInterface() 56 65 { 66 delete thePrecoModel; 67 delete theExcitationHandler; 68 57 69 delete hazard; 58 70 delete varntp; 59 delete inclInput;60 71 delete ws; 61 72 delete mat; … … 68 79 69 80 G4int particleI; 81 G4int bulletType = 0; 70 82 71 83 // Print diagnostic messages: 0 = silent, 1 and 2 = verbose … … 91 103 G4DynamicParticle *cascadeParticle = 0; 92 104 G4ParticleDefinition *aParticleDefinition = 0; 105 106 G4ReactionProductVector *thePrecoResult = 0; 107 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable(); 93 108 94 109 // INCL assumes the projectile particle is going in the direction of … … 102 117 G4LorentzRotation toLabFrame = toZ.inverse(); 103 118 104 inclInput = new G4InclInput(aTrack, theNucleus, false);105 106 119 theResult.Clear(); // Make sure the output data structure is clean. 120 121 calincl = new G4InclInput(aTrack, theNucleus, false); 122 incl->setInput(calincl); 123 124 // G4InclInput::printProjectileTargetInfo(aTrack, theNucleus); 125 // calincl->printInfo(); 107 126 108 127 #ifdef DEBUGINCL … … 117 136 G4double eKinSum = bulletE; 118 137 G4LorentzVector labv = G4LorentzVector(0.0, 0.0, std::sqrt(bulletE*(bulletE + 2.*mass)), bulletE + mass + amass); 138 G4LorentzVector labvA = G4LorentzVector(0.0, 0.0, 0.0, 0.0); 119 139 G4cout <<"Energy in the beginning = " << labv.e() / MeV << G4endl; 120 140 #endif 121 141 122 142 // Check wheter the input is acceptable. 123 if(( inclInput->bulletType() != 0) && ((inclInput->targetA() != 1) && (inclInput->targetZ() != 1))) {143 if((calincl->bulletType() != 0) && ((calincl->targetA() != 1) && (calincl->targetZ() != 1))) { 124 144 ws->nosurf = -2; // Nucleus surface, -2 = Woods-Saxon 125 145 ws->xfoisa = 8; … … 130 150 131 151 mat->nbmat = 1; 132 mat->amat[0] = int( inclInput->targetA());133 mat->zmat[0] = int( inclInput->targetZ());152 mat->amat[0] = int(calincl->targetA()); 153 mat->zmat[0] = int(calincl->targetZ()); 134 154 135 155 incl->initIncl(true); … … 140 160 G4cout <<"G4InclCascadeInterface: Try number = " << nTries << G4endl; 141 161 } 142 incl->processEventIncl( inclInput);162 incl->processEventIncl(calincl); 143 163 144 164 if(verboseLevel > 1) { … … 151 171 * Diagnostic output 152 172 */ 153 G4cout <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;154 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;155 156 G4cout <<"G4InclCascadeInterface: Target A: " << inclInput->targetA() << G4endl;157 G4cout <<"G4InclCascadeInterface: Target Z: " << inclInput->targetZ() << G4endl;173 G4cout <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl; 174 G4cout <<"G4Incl4AblaCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 175 176 G4cout <<"G4InclCascadeInterface: Target A: " << calincl->targetA() << G4endl; 177 G4cout <<"G4InclCascadeInterface: Target Z: " << calincl->targetZ() << G4endl; 158 178 159 179 if(verboseLevel > 3) { 160 diagdata <<"G4InclCascadeInterface: Bullet type: " << inclInput->bulletType() << G4endl;161 diagdata <<"G4InclCascadeInterface: Bullet energy: " << inclInput->bulletE() << " MeV" << G4endl;180 diagdata <<"G4InclCascadeInterface: Bullet type: " << calincl->bulletType() << G4endl; 181 diagdata <<"G4InclCascadeInterface: Bullet energy: " << calincl->bulletE() << " MeV" << G4endl; 162 182 163 diagdata <<"G4InclCascadeInterface: Target A: " << inclInput->targetA() << G4endl; 164 diagdata <<"G4InclCascadeInterface: Target Z: " << inclInput->targetZ() << G4endl; 165 } 166 167 // for(particleI = 0; particleI < varntp->ntrack; particleI++) { 168 // G4cout << n << " " << inclInput->f[6] << " " << inclInput->f[2] << " "; 169 // G4cout << varntp->massini << " " << varntp->mzini << " "; 170 // G4cout << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 171 // G4cout << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 172 // G4cout << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " " << varntp->itypcasc[particleI] << " "; 173 // G4cout << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 174 // G4cout << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 175 // // For diagnostic output 176 // if(verboseLevel > 3) { 177 // diagdata << n << " " << inclInput->f[6] << " " << inclInput->f[2] << " "; 178 // diagdata << varntp->massini << " " << varntp->mzini << " "; 179 // diagdata << varntp->exini << " " << varntp->mulncasc << " " << varntp->mulnevap << " " << varntp->mulntot << " "; 180 // diagdata << varntp->bimpact << " " << varntp->jremn << " " << varntp->kfis << " " << varntp->estfis << " "; 181 // diagdata << varntp->izfis << " " << varntp->iafis << " " << varntp->ntrack << " "; 182 // diagdata << varntp->itypcasc[particleI] << " "; 183 // diagdata << varntp->avv[particleI] << " " << varntp->zvv[particleI] << " " << varntp->enerj[particleI] << " "; 184 // diagdata << varntp->plab[particleI] << " " << varntp->tetlab[particleI] << " " << varntp->philab[particleI] << G4endl; 185 // } 186 // } 183 diagdata <<"G4InclCascadeInterface: Target A: " << calincl->targetA() << G4endl; 184 diagdata <<"G4InclCascadeInterface: Target Z: " << calincl->targetZ() << G4endl; 185 } 186 187 187 } 188 188 … … 196 196 197 197 theResult.SetStatusChange(stopAndKill); 198 aParticleDefinition = G4InclInput::getParticleDefinition(inclInput->bulletType()); 199 200 cascadeParticle = new G4DynamicParticle(); 201 cascadeParticle->SetDefinition(aParticleDefinition); 202 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 203 theResult.AddSecondary(cascadeParticle); 198 199 G4int bulletType = calincl->bulletType(); 200 aParticleDefinition = G4InclInput::getParticleDefinition(bulletType); 201 202 if(aParticleDefinition != 0) { 203 cascadeParticle = new G4DynamicParticle(); 204 cascadeParticle->SetDefinition(aParticleDefinition); 205 cascadeParticle->Set4Momentum(aTrack.Get4Momentum()); 206 theResult.AddSecondary(cascadeParticle); 207 } 204 208 } 205 209 … … 209 213 210 214 #ifdef DEBUGINCL 211 G4cout << "E [MeV]" << std::setw(12) << " Ekin [MeV]" << std::setw(12) << " E* [MeV]" << std::setw(12) << "Px [MeV]" << std::setw(12) << " Py [MeV]" << std::setw(12) << "Pz [MeV]" << std::setw(12) << "Pt [MeV]" << std::setw(12) << "A" << std::setw(12) << "Z" << G4endl; 215 G4cout << "E [MeV]" << std::setw(12) 216 << " Ekin [MeV]" << std::setw(12) 217 << "Px [MeV]" << std::setw(12) 218 << " Py [MeV]" << std::setw(12) 219 << "Pz [MeV]" << std::setw(12) 220 << "Pt [MeV]" << std::setw(12) 221 << "A" << std::setw(12) 222 << "Z" << G4endl; 212 223 #endif 213 224 … … 276 287 if((varntp->avv[particleI] > 1) && (varntp->zvv[particleI] >= 1)) { // Nucleus fragment 277 288 G4ParticleDefinition * aIonDef = 0; 278 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();279 289 280 290 G4int A = G4int(varntp->avv[particleI]); … … 315 325 G4double p = mom.mag(); 316 326 labv -= fm; 317 G4double px = mom.x() * MeV; 318 G4double py = mom.y() * MeV; 319 G4double pz = mom.z() * MeV; 327 if(varntp->avv[particleI] > 1) { 328 labvA += fm; 329 } 330 G4double px = mom.x() * MeV; 331 G4double py = mom.y() * MeV; 332 G4double pz = mom.z() * MeV; 320 333 G4double pt = std::sqrt(px*px+py*py); 321 334 G4double e = fm.e(); … … 330 343 G4cout << fm.e() / MeV 331 344 << std::setw(12) << cascadeParticle->GetKineticEnergy() / MeV 332 << std::setw(12) << exE / MeV333 345 << std::setw(12) << mom.x() / MeV 334 346 << std::setw(12) << mom.y() / MeV … … 351 363 } 352 364 } 365 366 G4double nuclearMass = G4NucleiProperties::GetNuclearMass(G4int(varntp->massini), G4int(varntp->mzini)) + varntp->exini * MeV; 367 G4LorentzVector fragmentMomentum(varntp->pxrem * MeV, varntp->pyrem * MeV, varntp->pzrem * MeV, 368 varntp->erecrem * MeV + nuclearMass); 369 G4double momentumScaling = G4InclUtils::calculate4MomentumScaling(G4int(varntp->massini), G4int(varntp->mzini), 370 varntp->exini, 371 varntp->erecrem, 372 varntp->pxrem, 373 varntp->pyrem, 374 varntp->pzrem); 375 G4LorentzVector p4(momentumScaling * varntp->pxrem * MeV, momentumScaling * varntp->pyrem * MeV, 376 momentumScaling * varntp->pzrem * MeV, 377 varntp->erecrem + nuclearMass); 378 379 // For four-momentum, baryon number and charge conservation check: 380 G4LorentzVector fourMomentumBalance = p4; 381 G4int baryonNumberBalance = G4int(varntp->massini); 382 G4int chargeBalance = G4int(varntp->mzini); 383 384 G4LorentzRotation toFragmentZ; 385 toFragmentZ.rotateZ(-p4.theta()); 386 toFragmentZ.rotateY(-p4.phi()); 387 G4LorentzRotation toFragmentLab = toFragmentZ.inverse(); 388 p4 *= toFragmentZ; 389 390 G4LorentzVector p4rest = p4; 391 p4rest.boost(-p4.boostVector()); 392 if(verboseLevel > 0) { 393 G4cout <<"Cascade remnant nucleus:" << G4endl; 394 G4cout <<"p4: " << G4endl; 395 G4cout <<" px: " << p4.px() <<" py: " << p4.py() <<" pz: " << p4.pz() << G4endl; 396 G4cout <<" E = " << p4.e() << G4endl; 397 398 G4cout <<"p4rest: " << G4endl; 399 G4cout <<" px: " << p4rest.px() <<" py: " << p4rest.py() <<" pz: " << p4rest.pz() << G4endl; 400 G4cout <<" E = " << p4rest.e() << G4endl; 401 } 402 403 G4Fragment theCascadeRemnant(G4int(varntp->massini), G4int(varntp->mzini), p4rest); 404 thePrecoResult = thePrecoModel->DeExcite(theCascadeRemnant); 405 if(thePrecoResult != 0) { 406 G4ReactionProductVector::iterator fragment; 407 for(fragment = thePrecoResult->begin(); fragment != thePrecoResult->end(); fragment++) { 408 G4ParticleDefinition *theFragmentDefinition = (*fragment)->GetDefinition(); 409 410 if(theFragmentDefinition != 0) { 411 G4DynamicParticle *theFragment = new G4DynamicParticle(theFragmentDefinition, (*fragment)->GetMomentum()); 412 G4LorentzVector labMomentum = theFragment->Get4Momentum(); 413 labMomentum.boost(p4.boostVector()); 414 labMomentum *= toFragmentLab; 415 labMomentum *= toLabFrame; 416 theFragment->Set4Momentum(labMomentum); 417 fourMomentumBalance -= theFragment->Get4Momentum(); 418 baryonNumberBalance -= theFragmentDefinition->GetAtomicMass(); 419 chargeBalance -= theFragmentDefinition->GetAtomicNumber(); 420 if(verboseLevel > 0) { 421 G4cout <<"Resulting fragment: " << G4endl; 422 G4cout <<" kinetic energy = " << theFragment->GetKineticEnergy() / MeV << " MeV" << G4endl; 423 G4cout <<" momentum = " << theFragment->GetMomentum().mag() / MeV << " MeV" << G4endl; 424 } 425 theResult.AddSecondary(theFragment); 426 } else { 427 G4cout <<"G4InclCascadeInterface: Error. Fragment produced by Fermi break-up does not exist." << G4endl; 428 G4cout <<"Resulting fragment: " << G4endl; 429 G4cout <<" momentum = " << (*fragment)->GetMomentum().mag() / MeV << " MeV" << G4endl; 430 } 431 } 432 delete thePrecoResult; 433 thePrecoResult = 0; 434 435 if(verboseLevel > 1 && std::abs(fourMomentumBalance.mag() / MeV) > 0.1 * MeV) { 436 G4cout <<"Four-momentum balance after remnant nucleus Fermi break-up:" << G4endl; 437 G4cout <<"Magnitude: " << fourMomentumBalance.mag() / MeV << " MeV" << G4endl; 438 G4cout <<"Vector components (px, py, pz, E) = (" 439 << fourMomentumBalance.px() << ", " 440 << fourMomentumBalance.py() << ", " 441 << fourMomentumBalance.pz() << ", " 442 << fourMomentumBalance.e() << ")" << G4endl; 443 } 444 if(baryonNumberBalance != 0 && verboseLevel > 1) { 445 G4cout <<"Baryon number balance after remnant nucleus Fermi break-up: " << baryonNumberBalance << G4endl; 446 } 447 if(chargeBalance != 0 && verboseLevel > 1) { 448 G4cout <<"Charge balance after remnant nucleus Fermi break-up: " << chargeBalance << G4endl; 449 } 450 } 451 // } // if(needsFermiBreakUp) 452 353 453 #ifdef DEBUGINCL 354 454 G4cout <<"--------------------------------------------------------------------------------" << G4endl; 355 455 G4double pt = std::sqrt(std::pow(labv.x(), 2) + std::pow(labv.y(), 2)); 356 G4cout << labv.e() / MeV << std::setw(12) << eKinSum / MeV << std::setw(12) << labv.x() << std::setw(12) << labv.y() << std::setw(12) << labv.z() << std::setw(12) << pt / MeV << std::setw(12) << baryonNumber << std::setw(12) << chargeNumber << " totals" << G4endl; 456 G4double ptA = std::sqrt(std::pow(labvA.x(), 2) + std::pow(labvA.y(), 2)); 457 G4cout << labv.e() / MeV << std::setw(12) 458 << eKinSum / MeV << std::setw(12) 459 << labv.x() / MeV << std::setw(12) 460 << labv.y() / MeV << std::setw(12) 461 << labv.z() / MeV << std::setw(12) 462 << pt / MeV << std::setw(12) 463 << baryonNumber << std::setw(12) 464 << chargeNumber << " totals" << G4endl; 465 G4cout << " - " << std::setw(12) 466 << " - " << std::setw(12) 467 << labvA.x() / MeV << std::setw(12) 468 << labvA.y() / MeV << std::setw(12) 469 << labvA.z() / MeV << std::setw(12) 470 << ptA / MeV << std::setw(12) 471 << " - " << std::setw(12) << " - " << " totals ABLA" << G4endl; 357 472 G4cout << G4endl; 358 473 359 474 if(verboseLevel > 3) { 360 475 if(baryonNumber != 0) { … … 369 484 } 370 485 #endif 371 486 372 487 varntp->ntrack = 0; // Clean up the number of generated particles in the event. 373 488 } … … 379 494 theResult.SetStatusChange(stopAndKill); 380 495 381 G4ParticleTable *theTableOfParticles = G4ParticleTable::GetParticleTable();382 496 cascadeParticle = new G4DynamicParticle(theTableOfParticles->FindParticle(aTrack.GetDefinition()), aTrack.Get4Momentum()); 383 497 … … 388 502 } 389 503 if(verboseLevel > 3) { 390 diagdata <<"ERROR G4InclCascadeInterface: Processing event number (internal) failed " << eventNumber << G4endl;391 } 392 393 if( inclInput->bulletType()== 0) {504 diagdata <<"ERROR G4InclCascadeInterface: Error processing event number (internal) failed " << eventNumber << G4endl; 505 } 506 507 if(bulletType == 0) { 394 508 if(verboseLevel > 1) { 395 509 G4cout <<"G4InclCascadeInterface: Unknown bullet type" << G4endl; … … 402 516 } 403 517 404 if(( inclInput->targetA() == 1) && (inclInput->targetZ() == 1)) { // Unsupported target518 if((calincl->targetA() == 1) && (calincl->targetZ() == 1)) { // Unsupported target 405 519 if(verboseLevel > 1) { 406 520 G4cout <<"Unsupported target: " << G4endl; 407 G4cout <<"Target A: " << inclInput->targetA() << G4endl;408 G4cout <<"TargetZ: " << inclInput->targetZ() << G4endl;521 G4cout <<"Target A: " << calincl->targetA() << G4endl; 522 G4cout <<"TargetZ: " << calincl->targetZ() << G4endl; 409 523 } 410 524 if(verboseLevel > 3) { 411 525 diagdata <<"Unsupported target: " << G4endl; 412 diagdata <<"Target A: " << inclInput->targetA() << G4endl;413 diagdata <<"TargetZ: " << inclInput->targetZ() << G4endl;414 } 415 } 416 417 if( inclInput->bulletE() < 100) { // INCL does not support E < 100 MeV.526 diagdata <<"Target A: " << calincl->targetA() << G4endl; 527 diagdata <<"TargetZ: " << calincl->targetZ() << G4endl; 528 } 529 } 530 531 if(calincl->bulletE() < 100) { // INCL does not support E < 100 MeV. 418 532 if(verboseLevel > 1) { 419 G4cout <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;533 G4cout <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 420 534 G4cout <<"WARNING: Returning the original bullet with original energy back to Geant4." << G4endl; 421 535 } 422 536 if(verboseLevel > 3) { 423 diagdata <<"Unsupported bullet energy: " << inclInput->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl;537 diagdata <<"Unsupported bullet energy: " << calincl->bulletE() << " MeV. (Lower limit is 100 MeV)." << G4endl; 424 538 } 425 539 } … … 430 544 } 431 545 546 delete calincl; 547 calincl = 0; 432 548 return &theResult; 433 549 }
Note: See TracChangeset
for help on using the changeset viewer.