| 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 // 081024 G4NucleiPropertiesTable:: to G4Nucle
27 //
28 // 230308 "samplingPosition_cluster" function
29 // 230308 Parameter "radam" shold be independe
30 // 230308 Sampling momentum of the particle do
31 // 230308 Binding energy evaluated by "GetTota
32
33 #include "G4LightIonQMDGroundStateNucleus.hh"
34
35 #include "G4NucleiProperties.hh"
36 #include "G4Proton.hh"
37 #include "G4Neutron.hh"
38 #include "G4Exp.hh"
39 #include "G4Pow.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "Randomize.hh"
42
43 G4LightIonQMDGroundStateNucleus::G4LightIonQMD
44 : maxTrial ( 1000 )
45 , r00 ( 1.124 ) // radius parameter for Woods
46 , r01 ( 0.5 ) // radius parameter for Woods
47 , saa ( 0.2 ) // diffuse parameter for init
48 , rada ( 0.9 ) // cutoff parameter
49 , radb ( 0.3 ) // cutoff parameter
50 , dsam ( 1.5 ) // minimum distance for same
51 , ddif ( 1.0 ) // minimum distance for diffe
52 , edepth ( 0.0 )
53 , epse ( 0.000001 ) // torelance for energy i
54 , meanfield ( NULL )
55 {
56
57 //std::cout << " G4LightIonQMDGroundStateNu
58
59 dsam2 = dsam*dsam;
60 ddif2 = ddif*ddif;
61
62 G4LightIonQMDParameters* parameters = G4Lig
63
64 hbc = parameters->Get_hbc();
65 gamm = parameters->Get_gamm();
66 cpw = parameters->Get_cpw();
67 cph = parameters->Get_cph();
68 epsx = parameters->Get_epsx();
69 cpc = parameters->Get_cpc();
70
71 cdp = parameters->Get_cdp();
72 c0p = parameters->Get_c0p();
73 c3p = parameters->Get_c3p();
74 csp = parameters->Get_csp();
75 clp = parameters->Get_clp();
76
77 ebini = 0.0;
78 // Following 10 lines should be here, right
79 // Otherwise, mass number cannot be conserv
80 // the target are nucleons.
81 //Nucleon primary or target case;
82 if ( z == 1 && a == 1 ) { // Hydrogen Cas
83 SetParticipant( new G4QMDParticipant( G4
84 // ebini = 0.0;
85 return;
86 }
87 else if ( z == 0 && a == 1 ) { // Neutron p
88 SetParticipant( new G4QMDParticipant( G4
89 // ebini = 0.0;
90 return;
91 }
92
93
94 //edepth = 0.0;
95
96 for ( int i = 0 ; i < a ; i++ )
97 {
98
99 G4ParticleDefinition* pd;
100
101 if ( i < z )
102 {
103 pd = G4Proton::Proton();
104 }
105 else
106 {
107 pd = G4Neutron::Neutron();
108 }
109
110 G4ThreeVector p( 0.0 );
111 G4ThreeVector r( 0.0 );
112 G4QMDParticipant* aParticipant = new G4Q
113 SetParticipant( aParticipant );
114
115 }
116
117 G4double radious = r00 * G4Pow::GetInstance
118
119 rt00 = radious - r01;
120 radm = radious; // JQMD, Skyrme, RMF, Y-H.
121 rmax = 1.0 / ( 1.0 + G4Exp ( -rt00/saa ) );
122
123 //maxTrial = 1000;
124
125
126 meanfield = new G4LightIonQMDMeanField();
127 meanfield->SetSystem( this );
128
129 //std::cout << "G4LightIonQMDGroundStateNuc
130 packNucleons();
131 //std::cout << "G4LightIonQMDGroundStateNuc
132
133 delete meanfield;
134
135 }
136
137
138
139 void G4LightIonQMDGroundStateNucleus::packNucl
140 {
141
142 //std::cout << "G4LightIonQMDGroundStateNuc
143
144 ebini = - G4NucleiProperties::GetBindingEne
145
146 G4double ebin00 = ebini * 0.001;
147
148 G4double ebin0 = 0.0;
149 G4double ebin1 = 0.0;
150
151 if ( GetMassNumber() != 4 )
152 {
153 ebin0 = ( ebini - 0.5 ) * 0.001;
154 ebin1 = ( ebini + 0.5 ) * 0.001;
155 }
156 else
157 {
158 ebin0 = ( ebini - 1.5 ) * 0.001;
159 ebin1 = ( ebini + 1.5 ) * 0.001;
160 }
161
162 G4int n0Try = 0;
163 G4bool isThisOK = false;
164 G4double energy_QMD = 0.0; // Y-H.S. and A.
165 G4int Nucle_Z = GetAtomicNumber(); // Y-H.S
166 G4int Nucle_A = GetMassNumber(); // Y-H.S.
167
168 while ( n0Try < maxTrial ) // Loop checking
169 {
170 n0Try++;
171 //std::cout << "TKDB packNucleons n0Try
172
173 // Sampling Position
174
175 //std::cout << "TKDB Sampling Position "
176
177 G4bool areThesePsOK = false;
178 G4int npTry = 0;
179 while ( npTry < maxTrial ) // Loop chec
180 {
181 //std::cout << "TKDB Sampling Position
182 if( (Nucle_Z == 6 && Nucle_A == 12)
183 {
184 //std::cout << "alpha cluster "
185 //npTry++;
186 //G4int i = 0;
187 if ( samplingPosition_cluster(
188 {
189 //G4double rsquare = meanfi
190 //G4cout << "R-square " <<
191 areThesePsOK = true;
192 break;
193 //exit(0);
194 /*
195 //std::cout << "packNucleo
196 for ( i = 1 ; i < GetMassN
197 {
198 //std::cout << "packNucleo
199 if ( !( samplingPosition_c
200 {
201 //std::cout << "packNucleo
202 break;
203 }
204 }
205 if ( i == GetMassNumber()
206 {
207 //std::cout << "packNucleo
208 areThesePsOK = true;
209 break;
210 }
211 */
212 }
213 } // Alpha cluster model added by Y
214 else
215 {
216 npTry++;
217 G4int i = 0;
218 if ( samplingPosition( i ) )
219 {
220 //std::cout << "packNucleon
221 for ( i = 1 ; i < GetMassNu
222 {
223 //std::cout << "packNuc
224 if ( !( samplingPositio
225 {
226 //std::cout << "pac
227 break;
228 }
229 }
230 if ( i == GetMassNumber() )
231 {
232 //std::cout << "packNuc
233 areThesePsOK = true;
234 break;
235 }
236 }
237 }
238 }
239 if ( areThesePsOK == false ) continue;
240
241 //std::cout << "TKDB Sampling Position E
242
243 // Calculate Two-body quantities
244
245 meanfield->Cal2BodyQuantities();
246 std::vector< G4double > rho_a ( GetMassN
247 std::vector< G4double > rho_s ( GetMassN
248 std::vector< G4double > rho_c ( GetMassN
249
250 rho_l.resize ( GetMassNumber() , 0.0 );
251 d_pot.resize ( GetMassNumber() , 0.0 );
252
253 for ( G4int i = 0 ; i < GetMassNumber()
254 {
255 for ( G4int j = 0 ; j < GetMassNumber
256 {
257
258 rho_a[ i ] += meanfield->GetRHA( j
259 G4int k = 0;
260 if ( participants[j]->GetDefinitio
261 {
262 k = 1;
263 }
264
265 rho_s[ i ] += meanfield->GetRHA( j
266
267 rho_c[ i ] += meanfield->GetRHE( j
268 }
269
270 }
271
272 for ( G4int i = 0 ; i < GetMassNumber()
273 {
274 rho_l[i] = cdp * ( rho_a[i] + 1.0 );
275 d_pot[i] = c0p * rho_a[i]
276 + c3p * G4Pow::GetInstance()
277 + csp * rho_s[i]
278 + clp * rho_c[i];
279
280 //std::cout << "d_pot[i] " << i << "
281 }
282
283
284 // Sampling Momentum
285
286 //std::cout << "TKDB Sampling Momentum "
287
288 phase_g.clear();
289 phase_g.resize( GetMassNumber() , 0.0 );
290
291 //std::cout << "TKDB Sampling Momentum 1
292 G4bool isThis1stMOK = false;
293 G4int nmTry = 0;
294 while ( nmTry < maxTrial ) // Loop check
295 {
296 nmTry++;
297 G4int i = 0;
298 if ( samplingMomentum( i ) )
299 {
300 isThis1stMOK = true;
301 break;
302 }
303 }
304 if ( isThis1stMOK == false ) continue;
305
306 //std::cout << "TKDB Sampling Momentum 2
307
308 G4bool areTheseMsOK = true;
309 nmTry = 0;
310 while ( nmTry < maxTrial ) // Loop check
311 {
312 nmTry++;
313 G4int i = 0;
314 for ( i = 1 ; i < GetMassNumber()
315 {
316 //std::cout << "TKDB packNucleo
317 if ( !( samplingMomentum( i ) )
318 {
319 //std::cout << "TKDB packNuc
320 areTheseMsOK = false;
321 break;
322 }
323 }
324 if ( i == GetMassNumber() )
325 {
326 areTheseMsOK = true;
327 }
328
329 break;
330 }
331 if ( areTheseMsOK == false ) continue;
332
333 // Kill Angluar Momentum
334
335 //std::cout << "TKDB Sampling Kill Anglu
336
337 killCMMotionAndAngularM();
338
339
340 // Check Binding Energy
341
342 //std::cout << "packNucleons Check Bindi
343 meanfield->Cal2BodyQuantities();
344 G4double etotal = meanfield->GetTotalEne
345 G4double totalmass = 0.0;
346 // G4double etotal = 0.0;
347 //G4double ekinal = 0.0;
348 for ( int i = 0 ; i < GetMassNumber() ;
349 {
350 // ekinal += participants[i]->GetKinet
351 // G4double emass = participants[i]->
352 // G4double ekinal2 = participants[i]
353 // ekinal2 = ekinal2 * ekinal2;
354 // //etotal += std::sqrt(ekinal2 + 2*
355 // etotal += meanfield->GetSingleEner
356 totalmass += participants[i]->GetMas
357 }
358
359 //G4double totalPotentialE = meanfield->
360 //G4double ebinal = ( totalPotentialE +
361 G4double ebinal = ( etotal-totalmass ) /
362
363 //std::cout << "packNucleons totalPotent
364 //std::cout << "packNucleons ebinal " <<
365 //std::cout << "packNucleons ekinal " <<
366
367 if ( ebinal < ebin0 || ebinal > ebin1 )
368 {
369 //std::cout << "packNucleons ebin0 "
370 //std::cout << "packNucleons ebin1 "
371 //std::cout << "packNucleons ebinal "
372 //std::cout << "packNucleons Check Bindi
373 continue;
374 }
375
376 //std::cout << "packNucleons Check Bindi
377
378
379 // Energy Adujstment
380
381 G4double dtc = 1.0;
382 G4double frg = -0.1;
383 G4double rdf0 = 0.5;
384
385 G4double edif0 = ebinal - ebin00;
386
387 G4double cfrc = 0.0;
388 if ( 0 < edif0 )
389 {
390 cfrc = frg;
391 }
392 else
393 {
394 cfrc = -frg;
395 }
396
397 G4int ifrc = 1;
398
399 G4int neaTry = 0;
400
401 G4bool isThisEAOK = false;
402 while ( neaTry < maxTrial ) // Loop che
403 {
404 neaTry++;
405
406 G4double edif = ebinal - ebin00;
407
408 //std::cout << "TKDB edif " << edif <
409 if ( std::abs ( edif ) < epse )
410 {
411
412 isThisEAOK = true;
413 //std::cout << "isThisEAOK " << is
414 break;
415 }
416
417 G4int jfrc = 0;
418 if ( edif < 0.0 )
419 {
420 jfrc = 1;
421 }
422 else
423 {
424 jfrc = -1;
425 }
426
427 if ( jfrc != ifrc )
428 {
429 cfrc = -rdf0 * cfrc;
430 dtc = rdf0 * dtc;
431 }
432
433 if ( jfrc == ifrc && std::abs( edif0
434 {
435 cfrc = -rdf0 * cfrc;
436 dtc = rdf0 * dtc;
437 }
438
439 ifrc = jfrc;
440 edif0 = edif;
441
442 meanfield->CalGraduate();
443
444 for ( int i = 0 ; i < GetMassNumber()
445 {
446 G4ThreeVector ri = participants[i]
447 G4ThreeVector p_i = participants[i
448
449 ri += dtc * ( meanfield->GetFFr(i)
450 p_i += dtc * ( meanfield->GetFFp(i
451
452 participants[i]->SetPosition( ri )
453 participants[i]->SetMomentum( p_i
454 }
455
456 //ekinal = 0.0;
457 etotal = 0.0;
458
459 meanfield->Cal2BodyQuantities();
460 etotal = meanfield->GetTotalEnergy()
461
462
463 //for ( int i = 0 ; i < GetMassNumber
464 //{
465 // ekinal += participants[i]->GetKi
466 //}
467
468 //meanfield->Cal2BodyQuantities();
469 //totalPotentialE = meanfield->GetTot
470
471 //ebinal = ( totalPotentialE + ekinal
472
473
474 energy_QMD = (etotal-totalmass)/ dou
475 ebinal = energy_QMD;
476
477
478 }
479 //std::cout << "isThisEAOK " << isThisEA
480 if ( isThisEAOK == false ) continue;
481
482 isThisOK = true;
483 //std::cout << "isThisOK " << isThisOK <
484
485 break;
486
487 } // while(n0Try < maxTrial)
488
489 if ( isThisOK == false )
490 {
491 G4cout << "GroundStateNucleus state cann
492 }
493
494 //std::cout << "packNucleons End" << std::e
495 return;
496 }
497
498
499 G4bool G4LightIonQMDGroundStateNucleus::sampli
500 {
501
502 G4bool result = false;
503
504 G4int nTry = 0;
505
506 while ( nTry < maxTrial ) // Loop checking
507 {
508
509 //std::cout << "samplingPosition i th pa
510
511 G4double rwod = -1.0;
512 G4double rrr = 0.0;
513
514 G4double rx = 0.0;
515 G4double ry = 0.0;
516 G4double rz = 0.0;
517
518 G4int icounter = 0;
519 G4int icounter_max = 1024;
520 while ( G4UniformRand() * rmax > rwod )
521 {
522 icounter++;
523 if ( icounter > icounter_max ) {
524 G4cout << "Loop-counter exceeded the thr
525 break;
526 }
527
528 G4double rsqr = 10.0;
529 G4int jcounter = 0;
530 G4int jcounter_max = 1024;
531 while ( rsqr > 1.0 ) // Loop checking
532 {
533 jcounter++;
534 if ( jcounter > jcounter_max ) {
535 G4cout << "Loop-counter exceeded the
536 break;
537 }
538 rx = 1.0 - 2.0 * G4UniformRand();
539 ry = 1.0 - 2.0 * G4UniformRand();
540 rz = 1.0 - 2.0 * G4UniformRand();
541 rsqr = rx*rx + ry*ry + rz*rz;
542 }
543 rrr = radm * std::sqrt ( rsqr );
544 rwod = 1.0 / ( 1.0 + G4Exp ( ( rrr -
545
546 }
547
548 participants[i]->SetPosition( G4ThreeVec
549
550 if ( i == 0 )
551 {
552 result = true;
553 return result;
554 }
555
556 // i > 1 ( Second Particle or later )
557 // Check Distance to others
558
559 G4bool isThisOK = true;
560 for ( G4int j = 0 ; j < i ; j++ )
561 {
562
563 G4double r2 = participants[j]->GetPo
564 G4double dmin2 = 0.0;
565
566 if ( participants[j]->GetDefinition()
567 {
568 dmin2 = dsam2;
569 }
570 else
571 {
572 dmin2 = ddif2;
573 }
574
575 //std::cout << "distance between j an
576 if ( r2 < dmin2 )
577 {
578 isThisOK = false;
579 break;
580 }
581
582 }
583
584 if ( isThisOK == true )
585 {
586 result = true;
587 return result;
588 }
589
590 nTry++;
591
592 }
593
594 // Here return "false"
595 return result;
596 }
597
598
599
600 G4bool G4LightIonQMDGroundStateNucleus::sampli
601 {
602
603 //std::cout << "TKDB samplingMomentum for "
604
605 G4bool result = false;
606
607 G4double pfm = hbc * G4Pow::GetInstance()->
608
609 if ( 10 < GetMassNumber() && -5.5 < ebini
610 {
611 pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std
612 }
613
614 //std::cout << "TKDB samplingMomentum pfm "
615
616 std::vector< G4double > phase;
617 phase.resize( i+1 ); // i start from 0
618
619 G4int ntry = 0;
620 // 710
621 while ( ntry < maxTrial ) // Loop checking
622 {
623
624 //std::cout << " TKDB ntry " << ntry <<
625 ntry++;
626
627 //G4double ke = DBL_MAX;
628
629 G4int tkdb_i =0;
630 // 700
631 //G4int icounter = 0;
632 //G4int icounter_max = 1024;
633
634
635 //while ( ke + d_pot [i] > edepth ) // L
636 //{
637 // icounter++;
638 // if ( icounter > icounter_max ) {
639 // G4cout << "Loop-counter exceeded the t
640 // break;
641 // }
642
643 G4double psqr = 10.0;
644 G4double px = 0.0;
645 G4double py = 0.0;
646 G4double pz = 0.0;
647
648 G4int jcounter = 0;
649 G4int jcounter_max = 1024;
650 while ( psqr > 1.0 ) // Loop checking
651 {
652 jcounter++;
653 if ( jcounter > jcounter_max ) {
654 G4cout << "Loop-counter exceeded the
655 break;
656 }
657 px = 1.0 - 2.0*G4UniformRand();
658 py = 1.0 - 2.0*G4UniformRand();
659 pz = 1.0 - 2.0*G4UniformRand();
660
661 psqr = px*px + py*py + pz*pz;
662 }
663
664 G4ThreeVector p ( px , py , pz );
665 p = pfm * p;
666 participants[i]->SetMomentum( p );
667 G4LorentzVector p4 = participants[i]-
668 //ke = p4.e() - p4.restMass();
669 //ke = participants[i]->GetKineticEne
670
671
672 tkdb_i++;
673 if ( tkdb_i > maxTrial ) return resul
674
675 //}
676
677 //std::cout << "TKDB ke d_pot[i] " << ke
678
679
680 if ( i == 0 )
681 {
682 result = true;
683 return result;
684 }
685
686 G4bool isThisOK = true;
687
688 // Check Pauli principle
689
690 phase[ i ] = 0.0;
691
692 //std::cout << "TKDB Check Pauli Princip
693
694 for ( G4int j = 0 ; j < i ; j++ )
695 {
696 phase[ j ] = 0.0;
697 //std::cout << "TKDB Check Pauli Prin
698 G4double expa = 0.0;
699 if ( participants[j]->GetDefinition()
700 {
701
702 expa = - meanfield->GetRR2(i,j) *
703
704 if ( expa > epsx )
705 {
706 G4ThreeVector p_i = participant
707 G4ThreeVector pj = participants
708 G4double dist2_p = p_i.diff2( p
709
710 dist2_p = dist2_p*cph;
711 expa = expa - dist2_p;
712
713 if ( expa > epsx )
714 {
715
716 phase[j] = G4Exp ( expa );
717
718 if ( phase[j] * cpc > 0.2 )
719 {
720 /*
721 G4cout << "TKDB Check Pauli Principle
722 G4cout << "TKDB Check Pauli Principle
723 G4cout << "TKDB Check Pauli Principle
724 */
725 isThisOK = false;
726 break;
727 }
728 if ( ( phase_g[j] + phase[j]
729 {
730 /*
731 G4cout << "TKDB Check Pauli Principle
732 G4cout << "TKDB Check Pauli Principle
733 G4cout << "TKDB Check Pauli Principle
734 G4cout << "TKDB Check Pauli Principle
735 */
736 isThisOK = false;
737 break;
738 }
739
740 phase[i] += phase[j];
741 if ( phase[i] * cpc > 0.3 )
742 {
743 /*
744 G4cout << "TKDB Check Pauli Principle
745 G4cout << "TKDB Check Pauli Principle
746 G4cout << "TKDB Check Pauli Principle
747 */
748 isThisOK = false;
749 break;
750 }
751
752 //std::cout << "TKDB Check P
753
754 }
755 else
756 {
757 //std::cout << "TKDB Check P
758 }
759
760 }
761 else
762 {
763 //std::cout << "TKDB Check Paul
764 }
765
766 }
767 else
768 {
769 //std::cout << "TKDB Check Pauli P
770 }
771
772 }
773
774 if ( isThisOK == true )
775 {
776
777 phase_g[i] = phase[i];
778
779 for ( int j = 0 ; j < i ; j++ )
780 {
781 phase_g[j] += phase[j];
782 }
783
784 result = true;
785 return result;
786 }
787
788 }
789
790 return result;
791
792 }
793
794
795
796 void G4LightIonQMDGroundStateNucleus::killCMMo
797 {
798
799 // CalEnergyAndAngularMomentumInCM();
800
801 //std::vector< G4ThreeVector > p ( GetMassN
802 //std::vector< G4ThreeVector > r ( GetMassN
803
804 // Move to cm system
805
806 G4ThreeVector pcm_tmp ( 0.0 );
807 G4ThreeVector rcm_tmp ( 0.0 );
808 G4double sumMass = 0.0;
809
810 for ( G4int i = 0 ; i < GetMassNumber() ; i
811 {
812 pcm_tmp += participants[i]->GetMomentum(
813 rcm_tmp += participants[i]->GetPosition(
814 sumMass += participants[i]->GetMass();
815 }
816
817 pcm_tmp = pcm_tmp/GetMassNumber();
818 rcm_tmp = rcm_tmp/sumMass;
819
820 for ( G4int i = 0 ; i < GetMassNumber() ; i
821 {
822 participants[i]->SetMomentum( participan
823 participants[i]->SetPosition( participan
824 }
825
826 // kill the angular momentum
827
828 G4ThreeVector ll ( 0.0 );
829 for ( G4int i = 0 ; i < GetMassNumber() ; i
830 {
831 ll += participants[i]->GetPosition().cro
832 }
833
834 G4double rr[3][3];
835 G4double ss[3][3];
836 for ( G4int i = 0 ; i < 3 ; i++ )
837 {
838 for ( G4int j = 0 ; j < 3 ; j++ )
839 {
840 rr [i][j] = 0.0;
841
842 if ( i == j )
843 {
844 ss [i][j] = 1.0;
845 }
846 else
847 {
848 ss [i][j] = 0.0;
849 }
850 }
851 }
852
853 for ( G4int i = 0 ; i < GetMassNumber() ; i
854 {
855 G4ThreeVector r = participants[i]->GetPo
856 rr[0][0] += r.y() * r.y() + r.z() * r.z(
857 rr[1][0] += - r.y() * r.x();
858 rr[2][0] += - r.z() * r.x();
859 rr[0][1] += - r.x() * r.y();
860 rr[1][1] += r.z() * r.z() + r.x() * r.x(
861 rr[2][1] += - r.z() * r.y();
862 rr[2][0] += - r.x() * r.z();
863 rr[2][1] += - r.y() * r.z();
864 rr[2][2] += r.x() * r.x() + r.y() * r.y(
865 }
866
867 for ( G4int i = 0 ; i < 3 ; i++ )
868 {
869 G4double x = rr [i][i];
870 for ( G4int j = 0 ; j < 3 ; j++ )
871 {
872 rr[i][j] = rr[i][j] / x;
873 ss[i][j] = ss[i][j] / x;
874 }
875 for ( G4int j = 0 ; j < 3 ; j++ )
876 {
877 if ( j != i )
878 {
879 G4double y = rr [j][i];
880 for ( G4int k = 0 ; k < 3 ; k++ )
881 {
882 rr[j][k] += -y * rr[i][k];
883 ss[j][k] += -y * ss[i][k];
884 }
885 }
886 }
887 }
888
889 G4double opl[3];
890 G4double rll[3];
891
892 rll[0] = ll.x();
893 rll[1] = ll.y();
894 rll[2] = ll.z();
895
896 for ( G4int i = 0 ; i < 3 ; i++ )
897 {
898 opl[i] = 0.0;
899
900 for ( G4int j = 0 ; j < 3 ; j++ )
901 {
902 opl[i] += ss[i][j]*rll[j];
903 }
904 }
905
906 for ( G4int i = 0 ; i < GetMassNumber() ; i
907 {
908 G4ThreeVector p_i = participants[i]->Get
909 G4ThreeVector ri = participants[i]->GetP
910 G4ThreeVector opl_v ( opl[0] , opl[1] ,
911
912 p_i += ri.cross(opl_v);
913
914 participants[i]->SetMomentum( p_i );
915 }
916
917 }
918
919 G4bool G4LightIonQMDGroundStateNucleus::sampli
920 {
921 std::vector < G4ThreeVector > base_x;
922 base_x.clear();
923 base_x.resize( 4 );
924
925 base_x[0] = G4ThreeVector( 1.0 , 1.0 , 1.0
926 base_x[1] = G4ThreeVector( -1.0 , -1.0 , 1
927 base_x[2] = G4ThreeVector( 1.0 , -1.0 , -1
928 base_x[3] = G4ThreeVector( -1.0 , 1.0 , -1
929
930 G4double al = 2 * pi * G4UniformRand();
931 G4double be = 2 * pi * G4UniformRand();
932 G4double ga = 2 * pi * G4UniformRand();
933 G4double sca = std::cos(al);
934 G4double ssa = std::sin(al);
935 G4double scb = std::cos(be);
936 G4double ssb = std::sin(be);
937 G4double scg = std::cos(ga);
938 G4double ssg = std::sin(ga);
939
940 if (z == 6)
941 {
942 std::vector < G4ThreeVector > sub_x;
943 sub_x.clear();
944 sub_x.resize( 3 );
945
946 sub_x[0] = G4ThreeVector( 1.0/std::sqr
947 sub_x[1] = G4ThreeVector( -0.5/std::sq
948 sub_x[2] = G4ThreeVector( -0.5/std::sq
949 G4double sbfac = 2.5 + 0.4*G4UniformRa
950 G4double safac = 0.5;
951
952 for ( G4int j = 0 ; j < 3 ; j++ )
953 {
954 G4double xx = (sca*scb*scg-ssa*ssg
955 G4double yy = (ssa*scb*scg+sca*ssg
956 G4double zz = (-ssb*scg)*sub_x[j].
957 G4ThreeVector sprime = G4ThreeVect
958
959 al = 2 * pi * G4UniformRand();
960 be = 2 * pi * G4UniformRand();
961 ga = 2 * pi * G4UniformRand();
962 G4double ca = std::cos(al);
963 G4double sa = std::sin(al);
964 G4double cb = std::cos(be);
965 G4double sb = std::sin(be);
966 G4double cg = std::cos(ga);
967 G4double sg = std::sin(ga);
968
969 for ( G4int i = 0 ; i < 4 ; i++ )
970 {
971 xx = (ca*cb*cg-sa*sg)*base_x[i
972 yy = (sa*cb*cg+ca*sg)*base_x[i
973 zz = (-sb*cg)*base_x[i].x() +
974 G4ThreeVector rprime = G4Three
975 //std::cout << "r base " << rp
976 //std::cout << "s base " << sp
977 G4ThreeVector pos = safac*rpri
978 //std::cout << GetParticipant(
979 participants[3*i+j]->SetPositi
980 }
981 }
982 }
983 else if(z == 8)
984 {
985 std::vector < G4ThreeVector > sub_x;
986 sub_x.clear();
987 sub_x.resize( 4 );
988
989 sub_x[0] = G4ThreeVector( 1.0 , 1.0 ,
990 sub_x[1] = G4ThreeVector( -1.0 , -1.0
991 sub_x[2] = G4ThreeVector( 1.0 , -1.0 ,
992 sub_x[3] = G4ThreeVector( -1.0 , 1.0 ,
993 G4double sbfac = 1.75 + 0.25*G4Uniform
994 G4double safac = 0.50;
995
996 for ( G4int j = 0 ; j < 4 ; j++ )
997 {
998 G4double xx = (sca*scb*scg-ssa*ssg
999 G4double yy = (ssa*scb*scg+sca*ssg
1000 G4double zz = (-ssb*scg)*sub_x[j]
1001 G4ThreeVector sprime = G4ThreeVec
1002
1003 al = 2 * pi * G4UniformRand();
1004 be = 2 * pi * G4UniformRand();
1005 ga = 2 * pi * G4UniformRand();
1006 G4double ca = std::cos(al);
1007 G4double sa = std::sin(al);
1008 G4double cb = std::cos(be);
1009 G4double sb = std::sin(be);
1010 G4double cg = std::cos(ga);
1011 G4double sg = std::sin(ga);
1012
1013 for ( G4int i = 0 ; i < 4 ; i++ )
1014 {
1015
1016 xx = (ca*cb*cg-sa*sg)*base_x[
1017 yy = (sa*cb*cg+ca*sg)*base_x[
1018 zz = (-sb*cg)*base_x[i].x() +
1019 G4ThreeVector rprime = G4Thre
1020
1021
1022 G4ThreeVector pos = safac*rpr
1023 //std::cout << GetParticipant
1024 participants[4*i+j]->SetPosit
1025 }
1026 }
1027 }
1028
1029 bool result = true;
1030 return result;
1031 }
1032