ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QMDReaction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4QMDReaction.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 // 080505 Fixed and changed sampling method of impact parameter by T. Koi
27 // 080602 Fix memory leaks by T. Koi
28 // 080612 Delete unnecessary dependency and unused functions
29 // Change criterion of reaction by T. Koi
30 // 081107 Add UnUseGEM (then use the default channel of G4Evaporation)
31 // UseFrag (chage criterion of a inelastic reaction)
32 // Fix bug in nucleon projectiles by T. Koi
33 // 090122 Be8 -> Alpha + Alpha
34 // 090331 Change member shenXS and genspaXS object to pointer
35 // 091119 Fix for incidence of neutral particles
36 //
37 #include "G4QMDReaction.hh"
38 #include "G4QMDNucleus.hh"
40 #include "G4Pow.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4NistManager.hh"
44 
46 
48 : G4HadronicInteraction("QMDModel")
49 , system ( NULL )
50 , deltaT ( 1 ) // in fsec (c=1)
51 , maxTime ( 100 ) // will have maxTime-th time step
52 , envelopF ( 1.05 ) // 10% for Peripheral reactions
53 , gem ( true )
54 , frag ( false )
55 {
56 
57  //090331
59  //genspaXS = new G4GeneralSpaceNNCrossSection();
61  meanField = new G4QMDMeanField();
62  collision = new G4QMDCollision();
63 
66  //fEvaporation - 8 default channels
67  //fCombined - 8 default + 60 GEM
68  //fGEM - 2 default + 66 GEM
72 
78 
84 
85 }
86 
87 
88 
90 {
91  delete evaporation;
92  delete excitationHandler;
93  delete collision;
94  delete meanField;
95 }
96 
97 
98 
100 {
101  //G4cout << "G4QMDReaction::ApplyYourself" << G4endl;
102 
104 
105  system = new G4QMDSystem;
106 
107  G4int proj_Z = 0;
108  G4int proj_A = 0;
109  const G4ParticleDefinition* proj_pd = ( const G4ParticleDefinition* ) projectile.GetDefinition();
110  if ( proj_pd->GetParticleType() == "nucleus" )
111  {
112  proj_Z = proj_pd->GetAtomicNumber();
113  proj_A = proj_pd->GetAtomicMass();
114  }
115  else
116  {
117  proj_Z = (int)( proj_pd->GetPDGCharge()/eplus );
118  proj_A = 1;
119  }
120  //G4int targ_Z = int ( target.GetZ() + 0.5 );
121  //G4int targ_A = int ( target.GetN() + 0.5 );
122  //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
123  G4int targ_Z = target.GetZ_asInt();
124  G4int targ_A = target.GetA_asInt();
125  const G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon( targ_Z , targ_A , 0.0 );
126 
127 
128  //G4NistManager* nistMan = G4NistManager::Instance();
129 // G4Element* G4NistManager::FindOrBuildElement( targ_Z );
130 
131  const G4DynamicParticle* proj_dp = new G4DynamicParticle ( proj_pd , projectile.Get4Momentum() );
132  //const G4Element* targ_ele = nistMan->FindOrBuildElement( targ_Z );
133  //G4double aTemp = projectile.GetMaterial()->GetTemperature();
134 
135  //090331
136 
138 
139  if ( proj_pd->GetParticleType() == "meson" ) theXS = piNucXS;
140 
141  //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
142  //G4double xs_0 = theXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
143  //110822
144  G4double xs_0 = theXS->GetIsoCrossSection ( proj_dp , targ_Z , targ_A );
145 
146  G4double bmax_0 = std::sqrt( xs_0 / pi );
147  //std::cout << "bmax_0 in fm (fermi) " << bmax_0/fermi << std::endl;
148 
149  //delete proj_dp;
150 
151  G4bool elastic = true;
152 
153  std::vector< G4QMDNucleus* > nucleuses; // Secondary nuceluses
154  G4ThreeVector boostToReac; // ReactionSystem (CM or NN);
155  G4ThreeVector boostBackToLAB; // Reaction System to LAB;
156 
157  G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ_pd->GetPDGMass()/GeV );
158  G4ThreeVector boostLABtoCM = targ4p.findBoostToCM( proj_dp->Get4Momentum()/GeV ); // CM of target and proj;
159 
160  G4double p1 = proj_dp->GetMomentum().mag()/GeV/proj_A;
161  G4double m1 = proj_dp->GetDefinition()->GetPDGMass()/GeV/proj_A;
162  G4double e1 = std::sqrt( p1*p1 + m1*m1 );
163  G4double e2 = targ_pd->GetPDGMass()/GeV/targ_A;
164  G4double beta_nn = -p1 / ( e1+e2 );
165 
166  G4ThreeVector boostLABtoNN ( 0. , 0. , beta_nn ); // CM of NN;
167 
168  G4double beta_nncm = ( - boostLABtoCM.beta() + boostLABtoNN.beta() ) / ( 1 - boostLABtoCM.beta() * boostLABtoNN.beta() ) ;
169 
170  //std::cout << targ4p << std::endl;
171  //std::cout << proj_dp->Get4Momentum()<< std::endl;
172  //std::cout << beta_nncm << std::endl;
173  G4ThreeVector boostNNtoCM( 0. , 0. , beta_nncm ); //
174  G4ThreeVector boostCMtoNN( 0. , 0. , -beta_nncm ); //
175 
176  boostToReac = boostLABtoNN;
177  boostBackToLAB = -boostLABtoNN;
178 
179  delete proj_dp;
180 
181  G4int icounter = 0;
182  G4int icounter_max = 1024;
183  while ( elastic ) // Loop checking, 11.03.2015, T. Koi
184  {
185  icounter++;
186  if ( icounter > icounter_max ) {
187  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
188  break;
189  }
190 
191 // impact parameter
192  //G4double bmax = 1.05*(bmax_0/fermi); // 10% for Peripheral reactions
193  G4double bmax = envelopF*(bmax_0/fermi);
194  G4double b = bmax * std::sqrt ( G4UniformRand() );
195 //071112
196  //G4double b = 0;
197  //G4double b = bmax;
198  //G4double b = bmax/1.05 * 0.7 * G4UniformRand();
199 
200  //G4cout << "G4QMDRESULT bmax_0 = " << bmax_0/fermi << " fm, bmax = " << bmax << " fm , b = " << b << " fm " << G4endl;
201 
202  G4double plab = projectile.GetTotalMomentum()/GeV;
203  G4double elab = ( projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV;
204 
205  calcOffSetOfCollision( b , proj_pd , targ_pd , plab , elab , bmax , boostCMtoNN );
206 
207 // Projectile
208  G4LorentzVector proj4pLAB = projectile.Get4Momentum()/GeV;
209 
211  if ( projectile.GetDefinition()->GetParticleType() == "nucleus"
212  || projectile.GetDefinition()->GetParticleName() == "proton"
213  || projectile.GetDefinition()->GetParticleName() == "neutron" )
214  {
215 
216  proj_Z = proj_pd->GetAtomicNumber();
217  proj_A = proj_pd->GetAtomicMass();
218 
219  proj = new G4QMDGroundStateNucleus( proj_Z , proj_A );
220  //proj->ShowParticipants();
221 
222 
223  meanField->SetSystem ( proj );
226 
227  }
228 
229 // Target
230  //G4int iz = int ( target.GetZ() );
231  //G4int ia = int ( target.GetN() );
232  //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
233  G4int iz = int ( target.GetZ_asInt() );
234  G4int ia = int ( target.GetA_asInt() );
235 
237 
238  meanField->SetSystem (targ );
239  targ->SetTotalPotential( meanField->GetTotalPotential() );
240  targ->CalEnergyAndAngularMomentumInCM();
241 
242  //G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ->GetNuclearMass()/GeV );
243 // Boost Vector to CM
244  //boostToCM = targ4p.findBoostToCM( proj4pLAB );
245 
246 // Target
247  for ( G4int i = 0 ; i < targ->GetTotalNumberOfParticipant() ; i++ )
248  {
249 
250  G4ThreeVector p0 = targ->GetParticipant( i )->GetMomentum();
251  G4ThreeVector r0 = targ->GetParticipant( i )->GetPosition();
252 
254  , p0.y()
256 
258  , r0.y()
260 
261  system->SetParticipant( new G4QMDParticipant( targ->GetParticipant( i )->GetDefinition() , p , r ) );
263 
264  }
265 
266  G4LorentzVector proj4pCM = CLHEP::boostOf ( proj4pLAB , boostToReac );
267  G4LorentzVector targ4pCM = CLHEP::boostOf ( targ4p , boostToReac );
268 
269 // Projectile
270  if ( proj != NULL )
271  {
272 
273 // projectile is nucleus
274 
275  for ( G4int i = 0 ; i < proj->GetTotalNumberOfParticipant() ; i++ )
276  {
277 
278  G4ThreeVector p0 = proj->GetParticipant( i )->GetMomentum();
279  G4ThreeVector r0 = proj->GetParticipant( i )->GetPosition();
280 
282  , p0.y()
284 
286  , r0.y()
288 
290  system->GetParticipant ( i + targ->GetTotalNumberOfParticipant() )->SetProjectile();
291  }
292 
293  }
294  else
295  {
296 
297 // projectile is particle
298 
299  // avoid multiple set in "elastic" loop
300  if ( system->GetTotalNumberOfParticipant() == targ->GetTotalNumberOfParticipant() )
301  {
302 
303  G4int i = targ->GetTotalNumberOfParticipant();
304 
305  G4ThreeVector p0( 0 );
306  G4ThreeVector r0( 0 );
307 
309  , p0.y()
311 
313  , r0.y()
315 
317  // This is not important becase only 1 projectile particle.
319  }
320 
321  }
322  //system->ShowParticipants();
323 
324  delete targ;
325  delete proj;
326 
329 
330 // Time Evolution
331  //std::cout << "Start time evolution " << std::endl;
332  //system->ShowParticipants();
333  for ( G4int i = 0 ; i < maxTime ; i++ )
334  {
335  //G4cout << " do Paropagate " << i << " th time step. " << G4endl;
337  //system->ShowParticipants();
339 
340  if ( i / 10 * 10 == i )
341  {
342  //G4cout << i << " th time step. " << G4endl;
343  //system->ShowParticipants();
344  }
345  //system->ShowParticipants();
346  }
347  //system->ShowParticipants();
348 
349 
350  //std::cout << "Doing Cluster Judgment " << std::endl;
351 
352  nucleuses = meanField->DoClusterJudgment();
353 
354 // Elastic Judgment
355 
356  G4int numberOfSecondary = int ( nucleuses.size() ) + system->GetTotalNumberOfParticipant();
357 
358  G4int sec_a_Z = 0;
359  G4int sec_a_A = 0;
360  const G4ParticleDefinition* sec_a_pd = NULL;
361  G4int sec_b_Z = 0;
362  G4int sec_b_A = 0;
363  const G4ParticleDefinition* sec_b_pd = NULL;
364 
365  if ( numberOfSecondary == 2 )
366  {
367 
368  G4bool elasticLike_system = false;
369  if ( nucleuses.size() == 2 )
370  {
371 
372  sec_a_Z = nucleuses[0]->GetAtomicNumber();
373  sec_a_A = nucleuses[0]->GetMassNumber();
374  sec_b_Z = nucleuses[1]->GetAtomicNumber();
375  sec_b_A = nucleuses[1]->GetMassNumber();
376 
377  if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_Z == targ_Z && sec_b_A == targ_A )
378  || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_Z == proj_Z && sec_b_A == proj_A ) )
379  {
380  elasticLike_system = true;
381  }
382 
383  }
384  else if ( nucleuses.size() == 1 )
385  {
386 
387  sec_a_Z = nucleuses[0]->GetAtomicNumber();
388  sec_a_A = nucleuses[0]->GetMassNumber();
389  sec_b_pd = system->GetParticipant( 0 )->GetDefinition();
390 
391  if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_pd == targ_pd )
392  || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_pd == proj_pd ) )
393  {
394  elasticLike_system = true;
395  }
396 
397  }
398  else
399  {
400 
401  sec_a_pd = system->GetParticipant( 0 )->GetDefinition();
402  sec_b_pd = system->GetParticipant( 1 )->GetDefinition();
403 
404  if ( ( sec_a_pd == proj_pd && sec_b_pd == targ_pd )
405  || ( sec_a_pd == targ_pd && sec_b_pd == proj_pd ) )
406  {
407  elasticLike_system = true;
408  }
409 
410  }
411 
412  if ( elasticLike_system == true )
413  {
414 
415  G4bool elasticLike_energy = true;
416 // Cal ExcitationEnergy
417  for ( G4int i = 0 ; i < int ( nucleuses.size() ) ; i++ )
418  {
419 
420  //meanField->SetSystem( nucleuses[i] );
421  meanField->SetNucleus( nucleuses[i] );
422  //nucleuses[i]->SetTotalPotential( meanField->GetTotalPotential() );
423  //nucleuses[i]->CalEnergyAndAngularMomentumInCM();
424 
425  if ( nucleuses[i]->GetExcitationEnergy()*GeV > 1.0*MeV ) elasticLike_energy = false;
426 
427  }
428 
429 // Check Collision
430  G4bool withCollision = true;
431  if ( system->GetNOCollision() == 0 ) withCollision = false;
432 
433 // Final judegement for Inelasitc or Elastic;
434 //
435 // ElasticLike without Collision
436  //if ( elasticLike_energy == true && withCollision == false ) elastic = true; // ielst = 0
437 // ElasticLike with Collision
438  //if ( elasticLike_energy == true && withCollision == true ) elastic = true; // ielst = 1
439 // InelasticLike without Collision
440  //if ( elasticLike_energy == false ) elastic = false; // ielst = 2
441  if ( frag == true )
442  if ( elasticLike_energy == false ) elastic = false;
443 // InelasticLike with Collision
444  if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3
445 
446  }
447 
448  }
449  else
450  {
451 
452 // numberOfSecondary != 2
453  elastic = false;
454 
455  }
456 
457 //071115
458  //G4cout << elastic << G4endl;
459  // if elastic is true try again from sampling of impact parameter
460 
461  if ( elastic == true )
462  {
463  // delete this nucleues
464  for ( std::vector< G4QMDNucleus* >::iterator
465  it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
466  {
467  delete *it;
468  }
469  nucleuses.clear();
470  }
471 
472  }
473 
474 
475 // Statical Decay Phase
476 
477  for ( std::vector< G4QMDNucleus* >::iterator it
478  = nucleuses.begin() ; it != nucleuses.end() ; it++ )
479  {
480 
481 /*
482  G4cout << "G4QMDRESULT "
483  << (*it)->GetAtomicNumber()
484  << " "
485  << (*it)->GetMassNumber()
486  << " "
487  << (*it)->Get4Momentum()
488  << " "
489  << (*it)->Get4Momentum().vect()
490  << " "
491  << (*it)->Get4Momentum().restMass()
492  << " "
493  << (*it)->GetNuclearMass()/GeV
494  << G4endl;
495 */
496 
497  meanField->SetNucleus ( *it );
498 
499  if ( (*it)->GetAtomicNumber() == 0 // neutron cluster
500  || (*it)->GetAtomicNumber() == (*it)->GetMassNumber() ) // proton cluster
501  {
502  // push back system
503  for ( G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ )
504  {
505  G4QMDParticipant* aP = new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() );
506  system->SetParticipant ( aP );
507  }
508  continue;
509  }
510 
511  G4double nucleus_e = std::sqrt ( G4Pow::GetInstance()->powN ( (*it)->GetNuclearMass()/GeV , 2 ) + G4Pow::GetInstance()->powN ( (*it)->Get4Momentum().vect().mag() , 2 ) );
512  G4LorentzVector nucleus_p4CM ( (*it)->Get4Momentum().vect() , nucleus_e );
513 
514 // std::cout << "G4QMDRESULT nucleus deltaQ " << deltaQ << std::endl;
515 
516  G4int ia = (*it)->GetMassNumber();
517  G4int iz = (*it)->GetAtomicNumber();
518 
519  G4LorentzVector lv ( G4ThreeVector( 0.0 ) , (*it)->GetExcitationEnergy()*GeV + G4IonTable::GetIonTable()->GetIonMass( iz , ia ) );
520 
521  G4Fragment* aFragment = new G4Fragment( ia , iz , lv );
522 
524  rv = excitationHandler->BreakItUp( *aFragment );
525  G4bool notBreak = true;
526  for ( G4ReactionProductVector::iterator itt
527  = rv->begin() ; itt != rv->end() ; itt++ )
528  {
529 
530  notBreak = false;
531  // Secondary from this nucleus (*it)
532  const G4ParticleDefinition* pd = (*itt)->GetDefinition();
533 
534  G4LorentzVector p4 ( (*itt)->GetMomentum()/GeV , (*itt)->GetTotalEnergy()/GeV ); //in nucleus(*it) rest system
535  G4LorentzVector p4_CM = CLHEP::boostOf( p4 , -nucleus_p4CM.findBoostToCM() ); // Back to CM
536  G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB
537 
538 
539 //090122
540  //theParticleChange.AddSecondary( dp );
541  if ( !( pd->GetAtomicNumber() == 4 && pd->GetAtomicMass() == 8 ) )
542  {
543  G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );
545  }
546  else
547  {
548  //Be8 -> Alpha + Alpha + Q
549  G4ThreeVector randomized_direction( G4UniformRand() , G4UniformRand() , G4UniformRand() );
550  randomized_direction = randomized_direction.unit();
551  G4double q_decay = (*itt)->GetMass() - 2*G4Alpha::Alpha()->GetPDGMass();
552  G4double p_decay = std::sqrt ( G4Pow::GetInstance()->powN(G4Alpha::Alpha()->GetPDGMass()+q_decay/2,2) - G4Pow::GetInstance()->powN(G4Alpha::Alpha()->GetPDGMass() , 2 ) );
553  G4LorentzVector p4_a1 ( p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 ); //in Be8 rest system
554 
555  G4LorentzVector p4_a1_Be8 = CLHEP::boostOf ( p4_a1/GeV , -p4.findBoostToCM() );
556  G4LorentzVector p4_a1_CM = CLHEP::boostOf ( p4_a1_Be8 , -nucleus_p4CM.findBoostToCM() );
557  G4LorentzVector p4_a1_LAB = CLHEP::boostOf ( p4_a1_CM , boostBackToLAB );
558 
559  G4LorentzVector p4_a2 ( -p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 ); //in Be8 rest system
560 
561  G4LorentzVector p4_a2_Be8 = CLHEP::boostOf ( p4_a2/GeV , -p4.findBoostToCM() );
562  G4LorentzVector p4_a2_CM = CLHEP::boostOf ( p4_a2_Be8 , -nucleus_p4CM.findBoostToCM() );
563  G4LorentzVector p4_a2_LAB = CLHEP::boostOf ( p4_a2_CM , boostBackToLAB );
564 
565  G4DynamicParticle* dp1 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a1_LAB*GeV );
566  G4DynamicParticle* dp2 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a2_LAB*GeV );
569  }
570 //090122
571 
572 /*
573  G4cout
574  << "Regist Secondary "
575  << (*itt)->GetDefinition()->GetParticleName()
576  << " "
577  << (*itt)->GetMomentum()/GeV
578  << " "
579  << (*itt)->GetKineticEnergy()/GeV
580  << " "
581  << (*itt)->GetMass()/GeV
582  << " "
583  << (*itt)->GetTotalEnergy()/GeV
584  << " "
585  << (*itt)->GetTotalEnergy()/GeV * (*itt)->GetTotalEnergy()/GeV
586  - (*itt)->GetMomentum()/GeV * (*itt)->GetMomentum()/GeV
587  << " "
588  << nucleus_p4CM.findBoostToCM()
589  << " "
590  << p4
591  << " "
592  << p4_CM
593  << " "
594  << p4_LAB
595  << G4endl;
596 */
597 
598  }
599  if ( notBreak == true )
600  {
601 
602  const G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( (*it)->GetAtomicNumber() , (*it)->GetMassNumber(), (*it)->GetExcitationEnergy()*GeV );
603  G4LorentzVector p4_CM = nucleus_p4CM;
604  G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB
605  G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );
607 
608  }
609 
610  for ( G4ReactionProductVector::iterator itt
611  = rv->begin() ; itt != rv->end() ; itt++ )
612  {
613  delete *itt;
614  }
615  delete rv;
616 
617  delete aFragment;
618 
619  }
620 
621 
622 
623  for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i++ )
624  {
625 
626  // Secondary particles
627 
630  G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB );
631  G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );
633 
634 /*
635  G4cout << "G4QMDRESULT "
636  << "r" << i << " " << system->GetParticipant ( i ) -> GetPosition() << " "
637  << "p" << i << " " << system->GetParticipant ( i ) -> Get4Momentum()
638  << G4endl;
639 */
640 
641  }
642 
643  for ( std::vector< G4QMDNucleus* >::iterator it
644  = nucleuses.begin() ; it != nucleuses.end() ; it++ )
645  {
646  delete *it; // delete nulceuse
647  }
648  nucleuses.clear();
649 
650  system->Clear();
651  delete system;
652 
654 
655  return &theParticleChange;
656 
657 }
658 
659 
660 
662 const G4ParticleDefinition* pd_proj ,
663 const G4ParticleDefinition* pd_targ ,
664 G4double ptot , G4double etot , G4double bmax , G4ThreeVector boostToCM )
665 {
666 
667  G4double mass_proj = pd_proj->GetPDGMass()/GeV;
668  G4double mass_targ = pd_targ->GetPDGMass()/GeV;
669 
670  G4double stot = std::sqrt ( etot*etot - ptot*ptot );
671 
672  G4double pstt = std::sqrt ( ( stot*stot - ( mass_proj + mass_targ ) * ( mass_proj + mass_targ )
673  ) * ( stot*stot - ( mass_proj - mass_targ ) * ( mass_proj - mass_targ ) ) )
674  / ( 2.0 * stot );
675 
676  G4double pzcc = pstt;
677  G4double eccm = stot - ( mass_proj + mass_targ );
678 
679  G4int zp = 1;
680  G4int ap = 1;
681  if ( pd_proj->GetParticleType() == "nucleus" )
682  {
683  zp = pd_proj->GetAtomicNumber();
684  ap = pd_proj->GetAtomicMass();
685  }
686  else
687  {
688  // proton, neutron, mesons
689  zp = int ( pd_proj->GetPDGCharge()/eplus + 0.5 );
690  // ap = 1;
691  }
692 
693 
694  G4int zt = pd_targ->GetAtomicNumber();
695  G4int at = pd_targ->GetAtomicMass();
696 
697 
698  //G4double rmax0 = 8.0; // T.K dicide parameter value // for low energy
699  G4double rmax0 = bmax + 4.0;
700  G4double rmax = std::sqrt( rmax0*rmax0 + b*b );
701 
702  G4double ccoul = 0.001439767;
703  G4double pcca = 1.0 - double ( zp * zt ) * ccoul / eccm / rmax - ( b / rmax )*( b / rmax );
704 
705  G4double pccf = std::sqrt( pcca );
706 
707  //Fix for neutral particles
708  G4double aas1 = 0.0;
709  G4double bbs = 0.0;
710 
711  if ( zp != 0 )
712  {
713  G4double aas = 2.0 * eccm * b / double ( zp * zt ) / ccoul;
714  bbs = 1.0 / std::sqrt ( 1.0 + aas*aas );
715  aas1 = ( 1.0 + aas * b / rmax ) * bbs;
716  }
717 
718  G4double cost = 0.0;
719  G4double sint = 0.0;
720  G4double thet1 = 0.0;
721  G4double thet2 = 0.0;
722  if ( 1.0 - aas1*aas1 <= 0 || 1.0 - bbs*bbs <= 0.0 )
723  {
724  cost = 1.0;
725  sint = 0.0;
726  }
727  else
728  {
729  G4double aat1 = aas1 / std::sqrt ( 1.0 - aas1*aas1 );
730  G4double aat2 = bbs / std::sqrt ( 1.0 - bbs*bbs );
731 
732  thet1 = std::atan ( aat1 );
733  thet2 = std::atan ( aat2 );
734 
735 // TK enter to else block
736  G4double theta = thet1 - thet2;
737  cost = std::cos( theta );
738  sint = std::sin( theta );
739  }
740 
741  G4double rzpr = -rmax * cost * ( mass_targ ) / ( mass_proj + mass_targ );
742  G4double rzta = rmax * cost * ( mass_proj ) / ( mass_proj + mass_targ );
743 
744  G4double rxpr = rmax / 2.0 * sint;
745 
746  G4double rxta = -rxpr;
747 
748 
749  G4double pzpc = pzcc * ( cost * pccf + sint * b / rmax );
750  G4double pxpr = pzcc * ( -sint * pccf + cost * b / rmax );
751 
752  G4double pztc = - pzpc;
753  G4double pxta = - pxpr;
754 
755  G4double epc = std::sqrt ( pzpc*pzpc + pxpr*pxpr + mass_proj*mass_proj );
756  G4double etc = std::sqrt ( pztc*pztc + pxta*pxta + mass_targ*mass_targ );
757 
758  G4double pzpr = pzpc;
759  G4double pzta = pztc;
760  G4double epr = epc;
761  G4double eta = etc;
762 
763 // CM -> NN
764  G4double gammacm = boostToCM.gamma();
765  //G4double betacm = -boostToCM.beta();
766  G4double betacm = boostToCM.z();
767  pzpr = pzpc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pzpc * betacm + epc );
768  pzta = pztc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pztc * betacm + etc );
769  epr = gammacm * ( epc + betacm * pzpc );
770  eta = gammacm * ( etc + betacm * pztc );
771 
772  //G4double betpr = pzpr / epr;
773  //G4double betta = pzta / eta;
774 
775  G4double gammpr = epr / ( mass_proj );
776  G4double gammta = eta / ( mass_targ );
777 
778  pzta = pzta / double ( at );
779  pxta = pxta / double ( at );
780 
781  pzpr = pzpr / double ( ap );
782  pxpr = pxpr / double ( ap );
783 
784  G4double zeroz = 0.0;
785 
786  rzpr = rzpr -zeroz;
787  rzta = rzta -zeroz;
788 
789  // Set results
795 
801 
802 }
803 
804 
805 
807 {
808 
809  if ( gem == true )
811  else
813 
814 }
815 
816 void G4QMDReaction::ModelDescription(std::ostream& outFile) const
817 {
818  outFile << "Lorentz covarianted Quantum Molecular Dynamics model for nucleus (particle) vs nucleus reactions\n";
819 }