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