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