ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QMDGroundStateNucleus.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4QMDGroundStateNucleus.cc
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 #include "G4QMDGroundStateNucleus.hh"
29 
30 #include "G4NucleiProperties.hh"
31 #include "G4Proton.hh"
32 #include "G4Neutron.hh"
33 #include "G4Exp.hh"
34 #include "G4Pow.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "Randomize.hh"
37 
38 G4QMDGroundStateNucleus::G4QMDGroundStateNucleus( G4int z , G4int a )
39 : maxTrial ( 1000 )
40 , r00 ( 1.124 ) // radius parameter for Woods-Saxon [fm]
41 , r01 ( 0.5 ) // radius parameter for Woods-Saxon
42 , saa ( 0.2 ) // diffuse parameter for initial Woods-Saxon shape
43 , rada ( 0.9 ) // cutoff parameter
44 , radb ( 0.3 ) // cutoff parameter
45 , dsam ( 1.5 ) // minimum distance for same particle [fm]
46 , ddif ( 1.0 ) // minimum distance for different particle
47 , edepth ( 0.0 )
48 , epse ( 0.000001 ) // torelance for energy in [GeV]
49 , meanfield ( NULL )
50 {
51 
52  //std::cout << " G4QMDGroundStateNucleus( G4int z , G4int a ) Begin " << z << " " << a << std::endl;
53 
54  dsam2 = dsam*dsam;
55  ddif2 = ddif*ddif;
56 
57  G4QMDParameters* parameters = G4QMDParameters::GetInstance();
58 
59  hbc = parameters->Get_hbc();
60  gamm = parameters->Get_gamm();
61  cpw = parameters->Get_cpw();
62  cph = parameters->Get_cph();
63  epsx = parameters->Get_epsx();
64  cpc = parameters->Get_cpc();
65 
66  cdp = parameters->Get_cdp();
67  c0p = parameters->Get_c0p();
68  c3p = parameters->Get_c3p();
69  csp = parameters->Get_csp();
70  clp = parameters->Get_clp();
71 
72  //edepth = 0.0;
73 
74  for ( int i = 0 ; i < a ; i++ )
75  {
76 
77  G4ParticleDefinition* pd;
78 
79  if ( i < z )
80  {
81  pd = G4Proton::Proton();
82  }
83  else
84  {
85  pd = G4Neutron::Neutron();
86  }
87 
88  G4ThreeVector p( 0.0 );
89  G4ThreeVector r( 0.0 );
90  G4QMDParticipant* aParticipant = new G4QMDParticipant( pd , p , r );
91  SetParticipant( aParticipant );
92 
93  }
94 
95  G4double radious = r00 * G4Pow::GetInstance()->A13( double ( GetMassNumber() ) );
96 
97  rt00 = radious - r01;
98  radm = radious - rada * ( gamm - 1.0 ) + radb;
99  rmax = 1.0 / ( 1.0 + G4Exp ( -rt00/saa ) );
100 
101  //maxTrial = 1000;
102 
103  //Nucleon primary or target case;
104  if ( z == 1 && a == 1 ) { // Hydrogen Case or proton primary
105  SetParticipant( new G4QMDParticipant( G4Proton::Proton() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
106  ebini = 0.0;
107  return;
108  }
109  else if ( z == 0 && a == 1 ) { // Neutron primary
110  SetParticipant( new G4QMDParticipant( G4Neutron::Neutron() , G4ThreeVector( 0.0 ) , G4ThreeVector( 0.0 ) ) );
111  ebini = 0.0;
112  return;
113  }
114 
115 
116  meanfield = new G4QMDMeanField();
117  meanfield->SetSystem( this );
118 
119  //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons Begin ( z , a ) ( " << z << ", " << a << " )" << std::endl;
120  packNucleons();
121  //std::cout << "G4QMDGroundStateNucleus( G4int z , G4int a ) packNucleons End" << std::endl;
122 
123  delete meanfield;
124 
125 }
126 
127 
128 
129 void G4QMDGroundStateNucleus::packNucleons()
130 {
131 
132  //std::cout << "G4QMDGroundStateNucleus::packNucleons" << std::endl;
133 
134  ebini = - G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() ) / GetMassNumber();
135 
136  G4double ebin00 = ebini * 0.001;
137 
138  G4double ebin0 = 0.0;
139  G4double ebin1 = 0.0;
140 
141  if ( GetMassNumber() != 4 )
142  {
143  ebin0 = ( ebini - 0.5 ) * 0.001;
144  ebin1 = ( ebini + 0.5 ) * 0.001;
145  }
146  else
147  {
148  ebin0 = ( ebini - 1.5 ) * 0.001;
149  ebin1 = ( ebini + 1.5 ) * 0.001;
150  }
151 
152  G4int n0Try = 0;
153  G4bool isThisOK = false;
154  while ( n0Try < maxTrial ) // Loop checking, 11.03.2015, T. Koi
155  {
156  n0Try++;
157  //std::cout << "TKDB packNucleons n0Try " << n0Try << std::endl;
158 
159 // Sampling Position
160 
161  //std::cout << "TKDB Sampling Position " << std::endl;
162 
163  G4bool areThesePsOK = false;
164  G4int npTry = 0;
165  while ( npTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
166  {
167  //std::cout << "TKDB Sampling Position npTry " << npTry << std::endl;
168  npTry++;
169  G4int i = 0;
170  if ( samplingPosition( i ) )
171  {
172  //std::cout << "packNucleons samplingPosition 0 succeed " << std::endl;
173  for ( i = 1 ; i < GetMassNumber() ; i++ )
174  {
175  //std::cout << "packNucleons samplingPosition " << i << " trying " << std::endl;
176  if ( !( samplingPosition( i ) ) )
177  {
178  //std::cout << "packNucleons samplingPosition " << i << " failed" << std::endl;
179  break;
180  }
181  }
182  if ( i == GetMassNumber() )
183  {
184  //std::cout << "packNucleons samplingPosition all scucceed " << std::endl;
185  areThesePsOK = true;
186  break;
187  }
188  }
189  }
190  if ( areThesePsOK == false ) continue;
191 
192  //std::cout << "TKDB Sampling Position End" << std::endl;
193 
194 // Calculate Two-body quantities
195 
196  meanfield->Cal2BodyQuantities();
197  std::vector< G4double > rho_a ( GetMassNumber() , 0.0 );
198  std::vector< G4double > rho_s ( GetMassNumber() , 0.0 );
199  std::vector< G4double > rho_c ( GetMassNumber() , 0.0 );
200 
201  rho_l.resize ( GetMassNumber() , 0.0 );
202  d_pot.resize ( GetMassNumber() , 0.0 );
203 
204  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
205  {
206  for ( G4int j = 0 ; j < GetMassNumber() ; j++ )
207  {
208 
209  rho_a[ i ] += meanfield->GetRHA( j , i );
210  G4int k = 0;
211  if ( participants[j]->GetDefinition() != participants[i]->GetDefinition() )
212  {
213  k = 1;
214  }
215 
216  rho_s[ i ] += meanfield->GetRHA( j , i )*( 1.0 - 2.0 * k ); // OK?
217 
218  rho_c[ i ] += meanfield->GetRHE( j , i );
219  }
220 
221  }
222 
223  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
224  {
225  rho_l[i] = cdp * ( rho_a[i] + 1.0 );
226  d_pot[i] = c0p * rho_a[i]
227  + c3p * G4Pow::GetInstance()->powA ( rho_a[i] , gamm )
228  + csp * rho_s[i]
229  + clp * rho_c[i];
230 
231  //std::cout << "d_pot[i] " << i << " " << d_pot[i] << std::endl;
232  }
233 
234 
235 // Sampling Momentum
236 
237  //std::cout << "TKDB Sampling Momentum " << std::endl;
238 
239  phase_g.clear();
240  phase_g.resize( GetMassNumber() , 0.0 );
241 
242  //std::cout << "TKDB Sampling Momentum 1st " << std::endl;
243  G4bool isThis1stMOK = false;
244  G4int nmTry = 0;
245  while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
246  {
247  nmTry++;
248  G4int i = 0;
249  if ( samplingMomentum( i ) )
250  {
251  isThis1stMOK = true;
252  break;
253  }
254  }
255  if ( isThis1stMOK == false ) continue;
256 
257  //std::cout << "TKDB Sampling Momentum 2nd so on" << std::endl;
258 
259  G4bool areTheseMsOK = true;
260  nmTry = 0;
261  while ( nmTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
262  {
263  nmTry++;
264  G4int i = 0;
265  for ( i = 1 ; i < GetMassNumber() ; i++ )
266  {
267  //std::cout << "TKDB packNucleons samplingMomentum try " << i << std::endl;
268  if ( !( samplingMomentum( i ) ) )
269  {
270  //std::cout << "TKDB packNucleons samplingMomentum " << i << " failed" << std::endl;
271  areTheseMsOK = false;
272  break;
273  }
274  }
275  if ( i == GetMassNumber() )
276  {
277  areTheseMsOK = true;
278  }
279 
280  break;
281  }
282  if ( areTheseMsOK == false ) continue;
283 
284 // Kill Angluar Momentum
285 
286  //std::cout << "TKDB Sampling Kill Angluar Momentum " << std::endl;
287 
288  killCMMotionAndAngularM();
289 
290 
291 // Check Binding Energy
292 
293  //std::cout << "packNucleons Check Binding Energy Begin " << std::endl;
294 
295  G4double ekinal = 0.0;
296  for ( int i = 0 ; i < GetMassNumber() ; i++ )
297  {
298  ekinal += participants[i]->GetKineticEnergy();
299  }
300 
301  meanfield->Cal2BodyQuantities();
302 
303  G4double totalPotentialE = meanfield->GetTotalPotential();
304  G4double ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
305 
306  //std::cout << "packNucleons totalPotentialE " << totalPotentialE << std::endl;
307  //std::cout << "packNucleons ebinal " << ebinal << std::endl;
308  //std::cout << "packNucleons ekinal " << ekinal << std::endl;
309 
310  if ( ebinal < ebin0 || ebinal > ebin1 )
311  {
312  //std::cout << "packNucleons ebin0 " << ebin0 << std::endl;
313  //std::cout << "packNucleons ebin1 " << ebin1 << std::endl;
314  //std::cout << "packNucleons ebinal " << ebinal << std::endl;
315  //std::cout << "packNucleons Check Binding Energy Failed " << std::endl;
316  continue;
317  }
318 
319  //std::cout << "packNucleons Check Binding Energy End = OK " << std::endl;
320 
321 
322 // Energy Adujstment
323 
324  G4double dtc = 1.0;
325  G4double frg = -0.1;
326  G4double rdf0 = 0.5;
327 
328  G4double edif0 = ebinal - ebin00;
329 
330  G4double cfrc = 0.0;
331  if ( 0 < edif0 )
332  {
333  cfrc = frg;
334  }
335  else
336  {
337  cfrc = -frg;
338  }
339 
340  G4int ifrc = 1;
341 
342  G4int neaTry = 0;
343 
344  G4bool isThisEAOK = false;
345  while ( neaTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
346  {
347  neaTry++;
348 
349  G4double edif = ebinal - ebin00;
350 
351  //std::cout << "TKDB edif " << edif << std::endl;
352  if ( std::abs ( edif ) < epse )
353  {
354 
355  isThisEAOK = true;
356  //std::cout << "isThisEAOK " << isThisEAOK << std::endl;
357  break;
358  }
359 
360  G4int jfrc = 0;
361  if ( edif < 0.0 )
362  {
363  jfrc = 1;
364  }
365  else
366  {
367  jfrc = -1;
368  }
369 
370  if ( jfrc != ifrc )
371  {
372  cfrc = -rdf0 * cfrc;
373  dtc = rdf0 * dtc;
374  }
375 
376  if ( jfrc == ifrc && std::abs( edif0 ) < std::abs( edif ) )
377  {
378  cfrc = -rdf0 * cfrc;
379  dtc = rdf0 * dtc;
380  }
381 
382  ifrc = jfrc;
383  edif0 = edif;
384 
385  meanfield->CalGraduate();
386 
387  for ( int i = 0 ; i < GetMassNumber() ; i++ )
388  {
389  G4ThreeVector ri = participants[i]->GetPosition();
390  G4ThreeVector p_i = participants[i]->GetMomentum();
391 
392  ri += dtc * ( meanfield->GetFFr(i) - cfrc * ( meanfield->GetFFp(i) ) );
393  p_i += dtc * ( meanfield->GetFFp(i) + cfrc * ( meanfield->GetFFr(i) ) );
394 
395  participants[i]->SetPosition( ri );
396  participants[i]->SetMomentum( p_i );
397  }
398 
399  ekinal = 0.0;
400 
401  for ( int i = 0 ; i < GetMassNumber() ; i++ )
402  {
403  ekinal += participants[i]->GetKineticEnergy();
404  }
405 
406  meanfield->Cal2BodyQuantities();
407  totalPotentialE = meanfield->GetTotalPotential();
408 
409  ebinal = ( totalPotentialE + ekinal ) / double ( GetMassNumber() );
410 
411  }
412  //std::cout << "isThisEAOK " << isThisEAOK << std::endl;
413  if ( isThisEAOK == false ) continue;
414 
415  isThisOK = true;
416  //std::cout << "isThisOK " << isThisOK << std::endl;
417  break;
418 
419  }
420 
421  if ( isThisOK == false )
422  {
423  G4cout << "GroundStateNucleus state cannot be created. Try again with another parameters." << G4endl;
424  }
425 
426  //std::cout << "packNucleons End" << std::endl;
427  return;
428 }
429 
430 
431 G4bool G4QMDGroundStateNucleus::samplingPosition( G4int i )
432 {
433 
434  G4bool result = false;
435 
436  G4int nTry = 0;
437 
438  while ( nTry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
439  {
440 
441  //std::cout << "samplingPosition i th particle, nTtry " << i << " " << nTry << std::endl;
442 
443  G4double rwod = -1.0;
444  G4double rrr = 0.0;
445 
446  G4double rx = 0.0;
447  G4double ry = 0.0;
448  G4double rz = 0.0;
449 
450  G4int icounter = 0;
451  G4int icounter_max = 1024;
452  while ( G4UniformRand() * rmax > rwod ) // Loop checking, 11.03.2015, T. Koi
453  {
454  icounter++;
455  if ( icounter > icounter_max ) {
456  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
457  break;
458  }
459 
460  G4double rsqr = 10.0;
461  G4int jcounter = 0;
462  G4int jcounter_max = 1024;
463  while ( rsqr > 1.0 ) // Loop checking, 11.03.2015, T. Koi
464  {
465  jcounter++;
466  if ( jcounter > jcounter_max ) {
467  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
468  break;
469  }
470  rx = 1.0 - 2.0 * G4UniformRand();
471  ry = 1.0 - 2.0 * G4UniformRand();
472  rz = 1.0 - 2.0 * G4UniformRand();
473  rsqr = rx*rx + ry*ry + rz*rz;
474  }
475  rrr = radm * std::sqrt ( rsqr );
476  rwod = 1.0 / ( 1.0 + G4Exp ( ( rrr - rt00 ) / saa ) );
477 
478  }
479 
480  participants[i]->SetPosition( G4ThreeVector( rx , ry , rz )*radm );
481 
482  if ( i == 0 )
483  {
484  result = true;
485  return result;
486  }
487 
488 // i > 1 ( Second Particle or later )
489 // Check Distance to others
490 
491  G4bool isThisOK = true;
492  for ( G4int j = 0 ; j < i ; j++ )
493  {
494 
495  G4double r2 = participants[j]->GetPosition().diff2( participants[i]->GetPosition() );
496  G4double dmin2 = 0.0;
497 
498  if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
499  {
500  dmin2 = dsam2;
501  }
502  else
503  {
504  dmin2 = ddif2;
505  }
506 
507  //std::cout << "distance between j and i " << j << " " << i << " " << r2 << " " << dmin2 << std::endl;
508  if ( r2 < dmin2 )
509  {
510  isThisOK = false;
511  break;
512  }
513 
514  }
515 
516  if ( isThisOK == true )
517  {
518  result = true;
519  return result;
520  }
521 
522  nTry++;
523 
524  }
525 
526 // Here return "false"
527  return result;
528 }
529 
530 
531 
532 G4bool G4QMDGroundStateNucleus::samplingMomentum( G4int i )
533 {
534 
535  //std::cout << "TKDB samplingMomentum for " << i << std::endl;
536 
537  G4bool result = false;
538 
539  G4double pfm = hbc * G4Pow::GetInstance()->A13 ( ( 3.0 / 2.0 * pi*pi * rho_l[i] ) );
540 
541  if ( 10 < GetMassNumber() && -5.5 < ebini )
542  {
543  pfm = pfm * ( 1.0 + 0.2 * std::sqrt( std::abs( 8.0 + ebini ) / 8.0 ) );
544  }
545 
546  //std::cout << "TKDB samplingMomentum pfm " << pfm << std::endl;
547 
548  std::vector< G4double > phase;
549  phase.resize( i+1 ); // i start from 0
550 
551  G4int ntry = 0;
552 // 710
553  while ( ntry < maxTrial ) // Loop checking, 11.03.2015, T. Koi
554  {
555 
556  //std::cout << " TKDB ntry " << ntry << std::endl;
557  ntry++;
558 
559  G4double ke = DBL_MAX;
560 
561  G4int tkdb_i =0;
562 // 700
563  G4int icounter = 0;
564  G4int icounter_max = 1024;
565  while ( ke + d_pot [i] > edepth ) // Loop checking, 11.03.2015, T. Koi
566  {
567  icounter++;
568  if ( icounter > icounter_max ) {
569  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
570  break;
571  }
572 
573  G4double psqr = 10.0;
574  G4double px = 0.0;
575  G4double py = 0.0;
576  G4double pz = 0.0;
577 
578  G4int jcounter = 0;
579  G4int jcounter_max = 1024;
580  while ( psqr > 1.0 ) // Loop checking, 11.03.2015, T. Koi
581  {
582  jcounter++;
583  if ( jcounter > jcounter_max ) {
584  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
585  break;
586  }
587  px = 1.0 - 2.0*G4UniformRand();
588  py = 1.0 - 2.0*G4UniformRand();
589  pz = 1.0 - 2.0*G4UniformRand();
590 
591  psqr = px*px + py*py + pz*pz;
592  }
593 
594  G4ThreeVector p ( px , py , pz );
595  p = pfm * p;
596  participants[i]->SetMomentum( p );
597  G4LorentzVector p4 = participants[i]->Get4Momentum();
598  //ke = p4.e() - p4.restMass();
599  ke = participants[i]->GetKineticEnergy();
600 
601 
602  tkdb_i++;
603  if ( tkdb_i > maxTrial ) return result; // return false
604 
605  }
606 
607  //std::cout << "TKDB ke d_pot[i] " << ke << " " << d_pot[i] << std::endl;
608 
609 
610  if ( i == 0 )
611  {
612  result = true;
613  return result;
614  }
615 
616  G4bool isThisOK = true;
617 
618  // Check Pauli principle
619 
620  phase[ i ] = 0.0;
621 
622  //std::cout << "TKDB Check Pauli Principle " << i << std::endl;
623 
624  for ( G4int j = 0 ; j < i ; j++ )
625  {
626  phase[ j ] = 0.0;
627  //std::cout << "TKDB Check Pauli Principle i , j " << i << " , " << j << std::endl;
628  G4double expa = 0.0;
629  if ( participants[j]->GetDefinition() == participants[i]->GetDefinition() )
630  {
631 
632  expa = - meanfield->GetRR2(i,j) * cpw;
633 
634  if ( expa > epsx )
635  {
636  G4ThreeVector p_i = participants[i]->GetMomentum();
637  G4ThreeVector pj = participants[j]->GetMomentum();
638  G4double dist2_p = p_i.diff2( pj );
639 
640  dist2_p = dist2_p*cph;
641  expa = expa - dist2_p;
642 
643  if ( expa > epsx )
644  {
645 
646  phase[j] = G4Exp ( expa );
647 
648  if ( phase[j] * cpc > 0.2 )
649  {
650 /*
651  G4cout << "TKDB Check Pauli Principle A i , j " << i << " , " << j << G4endl;
652  G4cout << "TKDB Check Pauli Principle phase[j] " << phase[j] << G4endl;
653  G4cout << "TKDB Check Pauli Principle phase[j]*cpc > 0.2 " << phase[j]*cpc << G4endl;
654 */
655  isThisOK = false;
656  break;
657  }
658  if ( ( phase_g[j] + phase[j] ) * cpc > 0.5 )
659  {
660 /*
661  G4cout << "TKDB Check Pauli Principle B i , j " << i << " , " << j << G4endl;
662  G4cout << "TKDB Check Pauli Principle B phase_g[j] " << phase_g[j] << G4endl;
663  G4cout << "TKDB Check Pauli Principle B phase[j] " << phase[j] << G4endl;
664  G4cout << "TKDB Check Pauli Principle B phase_g[j] + phase[j] ) * cpc > 0.5 " << ( phase_g[j] + phase[j] ) * cpc << G4endl;
665 */
666  isThisOK = false;
667  break;
668  }
669 
670  phase[i] += phase[j];
671  if ( phase[i] * cpc > 0.3 )
672  {
673 /*
674  G4cout << "TKDB Check Pauli Principle C i , j " << i << " , " << j << G4endl;
675  G4cout << "TKDB Check Pauli Principle C phase[i] " << phase[i] << G4endl;
676  G4cout << "TKDB Check Pauli Principle C phase[i] * cpc > 0.3 " << phase[i] * cpc << G4endl;
677 */
678  isThisOK = false;
679  break;
680  }
681 
682  //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
683 
684  }
685  else
686  {
687  //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
688  }
689 
690  }
691  else
692  {
693  //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
694  }
695 
696  }
697  else
698  {
699  //std::cout << "TKDB Check Pauli Principle OK i , j " << i << " , " << j << std::endl;
700  }
701 
702  }
703 
704  if ( isThisOK == true )
705  {
706 
707  phase_g[i] = phase[i];
708 
709  for ( int j = 0 ; j < i ; j++ )
710  {
711  phase_g[j] += phase[j];
712  }
713 
714  result = true;
715  return result;
716  }
717 
718  }
719 
720  return result;
721 
722 }
723 
724 
725 
726 void G4QMDGroundStateNucleus::killCMMotionAndAngularM()
727 {
728 
729 // CalEnergyAndAngularMomentumInCM();
730 
731  //std::vector< G4ThreeVector > p ( GetMassNumber() , 0.0 );
732  //std::vector< G4ThreeVector > r ( GetMassNumber() , 0.0 );
733 
734 // Move to cm system
735 
736  G4ThreeVector pcm_tmp ( 0.0 );
737  G4ThreeVector rcm_tmp ( 0.0 );
738  G4double sumMass = 0.0;
739 
740  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
741  {
742  pcm_tmp += participants[i]->GetMomentum();
743  rcm_tmp += participants[i]->GetPosition() * participants[i]->GetMass();
744  sumMass += participants[i]->GetMass();
745  }
746 
747  pcm_tmp = pcm_tmp/GetMassNumber();
748  rcm_tmp = rcm_tmp/sumMass;
749 
750  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
751  {
752  participants[i]->SetMomentum( participants[i]->GetMomentum() - pcm_tmp );
753  participants[i]->SetPosition( participants[i]->GetPosition() - rcm_tmp );
754  }
755 
756 // kill the angular momentum
757 
758  G4ThreeVector ll ( 0.0 );
759  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
760  {
761  ll += participants[i]->GetPosition().cross ( participants[i]->GetMomentum() );
762  }
763 
764  G4double rr[3][3];
765  G4double ss[3][3];
766  for ( G4int i = 0 ; i < 3 ; i++ )
767  {
768  for ( G4int j = 0 ; j < 3 ; j++ )
769  {
770  rr [i][j] = 0.0;
771 
772  if ( i == j )
773  {
774  ss [i][j] = 1.0;
775  }
776  else
777  {
778  ss [i][j] = 0.0;
779  }
780  }
781  }
782 
783  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
784  {
785  G4ThreeVector r = participants[i]->GetPosition();
786  rr[0][0] += r.y() * r.y() + r.z() * r.z();
787  rr[1][0] += - r.y() * r.x();
788  rr[2][0] += - r.z() * r.x();
789  rr[0][1] += - r.x() * r.y();
790  rr[1][1] += r.z() * r.z() + r.x() * r.x();
791  rr[2][1] += - r.z() * r.y();
792  rr[2][0] += - r.x() * r.z();
793  rr[2][1] += - r.y() * r.z();
794  rr[2][2] += r.x() * r.x() + r.y() * r.y();
795  }
796 
797  for ( G4int i = 0 ; i < 3 ; i++ )
798  {
799  G4double x = rr [i][i];
800  for ( G4int j = 0 ; j < 3 ; j++ )
801  {
802  rr[i][j] = rr[i][j] / x;
803  ss[i][j] = ss[i][j] / x;
804  }
805  for ( G4int j = 0 ; j < 3 ; j++ )
806  {
807  if ( j != i )
808  {
809  G4double y = rr [j][i];
810  for ( G4int k = 0 ; k < 3 ; k++ )
811  {
812  rr[j][k] += -y * rr[i][k];
813  ss[j][k] += -y * ss[i][k];
814  }
815  }
816  }
817  }
818 
819  G4double opl[3];
820  G4double rll[3];
821 
822  rll[0] = ll.x();
823  rll[1] = ll.y();
824  rll[2] = ll.z();
825 
826  for ( G4int i = 0 ; i < 3 ; i++ )
827  {
828  opl[i] = 0.0;
829 
830  for ( G4int j = 0 ; j < 3 ; j++ )
831  {
832  opl[i] += ss[i][j]*rll[j];
833  }
834  }
835 
836  for ( G4int i = 0 ; i < GetMassNumber() ; i++ )
837  {
838  G4ThreeVector p_i = participants[i]->GetMomentum() ;
839  G4ThreeVector ri = participants[i]->GetPosition() ;
840  G4ThreeVector opl_v ( opl[0] , opl[1] , opl[2] );
841 
842  p_i += ri.cross(opl_v);
843 
844  participants[i]->SetMomentum( p_i );
845  }
846 
847 }