Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/qmd/src/G4QMDGroundStateNucleus.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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