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

Last change on this file since 1253 was 1196, checked in by garnier, 15 years ago

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