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