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