00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "G4QMDCollision.hh"
00034 #include "G4Scatterer.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "Randomize.hh"
00038
00039 G4QMDCollision::G4QMDCollision()
00040 : deltar ( 4 )
00041 , bcmax0 ( 1.323142 )
00042 , bcmax1 ( 2.523 )
00043
00044
00045
00046 , epse ( 0.0001 )
00047 {
00048 theScatterer = new G4Scatterer();
00049 }
00050
00051
00052
00053 G4QMDCollision::~G4QMDCollision()
00054 {
00055 delete theScatterer;
00056 }
00057
00058
00059 void G4QMDCollision::CalKinematicsOfBinaryCollisions( G4double dt )
00060 {
00061 G4double deltaT = dt;
00062
00063 G4int n = theSystem->GetTotalNumberOfParticipant();
00064
00065
00066 for ( G4int i = 0 ; i < n ; i++ )
00067 {
00068 theSystem->GetParticipant( i )->UnsetHitMark();
00069 theSystem->GetParticipant( i )->UnsetHitMark();
00070
00071 }
00072
00073
00074
00075
00076 for ( G4int i = 0 ; i < n ; i++ )
00077 {
00078
00079
00080
00081 if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() )
00082 {
00083
00084 G4bool decayed = false;
00085
00086 G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition();
00087 G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum();
00088 G4ThreeVector r0 = theSystem->GetParticipant( i )->GetPosition();
00089
00090 G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum();
00091
00092 G4double eini = theMeanField->GetTotalPotential() + p40.e();
00093
00094 G4int n0 = theSystem->GetTotalNumberOfParticipant();
00095 G4int i0 = 0;
00096
00097 G4bool isThisEnergyOK = false;
00098
00099 for ( G4int ii = 0 ; ii < 4 ; ii++ )
00100 {
00101
00102
00103 G4LorentzVector p400 = p40;
00104
00105 p400 *= GeV;
00106
00107 G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
00108
00109 G4KineticTrackVector* secs = NULL;
00110 secs = kt.Decay();
00111 G4int id = 0;
00112 G4double et = 0;
00113 if ( secs )
00114 {
00115 for ( G4KineticTrackVector::iterator it
00116 = secs->begin() ; it != secs->end() ; it++ )
00117 {
00118
00119
00120
00121
00122
00123
00124
00125 if ( id == 0 )
00126 {
00127 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
00128 theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
00129 theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
00130
00131 et += (*it)->Get4Momentum().e()/GeV;
00132 }
00133 if ( id > 0 )
00134 {
00135
00136 theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
00137 et += (*it)->Get4Momentum().e()/GeV;
00138 if ( id > 1 )
00139 {
00140
00141
00142 }
00143 }
00144 id++;
00145
00146 delete *it;
00147 }
00148
00149 theMeanField->Update();
00150 i0 = id-1;
00151
00152 delete secs;
00153 }
00154
00155
00156
00157 G4double efin = theMeanField->GetTotalPotential() + et;
00158
00159
00160
00161 if ( std::abs ( eini - efin ) < epse*10 )
00162 {
00163
00164 isThisEnergyOK = true;
00165 break;
00166 }
00167 else
00168 {
00169
00170 theSystem->GetParticipant( i )->SetDefinition( pd0 );
00171 theSystem->GetParticipant( i )->SetPosition( r0 );
00172 theSystem->GetParticipant( i )->SetMomentum( p0 );
00173
00174 for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
00175 {
00176
00177
00178 theSystem->DeleteParticipant( i0i+n0 );
00179 }
00180
00181 theMeanField->Update();
00182 }
00183
00184 }
00185
00186
00187
00188 if ( isThisEnergyOK == true )
00189 {
00190 if ( theMeanField->IsPauliBlocked ( i ) != true )
00191 {
00192
00193 G4bool allOK = true;
00194 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
00195 {
00196 if ( theMeanField->IsPauliBlocked ( i0i+n0 ) == true )
00197 {
00198 allOK = false;
00199 break;
00200 }
00201 }
00202
00203 if ( allOK )
00204 {
00205 decayed = true;
00206 }
00207 }
00208
00209 }
00210
00211
00212 if ( decayed )
00213 {
00214
00215
00216 theSystem->GetParticipant( i )->SetHitMark();
00217 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
00218 {
00219 theSystem->GetParticipant( i0i+n0 )->SetHitMark();
00220 }
00221
00222 }
00223 else
00224 {
00225
00226
00227
00228 if ( isThisEnergyOK == true )
00229 {
00230
00231 theSystem->GetParticipant( i )->SetDefinition( pd0 );
00232 theSystem->GetParticipant( i )->SetPosition( r0 );
00233 theSystem->GetParticipant( i )->SetMomentum( p0 );
00234
00235 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
00236 {
00237
00238
00239 theSystem->DeleteParticipant( i0i+n0 );
00240 }
00241
00242 theMeanField->Update();
00243 }
00244
00245 }
00246
00247 }
00248 }
00249
00250
00251
00252 n = theSystem->GetTotalNumberOfParticipant();
00253
00254
00255
00256 for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ )
00257 {
00258
00259
00260 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
00261
00262 G4ThreeVector ri = theSystem->GetParticipant( i )->GetPosition();
00263 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
00264 G4double rmi = theSystem->GetParticipant( i )->GetMass();
00265 G4ParticleDefinition* pdi = theSystem->GetParticipant( i )->GetDefinition();
00266
00267 if ( pdi->GetPDGMass() == 0.0 ) continue;
00268
00269
00270 for ( G4int j = 0 ; j < i ; j++ )
00271 {
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
00284 if ( theSystem->GetParticipant( j )->IsThisHit() ) continue;
00285
00286
00287
00288
00289 if ( theSystem->GetParticipant( i )->IsThisProjectile() )
00290 {
00291 if ( theSystem->GetParticipant( j )->IsThisProjectile() ) continue;
00292 }
00293 else if ( theSystem->GetParticipant( i )->IsThisTarget() )
00294 {
00295 if ( theSystem->GetParticipant( j )->IsThisTarget() ) continue;
00296 }
00297
00298
00299 G4ThreeVector rj = theSystem->GetParticipant( j )->GetPosition();
00300 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
00301 G4double rmj = theSystem->GetParticipant( j )->GetMass();
00302 G4ParticleDefinition* pdj = theSystem->GetParticipant( j )->GetDefinition();
00303
00304 if ( pdj->GetPDGMass() == 0.0 ) continue;
00305
00306 G4double rr2 = theMeanField->GetRR2( i , j );
00307
00308
00309 if ( rr2 > deltar*deltar ) continue;
00310
00311
00312
00313
00314 G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) );
00315
00316 G4double cutoff = 0.0;
00317 G4double bcmax = 0.0;
00318
00319
00320
00321 if ( rmi < 0.94 && rmj < 0.94 )
00322 {
00323
00324 cutoff = rmi + rmj + 0.02;
00325 bcmax = bcmax0;
00326
00327
00328 }
00329 else
00330 {
00331 cutoff = rmi + rmj;
00332 bcmax = bcmax1;
00333
00334
00335 }
00336
00337
00338 if ( srt < cutoff ) continue;
00339
00340 G4ThreeVector dr = ri - rj;
00341 G4double rsq = dr*dr;
00342
00343 G4double pij = p4i*p4j;
00344 G4double pidr = p4i.vect()*dr;
00345 G4double pjdr = p4j.vect()*dr;
00346
00347 G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij );
00348 G4double bij = pidr / rmi - pjdr*rmi/pij;
00349 G4double cij = rsq + ( pidr / rmi ) * ( pidr / rmi );
00350 G4double brel = std::sqrt ( std::abs ( cij - bij*bij/aij ) );
00351
00352 if ( brel > bcmax ) continue;
00353
00354
00355 G4double bji = -pjdr/rmj + pidr * rmj /pij;
00356
00357 G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi;
00358 G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375 if ( std::abs ( ti + tj ) > deltaT ) continue;
00376
00377
00378 G4ThreeVector beta = ( p4i + p4j ).boostVector();
00379
00380 G4LorentzVector p = p4i;
00381 G4LorentzVector p4icm = p.boost( p.findBoostToCM ( p4j ) );
00382 G4ThreeVector pcm = p4icm.vect();
00383
00384 G4double prcm = pcm.mag();
00385
00386 if ( prcm <= 0.00001 ) continue;
00387
00388
00389 G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollision ( i , j ) );
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417 if ( energetically_forbidden == true )
00418 {
00419
00420
00421
00422
00423
00424 theSystem->GetParticipant( i )->SetMomentum( p4i.vect() );
00425 theSystem->GetParticipant( i )->SetDefinition( pdi );
00426 theSystem->GetParticipant( i )->SetPosition( ri );
00427
00428 theSystem->GetParticipant( j )->SetMomentum( p4j.vect() );
00429 theSystem->GetParticipant( j )->SetDefinition( pdj );
00430 theSystem->GetParticipant( j )->SetPosition( rj );
00431
00432 theMeanField->Cal2BodyQuantities( i );
00433 theMeanField->Cal2BodyQuantities( j );
00434
00435 }
00436 else
00437 {
00438
00439
00440 G4bool absorption = false;
00441 if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true;
00442 if ( absorption )
00443 {
00444
00445 i = i-1;
00446 n = n-1;
00447 }
00448
00449
00450
00451
00452 theSystem->GetParticipant( i )->UnsetInitialMark();
00453 if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark();
00454
00455 theSystem->GetParticipant( i )->SetHitMark();
00456 if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark();
00457
00458 theSystem->IncrementCollisionCounter();
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486 }
00487
00488 }
00489
00490 }
00491
00492
00493 }
00494
00495
00496
00497 G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollision( G4int i , G4int j )
00498 {
00499
00500
00501
00502
00503 G4bool result = false;
00504 G4bool energyOK = false;
00505 G4bool pauliOK = false;
00506 G4bool abs = false;
00507 G4QMDParticipant* absorbed = NULL;
00508
00509 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
00510 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
00511
00512
00513
00514 G4double epot = theMeanField->GetTotalPotential();
00515
00516 G4double eini = epot + p4i.e() + p4j.e();
00517
00518
00519
00520 G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition();
00521 G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition();
00522 G4LorentzVector p4i0 = p4i*GeV;
00523 G4LorentzVector p4j0 = p4j*GeV;
00524 G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi;
00525 G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi;
00526
00527 for ( G4int iitry = 0 ; iitry < 4 ; iitry++ )
00528 {
00529
00530 abs = false;
00531
00532 G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 );
00533 G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 );
00534
00535 G4LorentzVector p4ix_new;
00536 G4LorentzVector p4jx_new;
00537 G4KineticTrackVector* secs = NULL;
00538 secs = theScatterer->Scatter( kt1 , kt2 );
00539
00540
00541
00542
00543
00544
00545 if ( secs )
00546 {
00547 G4int iti = 0;
00548 if ( secs->size() == 2 )
00549 {
00550 for ( G4KineticTrackVector::iterator it
00551 = secs->begin() ; it != secs->end() ; it++ )
00552 {
00553 if ( iti == 0 )
00554 {
00555 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
00556 p4ix_new = (*it)->Get4Momentum()/GeV;
00557
00558 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
00559 }
00560 if ( iti == 1 )
00561 {
00562 theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() );
00563 p4jx_new = (*it)->Get4Momentum()/GeV;
00564
00565 theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() );
00566 }
00567
00568 iti++;
00569 }
00570 }
00571 else if ( secs->size() == 1 )
00572 {
00573
00574 abs = true;
00575
00576
00577 theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() );
00578 p4ix_new = secs->front()->Get4Momentum()/GeV;
00579 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
00580
00581 }
00582
00583
00584 if ( secs->size() > 2 )
00585 {
00586
00587 G4cout << "G4QMDCollision secs size > 2; " << secs->size() << G4endl;
00588
00589 for ( G4KineticTrackVector::iterator it
00590 = secs->begin() ; it != secs->end() ; it++ )
00591 {
00592 G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl;
00593 }
00594
00595 }
00596
00597
00598 for ( G4KineticTrackVector::iterator it
00599 = secs->begin() ; it != secs->end() ; it++ )
00600 {
00601 delete *it;
00602 }
00603
00604 delete secs;
00605 }
00606
00607
00608 if ( !abs )
00609 {
00610 theMeanField->Cal2BodyQuantities( i );
00611 theMeanField->Cal2BodyQuantities( j );
00612 }
00613 else
00614 {
00615 absorbed = theSystem->EraseParticipant( j );
00616 theMeanField->Update();
00617 }
00618
00619 epot = theMeanField->GetTotalPotential();
00620
00621 G4double efin = epot + p4ix_new.e() + p4jx_new.e();
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632 if ( std::abs ( eini - efin ) < epse )
00633 {
00634
00635
00636
00637
00638
00639
00640
00641 energyOK = true;
00642 break;
00643 }
00644 else
00645 {
00646
00647 if ( abs )
00648 {
00649
00650 theSystem->InsertParticipant( absorbed , j );
00651 theMeanField->Update();
00652 }
00653
00654 }
00655
00656 }
00657
00658
00659
00660 if ( energyOK )
00661 {
00662
00663
00664 if ( !abs )
00665 {
00666 if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) )
00667 {
00668
00669 pauliOK = true;
00670 }
00671 }
00672 else
00673 {
00674
00675
00676 if ( theMeanField->IsPauliBlocked ( i-1 ) == false )
00677 {
00678
00679 delete absorbed;
00680 pauliOK = true;
00681 }
00682 }
00683
00684
00685 if ( pauliOK )
00686 {
00687 result = true;
00688 }
00689 else
00690 {
00691
00692 if ( abs )
00693 {
00694
00695 theSystem->InsertParticipant( absorbed , j );
00696 theMeanField->Update();
00697 }
00698 }
00699 }
00700
00701 return result;
00702
00703 }
00704
00705
00706
00707 G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD( G4double sig , G4double cutoff , G4ThreeVector pcm , G4double prcm , G4double srt , G4ThreeVector beta , G4double gamma , G4int i , G4int j )
00708 {
00709
00710
00711
00712 G4bool result = true;
00713
00714 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
00715 G4double rmi = theSystem->GetParticipant( i )->GetMass();
00716 G4int zi = theSystem->GetParticipant( i )->GetChargeInUnitOfEplus();
00717
00718 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
00719 G4double rmj = theSystem->GetParticipant( j )->GetMass();
00720 G4int zj = theSystem->GetParticipant( j )->GetChargeInUnitOfEplus();
00721
00722 G4double pr = prcm;
00723
00724 G4double c2 = pcm.z()/pr;
00725
00726 G4double csrt = srt - cutoff;
00727
00728
00729
00730
00731 G4double asrt = srt - rmi - rmj;
00732 G4double pra = prcm;
00733
00734
00735
00736 G4double elastic = 0.0;
00737
00738 if ( zi == zj )
00739 {
00740 if ( csrt < 0.4286 )
00741 {
00742 elastic = 35.0 / ( 1. + csrt * 100.0 ) + 20.0;
00743 }
00744 else
00745 {
00746 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
00747 * 2. / pi + 1.0 ) * 9.65 + 7.0;
00748 }
00749 }
00750 else
00751 {
00752 if ( csrt < 0.4286 )
00753 {
00754 elastic = 28.0 / ( 1. + csrt * 100.0 ) + 27.0;
00755 }
00756 else
00757 {
00758 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
00759 * 2. / pi + 1.0 ) * 12.34 + 10.0;
00760 }
00761 }
00762
00763
00764
00765
00766
00767
00768 if ( G4UniformRand() > elastic / sig )
00769 {
00770
00771
00772 return result;
00773 }
00774 else
00775 {
00776
00777 }
00778
00779
00780
00781 G4double as = std::pow ( 3.65 * asrt , 6 );
00782 G4double a = 6.0 * as / (1.0 + as);
00783 G4double ta = -2.0 * pra*pra;
00784 G4double x = G4UniformRand();
00785 G4double t1 = std::log( (1-x) * std::exp(2.*a*ta) + x ) / a;
00786 G4double c1 = 1.0 - t1/ta;
00787
00788 if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0;
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798 t1 = 2.0*pi*G4UniformRand();
00799
00800 G4double t2 = 0.0;
00801 if ( pcm.x() == 0.0 && pcm.y() == 0 )
00802 {
00803 t2 = 0.0;
00804 }
00805 else
00806 {
00807 t2 = std::atan2( pcm.y() , pcm.x() );
00808 }
00809
00810
00811 G4double s1 = std::sqrt ( 1.0 - c1*c1 );
00812 G4double s2 = std::sqrt ( 1.0 - c2*c2 );
00813
00814 G4double ct1 = std::cos(t1);
00815 G4double st1 = std::sin(t1);
00816
00817 G4double ct2 = std::cos(t2);
00818 G4double st2 = std::sin(t2);
00819
00820 G4double ss = c2*s1*ct1 + s2*c1;
00821
00822 pcm.setX( pr * ( ss*ct2 - s1*st1*st2) );
00823 pcm.setY( pr * ( ss*st2 + s1*st1*ct2) );
00824 pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) );
00825
00826
00827
00828 G4double epot = theMeanField->GetTotalPotential();
00829
00830 G4double eini = epot + p4i.e() + p4j.e();
00831 G4double etwo = p4i.e() + p4j.e();
00832
00833
00834
00835
00836
00837
00838
00839
00840 for ( G4int itry = 0 ; itry < 4 ; itry++ )
00841 {
00842
00843 G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm );
00844 G4double pibeta = pcm*beta;
00845
00846 G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm );
00847
00848 G4ThreeVector pi_new = beta*trans + pcm;
00849
00850 G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm );
00851 trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm );
00852
00853 G4ThreeVector pj_new = beta*trans - pcm;
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863 theSystem->GetParticipant( i )->SetMomentum( pi_new );
00864 theSystem->GetParticipant( j )->SetMomentum( pj_new );
00865
00866 G4double pi_new_e = (theSystem->GetParticipant( i )->Get4Momentum()).e();
00867 G4double pj_new_e = (theSystem->GetParticipant( j )->Get4Momentum()).e();
00868
00869 theMeanField->Cal2BodyQuantities( i );
00870 theMeanField->Cal2BodyQuantities( j );
00871
00872 epot = theMeanField->GetTotalPotential();
00873
00874 G4double efin = epot + pi_new_e + pj_new_e ;
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884 if ( std::abs ( eini - efin ) < epse )
00885 {
00886
00887
00888
00889
00890
00891
00892
00893 }
00894
00895
00896 if ( std::abs ( eini - efin ) < epse ) return result;
00897
00898 G4double cona = ( eini - efin + etwo ) / gamma;
00899 G4double fac2 = 1.0 / ( 4.0 * cona*cona * pr*pr ) *
00900 ( ( cona*cona - ( rmi*rmi + rmj*rmj ) )*( cona*cona - ( rmi*rmi + rmj*rmj ) )
00901 - 4.0 * rmi*rmi * rmj*rmj );
00902
00903 if ( fac2 > 0 )
00904 {
00905 G4double fact = std::sqrt ( fac2 );
00906 pcm = fact*pcm;
00907 }
00908
00909
00910 }
00911
00912
00913 result = false;
00914
00915 return result;
00916
00917 }