Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FTFModel.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 //
27 // $Id: G4FTFModel.cc 100828 2016-11-02 15:25:59Z gcosmo $
28 // GEANT4 tag $Name: $
29 //
30 
31 // ------------------------------------------------------------
32 // GEANT 4 class implementation file
33 //
34 // ---------------- G4FTFModel ----------------
35 // by Gunter Folger, May 1998.
36 // class implementing the excitation in the FTF Parton String Model
37 //
38 // Vladimir Uzhinsky, November - December 2012
39 // simulation of nucleus-nucleus interactions was implemented.
40 // ------------------------------------------------------------
41 
42 #include <utility>
43 
44 #include "G4FTFModel.hh"
45 #include "G4ios.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4FTFParameters.hh"
49 #include "G4FTFParticipants.hh"
50 #include "G4DiffractiveSplitableHadron.hh"
51 #include "G4InteractionContent.hh"
52 #include "G4LorentzRotation.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4ParticleTable.hh"
55 #include "G4IonTable.hh"
56 #include "G4KineticTrack.hh"
57 
58 #include "G4Exp.hh"
59 #include "G4Log.hh"
60 
61 //============================================================================
62 
63 //#define debugFTFmodel
64 //#define debugReggeonCascade
65 //#define debugPutOnMassShell
66 //#define debugAdjust
67 //#define debugBuildString
68 
69 
70 //============================================================================
71 
72 G4FTFModel::G4FTFModel( const G4String& modelName ) :
73  G4VPartonStringModel( modelName ),
74  theExcitation( new G4DiffractiveExcitation() ),
75  theElastic( new G4ElasticHNScattering() ),
76  theAnnihilation( new G4FTFAnnihilation() )
77 {
78  G4VPartonStringModel::SetThisPointer( this );
79  theParameters = 0;
80  NumberOfInvolvedNucleonsOfTarget = 0;
81  NumberOfInvolvedNucleonsOfProjectile= 0;
82  for ( G4int i = 0; i < 250; i++ ) {
83  TheInvolvedNucleonsOfTarget[i] = 0;
84  TheInvolvedNucleonsOfProjectile[i] = 0;
85  }
86 
87  //LowEnergyLimit = 2000.0*MeV;
88  LowEnergyLimit = 1000.0*MeV;
89 
90  HighEnergyInter = true;
91 
92  G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 );
93  ProjectileResidual4Momentum = tmp;
94  ProjectileResidualMassNumber = 0;
95  ProjectileResidualCharge = 0;
96  ProjectileResidualExcitationEnergy = 0.0;
97 
98  TargetResidual4Momentum = tmp;
99  TargetResidualMassNumber = 0;
100  TargetResidualCharge = 0;
101  TargetResidualExcitationEnergy = 0.0;
102 
103  SetEnergyMomentumCheckLevels( 2.0*perCent, 150.0*MeV );
104 }
105 
106 
107 //============================================================================
108 
109 struct DeleteVSplitableHadron { void operator()( G4VSplitableHadron* aH ) { delete aH; } };
110 
111 
112 //============================================================================
113 
114 G4FTFModel::~G4FTFModel() {
115  // Because FTF model can be called for various particles
116  // theParameters must be erased at the end of each call.
117  // Thus the delete is also in G4FTFModel::GetStrings() method.
118  if ( theParameters != 0 ) delete theParameters;
119  if ( theExcitation != 0 ) delete theExcitation;
120  if ( theElastic != 0 ) delete theElastic;
121  if ( theAnnihilation != 0 ) delete theAnnihilation;
122 
123  // Erasing of strings created at annihilation.
124  if ( theAdditionalString.size() != 0 ) {
125  std::for_each( theAdditionalString.begin(), theAdditionalString.end(),
126  DeleteVSplitableHadron() );
127  }
128  theAdditionalString.clear();
129 
130  // Erasing of target involved nucleons.
131  if ( NumberOfInvolvedNucleonsOfTarget != 0 ) {
132  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
133  G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfTarget[i]->GetSplitableHadron();
134  if ( aNucleon ) delete aNucleon;
135  }
136  }
137 
138  // Erasing of projectile involved nucleons.
139  if ( NumberOfInvolvedNucleonsOfProjectile != 0 ) {
140  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
141  G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfProjectile[i]->GetSplitableHadron();
142  if ( aNucleon ) delete aNucleon;
143  }
144  }
145 }
146 
147 
148 //============================================================================
149 
150 void G4FTFModel::Init( const G4Nucleus& aNucleus, const G4DynamicParticle& aProjectile ) {
151 
152  theProjectile = aProjectile;
153 
154  G4double PlabPerParticle( 0.0 ); // Laboratory momentum Pz per particle/nucleon
155 
156  #ifdef debugFTFmodel
157  G4cout << "FTF init Proj Name " << theProjectile.GetDefinition()->GetParticleName() << G4endl
158  << "FTF init Proj Mass " << theProjectile.GetMass()
159  << " " << theProjectile.GetMomentum() << G4endl
160  << "FTF init Proj B Q " << theProjectile.GetDefinition()->GetBaryonNumber()
161  << " " << (G4int) theProjectile.GetDefinition()->GetPDGCharge() << G4endl
162  << "FTF init Target A Z " << aNucleus.GetA_asInt()
163  << " " << aNucleus.GetZ_asInt() << G4endl;
164  #endif
165 
166  theParticipants.Clean();
167 
168  theParticipants.SetProjectileNucleus( 0 );
169 
170  G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 );
171  ProjectileResidualMassNumber = 0;
172  ProjectileResidualCharge = 0;
173  ProjectileResidualExcitationEnergy = 0.0;
174  ProjectileResidual4Momentum = tmp;
175 
176  TargetResidualMassNumber = aNucleus.GetA_asInt();
177  TargetResidualCharge = aNucleus.GetZ_asInt();
178  TargetResidualExcitationEnergy = 0.0;
179  TargetResidual4Momentum = tmp;
180  G4double TargetResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()
181  ->GetIonMass( TargetResidualCharge, TargetResidualMassNumber );
182 
183  TargetResidual4Momentum.setE( TargetResidualMass );
184 
185  if ( std::abs( theProjectile.GetDefinition()->GetBaryonNumber() ) <= 1 ) {
186  // Projectile is a hadron : meson or baryon
187  PlabPerParticle = theProjectile.GetMomentum().z();
188  ProjectileResidualMassNumber = std::abs( theProjectile.GetDefinition()->GetBaryonNumber() );
189  ProjectileResidualCharge = G4int( theProjectile.GetDefinition()->GetPDGCharge() );
190  ProjectileResidualExcitationEnergy = 0.0;
191  //G4double ProjectileResidualMass = theProjectile.GetMass();
192  ProjectileResidual4Momentum.setVect( theProjectile.GetMomentum() );
193  ProjectileResidual4Momentum.setE( theProjectile.GetTotalEnergy() );
194  if ( PlabPerParticle < LowEnergyLimit ) {
195  HighEnergyInter = false;
196  } else {
197  HighEnergyInter = true;
198  }
199  } else {
200  if ( theProjectile.GetDefinition()->GetBaryonNumber() > 1 ) {
201  // Projectile is a nucleus
202  theParticipants.InitProjectileNucleus(theProjectile.GetDefinition()->GetBaryonNumber(),
203  G4int(theProjectile.GetDefinition()->GetPDGCharge()));
204  ProjectileResidualMassNumber = theProjectile.GetDefinition()->GetBaryonNumber();
205  ProjectileResidualCharge = G4int( theProjectile.GetDefinition()->GetPDGCharge() );
206  PlabPerParticle = theProjectile.GetMomentum().z() /
207  theProjectile.GetDefinition()->GetBaryonNumber();
208  if ( PlabPerParticle < LowEnergyLimit ) {
209  HighEnergyInter = false;
210  } else {
211  HighEnergyInter = true;
212  }
213  } else if ( theProjectile.GetDefinition()->GetBaryonNumber() < -1 ) {
214  // Projectile is an anti-nucleus
215  theParticipants.InitProjectileNucleus(
216  std::abs( theProjectile.GetDefinition()->GetBaryonNumber() ),
217  std::abs( G4int( theProjectile.GetDefinition()->GetPDGCharge() ) ) );
218  theParticipants.theProjectileNucleus->StartLoop();
219  G4Nucleon* aNucleon;
220  while ( ( aNucleon = theParticipants.theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
221  if ( aNucleon->GetDefinition() == G4Proton::Proton() ) {
222  aNucleon->SetParticleType( G4AntiProton::AntiProton() );
223  } else if ( aNucleon->GetDefinition() == G4Neutron::Neutron() ) {
224  aNucleon->SetParticleType( G4AntiNeutron::AntiNeutron() );
225  }
226  }
227  ProjectileResidualMassNumber = std::abs( theProjectile.GetDefinition()->GetBaryonNumber() );
228  ProjectileResidualCharge = std::abs( G4int(theProjectile.GetDefinition()->GetPDGCharge()) );
229  PlabPerParticle = theProjectile.GetMomentum().z() /
230  std::abs( theProjectile.GetDefinition()->GetBaryonNumber() );
231  if ( PlabPerParticle < LowEnergyLimit ) {
232  HighEnergyInter = false;
233  } else {
234  HighEnergyInter = true;
235  }
236  }
237 
238  G4ThreeVector BoostVector = theProjectile.GetMomentum() / theProjectile.GetTotalEnergy();
239  theParticipants.theProjectileNucleus->DoLorentzBoost( BoostVector );
240  theParticipants.theProjectileNucleus->DoLorentzContraction( BoostVector );
241  ProjectileResidualExcitationEnergy = 0.0;
242  //G4double ProjectileResidualMass = theProjectile.GetMass();
243  ProjectileResidual4Momentum.setVect( theProjectile.GetMomentum() );
244  ProjectileResidual4Momentum.setE( theProjectile.GetTotalEnergy() );
245  }
246 
247  // Init target nucleus
248  theParticipants.Init( aNucleus.GetA_asInt(), aNucleus.GetZ_asInt() );
249  //theParticipants.Init( aNucleus.GetA_asInt(), 0 ); // For h+neutron
250 
251  if ( theParameters != 0 ) delete theParameters;
252  theParameters = new G4FTFParameters( theProjectile.GetDefinition(), aNucleus.GetA_asInt(),
253  aNucleus.GetZ_asInt(), PlabPerParticle );
254 
255  if ( theAdditionalString.size() != 0 ) {
256  std::for_each( theAdditionalString.begin(), theAdditionalString.end(),
257  DeleteVSplitableHadron() );
258  }
259  theAdditionalString.clear();
260 
261  #ifdef debugFTFmodel
262  G4cout << "FTF end of Init" << G4endl << G4endl;
263  #endif
264 
265  //if ( std::abs( theProjectile.GetDefinition()->GetBaryonNumber() ) <= 1 &&
266  // aNucleus.GetA_asInt() < 2 ) theParameters->SetProbabilityOfElasticScatt( 0.0 );
267 }
268 
269 
270 //============================================================================
271 
272 G4ExcitedStringVector* G4FTFModel::GetStrings() {
273 
274  #ifdef debugFTFmodel
275  G4cout << "G4FTFModel::GetStrings() " << G4endl;
276  #endif
277 
278  G4ExcitedStringVector* theStrings( 0 );
279  theParticipants.GetList( theProjectile, theParameters );
280  StoreInvolvedNucleon();
281 
282  G4bool Success( true );
283 
284  if ( HighEnergyInter ) {
285  ReggeonCascade();
286 
287  #ifdef debugFTFmodel
288  G4cout << "FTF PutOnMassShell " << G4endl;
289  #endif
290 
291  Success = PutOnMassShell();
292 
293  #ifdef debugFTFmodel
294  G4cout << "FTF PutOnMassShell Success? " << Success << G4endl;
295  #endif
296 
297  }
298 
299  #ifdef debugFTFmodel
300  G4cout << "FTF ExciteParticipants " << G4endl;
301  #endif
302 
303  if ( Success ) Success = ExciteParticipants();
304 
305  #ifdef debugFTFmodel
306  G4cout << "FTF ExciteParticipants Success? " << Success << G4endl;
307  #endif
308 
309  if ( Success ) {
310 
311  #ifdef debugFTFmodel
312  G4cout << "FTF BuildStrings ";
313  #endif
314 
315  theStrings = BuildStrings();
316 
317  #ifdef debugFTFmodel
318  G4cout << "FTF BuildStrings " << theStrings << " OK" << G4endl
319  << "FTF GetResiduals of Nuclei " << G4endl;
320  #endif
321 
322  GetResiduals();
323 
324  if ( theParameters != 0 ) {
325  delete theParameters;
326  theParameters = 0;
327  }
328  } else if ( ! GetProjectileNucleus() ) {
329  // Erase the hadron projectile
330  std::vector< G4VSplitableHadron* > primaries;
331  theParticipants.StartLoop();
332  while ( theParticipants.Next() ) { /* Loop checking, 10.08.2015, A.Ribon */
333  const G4InteractionContent& interaction = theParticipants.GetInteraction();
334  // Do not allow for duplicates
335  if ( primaries.end() ==
336  std::find( primaries.begin(), primaries.end(), interaction.GetProjectile() ) ) {
337  primaries.push_back( interaction.GetProjectile() );
338  }
339  }
340  std::for_each( primaries.begin(), primaries.end(), DeleteVSplitableHadron() );
341  primaries.clear();
342  }
343 
344  // Cleaning of the memory
345  G4VSplitableHadron* aNucleon = 0;
346 
347  // Erase the projectile nucleons
348  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
349  aNucleon = TheInvolvedNucleonsOfProjectile[i]->GetSplitableHadron();
350  if ( aNucleon ) delete aNucleon;
351  }
352  NumberOfInvolvedNucleonsOfProjectile = 0;
353 
354  // Erase the target nucleons
355  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
356  aNucleon = TheInvolvedNucleonsOfTarget[i]->GetSplitableHadron();
357  if ( aNucleon ) delete aNucleon;
358  }
359  NumberOfInvolvedNucleonsOfTarget = 0;
360 
361  #ifdef debugFTFmodel
362  G4cout << "End of FTF. Go to fragmentation" << G4endl
363  << "To continue - enter 1, to stop - ^C" << G4endl;
364  //G4int Uzhi; G4cin >> Uzhi;
365  #endif
366 
367  theParticipants.Clean();
368 
369  return theStrings;
370 }
371 
372 
373 //============================================================================
374 
375 void G4FTFModel::StoreInvolvedNucleon() {
376  //To store nucleons involved in the interaction
377 
378  NumberOfInvolvedNucleonsOfTarget = 0;
379 
380  G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
381  theTargetNucleus->StartLoop();
382 
383  G4Nucleon* aNucleon;
384  while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
385  if ( aNucleon->AreYouHit() ) {
386  TheInvolvedNucleonsOfTarget[NumberOfInvolvedNucleonsOfTarget] = aNucleon;
387  NumberOfInvolvedNucleonsOfTarget++;
388  }
389  }
390 
391  #ifdef debugFTFmodel
392  G4cout << "G4FTFModel::StoreInvolvedNucleon -------------" << G4endl;
393  G4cout << "NumberOfInvolvedNucleonsOfTarget " << NumberOfInvolvedNucleonsOfTarget
394  << G4endl << G4endl;
395  #endif
396 
397  if ( ! GetProjectileNucleus() ) return; // The projectile is a hadron
398 
399  // The projectile is a nucleus or an anti-nucleus.
400 
401  NumberOfInvolvedNucleonsOfProjectile = 0;
402 
403  G4V3DNucleus* theProjectileNucleus = GetProjectileNucleus();
404  theProjectileNucleus->StartLoop();
405 
406  G4Nucleon* aProjectileNucleon;
407  while ( ( aProjectileNucleon = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
408  if ( aProjectileNucleon->AreYouHit() ) {
409  // Projectile nucleon was involved in the interaction.
410  TheInvolvedNucleonsOfProjectile[NumberOfInvolvedNucleonsOfProjectile] = aProjectileNucleon;
411  NumberOfInvolvedNucleonsOfProjectile++;
412  }
413  }
414 
415  #ifdef debugFTFmodel
416  G4cout << "NumberOfInvolvedNucleonsOfProjectile " << NumberOfInvolvedNucleonsOfProjectile
417  << G4endl << G4endl;
418  #endif
419  return;
420 }
421 
422 
423 //============================================================================
424 
425 void G4FTFModel::ReggeonCascade() {
426  // Implementation of the reggeon theory inspired model
427 
428  #ifdef debugReggeonCascade
429  G4cout << "G4FTFModel::ReggeonCascade -----------" << G4endl
430  << "theProjectile.GetTotalMomentum() " << theProjectile.GetTotalMomentum() << G4endl
431  << "theProjectile.GetTotalEnergy() " << theProjectile.GetTotalEnergy() << G4endl
432  << "ExcitationE/WN " << theParameters->GetExcitationEnergyPerWoundedNucleon() << G4endl;
433  #endif
434 
435  G4int InitNINt = NumberOfInvolvedNucleonsOfTarget;
436 
437  // Reggeon cascading in target nucleus
438  for ( G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
439  G4Nucleon* aTargetNucleon = TheInvolvedNucleonsOfTarget[ InvTN ];
440 
441  G4double CreationTime = aTargetNucleon->GetSplitableHadron()->GetTimeOfCreation();
442 
443  G4double XofWoundedNucleon = aTargetNucleon->GetPosition().x();
444  G4double YofWoundedNucleon = aTargetNucleon->GetPosition().y();
445 
446  G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
447  theTargetNucleus->StartLoop();
448 
449  G4Nucleon* Neighbour(0);
450  while ( ( Neighbour = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
451  if ( ! Neighbour->AreYouHit() ) {
452  G4double impact2 = sqr( XofWoundedNucleon - Neighbour->GetPosition().x() ) +
453  sqr( YofWoundedNucleon - Neighbour->GetPosition().y() );
454 
455  if ( G4UniformRand() < theParameters->GetCofNuclearDestruction() *
456  G4Exp( -impact2 / theParameters->GetR2ofNuclearDestruction() )
457  ) {
458  // The neighbour nucleon is involved in the reggeon cascade
459  TheInvolvedNucleonsOfTarget[ NumberOfInvolvedNucleonsOfTarget ] = Neighbour;
460  NumberOfInvolvedNucleonsOfTarget++;
461 
462  G4VSplitableHadron* targetSplitable;
463  targetSplitable = new G4DiffractiveSplitableHadron( *Neighbour );
464 
465  Neighbour->Hit( targetSplitable );
466  targetSplitable->SetTimeOfCreation( CreationTime );
467  targetSplitable->SetStatus( 3 ); // 2->3
468  }
469  }
470  }
471  }
472 
473  #ifdef debugReggeonCascade
474  G4cout << "Final NumberOfInvolvedNucleonsOfTarget "
475  << NumberOfInvolvedNucleonsOfTarget << G4endl << G4endl;
476  #endif
477 
478  if ( ! GetProjectileNucleus() ) return;
479 
480  // Nucleus-Nucleus Interaction : Destruction of Projectile
481  G4int InitNINp = NumberOfInvolvedNucleonsOfProjectile;
482 
483  //for ( G4int InvPN = 0; InvPN < NumberOfInvolvedNucleonsOfProjectile; InvPN++ ) {
484  for ( G4int InvPN = 0; InvPN < InitNINp; InvPN++ ) {
485  G4Nucleon* aProjectileNucleon = TheInvolvedNucleonsOfProjectile[ InvPN ];
486 
487  G4double CreationTime = aProjectileNucleon->GetSplitableHadron()->GetTimeOfCreation();
488 
489  G4double XofWoundedNucleon = aProjectileNucleon->GetPosition().x();
490  G4double YofWoundedNucleon = aProjectileNucleon->GetPosition().y();
491 
492  G4V3DNucleus* theProjectileNucleus = GetProjectileNucleus();
493  theProjectileNucleus->StartLoop();
494 
495  G4Nucleon* Neighbour( 0 );
496  while ( ( Neighbour = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
497  if ( ! Neighbour->AreYouHit() ) {
498  G4double impact2= sqr( XofWoundedNucleon - Neighbour->GetPosition().x() ) +
499  sqr( YofWoundedNucleon - Neighbour->GetPosition().y() );
500 
501  if ( G4UniformRand() < theParameters->GetCofNuclearDestructionPr() *
502  G4Exp( -impact2 / theParameters->GetR2ofNuclearDestruction() )
503  ) {
504  // The neighbour nucleon is involved in the reggeon cascade
505  TheInvolvedNucleonsOfProjectile[ NumberOfInvolvedNucleonsOfProjectile ] = Neighbour;
506  NumberOfInvolvedNucleonsOfProjectile++;
507 
508  G4VSplitableHadron* projectileSplitable;
509  projectileSplitable = new G4DiffractiveSplitableHadron( *Neighbour );
510 
511  Neighbour->Hit( projectileSplitable );
512  projectileSplitable->SetTimeOfCreation( CreationTime );
513  projectileSplitable->SetStatus( 3 );
514  }
515  }
516  }
517  }
518 
519  #ifdef debugReggeonCascade
520  G4cout << "NumberOfInvolvedNucleonsOfProjectile "
521  << NumberOfInvolvedNucleonsOfProjectile << G4endl << G4endl;
522  #endif
523 }
524 
525 
526 //============================================================================
527 
528 G4bool G4FTFModel::PutOnMassShell() {
529 
530  G4bool isProjectileNucleus = false;
531  if ( GetProjectileNucleus() ) {
532  isProjectileNucleus = true;
533  }
534 
535  #ifdef debugPutOnMassShell
536  G4cout << "PutOnMassShell start " << G4endl;
537  if ( isProjectileNucleus ) {
538  G4cout << "PutOnMassShell for Nucleus_Nucleus " << G4endl;
539  }
540  #endif
541 
542  G4LorentzVector Pprojectile( theProjectile.GetMomentum(), theProjectile.GetTotalEnergy() );
543  if ( Pprojectile.z() < 0.0 ) {
544  return false;
545  }
546 
547  G4bool isOk = true;
548 
549  G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 );
550  G4LorentzVector PtargetResidual( 0.0, 0.0, 0.0, 0.0 );
551  G4double SumMasses = 0.0;
552  G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
553  G4double TargetResidualMass = 0.0;
554 
555  #ifdef debugPutOnMassShell
556  G4cout << "Target : ";
557  #endif
558  isOk = ComputeNucleusProperties( theTargetNucleus, Ptarget, PtargetResidual, SumMasses,
559  TargetResidualExcitationEnergy, TargetResidualMass,
560  TargetResidualMassNumber, TargetResidualCharge );
561  if ( ! isOk ) return false;
562 
563  G4double Mprojectile = 0.0;
564  G4double M2projectile = 0.0;
565  G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 );
566  G4LorentzVector PprojResidual( 0.0, 0.0, 0.0, 0.0 );
567  G4V3DNucleus* thePrNucleus = GetProjectileNucleus();
568  G4double PrResidualMass = 0.0;
569 
570  if ( ! isProjectileNucleus ) { // hadron-nucleus collision
571  Mprojectile = Pprojectile.mag();
572  M2projectile = Pprojectile.mag2();
573  SumMasses += Mprojectile + 20.0*MeV;
574  } else { // nucleus-nucleus or antinucleus-nucleus collision
575  #ifdef debugPutOnMassShell
576  G4cout << "Projectile : ";
577  #endif
578  isOk = ComputeNucleusProperties( thePrNucleus, Pproj, PprojResidual, SumMasses,
579  ProjectileResidualExcitationEnergy, PrResidualMass,
580  ProjectileResidualMassNumber, ProjectileResidualCharge );
581  if ( ! isOk ) return false;
582  }
583 
584  G4LorentzVector Psum = Pprojectile + Ptarget;
585  G4double SqrtS = Psum.mag();
586  G4double S = Psum.mag2();
587 
588  #ifdef debugPutOnMassShell
589  G4cout << "Psum " << Psum/GeV << " GeV" << G4endl << "SqrtS " << SqrtS/GeV << " GeV" << G4endl
590  << "SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/GeV << " "
591  << PrResidualMass/GeV << " " << TargetResidualMass/GeV << " GeV" << G4endl;
592  #endif
593 
594  if ( SqrtS < SumMasses ) {
595  return false; // It is impossible to simulate after putting nuclear nucleons on mass-shell.
596  }
597 
598  // Try to consider also the excitation energy of the residual nucleus, if this is
599  // possible, with the available energy; otherwise, set the excitation energy to zero.
600  G4double savedSumMasses = SumMasses;
601  if ( isProjectileNucleus ) {
602  SumMasses -= std::sqrt( sqr( PrResidualMass ) + PprojResidual.perp2() );
603  SumMasses += std::sqrt( sqr( PrResidualMass + ProjectileResidualExcitationEnergy )
604  + PprojResidual.perp2() );
605  }
606  SumMasses -= std::sqrt( sqr( TargetResidualMass ) + PtargetResidual.perp2() );
607  SumMasses += std::sqrt( sqr( TargetResidualMass + TargetResidualExcitationEnergy )
608  + PtargetResidual.perp2() );
609 
610  if ( SqrtS < SumMasses ) {
611  SumMasses = savedSumMasses;
612  if ( isProjectileNucleus ) {
613  ProjectileResidualExcitationEnergy = 0.0;
614  }
615  TargetResidualExcitationEnergy = 0.0;
616  }
617 
618  TargetResidualMass += TargetResidualExcitationEnergy;
619  if ( isProjectileNucleus ) {
620  PrResidualMass += ProjectileResidualExcitationEnergy;
621  }
622 
623  #ifdef debugPutOnMassShell
624  if ( isProjectileNucleus ) {
625  G4cout << "PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/GeV << " "
626  << ProjectileResidualExcitationEnergy << " MeV" << G4endl;
627  }
628  G4cout << "TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/GeV << " "
629  << TargetResidualExcitationEnergy << " MeV" << G4endl
630  << "Sum masses " << SumMasses/GeV << G4endl;
631  #endif
632 
633  // Sampling of nucleons what can transfer to delta-isobars
634  if ( isProjectileNucleus && thePrNucleus->GetMassNumber() != 1 ) {
635  isOk = GenerateDeltaIsobar( SqrtS, NumberOfInvolvedNucleonsOfProjectile,
636  TheInvolvedNucleonsOfProjectile, SumMasses );
637  }
638  if ( theTargetNucleus->GetMassNumber() != 1 ) {
639  isOk = isOk &&
640  GenerateDeltaIsobar( SqrtS, NumberOfInvolvedNucleonsOfTarget,
641  TheInvolvedNucleonsOfTarget, SumMasses );
642  }
643  if ( ! isOk ) return false;
644 
645  // Now we know that it is kinematically possible to produce a final state made
646  // of the involved nucleons (or corresponding delta-isobars) and a residual nucleus.
647  // We have to sample the kinematical variables which will allow to define the 4-momenta
648  // of the final state. The sampled kinematical variables refer to the center-of-mass frame.
649  // Notice that the sampling of the transverse momentum corresponds to take into account
650  // Fermi motion.
651 
652  G4LorentzRotation toCms( -1*Psum.boostVector() );
653  G4LorentzVector Ptmp = toCms*Pprojectile;
654  if ( Ptmp.pz() <= 0.0 ) { // "String" moving backwards in c.m.s., abort collision!
655  return false;
656  }
657 
658  G4LorentzRotation toLab( toCms.inverse() );
659 
660  G4double YprojectileNucleus = 0.0;
661  if ( isProjectileNucleus ) {
662  Ptmp = toCms*Pproj;
663  YprojectileNucleus = Ptmp.rapidity();
664  }
665  Ptmp = toCms*Ptarget;
666  G4double YtargetNucleus = Ptmp.rapidity();
667 
668  // Ascribing of the involved nucleons Pt and Xminus
669  G4double DcorP = 0.0;
670  if ( isProjectileNucleus ) {
671  DcorP = theParameters->GetDofNuclearDestruction() / thePrNucleus->GetMassNumber();
672  }
673  G4double DcorT = theParameters->GetDofNuclearDestruction() / theTargetNucleus->GetMassNumber();
674  G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
675  G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
676 
677  #ifdef debugPutOnMassShell
678  if ( isProjectileNucleus ) {
679  G4cout << "Y projectileNucleus " << YprojectileNucleus << G4endl;
680  }
681  G4cout << "Y targetNucleus " << YtargetNucleus << G4endl
682  << "Dcor " << theParameters->GetDofNuclearDestruction()
683  << " DcorP DcorT " << DcorP << " " << DcorT << " AveragePt2 " << AveragePt2 << G4endl;
684  #endif
685 
686  G4double M2proj = M2projectile; // Initialization needed only for hadron-nucleus collisions
687  G4double WplusProjectile = 0.0;
688  G4double M2target = 0.0;
689  G4double WminusTarget = 0.0;
690  G4int NumberOfTries = 0;
691  G4double ScaleFactor = 1.0;
692  G4bool OuterSuccess = true;
693 
694  const G4int maxNumberOfLoops = 1000;
695  G4int loopCounter = 0;
696  do { // while ( ! OuterSuccess )
697  OuterSuccess = true;
698  const G4int maxNumberOfInnerLoops = 10000;
699  do { // while ( SqrtS < Mprojectile + std::sqrt( M2target ) )
700  NumberOfTries++;
701  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
702  // After many tries, it is convenient to reduce the values of DcorP, DcorT and
703  // AveragePt2, so that the sampled momenta (respectively, pz, and pt) of the
704  // involved nucleons (or corresponding delta-isomers) are smaller, and therefore
705  // it is more likely to satisfy the momentum conservation.
706  ScaleFactor /= 2.0;
707  DcorP *= ScaleFactor;
708  DcorT *= ScaleFactor;
709  AveragePt2 *= ScaleFactor;
710  }
711  if ( isProjectileNucleus ) {
712  // Sampling of kinematical properties of projectile nucleons
713  isOk = SamplingNucleonKinematics( AveragePt2, maxPtSquare, DcorP,
714  thePrNucleus, PprojResidual,
715  PrResidualMass, ProjectileResidualMassNumber,
716  NumberOfInvolvedNucleonsOfProjectile,
717  TheInvolvedNucleonsOfProjectile, M2proj );
718  }
719  // Sampling of kinematical properties of target nucleons
720  isOk = isOk &&
721  SamplingNucleonKinematics( AveragePt2, maxPtSquare, DcorT,
722  theTargetNucleus, PtargetResidual,
723  TargetResidualMass, TargetResidualMassNumber,
724  NumberOfInvolvedNucleonsOfTarget,
725  TheInvolvedNucleonsOfTarget, M2target );
726 
727  #ifdef debugPutOnMassShell
728  G4cout << "SqrtS, Mp+Mt, Mp, Mt " << SqrtS/GeV << " "
729  << ( std::sqrt( M2proj ) + std::sqrt( M2target) )/GeV << " "
730  << std::sqrt( M2proj )/GeV << " " << std::sqrt( M2target )/GeV << G4endl;
731  #endif
732 
733  if ( ! isOk ) return false;
734  } while ( ( SqrtS < std::sqrt( M2proj ) + std::sqrt( M2target ) ) &&
735  NumberOfTries < maxNumberOfInnerLoops ); /* Loop checking, 10.08.2015, A.Ribon */
736  if ( NumberOfTries >= maxNumberOfInnerLoops ) {
737  #ifdef debugPutOnMassShell
738  G4cout << "BAD situation: forced exit of the inner while loop!" << G4endl;
739  #endif
740  return false;
741  }
742  if ( isProjectileNucleus ) {
743  isOk = CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus, true,
744  NumberOfInvolvedNucleonsOfProjectile,
745  TheInvolvedNucleonsOfProjectile,
746  WminusTarget, WplusProjectile, OuterSuccess );
747  }
748  isOk = isOk &&
749  CheckKinematics( S, SqrtS, M2proj, M2target, YtargetNucleus, false,
750  NumberOfInvolvedNucleonsOfTarget, TheInvolvedNucleonsOfTarget,
751  WminusTarget, WplusProjectile, OuterSuccess );
752  if ( ! isOk ) return false;
753  } while ( ( ! OuterSuccess ) &&
754  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
755  if ( loopCounter >= maxNumberOfLoops ) {
756  #ifdef debugPutOnMassShell
757  G4cout << "BAD situation: forced exit of the while loop!" << G4endl;
758  #endif
759  return false;
760  }
761 
762  // Now the sampling is completed, and we can determine the kinematics of the
763  // whole system. This is done first in the center-of-mass frame, and then it is boosted
764  // to the lab frame. The transverse momentum of the residual nucleus is determined as
765  // the recoil of each hadron (nucleon or delta) which is emitted, i.e. in such a way
766  // to conserve (by construction) the transverse momentum.
767 
768  if ( ! isProjectileNucleus ) { // hadron-nucleus collision
769 
770  G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
771  G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
772  Pprojectile.setPz( Pzprojectile );
773  Pprojectile.setE( Eprojectile );
774 
775  #ifdef debugPutOnMassShell
776  G4cout << "Proj after in CMS " << Pprojectile << G4endl;
777  #endif
778 
779  Pprojectile.transform( toLab );
780  theProjectile.SetMomentum( Pprojectile.vect() );
781  theProjectile.SetTotalEnergy( Pprojectile.e() );
782 
783  theParticipants.StartLoop();
784  theParticipants.Next();
785  G4VSplitableHadron* primary = theParticipants.GetInteraction().GetProjectile();
786  primary->Set4Momentum( Pprojectile );
787 
788  #ifdef debugPutOnMassShell
789  G4cout << "Final proj. mom in Lab. " << primary->Get4Momentum() << G4endl;
790  #endif
791 
792  } else { // nucleus-nucleus or antinucleus-nucleus collision
793 
794  isOk = FinalizeKinematics( WplusProjectile, true, toLab, PrResidualMass,
795  ProjectileResidualMassNumber, NumberOfInvolvedNucleonsOfProjectile,
796  TheInvolvedNucleonsOfProjectile, ProjectileResidual4Momentum );
797 
798  #ifdef debugPutOnMassShell
799  G4cout << "Projectile Residual4Momentum in CMS " << ProjectileResidual4Momentum << G4endl;
800  #endif
801 
802  if ( ! isOk ) return false;
803 
804  ProjectileResidual4Momentum.transform( toLab );
805 
806  #ifdef debugPutOnMassShell
807  G4cout << "Projectile Residual4Momentum in Lab " << ProjectileResidual4Momentum << G4endl;
808  #endif
809 
810  }
811 
812  isOk = FinalizeKinematics( WminusTarget, false, toLab, TargetResidualMass,
813  TargetResidualMassNumber, NumberOfInvolvedNucleonsOfTarget,
814  TheInvolvedNucleonsOfTarget, TargetResidual4Momentum );
815 
816  #ifdef debugPutOnMassShell
817  G4cout << "Target Residual4Momentum in CMS " << TargetResidual4Momentum << G4endl;
818  #endif
819 
820  if ( ! isOk ) return false;
821 
822  TargetResidual4Momentum.transform( toLab );
823 
824  #ifdef debugPutOnMassShell
825  G4cout << "Target Residual4Momentum in Lab " << TargetResidual4Momentum << G4endl;
826  #endif
827 
828  return true;
829 
830 }
831 
832 
833 //============================================================================
834 
835 G4bool G4FTFModel::ExciteParticipants() {
836 
837  #ifdef debugBuildString
838  G4cout << "G4FTFModel::ExciteParticipants() " << G4endl;
839  #endif
840 
841  G4bool Successfull( true );
842  G4int MaxNumOfInelCollisions = G4int( theParameters->GetMaxNumberOfCollisions() );
843  if ( MaxNumOfInelCollisions > 0 ) { // Plab > Pbound, normal application of FTF is possible
844  G4double ProbMaxNumber = theParameters->GetMaxNumberOfCollisions() - MaxNumOfInelCollisions;
845  if ( G4UniformRand() < ProbMaxNumber ) MaxNumOfInelCollisions++;
846  } else {
847  // Plab < Pbound, normal application of FTF is impossible,low energy corrections applied
848  MaxNumOfInelCollisions = 1;
849  }
850 
851  #ifdef debugBuildString
852  G4cout << "MaxNumOfInelCollisions MaxNumOfInelCollisions " << MaxNumOfInelCollisions << G4endl;
853  #endif
854 
855  G4int CurrentInteraction( 0 );
856  theParticipants.StartLoop();
857 
858  while ( theParticipants.Next() ) { /* Loop checking, 10.08.2015, A.Ribon */
859 
860  CurrentInteraction++;
861  const G4InteractionContent& collision = theParticipants.GetInteraction();
862  G4VSplitableHadron* projectile = collision.GetProjectile();
863  G4Nucleon* ProjectileNucleon = collision.GetProjectileNucleon();
864  G4VSplitableHadron* target = collision.GetTarget();
865  G4Nucleon* TargetNucleon = collision.GetTargetNucleon();
866 
867  #ifdef debugBuildString
868  G4cout << G4endl << "Interaction # Status " << CurrentInteraction << " "
869  << collision.GetStatus() << G4endl << "Pr* Tr* " << projectile << " "
870  << target << G4endl << "projectile->GetStatus target->GetStatus "
871  << projectile->GetStatus() << " " << target->GetStatus() << G4endl
872  << "projectile->GetSoftC target->GetSoftC " << projectile->GetSoftCollisionCount()
873  << " " << target->GetSoftCollisionCount() << G4endl;
874  #endif
875 
876  if ( collision.GetStatus() ) {
877  if ( G4UniformRand() < theParameters->GetProbabilityOfElasticScatt() ) {
878  // Elastic scattering
879 
880  #ifdef debugBuildString
881  G4cout << "Elastic scattering" << G4endl;
882  #endif
883 
884  if ( ! HighEnergyInter ) {
885  G4bool Annihilation = false;
886  G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
887  TargetNucleon, Annihilation );
888  if ( ! Result ) continue;
889  }
890  Successfull = theElastic->ElasticScattering( projectile, target, theParameters )
891  || Successfull;
892  } else if ( G4UniformRand() > theParameters->GetProbabilityOfAnnihilation() ) {
893  // Inelastic scattering
894 
895  #ifdef debugBuildString
896  G4cout << "Inelastic interaction" << G4endl
897  << "MaxNumOfInelCollisions " << MaxNumOfInelCollisions << G4endl;
898  #endif
899 
900  if ( ! HighEnergyInter ) {
901  G4bool Annihilation = false;
902  G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
903  TargetNucleon, Annihilation );
904  if ( ! Result ) continue;
905  }
906  if ( G4UniformRand() <
907  ( 1.0 - target->GetSoftCollisionCount() / MaxNumOfInelCollisions ) *
908  ( 1.0 - projectile->GetSoftCollisionCount() / MaxNumOfInelCollisions ) ) {
909  //if ( ! HighEnergyInter ) {
910  // G4bool Annihilation = false;
911  // G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
912  // TargetNucleon, Annihilation );
913  // if ( ! Result ) continue;
914  //}
915  if (theExcitation->ExciteParticipants( projectile, target, theParameters, theElastic )){
916 
917  #ifdef debugBuildString
918  G4cout << "FTF excitation Successfull " << G4endl;
919  // G4cout << "After pro " << projectile->Get4Momentum() << " "
920  // << projectile->Get4Momentum().mag() << G4endl
921  // << "After tar " << target->Get4Momentum() << " "
922  // << target->Get4Momentum().mag() << G4endl;
923  #endif
924 
925  } else {
926 
927  Successfull = theElastic->ElasticScattering( projectile, target, theParameters )
928  && Successfull;
929 
930  #ifdef debugBuildString
931  G4cout << "FTF excitation Non Successfull -> Elastic scattering "
932  << Successfull << G4endl;
933  #endif
934  }
935  } else { // The inelastic interactition was rejected -> elastic scattering
936 
937  #ifdef debugBuildString
938  G4cout << "Elastic scat. at rejection inelastic scattering" << G4endl;
939  #endif
940 
941  //if ( ! HighEnergyInter ) {
942  // G4bool Annihilation = false;
943  // G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
944  // TargetNucleon, Annihilation );
945  // if ( ! Result) continue;
946  //}
947  Successfull = theElastic->ElasticScattering( projectile, target, theParameters )
948  || Successfull;
949  }
950  } else { // Annihilation
951 
952  #ifdef debugBuildString
953  G4cout << "Annihilation" << G4endl;
954  #endif
955 
956  // Skipping possible interactions of the annihilated nucleons
957  while ( theParticipants.Next() ) { /* Loop checking, 10.08.2015, A.Ribon */
958  G4InteractionContent& acollision = theParticipants.GetInteraction();
959  G4VSplitableHadron* NextProjectileNucleon = acollision.GetProjectile();
960  G4VSplitableHadron* NextTargetNucleon = acollision.GetTarget();
961  if ( projectile == NextProjectileNucleon || target == NextTargetNucleon ) {
962  acollision.SetStatus( 0 );
963  }
964  }
965 
966  // Return to the annihilation
967  theParticipants.StartLoop();
968  for ( G4int I = 0; I < CurrentInteraction; I++ ) theParticipants.Next();
969 
970  // At last, annihilation
971  if ( ! HighEnergyInter ) {
972  G4bool Annihilation = true;
973  G4bool Result = AdjustNucleons( projectile, ProjectileNucleon, target,
974  TargetNucleon, Annihilation );
975  if ( ! Result ) continue;
976  }
977  G4VSplitableHadron* AdditionalString = 0;
978  if ( theAnnihilation->Annihilate( projectile, target, AdditionalString, theParameters ) ){
979  Successfull = Successfull || true;
980 
981  #ifdef debugBuildString
982  G4cout << "Annihilation successfull. " << "*AdditionalString "
983  << AdditionalString << G4endl;
984  //G4cout << "After pro " << projectile->Get4Momentum() << G4endl;
985  //G4cout << "After tar " << target->Get4Momentum() << G4endl;
986  #endif
987 
988  if ( AdditionalString != 0 ) theAdditionalString.push_back( AdditionalString );
989 
990  /*
991  if ( target->GetStatus() == 4 ) {
992  // Skipping possible interactions of the annihilated nucleons
993  while ( theParticipants.Next() ) {
994  G4InteractionContent& acollision = theParticipants.GetInteraction();
995  G4VSplitableHadron* NextProjectileNucleon = acollision.GetProjectile();
996  G4VSplitableHadron* NextTargetNucleon = acollision.GetTarget();
997  if ( target == NextTargetNucleon ) { acollision.SetStatus( 0 ); }
998  }
999  }
1000  theParticipants.StartLoop();
1001  for ( G4int I = 0; I < CurrentInteraction; I++ ) theParticipants.Next();
1002  */
1003 
1004  }
1005  }
1006  }
1007 
1008  #ifdef debugBuildString
1009  G4cout << "----------------------------- Final properties " << G4endl
1010  << "projectile->GetStatus target->GetStatus " << projectile->GetStatus()
1011  << " " << target->GetStatus() << G4endl << "projectile->GetSoftC target->GetSoftC "
1012  << projectile->GetSoftCollisionCount() << " " << target->GetSoftCollisionCount()
1013  << G4endl << "ExciteParticipants() Successfull? " << Successfull << G4endl;
1014  #endif
1015 
1016  } // end of while ( theParticipants.Next() )
1017 
1018  return Successfull;
1019 }
1020 
1021 
1022 //============================================================================
1023 
1024 G4bool G4FTFModel::AdjustNucleons( G4VSplitableHadron* SelectedAntiBaryon,
1025  G4Nucleon* ProjectileNucleon,
1026  G4VSplitableHadron* SelectedTargetNucleon,
1027  G4Nucleon* TargetNucleon,
1028  G4bool Annihilation ) {
1029 
1030  #ifdef debugAdjust
1031  G4cout << "AdjustNucleons ---------------------------------------" << G4endl
1032  << "Proj is nucleus? " << GetProjectileNucleus() << G4endl
1033  << "Proj 4mom " << SelectedAntiBaryon->Get4Momentum() << G4endl
1034  << "Targ 4mom " << SelectedTargetNucleon->Get4Momentum() << G4endl
1035  << "Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
1036  << ProjectileResidualMassNumber << " " << ProjectileResidualCharge << " "
1037  << ProjectileResidualExcitationEnergy << G4endl
1038  << "Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
1039  << TargetResidualMassNumber << " " << TargetResidualCharge << " "
1040  << TargetResidualExcitationEnergy << G4endl
1041  << "Collis. pr tr " << SelectedAntiBaryon->GetSoftCollisionCount()
1042  << SelectedTargetNucleon->GetSoftCollisionCount() << G4endl;
1043  #endif
1044 
1045  if ( SelectedAntiBaryon->GetSoftCollisionCount() != 0 &&
1046  SelectedTargetNucleon->GetSoftCollisionCount() != 0 ) {
1047  return true; // Selected hadrons were adjusted before.
1048  }
1049 
1050  // Ascribing of the involved nucleons Pt and X
1051  G4double Dcor = theParameters->GetDofNuclearDestruction();
1052 
1053  G4double DcorP( 0.0 ), DcorT( 0.0 );
1054  if ( ProjectileResidualMassNumber != 0 ) DcorP = Dcor / G4double(ProjectileResidualMassNumber);
1055  if ( TargetResidualMassNumber != 0 ) DcorT = Dcor / G4double(TargetResidualMassNumber);
1056 
1057  G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
1058  G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
1059  G4double ExcitationEnergyPerWoundedNucleon =
1060  theParameters->GetExcitationEnergyPerWoundedNucleon();
1061 
1062  if ( ( ! GetProjectileNucleus() &&
1063  SelectedAntiBaryon->GetSoftCollisionCount() == 0 &&
1064  SelectedTargetNucleon->GetSoftCollisionCount() == 0 )
1065  ||
1066  ( SelectedAntiBaryon->GetSoftCollisionCount() != 0 &&
1067  SelectedTargetNucleon->GetSoftCollisionCount() == 0 ) ) {
1068  // The case of hadron-nucleus interactions, or
1069  // the case when projectile nuclear nucleon participated in
1070  // a collision, but target nucleon did not participate.
1071 
1072  #ifdef debugAdjust
1073  G4cout << "case 1, hA prcol=0 trcol=0, AA prcol#0 trcol=0" << G4endl;
1074  #endif
1075 
1076  if ( TargetResidualMassNumber < 1 ) {
1077  return false;
1078  }
1079 
1080  if ( SelectedAntiBaryon->Get4Momentum().rapidity() < TargetResidual4Momentum.rapidity() ) {
1081  return false;
1082  }
1083 
1084  if ( TargetResidualMassNumber == 1 ) {
1085  TargetResidualMassNumber = 0;
1086  TargetResidualCharge = 0;
1087  TargetResidualExcitationEnergy = 0.0;
1088  SelectedTargetNucleon->Set4Momentum( TargetResidual4Momentum );
1089  TargetResidual4Momentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1090  return true;
1091  }
1092 
1093  G4LorentzVector Psum = SelectedAntiBaryon->Get4Momentum() + TargetResidual4Momentum;
1094 
1095  #ifdef debugAdjust
1096  G4cout << "Targ res Init " << TargetResidual4Momentum << G4endl;
1097  #endif
1098 
1099  // Transform momenta to cms and then rotate parallel to z axis;
1100  G4LorentzRotation toCms( -1*Psum.boostVector() );
1101  G4LorentzVector Pprojectile = SelectedAntiBaryon->Get4Momentum();
1102  G4LorentzVector Ptmp = toCms * Pprojectile;
1103  toCms.rotateZ( -1*Ptmp.phi() );
1104  toCms.rotateY( -1*Ptmp.theta() );
1105  Pprojectile.transform( toCms );
1106  G4LorentzRotation toLab( toCms.inverse() );
1107 
1108  G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 );
1109 
1110  G4double SqrtS = Psum.mag();
1111  G4double S = sqr( SqrtS );
1112 
1113  G4int TResidualMassNumber = TargetResidualMassNumber - 1;
1114  G4int TResidualCharge = TargetResidualCharge -
1115  G4int( TargetNucleon->GetDefinition()->GetPDGCharge() );
1116 
1117  //AR-30May2016 : brought back residual excitation energy as it was in G4 10.1
1118  G4double TResidualExcitationEnergy = TargetResidualExcitationEnergy +
1119  ExcitationEnergyPerWoundedNucleon;
1120  //AR-30May2016 G4double TResidualExcitationEnergy = TargetResidualExcitationEnergy -
1121  //AR-30May2016 ExcitationEnergyPerWoundedNucleon*G4Log( G4UniformRand());
1122 
1123  if ( TResidualMassNumber <= 1 ) {
1124  TResidualExcitationEnergy = 0.0;
1125  }
1126 
1127  G4double TResidualMass( 0.0 );
1128  if ( TResidualMassNumber != 0 ) {
1129  TResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()
1130  ->GetIonMass( TResidualCharge, TResidualMassNumber );
1131  }
1132 
1133  G4double TNucleonMass = TargetNucleon->GetDefinition()->GetPDGMass();
1134  G4double SumMasses = SelectedAntiBaryon->Get4Momentum().mag() + TNucleonMass + TResidualMass;
1135 
1136  G4bool Stopping = false;
1137 
1138  #ifdef debugAdjust
1139  G4cout << "Annihilation " << Annihilation << G4endl;
1140  #endif
1141 
1142  if ( ! Annihilation ) {
1143  if ( SqrtS < SumMasses ) {
1144  return false;
1145  }
1146  if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1147 
1148  #ifdef debugAdjust
1149  G4cout << "TResidualExcitationEnergy " << TResidualExcitationEnergy << G4endl;
1150  #endif
1151 
1152  TResidualExcitationEnergy = SqrtS - SumMasses;
1153 
1154  #ifdef debugAdjust
1155  G4cout << "TResidualExcitationEnergy " << TResidualExcitationEnergy << G4endl;
1156  #endif
1157 
1158  Stopping = true;
1159  return false;
1160  }
1161  }
1162 
1163  if ( Annihilation ) {
1164 
1165  #ifdef debugAdjust
1166  G4cout << "SqrtS < SumMasses - TNucleonMass " << SqrtS << " "
1167  << SumMasses - TNucleonMass << G4endl;
1168  #endif
1169 
1170  if ( SqrtS < SumMasses - TNucleonMass ) {
1171  return false;
1172  }
1173 
1174  #ifdef debugAdjust
1175  G4cout << "SqrtS < SumMasses " << SqrtS << " " << SumMasses << G4endl;
1176  #endif
1177 
1178  if ( SqrtS < SumMasses ) {
1179  TNucleonMass = SqrtS - (SumMasses - TNucleonMass) - TResidualExcitationEnergy;
1180 
1181  #ifdef debugAdjust
1182  G4cout << "TNucleonMass " << TNucleonMass << G4endl;
1183  #endif
1184 
1185  SumMasses = SqrtS - TResidualExcitationEnergy;
1186  //TResidualExcitationEnergy =0.0;
1187  Stopping = true;
1188  }
1189 
1190  #ifdef debugAdjust
1191  G4cout << "SqrtS < SumMasses " << SqrtS << " " << SumMasses << G4endl;
1192  #endif
1193 
1194  if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1195  TResidualExcitationEnergy = SqrtS - SumMasses;
1196  Stopping = true;
1197  }
1198  }
1199 
1200  #ifdef debugAdjust
1201  G4cout << "Stopping " << Stopping << G4endl;
1202  #endif
1203 
1204  if ( Stopping ) {
1205  // All 3-momenta of particles = 0
1206  // New projectile
1207  Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1208  Ptmp.setE( SelectedAntiBaryon->Get4Momentum().mag() );
1209 
1210  #ifdef debugAdjust
1211  G4cout << "Proj stop " << Ptmp << G4endl;
1212  #endif
1213 
1214  Pprojectile = Ptmp; Pprojectile.transform( toLab );
1215  SelectedAntiBaryon->Set4Momentum( Pprojectile );
1216 
1217  // New target nucleon
1218  Ptmp.setE( TNucleonMass );
1219 
1220  #ifdef debugAdjust
1221  G4cout << "Targ stop " << Ptmp << G4endl;
1222  #endif
1223 
1224  Ptarget = Ptmp; Ptarget.transform( toLab );
1225  SelectedTargetNucleon->Set4Momentum( Ptarget );
1226 
1227  // New target residual
1228  TargetResidualMassNumber = TResidualMassNumber;
1229  TargetResidualCharge = TResidualCharge;
1230  TargetResidualExcitationEnergy = TResidualExcitationEnergy;
1231 
1232  Ptmp.setE( TResidualMass + TargetResidualExcitationEnergy );
1233 
1234  #ifdef debugAdjust
1235  G4cout << "Resi stop " << Ptmp << G4endl;
1236  #endif
1237 
1238  Ptmp.transform( toLab );
1239  TargetResidual4Momentum = Ptmp;
1240 
1241  #ifdef debugAdjust
1242  G4cout << Pprojectile << G4endl << Ptarget << G4endl << TargetResidual4Momentum << G4endl;
1243  #endif
1244 
1245  return true;
1246  }
1247 
1248  G4double Mprojectile = Pprojectile.mag();
1249  G4double M2projectile = Pprojectile.mag2();
1250  G4double WplusProjectile( 0.0 );
1251 
1252  G4LorentzVector TResidual4Momentum = toCms * TargetResidual4Momentum;
1253  G4double YtargetNucleus = TResidual4Momentum.rapidity();
1254 
1255  TResidualMass += TResidualExcitationEnergy;
1256  G4double M2target( 0.0 );
1257  G4double WminusTarget( 0.0 );
1258 
1259  G4ThreeVector PtNucleon( 0.0, 0.0, 0.0 );
1260  G4double XminusNucleon( 0.0 );
1261  G4ThreeVector PtResidual( 0.0, 0.0, 0.0 );
1262  G4double XminusResidual( 0.0 );
1263 
1264  G4int NumberOfTries( 0 );
1265  G4double ScaleFactor( 1.0 );
1266  G4bool OuterSuccess( true );
1267 
1268  const G4int maxNumberOfLoops = 1000;
1269  G4int loopCounter = 0;
1270  do {
1271  OuterSuccess = true;
1272 
1273  const G4int maxNumberOfTries = 10000;
1274  do {
1275 
1276  NumberOfTries++;
1277 
1278  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1279  // At large number of tries it would be better to reduce the values
1280  ScaleFactor /= 2.0;
1281  DcorT *= ScaleFactor;
1282  AveragePt2 *= ScaleFactor;
1283  }
1284 
1285  //if ( TargetResidualMassNumber > 1 ) {
1286  // PtNucleon = GaussianPt( AveragePt2, maxPtSquare );
1287  //} else {
1288  // PtNucleon = G4ThreeVector( 0.0, 0.0, 0.0 );
1289  //}
1290  //PtResidual = -PtNucleon;
1291 
1292  G4bool InerSuccess = true;
1293  if ( TargetResidualMassNumber > 1 ) {
1294  const G4int maxNumberOfInnerLoops = 1000;
1295  G4int innerLoopCounter = 0;
1296  do {
1297  InerSuccess = true;
1298 
1299  PtNucleon = GaussianPt( AveragePt2, maxPtSquare );
1300  PtResidual = -PtNucleon;
1301 
1302  G4double Mtarget = std::sqrt( sqr( TNucleonMass ) + PtNucleon.mag2() ) +
1303  std::sqrt( sqr( TResidualMass ) + PtResidual.mag2() );
1304  if ( SqrtS < Mprojectile + Mtarget ) {
1305  InerSuccess = false;
1306  continue;
1307  }
1308 
1309  G4ThreeVector tmpX = GaussianPt( DcorT*DcorT, 1.0 );
1310  G4double Xcenter = std::sqrt( sqr( TNucleonMass ) + PtNucleon.mag2() ) / Mtarget;
1311  XminusNucleon = Xcenter + tmpX.x();
1312  if ( XminusNucleon <= 0.0 || XminusNucleon >= 1.0 ) {
1313  InerSuccess = false;
1314  continue;
1315  }
1316 
1317  XminusResidual = 1.0 - XminusNucleon;
1318  } while ( ( ! InerSuccess ) &&
1319  ++innerLoopCounter < maxNumberOfInnerLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1320  if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
1321  #ifdef debugAdjust
1322  G4cout << "BAD situation: forced exit of the inner while loop!" << G4endl;
1323  #endif
1324  return false;
1325  }
1326 
1327  } else {
1328  XminusNucleon = 1.0;
1329  XminusResidual = 1.0; // It must be 0, but in the case calculation of Pz, E is problematic
1330  }
1331 
1332  M2target = ( sqr( TNucleonMass ) + PtNucleon.mag2() ) / XminusNucleon +
1333  ( sqr( TResidualMass ) + PtResidual.mag2() ) / XminusResidual;
1334 
1335  } while ( ( SqrtS < Mprojectile + std::sqrt( M2target) ) &&
1336  ++NumberOfTries < maxNumberOfTries ); /* Loop checking, 10.08.2015, A.Ribon */
1337  if ( NumberOfTries >= maxNumberOfTries ) {
1338  #ifdef debugAdjust
1339  G4cout << "BAD situation: forced exit of the intermediate while loop!" << G4endl;
1340  #endif
1341  return false;
1342  }
1343 
1344  G4double DecayMomentum2 = sqr( S ) + sqr( M2projectile ) + sqr( M2target )
1345  - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
1346 
1347  WminusTarget = ( S - M2projectile + M2target + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS;
1348  WplusProjectile = SqrtS - M2target / WminusTarget;
1349 
1350  G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
1351  G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
1352  G4double Yprojectile = 0.5 * G4Log( (Eprojectile + Pzprojectile) /
1353  (Eprojectile - Pzprojectile) );
1354 
1355  #ifdef debugAdjust
1356  G4cout << "DecayMomentum2 " << DecayMomentum2 << G4endl
1357  << "WminusTarget WplusProjectile " << WminusTarget << " " << WplusProjectile
1358  << G4endl << "Yprojectile " << Yprojectile << G4endl;
1359  #endif
1360 
1361  G4double Mt2 = sqr( TNucleonMass ) + PtNucleon.mag2();
1362  G4double Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1363  G4double E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1364  G4double YtargetNucleon = 0.5 * G4Log( (E + Pz)/(E - Pz) );
1365 
1366  #ifdef debugAdjust
1367  G4cout << "YtN Ytr YtN-Ytr " << " " << YtargetNucleon << " " << YtargetNucleus << " "
1368  << YtargetNucleon - YtargetNucleus << G4endl
1369  << "YtN Ypr YtN-Ypr " << " " << YtargetNucleon << " " << Yprojectile
1370  << " " << YtargetNucleon - Yprojectile << G4endl;
1371  #endif
1372 
1373  if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 || Yprojectile < YtargetNucleon ) {
1374  OuterSuccess = false;
1375  continue;
1376  }
1377 
1378  } while ( ( ! OuterSuccess ) &&
1379  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1380  if ( loopCounter >= maxNumberOfLoops ) {
1381  #ifdef debugAdjust
1382  G4cout << "BAD situation: forced exit of the while loop!" << G4endl;
1383  #endif
1384  return false;
1385  }
1386 
1387  G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
1388  G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
1389  Pprojectile.setPz( Pzprojectile ); Pprojectile.setE( Eprojectile );
1390 
1391  #ifdef debugAdjust
1392  G4cout << "Proj after in CMS " << Pprojectile << G4endl;
1393  #endif
1394 
1395  Pprojectile.transform( toLab ); // The work with the projectile is finished at the moment
1396 
1397  SelectedAntiBaryon->Set4Momentum( Pprojectile );
1398 
1399  #ifdef debugAdjust
1400  G4cout << "New proj4M " << Pprojectile << G4endl;
1401  #endif
1402 
1403  G4double Mt2 = sqr( TNucleonMass ) + PtNucleon.mag2();
1404  G4double Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1405  G4double E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
1406 
1407  Ptarget.setPx( PtNucleon.x() ); Ptarget.setPy( PtNucleon.y() );
1408  Ptarget.setPz( Pz ); Ptarget.setE( E );
1409  Ptarget.transform( toLab );
1410  SelectedTargetNucleon->Set4Momentum( Ptarget );
1411 
1412  #ifdef debugAdjust
1413  G4cout << "New targ4M " << Ptarget << G4endl;
1414  #endif
1415 
1416  // New target residual
1417  TargetResidualMassNumber = TResidualMassNumber;
1418  TargetResidualCharge = TResidualCharge;
1419  TargetResidualExcitationEnergy = TResidualExcitationEnergy;
1420 
1421  #ifdef debugAdjust
1422  G4cout << "TargetResidualMassNumber TargetResidualCharge TargetResidualExcitationEnergy "
1423  << TargetResidualMassNumber << " " << TargetResidualCharge << " "
1424  << TargetResidualExcitationEnergy << G4endl;
1425  #endif
1426 
1427  if ( TargetResidualMassNumber != 0 ) {
1428  Mt2 = sqr( TResidualMass ) + PtResidual.mag2();
1429  Pz = -WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
1430  E = WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
1431 
1432  TargetResidual4Momentum.setPx( PtResidual.x() );
1433  TargetResidual4Momentum.setPy( PtResidual.y() );
1434  TargetResidual4Momentum.setPz( Pz );
1435  TargetResidual4Momentum.setE( E );
1436 
1437  #ifdef debugAdjust
1438  G4cout << "New Residu " << TargetResidual4Momentum << " CMS" << G4endl;
1439  #endif
1440 
1441  TargetResidual4Momentum.transform( toLab );
1442 
1443  #ifdef debugAdjust
1444  G4cout << "New Residu " << TargetResidual4Momentum << " Lab" << G4endl;
1445  #endif
1446 
1447  } else {
1448  TargetResidual4Momentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1449  }
1450  return true;
1451 
1452  } else if ( SelectedAntiBaryon->GetSoftCollisionCount() == 0 &&
1453  SelectedTargetNucleon->GetSoftCollisionCount() != 0 ) {
1454  // It is assumed that in the case there is ProjectileResidualNucleus
1455 
1456  #ifdef debugAdjust
1457  G4cout << "case 2, prcol=0 trcol#0" << G4endl;
1458  #endif
1459 
1460  if ( ProjectileResidualMassNumber < 1 ) return false;
1461 
1462  if ( ProjectileResidual4Momentum.rapidity() <=
1463  SelectedTargetNucleon->Get4Momentum().rapidity() ) {
1464  return false;
1465  }
1466 
1467  if ( ProjectileResidualMassNumber == 1 ) {
1468  ProjectileResidualMassNumber = 0;
1469  ProjectileResidualCharge = 0;
1470  ProjectileResidualExcitationEnergy = 0.0;
1471  SelectedAntiBaryon->Set4Momentum( ProjectileResidual4Momentum );
1472  ProjectileResidual4Momentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1473  return true;
1474  }
1475 
1476  G4LorentzVector Psum = ProjectileResidual4Momentum + SelectedTargetNucleon->Get4Momentum();
1477 
1478  // Transform momenta to cms and then rotate parallel to z axis;
1479  G4LorentzRotation toCms( -1*Psum.boostVector() );
1480  G4LorentzVector Pprojectile = ProjectileResidual4Momentum;
1481  G4LorentzVector Ptmp = toCms * Pprojectile;
1482  toCms.rotateZ( -1*Ptmp.phi() );
1483  toCms.rotateY( -1*Ptmp.theta() );
1484  G4LorentzRotation toLab( toCms.inverse() );
1485  G4LorentzVector Ptarget = toCms * SelectedTargetNucleon->Get4Momentum();
1486  Pprojectile.transform( toCms );
1487 
1488  G4double SqrtS = Psum.mag();
1489  G4double S = sqr( SqrtS );
1490 
1491  G4int TResidualMassNumber = ProjectileResidualMassNumber - 1;
1492  G4int TResidualCharge = ProjectileResidualCharge
1493  - std::abs( G4int(ProjectileNucleon->GetDefinition()->GetPDGCharge()) );
1494  //AR-11Oct2016 : brought back residual excitation energy as it was in G4 10.1
1495  G4double TResidualExcitationEnergy = ProjectileResidualExcitationEnergy +
1496  ExcitationEnergyPerWoundedNucleon;
1497  //AR-11Oct2016 G4double TResidualExcitationEnergy = ProjectileResidualExcitationEnergy -
1498  //AR-11Oct2016 ExcitationEnergyPerWoundedNucleon*G4Log( G4UniformRand() );
1499  if ( TResidualMassNumber <= 1 ) {
1500  TResidualExcitationEnergy = 0.0;
1501  }
1502 
1503  G4double TResidualMass( 0.0 );
1504  if ( TResidualMassNumber != 0 ) {
1505  TResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()
1506  ->GetIonMass( TResidualCharge , TResidualMassNumber );
1507  }
1508 
1509  G4double TNucleonMass = ProjectileNucleon->GetDefinition()->GetPDGMass();
1510 
1511  G4double SumMasses = SelectedTargetNucleon->Get4Momentum().mag() +
1512  TNucleonMass + TResidualMass;
1513 
1514  #ifdef debugAdjust
1515  G4cout << "SelectedTN.mag() PNMass + PResidualMass "
1516  << SelectedTargetNucleon->Get4Momentum().mag() << " "
1517  << TNucleonMass << " " << TResidualMass << G4endl;
1518  #endif
1519 
1520  G4bool Stopping = false;
1521 
1522  if ( ! Annihilation ) {
1523  if ( SqrtS < SumMasses ) {
1524  return false;
1525  }
1526  if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1527  TResidualExcitationEnergy = SqrtS - SumMasses;
1528  Stopping = true;
1529  return false;
1530  }
1531  }
1532 
1533  if ( Annihilation ) {
1534  if ( SqrtS < SumMasses - TNucleonMass ) {
1535  return false;
1536  }
1537  if ( SqrtS < SumMasses ) {
1538  TNucleonMass = SqrtS - (SumMasses - TNucleonMass);
1539  SumMasses = SqrtS;
1540  TResidualExcitationEnergy = 0.0;
1541  Stopping = true;
1542  }
1543 
1544  if ( SqrtS < SumMasses + TResidualExcitationEnergy ) {
1545  TResidualExcitationEnergy = SqrtS - SumMasses;
1546  Stopping=true;
1547  }
1548  }
1549 
1550  #ifdef debugAdjust
1551  G4cout << "Stopping " << Stopping << G4endl;
1552  #endif
1553 
1554  if ( Stopping ) {
1555  // All 3-momenta of particles = 0
1556  // New target nucleon
1557  Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1558  Ptmp.setE( SelectedTargetNucleon->Get4Momentum().mag() );
1559  Ptarget = Ptmp; Ptarget.transform( toLab );
1560  SelectedTargetNucleon->Set4Momentum( Ptarget );
1561 
1562  // New projectile nucleon
1563  Ptmp.setE( TNucleonMass );
1564  Pprojectile = Ptmp; Pprojectile.transform( toLab );
1565  SelectedAntiBaryon->Set4Momentum( Pprojectile );
1566 
1567  // New projectile residual
1568  ProjectileResidualMassNumber = TResidualMassNumber;
1569  ProjectileResidualCharge = TResidualCharge;
1570  ProjectileResidualExcitationEnergy = TResidualExcitationEnergy;
1571 
1572  Ptmp.setE( TResidualMass + ProjectileResidualExcitationEnergy );
1573  Ptmp.transform( toLab );
1574  ProjectileResidual4Momentum = Ptmp;
1575 
1576  return true;
1577  }
1578 
1579  G4double Mtarget = Ptarget.mag();
1580  G4double M2target = Ptarget.mag2();
1581 
1582  G4LorentzVector TResidual4Momentum = toCms * ProjectileResidual4Momentum;
1583  G4double YprojectileNucleus = TResidual4Momentum.rapidity();
1584 
1585  TResidualMass += TResidualExcitationEnergy;
1586 
1587  G4double M2projectile( 0.0 );
1588  G4double WminusTarget( 0.0 );
1589  G4double WplusProjectile( 0.0 );
1590  G4ThreeVector PtNucleon( 0.0, 0.0, 0.0 );
1591  G4double XplusNucleon( 0.0 );
1592  G4ThreeVector PtResidual( 0.0, 0.0, 0.0 );
1593  G4double XplusResidual( 0.0 );
1594  G4int NumberOfTries( 0 );
1595  G4double ScaleFactor( 1.0 );
1596  G4bool OuterSuccess( true );
1597 
1598  const G4int maxNumberOfLoops = 1000;
1599  G4int loopCounter = 0;
1600  do { // while ( ! OuterSuccess )
1601 
1602  OuterSuccess = true;
1603  const G4int maxNumberOfTries = 10000;
1604  do { // while ( SqrtS < Mtarget + std::sqrt( M2projectile ) )
1605 
1606  NumberOfTries++;
1607 
1608  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
1609  // At large number of tries it would be better to reduce the values
1610  ScaleFactor /= 2.0;
1611  DcorP *= ScaleFactor;
1612  AveragePt2 *= ScaleFactor;
1613  }
1614 
1615  #ifdef debugAdjust
1616  G4cout << "ProjectileResidualMassNumber " << ProjectileResidualMassNumber << G4endl;
1617  #endif
1618 
1619  if ( ProjectileResidualMassNumber > 1 ) {
1620  PtNucleon = GaussianPt( AveragePt2, maxPtSquare );
1621  } else {
1622  PtNucleon = G4ThreeVector( 0.0, 0.0, 0.0 );
1623  }
1624  PtResidual = -PtNucleon;
1625 
1626  G4double Mprojectile = std::sqrt( sqr( TNucleonMass ) + PtNucleon.mag2() ) +
1627  std::sqrt( sqr( TResidualMass ) + PtResidual.mag2() );
1628 
1629  #ifdef debugAdjust
1630  G4cout << "SqrtS < Mtarget + Mprojectile " << SqrtS << " " << Mtarget
1631  << " " << Mprojectile << " " << Mtarget + Mprojectile << G4endl;
1632  #endif
1633 
1634  M2projectile = sqr( Mprojectile );
1635  if ( SqrtS < Mtarget + Mprojectile ) {
1636  OuterSuccess = false;
1637  continue;
1638  }
1639 
1640  G4double Xcenter = std::sqrt( sqr( TNucleonMass ) + PtNucleon.mag2() ) / Mprojectile;
1641 
1642  G4bool InerSuccess = true;
1643  if ( ProjectileResidualMassNumber > 1 ) {
1644  const G4int maxNumberOfInnerLoops = 1000;
1645  G4int innerLoopCounter = 0;
1646  do {
1647  InerSuccess = true;
1648  G4ThreeVector tmpX = GaussianPt( DcorP*DcorP, 1.0 );
1649  XplusNucleon = Xcenter + tmpX.x();
1650  if ( XplusNucleon <= 0.0 || XplusNucleon >= 1.0 ) {
1651  InerSuccess = false;
1652  continue;
1653  }
1654  XplusResidual = 1.0 - XplusNucleon;
1655  } while ( ( ! InerSuccess ) &&
1656  ++innerLoopCounter < maxNumberOfInnerLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1657  if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
1658  #ifdef debugAdjust
1659  G4cout << "BAD situation: forced exit of the inner while loop!" << G4endl;
1660  #endif
1661  return false;
1662  }
1663 
1664  } else {
1665  XplusNucleon = 1.0;
1666  XplusResidual = 1.0; // It must be 0, but in the case determination
1667  // of Pz and E will be problematic.
1668  }
1669 
1670  #ifdef debugAdjust
1671  G4cout << "TNucleonMass PtNucleon XplusNucleon " << TNucleonMass << " " << PtNucleon
1672  << " " << XplusNucleon << G4endl
1673  << "TResidualMass PtResidual XplusResidual " << TResidualMass << " " << PtResidual
1674  << " " << XplusResidual << G4endl;
1675  #endif
1676 
1677  M2projectile = ( sqr( TNucleonMass ) + PtNucleon.mag2() ) / XplusNucleon +
1678  ( sqr( TResidualMass ) + PtResidual.mag2() ) / XplusResidual;
1679 
1680  #ifdef debugAdjust
1681  G4cout << "SqrtS < Mtarget + std::sqrt(M2projectile) " << SqrtS << " " << Mtarget
1682  << " " << std::sqrt( M2projectile ) << " " << Mtarget + std::sqrt( M2projectile )
1683  << G4endl;
1684  #endif
1685 
1686  } while ( ( SqrtS < Mtarget + std::sqrt( M2projectile ) ) &&
1687  ++NumberOfTries < maxNumberOfTries ); /* Loop checking, 10.08.2015, A.Ribon */
1688  if ( NumberOfTries >= maxNumberOfTries ) {
1689  #ifdef debugAdjust
1690  G4cout << "BAD situation: forced exit of the intermediate while loop!" << G4endl;
1691  #endif
1692  return false;
1693  }
1694 
1695  G4double DecayMomentum2 = sqr( S ) + sqr( M2projectile ) + sqr( M2target )
1696  - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
1697 
1698  WplusProjectile = ( S + M2projectile - M2target + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
1699  WminusTarget = SqrtS - M2projectile/WplusProjectile;
1700 
1701  G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
1702  G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
1703  G4double Ytarget = 0.5 * G4Log( (Etarget + Pztarget)/(Etarget - Pztarget) );
1704 
1705  #ifdef debugAdjust
1706  G4cout << "DecayMomentum2 " << DecayMomentum2 << G4endl
1707  << "WminusTarget WplusProjectile " << WminusTarget << " " << WplusProjectile
1708  << G4endl << "YtargetNucleon " << Ytarget << G4endl;
1709  #endif
1710 
1711  G4double Mt2 = sqr( TNucleonMass ) + PtNucleon.mag2();
1712  G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
1713  G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
1714  G4double YprojectileNucleon = 0.5 * G4Log( (E + Pz)/(E - Pz) );
1715 
1716  #ifdef debugAdjust
1717  G4cout << "YpN Ypr YpN-Ypr " << " " << YprojectileNucleon << " " << YprojectileNucleus
1718  << " " << YprojectileNucleon - YprojectileNucleus << G4endl
1719  << "YpN Ytr YpN-Ytr " << " " << YprojectileNucleon << " " << Ytarget
1720  << " " << YprojectileNucleon - Ytarget << G4endl;
1721  #endif
1722 
1723  if ( std::abs( YprojectileNucleon - YprojectileNucleus ) > 2 ||
1724  Ytarget > YprojectileNucleon ) {
1725  OuterSuccess = false;
1726  continue;
1727  }
1728 
1729  } while ( ( ! OuterSuccess ) &&
1730  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
1731  if ( loopCounter >= maxNumberOfLoops ) {
1732  #ifdef debugAdjust
1733  G4cout << "BAD situation: forced exit of the while loop!" << G4endl;
1734  #endif
1735  return false;
1736  }
1737 
1738  // New target
1739  G4double Pztarget = -WminusTarget/2.0 + M2target/2.0/WminusTarget;
1740  G4double Etarget = WminusTarget/2.0 + M2target/2.0/WminusTarget;
1741  Ptarget.setPz( Pztarget ); Ptarget.setE( Etarget );
1742  Ptarget.transform( toLab ); // The work with the target nucleon is finished at the moment.
1743  SelectedTargetNucleon->Set4Momentum( Ptarget );
1744 
1745  #ifdef debugAdjust
1746  G4cout << "Targ after in Lab " << Ptarget << G4endl;
1747  #endif
1748 
1749  // New projectile
1750  G4double Mt2 = sqr( TNucleonMass ) + PtNucleon.mag2();
1751  G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
1752  G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
1753  Pprojectile.setPx( PtNucleon.x() ); Pprojectile.setPy( PtNucleon.y() );
1754  Pprojectile.setPz( Pz ); Pprojectile.setE( E );
1755  Pprojectile.transform( toLab );
1756  SelectedAntiBaryon->Set4Momentum( Pprojectile );
1757 
1758  #ifdef debugAdjust
1759  G4cout << "Proj after in Lab " << Pprojectile << G4endl;
1760  #endif
1761 
1762  // New projectile residual
1763  ProjectileResidualMassNumber = TResidualMassNumber;
1764  ProjectileResidualCharge = TResidualCharge;
1765  ProjectileResidualExcitationEnergy = TResidualExcitationEnergy;
1766 
1767  if ( ProjectileResidualMassNumber != 0 ) {
1768  Mt2 = sqr( TResidualMass ) + PtResidual.mag2();
1769  Pz = WplusProjectile*XplusResidual/2.0 - Mt2/(2.0*WplusProjectile*XplusResidual);
1770  E = WplusProjectile*XplusResidual/2.0 + Mt2/(2.0*WplusProjectile*XplusResidual);
1771  ProjectileResidual4Momentum.setPx( PtResidual.x() );
1772  ProjectileResidual4Momentum.setPy( PtResidual.y() );
1773  ProjectileResidual4Momentum.setPz( Pz );
1774  ProjectileResidual4Momentum.setE( E );
1775  ProjectileResidual4Momentum.transform( toLab );
1776  } else {
1777  ProjectileResidual4Momentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
1778  }
1779  return true;
1780 
1781  } else { // if ( SelectedAntiBaryon->GetSoftCollisionCount() == 0 &&
1782  // SelectedTargetNucleon->GetSoftCollisionCount() == 0 )
1783 
1784  // It can be in the case of nucleus-nucleus interaction only!
1785 
1786  #ifdef debugAdjust
1787  G4cout << "case 3, prcol=0 trcol=0" << G4endl;
1788  #endif
1789 
1790  if ( ! GetProjectileNucleus() ) return false;
1791 
1792  #ifdef debugAdjust
1793  G4cout << "Proj res Init " << ProjectileResidual4Momentum << G4endl
1794  << "Targ res Init " << TargetResidual4Momentum << G4endl
1795  << "ProjectileResidualMassNumber ProjectileResidualCharge "
1796  << ProjectileResidualMassNumber << " " << ProjectileResidualCharge << G4endl
1797  << "TargetResidualMassNumber TargetResidualCharge " << TargetResidualMassNumber
1798  << " " << TargetResidualCharge << G4endl;
1799  #endif
1800 
1801  G4LorentzVector Psum = ProjectileResidual4Momentum + TargetResidual4Momentum;
1802 
1803  // Transform momenta to cms and then rotate parallel to z axis;
1804  G4LorentzRotation toCms( -1*Psum.boostVector() );
1805  G4LorentzVector Pprojectile = ProjectileResidual4Momentum;
1806  G4LorentzVector Ptmp = toCms * Pprojectile;
1807  toCms.rotateZ( -1*Ptmp.phi() );
1808  toCms.rotateY( -1*Ptmp.theta() );
1809  G4LorentzRotation toLab( toCms.inverse() );
1810  Pprojectile.transform( toCms );
1811  G4LorentzVector Ptarget = toCms * TargetResidual4Momentum;
1812 
1813  G4double SqrtS = Psum.mag();
1814  G4double S = sqr( SqrtS );
1815 
1816  G4int PResidualMassNumber = ProjectileResidualMassNumber - 1;
1817  G4int PResidualCharge = ProjectileResidualCharge -
1818  std::abs( G4int(ProjectileNucleon->GetDefinition()->GetPDGCharge()) );
1819  //AR-11Oct2016 : brought back residual excitation energy as it was in G4 10.1
1820  G4double PResidualExcitationEnergy = ProjectileResidualExcitationEnergy +
1821  ExcitationEnergyPerWoundedNucleon;
1822  //AR-11Oct2016 G4double PResidualExcitationEnergy = ProjectileResidualExcitationEnergy -
1823  //AR-11Oct2016 ExcitationEnergyPerWoundedNucleon*G4Log( G4UniformRand() );
1824  if ( PResidualMassNumber <= 1 ) {
1825  PResidualExcitationEnergy = 0.0;
1826  }
1827 
1828  G4double PResidualMass( 0.0 );
1829  if ( PResidualMassNumber != 0 ) {
1830  PResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()
1831  ->GetIonMass( PResidualCharge, PResidualMassNumber );
1832  }
1833 
1834  G4double PNucleonMass = ProjectileNucleon->GetDefinition()->GetPDGMass();
1835 
1836  G4int TResidualMassNumber = TargetResidualMassNumber - 1;
1837  G4int TResidualCharge = TargetResidualCharge -
1838  G4int( TargetNucleon->GetDefinition()->GetPDGCharge() );
1839  //AR-11Oct2016 : brought back residual excitation energy as it was in G4 10.1
1840  G4double TResidualExcitationEnergy = TargetResidualExcitationEnergy +
1841  ExcitationEnergyPerWoundedNucleon;
1842  //AR-11Oct2016 G4double TResidualExcitationEnergy = TargetResidualExcitationEnergy -
1843  //AR-11Oct2016 ExcitationEnergyPerWoundedNucleon*G4Log( G4UniformRand() );
1844  if ( TResidualMassNumber <= 1 ) {
1845  TResidualExcitationEnergy = 0.0;
1846  }
1847  G4double TResidualMass( 0.0 );
1848  if ( TResidualMassNumber != 0 ) {
1849  TResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()
1850  ->GetIonMass( TResidualCharge, TResidualMassNumber );
1851  }
1852 
1853  G4double TNucleonMass = TargetNucleon->GetDefinition()->GetPDGMass();
1854 
1855  G4double SumMasses = PNucleonMass + PResidualMass + TNucleonMass + TResidualMass;
1856 
1857  #ifdef debugAdjust
1858  G4cout << "PNucleonMass PResidualMass TNucleonMass TResidualMass " << PNucleonMass
1859  << " " << PResidualMass << " " << TNucleonMass << " " << TResidualMass << G4endl
1860  << "PResidualExcitationEnergy " << PResidualExcitationEnergy << G4endl
1861  << "TResidualExcitationEnergy " << TResidualExcitationEnergy << G4endl;
1862  #endif
1863 
1864  G4bool Stopping = false;
1865 
1866  if ( ! Annihilation ) {
1867 
1868  #ifdef debugAdjust
1869  G4cout << "SqrtS < SumMasses " << SqrtS << " " << SumMasses << G4endl;
1870  #endif
1871 
1872  if ( SqrtS < SumMasses ) {
1873  return false;
1874  }
1875 
1876  #ifdef debugAdjust
1877  G4cout << "SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy "
1878  << SqrtS << " " << SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy
1879  << G4endl;
1880  #endif
1881 
1882  if ( SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy ) {
1883  Stopping = true;
1884  if ( PResidualExcitationEnergy <= 0.0 ) {
1885  TResidualExcitationEnergy = SqrtS - SumMasses;
1886  } else if ( TResidualExcitationEnergy <= 0.0 ) {
1887  PResidualExcitationEnergy = SqrtS - SumMasses;
1888  } else {
1889  G4double Fraction = (SqrtS - SumMasses) /
1890  (PResidualExcitationEnergy + TResidualExcitationEnergy);
1891  PResidualExcitationEnergy *= Fraction;
1892  TResidualExcitationEnergy *= Fraction;
1893  }
1894  }
1895  }
1896 
1897  #ifdef debugAdjust
1898  G4cout << "Stopping " << Stopping << G4endl;
1899  #endif
1900 
1901  if ( Annihilation ) {
1902  if ( SqrtS < SumMasses - TNucleonMass ) {
1903  return false;
1904  }
1905  if ( SqrtS < SumMasses ) {
1906  Stopping = true;
1907  TNucleonMass = SqrtS - (SumMasses - TNucleonMass);
1908  SumMasses = SqrtS;
1909  TResidualExcitationEnergy = 0.0;
1910  }
1911  if ( SqrtS < SumMasses + PResidualExcitationEnergy + TResidualExcitationEnergy ) {
1912  Stopping = true;
1913  if ( PResidualExcitationEnergy <= 0.0 ) {
1914  TResidualExcitationEnergy = SqrtS - SumMasses;
1915  } else if ( TResidualExcitationEnergy <= 0.0 ) {
1916  PResidualExcitationEnergy = SqrtS - SumMasses;
1917  } else {
1918  G4double Fraction = (SqrtS - SumMasses) /
1919  (PResidualExcitationEnergy + TResidualExcitationEnergy);
1920  PResidualExcitationEnergy *= Fraction;
1921  TResidualExcitationEnergy *= Fraction;
1922  }
1923  }
1924  }
1925 
1926  if ( Stopping ) {
1927  // All 3-momenta of particles = 0
1928  // New projectile
1929  Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1930  Ptmp.setE( PNucleonMass );
1931  Pprojectile = Ptmp; Pprojectile.transform( toLab );
1932  SelectedAntiBaryon->Set4Momentum( Pprojectile );
1933 
1934  // New projectile residual
1935  ProjectileResidualMassNumber = PResidualMassNumber;
1936  ProjectileResidualCharge = PResidualCharge;
1937  ProjectileResidualExcitationEnergy = PResidualExcitationEnergy;
1938 
1939  Ptmp.setE( PResidualMass + ProjectileResidualExcitationEnergy );
1940  Ptmp.transform( toLab );
1941  ProjectileResidual4Momentum = Ptmp;
1942 
1943  // New target nucleon
1944  Ptmp.setPx( 0.0 ); Ptmp.setPy( 0.0 ); Ptmp.setPz( 0.0 );
1945  Ptmp.setE( TNucleonMass );
1946  Ptarget = Ptmp; Ptarget.transform( toLab );
1947  SelectedTargetNucleon->Set4Momentum( Ptarget );
1948 
1949  // New target residual
1950  TargetResidualMassNumber = TResidualMassNumber;
1951  TargetResidualCharge = TResidualCharge;
1952  TargetResidualExcitationEnergy = TResidualExcitationEnergy;
1953 
1954  Ptmp.setE( TResidualMass + TargetResidualExcitationEnergy );
1955  Ptmp.transform( toLab );
1956  TargetResidual4Momentum = Ptmp;
1957 
1958  return true;
1959  }
1960 
1961  G4LorentzVector PResidual4Momentum = toCms * ProjectileResidual4Momentum;
1962  G4double YprojectileNucleus = PResidual4Momentum.rapidity();
1963 
1964  #ifdef debugAdjust
1965  G4cout << "YprojectileNucleus XcenterP " << YprojectileNucleus << G4endl;
1966  #endif
1967 
1968  G4LorentzVector TResidual4Momentum = toCms*TargetResidual4Momentum;
1969  G4double YtargetNucleus = TResidual4Momentum.rapidity();
1970 
1971  PResidualMass += PResidualExcitationEnergy;
1972  TResidualMass += TResidualExcitationEnergy;
1973 
1974  G4double M2projectile( 0.0 );
1975  G4double M2target( 0.0 );
1976  G4double WminusTarget( 0.0 );
1977  G4double WplusProjectile( 0.0 );
1978 
1979  G4ThreeVector PtNucleonP( 0.0, 0.0, 0.0 );
1980  G4double XplusNucleon( 0.0 );
1981  G4ThreeVector PtResidualP( 0.0, 0.0, 0.0 );
1982  G4double XplusResidual( 0.0 );
1983 
1984  G4ThreeVector PtNucleonT( 0.0, 0.0, 0.0 );
1985  G4double XminusNucleon( 0.0 );
1986  G4ThreeVector PtResidualT( 0.0, 0.0, 0.0 );
1987  G4double XminusResidual( 0.0 );
1988 
1989  G4int NumberOfTries( 0 );
1990  G4double ScaleFactor( 1.0 );
1991  G4bool OuterSuccess( true );
1992 
1993  const G4int maxNumberOfLoops = 1000;
1994  G4int loopCounter = 0;
1995  do { // while ( ! OuterSuccess )
1996 
1997  OuterSuccess = true;
1998  const G4int maxNumberOfTries = 10000;
1999  do { // while ( SqrtS < std::sqrt( M2projectile ) + std::sqrt( M2target ) )
2000 
2001  NumberOfTries++;
2002 
2003  if ( NumberOfTries == 100*(NumberOfTries/100) ) {
2004  // At large number of tries it would be better to reduce the values
2005  ScaleFactor /= 2.0;
2006  DcorP *= ScaleFactor;
2007  DcorT *= ScaleFactor;
2008  AveragePt2 *= ScaleFactor;
2009  }
2010 
2011  #ifdef debugAdjust
2012  //G4cout << "NumberOfTries ScaleFactor " << NumberOfTries << " " << ScaleFactor << G4endl;
2013  #endif
2014 
2015  if ( ProjectileResidualMassNumber > 1 ) {
2016  PtNucleonP = GaussianPt( AveragePt2, maxPtSquare );
2017  } else {
2018  PtNucleonP = G4ThreeVector( 0.0, 0.0, 0.0 );
2019  }
2020  PtResidualP = -PtNucleonP;
2021 
2022  if ( TargetResidualMassNumber > 1 ) {
2023  PtNucleonT = GaussianPt( AveragePt2, maxPtSquare );
2024  } else {
2025  PtNucleonT = G4ThreeVector( 0.0, 0.0, 0.0 );
2026  }
2027  PtResidualT = -PtNucleonT;
2028 
2029  G4double Mprojectile = std::sqrt( sqr( PNucleonMass ) + PtNucleonP.mag2() ) +
2030  std::sqrt( sqr( PResidualMass ) + PtResidualP.mag2() );
2031  M2projectile = sqr( Mprojectile );
2032 
2033  G4double Mtarget = std::sqrt( sqr( TNucleonMass ) + PtNucleonT.mag2() ) +
2034  std::sqrt( sqr( TResidualMass ) + PtResidualT.mag2() );
2035  M2target = sqr( Mtarget );
2036 
2037  if ( SqrtS < Mprojectile + Mtarget ) {
2038  OuterSuccess = false;
2039  continue;
2040  }
2041 
2042  G4bool InerSuccess = true;
2043 
2044  if ( ProjectileResidualMassNumber > 1 ) {
2045  const G4int maxNumberOfInnerLoops = 1000;
2046  G4int innerLoopCounter = 0;
2047  do {
2048  InerSuccess = true;
2049  G4ThreeVector tmpX = GaussianPt( DcorP*DcorP, 1.0 );
2050  G4double XcenterP = std::sqrt( sqr( PNucleonMass ) + PtNucleonP.mag2() ) / Mprojectile;
2051  XplusNucleon = XcenterP + tmpX.x();
2052 
2053  #ifdef debugAdjust
2054  //G4cout << "XplusNucleon 1 " << XplusNucleon << G4endl;
2055  //G4int Uzhi; G4cin >> Uzhi;
2056  #endif
2057 
2058  if ( XplusNucleon <= 0.0 || XplusNucleon >= 1.0 ) {
2059  InerSuccess = false;
2060  continue;
2061  }
2062  XplusResidual = 1.0 - XplusNucleon;
2063  } while ( ( ! InerSuccess ) &&
2064  ++innerLoopCounter < maxNumberOfInnerLoops ); /* Loop checking, 10.08.2015, A.Ribon */
2065  if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
2066  #ifdef debugAdjust
2067  G4cout << "BAD situation: forced exit of the first inner while loop!" << G4endl;
2068  #endif
2069  return false;
2070  }
2071 
2072  #ifdef debugAdjust
2073  //G4cout << "XplusNucleon XplusResidual 2 " << XplusNucleon
2074  // << " " << XplusResidual << G4endl;
2075  //{ G4int Uzhi; G4cin >> Uzhi; }
2076  #endif
2077 
2078  } else {
2079  XplusNucleon = 1.0;
2080  XplusResidual = 1.0; // It must be 0
2081  }
2082 
2083  if ( TargetResidualMassNumber > 1 ) {
2084 
2085  const G4int maxNumberOfInnerLoops = 1000;
2086  G4int innerLoopCounter = 0;
2087  do {
2088  InerSuccess = true;
2089 
2090  G4ThreeVector tmpX = GaussianPt( DcorT*DcorT, 1.0 );
2091  G4double XcenterT = std::sqrt( sqr( TNucleonMass ) + PtNucleonT.mag2() ) / Mtarget;
2092  XminusNucleon = XcenterT + tmpX.x();
2093  if ( XminusNucleon <= 0.0 || XminusNucleon >= 1.0 ) {
2094  InerSuccess = false;
2095  continue;
2096  }
2097  XminusResidual = 1.0 - XminusNucleon;
2098  } while ( ( ! InerSuccess ) &&
2099  ++innerLoopCounter < maxNumberOfInnerLoops ); /* Loop checking, 10.08.2015, A.Ribon */
2100  if ( innerLoopCounter >= maxNumberOfInnerLoops ) {
2101  #ifdef debugAdjust
2102  G4cout << "BAD situation: forced exit of the second inner while loop!" << G4endl;
2103  #endif
2104  return false;
2105  }
2106  } else {
2107  XminusNucleon = 1.0;
2108  XminusResidual = 1.0; // It must be 0
2109  }
2110 
2111  #ifdef debugAdjust
2112  G4cout << "PtNucleonP " << PtNucleonP << " " << PtResidualP << G4endl
2113  << "XplusNucleon XplusResidual " << XplusNucleon << " " << XplusResidual << G4endl
2114  << "PtNucleonT " << PtNucleonT << " " << PtResidualT << G4endl
2115  << "XminusNucleon XminusResidual " << XminusNucleon << " " << XminusResidual
2116  << G4endl;
2117  #endif
2118 
2119  M2projectile = ( sqr( PNucleonMass ) + PtNucleonP.mag2() ) / XplusNucleon +
2120  ( sqr( PResidualMass) + PtResidualP.mag2() ) / XplusResidual;
2121  M2target = ( sqr( TNucleonMass ) + PtNucleonT.mag2() ) / XminusNucleon +
2122  ( sqr( TResidualMass ) + PtResidualT.mag2() ) / XminusResidual;
2123 
2124  } while ( ( SqrtS < std::sqrt( M2projectile ) + std::sqrt( M2target ) ) &&
2125  ++NumberOfTries < maxNumberOfTries ); /* Loop checking, 10.08.2015, A.Ribon */
2126  if ( NumberOfTries >= maxNumberOfTries ) {
2127  #ifdef debugAdjust
2128  G4cout << "BAD situation: forced exit of the intermediate while loop!" << G4endl;
2129  #endif
2130  return false;
2131  }
2132 
2133  G4double DecayMomentum2 = sqr( S ) + sqr( M2projectile ) + sqr( M2target )
2134  - 2.0*S*M2projectile - 2.0*S*M2target - 2.0*M2projectile*M2target;
2135 
2136  WplusProjectile = ( S + M2projectile - M2target + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS;
2137  WminusTarget = SqrtS - M2projectile/WplusProjectile;
2138 
2139  G4double Mt2 = sqr( PNucleonMass ) + PtNucleonP.mag2();
2140  G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2141  G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2142  G4double YprojectileNucleon = 0.5 * G4Log( (E + Pz)/(E - Pz) );
2143 
2144  Mt2 = sqr( TNucleonMass ) + PtNucleonT.mag2();
2145  Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2146  E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2147  G4double YtargetNucleon = 0.5 * G4Log( (E + Pz)/(E - Pz) );
2148 
2149  if ( std::abs( YtargetNucleon - YtargetNucleus ) > 2 ||
2150  std::abs( YprojectileNucleon - YprojectileNucleus ) > 2 ||
2151  YprojectileNucleon < YtargetNucleon ) {
2152  OuterSuccess = false;
2153  continue;
2154  }
2155 
2156  } while ( ( ! OuterSuccess ) &&
2157  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
2158  if ( loopCounter >= maxNumberOfLoops ) {
2159  #ifdef debugAdjust
2160  G4cout << "BAD situation: forced exit of the while loop!" << G4endl;
2161  #endif
2162  return false;
2163  }
2164 
2165  #ifdef debugAdjust
2166  G4cout << "PtNucleonP " << PtNucleonP << G4endl;
2167  #endif
2168 
2169  G4double Mt2 = sqr( PNucleonMass ) + PtNucleonP.mag2();
2170  G4double Pz = WplusProjectile*XplusNucleon/2.0 - Mt2/(2.0*WplusProjectile*XplusNucleon);
2171  G4double E = WplusProjectile*XplusNucleon/2.0 + Mt2/(2.0*WplusProjectile*XplusNucleon);
2172 
2173  Pprojectile.setPx( PtNucleonP.x() ); Pprojectile.setPy( PtNucleonP.y() );
2174  Pprojectile.setPz( Pz ); Pprojectile.setE( E );
2175  Pprojectile.transform( toLab );
2176  SelectedAntiBaryon->Set4Momentum( Pprojectile );
2177 
2178  // New projectile residual
2179  ProjectileResidualMassNumber = PResidualMassNumber;
2180  ProjectileResidualCharge = PResidualCharge;
2181  ProjectileResidualExcitationEnergy = PResidualExcitationEnergy;
2182 
2183  #ifdef debugAdjust
2184  G4cout << "PResidualMass PtResidualP " << PResidualMass << " " << PtResidualP << G4endl;
2185  #endif
2186 
2187  if ( ProjectileResidualMassNumber != 0 ) {
2188  Mt2 = sqr( PResidualMass ) + PtResidualP.mag2();
2189  Pz = WplusProjectile*XplusResidual/2.0 - Mt2/(2.0*WplusProjectile*XplusResidual);
2190  E = WplusProjectile*XplusResidual/2.0 + Mt2/(2.0*WplusProjectile*XplusResidual);
2191  ProjectileResidual4Momentum.setPx( PtResidualP.x() );
2192  ProjectileResidual4Momentum.setPy( PtResidualP.y() );
2193  ProjectileResidual4Momentum.setPz( Pz );
2194  ProjectileResidual4Momentum.setE( E );
2195  ProjectileResidual4Momentum.transform( toLab );
2196  } else {
2197  ProjectileResidual4Momentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
2198  }
2199 
2200  #ifdef debugAdjust
2201  G4cout << "Pr N R " << Pprojectile << G4endl << " "
2202  << ProjectileResidual4Momentum << G4endl;
2203  #endif
2204 
2205  Mt2 = sqr( TNucleonMass ) + PtNucleonT.mag2();
2206  Pz = -WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2207  E = WminusTarget*XminusNucleon/2.0 + Mt2/(2.0*WminusTarget*XminusNucleon);
2208 
2209  Ptarget.setPx( PtNucleonT.x() ); Ptarget.setPy( PtNucleonT.y() );
2210  Ptarget.setPz( Pz ); Ptarget.setE( E );
2211  Ptarget.transform( toLab );
2212  SelectedTargetNucleon->Set4Momentum( Ptarget );
2213 
2214  // New target residual
2215  TargetResidualMassNumber = TResidualMassNumber;
2216  TargetResidualCharge = TResidualCharge;
2217  TargetResidualExcitationEnergy = TResidualExcitationEnergy;
2218 
2219  if ( TargetResidualMassNumber != 0 ) {
2220  Mt2 = sqr( TResidualMass ) + PtResidualT.mag2();
2221  Pz = -WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2222  E = WminusTarget*XminusResidual/2.0 + Mt2/(2.0*WminusTarget*XminusResidual);
2223 
2224  TargetResidual4Momentum.setPx( PtResidualT.x() );
2225  TargetResidual4Momentum.setPy( PtResidualT.y() );
2226  TargetResidual4Momentum.setPz( Pz );
2227  TargetResidual4Momentum.setE( E) ;
2228  TargetResidual4Momentum.transform( toLab );
2229  } else {
2230  TargetResidual4Momentum = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
2231  }
2232 
2233  #ifdef debugAdjust
2234  G4cout << "Tr N R " << Ptarget << G4endl << " " << TargetResidual4Momentum << G4endl;
2235  #endif
2236 
2237  return true;
2238 
2239  }
2240 
2241 }
2242 
2243 
2244 //============================================================================
2245 
2246 G4ExcitedStringVector* G4FTFModel::BuildStrings() {
2247  // Loop over all collisions; find all primaries, and all targets
2248  // (targets may be duplicate in the List (to unique G4VSplitableHadrons) ).
2249 
2250  G4ExcitedStringVector* strings = new G4ExcitedStringVector();
2251  G4ExcitedString* FirstString( 0 ); // If there will be a kink,
2252  G4ExcitedString* SecondString( 0 ); // two strings will be produced.
2253 
2254  if ( ! GetProjectileNucleus() ) {
2255 
2256  std::vector< G4VSplitableHadron* > primaries;
2257  theParticipants.StartLoop();
2258  while ( theParticipants.Next() ) { /* Loop checking, 10.08.2015, A.Ribon */
2259  const G4InteractionContent& interaction = theParticipants.GetInteraction();
2260  // do not allow for duplicates ...
2261  if ( interaction.GetStatus() ) {
2262  if ( primaries.end() == std::find( primaries.begin(), primaries.end(),
2263  interaction.GetProjectile() ) ) {
2264  primaries.push_back( interaction.GetProjectile() );
2265  }
2266  }
2267  }
2268 
2269  #ifdef debugBuildString
2270  G4cout << "G4FTFModel::BuildStrings()" << G4endl
2271  << "Number of projectile strings " << primaries.size() << G4endl;
2272  #endif
2273 
2274  for ( unsigned int ahadron = 0; ahadron < primaries.size(); ahadron++ ) {
2275  G4bool isProjectile( true );
2276  //G4cout << "primaries[ ahadron ] " << primaries[ ahadron ] << G4endl;
2277  //if ( primaries[ ahadron ]->GetStatus() <= 1 ) isProjectile = true;
2278  FirstString = 0; SecondString = 0;
2279  if ( primaries[ ahadron ]->GetStatus() <= 1 ) {
2280  theExcitation->CreateStrings( primaries[ ahadron ], isProjectile,
2281  FirstString, SecondString, theParameters );
2282  } else if ( primaries[ ahadron ]->GetStatus() == 2 ) {
2283  G4LorentzVector ParticleMomentum = primaries[ ahadron ]->Get4Momentum();
2284  G4KineticTrack* aTrack = new G4KineticTrack(
2285  primaries[ ahadron ]->GetDefinition(),
2286  primaries[ ahadron ]->GetTimeOfCreation(),
2287  primaries[ ahadron ]->GetPosition(),
2288  ParticleMomentum );
2289  FirstString = new G4ExcitedString( aTrack );
2290  } else {
2291  G4cout << "Something wrong in FTF Model Build String" << G4endl;
2292  }
2293 
2294  if ( FirstString != 0 ) strings->push_back( FirstString );
2295  if ( SecondString != 0 ) strings->push_back( SecondString );
2296 
2297  #ifdef debugBuildString
2298  G4cout << "FirstString & SecondString? " << FirstString << " " << SecondString << G4endl;
2299  if ( FirstString->IsExcited() ) {
2300  G4cout << "Quarks on the FirstString ends " << FirstString->GetRightParton()->GetPDGcode()
2301  << " " << FirstString->GetLeftParton()->GetPDGcode() << G4endl;
2302  } else {
2303  G4cout << "Kinetic track is stored" << G4endl;
2304  }
2305  #endif
2306 
2307  }
2308 
2309  #ifdef debugBuildString
2310  if ( FirstString->IsExcited() ) {
2311  G4cout << "Check 1 string " << strings->operator[](0)->GetRightParton()->GetPDGcode()
2312  << " " << strings->operator[](0)->GetLeftParton()->GetPDGcode() << G4endl << G4endl;
2313  }
2314  #endif
2315 
2316  std::for_each( primaries.begin(), primaries.end(), DeleteVSplitableHadron() );
2317  primaries.clear();
2318 
2319  } else { // Projectile is a nucleus
2320 
2321  #ifdef debugBuildString
2322  G4cout << "Building of projectile-like strings" << G4endl;
2323  #endif
2324 
2325  G4bool isProjectile = true;
2326  for ( G4int ahadron = 0; ahadron < NumberOfInvolvedNucleonsOfProjectile; ahadron++ ) {
2327 
2328  #ifdef debugBuildString
2329  G4cout << "Nucleon #, status, intCount " << ahadron << " "
2330  << TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron()->GetStatus()
2331  << " " << TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron()
2332  ->GetSoftCollisionCount()<<G4endl;
2333  #endif
2334 
2335  G4VSplitableHadron* aProjectile =
2336  TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron();
2337 
2338  #ifdef debugBuildString
2339  G4cout << G4endl << "ahadron aProjectile Status " << ahadron << " " << aProjectile
2340  << " " << aProjectile->GetStatus() << G4endl;
2341  #endif
2342 
2343  if ( aProjectile->GetStatus() == 0 ) { // A nucleon took part in non-diffractive interaction
2344 
2345  #ifdef debugBuildString
2346  G4cout << "Case1 aProjectile->GetStatus() == 0 " << G4endl;
2347  #endif
2348 
2349  FirstString = 0; SecondString = 0;
2350  theExcitation->CreateStrings(
2351  TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2352  isProjectile, FirstString, SecondString, theParameters );
2353  if ( FirstString != 0 ) strings->push_back( FirstString );
2354  if ( SecondString != 0 ) strings->push_back( SecondString );
2355  } else if ( aProjectile->GetStatus() == 1 && aProjectile->GetSoftCollisionCount() != 0 ) {
2356  // Nucleon took part in diffractive interaction
2357 
2358  #ifdef debugBuildString
2359  G4cout << "Case2 aProjectile->GetStatus() !=0 St==1 SoftCol!=0" << G4endl;
2360  #endif
2361 
2362  FirstString = 0; SecondString = 0;
2363  theExcitation->CreateStrings(
2364  TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2365  isProjectile, FirstString, SecondString, theParameters );
2366  if ( FirstString != 0 ) strings->push_back( FirstString );
2367  if ( SecondString != 0 ) strings->push_back( SecondString );
2368  } else if ( aProjectile->GetStatus() == 1 && aProjectile->GetSoftCollisionCount() == 0 &&
2369  HighEnergyInter ) {
2370  // Nucleon was considered as a paricipant of an interaction,
2371  // but the interaction was skipped due to annihilation.
2372  // It is now considered as an involved nucleon at high energies.
2373 
2374  #ifdef debugBuildString
2375  G4cout << "Case3 aProjectile->GetStatus() !=0 St==1 SoftCol==0" << G4endl;
2376  #endif
2377 
2378  FirstString = 0; SecondString = 0;
2379  theExcitation->CreateStrings(
2380  TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2381  isProjectile, FirstString, SecondString, theParameters );
2382  if ( FirstString != 0 ) strings->push_back( FirstString );
2383  if ( SecondString != 0 ) strings->push_back( SecondString );
2384 
2385  #ifdef debugBuildString
2386  G4cout << " Strings are built for nucleon marked for an interaction, but"
2387  << " the interaction was skipped." << G4endl;
2388  #endif
2389 
2390  } else if ( aProjectile->GetStatus() == 2 || aProjectile->GetStatus() == 3 ) {
2391  // Nucleon which was involved in the Reggeon cascading
2392 
2393  #ifdef debugBuildString
2394  G4cout << "Case4 aProjectile->GetStatus() !=0 St==2 " << G4endl;
2395  #endif
2396 
2397  FirstString = 0; SecondString = 0;
2398  theExcitation->CreateStrings(
2399  TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron(),
2400  isProjectile, FirstString, SecondString, theParameters );
2401  if ( FirstString != 0 ) strings->push_back( FirstString );
2402  if ( SecondString != 0 ) strings->push_back( SecondString );
2403 
2404  #ifdef debugBuildString
2405  G4cout << " Strings are build for involved nucleon." << G4endl;
2406  #endif
2407 
2408  } else {
2409 
2410  #ifdef debugBuildString
2411  G4cout << "Case5 " << G4endl;
2412  #endif
2413 
2414  //TheInvolvedNucleonsOfProjectile[ ahadron ]->Hit( 0 );
2415  //G4cout << TheInvolvedNucleonsOfProjectile[ ahadron ]->GetSplitableHadron() << G4endl;
2416 
2417  #ifdef debugBuildString
2418  G4cout << " No string" << G4endl;
2419  #endif
2420 
2421  }
2422  }
2423  }
2424 
2425  #ifdef debugBuildString
2426  G4cout << "Building of target-like strings" << G4endl;
2427  #endif
2428 
2429  G4bool isProjectile = false;
2430  for ( G4int ahadron = 0; ahadron < NumberOfInvolvedNucleonsOfTarget; ahadron++ ) {
2431  G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfTarget[ ahadron ]->GetSplitableHadron();
2432 
2433  #ifdef debugBuildString
2434  G4cout << "Nucleon #, status, intCount " << aNucleon << " " << ahadron << " "
2435  << aNucleon->GetStatus() << " " << aNucleon->GetSoftCollisionCount()<<G4endl;;
2436  #endif
2437 
2438  if ( aNucleon->GetStatus() == 0 ) { // A nucleon took part in non-diffractive interaction
2439  FirstString = 0 ; SecondString = 0;
2440  theExcitation->CreateStrings( aNucleon, isProjectile,
2441  FirstString, SecondString, theParameters );
2442  if ( FirstString != 0 ) strings->push_back( FirstString );
2443  if ( SecondString != 0 ) strings->push_back( SecondString );
2444 
2445  #ifdef debugBuildString
2446  G4cout << " 1 case A string is build" << G4endl;
2447  #endif
2448 
2449  } else if ( aNucleon->GetStatus() == 1 && aNucleon->GetSoftCollisionCount() != 0 ) {
2450  // A nucleon took part in diffractive interaction
2451  FirstString = 0; SecondString = 0;
2452  theExcitation->CreateStrings( aNucleon, isProjectile,
2453  FirstString, SecondString, theParameters );
2454  if ( FirstString != 0 ) strings->push_back( FirstString );
2455  if ( SecondString != 0 ) strings->push_back( SecondString );
2456 
2457  #ifdef debugBuildString
2458  G4cout << " 2 case A string is build, nucleon was excited." << G4endl;
2459  #endif
2460 
2461  } else if ( aNucleon->GetStatus() == 1 && aNucleon->GetSoftCollisionCount() == 0 &&
2462  HighEnergyInter ) {
2463  // A nucleon was considered as a participant but due to annihilation
2464  // its interactions were skipped. It will be considered as involved one
2465  // at high energies.
2466  FirstString = 0; SecondString = 0;
2467  theExcitation->CreateStrings( aNucleon, isProjectile,
2468  FirstString, SecondString, theParameters );
2469 
2470  if ( SecondString == 0 ) {
2471  G4LorentzVector ParticleMomentum = aNucleon->Get4Momentum();
2472  G4KineticTrack* aTrack = new G4KineticTrack(
2473  aNucleon->GetDefinition(),
2474  aNucleon->GetTimeOfCreation(),
2475  FirstString->GetPosition(),
2476  ParticleMomentum );
2477  delete FirstString;
2478  FirstString = new G4ExcitedString( aTrack );
2479  };
2480 
2481  if ( FirstString != 0 ) strings->push_back( FirstString );
2482  if ( SecondString != 0 ) strings->push_back( SecondString );
2483 
2484  #ifdef debugBuildString
2485  G4cout << "3 case A string is build" << G4endl;
2486  #endif
2487 
2488  } else if ( aNucleon->GetStatus() == 1 && aNucleon->GetSoftCollisionCount() == 0 &&
2489  ! HighEnergyInter ) {
2490  // A nucleon was considered as a participant but due to annihilation
2491  // its interactions were skipped. It will be returned to nucleus
2492  // at low energies energies.
2493  aNucleon->SetStatus( 5 ); // 4->5
2494  // ??? delete aNucleon;
2495 
2496  #ifdef debugBuildString
2497  G4cout << "4 case A string is not build" << G4endl;
2498  #endif
2499 
2500  } else if ( aNucleon->GetStatus() == 2 || // A nucleon took part in quark exchange
2501  aNucleon->GetStatus() == 3 ) { // A nucleon was involved in Reggeon cascading
2502  FirstString = 0; SecondString = 0;
2503  theExcitation->CreateStrings( aNucleon, isProjectile,
2504  FirstString, SecondString, theParameters );
2505 
2506  if ( SecondString == 0 ) {
2507  G4LorentzVector ParticleMomentum = aNucleon->Get4Momentum();
2508  G4KineticTrack* aTrack = new G4KineticTrack(
2509  aNucleon->GetDefinition(),
2510  aNucleon->GetTimeOfCreation(),
2511  aNucleon->GetPosition(), //FirstString->GetPosition(),
2512  ParticleMomentum );
2513  delete FirstString;
2514  FirstString = new G4ExcitedString( aTrack );
2515  };
2516 
2517  if ( FirstString != 0 ) strings->push_back( FirstString );
2518  if ( SecondString != 0 ) strings->push_back( SecondString );
2519 
2520  #ifdef debugBuildString
2521  G4cout << "5 case A string is build" << G4endl;
2522  #endif
2523 
2524  } else {
2525 
2526  #ifdef debugBuildString
2527  G4cout << "6 case No string" << G4endl;
2528  #endif
2529 
2530  }
2531  }
2532 
2533  #ifdef debugBuildString
2534  G4cout << G4endl << "theAdditionalString.size() " << theAdditionalString.size()
2535  << G4endl << G4endl;
2536  #endif
2537 
2538  isProjectile = true;
2539  if ( theAdditionalString.size() != 0 ) {
2540  for ( unsigned int ahadron = 0; ahadron < theAdditionalString.size(); ahadron++ ) {
2541  //if ( theAdditionalString[ ahadron ]->GetStatus() <= 1 ) isProjectile = true;
2542  FirstString = 0; SecondString = 0;
2543  theExcitation->CreateStrings( theAdditionalString[ ahadron ], isProjectile,
2544  FirstString, SecondString, theParameters );
2545  if ( FirstString != 0 ) strings->push_back( FirstString );
2546  if ( SecondString != 0 ) strings->push_back( SecondString );
2547  }
2548  }
2549 
2550  //for ( unsigned int ahadron = 0; ahadron < strings->size(); ahadron++ ) {
2551  // G4cout << ahadron << " " << strings->operator[]( ahadron )->GetRightParton()->GetPDGcode()
2552  // << " " << strings->operator[]( ahadron )->GetLeftParton()->GetPDGcode() << G4endl;
2553  //}
2554  //G4cout << "------------------------" << G4endl;
2555 
2556  return strings;
2557 }
2558 
2559 
2560 //============================================================================
2561 
2562 void G4FTFModel::GetResiduals() {
2563  // This method is needed for the correct application of G4PrecompoundModelInterface
2564 
2565  #ifdef debugFTFmodel
2566  G4cout << "GetResiduals(): HighEnergyInter? GetProjectileNucleus()?"
2567  << HighEnergyInter << " " << GetProjectileNucleus() << G4endl;
2568  #endif
2569 
2570  if ( HighEnergyInter ) {
2571 
2572  #ifdef debugFTFmodel
2573  G4cout << "NumberOfInvolvedNucleonsOfTarget "<< NumberOfInvolvedNucleonsOfTarget << G4endl;
2574  #endif
2575 
2576  G4double DeltaExcitationE = TargetResidualExcitationEnergy /
2577  G4double( NumberOfInvolvedNucleonsOfTarget );
2578  G4LorentzVector DeltaPResidualNucleus = TargetResidual4Momentum /
2579  G4double( NumberOfInvolvedNucleonsOfTarget );
2580 
2581  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
2582  G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
2583 
2584  #ifdef debugFTFmodel
2585  G4VSplitableHadron* targetSplitable = aNucleon->GetSplitableHadron();
2586  G4cout << i << " Hit? " << aNucleon->AreYouHit() << " " << targetSplitable << G4endl;
2587  if ( targetSplitable ) G4cout << i << "Status " << targetSplitable->GetStatus() << G4endl;
2588  #endif
2589 
2590  G4LorentzVector tmp = -DeltaPResidualNucleus;
2591  aNucleon->SetMomentum( tmp );
2592  aNucleon->SetBindingEnergy( DeltaExcitationE );
2593  }
2594 
2595  if ( TargetResidualMassNumber != 0 ) {
2596  G4ThreeVector bstToCM = TargetResidual4Momentum.findBoostToCM();
2597 
2598  G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
2599  G4LorentzVector residualMomentum( 0.0, 0.0, 0.0, 0.0 );
2600  G4Nucleon* aNucleon = 0;
2601  theTargetNucleus->StartLoop();
2602  while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2603  if ( ! aNucleon->AreYouHit() ) {
2604  G4LorentzVector tmp = aNucleon->Get4Momentum(); tmp.boost( bstToCM );
2605  aNucleon->SetMomentum( tmp );
2606  residualMomentum += tmp;
2607  }
2608  }
2609 
2610  residualMomentum /= TargetResidualMassNumber;
2611 
2612  G4double Mass = TargetResidual4Momentum.mag();
2613  G4double SumMasses = 0.0;
2614 
2615  aNucleon = 0;
2616  theTargetNucleus->StartLoop();
2617  while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2618  if ( ! aNucleon->AreYouHit() ) {
2619  G4LorentzVector tmp = aNucleon->Get4Momentum() - residualMomentum;
2620  G4double E = std::sqrt( tmp.vect().mag2() +
2621  sqr( aNucleon->GetDefinition()->GetPDGMass() - aNucleon->GetBindingEnergy() ) );
2622  tmp.setE( E ); aNucleon->SetMomentum( tmp );
2623  SumMasses += E;
2624  }
2625  }
2626 
2627  G4double Chigh = Mass / SumMasses; G4double Clow = 0.0; G4double C;
2628  const G4int maxNumberOfLoops = 1000;
2629  G4int loopCounter = 0;
2630  do {
2631  C = ( Chigh + Clow ) / 2.0;
2632  SumMasses = 0.0;
2633  aNucleon = 0;
2634  theTargetNucleus->StartLoop();
2635  while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2636  if ( ! aNucleon->AreYouHit() ) {
2637  G4LorentzVector tmp = aNucleon->Get4Momentum();
2638  G4double E = std::sqrt( tmp.vect().mag2()*sqr(C) +
2639  sqr( aNucleon->GetDefinition()->GetPDGMass() - aNucleon->GetBindingEnergy() ) );
2640  SumMasses += E;
2641  }
2642  }
2643 
2644  if ( SumMasses > Mass ) Chigh = C;
2645  else Clow = C;
2646 
2647  } while ( Chigh - Clow > 0.01 &&
2648  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
2649  if ( loopCounter >= maxNumberOfLoops ) {
2650  #ifdef debugFTFmodel
2651  G4cout << "BAD situation: forced exit of the first while loop in G4FTFModel::GetResidual" << G4endl
2652  << "\t return immediately from the method!" << G4endl;
2653  #endif
2654  return;
2655  }
2656 
2657  aNucleon = 0;
2658  theTargetNucleus->StartLoop();
2659  while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2660  if ( !aNucleon->AreYouHit() ) {
2661  G4LorentzVector tmp = aNucleon->Get4Momentum()*C;
2662  G4double E = std::sqrt( tmp.vect().mag2()+
2663  sqr( aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy() ) );
2664  tmp.setE( E ); tmp.boost( -bstToCM );
2665  aNucleon->SetMomentum( tmp );
2666  }
2667  }
2668  }
2669 
2670  if ( ! GetProjectileNucleus() ) return; // The projectile is a hadron
2671 
2672  #ifdef debugFTFmodel
2673  G4cout << "NumberOfInvolvedNucleonsOfProjectile " << NumberOfInvolvedNucleonsOfProjectile
2674  << G4endl << "ProjectileResidualExcitationEnergy ProjectileResidual4Momentum "
2675  << ProjectileResidualExcitationEnergy << " " << ProjectileResidual4Momentum << G4endl;
2676  #endif
2677 
2678  DeltaExcitationE = ProjectileResidualExcitationEnergy /
2679  G4double( NumberOfInvolvedNucleonsOfProjectile );
2680  DeltaPResidualNucleus = ProjectileResidual4Momentum /
2681  G4double( NumberOfInvolvedNucleonsOfProjectile );
2682 
2683  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
2684  G4Nucleon* aNucleon = TheInvolvedNucleonsOfProjectile[i];
2685 
2686  #ifdef debugFTFmodel
2687  G4VSplitableHadron* projSplitable = aNucleon->GetSplitableHadron();
2688  G4cout << i << " Hit? " << aNucleon->AreYouHit() << " " << projSplitable << G4endl;
2689  if ( projSplitable ) G4cout << i << "Status " << projSplitable->GetStatus() << G4endl;
2690  #endif
2691 
2692  G4LorentzVector tmp = -DeltaPResidualNucleus;
2693  aNucleon->SetMomentum( tmp );
2694  aNucleon->SetBindingEnergy( DeltaExcitationE );
2695  }
2696 
2697  if ( ProjectileResidualMassNumber != 0 ) {
2698  G4ThreeVector bstToCM = ProjectileResidual4Momentum.findBoostToCM();
2699 
2700  G4V3DNucleus* theProjectileNucleus = GetProjectileNucleus();
2701  G4LorentzVector residualMomentum( 0.0, 0.0, 0.0, 0.0);
2702  G4Nucleon* aNucleon = 0;
2703  theProjectileNucleus->StartLoop();
2704  while ( ( aNucleon = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2705  if ( ! aNucleon->AreYouHit() ) {
2706  G4LorentzVector tmp = aNucleon->Get4Momentum(); tmp.boost( bstToCM );
2707  aNucleon->SetMomentum( tmp );
2708  residualMomentum += tmp;
2709  }
2710  }
2711 
2712  residualMomentum /= ProjectileResidualMassNumber;
2713 
2714  G4double Mass = ProjectileResidual4Momentum.mag();
2715  G4double SumMasses= 0.0;
2716 
2717  aNucleon = 0;
2718  theProjectileNucleus->StartLoop();
2719  while ( ( aNucleon = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2720  if ( ! aNucleon->AreYouHit() ) {
2721  G4LorentzVector tmp = aNucleon->Get4Momentum() - residualMomentum;
2722  G4double E=std::sqrt( tmp.vect().mag2() +
2723  sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy() ) );
2724  tmp.setE( E ); aNucleon->SetMomentum( tmp );
2725  SumMasses += E;
2726  }
2727  }
2728 
2729  G4double Chigh = Mass / SumMasses; G4double Clow = 0.0; G4double C;
2730  const G4int maxNumberOfLoops = 1000;
2731  G4int loopCounter = 0;
2732  do {
2733  C = ( Chigh + Clow ) / 2.0;
2734 
2735  SumMasses = 0.0;
2736  aNucleon = 0;
2737  theProjectileNucleus->StartLoop();
2738  while ( ( aNucleon = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2739  if ( ! aNucleon->AreYouHit() ) {
2740  G4LorentzVector tmp = aNucleon->Get4Momentum();
2741  G4double E = std::sqrt( tmp.vect().mag2()*sqr(C) +
2742  sqr( aNucleon->GetDefinition()->GetPDGMass() - aNucleon->GetBindingEnergy() ) );
2743  SumMasses += E;
2744  }
2745  }
2746 
2747  if ( SumMasses > Mass) Chigh = C;
2748  else Clow = C;
2749 
2750  } while ( Chigh - Clow > 0.01 &&
2751  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
2752  if ( loopCounter >= maxNumberOfLoops ) {
2753  #ifdef debugFTFmodel
2754  G4cout << "BAD situation: forced exit of the second while loop in G4FTFModel::GetResidual" << G4endl
2755  << "\t return immediately from the method!" << G4endl;
2756  #endif
2757  return;
2758  }
2759 
2760  aNucleon = 0;
2761  theProjectileNucleus->StartLoop();
2762  while ( ( aNucleon = theProjectileNucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2763  if ( ! aNucleon->AreYouHit() ) {
2764  G4LorentzVector tmp = aNucleon->Get4Momentum()*C;
2765  G4double E = std::sqrt( tmp.vect().mag2() +
2766  sqr( aNucleon->GetDefinition()->GetPDGMass() - aNucleon->GetBindingEnergy() ) );
2767  tmp.setE( E ); tmp.boost( -bstToCM );
2768  aNucleon->SetMomentum( tmp );
2769  }
2770  }
2771  } // End of if ( ProjectileResidualMassNumber != 0 )
2772 
2773  #ifdef debugFTFmodel
2774  G4cout << "End projectile" << G4endl;
2775  #endif
2776 
2777  } else {
2778 
2779  #ifdef debugFTFmodel
2780  G4cout << "Low energy interaction: Target nucleus --------------" << G4endl
2781  << "Tr ResidualMassNumber Tr ResidualCharge Tr ResidualExcitationEnergy "
2782  << TargetResidualMassNumber << " " << TargetResidualCharge << " "
2783  << TargetResidualExcitationEnergy << G4endl;
2784  #endif
2785 
2786  G4int NumberOfTargetParticipant( 0 );
2787  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
2788  G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
2789  G4VSplitableHadron* targetSplitable = aNucleon->GetSplitableHadron();
2790  if ( targetSplitable->GetSoftCollisionCount() != 0 ) NumberOfTargetParticipant++;
2791  }
2792 
2793  G4double DeltaExcitationE( 0.0 );
2794  G4LorentzVector DeltaPResidualNucleus = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
2795 
2796  if ( NumberOfTargetParticipant != 0 ) {
2797  DeltaExcitationE = TargetResidualExcitationEnergy / G4double( NumberOfTargetParticipant );
2798  DeltaPResidualNucleus = TargetResidual4Momentum / G4double( NumberOfTargetParticipant );
2799  }
2800 
2801  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
2802  G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
2803  G4VSplitableHadron* targetSplitable = aNucleon->GetSplitableHadron();
2804  if ( targetSplitable->GetSoftCollisionCount() != 0 ) {
2805  G4LorentzVector tmp = -DeltaPResidualNucleus;
2806  aNucleon->SetMomentum( tmp );
2807  aNucleon->SetBindingEnergy( DeltaExcitationE );
2808  } else {
2809  delete targetSplitable;
2810  targetSplitable = 0;
2811  aNucleon->Hit( targetSplitable );
2812  aNucleon->SetBindingEnergy( 0.0 );
2813  }
2814  }
2815 
2816  #ifdef debugFTFmodel
2817  G4cout << "NumberOfTargetParticipant " << NumberOfTargetParticipant << G4endl
2818  << "TargetResidual4Momentum " << TargetResidual4Momentum << G4endl;
2819  #endif
2820 
2821  if ( ! GetProjectileNucleus() ) return; // The projectile is a hadron
2822 
2823  #ifdef debugFTFmodel
2824  G4cout << "Low energy interaction: Projectile nucleus --------------" << G4endl
2825  << "Pr ResidualMassNumber Pr ResidualCharge Pr ResidualExcitationEnergy "
2826  << ProjectileResidualMassNumber << " " << ProjectileResidualCharge << " "
2827  << ProjectileResidualExcitationEnergy << G4endl;
2828  #endif
2829 
2830  G4int NumberOfProjectileParticipant( 0 );
2831  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
2832  G4Nucleon* aNucleon = TheInvolvedNucleonsOfProjectile[i];
2833  G4VSplitableHadron* projectileSplitable = aNucleon->GetSplitableHadron();
2834  if ( projectileSplitable->GetSoftCollisionCount() != 0 )
2835  NumberOfProjectileParticipant++;
2836  }
2837 
2838  #ifdef debugFTFmodel
2839  G4cout << "NumberOfProjectileParticipant" << G4endl;
2840  #endif
2841 
2842  DeltaExcitationE = 0.0;
2843  DeltaPResidualNucleus = G4LorentzVector( 0.0, 0.0, 0.0, 0.0 );
2844 
2845  if ( NumberOfProjectileParticipant != 0 ) {
2846  DeltaExcitationE = ProjectileResidualExcitationEnergy /
2847  G4double( NumberOfProjectileParticipant );
2848  DeltaPResidualNucleus = ProjectileResidual4Momentum /
2849  G4double( NumberOfProjectileParticipant );
2850  }
2851  //G4cout << "DeltaExcitationE DeltaPResidualNucleus " << DeltaExcitationE
2852  // << " " << DeltaPResidualNucleus << G4endl;
2853  for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
2854  G4Nucleon* aNucleon = TheInvolvedNucleonsOfProjectile[i];
2855  G4VSplitableHadron* projectileSplitable = aNucleon->GetSplitableHadron();
2856  if ( projectileSplitable->GetSoftCollisionCount() != 0 ) {
2857  G4LorentzVector tmp = -DeltaPResidualNucleus;
2858  aNucleon->SetMomentum( tmp );
2859  aNucleon->SetBindingEnergy( DeltaExcitationE );
2860  } else {
2861  delete projectileSplitable;
2862  projectileSplitable = 0;
2863  aNucleon->Hit( projectileSplitable );
2864  aNucleon->SetBindingEnergy( 0.0 );
2865  }
2866  }
2867 
2868  #ifdef debugFTFmodel
2869  G4cout << "NumberOfProjectileParticipant " << NumberOfProjectileParticipant << G4endl
2870  << "ProjectileResidual4Momentum " << ProjectileResidual4Momentum << G4endl;
2871  #endif
2872 
2873  }
2874 
2875  #ifdef debugFTFmodel
2876  G4cout << "End GetResiduals -----------------" << G4endl;
2877  #endif
2878 
2879 }
2880 
2881 
2882 //============================================================================
2883 
2884 G4ThreeVector G4FTFModel::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
2885  // @@ this method is used in FTFModel as well. Should go somewhere common!
2886 
2887  G4double Pt2( 0.0 );
2888  if ( AveragePt2 <= 0.0 ) {
2889  Pt2 = 0.0;
2890  } else {
2891  Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() *
2892  ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
2893  }
2894  G4double Pt = std::sqrt( Pt2 );
2895  G4double phi = G4UniformRand() * twopi;
2896 
2897  return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
2898 }
2899 
2900 
2901 //============================================================================
2902 
2903 G4bool G4FTFModel::
2904 ComputeNucleusProperties( G4V3DNucleus* nucleus, // input parameter
2905  G4LorentzVector& nucleusMomentum, // input & output parameter
2906  G4LorentzVector& residualMomentum, // input & output parameter
2907  G4double& sumMasses, // input & output parameter
2908  G4double& residualExcitationEnergy, // input & output parameter
2909  G4double& residualMass, // input & output parameter
2910  G4int& residualMassNumber, // input & output parameter
2911  G4int& residualCharge ) { // input & output parameter
2912 
2913  // This method, which is called only by PutOnMassShell, computes some nucleus properties for:
2914  // - either the target nucleus (which is never an antinucleus): this for any kind
2915  // of hadronic interaction (hadron-nucleus, nucleus-nucleus, antinucleus-nucleus);
2916  // - or the projectile nucleus or antinucleus: this only in the case of nucleus-nucleus
2917  // or antinucleus-nucleus interaction.
2918  // This method assumes that the all the parameters have been initialized by the caller;
2919  // the action of this method consists in modifying all these parameters, except the
2920  // first one. The return value is "false" only in the case the pointer to the nucleus
2921  // is null.
2922 
2923  if ( ! nucleus ) return false;
2924 
2925  G4double ExcitationEnergyPerWoundedNucleon =
2926  theParameters->GetExcitationEnergyPerWoundedNucleon();
2927 
2928  // Loop over the nucleons of the nucleus.
2929  // The nucleons that have been involved in the interaction (either from Glauber or
2930  // Reggeon Cascading) will be candidate to be emitted.
2931  // All the remaining nucleons will be the nucleons of the candidate residual nucleus.
2932  // The variable sumMasses is the amount of energy corresponding to:
2933  // 1. transverse mass of each involved nucleon
2934  // 2. 20.0*MeV separation energy for each involved nucleon
2935  // 3. transverse mass of the residual nucleus
2936  // In this first evaluation of sumMasses, the excitation energy of the residual nucleus
2937  // (residualExcitationEnergy, estimated by adding a constant value to each involved
2938  // nucleon) is not taken into account.
2939  G4Nucleon* aNucleon = 0;
2940  nucleus->StartLoop();
2941  while ( ( aNucleon = nucleus->GetNextNucleon() ) ) { /* Loop checking, 10.08.2015, A.Ribon */
2942  nucleusMomentum += aNucleon->Get4Momentum();
2943  if ( aNucleon->AreYouHit() ) { // Involved nucleons
2944  // Consider in sumMasses the nominal, i.e. on-shell, masses of the nucleons
2945  // (not the current masses, which could be different because the nucleons are off-shell).
2946  sumMasses += std::sqrt( sqr( aNucleon->GetDefinition()->GetPDGMass() )
2947  + aNucleon->Get4Momentum().perp2() );
2948  sumMasses += 20.0*MeV; // Separation energy for a nucleon
2949 
2950  //AR-11Oct2016 : brought back residual excitation energy as it was in G4 10.1
2951  residualExcitationEnergy += ExcitationEnergyPerWoundedNucleon;
2952  //AR-11Oct2016 residualExcitationEnergy += -ExcitationEnergyPerWoundedNucleon*
2953  //AR-11Oct2016 G4Log( G4UniformRand());
2954 
2955  residualMassNumber--;
2956  // The absolute value below is needed only in the case of anti-nucleus.
2957  residualCharge -= std::abs( G4int( aNucleon->GetDefinition()->GetPDGCharge() ) );
2958  } else { // Spectator nucleons
2959  residualMomentum += aNucleon->Get4Momentum();
2960  }
2961  }
2962  #ifdef debugPutOnMassShell
2963  G4cout << "ExcitationEnergyPerWoundedNucleon " << ExcitationEnergyPerWoundedNucleon << G4endl
2964  << "\t Residual Charge, MassNumber " << residualCharge << " " << residualMassNumber
2965  << G4endl << "\t Initial Momentum " << nucleusMomentum
2966  << G4endl << "\t Residual Momentum " << residualMomentum << G4endl;
2967  #endif
2968  residualMomentum.setPz( 0.0 );
2969  residualMomentum.setE( 0.0 );
2970  if ( residualMassNumber == 0 ) {
2971  residualMass = 0.0;
2972  residualExcitationEnergy = 0.0;
2973  } else {
2974  residualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
2975  GetIonMass( residualCharge, residualMassNumber );
2976  if ( residualMassNumber == 1 ) {
2977  residualExcitationEnergy = 0.0;
2978  }
2979  residualMass += residualExcitationEnergy; // Uzhi March 2016 ?
2980  }
2981  sumMasses += std::sqrt( sqr( residualMass ) + residualMomentum.perp2() );
2982  return true;
2983 }
2984 
2985 
2986 //============================================================================
2987 
2988 G4bool G4FTFModel::
2989 GenerateDeltaIsobar( const G4double sqrtS, // input parameter
2990  const G4int numberOfInvolvedNucleons, // input parameter
2991  G4Nucleon* involvedNucleons[], // input & output parameter
2992  G4double& sumMasses ) { // input & output parameter
2993 
2994  // This method, which is called only by PutOnMassShell, check whether is possible to
2995  // re-interpret some of the involved nucleons as delta-isobars:
2996  // - either by replacing a proton (2212) with a Delta+ (2214),
2997  // - or by replacing a neutron (2112) with a Delta0 (2114).
2998  // The on-shell mass of these delta-isobars is ~1232 MeV, so ~292-294 MeV heavier than
2999  // the corresponding nucleon on-shell mass. However 400.0*MeV is considered to estimate
3000  // the max number of deltas compatible with the available energy.
3001  // The delta-isobars are considered with the same transverse momentum as their
3002  // corresponding nucleons.
3003  // This method assumes that all the parameters have been initialized by the caller;
3004  // the action of this method consists in modifying (eventually) involveNucleons and
3005  // sumMasses. The return value is "false" only in the case that the input parameters
3006  // have unphysical values.
3007 
3008  if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 ) return false;
3009 
3010  //const G4double ProbDeltaIsobar = 0.05;
3011  //const G4double ProbDeltaIsobar = 0.25;
3012  const G4double probDeltaIsobar = 0.05; // A.R. 07.08.2013 0.10 -> 0.05 Uzhi March 2016
3013 
3014  G4int maxNumberOfDeltas = G4int( (sqrtS - sumMasses)/(400.0*MeV) );
3015  G4int numberOfDeltas = 0;
3016 
3017  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3018  //G4cout << "i maxNumberOfDeltas probDeltaIsobar " << i << " " << maxNumberOfDeltas
3019  // << " " << probDeltaIsobar << G4endl;
3020  if ( G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
3021  numberOfDeltas++;
3022  if ( ! involvedNucleons[i] ) continue;
3023  G4VSplitableHadron* splitableHadron = involvedNucleons[i]->GetSplitableHadron();
3024  G4double massNuc = std::sqrt( sqr( splitableHadron->GetDefinition()->GetPDGMass() )
3025  + splitableHadron->Get4Momentum().perp2() );
3026  //AR The absolute value below is needed in the case of an antinucleus.
3027  G4int pdgCode = std::abs( splitableHadron->GetDefinition()->GetPDGEncoding() );
3028  const G4ParticleDefinition* old_def = splitableHadron->GetDefinition();
3029  G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4; // Delta
3030  if ( splitableHadron->GetDefinition()->GetPDGEncoding() < 0 ) newPdgCode *= -1;
3031  const G4ParticleDefinition* ptr =
3032  G4ParticleTable::GetParticleTable()->FindParticle( newPdgCode );
3033  splitableHadron->SetDefinition( ptr );
3034  G4double massDelta = std::sqrt( sqr( splitableHadron->GetDefinition()->GetPDGMass() )
3035  + splitableHadron->Get4Momentum().perp2() );
3036  //G4cout << i << " " << sqrtS/GeV << " " << sumMasses/GeV << " " << massDelta/GeV
3037  // << " " << massNuc << G4endl;
3038  if ( sqrtS < sumMasses + massDelta - massNuc ) { // Change cannot be accepted!
3039  splitableHadron->SetDefinition( old_def );
3040  break;
3041  } else { // Change is accepted
3042  sumMasses += ( massDelta - massNuc ); // Uzhi March 2016 ?
3043  }
3044  }
3045  }
3046  //G4cout << "maxNumberOfDeltas numberOfDeltas " << maxNumberOfDeltas << " "
3047  // << numberOfDeltas << G4endl;
3048  return true;
3049 }
3050 
3051 
3052 //============================================================================
3053 
3054 G4bool G4FTFModel::
3055 SamplingNucleonKinematics( G4double averagePt2, // input parameter
3056  const G4double maxPt2, // input parameter
3057  G4double dCor, // input parameter
3058  G4V3DNucleus* nucleus, // input parameter
3059  const G4LorentzVector& pResidual, // input parameter
3060  const G4double residualMass, // input parameter
3061  const G4int residualMassNumber, // input parameter
3062  const G4int numberOfInvolvedNucleons, // input parameter
3063  G4Nucleon* involvedNucleons[], // input & output parameter
3064  G4double& mass2 ) { // output parameter
3065 
3066  // This method, which is called only by PutOnMassShell, does the sampling of:
3067  // - either the target nucleons: this for any kind of hadronic interactions
3068  // (hadron-nucleus, nucleus-nucleus, antinucleus-nucleus);
3069  // - or the projectile nucleons or antinucleons: this only in the case of
3070  // nucleus-nucleus or antinucleus-nucleus interactions, respectively.
3071  // This method assumes that all the parameters have been initialized by the caller;
3072  // the action of this method consists in changing the properties of the nucleons
3073  // whose pointers are in the vector involvedNucleons, as well as changing the
3074  // variable mass2.
3075 
3076  if ( ! nucleus ) return false;
3077 
3078  if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
3079  dCor = 0.0;
3080  averagePt2 = 0.0;
3081  }
3082 
3083  G4bool success = true;
3084 
3085  G4double SumMasses = residualMass;
3086 
3087  /*
3088  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3089  G4Nucleon* aNucleon = involvedNucleons[i];
3090  if ( ! aNucleon ) continue;
3091  SumMasses += aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
3092  }
3093  */
3094 
3095  const G4int maxNumberOfLoops = 1000;
3096  G4int loopCounter = 0;
3097  do {
3098 
3099  success = true;
3100 
3101  // Sampling of nucleon Pt
3102  G4ThreeVector ptSum( 0.0, 0.0, 0.0 );
3103 
3104  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3105  G4Nucleon* aNucleon = involvedNucleons[i];
3106  if ( ! aNucleon ) continue;
3107  G4ThreeVector tmpPt = GaussianPt( averagePt2, maxPt2 );
3108  ptSum += tmpPt;
3109  G4LorentzVector tmp( tmpPt.x(), tmpPt.y(), 0.0, 0.0 );
3110  aNucleon->SetMomentum( tmp );
3111  }
3112 
3113  G4double deltaPx = ( ptSum.x() - pResidual.x() ) / numberOfInvolvedNucleons;
3114  G4double deltaPy = ( ptSum.y() - pResidual.y() ) / numberOfInvolvedNucleons;
3115 
3116  SumMasses = residualMass;
3117  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3118  G4Nucleon* aNucleon = involvedNucleons[i];
3119  if ( ! aNucleon ) continue;
3120  G4double px = aNucleon->Get4Momentum().px() - deltaPx;
3121  G4double py = aNucleon->Get4Momentum().py() - deltaPy;
3122  G4double MtN = std::sqrt( sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() )
3123  + sqr( px ) + sqr( py ) );
3124  SumMasses += MtN;
3125  G4LorentzVector tmp( px, py, 0.0, MtN);
3126  aNucleon->SetMomentum( tmp );
3127  }
3128 
3129  // Sampling X of nucleon
3130  G4double xSum = 0.0;
3131 
3132  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3133  G4Nucleon* aNucleon = involvedNucleons[i];
3134  if ( ! aNucleon ) continue;
3135 
3136  G4ThreeVector tmpX = GaussianPt( dCor*dCor, 1.0 );
3137  //G4double x = tmpX.x() + aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()/SumMasses;
3138  G4double x = tmpX.x() + aNucleon->Get4Momentum().e()/SumMasses;
3139  if ( x < 0.0 || x > 1.0 ) {
3140  success = false;
3141  break;
3142  }
3143  xSum += x;
3144  //AR The energy is in the lab (instead of cms) frame but it will not be used.
3145  //G4LorentzVector tmp( tmpPt.x(), tmpPt.y(), x, aNucleon->Get4Momentum().e() );
3146  G4LorentzVector tmp( aNucleon->Get4Momentum().x(), aNucleon->Get4Momentum().y(),
3147  x, aNucleon->Get4Momentum().e() );
3148  aNucleon->SetMomentum( tmp );
3149  }
3150 
3151  if ( xSum < 0.0 || xSum > 1.0 ) success = false;
3152 
3153  if ( ! success ) continue;
3154 
3155  //G4double deltaPx = ( ptSum.x() - pResidual.x() ) / numberOfInvolvedNucleons;
3156  //G4double deltaPy = ( ptSum.y() - pResidual.y() ) / numberOfInvolvedNucleons;
3157  G4double delta = 0.0;
3158  if ( residualMassNumber == 0 ) {
3159  delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
3160  } else {
3161  delta = 0.0;
3162  }
3163 
3164  xSum = 1.0;
3165  mass2 = 0.0;
3166  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3167  G4Nucleon* aNucleon = involvedNucleons[i];
3168  if ( ! aNucleon ) continue;
3169  G4double x = aNucleon->Get4Momentum().pz() - delta;
3170  xSum -= x;
3171  if ( residualMassNumber == 0 ) {
3172  if ( x <= 0.0 || x > 1.0 ) {
3173  success = false;
3174  break;
3175  }
3176  } else {
3177  if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
3178  success = false;
3179  break;
3180  }
3181  }
3182 
3183  /*
3184  G4double px = aNucleon->Get4Momentum().px() - deltaPx;
3185  G4double py = aNucleon->Get4Momentum().py() - deltaPy;
3186  mass2 += ( sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() )
3187  + sqr( px ) + sqr( py ) ) / x;
3188  G4LorentzVector tmp( px, py, x, aNucleon->Get4Momentum().e() );
3189  */
3190 
3191  mass2 += sqr( aNucleon->Get4Momentum().e() ) / x;
3192  G4LorentzVector tmp( aNucleon->Get4Momentum().px(), aNucleon->Get4Momentum().py(),
3193  x, aNucleon->Get4Momentum().e() );
3194  aNucleon->SetMomentum( tmp );
3195  }
3196  if ( ! success ) continue;
3197 
3198  if ( success && residualMassNumber != 0 ) {
3199  //mass2 += ( sqr( residualMass ) + pResidual.perp2() ) / xSum;
3200  mass2 += sqr( residualMass ) / xSum;
3201  }
3202 
3203  #ifdef debugPutOnMassShell
3204  G4cout << "success " << success << G4endl << " Mt " << std::sqrt( mass2 )/GeV << G4endl;
3205  #endif
3206 
3207  } while ( ( ! success ) &&
3208  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
3209  if ( loopCounter >= maxNumberOfLoops ) {
3210  return false;
3211  }
3212 
3213  return true;
3214 }
3215 
3216 
3217 //============================================================================
3218 
3219 G4bool G4FTFModel::
3220 CheckKinematics( const G4double sValue, // input parameter
3221  const G4double sqrtS, // input parameter
3222  const G4double projectileMass2, // input parameter
3223  const G4double targetMass2, // input parameter
3224  const G4double nucleusY, // input parameter
3225  const G4bool isProjectileNucleus, // input parameter
3226  const G4int numberOfInvolvedNucleons, // input parameter
3227  G4Nucleon* involvedNucleons[], // input parameter
3228  G4double& targetWminus, // output parameter
3229  G4double& projectileWplus, // output parameter
3230  G4bool& success ) { // input & output parameter
3231 
3232  // This method, which is called only by PutOnMassShell, checks whether the
3233  // kinematics is acceptable or not.
3234  // This method assumes that all the parameters have been initialized by the caller;
3235  // notice that the input boolean parameter isProjectileNucleus is meant to be true
3236  // only in the case of nucleus or antinucleus projectile.
3237  // The action of this method consists in computing targetWminus and projectileWplus
3238  // and setting the parameter success to false in the case that the kinematics should
3239  // be rejeted.
3240 
3241  G4double decayMomentum2 = sqr( sValue ) + sqr( projectileMass2 ) + sqr( targetMass2 )
3242  - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
3243  - 2.0*projectileMass2*targetMass2;
3244  targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
3245  / 2.0 / sqrtS;
3246  projectileWplus = sqrtS - targetMass2/targetWminus;
3247  G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
3248  G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
3249  G4double projectileY = 0.5 * G4Log( (projectileE + projectilePz)/
3250  (projectileE - projectilePz) );
3251  G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
3252  G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
3253  G4double targetY = 0.5 * G4Log( (targetE + targetPz)/(targetE - targetPz) );
3254 
3255  #ifdef debugPutOnMassShell
3256  G4cout << "decayMomentum2 " << decayMomentum2 << G4endl
3257  << "\t targetWminus projectileWplus " << targetWminus << " " << projectileWplus << G4endl
3258  << "\t projectileY targetY " << projectileY << " " << targetY << G4endl;
3259  #endif
3260 
3261  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3262  G4Nucleon* aNucleon = involvedNucleons[i];
3263  if ( ! aNucleon ) continue;
3264  G4LorentzVector tmp = aNucleon->Get4Momentum();
3265  G4double mt2 = sqr( tmp.x() ) + sqr( tmp.y() ) +
3266  sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() );
3267  G4double x = tmp.z();
3268  G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
3269  G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
3270  if ( isProjectileNucleus ) {
3271  pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*x);
3272  e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*x);
3273  }
3274  G4double nucleonY = 0.5 * G4Log( (e + pz)/(e - pz) );
3275 
3276  #ifdef debugPutOnMassShell
3277  G4cout << "i nY pY nY-AY AY " << i << " " << nucleonY << " " << projectileY <<G4endl;
3278  #endif
3279 
3280  if ( std::abs( nucleonY - nucleusY ) > 2 ||
3281  ( isProjectileNucleus && targetY > nucleonY ) ||
3282  ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
3283  success = false;
3284  break;
3285  }
3286  }
3287  return true;
3288 }
3289 
3290 
3291 //============================================================================
3292 
3293 G4bool G4FTFModel::
3294 FinalizeKinematics( const G4double w, // input parameter
3295  const G4bool isProjectileNucleus, // input parameter
3296  const G4LorentzRotation& boostFromCmsToLab, // input parameter
3297  const G4double residualMass, // input parameter
3298  const G4int residualMassNumber, // input parameter
3299  const G4int numberOfInvolvedNucleons, // input parameter
3300  G4Nucleon* involvedNucleons[], // input & output parameter
3301  G4LorentzVector& residual4Momentum ) { // output parameter
3302 
3303  // This method, which is called only by PutOnMassShell, finalizes the kinematics:
3304  // this method is called when we are sure that the sampling of the kinematics is
3305  // acceptable.
3306  // This method assumes that all the parameters have been initialized by the caller;
3307  // notice that the input boolean parameter isProjectileNucleus is meant to be true
3308  // only in the case of nucleus or antinucleus projectile: this information is needed
3309  // because the sign of pz (in the center-of-mass frame) in this case is opposite
3310  // with respect to the case of a normal hadron projectile.
3311  // The action of this method consists in modifying the momenta of the nucleons
3312  // (in the lab frame) and computing the residual 4-momentum (in the center-of-mass
3313  // frame).
3314 
3315  G4ThreeVector residual3Momentum( 0.0, 0.0, 1.0 );
3316 
3317  for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
3318  G4Nucleon* aNucleon = involvedNucleons[i];
3319  if ( ! aNucleon ) continue;
3320  G4LorentzVector tmp = aNucleon->Get4Momentum();
3321  residual3Momentum -= tmp.vect();
3322  G4double mt2 = sqr( tmp.x() ) + sqr( tmp.y() ) +
3323  sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() );
3324  G4double x = tmp.z();
3325  G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w * x );
3326  G4double e = w * x / 2.0 + mt2 / ( 2.0 * w * x );
3327  // Reverse the sign of pz in the case of nucleus or antinucleus projectile
3328  if ( isProjectileNucleus ) pz *= -1.0;
3329  tmp.setPz( pz );
3330  tmp.setE( e );
3331  tmp.transform( boostFromCmsToLab );
3332  aNucleon->SetMomentum( tmp );
3333  G4VSplitableHadron* splitableHadron = aNucleon->GetSplitableHadron();
3334  splitableHadron->Set4Momentum( tmp );
3335  }
3336 
3337  G4double residualMt2 = sqr( residualMass ) + sqr( residual3Momentum.x() )
3338  + sqr( residual3Momentum.y() );
3339 
3340  #ifdef debugPutOnMassShell
3341  G4cout << "w residual3Momentum.z() " << w << " " << residual3Momentum.z() << G4endl;
3342  #endif
3343 
3344  G4double residualPz = 0.0;
3345  G4double residualE = 0.0;
3346  if ( residualMassNumber != 0 ) {
3347  residualPz = -w * residual3Momentum.z() / 2.0 +
3348  residualMt2 / ( 2.0 * w * residual3Momentum.z() );
3349  residualE = w * residual3Momentum.z() / 2.0 +
3350  residualMt2 / ( 2.0 * w * residual3Momentum.z() );
3351  // Reverse the sign of residualPz in the case of nucleus or antinucleus projectile
3352  if ( isProjectileNucleus ) residualPz *= -1.0;
3353  }
3354 
3355  residual4Momentum.setPx( residual3Momentum.x() );
3356  residual4Momentum.setPy( residual3Momentum.y() );
3357  residual4Momentum.setPz( residualPz );
3358  residual4Momentum.setE( residualE );
3359 
3360  return true;
3361 }
3362 
3363 
3364 //============================================================================
3365 
3366 void G4FTFModel::ModelDescription( std::ostream& desc ) const {
3367  desc << " FTF (Fritiof) Model \n"
3368  << "The FTF model is based on the well-known FRITIOF \n"
3369  << "model (B. Andersson et al., Nucl. Phys. B281, 289 \n"
3370  << "(1987)). Its first program implementation was given\n"
3371  << "by B. Nilsson-Almquist and E. Stenlund (Comp. Phys.\n"
3372  << "Comm. 43, 387 (1987)). The Fritiof model assumes \n"
3373  << "that all hadron-hadron interactions are binary \n"
3374  << "reactions, h_1+h_2->h_1'+h_2' where h_1' and h_2' \n"
3375  << "are excited states of the hadrons with continuous \n"
3376  << "mass spectra. The excited hadrons are considered as\n"
3377  << "QCD-strings, and the corresponding LUND-string \n"
3378  << "fragmentation model is applied for a simulation of \n"
3379  << "their decays. \n"
3380  << " The Fritiof model assumes that in the course of \n"
3381  << "a hadron-nucleus interaction a string originated \n"
3382  << "from the projectile can interact with various intra\n"
3383  << "nuclear nucleons and becomes into highly excited \n"
3384  << "states. The probability of multiple interactions is\n"
3385  << "calculated in the Glauber approximation. A cascading\n"
3386  << "of secondary particles was neglected as a rule. Due\n"
3387  << "to these, the original Fritiof model fails to des- \n"
3388  << "cribe a nuclear destruction and slow particle spectra.\n"
3389  << " In order to overcome the difficulties we enlarge\n"
3390  << "the model by the reggeon theory inspired model of \n"
3391  << "nuclear desctruction (Kh. Abdel-Waged and V.V. Uzhi-\n"
3392  << "nsky, Phys. Atom. Nucl. 60, 828 (1997); Yad. Fiz. 60, 925\n"
3393  << "(1997)). Momenta of the nucleons ejected from a nuc-\n"
3394  << "leus in the reggeon cascading are sampled according\n"
3395  << "to a Fermi motion algorithm presented in (EMU-01 \n"
3396  << "Collaboration (M.I. Adamovich et al.) Zeit. fur Phys.\n"
3397  << "A358, 337 (1997)). \n"
3398  << " New features were also added to the Fritiof model\n"
3399  << "implemented in Geant4: a simulation of elastic had-\n"
3400  << "ron-nucleon scatterings, a simulation of binary \n"
3401  << "reactions like NN>NN* in hadron-nucleon interactions,\n"
3402  << "a separate simulation of single diffractive and non-\n"
3403  << " diffractive events. These allowed to describe after\n"
3404  << "model parameter tuning a wide set of experimental \n"
3405  << "data. \n";
3406 }
3407 
G4Exp.hh
G4ElasticHNScattering::ElasticScattering
virtual G4bool ElasticScattering(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters) const
Definition: G4ElasticHNScattering.cc:64
G4FTFParticipants::GetList
void GetList(const G4ReactionProduct &thePrimary, G4FTFParameters *theParameters)
Definition: G4FTFParticipants.cc:83
G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:177
G4FTFParameters.hh
G4ExcitedStringVector
std::vector< G4ExcitedString * > G4ExcitedStringVector
Definition: G4ExcitedStringVector.hh:41
G4InteractionContent.hh
G4InteractionContent::GetTargetNucleon
G4Nucleon * GetTargetNucleon() const
Definition: G4InteractionContent.hh:149
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:566
G4InteractionContent::GetProjectileNucleon
G4Nucleon * GetProjectileNucleon() const
Definition: G4InteractionContent.hh:139
G4Nucleon::SetParticleType
void SetParticleType(G4Proton *aProton)
Definition: G4Nucleon.hh:77
target
const XML_Char * target
Definition: expat.h:268
G4FTFParameters::GetProbabilityOfElasticScatt
G4double GetProbabilityOfElasticScatt()
Definition: G4FTFParameters.hh:401
G4ReactionProduct::GetTotalMomentum
G4double GetTotalMomentum() const
Definition: G4ReactionProduct.hh:126
G4FTFModel::G4FTFModel
G4FTFModel(const G4String &modelName="FTF")
Definition: G4FTFModel.cc:72
S
double S(double temp)
Definition: G4DNAElectronHoleRecombination.cc:77
G4Nucleon::SetMomentum
void SetMomentum(G4LorentzVector &aMomentum)
Definition: G4Nucleon.hh:70
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
DeleteVSplitableHadron::operator()
void operator()(G4VSplitableHadron *aH)
Definition: G4FTFModel.cc:109
CLHEP::Hep3Vector::x
double x() const
G4DiffractiveExcitation::ExciteParticipants
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
Definition: G4DiffractiveExcitation.cc:94
G4V3DNucleus::StartLoop
virtual G4bool StartLoop()=0
perCent
static constexpr double perCent
Definition: G4SIunits.hh:332
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:171
G4InteractionContent
Definition: G4InteractionContent.hh:55
G4V3DNucleus::DoLorentzContraction
virtual void DoLorentzContraction(const G4LorentzVector &theBoost)=0
G4V3DNucleus::GetMassNumber
virtual G4int GetMassNumber()=0
G4VSplitableHadron::SetTimeOfCreation
void SetTimeOfCreation(G4double aTime)
Definition: G4VSplitableHadron.hh:148
G4Nucleon::Get4Momentum
virtual const G4LorentzVector & Get4Momentum() const
Definition: G4Nucleon.hh:72
test.x
tuple x
Definition: test.py:50
C
double C(double temp)
Definition: G4DNAElectronHoleRecombination.cc:86
G4FTFModel::GetTargetNucleus
G4V3DNucleus * GetTargetNucleus() const
Definition: G4FTFModel.hh:165
G4Nucleon::GetPosition
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
G4VSplitableHadron::SetDefinition
void SetDefinition(const G4ParticleDefinition *aDefinition)
Definition: G4VSplitableHadron.hh:133
G4VParticipants::SetProjectileNucleus
virtual void SetProjectileNucleus(G4V3DNucleus *aNucleus)
Definition: G4VParticipants.cc:79
CLHEP::HepLorentzVector::rotateZ
HepLorentzVector & rotateZ(double)
Definition: LorentzVector.cc:80
G4int
int G4int
Definition: G4Types.hh:78
G4VSplitableHadron::SetStatus
void SetStatus(const G4int aStatus)
Definition: G4VSplitableHadron.hh:168
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::Hep3Vector::z
double z() const
G4IonTable.hh
CLHEP::HepLorentzVector::phi
double phi() const
twopi
static constexpr double twopi
Definition: G4SIunits.hh:76
G4Nucleon::GetDefinition
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
G4VPartonStringModel::SetThisPointer
void SetThisPointer(G4VPartonStringModel *aPointer)
Definition: G4VPartonStringModel.hh:87
G4Nucleon::GetSplitableHadron
G4VSplitableHadron * GetSplitableHadron() const
Definition: G4Nucleon.hh:96
G4InteractionContent::GetStatus
G4int GetStatus() const
Definition: G4InteractionContent.cc:70
G4FTFParameters::GetCofNuclearDestructionPr
G4double GetCofNuclearDestructionPr()
Definition: G4FTFParameters.hh:477
G4VSplitableHadron::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4VSplitableHadron.hh:138
CLHEP::HepLorentzVector::rapidity
double rapidity() const
Definition: LorentzVectorK.cc:128
G4FTFParticipants::GetInteraction
G4InteractionContent & GetInteraction()
Definition: G4FTFParticipants.hh:84
G4ReactionProduct::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:107
G4DiffractiveSplitableHadron.hh
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:658
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4FTFParameters::GetMaxNumberOfCollisions
G4double GetMaxNumberOfCollisions()
Definition: G4FTFParameters.hh:469
G4FTFModel.hh
G4ParticleTable.hh
CLHEP::HepLorentzVector::theta
double theta() const
G4VParticipants::Init
virtual void Init(G4int theZ, G4int theA)
Definition: G4VParticipants.cc:57
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4VHighEnergyGenerator::SetEnergyMomentumCheckLevels
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double AbsoluteLevel)
Definition: G4VHighEnergyGenerator.cc:51
CLHEP::HepLorentzVector::mag
double mag() const
CLHEP::HepLorentzVector::py
double py() const
G4FTFModel::GetProjectileNucleus
G4V3DNucleus * GetProjectileNucleus() const
Definition: G4FTFModel.hh:170
G4bool
bool G4bool
Definition: G4Types.hh:79
G4IonTable::GetIonMass
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1335
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
G4Nucleon::AreYouHit
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
CLHEP::HepLorentzVector::px
double px() const
G4VParticipants::theProjectileNucleus
G4V3DNucleus * theProjectileNucleus
Definition: G4VParticipants.hh:68
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
CLHEP::HepLorentzRotation::transform
HepLorentzRotation & transform(const HepBoost &b)
G4FTFParameters::GetMaxPt2ofNuclearDestruction
G4double GetMaxPt2ofNuclearDestruction()
Definition: G4FTFParameters.hh:501
G4ParticleDefinition.hh
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Nucleon::SetBindingEnergy
void SetBindingEnergy(G4double anEnergy)
Definition: G4Nucleon.hh:74
G4Log.hh
G4KineticTrack.hh
G4FTFModel::Init
void Init(const G4Nucleus &aNucleus, const G4DynamicParticle &aProjectile)
Definition: G4FTFModel.cc:150
G4PhysicalConstants.hh
G4VSplitableHadron::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4VSplitableHadron.hh:128
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4FTFModel::ModelDescription
virtual void ModelDescription(std::ostream &) const
Definition: G4FTFModel.cc:3366
G4DiffractiveExcitation::CreateStrings
virtual void CreateStrings(G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const
Definition: G4DiffractiveExcitation.cc:1193
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4VParticipants::InitProjectileNucleus
virtual void InitProjectileNucleus(G4int theZ, G4int theA)
Definition: G4VParticipants.cc:71
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4V3DNucleus::DoLorentzBoost
virtual void DoLorentzBoost(const G4LorentzVector &theBoost)=0
G4InteractionContent::SetStatus
void SetStatus(G4int aValue)
Definition: G4InteractionContent.cc:67
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4FTFParameters::GetPt2ofNuclearDestruction
G4double GetPt2ofNuclearDestruction()
Definition: G4FTFParameters.hh:497
G4FTFParticipants.hh
G4FTFParameters::GetR2ofNuclearDestruction
G4double GetR2ofNuclearDestruction()
Definition: G4FTFParameters.hh:485
CLHEP::HepLorentzVector::mag2
double mag2() const
G4ios.hh
G4FTFParameters::GetExcitationEnergyPerWoundedNucleon
G4double GetExcitationEnergyPerWoundedNucleon()
Definition: G4FTFParameters.hh:489
CLHEP::Hep3Vector::y
double y() const
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
CLHEP::HepLorentzVector::rotateY
HepLorentzVector & rotateY(double)
Definition: LorentzVector.cc:76
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
CLHEP::Hep3Vector::mag2
double mag2() const
G4InteractionContent::GetTarget
G4VSplitableHadron * GetTarget() const
Definition: G4InteractionContent.hh:124
G4InteractionContent::GetProjectile
G4VSplitableHadron * GetProjectile() const
Definition: G4InteractionContent.hh:119
CLHEP::HepLorentzVector::perp2
double perp2() const
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
CLHEP::HepLorentzVector::pz
double pz() const
G4VSplitableHadron::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4VSplitableHadron.hh:163
G4FTFAnnihilation::Annihilate
virtual G4bool Annihilate(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4VSplitableHadron *&AdditionalString, G4FTFParameters *theParameters) const
Definition: G4FTFAnnihilation.cc:88
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4FTFParameters::GetCofNuclearDestruction
G4double GetCofNuclearDestruction()
Definition: G4FTFParameters.hh:481
CLHEP::HepLorentzVector::transform
HepLorentzVector & transform(const HepRotation &)
G4FTFParameters::GetProbabilityOfAnnihilation
G4double GetProbabilityOfAnnihilation()
Definition: G4FTFParameters.hh:410
G4V3DNucleus::GetNextNucleon
virtual G4Nucleon * GetNextNucleon()=0
G4VSplitableHadron::Set4Momentum
void Set4Momentum(const G4LorentzVector &a4Momentum)
Definition: G4VSplitableHadron.hh:123
G4Nucleon::Hit
void Hit(G4VSplitableHadron *aHit)
Definition: G4Nucleon.hh:90
sqr
T sqr(const T &x)
Definition: templates.hh:145
CLHEP::HepLorentzVector::setVect
void setVect(const Hep3Vector &)
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4Nucleon::GetBindingEnergy
G4double GetBindingEnergy() const
Definition: G4Nucleon.hh:75
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4FTFModel::GetStrings
G4ExcitedStringVector * GetStrings()
Definition: G4FTFModel.cc:272
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4LorentzVector
CLHEP::HepLorentzVector G4LorentzVector
Definition: G4LorentzVector.hh:35
G4FTFParameters::GetDofNuclearDestruction
G4double GetDofNuclearDestruction()
Definition: G4FTFParameters.hh:493
G4LorentzRotation.hh