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

Last change on this file since 1055 was 962, checked in by garnier, 15 years ago

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