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