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