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