source: trunk/source/processes/hadronic/models/qmd/src/G4QMDGroundStateNucleus.cc @ 869

Last change on this file since 869 was 819, checked in by garnier, 16 years ago

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