Geant4  10.02.p01
G4GeneratorPrecompoundInterface.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 // $Id: G4GeneratorPrecompoundInterface.cc 92692 2015-09-14 07:06:19Z gcosmo $
27 //
28 // -----------------------------------------------------------------------------
29 // GEANT 4 class file
30 //
31 // History: first implementation
32 // HPW, 10DEC 98, the decay part originally written by Gunter Folger
33 // in his FTF-test-program.
34 //
35 // M.Kelsey, 28 Jul 2011 -- Replace loop to decay input secondaries
36 // with new utility class, simplify cleanup loops
37 // -----------------------------------------------------------------------------
38 
39 #include <algorithm>
40 #include <vector>
41 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
46 #include "G4KineticTrackVector.hh"
47 #include "G4Proton.hh"
48 #include "G4Neutron.hh"
49 
50 #include "G4Deuteron.hh"
51 #include "G4Triton.hh"
52 #include "G4He3.hh"
53 #include "G4Alpha.hh"
54 
55 #include "G4V3DNucleus.hh"
56 #include "G4Nucleon.hh"
57 
58 #include "G4AntiProton.hh"
59 #include "G4AntiNeutron.hh"
60 #include "G4AntiDeuteron.hh"
61 #include "G4AntiTriton.hh"
62 #include "G4AntiHe3.hh"
63 #include "G4AntiAlpha.hh"
64 
65 #include "G4FragmentVector.hh"
66 #include "G4ReactionProduct.hh"
68 #include "G4PreCompoundModel.hh"
69 #include "G4ExcitationHandler.hh"
70 #include "G4DecayKineticTracks.hh"
72 
73 //---------------------------------------------------------------------
74 #include "Randomize.hh"
75 #include "G4Log.hh"
76 
77 //#define debugPrecoInt
78 
80 : CaptureThreshold(70*MeV) // Uzhi 1.05.2015 10 ->70
81 {
84 
87  He3 =G4He3::He3();
89 
92 
97 
98  if(preModel) { SetDeExcitation(preModel); }
99  else {
100  G4HadronicInteraction* hadi =
102  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
103  if(!pre) { pre = new G4PreCompoundModel(); }
104  SetDeExcitation(pre);
105  }
106 }
107 
109 {
110 }
111 
113 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
114 {
115  #ifdef debugPrecoInt
116  G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::Propagate"<<G4endl;
117  G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
118  G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
119  #endif
120 
121  G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
122 
123  // decay the strong resonances
124  G4DecayKineticTracks decay(theSecondaries);
125  #ifdef debugPrecoInt
126  G4cout<<"Final stable particles number "<<theSecondaries->size()<<G4endl;
127  #endif
128 
129  // prepare the fragment
130  G4int anA=theNucleus->GetMassNumber();
131  G4int aZ=theNucleus->GetCharge();
132 // G4double TargetNucleusMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
133 
134  G4int numberOfEx = 0;
135  G4int numberOfCh = 0;
136  G4int numberOfHoles = 0;
137 
138  G4double R = theNucleus->GetNuclearRadius();
139 
140  G4LorentzVector captured4Momentum(0.,0.,0.,0.);
141  G4LorentzVector Residual4Momentum(0.,0.,0.,0.); // TargetNucleusMass is not need at the moment
142  G4LorentzVector Secondary4Momentum(0.,0.,0.,0.);
143 
144  // loop over secondaries
145  G4KineticTrackVector::iterator iter;
146  for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
147  {
148  const G4ParticleDefinition* part = (*iter)->GetDefinition();
149  G4double e = (*iter)->Get4Momentum().e();
150  G4double mass = (*iter)->Get4Momentum().mag();
151  G4ThreeVector mom = (*iter)->Get4Momentum().vect();
152  if((part != proton && part != neutron) ||
153 // Uzhi 2.05.2015 (e > mass + CaptureThreshold) ||
154  ((*iter)->GetPosition().mag() > R)) {
155  G4ReactionProduct * theNew = new G4ReactionProduct(part);
156  theNew->SetMomentum(mom);
157  theNew->SetTotalEnergy(e);
158  theTotalResult->push_back(theNew);
159  Secondary4Momentum += (*iter)->Get4Momentum(); // Uzhi 29 April
160  #ifdef debugPrecoInt
161  G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
162  <<(*iter)->Get4Momentum().mag()<<G4endl;
163  #endif
164  } else {
165  if( e-mass > -CaptureThreshold*G4Log( G4UniformRand()) ) { // Added by Uzhi 2.05.2015
166  G4ReactionProduct * theNew = new G4ReactionProduct(part);
167  theNew->SetMomentum(mom);
168  theNew->SetTotalEnergy(e);
169  theTotalResult->push_back(theNew);
170  Secondary4Momentum += (*iter)->Get4Momentum(); // Uzhi 29 April
171  #ifdef debugPrecoInt
172  G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
173  <<(*iter)->Get4Momentum().mag()<<G4endl;
174  #endif
175  } else {
176  // within the nucleus, neutron or proton
177  // now calculate A, Z of the fragment, momentum, number of exciton states
178  ++anA;
179  ++numberOfEx;
180  G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
181  aZ += Z;
182  numberOfCh += Z;
183  captured4Momentum += (*iter)->Get4Momentum();
184  #ifdef debugPrecoInt
185  G4cout<<"Captured 4Mom "<<part->GetParticleName()<<(*iter)->Get4Momentum()<<G4endl;
186  #endif
187  }
188  }
189  delete (*iter);
190  }
191  delete theSecondaries;
192 
193  // loop over wounded nucleus
194  G4Nucleon * theCurrentNucleon =
195  theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
196  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
197  {
198  if(theCurrentNucleon->AreYouHit()) {
199  ++numberOfHoles;
200  ++numberOfEx;
201  --anA;
202  aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
203 
204  Residual4Momentum -= theCurrentNucleon->Get4Momentum();
205  }
206  theCurrentNucleon = theNucleus->GetNextNucleon();
207  }
208 
209  #ifdef debugPrecoInt
210  G4cout<<G4endl;
211  G4cout<<"Secondary 4Mom "<<Secondary4Momentum<<G4endl;
212  G4cout<<"Captured 4Mom "<<captured4Momentum<<G4endl;
213  G4cout<<"Sec + Captured "<<Secondary4Momentum+captured4Momentum<<G4endl;
214  G4cout<<"Residual4Mom "<<Residual4Momentum<<G4endl;
215  G4cout<<"Sum 4 momenta "
216  <<Secondary4Momentum + captured4Momentum + Residual4Momentum <<G4endl;
217  #endif
218 
219  // Check that we use QGS model; loop over wounded nucleus
220  G4bool QGSM(false);
221  theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
222  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
223  {
224  if(theCurrentNucleon->AreYouHit())
225  {
226  if(theCurrentNucleon->Get4Momentum().mag() <
227  theCurrentNucleon->GetDefinition()->GetPDGMass()) QGSM=true;
228  }
229  theCurrentNucleon = theNucleus->GetNextNucleon();
230  }
231 
232  #ifdef debugPrecoInt
233  if(!QGSM){
234  G4cout<<G4endl;
235  G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
236  G4cout<<"Residual 4Mom "<<Residual4Momentum<<G4endl;
237  if(numberOfEx == 0)
238  {G4cout<<"Residual 4Mom = 0 means that there were not wounded and captured nucleons"<<G4endl;}
239  }
240  #endif
241 
242  if(anA == 0) return theTotalResult;
243 
244  G4LorentzVector exciton4Momentum(0.,0.,0.,0.); // Uzhi 29 April
245  if(anA >= aZ)
246  {
247  if(!QGSM)
248  { // FTF model was used
250 
251 // G4LorentzVector exciton4Momentum = Residual4Momentum + captured4Momentum;
252  exciton4Momentum = Residual4Momentum + captured4Momentum;
253 //exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
254  G4double ActualMass = exciton4Momentum.mag();
255  if(ActualMass <= fMass ) { //E*<=0, Uzhi 5.05.2015
256  exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass))); // Uzhi 13.05.2015
257  }
258 
259  #ifdef debugPrecoInt
260  G4double exEnergy = 0.0;
261  if(ActualMass <= fMass ) {exEnergy = 0.;} // Uzhi 5.05.2015
262  else {exEnergy = ActualMass - fMass;}
263  G4cout<<"Ground state residual Mass "<<fMass<<" E* "<<exEnergy<<G4endl;
264  #endif
265  }
266  else
267  { // QGS model was used
268  G4double InitialTargetMass =
269  G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(), theNucleus->GetCharge());
270 
271  exciton4Momentum =
272  GetPrimaryProjectile()->Get4Momentum() + G4LorentzVector(0.,0.,0.,InitialTargetMass)
273  -Secondary4Momentum;
274 
276  G4double ActualMass = exciton4Momentum.mag();
277 
278  #ifdef debugPrecoInt
279  G4cout<<G4endl;
280  G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
281  G4cout<<"Residual4Momentum "<<exciton4Momentum<<G4endl;
282  G4cout<<"ResidualMass, GroundStateMass and E* "<<ActualMass<<" "<<fMass<<" "
283  <<ActualMass - fMass<<G4endl;
284  #endif
285 
286  if(ActualMass - fMass < 0.)
287  {
288  G4double ResE = std::sqrt(exciton4Momentum.vect().mag2() + sqr(fMass+10*MeV));
289  exciton4Momentum.setE(ResE);
290  #ifdef debugPrecoInt
291  G4cout<<"ActualMass - fMass < 0. "<<ActualMass<<" "<<fMass<<" "<<ActualMass - fMass<<G4endl;
292  G4int Uzhi; G4cin>>Uzhi;
293  #endif
294  }
295  }
296  // Need to de-excite the remnant nucleus only if excitation energy > 0.
297  G4Fragment anInitialState(anA, aZ, exciton4Momentum);
298  anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
299  anInitialState.SetNumberOfCharged(numberOfCh);
300  anInitialState.SetNumberOfHoles(numberOfHoles);
301 
302  G4ReactionProductVector * aPrecoResult =
303  theDeExcitation->DeExcite(anInitialState);
304  // fill pre-compound part into the result, and return
305  #ifdef debugPrecoInt
306  G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
307  #endif
308  for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
309  {
310  theTotalResult->push_back(aPrecoResult->operator[](ll));
311  #ifdef debugPrecoInt
312  G4cout<<"Fragment "<<ll<<" "
313  <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
314  <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
315  <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
316  <<aPrecoResult->operator[](ll)->GetDefinition()->GetPDGMass()<<G4endl;
317  #endif
318  }
319  delete aPrecoResult;
320  }
321 
322  return theTotalResult;
323 }
324 
327 {
328  G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
329  << G4endl;
330  G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
331  G4cout << "Please remove from your physics list."<<G4endl;
332  throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
333  return new G4HadFinalState;
334 }
336 {
337  outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
338  << "energy model through the wounded nucleus to precompound de-excition.\n"
339  << "Low energy protons and neutron present among secondaries produced by \n"
340  << "the high energy generator and within the nucleus are captured. The wounded\n"
341  << "nucleus and the captured particles form an excited nuclear fragment. This\n"
342  << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
343  << "Nuclear de-excitation:\n";
344  // preco
345 
346 }
347 
348 // Uzhi Nov. 2012 ------------------------------------------------
351  G4V3DNucleus* theProjectileNucleus)
352 {
353 #ifdef debugPrecoInt
354  G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
355  G4cout<<"Projectile A and Z "<<theProjectileNucleus->GetMassNumber()<<" "
356  <<theProjectileNucleus->GetCharge()<<G4endl;
357  G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "
358  <<theNucleus->GetCharge()<<G4endl;
359  G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
360  G4cout<<"Projectile 4Mom and mass "<<GetPrimaryProjectile()->Get4Momentum()<<" "
361  <<GetPrimaryProjectile()->Get4Momentum().mag()<<G4endl<<G4endl;
362 #endif
363 
364  // prepare the target residual
365  G4int anA=theNucleus->GetMassNumber();
366  G4int aZ=theNucleus->GetCharge();
367  G4int numberOfEx = 0;
368  G4int numberOfCh = 0;
369  G4int numberOfHoles = 0;
370  G4double exEnergy = 0.0;
371  G4double R = theNucleus->GetNuclearRadius();
372  G4LorentzVector Target4Momentum(0.,0.,0.,0.);
373 
374  // loop over wounded target nucleus
375  G4Nucleon * theCurrentNucleon =
376  theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
377  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
378  {
379  if(theCurrentNucleon->AreYouHit()) {
380  ++numberOfHoles;
381  ++numberOfEx;
382  --anA;
383  aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
384  eplus + 0.1);
385  exEnergy += theCurrentNucleon->GetBindingEnergy();
386  Target4Momentum -=theCurrentNucleon->Get4Momentum();
387  }
388  theCurrentNucleon = theNucleus->GetNextNucleon();
389  }
390 
391 #ifdef debugPrecoInt
392  G4cout<<"Residual Target A Z E* 4mom "<<anA<<" "<<aZ<<" "<<exEnergy<<" "
393  <<Target4Momentum<<G4endl;
394 #endif
395 
396  // prepare the projectile residual
397 
398  G4bool ProjectileIsAntiNucleus=
400 
401  G4ThreeVector bst = GetPrimaryProjectile()->Get4Momentum().boostVector();
402 
403  G4int anAb=theProjectileNucleus->GetMassNumber();
404  G4int aZb=theProjectileNucleus->GetCharge();
405  G4int numberOfExB = 0;
406  G4int numberOfChB = 0;
407  G4int numberOfHolesB = 0;
408  G4double exEnergyB = 0.0;
409  G4double Rb = theProjectileNucleus->GetNuclearRadius();
410  G4LorentzVector Projectile4Momentum(0.,0.,0.,0.);
411 
412  // loop over wounded projectile nucleus
413  theCurrentNucleon =
414  theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
415  while(theCurrentNucleon) /* Loop checking, 31.08.2015, G.Folger */
416  {
417  if(theCurrentNucleon->AreYouHit()) {
418  ++numberOfHolesB;
419  ++numberOfExB;
420  --anAb;
421  if(!ProjectileIsAntiNucleus) {
422  aZb -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
423  eplus + 0.1);
424  } else {
425  aZb += G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
426  eplus - 0.1);
427  }
428  exEnergyB += theCurrentNucleon->GetBindingEnergy();
429  Projectile4Momentum -=theCurrentNucleon->Get4Momentum();
430  }
431  theCurrentNucleon = theProjectileNucleus->GetNextNucleon();
432  }
433 
434  G4bool ExistTargetRemnant = G4double (numberOfHoles) <
435  0.3* G4double (numberOfHoles + anA);
436  G4bool ExistProjectileRemnant= G4double (numberOfHolesB) <
437  0.3*G4double (numberOfHolesB + anAb);
438 
439 #ifdef debugPrecoInt
440  G4cout<<"Projectile residual A Z E* 4mom "<<anAb<<" "<<aZb<<" "<<exEnergyB<<" "
441  <<Projectile4Momentum<<G4endl;
442  G4cout<<" ExistTargetRemnant ExistProjectileRemnant "
443  <<ExistTargetRemnant<<" "<< ExistProjectileRemnant<<G4endl;
444 #endif
445  //-----------------------------------------------------------------------------
446  // decay the strong resonances
447  G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
448  G4DecayKineticTracks decay(theSecondaries);
449  #ifdef debugPrecoInt
450  G4cout<<"Secondary stable particles number "<<theSecondaries->size()<<G4endl;
451  #endif
452 
453 #ifdef debugPrecoInt
454  G4LorentzVector secondary4Momemtum(0,0,0,0);
455  G4int SecondrNum(0);
456 #endif
457 
458  // loop over secondaries
459  G4KineticTrackVector::iterator iter;
460  for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
461  {
462  const G4ParticleDefinition* part = (*iter)->GetDefinition();
463  G4LorentzVector aTrack4Momentum=(*iter)->Get4Momentum();
464 
465  if( part != proton && part != neutron &&
466  (part != ANTIproton && ProjectileIsAntiNucleus) &&
467  (part != ANTIneutron && ProjectileIsAntiNucleus) )
468  {
469  G4ReactionProduct * theNew = new G4ReactionProduct(part);
470  theNew->SetMomentum(aTrack4Momentum.vect());
471  theNew->SetTotalEnergy(aTrack4Momentum.e());
472  theTotalResult->push_back(theNew);
473 #ifdef debugPrecoInt
474  SecondrNum++;
475  secondary4Momemtum += (*iter)->Get4Momentum();
476  G4cout<<"Secondary "<<SecondrNum<<" "
477  <<theNew->GetDefinition()->GetParticleName()<<" "
478  <<theNew->GetMomentum()<<" "<<theNew->GetTotalEnergy()<<G4endl;
479 
480 #endif
481  delete (*iter);
482  continue;
483  }
484 
485  G4bool CanBeCapturedByTarget = false;
486  if( part == proton || part == neutron)
487  {
488  CanBeCapturedByTarget = ExistTargetRemnant &&
490  (aTrack4Momentum + Target4Momentum).mag() -
491  aTrack4Momentum.mag() - Target4Momentum.mag()) &&
492  ((*iter)->GetPosition().mag() < R);
493  }
494  // ---------------------------
495  G4LorentzVector Position((*iter)->GetPosition(), (*iter)->GetFormationTime());
496  Position.boost(bst);
497 
498  G4bool CanBeCapturedByProjectile = false;
499 
500  if( !ProjectileIsAntiNucleus &&
501  ( part == proton || part == neutron))
502  {
503  CanBeCapturedByProjectile = ExistProjectileRemnant &&
505  (aTrack4Momentum + Projectile4Momentum).mag() -
506  aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
507  (Position.vect().mag() < Rb);
508  }
509 
510  if( ProjectileIsAntiNucleus &&
511  ( part == ANTIproton || part == ANTIneutron))
512  {
513  CanBeCapturedByProjectile = ExistProjectileRemnant &&
515  (aTrack4Momentum + Projectile4Momentum).mag() -
516  aTrack4Momentum.mag() - Projectile4Momentum.mag()) &&
517  (Position.vect().mag() < Rb);
518  }
519 
520  if(CanBeCapturedByTarget && CanBeCapturedByProjectile)
521  {
522  if(G4UniformRand() < 0.5)
523  { CanBeCapturedByTarget = true; CanBeCapturedByProjectile = false;}
524  else
525  { CanBeCapturedByTarget = false; CanBeCapturedByProjectile = true;}
526  }
527 
528  if(CanBeCapturedByTarget)
529  {
530  // within the target nucleus, neutron or proton
531  // now calculate A, Z of the fragment, momentum,
532  // number of exciton states
533 #ifdef debugPrecoInt
534  G4cout<<"Track is CapturedByTarget "<<" "<<part->GetParticleName()<<" "
535  <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
536 #endif
537  ++anA;
538  ++numberOfEx;
539  G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
540  aZ += Z;
541  numberOfCh += Z;
542  Target4Momentum +=aTrack4Momentum;
543  delete (*iter);
544  } else if(CanBeCapturedByProjectile)
545  {
546  // within the projectile nucleus, neutron or proton
547  // now calculate A, Z of the fragment, momentum,
548  // number of exciton states
549 #ifdef debugPrecoInt
550  G4cout<<"Track is CapturedByProjectile"<<" "<<part->GetParticleName()<<" "
551  <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
552 #endif
553  ++anAb;
554  ++numberOfExB;
555  G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
556  if( ProjectileIsAntiNucleus ) Z=-Z;
557  aZb += Z;
558  numberOfChB += Z;
559  Projectile4Momentum +=aTrack4Momentum;
560  delete (*iter);
561  } else
562  { // the track is not captured
563  G4ReactionProduct * theNew = new G4ReactionProduct(part);
564  theNew->SetMomentum(aTrack4Momentum.vect());
565  theNew->SetTotalEnergy(aTrack4Momentum.e());
566  theTotalResult->push_back(theNew);
567 
568 #ifdef debugPrecoInt
569  SecondrNum++;
570  secondary4Momemtum += (*iter)->Get4Momentum();
571 /*
572  G4cout<<"Secondary "<<SecondrNum<<" "
573  <<theNew->GetDefinition()->GetParticleName()<<" "
574  <<secondary4Momemtum<<G4endl;
575 */
576 #endif
577  delete (*iter);
578  continue;
579  }
580  }
581  delete theSecondaries;
582  //-----------------------------------------------------
583 
584  #ifdef debugPrecoInt
585  G4cout<<"Final target residual A Z E* 4mom "<<anA<<" "<<aZ<<" "
586  <<exEnergy<<" "<<Target4Momentum<<G4endl;
587  #endif
588 
589  if(0!=anA )
590  {
592 
593  if((anA == theNucleus->GetMassNumber()) && (exEnergy <= 0.))
594  {Target4Momentum.setE(fMass);}
595 
596  G4double RemnMass=Target4Momentum.mag();
597 
598  if(RemnMass < fMass)
599  {
600  RemnMass=fMass + exEnergy;
601  Target4Momentum.setE(std::sqrt(Target4Momentum.vect().mag2() +
602  RemnMass*RemnMass));
603  } else
604  { exEnergy=RemnMass-fMass;}
605 
606  if( exEnergy < 0.) exEnergy=0.;
607 
608  // Need to de-excite the remnant nucleus
609  G4Fragment anInitialState(anA, aZ, Target4Momentum);
610  anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
611  anInitialState.SetNumberOfCharged(numberOfCh);
612  anInitialState.SetNumberOfHoles(numberOfHoles);
613 
614  G4ReactionProductVector * aPrecoResult =
615  theDeExcitation->DeExcite(anInitialState);
616 
617  #ifdef debugPrecoInt
618  G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
619  #endif
620 
621  // fill pre-compound part into the result, and return
622  for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
623  {
624  theTotalResult->push_back(aPrecoResult->operator[](ll));
625  #ifdef debugPrecoInt
626  G4cout<<"Target fragment "<<ll<<" "
627  <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
628  <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
629  <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
630  <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
631  #endif
632  }
633  delete aPrecoResult;
634  }
635 
636  //-----------------------------------------------------
637  if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
638  {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}
639 
640  #ifdef debugPrecoInt
641  G4cout<<"Final projectile residual A Z E* Pmom Pmag2 "<<anAb<<" "<<aZb<<" "
642  <<exEnergyB<<" "<<Projectile4Momentum<<" "
643  <<Projectile4Momentum.mag2()<<G4endl;
644  #endif
645 
646  if(0!=anAb)
647  {
648  // G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM(); // Uzhi Apr. 2015
649  // Projectile4Momentum.boost(bstToCM); // Uzhi Apr. 2015
650 
651  G4double fMass = G4NucleiProperties::GetNuclearMass(anAb, aZb);
652  G4double RemnMass=Projectile4Momentum.mag();
653 
654  if(RemnMass < fMass)
655  {
656  RemnMass=fMass + exEnergyB;
657  Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() + // Uzhi 8.05.2015
658  RemnMass*RemnMass)); // Uzhi 8.05.2015
659  } else
660  { exEnergyB=RemnMass-fMass;}
661 
662  if( exEnergyB < 0.) exEnergyB=0.;
663 
664  G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM(); // Uzhi Apr. 2015
665  Projectile4Momentum.boost(bstToCM); // Uzhi Apr. 2015
666 
667  // Need to de-excite the remnant nucleus
668  G4Fragment anInitialState(anAb, aZb, Projectile4Momentum);
669  anInitialState.SetNumberOfParticles(numberOfExB-numberOfHolesB);
670  anInitialState.SetNumberOfCharged(numberOfChB);
671  anInitialState.SetNumberOfHoles(numberOfHolesB);
672 
673  G4ReactionProductVector * aPrecoResult =
674  theDeExcitation->DeExcite(anInitialState);
675 
676  #ifdef debugPrecoInt
677  G4cout<<"Projectile fragment number "<<aPrecoResult->size()<<G4endl;
678  #endif
679 
680  // fill pre-compound part into the result, and return
681  for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
682  {
683  G4LorentzVector tmp=G4LorentzVector(aPrecoResult->operator[](ll)->GetMomentum(), // Uzhi 2015
684  aPrecoResult->operator[](ll)->GetTotalEnergy());// Uzhi 2015
685  tmp.boost(-bstToCM); // Transformation to the system of original remnant // Uzhi 2015
686  aPrecoResult->operator[](ll)->SetMomentum(tmp.vect()); // Uzhi 2015
687  aPrecoResult->operator[](ll)->SetTotalEnergy(tmp.e()); // Uzhi 2015
688 
689  if(ProjectileIsAntiNucleus)
690  {
691  const G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
692  const G4ParticleDefinition * LastFragment=aFragment;
693  if (aFragment == proton) {LastFragment=G4AntiProton::AntiProtonDefinition();}
694  else if(aFragment == neutron) {LastFragment=G4AntiNeutron::AntiNeutronDefinition();}
695  else if(aFragment == deuteron){LastFragment=G4AntiDeuteron::AntiDeuteronDefinition();}
696  else if(aFragment == triton) {LastFragment=G4AntiTriton::AntiTritonDefinition();}
697  else if(aFragment == He3) {LastFragment=G4AntiHe3::AntiHe3Definition();}
698  else if(aFragment == He4) {LastFragment=G4AntiAlpha::AntiAlphaDefinition();}
699  else {}
700 
701  aPrecoResult->operator[](ll)->SetDefinitionAndUpdateE(LastFragment);
702  }
703 
704  #ifdef debugPrecoInt
705  G4cout<<"Projectile fragment "<<ll<<" "
706  <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
707  <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
708  <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
709  <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
710  #endif
711 //Uzhi
712  theTotalResult->push_back(aPrecoResult->operator[](ll));
713  }
714 
715  delete aPrecoResult;
716  }
717 
718  return theTotalResult;
719 }
720 
721 
static G4AntiTriton * AntiTritonDefinition()
Definition: G4AntiTriton.cc:89
static G4AntiHe3 * AntiHe3()
Definition: G4AntiHe3.cc:94
static const double MeV
Definition: G4SIunits.hh:211
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4GeneratorPrecompoundInterface(G4VPreCompoundModel *p=0)
virtual G4int GetCharge()=0
CLHEP::Hep3Vector G4ThreeVector
const G4HadProjectile * GetPrimaryProjectile() const
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
virtual G4double GetNuclearRadius()=0
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
virtual G4bool StartLoop()=0
static G4AntiDeuteron * AntiDeuteron()
void SetMomentum(const G4double x, const G4double y, const G4double z)
virtual G4int GetMassNumber()=0
virtual void PropagateModelDescription(std::ostream &) const
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:352
static G4AntiDeuteron * AntiDeuteronDefinition()
virtual const G4LorentzVector & Get4Momentum() const
Definition: G4Nucleon.hh:72
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
static G4AntiAlpha * AntiAlpha()
Definition: G4AntiAlpha.cc:89
static G4AntiProton * AntiProtonDefinition()
Definition: G4AntiProton.cc:88
int G4int
Definition: G4Types.hh:78
static G4AntiNeutron * AntiNeutronDefinition()
const G4String & GetParticleName() const
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
std::vector< G4ReactionProduct * > G4ReactionProductVector
#define G4cin
Definition: G4ios.hh:60
const G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
const G4ParticleDefinition * GetDefinition() const
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
bool G4bool
Definition: G4Types.hh:79
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
void SetTotalEnergy(const G4double en)
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
void SetNumberOfParticles(G4int value)
Definition: G4Fragment.hh:361
virtual G4ReactionProductVector * PropagateNuclNucl(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4V3DNucleus *theProjectileNucleus)
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
static G4AntiHe3 * AntiHe3Definition()
Definition: G4AntiHe3.cc:89
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4HadronicInteraction * FindModel(const G4String &name)
G4double GetTotalEnergy() const
G4double GetPDGMass() const
static G4HadronicInteractionRegistry * Instance()
G4ThreeVector GetMomentum() const
void SetDeExcitation(G4VPreCompoundModel *ptr)
static G4AntiAlpha * AntiAlphaDefinition()
Definition: G4AntiAlpha.cc:84
#define G4endl
Definition: G4ios.hh:61
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
virtual G4Nucleon * GetNextNucleon()=0
void SetNumberOfCharged(G4int value)
Definition: G4Fragment.hh:366
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
static G4AntiTriton * AntiTriton()
Definition: G4AntiTriton.cc:94
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
static G4He3 * He3()
Definition: G4He3.cc:94
G4double GetBindingEnergy() const
Definition: G4Nucleon.hh:75
static G4AntiNeutron * AntiNeutron()
CLHEP::HepLorentzVector G4LorentzVector