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