Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QMDGroundStateNucleus.cc
Go to the documentation of this file.
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 //
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 
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 
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 
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
106  ebini = 0.0;
107  return;
108  }
109  else if ( z == 0 && a == 1 ) { // Neutron primary
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 
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 }
G4double G4ParticleHPJENDLHEData::G4double result
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double GetRHA(G4int i, G4int j)
G4int GetAtomicNumber()
Definition: G4QMDNucleus.cc:89
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4int GetMassNumber()
Definition: G4QMDNucleus.cc:68
void SetParticipant(G4QMDParticipant *particle)
Definition: G4QMDSystem.hh:51
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetFFr(G4int i)
double x() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double Get_hbc()
const char * p
Definition: xmltok.h:285
G4double Get_c0p()
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
void Cal2BodyQuantities()
double z() const
G4double Get_cpc()
double diff2(const Hep3Vector &v) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static G4QMDParameters * GetInstance()
bool G4bool
Definition: G4Types.hh:79
G4double Get_csp()
G4double Get_cdp()
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double Get_epsx()
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double Get_gamm()
G4double GetTotalPotential()
std::vector< G4QMDParticipant * > participants
Definition: G4QMDSystem.hh:72
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetRHE(G4int i, G4int j)
G4double A13(G4double A) const
Definition: G4Pow.hh:132
static G4double GetBindingEnergy(const G4int A, const G4int Z)
G4double GetRR2(G4int i, G4int j)
G4double Get_clp()
void SetSystem(G4QMDSystem *aSystem)
double y() const
G4double Get_c3p()
tuple z
Definition: test.py:28
G4ThreeVector GetFFp(G4int i)
#define G4endl
Definition: G4ios.hh:61
G4double Get_cpw()
static constexpr double pi
Definition: G4SIunits.hh:75
Hep3Vector cross(const Hep3Vector &) const
G4double Get_cph()
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83