Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QMDCollision.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // 080602 Fix memory leaks by T. Koi
27 // 081120 Add deltaT in signature of CalKinematicsOfBinaryCollisions
28 // Add several required updating of Mean Filed
29 // Modified handling of absorption case by T. Koi
30 // 090126 Fix in absorption case by T. Koi
31 // 090331 Fix for gamma participant by T. Koi
32 //
33 #include "G4QMDCollision.hh"
34 #include "G4Scatterer.hh"
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "Randomize.hh"
38 
39 G4QMDCollision::G4QMDCollision()
40 : deltar ( 4 )
41 , bcmax0 ( 1.323142 ) // NN maximum impact parameter
42 , bcmax1 ( 2.523 ) // others maximum impact parameter
43 // , sig0 ( 55 ) // NN cross section
44 //110617 fix for gcc 4.6 compilation warnings
45 //, sig1 ( 200 ) // others cross section
46 , epse ( 0.0001 )
47 {
48  theScatterer = new G4Scatterer();
49 }
50 
51 
52 
53 G4QMDCollision::~G4QMDCollision()
54 {
55  delete theScatterer;
56 }
57 
58 
59 void G4QMDCollision::CalKinematicsOfBinaryCollisions( G4double dt )
60 {
61  G4double deltaT = dt;
62 
63  G4int n = theSystem->GetTotalNumberOfParticipant();
64 //081118
65  //G4int nb = 0;
66  for ( G4int i = 0 ; i < n ; i++ )
67  {
68  theSystem->GetParticipant( i )->UnsetHitMark();
69  theSystem->GetParticipant( i )->UnsetHitMark();
70  //nb += theSystem->GetParticipant( i )->GetBaryonNumber();
71  }
72  //G4cout << "nb = " << nb << " n = " << n << G4endl;
73 
74 
75 //071101
76  for ( G4int i = 0 ; i < n ; i++ )
77  {
78 
79  //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl;
80 
81  if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() )
82  {
83 
84  G4bool decayed = false;
85 
86  G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition();
87  G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum();
88  G4ThreeVector r0 = theSystem->GetParticipant( i )->GetPosition();
89 
90  G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum();
91 
92  G4double eini = theMeanField->GetTotalPotential() + p40.e();
93 
94  G4int n0 = theSystem->GetTotalNumberOfParticipant();
95  G4int i0 = 0;
96 
97  G4bool isThisEnergyOK = false;
98 
99  for ( G4int ii = 0 ; ii < 4 ; ii++ )
100  {
101 
102  //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum();
103  G4LorentzVector p400 = p40;
104 
105  p400 *= GeV;
106  //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 );
107  G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 );
108  //std::cout << "G4KineticTrack " << i << " " << kt.GetDefinition()->GetParticleName() << kt.GetPosition() << std::endl;
109  G4KineticTrackVector* secs = NULL;
110  secs = kt.Decay();
111  G4int id = 0;
112  G4double et = 0;
113  if ( secs )
114  {
115  for ( G4KineticTrackVector::iterator it
116  = secs->begin() ; it != secs->end() ; it++ )
117  {
118 /*
119  G4cout << "G4KineticTrack"
120  << " " << (*it)->GetDefinition()->GetParticleName()
121  << " " << (*it)->Get4Momentum()
122  << " " << (*it)->GetPosition()/fermi
123  << G4endl;
124 */
125  if ( id == 0 )
126  {
127  theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
128  theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV );
129  theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi );
130  //theMeanField->Cal2BodyQuantities( i );
131  et += (*it)->Get4Momentum().e()/GeV;
132  }
133  if ( id > 0 )
134  {
135  // Append end;
136  theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) );
137  et += (*it)->Get4Momentum().e()/GeV;
138  if ( id > 1 )
139  {
140  //081118
141  //G4cout << "G4QMDCollision id >2; id= " << id << G4endl;
142  }
143  }
144  id++; // number of daughter particles
145 
146  delete *it;
147  }
148 
149  theMeanField->Update();
150  i0 = id-1; // 0 enter to i
151 
152  delete secs;
153  }
154 
155 // EnergyCheck
156 
157  G4double efin = theMeanField->GetTotalPotential() + et;
158  //std::cout << std::abs ( eini - efin ) - epse << std::endl;
159 // std::cout << std::abs ( eini - efin ) - epse*10 << std::endl;
160 // *10 TK
161  if ( std::abs ( eini - efin ) < epse*10 )
162  {
163  // Energy OK
164  isThisEnergyOK = true;
165  break;
166  }
167  else
168  {
169 
170  theSystem->GetParticipant( i )->SetDefinition( pd0 );
171  theSystem->GetParticipant( i )->SetPosition( r0 );
172  theSystem->GetParticipant( i )->SetMomentum( p0 );
173 
174  for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ )
175  {
176  //081118
177  //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0 << std::endl;
178  theSystem->DeleteParticipant( i0i+n0 );
179  }
180  //081103
181  theMeanField->Update();
182  }
183 
184  }
185 
186 
187 // Pauli Check
188  if ( isThisEnergyOK == true )
189  {
190  if ( theMeanField->IsPauliBlocked ( i ) != true )
191  {
192 
193  G4bool allOK = true;
194  for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
195  {
196  if ( theMeanField->IsPauliBlocked ( i0i+n0 ) == true )
197  {
198  allOK = false;
199  break;
200  }
201  }
202 
203  if ( allOK )
204  {
205  decayed = true; //Decay Succeeded
206  }
207  }
208 
209  }
210 //
211 
212  if ( decayed )
213  {
214  //081119
215  //G4cout << "Decay Suceeded! " << std::endl;
216  theSystem->GetParticipant( i )->SetHitMark();
217  for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
218  {
219  theSystem->GetParticipant( i0i+n0 )->SetHitMark();
220  }
221 
222  }
223  else
224  {
225 
226 // Decay Blocked and re-enter orginal participant;
227 
228  if ( isThisEnergyOK == true ) // for false case already done
229  {
230 
231  theSystem->GetParticipant( i )->SetDefinition( pd0 );
232  theSystem->GetParticipant( i )->SetPosition( r0 );
233  theSystem->GetParticipant( i )->SetMomentum( p0 );
234 
235  for ( G4int i0i = 0 ; i0i < i0 ; i0i++ )
236  {
237  //081118
238  //std::cout << "Decay Blocked deleteing " << i0i+n0 << std::endl;
239  theSystem->DeleteParticipant( i0i+n0 );
240  }
241  //081103
242  theMeanField->Update();
243  }
244 
245  }
246 
247  } //shortlive
248  } // go next participant
249 //071101
250 
251 
252  n = theSystem->GetTotalNumberOfParticipant();
253 
254 //081118
255  //for ( G4int i = 1 ; i < n ; i++ )
256  for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ )
257  {
258 
259  //std::cout << "Collision i " << i << std::endl;
260  if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
261 
262  G4ThreeVector ri = theSystem->GetParticipant( i )->GetPosition();
263  G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
264  G4double rmi = theSystem->GetParticipant( i )->GetMass();
265  G4ParticleDefinition* pdi = theSystem->GetParticipant( i )->GetDefinition();
266 //090331 gamma
267  if ( pdi->GetPDGMass() == 0.0 ) continue;
268 
269  //std::cout << " p4i00 " << p4i << std::endl;
270  for ( G4int j = 0 ; j < i ; j++ )
271  {
272 
273 
274 /*
275  std::cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisProjectile() << std::endl;
276  std::cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisProjectile() << std::endl;
277  std::cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisTarget() << std::endl;
278  std::cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisTarget() << std::endl;
279 */
280 
281  // Only 1 Collision allowed for each particle in a time step.
282  //081119
283  if ( theSystem->GetParticipant( i )->IsThisHit() ) continue;
284  if ( theSystem->GetParticipant( j )->IsThisHit() ) continue;
285 
286  //std::cout << "Collision " << i << " " << j << std::endl;
287 
288  // Do not allow collision between nucleons in target/projectile til its first collision.
289  if ( theSystem->GetParticipant( i )->IsThisProjectile() )
290  {
291  if ( theSystem->GetParticipant( j )->IsThisProjectile() ) continue;
292  }
293  else if ( theSystem->GetParticipant( i )->IsThisTarget() )
294  {
295  if ( theSystem->GetParticipant( j )->IsThisTarget() ) continue;
296  }
297 
298 
299  G4ThreeVector rj = theSystem->GetParticipant( j )->GetPosition();
300  G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
301  G4double rmj = theSystem->GetParticipant( j )->GetMass();
302  G4ParticleDefinition* pdj = theSystem->GetParticipant( j )->GetDefinition();
303 //090331 gamma
304  if ( pdj->GetPDGMass() == 0.0 ) continue;
305 
306  G4double rr2 = theMeanField->GetRR2( i , j );
307 
308 // Here we assume elab (beam momentum less than 5 GeV/n )
309  if ( rr2 > deltar*deltar ) continue;
310 
311  //G4double s = (p4i+p4j)*(p4i+p4j);
312  //G4double srt = std::sqrt ( s );
313 
314  G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) );
315 
316  G4double cutoff = 0.0;
317  G4double bcmax = 0.0;
318  //110617 fix for gcc 4.6 compilation warnings
319  //G4double sig = 0.0;
320 
321  if ( rmi < 0.94 && rmj < 0.94 )
322  {
323 // nucleon or pion case
324  cutoff = rmi + rmj + 0.02;
325  bcmax = bcmax0;
326  //110617 fix for gcc 4.6 compilation warnings
327  //sig = sig0;
328  }
329  else
330  {
331  cutoff = rmi + rmj;
332  bcmax = bcmax1;
333  //110617 fix for gcc compilation warnings
334  //sig = sig1;
335  }
336 
337  //std::cout << "Collision cutoff " << i << " " << j << " " << cutoff << std::endl;
338  if ( srt < cutoff ) continue;
339 
340  G4ThreeVector dr = ri - rj;
341  G4double rsq = dr*dr;
342 
343  G4double pij = p4i*p4j;
344  G4double pidr = p4i.vect()*dr;
345  G4double pjdr = p4j.vect()*dr;
346 
347  G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij );
348  G4double bij = pidr / rmi - pjdr*rmi/pij;
349  G4double cij = rsq + ( pidr / rmi ) * ( pidr / rmi );
350  G4double brel = std::sqrt ( std::abs ( cij - bij*bij/aij ) );
351 
352  if ( brel > bcmax ) continue;
353  //std::cout << "collisions3 " << std::endl;
354 
355  G4double bji = -pjdr/rmj + pidr * rmj /pij;
356 
357  G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi;
358  G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj;
359 
360 
361 /*
362  std::cout << "collisions4 p4i " << p4i << std::endl;
363  std::cout << "collisions4 ri " << ri << std::endl;
364  std::cout << "collisions4 p4j " << p4j << std::endl;
365  std::cout << "collisions4 rj " << rj << std::endl;
366  std::cout << "collisions4 dr " << dr << std::endl;
367  std::cout << "collisions4 pij " << pij << std::endl;
368  std::cout << "collisions4 aij " << aij << std::endl;
369  std::cout << "collisions4 bij bji " << bij << " " << bji << std::endl;
370  std::cout << "collisions4 pidr pjdr " << pidr << " " << pjdr << std::endl;
371  std::cout << "collisions4 p4i.e() p4j.e() " << p4i.e() << " " << p4j.e() << std::endl;
372  std::cout << "collisions4 rmi rmj " << rmi << " " << rmj << std::endl;
373  std::cout << "collisions4 " << ti << " " << tj << std::endl;
374 */
375  if ( std::abs ( ti + tj ) > deltaT ) continue;
376  //std::cout << "collisions4 " << std::endl;
377 
378  G4ThreeVector beta = ( p4i + p4j ).boostVector();
379 
380  G4LorentzVector p = p4i;
381  G4LorentzVector p4icm = p.boost( p.findBoostToCM ( p4j ) );
382  G4ThreeVector pcm = p4icm.vect();
383 
384  G4double prcm = pcm.mag();
385 
386  if ( prcm <= 0.00001 ) continue;
387  //std::cout << "collisions5 " << std::endl;
388 
389  G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollision ( i , j ) ); // Use Geant4 Collision Library
390  //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic
391 
392 /*
393  G4bool pauli_blocked = false;
394  if ( energetically_forbidden == false ) // result true
395  {
396  if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true )
397  {
398  pauli_blocked = true;
399  //std::cout << "G4QMDRESULT Collsion Pauli Blocked " << std::endl;
400  }
401  }
402  else
403  {
404  if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true )
405  pauli_blocked = false;
406  //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl;
407  }
408 */
409 
410 /*
411  std::cout << "G4QMDRESULT Collsion initial p4 i and j "
412  << p4i << " " << p4j
413  << std::endl;
414 */
415 // 081118
416  //if ( energetically_forbidden == true || pauli_blocked == true )
417  if ( energetically_forbidden == true )
418  {
419 
420  //G4cout << " energetically_forbidden " << G4endl;
421 // Collsion not allowed then re enter orginal participants
422 // Now only momentum, becasuse we only consider elastic scattering of nucleons
423 
424  theSystem->GetParticipant( i )->SetMomentum( p4i.vect() );
425  theSystem->GetParticipant( i )->SetDefinition( pdi );
426  theSystem->GetParticipant( i )->SetPosition( ri );
427 
428  theSystem->GetParticipant( j )->SetMomentum( p4j.vect() );
429  theSystem->GetParticipant( j )->SetDefinition( pdj );
430  theSystem->GetParticipant( j )->SetPosition( rj );
431 
432  theMeanField->Cal2BodyQuantities( i );
433  theMeanField->Cal2BodyQuantities( j );
434 
435  }
436  else
437  {
438 
439 
440  G4bool absorption = false;
441  if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true;
442  if ( absorption )
443  {
444  //G4cout << "Absorption happend " << G4endl;
445  i = i-1;
446  n = n-1;
447  }
448 
449 // Collsion allowed (really happened)
450 
451  // Unset Projectile/Target flag
452  theSystem->GetParticipant( i )->UnsetInitialMark();
453  if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark();
454 
455  theSystem->GetParticipant( i )->SetHitMark();
456  if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark();
457 
458  theSystem->IncrementCollisionCounter();
459 
460 /*
461  std::cout << "G4QMDRESULT Collsion Really Happened between "
462  << i << " and " << j
463  << std::endl;
464  std::cout << "G4QMDRESULT Collsion initial p4 i and j "
465  << p4i << " " << p4j
466  << std::endl;
467  std::cout << "G4QMDRESULT Collsion after p4 i and j "
468  << theSystem->GetParticipant( i )->Get4Momentum()
469  << " "
470  << theSystem->GetParticipant( j )->Get4Momentum()
471  << std::endl;
472  std::cout << "G4QMDRESULT Collsion Diff "
473  << p4i + p4j - theSystem->GetParticipant( i )->Get4Momentum() - theSystem->GetParticipant( j )->Get4Momentum()
474  << std::endl;
475  std::cout << "G4QMDRESULT Collsion initial r i and j "
476  << ri << " " << rj
477  << std::endl;
478  std::cout << "G4QMDRESULT Collsion after r i and j "
479  << theSystem->GetParticipant( i )->GetPosition()
480  << " "
481  << theSystem->GetParticipant( j )->GetPosition()
482  << std::endl;
483 */
484 
485 
486  }
487 
488  }
489 
490  }
491 
492 
493 }
494 
495 
496 
497 G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollision( G4int i , G4int j )
498 {
499 
500 //081103
501  //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl;
502 
503  G4bool result = false;
504  G4bool energyOK = false;
505  G4bool pauliOK = false;
506  G4bool abs = false;
507  G4QMDParticipant* absorbed = NULL;
508 
509  G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
510  G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
511 
512 //071031
513 
514  G4double epot = theMeanField->GetTotalPotential();
515 
516  G4double eini = epot + p4i.e() + p4j.e();
517 
518 //071031
519  // will use KineticTrack
520  G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition();
521  G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition();
522  G4LorentzVector p4i0 = p4i*GeV;
523  G4LorentzVector p4j0 = p4j*GeV;
524  G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi;
525  G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi;
526 
527  for ( G4int iitry = 0 ; iitry < 4 ; iitry++ )
528  {
529 
530  abs = false;
531 
532  G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 );
533  G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 );
534 
535  G4LorentzVector p4ix_new;
536  G4LorentzVector p4jx_new;
537  G4KineticTrackVector* secs = NULL;
538  secs = theScatterer->Scatter( kt1 , kt2 );
539 
540  //std::cout << "G4QMDSCATTERER BEFORE " << kt1.GetDefinition()->GetParticleName() << " " << kt1.Get4Momentum()/GeV << " " << kt1.GetPosition()/fermi << std::endl;
541  //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl;
542  //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl;
543 
544 
545  if ( secs )
546  {
547  G4int iti = 0;
548  if ( secs->size() == 2 )
549  {
550  for ( G4KineticTrackVector::iterator it
551  = secs->begin() ; it != secs->end() ; it++ )
552  {
553  if ( iti == 0 )
554  {
555  theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() );
556  p4ix_new = (*it)->Get4Momentum()/GeV;
557  //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl;
558  theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
559  }
560  if ( iti == 1 )
561  {
562  theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() );
563  p4jx_new = (*it)->Get4Momentum()/GeV;
564  //std::cout << "THESCATTERER " << p4jx_new.e()-p4jx_new.m() << std::endl;
565  theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() );
566  }
567  //std::cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << std::endl;
568  iti++;
569  }
570  }
571  else if ( secs->size() == 1 )
572  {
573 //081118
574  abs = true;
575  //G4cout << "G4QMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl;
576  //secs->front()->Decay();
577  theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() );
578  p4ix_new = secs->front()->Get4Momentum()/GeV;
579  theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() );
580 
581  }
582 
583 //081118
584  if ( secs->size() > 2 )
585  {
586 
587  G4cout << "G4QMDCollision secs size > 2; " << secs->size() << G4endl;
588 
589  for ( G4KineticTrackVector::iterator it
590  = secs->begin() ; it != secs->end() ; it++ )
591  {
592  G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl;
593  }
594 
595  }
596 
597  // deleteing KineticTrack
598  for ( G4KineticTrackVector::iterator it
599  = secs->begin() ; it != secs->end() ; it++ )
600  {
601  delete *it;
602  }
603 
604  delete secs;
605  }
606 //071031
607 
608  if ( !abs )
609  {
610  theMeanField->Cal2BodyQuantities( i );
611  theMeanField->Cal2BodyQuantities( j );
612  }
613  else
614  {
615  absorbed = theSystem->EraseParticipant( j );
616  theMeanField->Update();
617  }
618 
619  epot = theMeanField->GetTotalPotential();
620 
621  G4double efin = epot + p4ix_new.e() + p4jx_new.e();
622 
623  //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - epse << std::endl;
624 
625 /*
626  std::cout << "Collision efin " << i << " " << j << " " << efin << std::endl;
627  std::cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << epse << std::endl;
628  std::cout << "Collision " << std::abs ( eini - efin ) << " " << epse << std::endl;
629 */
630 
631 //071031
632  if ( std::abs ( eini - efin ) < epse )
633  {
634  // Collison OK
635  //std::cout << "collisions6" << std::endl;
636  //std::cout << "collisions before " << p4i << " " << p4j << std::endl;
637  //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
638  //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
639  //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl;
640  //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
641  energyOK = true;
642  break;
643  }
644  else
645  {
646  //G4cout << "Energy Not OK " << G4endl;
647  if ( abs )
648  {
649  //G4cout << "TKDB reinsert j " << G4endl;
650  theSystem->InsertParticipant( absorbed , j );
651  theMeanField->Update();
652  }
653  // do not need reinsert in no absroption case
654  }
655 //071031
656  }
657 
658 // Energetically forbidden collision
659 
660  if ( energyOK )
661  {
662  // Pauli Check
663  //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl;
664  if ( !abs )
665  {
666  if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) )
667  {
668  //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
669  pauliOK = true;
670  }
671  }
672  else
673  {
674  //if ( theMeanField->IsPauliBlocked ( i ) == false )
675  //090126 i-1 cause jth is erased
676  if ( theMeanField->IsPauliBlocked ( i-1 ) == false )
677  {
678  //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl;
679  delete absorbed;
680  pauliOK = true;
681  }
682  }
683 
684 
685  if ( pauliOK )
686  {
687  result = true;
688  }
689  else
690  {
691  //G4cout << "Pauli Blocked" << G4endl;
692  if ( abs )
693  {
694  //G4cout << "TKDB reinsert j pauli block" << G4endl;
695  theSystem->InsertParticipant( absorbed , j );
696  theMeanField->Update();
697  }
698  }
699  }
700 
701  return result;
702 
703 }
704 
705 
706 
707 G4bool G4QMDCollision::CalFinalStateOfTheBinaryCollisionJQMD( G4double sig , G4double cutoff , G4ThreeVector pcm , G4double prcm , G4double srt , G4ThreeVector beta , G4double gamma , G4int i , G4int j )
708 {
709 
710  //G4cout << "CalFinalStateOfTheBinaryCollisionJQMD" << G4endl;
711 
712  G4bool result = true;
713 
714  G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum();
715  G4double rmi = theSystem->GetParticipant( i )->GetMass();
716  G4int zi = theSystem->GetParticipant( i )->GetChargeInUnitOfEplus();
717 
718  G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum();
719  G4double rmj = theSystem->GetParticipant( j )->GetMass();
720  G4int zj = theSystem->GetParticipant( j )->GetChargeInUnitOfEplus();
721 
722  G4double pr = prcm;
723 
724  G4double c2 = pcm.z()/pr;
725 
726  G4double csrt = srt - cutoff;
727 
728  //G4double pri = prcm;
729  //G4double prf = sqrt ( 0.25 * srt*srt -rm2 );
730 
731  G4double asrt = srt - rmi - rmj;
732  G4double pra = prcm;
733 
734 
735 
736  G4double elastic = 0.0;
737 
738  if ( zi == zj )
739  {
740  if ( csrt < 0.4286 )
741  {
742  elastic = 35.0 / ( 1. + csrt * 100.0 ) + 20.0;
743  }
744  else
745  {
746  elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
747  * 2. / pi + 1.0 ) * 9.65 + 7.0;
748  }
749  }
750  else
751  {
752  if ( csrt < 0.4286 )
753  {
754  elastic = 28.0 / ( 1. + csrt * 100.0 ) + 27.0;
755  }
756  else
757  {
758  elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 )
759  * 2. / pi + 1.0 ) * 12.34 + 10.0;
760  }
761  }
762 
763 // std::cout << "Collision csrt " << i << " " << j << " " << csrt << std::endl;
764 // std::cout << "Collision elstic " << i << " " << j << " " << elastic << std::endl;
765 
766 
767 // std::cout << "Collision sig " << i << " " << j << " " << sig << std::endl;
768  if ( G4UniformRand() > elastic / sig )
769  {
770  //std::cout << "Inelastic " << std::endl;
771  //std::cout << "elastic/sig " << elastic/sig << std::endl;
772  return result;
773  }
774  else
775  {
776  //std::cout << "Elastic " << std::endl;
777  }
778 // std::cout << "Collision ELSTIC " << i << " " << j << std::endl;
779 
780 
781  G4double as = std::pow ( 3.65 * asrt , 6 );
782  G4double a = 6.0 * as / (1.0 + as);
783  G4double ta = -2.0 * pra*pra;
784  G4double x = G4UniformRand();
785  G4double t1 = std::log( (1-x) * std::exp(2.*a*ta) + x ) / a;
786  G4double c1 = 1.0 - t1/ta;
787 
788  if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0;
789 
790 /*
791  std::cout << "Collision as " << i << " " << j << " " << as << std::endl;
792  std::cout << "Collision a " << i << " " << j << " " << a << std::endl;
793  std::cout << "Collision ta " << i << " " << j << " " << ta << std::endl;
794  std::cout << "Collision x " << i << " " << j << " " << x << std::endl;
795  std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl;
796  std::cout << "Collision c1 " << i << " " << j << " " << c1 << std::endl;
797 */
798  t1 = 2.0*pi*G4UniformRand();
799 // std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl;
800  G4double t2 = 0.0;
801  if ( pcm.x() == 0.0 && pcm.y() == 0 )
802  {
803  t2 = 0.0;
804  }
805  else
806  {
807  t2 = std::atan2( pcm.y() , pcm.x() );
808  }
809 // std::cout << "Collision t2 " << i << " " << j << " " << t2 << std::endl;
810 
811  G4double s1 = std::sqrt ( 1.0 - c1*c1 );
812  G4double s2 = std::sqrt ( 1.0 - c2*c2 );
813 
814  G4double ct1 = std::cos(t1);
815  G4double st1 = std::sin(t1);
816 
817  G4double ct2 = std::cos(t2);
818  G4double st2 = std::sin(t2);
819 
820  G4double ss = c2*s1*ct1 + s2*c1;
821 
822  pcm.setX( pr * ( ss*ct2 - s1*st1*st2) );
823  pcm.setY( pr * ( ss*st2 + s1*st1*ct2) );
824  pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) );
825 
826 // std::cout << "Collision pcm " << i << " " << j << " " << pcm << std::endl;
827 
828  G4double epot = theMeanField->GetTotalPotential();
829 
830  G4double eini = epot + p4i.e() + p4j.e();
831  G4double etwo = p4i.e() + p4j.e();
832 
833 /*
834  std::cout << "Collision epot " << i << " " << j << " " << epot << std::endl;
835  std::cout << "Collision eini " << i << " " << j << " " << eini << std::endl;
836  std::cout << "Collision etwo " << i << " " << j << " " << etwo << std::endl;
837 */
838 
839 
840  for ( G4int itry = 0 ; itry < 4 ; itry++ )
841  {
842 
843  G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm );
844  G4double pibeta = pcm*beta;
845 
846  G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm );
847 
848  G4ThreeVector pi_new = beta*trans + pcm;
849 
850  G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm );
851  trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm );
852 
853  G4ThreeVector pj_new = beta*trans - pcm;
854 
855 //
856 // Delete old
857 // Add new Particitipants
858 //
859 // Now only change momentum ( Beacuse we only have elastic sctter of nucleon
860 // In future Definition also will be change
861 //
862 
863  theSystem->GetParticipant( i )->SetMomentum( pi_new );
864  theSystem->GetParticipant( j )->SetMomentum( pj_new );
865 
866  G4double pi_new_e = (theSystem->GetParticipant( i )->Get4Momentum()).e();
867  G4double pj_new_e = (theSystem->GetParticipant( j )->Get4Momentum()).e();
868 
869  theMeanField->Cal2BodyQuantities( i );
870  theMeanField->Cal2BodyQuantities( j );
871 
872  epot = theMeanField->GetTotalPotential();
873 
874  G4double efin = epot + pi_new_e + pj_new_e ;
875 
876  //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - epse << std::endl;
877 /*
878  std::cout << "Collision efin " << i << " " << j << " " << efin << std::endl;
879  std::cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << epse << std::endl;
880  std::cout << "Collision " << std::abs ( eini - efin ) << " " << epse << std::endl;
881 */
882 
883 //071031
884  if ( std::abs ( eini - efin ) < epse )
885  {
886  // Collison OK
887  //std::cout << "collisions6" << std::endl;
888  //std::cout << "collisions before " << p4i << " " << p4j << std::endl;
889  //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl;
890  //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl;
891  //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl;
892  //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl;
893  }
894 //071031
895 
896  if ( std::abs ( eini - efin ) < epse ) return result; // Collison OK
897 
898  G4double cona = ( eini - efin + etwo ) / gamma;
899  G4double fac2 = 1.0 / ( 4.0 * cona*cona * pr*pr ) *
900  ( ( cona*cona - ( rmi*rmi + rmj*rmj ) )*( cona*cona - ( rmi*rmi + rmj*rmj ) )
901  - 4.0 * rmi*rmi * rmj*rmj );
902 
903  if ( fac2 > 0 )
904  {
905  G4double fact = std::sqrt ( fac2 );
906  pcm = fact*pcm;
907  }
908 
909 
910  }
911 
912 // Energetically forbidden collision
913  result = false;
914 
915  return result;
916 
917 }
Generated on Sat May 25 2013 14:34:06 for Geant4 by   doxygen 1.8.4