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