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 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 #include <algorithm> 27 #include <vector> 28 #include <cmath> 29 #include <numeric> 30 31 #include "G4BinaryLightIonReaction.hh" 32 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4LorentzVector.hh" 35 #include "G4LorentzRotation.hh" 36 #include "G4ReactionProductVector.hh" 37 #include "G4ping.hh" 38 #include "G4Delete.hh" 39 #include "G4Neutron.hh" 40 #include "G4VNuclearDensity.hh" 41 #include "G4FermiMomentum.hh" 42 #include "G4HadTmpUtil.hh" 43 #include "G4PreCompoundModel.hh" 44 #include "G4HadronicInteractionRegistry.hh" 45 #include "G4Log.hh" 46 #include "G4PhysicsModelCatalog.hh" 47 #include "G4HadronicParameters.hh" 48 49 G4int G4BinaryLightIonReaction::theBLIR_ID = - 50 51 //#define debug_G4BinaryLightIonReaction 52 //#define debug_BLIR_finalstate 53 //#define debug_BLIR_result 54 55 G4BinaryLightIonReaction::G4BinaryLightIonReac 56 : G4HadronicInteraction("Binary Light Ion Casc 57 theProjectileFragmentation(ptr), 58 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spect 59 projectile3dNucleus(0),target3dNucleus(0) 60 { 61 if(!ptr) { 62 G4HadronicInteraction* p = 63 G4HadronicInteractionRegistry::Instance( 64 G4VPreCompoundModel* pre = static_cast<G4V 65 if(!pre) { pre = new G4PreCompoundModel(); 66 theProjectileFragmentation = pre; 67 } 68 theModel = new G4BinaryCascade(theProjectile 69 theHandler = theProjectileFragmentation->Get 70 theBLIR_ID = G4PhysicsModelCatalog::GetM 71 debug_G4BinaryLightIonReactionResults = G4Ha 72 } 73 74 G4BinaryLightIonReaction::~G4BinaryLightIonRea 75 {} 76 77 void G4BinaryLightIonReaction::ModelDescriptio 78 { 79 outFile << "G4Binary Light Ion Cascade is an 80 << "using G4BinaryCasacde to model the i 81 << "nucleus with a nucleus.\n" 82 << "The lighter of the two nuclei is tre 83 << "which are transported simultaneously 84 } 85 86 //-------------------------------------------- 87 struct ReactionProduct4Mom 88 { 89 G4LorentzVector operator()(G4LorentzVector 90 }; 91 92 G4HadFinalState *G4BinaryLightIonReaction:: 93 ApplyYourself(const G4HadProjectile &aTrack, G 94 { 95 if(debug_G4BinaryLightIonReactionResults) G4 96 G4ping debug("debug_G4BinaryLightIonReaction 97 pA=aTrack.GetDefinition()->GetBaryonNumber() 98 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCha 99 tA=targetNucleus.GetA_asInt(); 100 tZ=targetNucleus.GetZ_asInt(); 101 G4double timePrimary = aTrack.GetGlobalTime( 102 G4LorentzVector mom(aTrack.Get4Momentum()); 103 //G4cout << "proj mom : " << mom << G4endl; 104 G4LorentzRotation toBreit(mom.boostVector()) 105 106 G4bool swapped=SetLighterAsProjectile(mom, t 107 //G4cout << "after swap, swapped? / mom " < 108 G4ReactionProductVector * result = 0; 109 G4ReactionProductVector * cascaders=0; //new 110 // G4double m_nucl(0); // to check energ 111 112 113 // G4double m1=G4ParticleTable::GetPartic 114 // G4cout << "Entering the decision point 115 // << (mom.t()-mom.mag())/pA << " 116 // << pA<<" "<< pZ<<" " 117 // << tA<<" "<< tZ<<G4endl 118 // << " "<<mom.t()-mom.mag()<<" " 119 // << mom.t()- m1<<G4endl; 120 if( (mom.t()-mom.mag())/pA < 50*MeV ) 121 { 122 // G4cout << "Using pre-compound only 123 // m_nucl = mom.mag(); 124 cascaders=FuseNucleiAndPrompound(mom); 125 if( !cascaders ) 126 { 127 128 // abort!! happens for too low e 129 130 theResult.Clear(); 131 theResult.SetStatusChange(isAliv 132 theResult.SetEnergyChange(aTrack 133 theResult.SetMomentumChange(aTra 134 return &theResult; 135 } 136 } 137 else 138 { 139 result=Interact(mom,toBreit); 140 141 if(! result ) 142 { 143 // abort!! 144 145 G4cerr << "G4BinaryLightIonReaction 146 G4cerr << " Primary " << aTrack.Get 147 << ", (A,Z)=(" << aTrack.GetDefi 148 << "," << aTrack.GetDefinition() 149 << ", kinetic energy " << aTrack 150 << G4endl; 151 G4cerr << " Target nucleus (A,Z)=(" 152 << (swapped?pA:tA) << "," 153 << (swapped?pZ:tZ) << ")" << 154 G4cerr << " if frequent, please sub 155 << G4endl << G4endl; 156 157 theResult.Clear(); 158 theResult.SetStatusChange(isAlive); 159 theResult.SetEnergyChange(aTrack.Ge 160 theResult.SetMomentumChange(aTrack. 161 return &theResult; 162 } 163 164 // Calculate excitation energy, 165 G4double theStatisticalExEnergy = GetProj 166 167 168 pInitialState = mom; 169 //G4cout << "BLIC: pInitialState from 170 pInitialState.setT(pInitialState.getT() + 171 G4ParticleTable::GetParticleTable()->GetIo 172 //G4cout << "BLIC: target nucleus added 173 174 delete target3dNucleus;target3dNucleus=0; 175 delete projectile3dNucleus;projectile3dNu 176 177 G4ReactionProductVector * spectators= new 178 179 cascaders = new G4ReactionProductVector; 180 181 G4LorentzVector pspectators=SortResult(re 182 // this also sets spectatorA and 183 184 // pFinalState=std::accumulate(casca 185 186 std::vector<G4ReactionProduct *>::iterato 187 188 // G4cout << "pInitialState, pFin 189 // if ( spectA-spectatorA !=0 || spec 190 // { 191 // G4cout << "spect Nucl != spect 192 // spectatorA <<" "<< spectatorZ << 193 // } 194 delete result; 195 result=0; 196 G4LorentzVector momentum(pInitialState-pF 197 G4int loopcount(0); 198 //G4cout << "BLIC: momentum, pspectato 199 while (std::abs(momentum.e()-pspectators. 200 201 { 202 G4LorentzVector pCorrect(pInitialState- 203 //G4cout << "BLIC:: BIC nonconservatio 204 // Correct outgoing casacde particles.. 205 G4bool EnergyIsCorrect=EnergyAndMomentu 206 if ( ! EnergyIsCorrect && debug_G4Binar 207 { 208 G4cout << "Warning - G4BinaryLightIon 209 } 210 pFinalState=G4LorentzVector(0,0,0,0); 211 for(iter=cascaders->begin(); iter!=casc 212 { 213 pFinalState += G4LorentzVector( (*ite 214 } 215 momentum=pInitialState-pFinalState; 216 if (++loopcount > 10 ) 217 { 218 break; 219 } 220 } 221 222 // Check if Energy/Momemtum is now ok, if 223 if ( std::abs(momentum.e()-pspectators.e( 224 { 225 for (iter=spectators->begin();iter!=specta 226 { 227 delete *iter; 228 } 229 delete spectators; 230 for(iter=cascaders->begin(); iter!=cascad 231 { 232 delete *iter; 233 } 234 delete cascaders; 235 236 G4cout << "G4BinaryLightIonReaction.cc: m 237 << " initial - final " << momentum 238 << " .. pInitialState/pFinalState/s 239 << pInitialState << G4endl 240 << pFinalState << G4endl 241 << pspectators << G4endl 242 << " .. A,Z " << spectatorA <<" "<< 243 G4cout << "G4BinaryLightIonReaction inval 244 G4cout << " Primary " << aTrack.GetDefini 245 << ", (A,Z)=(" << aTrack.GetDefin 246 << "," << aTrack.GetDefinition()- 247 << ", kinetic energy " << aTrack. 248 << G4endl; 249 G4cout << " Target nucleus (A,Z)=(" << t 250 << "," << targetNucleus.GetZ_ 251 G4cout << " if frequent, please submit ab 252 << G4endl << G4endl; 253 #ifdef debug_G4BinaryLightIonReaction 254 G4ExceptionDescription ed; 255 ed << "G4BinaryLightIonreaction: Ter 256 G4Exception("G4BinaryLightIonreactio 257 ed); 258 259 #endif 260 theResult.Clear(); 261 theResult.SetStatusChange(isAlive); 262 theResult.SetEnergyChange(aTrack.GetKinet 263 theResult.SetMomentumChange(aTrack.Get4Mo 264 return &theResult; 265 266 } 267 if (spectatorA > 0 ) 268 { 269 // DeExciteSpectatorNucleus() also 270 DeExciteSpectatorNucleus(specta 271 } else { // no spectators 272 delete spectators; 273 } 274 } 275 // Rotate to lab 276 G4LorentzRotation toZ; 277 toZ.rotateZ(-1*mom.phi()); 278 toZ.rotateY(-1*mom.theta()); 279 G4LorentzRotation toLab(toZ.inverse()); 280 281 // Fill the particle change, while rotating. 282 // theResult.Clear(); 283 theResult.Clear(); 284 theResult.SetStatusChange(stopAndKill); 285 G4LorentzVector ptot(0); 286 #ifdef debug_BLIR_result 287 G4LorentzVector p_raw; 288 #endif 289 //G4int i=0; 290 291 G4ReactionProductVector::iterator iter 292 for(iter=cascaders->begin(); iter!=cascaders 293 { 294 if((*iter)->GetNewlyAdded()) 295 { 296 G4DynamicParticle * aNewDP = 297 new G4DynamicParticle((*iter)->GetDe 298 (*iter)->GetTotalEnergy(), 299 (*iter)->GetMomentum() ); 300 G4LorentzVector tmp = aNewDP->Get4Moment 301 #ifdef debug_BLIR_result 302 p_raw+= tmp; 303 #endif 304 if(swapped) 305 { 306 tmp*=toBreit.inverse(); 307 tmp.setVect(-tmp.vect()); 308 } 309 tmp *= toLab; 310 aNewDP->Set4Momentum(tmp); 311 G4HadSecondary aNew = G4HadSecondary(aNe 312 G4double time = 0; 313 //if(time < 0.0) { time = 0.0; } 314 aNew.SetTime(timePrimary + time); 315 //aNew.SetCreatorModelID((*iter)-> 316 aNew.SetCreatorModelID(theBLIR_ID) 317 318 theResult.AddSecondary(aNew); 319 ptot += tmp; 320 //G4cout << "BLIC: Secondary " < 321 // <<" "<< aNew->GetMomen 322 } 323 delete *iter; 324 } 325 delete cascaders; 326 327 #ifdef debug_BLIR_result 328 //G4cout << "Result analysis, secondaries " 329 //G4cout << "p_tot_raw " << p_raw << " sum p 330 G4double m_nucl= G4ParticleTable::GetPartic 331 GetIonMass(targetNucleus.GetZ_asInt( 332 // delete? tZ=targetNucleus.GetZ_asInt(); 333 334 //G4cout << "BLIC Energy conservation initia 335 // << aTrack.GetTotalEnergy() + m_nuc 336 // <<" "<< aTrack.GetTotalEnergy() + m_ 337 G4cout << "BLIC momentum conservation " << a 338 << " ptot " << ptot << " delta " << aTra 339 << " 3mom.mag() " << (aTrack.Get4 340 #endif 341 342 if(debug_G4BinaryLightIonReactionResults) G4 343 344 return &theResult; 345 } 346 347 //-------------------------------------------- 348 349 //******************************************** 350 G4bool G4BinaryLightIonReaction::EnergyAndMome 351 G4ReactionProductVector* Output, G4Lorentz 352 //******************************************** 353 { 354 const int nAttemptScale = 2500; 355 const double ErrLimit = 1.E-6; 356 if (Output->empty()) 357 return TRUE; 358 G4LorentzVector SumMom(0,0,0,0); 359 G4double SumMass = 0; 360 G4double TotalCollisionMass = TotalCo 361 size_t i = 0; 362 // Calculate sum hadron 4-momenta and summin 363 for(i = 0; i < Output->size(); i++) 364 { 365 SumMom += G4LorentzVector((*Output)[i]->G 366 SumMass += (*Output)[i]->GetDefinition()-> 367 } 368 // G4cout << " E/P corrector, SumMass, Sum 369 // << SumMass <<" "<< SumMom.m2() << 370 if (SumMass > TotalCollisionMass) return FAL 371 SumMass = SumMom.m2(); 372 if (SumMass < 0) return FALSE; 373 SumMass = std::sqrt(SumMass); 374 375 // Compute c.m.s. hadron velocity and boost 376 G4ThreeVector Beta = -SumMom.boostVector(); 377 //G4cout << " == pre boost 2 "<< SumMo 378 //--old Output->Boost(Beta); 379 for(i = 0; i < Output->size(); i++) 380 { 381 G4LorentzVector mom = G4LorentzVector((*Ou 382 mom *= Beta; 383 (*Output)[i]->SetMomentum(mom.vect()); 384 (*Output)[i]->SetTotalEnergy(mom.e()); 385 } 386 387 // Scale total c.m.s. hadron energy (hadron 388 // It should be equal interaction mass 389 G4double Scale = 0,OldScale=0; 390 G4double factor = 1.; 391 G4int cAttempt = 0; 392 G4double Sum = 0; 393 G4bool success = false; 394 for(cAttempt = 0; cAttempt < nAttemptScale; 395 { 396 Sum = 0; 397 for(i = 0; i < Output->size(); i++) 398 { 399 G4LorentzVector HadronMom = G4LorentzVec 400 HadronMom.setVect(HadronMom.vect()+ fact 401 G4double E = std::sqrt(HadronMom.vect(). 402 HadronMom.setE(E); 403 (*Output)[i]->SetMomentum(HadronMom.vect 404 (*Output)[i]->SetTotalEnergy(HadronMom.e 405 Sum += E; 406 } 407 OldScale=Scale; 408 Scale = TotalCollisionMass/Sum - 1; 409 // G4cout << "E/P corr - " << cAttempt << 410 if (std::abs(Scale) <= ErrLimit 411 || OldScale == Scale) // protect ' 412 { 413 if (debug_G4BinaryLightIonReactionResult 414 success = true; 415 break; 416 } 417 if ( cAttempt > 10 ) 418 { 419 // G4cout << " speed it up? " << 420 factor=std::max(1.,G4Log(std::abs(OldSca 421 // G4cout << " ? factor ? " << factor 422 } 423 } 424 425 if( (!success) && debug_G4BinaryLightIonRea 426 { 427 G4cout << "G4G4BinaryLightIonReaction::Ene 428 G4cout << " Scale not unity at end of it 429 G4cout << " Increase number of attempts 430 } 431 432 // Compute c.m.s. interaction velocity and K 433 Beta = TotalCollisionMom.boostVector(); 434 //--old Output->Boost(Beta); 435 for(i = 0; i < Output->size(); i++) 436 { 437 G4LorentzVector mom = G4LorentzVector((*Ou 438 mom *= Beta; 439 (*Output)[i]->SetMomentum(mom.vect()); 440 (*Output)[i]->SetTotalEnergy(mom.e()); 441 } 442 return TRUE; 443 } 444 G4bool G4BinaryLightIonReaction::SetLighterAsP 445 { 446 G4bool swapped = false; 447 if(tA<pA) 448 { 449 swapped = true; 450 G4int tmp(0); 451 tmp = tA; tA=pA; pA=tmp; 452 tmp = tZ; tZ=pZ; pZ=tmp; 453 G4double m1=G4ParticleTable::GetParticle 454 G4LorentzVector it(m1, G4ThreeVector(0,0 455 mom = toBreit*it; 456 } 457 return swapped; 458 } 459 G4ReactionProductVector * G4BinaryLightIonReac 460 { 461 // Check if kinematically nuclei can fuse. 462 G4double mFused=G4ParticleTable::GetParticl 463 G4double mTarget=G4ParticleTable::GetPartic 464 G4LorentzVector pCompound(mom.e()+mTarget,m 465 G4double m2Compound=pCompound.m2(); 466 if (m2Compound < sqr(mFused) ) { 467 //G4cout << "G4BLIC: projectile p, mTarge 468 // << " " << sqrt(m2Compound)<< " " 469 return 0; 470 } 471 472 G4Fragment aPreFrag; 473 aPreFrag.SetZandA_asInt(pZ+tZ, pA+tA); 474 aPreFrag.SetNumberOfParticles(pA); 475 aPreFrag.SetNumberOfCharged(pZ); 476 aPreFrag.SetNumberOfHoles(0); 477 //GF FIXME: whyusing plop in z direction? t 478 //G4ThreeVector plop(0.,0., mom.vect().mag( 479 //G4LorentzVector aL(mom.t()+mTarget, plop) 480 G4LorentzVector aL(mom.t()+mTarget,mom.vect 481 aPreFrag.SetMomentum(aL); 482 483 484 //G4cout << "Fragment INFO "<< pA+tA 485 // << aL <<" "<<G4endl << aPreF 486 G4ReactionProductVector * cascaders = thePr 487 //G4double tSum = 0; 488 for(size_t count = 0; count<cascaders->size 489 { 490 cascaders->operator[](count)->SetNewlyAd 491 //tSum += cascaders->operator[](count)-> 492 } 493 // G4cout << "Exiting pre-compound on 494 return cascaders; 495 } 496 G4ReactionProductVector * G4BinaryLightIonReac 497 { 498 G4ReactionProductVector * result = 0; 499 G4double projectileMass(0); 500 G4LorentzVector it; 501 502 G4int tryCount(0); 503 do 504 { 505 ++tryCount; 506 projectile3dNucleus = new G4Fancy3DNu 507 projectile3dNucleus->Init(pA, pZ); 508 projectile3dNucleus->CenterNucleons() 509 projectileMass=G4ParticleTable::GetPa 510 projectile3dNucleus->GetCharge( 511 it=toBreit * G4LorentzVector(projecti 512 513 target3dNucleus = new G4Fancy3DNucleu 514 target3dNucleus->Init(tA, tZ); 515 G4double impactMax = target3dNucleus- 516 // G4cout << "out radius - nuc 517 G4double aX=(2.*G4UniformRand()-1.)*i 518 G4double aY=(2.*G4UniformRand()-1.)*i 519 G4ThreeVector pos(aX, aY, -2.*impactM 520 521 G4KineticTrackVector * initalState = 522 projectile3dNucleus->StartLoop(); 523 G4Nucleon * aNuc; 524 G4LorentzVector tmpV(0,0,0,0); 525 #ifdef debug_BLIR_finalstate 526 G4LorentzVector pinitial; 527 #endif 528 G4LorentzVector nucleonMom(1./pA*mom) 529 nucleonMom.setZ(nucleonMom.vect().mag 530 nucleonMom.setX(0); 531 nucleonMom.setY(0); 532 theFermi.Init(pA,pZ); 533 while( (aNuc=projectile3dNucleus->Get 534 { 535 G4LorentzVector p4 = aNuc->GetMome 536 tmpV+=p4; 537 G4ThreeVector nucleonPosition(aNuc 538 G4double density=(projectile3dNucl 539 nucleonPosition += pos; 540 G4KineticTrack * it1 = new G4Kinet 541 it1->SetState(G4KineticTrack::outs 542 G4double pfermi= theFermi.GetFermi 543 G4double mass = aNuc->GetDefinitio 544 G4double Efermi= std::sqrt( sqr(ma 545 it1->SetProjectilePotential(-Eferm 546 initalState->push_back(it1); 547 #ifdef debug_BLIR_finalstate 548 pinitial += it1->Get4Momentum() 549 #endif 550 } 551 552 result=theModel->Propagate(initalStat 553 #ifdef debug_BLIR_finalstate 554 if( result && result->size()>0) 555 { 556 G4cout << " Cascade result " << G4endl; 557 G4LorentzVector presult; 558 G4ReactionProductVector::iterator 559 G4ReactionProduct xp; 560 for (iter=result->begin(); iter ! 561 { 562 presult += G4LorentzVector((*ite 563 G4cout << (*iter)->GetDefinition()->Ge 564 << "("<< (*iter)->GetMomentum().x()<<" 565 << (*iter)->GetMomentum().y()<<"," 566 << (*iter)->GetMomentum().z()<<";" 567 << (*iter)->GetTotalEnergy() <<")"< 568 } 569 570 G4cout << "BLIC check result : in 571 << " final " << presult 572 << " IF - FF " << pinitial +G 573 574 } 575 #endif 576 if( result && result->size()==0) 577 { 578 delete result; 579 result=0; 580 } 581 if ( ! result ) 582 { 583 delete target3dNucleus; 584 delete projectile3dNucleus; 585 } 586 587 // std::for_each(initalState->begin() 588 // delete initalState; 589 590 } while (! result && tryCount< 150); / 591 return result; 592 } 593 G4double G4BinaryLightIonReaction::GetProjecti 594 { 595 596 G4Nucleon * aNuc; 597 // the projectileNucleus excitation ener 598 G4double theStatisticalExEnergy = 0; 599 projectile3dNucleus->StartLoop(); 600 while( (aNuc=projectile3dNucleus->GetNex 601 { 602 //G4cout << " Nucleon : " << a 603 if(aNuc->AreYouHit()) { 604 G4ThreeVector aPosition(aNuc->GetP 605 G4double localDensity = projectile 606 G4double localPfermi = theFermi.Ge 607 G4double nucMass = aNuc->GetDefini 608 G4double localFermiEnergy = std::s 609 G4double deltaE = localFermiEnergy 610 theStatisticalExEnergy += deltaE; 611 } 612 } 613 return theStatisticalExEnergy; 614 } 615 616 G4LorentzVector G4BinaryLightIonReaction::Sort 617 { 618 unsigned int i(0); 619 spectatorA=spectatorZ=0; 620 G4LorentzVector pspectators(0,0,0,0); 621 pFinalState=G4LorentzVector(0,0,0,0); 622 for(i=0; i<result->size(); i++) 623 { 624 if( (*result)[i]->GetNewlyAdded() ) 625 { 626 pFinalState += G4LorentzVector( (*res 627 cascaders->push_back((*result)[i]); 628 } 629 else { 630 // G4cout <<" spectator ... 631 pspectators += G4LorentzVector( (*res 632 spectators->push_back((*result)[i]); 633 spectatorA++; 634 spectatorZ+= G4lrint((*result)[i]->Ge 635 } 636 637 // G4cout << (*result)[i]<< " " 638 // << (*result)[i]->GetDefinition 639 // << (*result)[i]->GetMomentum() 640 // << (*result)[i]->GetTotalEnerg 641 } 642 //G4cout << "pFinalState / pspectators, 643 // << " (" << spectatorA << ", "<< sp 644 645 return pspectators; 646 } 647 648 void G4BinaryLightIonReaction::DeExciteSpectat 649 650 { 651 // call precompound model 652 G4ReactionProductVector * proFrag = 0; 653 G4LorentzVector pFragment(0.,0.,0.,0.); 654 // G4cout << " == pre boost 1 "<< mome 655 G4LorentzRotation boost_fragments; 656 // G4cout << " == post boost 1 "<< mom 657 // G4LorentzRotation boost_spectator_mom 658 // G4cout << "- momentum " << boost_spe 659 G4LorentzVector pFragments(0,0,0,0); 660 661 if(spectatorZ>0 && spectatorA>1) 662 { 663 // Make the fragment 664 G4Fragment aProRes; 665 aProRes.SetZandA_asInt(spectatorZ, spect 666 aProRes.SetNumberOfParticles(0); 667 aProRes.SetNumberOfCharged(0); 668 aProRes.SetNumberOfHoles(pA-spectatorA); 669 G4double mFragment=G4ParticleTable::GetP 670 pFragment=G4LorentzVector(0,0,0,mFragmen 671 aProRes.SetMomentum(pFragment); 672 673 proFrag = theHandler->BreakItUp(aProRes) 674 675 boost_fragments = G4LorentzRotation(pSpe 676 677 // G4cout << " Fragment a,z, Mass Fr 678 // << spectatorA <<" "<< spectator 679 // << momentum.mag() <<" "<< momen 680 // << " "<<theStatisticalExEnergy 681 // << " "<< boost_fragments*pFragm 682 G4ReactionProductVector::iterator ispect 683 for (ispectator=spectators->begin();ispe 684 { 685 delete *ispectator; 686 } 687 } 688 else if(spectatorA!=0) 689 { 690 G4ReactionProductVector::iterator ispecta 691 for (ispectator=spectators->begin();ispec 692 { 693 (*ispectator)->SetNewlyAdded(true); 694 cascaders->push_back(*ispectator); 695 pFinalState+=G4LorentzVector((*ispect 696 //G4cout << "BLIC: spectator 697 // << (*ispectator)->GetDefi 698 // << (*ispectator)->GetMome 699 // << (*ispectator)->GetTota 700 } 701 702 } 703 // / if (spectators) 704 delete spectators; 705 706 // collect the evaporation part and boost t 707 G4ReactionProductVector::iterator ii; 708 if(proFrag) 709 { 710 for(ii=proFrag->begin(); ii!=proFrag->en 711 { 712 (*ii)->SetNewlyAdded(true); 713 G4LorentzVector tmp((*ii)->GetMomentu 714 tmp *= boost_fragments; 715 (*ii)->SetMomentum(tmp.vect()); 716 (*ii)->SetTotalEnergy(tmp.e()); 717 // result->push_back(*ii); 718 pFragments += tmp; 719 } 720 } 721 722 // G4cout << "Fragmented p, momentum, de 723 // <<" "<< pFragments-momentum < 724 725 // correct p/E of Cascade secondaries 726 G4LorentzVector pCas=pInitialState - pFragm 727 728 //G4cout <<"BLIC: Going to correct from 729 // the creation of excited fragment did vi 730 G4bool EnergyIsCorrect=EnergyAndMomentumCor 731 if ( ! EnergyIsCorrect && debug_G4BinaryLig 732 { 733 G4cout << "G4BinaryLightIonReaction E/P 734 } 735 736 // Add deexcitation secondaries 737 if(proFrag) 738 { 739 for(ii=proFrag->begin(); ii!=proFrag->en 740 { 741 cascaders->push_back(*ii); 742 } 743 delete proFrag; 744 } 745 //G4cout << "EnergyIsCorrect? " << Energ 746 if ( ! EnergyIsCorrect ) 747 { 748 // G4cout <<" ! EnergyIsCorrect " << 749 if (! EnergyAndMomentumCorrector(cascade 750 { 751 if(debug_G4BinaryLightIonReactionResu 752 G4cout << "G4BinaryLightIonReactio 753 } 754 } 755 756 } 757 758