| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
27 //
28 // 230308 "samplingPosition_cluster" function added by S-H. Sato and A. Haga
29 // 230308 Parameter "radam" shold be independent of parameter "gamm" (S-H. Sato and A. Haga)
30 // 230308 Sampling momentum of the particle does not consider its binding energy (S-H. Sato and A. Haga)
31 // 230308 Binding energy evaluated by "GetTotalEnergy" calculated in Mean Field (S-H. Sato and A. Haga)
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::G4LightIonQMDGroundStateNucleus( G4int z , G4int a )
44 : maxTrial ( 1000 )
45 , r00 ( 1.124 ) // radius parameter for Woods-Saxon [fm]
46 , r01 ( 0.5 ) // radius parameter for Woods-Saxon
47 , saa ( 0.2 ) // diffuse parameter for initial Woods-Saxon shape
48 , rada ( 0.9 ) // cutoff parameter
49 , radb ( 0.3 ) // cutoff parameter
50 , dsam ( 1.5 ) // minimum distance for same particle [fm]
51 , ddif ( 1.0 ) // minimum distance for different particle
52 , edepth ( 0.0 )
53 , epse ( 0.000001 ) // torelance for energy in [GeV]
54 , meanfield ( NULL )
55 {
56
57 //std::cout << " G4LightIonQMDGroundStateNucleus( G4int z , G4int a ) Begin " << z << " " << a << std::endl;
58
59 dsam2 = dsam*dsam;
60 ddif2 = ddif*ddif;
61
62 G4LightIonQMDParameters* parameters = G4LightIonQMDParameters::GetInstance();
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 before the line 90.
79 // Otherwise, mass number cannot be conserved if the projectile or
80 // the target are nucleons.
81 //Nucleon primary or target case;
82 if ( z == 1 && a == 1 ) { // Hydrogen Case or proton primary
83 SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
84 // ebini = 0.0;
85 return;
86 }
87 else if ( z == 0 && a == 1 ) { // Neutron primary
88 SetParticipant( new G4QMDParticipant( G4Neutron::Neutron() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
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 G4QMDParticipant( pd , p , r );
113 SetParticipant( aParticipant );
114
115 }
116
117 G4double radious = r00 * G4Pow::GetInstance()->A13( double ( GetMassNumber() ) );
118
119 rt00 = radious - r01;
120 radm = radious; // JQMD, Skyrme, RMF, Y-H. Sato and A. Haga, 20230308
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 << "G4LightIonQMDGroundStateNucleus( G4int z , G4int a ) packNucleons Begin ( z , a ) ( " << z << ", " << a << " )" << std::endl;
130 packNucleons();
131 //std::cout << "G4LightIonQMDGroundStateNucleus( G4int z , G4int a ) packNucleons End" << std::endl;
132
133 delete meanfield;
134
135 }
136
137
138
139 void G4LightIonQMDGroundStateNucleus::packNucleons()
140 {
141
142 //std::cout << "G4LightIonQMDGroundStateNucleus::packNucleons" << std::endl;
143
144 ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber();
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.H. 20230308
165 G4int Nucle_Z = GetAtomicNumber(); // Y-H.S. and A.H. 20230308
166 G4int Nucle_A = GetMassNumber(); // Y-H.S. and A.H. 20230308
167
168 while ( n0Try < maxTrial ) // Loop checking, 11.03.2015, T. Koi
169 {
170 n0Try++;
171 //std::cout << "TKDB packNucleons n0Try " << n0Try << std::endl;
172
173 // Sampling Position
174
175 //std::cout << "TKDB Sampling Position " << std::endl;
176
177 G4bool areThesePsOK = false;
178 G4int npTry = 0;
179 while ( npTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
180 {
181 //std::cout << "TKDB Sampling Position npTry " << npTry << std::endl;
182 if( (Nucle_Z == 6 && Nucle_A == 12) || (Nucle_Z == 8 && Nucle_A == 16))
183 {
184 //std::cout << "alpha cluster " << Nucle_Z << " ," << Nucle_A <<std::endl;
185 //npTry++;
186 //G4int i = 0;
187 if ( samplingPosition_cluster( Nucle_Z ) )
188 {
189 //G4double rsquare = meanfield->CalDensityProfile();
190 //G4cout << "R-square " << std::sqrt(rsquare) << G4endl;
191 areThesePsOK = true;
192 break;
193 //exit(0);
194 /*
195 //std::cout << "packNucleons samplingPosition 0 succeed " << std::endl;
196 for ( i = 1 ; i < GetMassNumber() ; i++ )
197 {
198 //std::cout << "packNucleons samplingPosition " << i << " trying " << std::endl;
199 if ( !( samplingPosition_cluster( Nucle_Z, Nucle_A ) ) )
200 {
201 //std::cout << "packNucleons samplingPosition " << i << " failed" << std::endl;
202 break;
203 }
204 }
205 if ( i == GetMassNumber() )
206 {
207 //std::cout << "packNucleons samplingPosition all scucceed " << std::endl;
208 areThesePsOK = true;
209 break;
210 }
211 */
212 }
213 } // Alpha cluster model added by Y-H.S. and A.H. 20230308
214 else
215 {
216 npTry++;
217 G4int i = 0;
218 if ( samplingPosition( i ) )
219 {
220 //std::cout << "packNucleons samplingPosition 0 succeed " << std::endl;
221 for ( i = 1 ; i < GetMassNumber() ; i++ )
222 {
223 //std::cout << "packNucleons samplingPosition " << i << " trying " << std::endl;
224 if ( !( samplingPosition( i ) ) )
225 {
226 //std::cout << "packNucleons samplingPosition " << i << " failed" << std::endl;
227 break;
228 }
229 }
230 if ( i == GetMassNumber() )
231 {
232 //std::cout << "packNucleons samplingPosition all scucceed " << std::endl;
233 areThesePsOK = true;
234 break;
235 }
236 }
237 }
238 }
239 if ( areThesePsOK == false ) continue;
240
241 //std::cout << "TKDB Sampling Position End" << std::endl;
242
243 // Calculate Two-body quantities
244
245 meanfield->Cal2BodyQuantities();
246 std::vector< G4double > rho_a ( GetMassNumber() , 0.0 );
247 std::vector< G4double > rho_s ( GetMassNumber() , 0.0 );
248 std::vector< G4double > rho_c ( GetMassNumber() , 0.0 );
249
250 rho_l.resize ( GetMassNumber() , 0.0 );
251 d_pot.resize ( GetMassNumber() , 0.0 );
252
253 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
254 {
255 for ( G4int j = 0 ; j < GetMassNumber() ; j++ )
256 {
257
258 rho_a[ i ] += meanfield->GetRHA( j , i );
259 G4int k = 0;
260 if ( participants[j]->GetDefinition() != participants[i]->GetDefinition() )
261 {
262 k = 1;
263 }
264
265 rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k ); // OK?
266
267 rho_c[ i ] += meanfield->GetRHE( j , i );
268 }
269
270 }
271
272 for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
273 {
274 rho_l[i] = cdp * ( rho_a[i] + 1.0 );
275 d_pot[i] = c0p * rho_a[i]
276 + c3p * G4Pow::GetInstance()->powA ( rho_a[i] , gamm )
277 + csp * rho_s[i]
278 + clp * rho_c[i];
279
280 //std::cout << "d_pot[i] " << i << " " << d_pot[i] << std::endl;
281 }
282
283
284 // Sampling Momentum
285
286 //std::cout << "TKDB Sampling Momentum " << std::endl;
287
288 phase_g.clear();
289 phase_g.resize( GetMassNumber() , 0.0 );
290
291 //std::cout << "TKDB Sampling Momentum 1st " << std::endl;
292 G4bool isThis1stMOK = false;
293 G4int nmTry = 0;
294 while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
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 2nd so on" << std::endl;
307
308 G4bool areTheseMsOK = true;
309 nmTry = 0;
310 while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
311 {
312 nmTry++;
313 G4int i = 0;
314 for ( i = 1 ; i < GetMassNumber() ; i++ )
315 {
316 //std::cout << "TKDB packNucleons samplingMomentum try " << i << std::endl;
317 if ( !( samplingMomentum( i ) ) )
318 {
319 //std::cout << "TKDB packNucleons samplingMomentum " << i << " failed" << std::endl;
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 Angluar Momentum " << std::endl;
336
337 killCMMotionAndAngularM();
338
339
340 // Check Binding Energy
341
342 //std::cout << "packNucleons Check Binding Energy Begin " << std::endl;
343 meanfield->Cal2BodyQuantities();
344 G4double etotal = meanfield->GetTotalEnergy();
345 G4double totalmass = 0.0;
346 // G4double etotal = 0.0;
347 //G4double ekinal = 0.0;
348 for ( int i = 0 ; i < GetMassNumber() ; i++ )
349 {
350 // ekinal += participants[i]->GetKineticEnergy();
351 // G4double emass = participants[i]->GetMass();
352 // G4double ekinal2 = participants[i]->GetKineticEnergy() + emass;
353 // ekinal2 = ekinal2 * ekinal2;
354 // //etotal += std::sqrt(ekinal2 + 2*emass* meanfield->GetPotential(i)) - emass;
355 // etotal += meanfield->GetSingleEnergy(i) -emass;
356 totalmass += participants[i]->GetMass();
357 }
358
359 //G4double totalPotentialE = meanfield->GetTotalPotential();
360 //G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
361 G4double ebinal = ( etotal-totalmass ) / double ( GetMassNumber() );
362
363 //std::cout << "packNucleons totalPotentialE " << totalPotentialE << std::endl;
364 //std::cout << "packNucleons ebinal " << ebinal << std::endl;
365 //std::cout << "packNucleons ekinal " << ekinal << std::endl;
366
367 if ( ebinal < ebin0 || ebinal > ebin1 )
368 {
369 //std::cout << "packNucleons ebin0 " << ebin0 << std::endl;
370 //std::cout << "packNucleons ebin1 " << ebin1 << std::endl;
371 //std::cout << "packNucleons ebinal " << ebinal << std::endl;
372 //std::cout << "packNucleons Check Binding Energy Failed " << std::endl;
373 continue;
374 }
375
376 //std::cout << "packNucleons Check Binding Energy End = OK " << std::endl;
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 checking, 11.03.2015, T. Koi
403 {
404 neaTry++;
405
406 G4double edif = ebinal - ebin00;
407
408 //std::cout << "TKDB edif " << edif << std::endl;
409 if ( std::abs ( edif ) < epse )
410 {
411
412 isThisEAOK = true;
413 //std::cout << "isThisEAOK " << isThisEAOK << std::endl;
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 ) < std::abs( edif ) )
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() ; i++ )
445 {
446 G4ThreeVector ri = participants[i]->GetPosition();
447 G4ThreeVector p_i = participants[i]->GetMomentum();
448
449 ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) );
450 p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(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() ; i++ )
464 //{
465 // ekinal += participants[i]->GetKineticEnergy();
466 //}
467
468 //meanfield->Cal2BodyQuantities();
469 //totalPotentialE = meanfield->GetTotalPotential();
470
471 //ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
472
473
474 energy_QMD = (etotal-totalmass)/ double ( GetMassNumber() );
475 ebinal = energy_QMD;
476
477
478 }
479 //std::cout << "isThisEAOK " << isThisEAOK << std::endl;
480 if ( isThisEAOK == false ) continue;
481
482 isThisOK = true;
483 //std::cout << "isThisOK " << isThisOK << std::endl;
484
485 break;
486
487 } // while(n0Try < maxTrial)
488
489 if ( isThisOK == false )
490 {
491 G4cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << G4endl;
492 }
493
494 //std::cout << "packNucleons End" << std::endl;
495 return;
496 }
497
498
499 G4bool G4LightIonQMDGroundStateNucleus::samplingPosition( G4int i )
500 {
501
502 G4bool result = false;
503
504 G4int nTry = 0;
505
506 while ( nTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
507 {
508
509 //std::cout << "samplingPosition i th particle, nTtry " << i << " " << nTry << std::endl;
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 ) // Loop checking, 11.03.2015, T. Koi
521 {
522 icounter++;
523 if ( icounter > icounter_max ) {
524 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
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, 11.03.2015, T. Koi
532 {
533 jcounter++;
534 if ( jcounter > jcounter_max ) {
535 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
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 - rt00 ) / saa ) );
545
546 }
547
548 participants[i]->SetPosition( G4ThreeVector( rx , ry , rz )*radm );
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]->GetPosition().diff2( participants[i]->GetPosition() );
564 G4double dmin2 = 0.0;
565
566 if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
567 {
568 dmin2 = dsam2;
569 }
570 else
571 {
572 dmin2 = ddif2;
573 }
574
575 //std::cout << "distance between j and i " << j << " " << i << " " << r2 << " " << dmin2 << std::endl;
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::samplingMomentum( G4int i )
601 {
602
603 //std::cout << "TKDB samplingMomentum for " << i << std::endl;
604
605 G4bool result = false;
606
607 G4double pfm = hbc * G4Pow::GetInstance()->A13 ( ( 3.0 / 2.0 * pi*pi * rho_l[i] ) );
608
609 if ( 10 < GetMassNumber() && -5.5 < ebini )
610 {
611 pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std::abs( 8.0 + ebini ) / 8.0 ) );
612 }
613
614 //std::cout << "TKDB samplingMomentum pfm " << pfm << std::endl;
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, 11.03.2015, T. Koi
622 {
623
624 //std::cout << " TKDB ntry " << ntry << std::endl;
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 ) // Loop checking, 11.03.2015, T. Koi
636 //{
637 // icounter++;
638 // if ( icounter > icounter_max ) {
639 // G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
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, 11.03.2015, T. Koi
651 {
652 jcounter++;
653 if ( jcounter > jcounter_max ) {
654 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
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]->Get4Momentum();
668 //ke = p4.e() - p4.restMass();
669 //ke = participants[i]->GetKineticEnergy();
670
671
672 tkdb_i++;
673 if ( tkdb_i > maxTrial ) return result; // return false
674
675 //}
676
677 //std::cout << "TKDB ke d_pot[i] " << ke << " " << d_pot[i] << std::endl;
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 Principle " << i << std::endl;
693
694 for ( G4int j = 0 ; j < i ; j++ )
695 {
696 phase[ j ] = 0.0;
697 //std::cout << "TKDB Check Pauli Principle i , j " << i << " , " << j << std::endl;
698 G4double expa = 0.0;
699 if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
700 {
701
702 expa = - meanfield->GetRR2(i,j) * cpw;
703
704 if ( expa > epsx )
705 {
706 G4ThreeVector p_i = participants[i]->GetMomentum();
707 G4ThreeVector pj = participants[j]->GetMomentum();
708 G4double dist2_p = p_i.diff2( pj );
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 A i , j " << i << " , " << j << G4endl;
722 G4cout << "TKDB Check Pauli Principle phase[j] " << phase[j] << G4endl;
723 G4cout << "TKDB Check Pauli Principle phase[j]*cpc > 0.2 " << phase[j]*cpc << G4endl;
724 */
725 isThisOK = false;
726 break;
727 }
728 if ( ( phase_g[j] + phase[j] ) * cpc > 0.5 )
729 {
730 /*
731 G4cout << "TKDB Check Pauli Principle B i , j " << i << " , " << j << G4endl;
732 G4cout << "TKDB Check Pauli Principle B phase_g[j] " << phase_g[j] << G4endl;
733 G4cout << "TKDB Check Pauli Principle B phase[j] " << phase[j] << G4endl;
734 G4cout << "TKDB Check Pauli Principle B phase_g[j] + phase[j] ) * cpc > 0.5 " << ( phase_g[j] + phase[j] ) * cpc << G4endl;
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 C i , j " << i << " , " << j << G4endl;
745 G4cout << "TKDB Check Pauli Principle C phase[i] " << phase[i] << G4endl;
746 G4cout << "TKDB Check Pauli Principle C phase[i] * cpc > 0.3 " << phase[i] * cpc << G4endl;
747 */
748 isThisOK = false;
749 break;
750 }
751
752 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
753
754 }
755 else
756 {
757 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
758 }
759
760 }
761 else
762 {
763 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
764 }
765
766 }
767 else
768 {
769 //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
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::killCMMotionAndAngularM()
797 {
798
799 // CalEnergyAndAngularMomentumInCM();
800
801 //std::vector< G4ThreeVector > p ( GetMassNumber() , 0.0 );
802 //std::vector< G4ThreeVector > r ( GetMassNumber() , 0.0 );
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() * participants[i]->GetMass();
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( participants[i]->GetMomentum() - pcm_tmp );
823 participants[i]->SetPosition( participants[i]->GetPosition() - rcm_tmp );
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().cross ( participants[i]->GetMomentum() );
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]->GetPosition();
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]->GetMomentum() ;
909 G4ThreeVector ri = participants[i]->GetPosition() ;
910 G4ThreeVector opl_v ( opl[0] , opl[1] , opl[2] );
911
912 p_i += ri.cross(opl_v);
913
914 participants[i]->SetMomentum( p_i );
915 }
916
917 }
918
919 G4bool G4LightIonQMDGroundStateNucleus::samplingPosition_cluster( G4int z )
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 )/std::sqrt( 3.0 );
926 base_x[1] = G4ThreeVector( -1.0 , -1.0 , 1.0 )/std::sqrt( 3.0 );
927 base_x[2] = G4ThreeVector( 1.0 , -1.0 , -1.0 )/std::sqrt( 3.0 );
928 base_x[3] = G4ThreeVector( -1.0 , 1.0 , -1.0 )/std::sqrt( 3.0 );
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::sqrt( 3.0 ), 0.0 , 0.0 );
947 sub_x[1] = G4ThreeVector( -0.5/std::sqrt( 3.0 ) , 0.5 , 0.0 );
948 sub_x[2] = G4ThreeVector( -0.5/std::sqrt( 3.0 ) , -0.5 , 0.0 );
949 G4double sbfac = 2.5 + 0.4*G4UniformRand();
950 G4double safac = 0.5;
951
952 for ( G4int j = 0 ; j < 3 ; j++ )
953 {
954 G4double xx = (sca*scb*scg-ssa*ssg)*sub_x[j].x() + (-sca*scb*ssg-ssa*scg)*sub_x[j].y() + sca*ssb*sub_x[j].z();
955 G4double yy = (ssa*scb*scg+sca*ssg)*sub_x[j].x() + (-ssa*scb*ssg+sca*scg)*sub_x[j].y() + ssa*ssb*sub_x[j].z();
956 G4double zz = (-ssb*scg)*sub_x[j].x() + (ssb*ssg)*sub_x[j].y() + scb*sub_x[j].z();
957 G4ThreeVector sprime = G4ThreeVector(xx, yy, zz);
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].x() + (-ca*cb*sg-sa*cg)*base_x[i].y() + ca*sb*base_x[i].z();
972 yy = (sa*cb*cg+ca*sg)*base_x[i].x() + (-sa*cb*sg+ca*cg)*base_x[i].y() + sa*sb*base_x[i].z();
973 zz = (-sb*cg)*base_x[i].x() + (sb*sg)*base_x[i].y() + cb*base_x[i].z();
974 G4ThreeVector rprime = G4ThreeVector(xx, yy, zz);
975 //std::cout << "r base " << rprime << " r2 " << rprime*rprime << std::endl;
976 //std::cout << "s base " << sprime << " s2 " << sprime*sprime << std::endl;
977 G4ThreeVector pos = safac*rprime + sbfac*sprime;
978 //std::cout << GetParticipant(k)->GetChargeInUnitOfEplus() << " samplingPosition " << pos << std::endl;
979 participants[3*i+j]->SetPosition( pos );
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 , 1.0 )/std::sqrt( 3.0 );
990 sub_x[1] = G4ThreeVector( -1.0 , -1.0 , 1.0 )/std::sqrt( 3.0 );
991 sub_x[2] = G4ThreeVector( 1.0 , -1.0 , -1.0 )/std::sqrt( 3.0 );
992 sub_x[3] = G4ThreeVector( -1.0 , 1.0 , -1.0 )/std::sqrt( 3.0 );
993 G4double sbfac = 1.75 + 0.25*G4UniformRand();
994 G4double safac = 0.50;
995
996 for ( G4int j = 0 ; j < 4 ; j++ )
997 {
998 G4double xx = (sca*scb*scg-ssa*ssg)*sub_x[j].x() + (-sca*scb*ssg-ssa*scg)*sub_x[j].y() + sca*ssb*sub_x[j].z();
999 G4double yy = (ssa*scb*scg+sca*ssg)*sub_x[j].x() + (-ssa*scb*ssg+sca*scg)*sub_x[j].y() + ssa*ssb*sub_x[j].z();
1000 G4double zz = (-ssb*scg)*sub_x[j].x() + (ssb*ssg)*sub_x[j].y() + scb*sub_x[j].z();
1001 G4ThreeVector sprime = G4ThreeVector(xx, yy, zz);
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[i].x() + (-ca*cb*sg-sa*cg)*base_x[i].y() + ca*sb*base_x[i].z();
1017 yy = (sa*cb*cg+ca*sg)*base_x[i].x() + (-sa*cb*sg+ca*cg)*base_x[i].y() + sa*sb*base_x[i].z();
1018 zz = (-sb*cg)*base_x[i].x() + (sb*sg)*base_x[i].y() + cb*base_x[i].z();
1019 G4ThreeVector rprime = G4ThreeVector(xx, yy, zz);
1020
1021
1022 G4ThreeVector pos = safac*rprime + sbfac*sprime;
1023 //std::cout << GetParticipant(k)->GetChargeInUnitOfEplus() << " samplingPosition " << pos << std::endl;
1024 participants[4*i+j]->SetPosition( pos );
1025 }
1026 }
1027 }
1028
1029 bool result = true;
1030 return result;
1031 }
1032