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