| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
35 //
36 // Vladimir Uzhinsky, November
37 // simulation of nucleus-nucleus interac
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& modelN
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, 1
112 }
113
114
115 //============================================
116
117 struct DeleteVSplitableHadron { void operator(
118
119
120 //============================================
121
122 G4FTFModel::~G4FTFModel() {
123 // Because FTF model can be called for vari
124 //
125 // ---> NOTE (JVY): This statement below is
126 // theParameters must be erased at the end
127 // Thus the delete is also in G4FTFModel::G
128 // ---> JVY
129 //
130 if ( theParameters != 0 ) delete theParam
131 if ( theExcitation != 0 ) delete theExcit
132 if ( theElastic != 0 ) delete theElast
133 if ( theAnnihilation != 0 ) delete theAnnih
134
135 // Erasing of strings created at annihilati
136 if ( theAdditionalString.size() != 0 ) {
137 std::for_each( theAdditionalString.begin(
138 DeleteVSplitableHadron() )
139 }
140 theAdditionalString.clear();
141
142 // Erasing of target involved nucleons.
143 if ( NumberOfInvolvedNucleonsOfTarget != 0
144 for ( G4int i = 0; i < NumberOfInvolvedNu
145 G4VSplitableHadron* aNucleon = TheInvol
146 if ( aNucleon ) delete aNucleon;
147 }
148 }
149
150 // Erasing of projectile involved nucleons.
151 if ( NumberOfInvolvedNucleonsOfProjectile !
152 for ( G4int i = 0; i < NumberOfInvolvedNu
153 G4VSplitableHadron* aNucleon = TheInvol
154 if ( aNucleon ) delete aNucleon;
155 }
156 }
157 }
158
159
160 //============================================
161
162 void G4FTFModel::Init( const G4Nucleus& aNucle
163
164 theProjectile = aProjectile;
165
166 G4double PlabPerParticle( 0.0 ); // Laborat
167
168 #ifdef debugFTFmodel
169 G4cout << "FTF init Proj Name " << theProjec
170 << "FTF init Proj Mass " << theProjec
171 << " " << theProjectile.GetMomentum()
172 << "FTF init Proj B Q " << theProjec
173 << " " << (G4int) theProjectile.GetDe
174 << "FTF init Target A Z " << aNucleus
175 << " " << aNucleus.GetZ_asInt() << G4
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.Ge
190 TargetResidualCharge = aNucleus.Ge
191 TargetResidualExcitationEnergy = 0.0;
192 TargetResidual4Momentum = tmp;
193 G4double TargetResidualMass = G4ParticleTabl
194 ->GetIonMass(
195
196 TargetResidual4Momentum.setE( TargetResidual
197
198 if ( std::abs( theProjectile.GetDefinition()
199 // Projectile is a hadron : meson or baryo
200 ProjectileResidualMassNumber = std::abs( t
201 ProjectileResidualCharge = G4int( theProje
202 PlabPerParticle = theProjectile.GetMomentu
203 ProjectileResidualExcitationEnergy = 0.0;
204 //G4double ProjectileResidualMass = thePro
205 ProjectileResidual4Momentum.setVect( thePr
206 ProjectileResidual4Momentum.setE( theProje
207 if ( PlabPerParticle < LowEnergyLimit ) {
208 HighEnergyInter = false;
209 } else {
210 HighEnergyInter = true;
211 }
212 } else {
213 if ( theProjectile.GetDefinition()->GetBar
214 // Projectile is a nucleus
215 ProjectileResidualMassNumber = theProjec
216 ProjectileResidualCharge = G4int( thePro
217 ProjectileResidualLambdaNumber = theProj
218 PlabPerParticle = theProjectile.GetMomen
219 if ( PlabPerParticle < LowEnergyLimit )
220 HighEnergyInter = false;
221 } else {
222 HighEnergyInter = true;
223 }
224 theParticipants.InitProjectileNucleus( P
225 ProjectileResidualLambdaNumber
226 } else if ( theProjectile.GetDefinition()-
227 // Projectile is an anti-nucleus
228 ProjectileResidualMassNumber = std::abs(
229 ProjectileResidualCharge = std::abs( G4i
230 ProjectileResidualLambdaNumber = theProj
231 PlabPerParticle = theProjectile.GetMomen
232 if ( PlabPerParticle < LowEnergyLimit )
233 HighEnergyInter = false;
234 } else {
235 HighEnergyInter = true;
236 }
237 theParticipants.InitProjectileNucleus( P
238 P
239 theParticipants.GetProjectileNucleus()->
240 G4Nucleon* aNucleon;
241 while ( ( aNucleon = theParticipants.Get
242 if ( aNucleon->GetDefinition() == G4Pr
243 aNucleon->SetParticleType( G4AntiPro
244 } else if ( aNucleon->GetDefinition()
245 aNucleon->SetParticleType( G4AntiNeu
246 } else if ( aNucleon->GetDefinition()
247 aNucleon->SetParticleType( G4AntiLambda::D
248 }
249 }
250 }
251
252 G4ThreeVector BoostVector = theProjectile.
253 theParticipants.GetProjectileNucleus()->Do
254 theParticipants.GetProjectileNucleus()->Do
255 ProjectileResidualExcitationEnergy = 0.0;
256 //G4double ProjectileResidualMass = thePro
257 ProjectileResidual4Momentum.setVect( thePr
258 ProjectileResidual4Momentum.setE( theProje
259 }
260
261 // Init target nucleus (assumed to be never
262 theParticipants.Init( aNucleus.GetA_asInt(),
263
264 NumberOfProjectileSpectatorNucleons = std::a
265 NumberOfTargetSpectatorNucleons = aNucleus.G
266 NumberOfNNcollisions = 0;
267
268 // reset/recalculate everything for the new
269 theParameters->InitForInteraction( theProjec
270 aNucleus.
271
272 if ( theAdditionalString.size() != 0 ) {
273 std::for_each( theAdditionalString.begin()
274 DeleteVSplitableHadron() );
275 }
276 theAdditionalString.clear();
277
278 #ifdef debugFTFmodel
279 G4cout << "FTF end of Init" << G4endl << G4e
280 #endif
281
282 // In the case of Hydrogen target, for non-i
283 // do NOT simulate quasi-elastic (by forcing
284 // elastic scatering in theParameters - whic
285 // This is necessary because in this case qu
286 // with only one nucleon would be identical
287 // and the latter is already included in the
288 // (i.e. G4HadronElasticProcess).
289 if ( std::abs( theProjectile.GetDefinition()
290 aNucleus.GetA_asInt() < 2 ) theParamete
291
292 if ( SampleBinInterval() ) theParticipants.S
293 }
294
295
296 //============================================
297
298 G4ExcitedStringVector* G4FTFModel::GetStrings(
299
300 #ifdef debugFTFmodel
301 G4cout << "G4FTFModel::GetStrings() " << G4e
302 #endif
303
304 G4ExcitedStringVector* theStrings = new G4Ex
305 theParticipants.GetList( theProjectile, theP
306
307 SetImpactParameter( theParticipants.GetImpac
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? "
324 #endif
325
326 }
327
328 #ifdef debugFTFmodel
329 G4cout << "FTF ExciteParticipants " << G4end
330 #endif
331
332 if ( Success ) Success = ExciteParticipants(
333
334 #ifdef debugFTFmodel
335 G4cout << "FTF ExciteParticipants Success? "
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 " << theString
348 << "FTF GetResiduals of Nuclei " <<
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* > primari
362 theParticipants.StartLoop();
363 while ( theParticipants.Next() ) { /* Loo
364 const G4InteractionContent& interaction
365 // Do not allow for duplicates
366 if ( primaries.end() ==
367 std::find( primaries.begin(), prima
368 primaries.push_back( interaction.GetPr
369 }
370 }
371 std::for_each( primaries.begin(), primarie
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 < NumberOfInvolvedNucle
380 aNucleon = TheInvolvedNucleonsOfProjectile
381 if ( aNucleon ) delete aNucleon;
382 }
383 NumberOfInvolvedNucleonsOfProjectile = 0;
384
385 // Erase the target nucleons
386 for ( G4int i = 0; i < NumberOfInvolvedNucle
387 aNucleon = TheInvolvedNucleonsOfTarget[i]-
388 if ( aNucleon ) delete aNucleon;
389 }
390 NumberOfInvolvedNucleonsOfTarget = 0;
391
392 #ifdef debugFTFmodel
393 G4cout << "End of FTF. Go to fragmentation"
394 << "To continue - enter 1, to stop -
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 interact
407
408 NumberOfInvolvedNucleonsOfTarget = 0;
409
410 G4V3DNucleus* theTargetNucleus = GetTargetNu
411 theTargetNucleus->StartLoop();
412
413 G4Nucleon* aNucleon;
414 while ( ( aNucleon = theTargetNucleus->GetNe
415 if ( aNucleon->AreYouHit() ) {
416 TheInvolvedNucleonsOfTarget[NumberOfInvo
417 NumberOfInvolvedNucleonsOfTarget++;
418 }
419 }
420
421 #ifdef debugFTFmodel
422 G4cout << "G4FTFModel::StoreInvolvedNucleon
423 G4cout << "NumberOfInvolvedNucleonsOfTarget
424 << G4endl << G4endl;
425 #endif
426
427
428 if ( ! GetProjectileNucleus() ) return; //
429
430 // The projectile is a nucleus or an anti-nu
431
432 NumberOfInvolvedNucleonsOfProjectile = 0;
433
434 G4V3DNucleus* theProjectileNucleus = GetProj
435 theProjectileNucleus->StartLoop();
436
437 G4Nucleon* aProjectileNucleon;
438 while ( ( aProjectileNucleon = theProjectile
439 if ( aProjectileNucleon->AreYouHit() ) {
440 // Projectile nucleon was involved in th
441 TheInvolvedNucleonsOfProjectile[NumberOf
442 NumberOfInvolvedNucleonsOfProjectile++;
443 }
444 }
445
446 #ifdef debugFTFmodel
447 G4cout << "NumberOfInvolvedNucleonsOfProject
448 << G4endl << G4endl;
449 #endif
450 return;
451 }
452
453
454 //============================================
455
456 void G4FTFModel::ReggeonCascade() {
457 // Implementation of the reggeon theory insp
458
459 #ifdef debugReggeonCascade
460 G4cout << "G4FTFModel::ReggeonCascade ------
461 << "theProjectile.GetTotalMomentum()
462 << "theProjectile.GetTotalEnergy() "
463 << "ExcitationE/WN " << theParameters
464 #endif
465
466 G4int InitNINt = NumberOfInvolvedNucleonsOfT
467
468 // Reggeon cascading in target nucleus
469 for ( G4int InvTN = 0; InvTN < InitNINt; Inv
470 G4Nucleon* aTargetNucleon = TheInvolvedNuc
471
472 G4double CreationTime = aTargetNucleon->Ge
473
474 G4double XofWoundedNucleon = aTargetNucleo
475 G4double YofWoundedNucleon = aTargetNucleo
476
477 G4V3DNucleus* theTargetNucleus = GetTarget
478 theTargetNucleus->StartLoop();
479
480 G4Nucleon* Neighbour(0);
481 while ( ( Neighbour = theTargetNucleus->Ge
482 if ( ! Neighbour->AreYouHit() ) {
483 G4double impact2 = sqr( XofWoundedNucl
484 sqr( YofWoundedNucl
485
486 if ( G4UniformRand() < theParameters->
487 G4Exp( -impact2
488 ) {
489 // The neighbour nucleon is involved
490 TheInvolvedNucleonsOfTarget[ NumberO
491 NumberOfInvolvedNucleonsOfTarget++;
492
493 G4VSplitableHadron* targetSplitable;
494 targetSplitable = new G4DiffractiveS
495
496 Neighbour->Hit( targetSplitable );
497 targetSplitable->SetTimeOfCreation(
498 targetSplitable->SetStatus( 3 ); //
499 }
500 }
501 }
502 }
503
504 #ifdef debugReggeonCascade
505 G4cout << "Final NumberOfInvolvedNucleonsOfT
506 << NumberOfInvolvedNucleonsOfTarget <
507 #endif
508
509 if ( ! GetProjectileNucleus() ) return;
510
511 // Nucleus-Nucleus Interaction : Destruction
512 G4int InitNINp = NumberOfInvolvedNucleonsOfP
513
514 //for ( G4int InvPN = 0; InvPN < NumberOfInv
515 for ( G4int InvPN = 0; InvPN < InitNINp; Inv
516 G4Nucleon* aProjectileNucleon = TheInvolve
517
518 G4double CreationTime = aProjectileNucleon
519
520 G4double XofWoundedNucleon = aProjectileNu
521 G4double YofWoundedNucleon = aProjectileNu
522
523 G4V3DNucleus* theProjectileNucleus = GetPr
524 theProjectileNucleus->StartLoop();
525
526 G4Nucleon* Neighbour( 0 );
527 while ( ( Neighbour = theProjectileNucleus
528 if ( ! Neighbour->AreYouHit() ) {
529 G4double impact2= sqr( XofWoundedNucle
530 sqr( YofWoundedNucle
531
532 if ( G4UniformRand() < theParameters->
533 G4Exp( -impact2
534 ) {
535 // The neighbour nucleon is involved
536 TheInvolvedNucleonsOfProjectile[ Num
537 NumberOfInvolvedNucleonsOfProjectile
538
539 G4VSplitableHadron* projectileSplita
540 projectileSplitable = new G4Diffract
541
542 Neighbour->Hit( projectileSplitable
543 projectileSplitable->SetTimeOfCreati
544 projectileSplitable->SetStatus( 3 );
545 }
546 }
547 }
548 }
549
550 #ifdef debugReggeonCascade
551 G4cout << "NumberOfInvolvedNucleonsOfProject
552 << NumberOfInvolvedNucleonsOfProjecti
553 #endif
554 }
555
556
557 //============================================
558
559 G4bool G4FTFModel::PutOnMassShell() {
560
561 G4bool isProjectileNucleus = false;
562 if ( GetProjectileNucleus() ) isProjectileNu
563
564 #ifdef debugPutOnMassShell
565 G4cout << "PutOnMassShell start " << G4endl;
566 if ( isProjectileNucleus ) {
567 G4cout << "PutOnMassShell for Nucleus_Nucl
568 }
569 #endif
570
571 G4LorentzVector Pprojectile( theProjectile.G
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
578 G4double SumMasses = 0.0;
579 G4V3DNucleus* theTargetNucleus = GetTargetNu
580 G4double TargetResidualMass = 0.0;
581
582 #ifdef debugPutOnMassShell
583 G4cout << "Target : ";
584 #endif
585 isOk = ComputeNucleusProperties( theTargetNu
586 TargetResid
587 TargetResid
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
594 G4V3DNucleus* thePrNucleus = GetProjectileNu
595 G4double PrResidualMass = 0.0;
596
597 if ( ! isProjectileNucleus ) { // hadron-nu
598 Mprojectile = Pprojectile.mag();
599 M2projectile = Pprojectile.mag2();
600 SumMasses += Mprojectile + 20.0*MeV;
601 } else { // nucleus-nucleus or antinucleus-
602 #ifdef debugPutOnMassShell
603 G4cout << "Projectile : ";
604 #endif
605 isOk = ComputeNucleusProperties( thePrNucl
606 Projectil
607 Projectil
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" << G
617 << "SumMasses, PrResidualMass and Tar
618 << PrResidualMass/GeV << " " << Targe
619 #endif
620
621 if ( SqrtS < SumMasses ) return false; // I
622
623 // Try to consider also the excitation energ
624 // possible, with the available energy; othe
625 G4double savedSumMasses = SumMasses;
626 if ( isProjectileNucleus ) {
627 SumMasses -= std::sqrt( sqr( PrResidualMas
628 SumMasses += std::sqrt( sqr( PrResidualMas
629 + PprojResidual.pe
630 }
631 SumMasses -= std::sqrt( sqr( TargetResidualM
632 SumMasses += std::sqrt( sqr( TargetResidualM
633 + PtargetResidual.pe
634
635 if ( SqrtS < SumMasses ) {
636 SumMasses = savedSumMasses;
637 if ( isProjectileNucleus ) ProjectileResid
638 TargetResidualExcitationEnergy = 0.0;
639 }
640
641 TargetResidualMass += TargetResidualExcitati
642 if ( isProjectileNucleus ) PrResidualMass +=
643
644 #ifdef debugPutOnMassShell
645 if ( isProjectileNucleus ) {
646 G4cout << "PrResidualMass ProjResidualExci
647 << ProjectileResidualExcitationEnergy <<
648 }
649 G4cout << "TargetResidualMass TargetResidual
650 << TargetResidualExcitationEnergy <<
651 << "Sum masses " << SumMasses/GeV <<
652 #endif
653
654 // Sampling of nucleons what can transfer to
655 if ( isProjectileNucleus && thePrNucleus->
656 isOk = GenerateDeltaIsobar( SqrtS, Numbe
657 TheInvolvedN
658 }
659 if ( theTargetNucleus->GetMassNumber() != 1
660 isOk = isOk && GenerateDeltaIsobar( Sqrt
661 TheI
662 }
663 if ( ! isOk ) return false;
664
665 // Now we know that it is kinematically poss
666 // of the involved nucleons (or correspondin
667 // We have to sample the kinematical variabl
668 // of the final state. The sampled kinematic
669 // Notice that the sampling of the transvers
670 // Fermi motion.
671
672 G4LorentzRotation toCms( -1*Psum.boostVector
673 G4LorentzVector Ptmp = toCms*Pprojectile;
674 if ( Ptmp.pz() <= 0.0 ) return false; // "S
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
687 G4double DcorP = 0.0;
688 if ( isProjectileNucleus ) DcorP = theParame
689 G4double DcorT = theParameters->GetDof
690 G4double AveragePt2 = theParameters->GetPt2
691 G4double maxPtSquare = theParameters->GetMax
692
693 #ifdef debugPutOnMassShell
694 if ( isProjectileNucleus ) {
695 G4cout << "Y projectileNucleus " << Yproje
696 }
697 G4cout << "Y targetNucleus " << YtargetN
698 << "Dcor " << theParameters->GetDofNu
699 << " DcorP DcorT " << DcorP << " " <<
700 #endif
701
702 G4double M2proj = M2projectile; // Initiali
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
716 NumberOfTries++;
717 if ( NumberOfTries == 100*(NumberOfTries
718 // After many tries, it is convenient
719 // AveragePt2, so that the sampled mom
720 // involved nucleons (or corresponding delta
721 // it is more likely to satisfy the mo
722 ScaleFactor /= 2.0;
723 DcorP *= ScaleFactor;
724 DcorT *= ScaleFactor;
725 AveragePt2 *= ScaleFactor;
726 }
727 if ( isProjectileNucleus ) {
728 // Sampling of kinematical properties
729 isOk = SamplingNucleonKinematics( Aver
730 theP
731 PrRe
732 Numb
733 TheI
734 }
735 // Sampling of kinematical properties of
736 isOk = isOk && SamplingNucleonKinemati
737
738
739
740
741 #ifdef debugPutOnMassShell
742 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << Sqr
743 << ( std::sqrt( M2proj ) + std::s
744 << std::sqrt( M2proj )/GeV << " "
745 #endif
746 if ( ! isOk ) return false;
747 } while ( ( SqrtS < std::sqrt( M2proj ) +
748 NumberOfTries < maxNumberOfInner
749 if ( NumberOfTries >= maxNumberOfInnerLoop
750 #ifdef debugPutOnMassShell
751 G4cout << "BAD situation: forced exit of
752 #endif
753 return false;
754 }
755 if ( isProjectileNucleus ) {
756 isOk = CheckKinematics( S, SqrtS, M2proj
757 NumberOfInvolved
758 TheInvolvedNucle
759 WminusTarget, Wp
760 }
761 isOk = isOk && CheckKinematics( S, SqrtS
762 NumberOf
763 WminusTa
764 if ( ! isOk ) return false;
765 } while ( ( ! OuterSuccess ) &&
766 ++loopCounter < maxNumberOfLoops )
767 if ( loopCounter >= maxNumberOfLoops ) {
768 #ifdef debugPutOnMassShell
769 G4cout << "BAD situation: forced exit of t
770 #endif
771 return false;
772 }
773
774 // Now the sampling is completed, and we can
775 // whole system. This is done first in the c
776 // to the lab frame. The transverse momentum
777 // the recoil of each hadron (nucleon or del
778 // to conserve (by construction) the transve
779
780 if ( ! isProjectileNucleus ) { // hadron-nu
781
782 G4double Pzprojectile = WplusProjectile/2.
783 G4double Eprojectile = WplusProjectile/2.
784 Pprojectile.setPz( Pzprojectile );
785 Pprojectile.setE( Eprojectile );
786
787 #ifdef debugPutOnMassShell
788 G4cout << "Proj after in CMS " << Pproject
789 #endif
790
791 Pprojectile.transform( toLab );
792 theProjectile.SetMomentum( Pprojectile.vec
793 theProjectile.SetTotalEnergy( Pprojectile.
794
795 theParticipants.StartLoop();
796 theParticipants.Next();
797 G4VSplitableHadron* primary = theParticipa
798 primary->Set4Momentum( Pprojectile );
799
800 #ifdef debugPutOnMassShell
801 G4cout << "Final proj. mom in Lab. " << pr
802 #endif
803
804 } else { // nucleus-nucleus or antinucleus-
805
806 isOk = FinalizeKinematics( WplusProjectile
807 ProjectileResid
808 TheInvolvedNucl
809
810 #ifdef debugPutOnMassShell
811 G4cout << "Projectile Residual4Momentum in
812 #endif
813
814 if ( ! isOk ) return false;
815
816 ProjectileResidual4Momentum.transform( toL
817
818 #ifdef debugPutOnMassShell
819 G4cout << "Projectile Residual4Momentum in
820 #endif
821
822 }
823
824 isOk = FinalizeKinematics( WminusTarget, fal
825 TargetResidualMas
826 TheInvolvedNucleo
827
828 #ifdef debugPutOnMassShell
829 G4cout << "Target Residual4Momentum in CMS "
830 #endif
831
832 if ( ! isOk ) return false;
833
834 TargetResidual4Momentum.transform( toLab );
835
836 #ifdef debugPutOnMassShell
837 G4cout << "Target Residual4Momentum in Lab "
838 #endif
839
840 return true;
841
842 }
843
844
845 //============================================
846
847 G4bool G4FTFModel::ExciteParticipants() {
848
849 #ifdef debugBuildString
850 G4cout << "G4FTFModel::ExciteParticipants()
851 #endif
852
853 G4bool Success( false );
854 G4int MaxNumOfInelCollisions = G4int( thePar
855 if ( MaxNumOfInelCollisions > 0 ) { // Pla
856 G4double ProbMaxNumber = theParameters->Ge
857 if ( G4UniformRand() < ProbMaxNumber ) Max
858 } else {
859 // Plab < Pbound, normal application of FT
860 MaxNumOfInelCollisions = 1;
861 }
862
863 #ifdef debugBuildString
864 G4cout << "MaxNumOfInelCollisions per hadron
865 #endif
866
867 G4int CurrentInteraction( 0 );
868 theParticipants.StartLoop();
869
870 G4bool InnerSuccess( true );
871 while ( theParticipants.Next() ) { /* Loop
872 CurrentInteraction++;
873 const G4InteractionContent& collision = th
874 G4VSplitableHadron* projectile = collision
875 G4Nucleon* ProjectileNucleon = collision.G
876 G4VSplitableHadron* target = collision.Get
877 G4Nucleon* TargetNucleon = collision.GetTa
878
879 #ifdef debugBuildString
880 G4cout << G4endl << "Interaction # Status
881 << collision.GetStatus() << G4endl
882 << target << G4endl << "projectile-
883 << projectile->GetStatus() << " " <
884 << "projectile->GetSoftC target->G
885 << " " << target->GetSoftCollisionC
886 #endif
887
888 if ( collision.GetStatus() ) {
889 if ( G4UniformRand() < theParameters->Ge
890 // Elastic scattering
891
892 #ifdef debugBuildString
893 G4cout << "Elastic scattering" << G4en
894 #endif
895
896 if ( ! HighEnergyInter ) {
897 G4bool Annihilation = false;
898 G4bool Result = AdjustNucleons( proj
899 Targ
900 if ( ! Result ) continue;
901 }
902 InnerSuccess = theElastic->ElasticSca
903 } else if ( G4UniformRand() > theParamet
904 // Inelastic scattering
905
906 #ifdef debugBuildString
907 G4cout << "Inelastic interaction" << G
908 << "MaxNumOfInelCollisions per
909 #endif
910
911 if ( ! HighEnergyInter ) {
912 G4bool Annihilation = false;
913 G4bool Result = AdjustNucleons( proj
914 Targ
915 if ( ! Result ) continue;
916 }
917 if ( G4UniformRand() <
918 ( 1.0 - target->GetSoftCollisionC
919 ( 1.0 - projectile->GetSoftCollis
920 //if ( ! HighEnergyInter ) {
921 // G4bool Annihilation = false;
922 // G4bool Result = AdjustNucleons(
923 //
924 // if ( ! Result ) continue;
925 //}
926 if ( theExcitation->ExciteParticipan
927 InnerSuccess = true;
928 NumberOfNNcollisions++;
929 #ifdef debugBuildString
930 G4cout << "FTF excitation Successf
931 // G4cout << "After pro " << proj
932 // << projectile->Get4Momen
933 // << "After tar " << targ
934 // << target->Get4Momentum(
935 #endif
936 } else {
937 InnerSuccess = theElastic->Elastic
938 #ifdef debugBuildString
939 G4cout << "FTF excitation Non Inne
940 << InnerSuccess << G4endl;
941 #endif
942 }
943 } else { // The inelastic interactiti
944 #ifdef debugBuildString
945 G4cout << "Elastic scat. at rejectio
946 #endif
947 //if ( ! HighEnergyInter ) {
948 // G4bool Annihilation = false;
949 // G4bool Result = AdjustNucleons(
950 //
951 // if ( ! Result) continue;
952 //}
953 InnerSuccess = theElastic->ElasticSc
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( proj
965 Targ
966 if ( ! Result ) continue;
967 }
968
969 G4VSplitableHadron* AdditionalString =
970 if ( theAnnihilation->Annihilate( proj
971 InnerSuccess = true;
972 #ifdef debugBuildString
973 G4cout << "Annihilation successfull.
974 << AdditionalString << G4endl
975 //G4cout << "After pro " << project
976 //G4cout << "After tar " << target-
977 #endif
978
979 if ( AdditionalString != 0 ) theAddi
980
981 NumberOfNNcollisions++;
982
983 // Skipping possible interactions of
984 while ( theParticipants.Next() ) {
985 G4InteractionContent& acollision =
986 G4VSplitableHadron* NextProjectile
987 G4VSplitableHadron* NextTargetNucl
988 if ( projectile == NextProjectileN
989 acollision.SetStatus( 0 );
990 }
991 }
992
993 // Continue the interactions
994 theParticipants.StartLoop();
995 for ( G4int i = 0; i < CurrentIntera
996
997 /*
998 if ( target->GetStatus() == 4 ) {
999 // Skipping possible interactions
1000 while ( theParticipants.Next() )
1001 G4InteractionContent& acollisio
1002 G4VSplitableHadron* NextProject
1003 G4VSplitableHadron* NextTargetN
1004 if ( target == NextTargetNucleo
1005 }
1006 }
1007 theParticipants.StartLoop();
1008 for ( G4int I = 0; I < CurrentInter
1009 */
1010
1011 }
1012 }
1013 }
1014
1015 if( InnerSuccess ) Success = true;
1016
1017 #ifdef debugBuildString
1018 G4cout << "-----------------------------
1019 << "projectile->GetStatus target->
1020 << " " << target->GetStatus() << G
1021 << projectile->GetSoftCollisionCou
1022 << G4endl << "ExciteParticipants()
1023 #endif
1024
1025 } // end of while ( theParticipants.Next()
1026
1027 return Success;
1028 }
1029
1030
1031 //===========================================
1032
1033 G4bool G4FTFModel::AdjustNucleons( G4VSplitab
1034 G4Nucleon*
1035 G4VSplitab
1036 G4Nucleon*
1037 G4bool Ann
1038
1039 #ifdef debugAdjust
1040 G4cout << "AdjustNucleons -----------------
1041 << "Proj is nucleus? " << GetProject
1042 << "Proj 4mom " << SelectedAntiBaryo
1043 << "Targ 4mom " << SelectedTargetNuc
1044 << "Pr ResidualMassNumber Pr Residua
1045 << ProjectileResidualMassNumber << "
1046 << ProjectileResidualExcitationEnerg
1047 << "Tr ResidualMassNumber Tr Residua
1048 << TargetResidualMassNumber << " " <
1049 << TargetResidualExcitationEnergy <<
1050 << "Collis. pr tr " << SelectedAntiB
1051 << SelectedTargetNucleon->GetSoftCol
1052 #endif
1053
1054 if ( SelectedAntiBaryon->GetSoftCollisionCo
1055 SelectedTargetNucleon->GetSoftCollisio
1056 return true; // Selected hadrons were adj
1057 }
1058
1059 G4int interactionCase = 0;
1060 if ( ( ! GetProjectileNucleus() &&
1061 SelectedAntiBaryon->GetSoftCollis
1062 SelectedTargetNucleon->GetSoftCol
1063 ||
1064 ( SelectedAntiBaryon->GetSoftCollis
1065 SelectedTargetNucleon->GetSoftCol
1066 // The case of hadron-nucleus interaction
1067 // the case when projectile nuclear nucle
1068 // a collision, but target nucleon did no
1069 interactionCase = 1;
1070 #ifdef debugAdjust
1071 G4cout << "case 1, hA prcol=0 trcol=0, AA
1072 #endif
1073 if ( TargetResidualMassNumber < 1 ) {
1074 return false;
1075 }
1076 if ( SelectedAntiBaryon->Get4Momentum().r
1077 return false;
1078 }
1079 if ( TargetResidualMassNumber == 1 ) {
1080 TargetResidualMassNumber = 0;
1081 TargetResidualCharge = 0;
1082 TargetResidualExcitationEnergy = 0.0;
1083 SelectedTargetNucleon->Set4Momentum( Ta
1084 TargetResidual4Momentum = G4LorentzVect
1085 return true;
1086 }
1087 } else if ( SelectedAntiBaryon->GetSoftColl
1088 SelectedTargetNucleon->GetSoftC
1089 // It is assumed that in the case there i
1090 interactionCase = 2;
1091 #ifdef debugAdjust
1092 G4cout << "case 2, prcol=0 trcol#0" << G
1093 #endif
1094 if ( ProjectileResidualMassNumber < 1 ) {
1095 return false;
1096 }
1097 if ( ProjectileResidual4Momentum.rapidity
1098 SelectedTargetNucleon->Get4Momentum(
1099 return false;
1100 }
1101 if ( ProjectileResidualMassNumber == 1 )
1102 ProjectileResidualMassNumber = 0;
1103 ProjectileResidualCharge = 0;
1104 ProjectileResidualExcitationEnergy = 0.
1105 SelectedAntiBaryon->Set4Momentum( Proje
1106 ProjectileResidual4Momentum = G4Lorentz
1107 return true;
1108 }
1109 } else { // It has to be a nucleus-nucleus
1110 interactionCase = 3;
1111 #ifdef debugAdjust
1112 G4cout << "case 3, prcol=0 trcol=0" << G
1113 #endif
1114 if ( ! GetProjectileNucleus() ) {
1115 return false;
1116 }
1117 #ifdef debugAdjust
1118 G4cout << "Proj res Init " << ProjectileR
1119 << "Targ res Init " << TargetResid
1120 << "ProjectileResidualMassNumber P
1121 << ProjectileResidualMassNumber <<
1122 << " (" << ProjectileResidualLambd
1123 << "TargetResidualMassNumber Targe
1124 << " " << TargetResidualCharge <<
1125 #endif
1126 }
1127
1128 CommonVariables common;
1129 G4int returnCode = AdjustNucleonsAlgorithm_
1130
1131
1132 G4bool returnResult = false;
1133 if ( returnCode == 0 ) {
1134 returnResult = true; // Successfully end
1135 } else if ( returnCode == 1 ) {
1136 // The part before sampling has been succ
1137 returnResult = AdjustNucleonsAlgorithm_Sa
1138 if ( returnResult ) { // The sampling ha
1139 AdjustNucleonsAlgorithm_afterSampling(
1140
1141 }
1142 }
1143
1144 return returnResult;
1145 }
1146
1147 //-------------------------------------------
1148
1149 G4int G4FTFModel::AdjustNucleonsAlgorithm_bef
1150
1151
1152
1153
1154
1155
1156 // First of the three utility methods used
1157 // This method returns an integer code - in
1158 // "0" : successfully ended and nothing e
1159 // "1" : successfully completed, but the
1160 // "99" : unsuccessfully ended, nothing el
1161 G4int returnCode = 99;
1162
1163 G4double ExcitationEnergyPerWoundedNucleon
1164
1165 // some checks and initializations
1166 if ( interactionCase == 1 ) {
1167 common.Psum = SelectedAntiBaryon->Get4Mom
1168 #ifdef debugAdjust
1169 G4cout << "Targ res Init " << TargetResid
1170 #endif
1171 common.Pprojectile = SelectedAntiBaryon->
1172 } else if ( interactionCase == 2 ) {
1173 common.Psum = ProjectileResidual4Momentum
1174 common.Pprojectile = ProjectileResidual4M
1175 } else if ( interactionCase == 3 ) {
1176 common.Psum = ProjectileResidual4Momentum
1177 common.Pprojectile = ProjectileResidual4M
1178 }
1179
1180 // transform momenta to cms and then rotate
1181 common.toCms = G4LorentzRotation( -1*common
1182 common.Ptmp = common.toCms * common.Pprojec
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 an
1191 G4bool Stopping = false;
1192 if ( interactionCase == 1 ) {
1193 common.TResidualMassNumber = TargetResidu
1194 common.TResidualCharge = TargetResidual
1195 - G4int( TargetN
1196 common.TResidualExcitationEnergy = Targ
1197 - Exci
1198 if ( common.TResidualMassNumber <= 1 ) {
1199 common.TResidualExcitationEnergy = 0.0;
1200 }
1201 if ( common.TResidualMassNumber != 0 ) {
1202 common.TResidualMass = G4ParticleTable:
1203 ->GetIonMass( co
1204 }
1205 common.TNucleonMass = TargetNucleon->GetD
1206 common.SumMasses = SelectedAntiBaryon->
1207 + common.TResidualMass
1208 #ifdef debugAdjust
1209 G4cout << "Annihilation " << Annihilation
1210 #endif
1211 } else if ( interactionCase == 2 ) {
1212 common.Ptarget = common.toCms * SelectedT
1213 common.TResidualMassNumber = ProjectileRe
1214 common.TResidualCharge = ProjectileResi
1215 - std::abs( G4in
1216 common.TResidualExcitationEnergy = Proj
1217 - Exci
1218 if ( common.TResidualMassNumber <= 1 ) {
1219 common.TResidualExcitationEnergy = 0.0;
1220 }
1221 if ( common.TResidualMassNumber != 0 ) {
1222 common.TResidualMass = G4ParticleTable:
1223 ->GetIonMass( co
1224 }
1225 common.TNucleonMass = ProjectileNucleon->
1226 common.SumMasses = SelectedTargetNucleo
1227 + common.TResidualMass
1228 #ifdef debugAdjust
1229 G4cout << "SelectedTN.mag() PNMass + PRes
1230 << SelectedTargetNucleon->Get4Mome
1231 << common.TNucleonMass << " " << c
1232 #endif
1233 } else if ( interactionCase == 3 ) {
1234 common.Ptarget = common.toCms * TargetRes
1235 common.PResidualMassNumber = ProjectileRe
1236 common.PResidualCharge = ProjectileResi
1237 - std::abs( G4in
1238 common.PResidualLambdaNumber = Projectile
1239 if ( ProjectileNucleon->GetDefinition() =
1240 ProjectileNucleon->GetDefinition() == G4An
1241 --common.PResidualLambdaNumber;
1242 }
1243 common.PResidualExcitationEnergy = Proj
1244 - Exci
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
1251 common.PResidualMass = G4Proton::De
1252 } else if ( common.PResidualLambdaNum
1253 common.PResidualMass = G4Lambda::De
1254 } else {
1255 common.PResidualMass = G4Neutron::D
1256 }
1257 } else {
1258 if ( common.PResidualLambdaNumber > 0
1259 if ( common.PResidualMassNumber ==
1260 common.PResidualMass = G4Lambda::
1261 if ( std::abs( common.PResidualCharge )
1262 common.PResidualMass += G4Proton::Def
1263 } else if ( common.PResidualLambdaNumbe
1264 common.PResidualMass += G4Neutron::De
1265 } else {
1266 common.PResidualMass += G4Lambda::Def
1267 }
1268 } else {
1269 common.PResidualMass = G4HyperNucleiPro
1270 std::abs( common.PRes
1271 common.PResidualLambdaN
1272 }
1273 } else {
1274 common.PResidualMass = G4ParticleTable::G
1275 GetIonMass( std::a
1276 }
1277 }
1278 }
1279 common.PNucleonMass = ProjectileNucleon->
1280 common.TResidualMassNumber = TargetResidu
1281 common.TResidualCharge = TargetResidual
1282 - G4int( TargetN
1283 common.TResidualExcitationEnergy = Targ
1284 - Exci
1285 if ( common.TResidualMassNumber <= 1 ) {
1286 common.TResidualExcitationEnergy = 0.0;
1287 }
1288 if ( common.TResidualMassNumber != 0 ) {
1289 common.TResidualMass = G4ParticleTable:
1290 GetIonMass( comm
1291 }
1292 common.TNucleonMass = TargetNucleon->GetD
1293 common.SumMasses = common.PNucleonMass
1294 + common.TResidualMass
1295 #ifdef debugAdjust
1296 G4cout << "PNucleonMass PResidualMass TNu
1297 << " " << common.PResidualMass <<
1298 << common.TResidualMass << G4endl
1299 << "PResidualExcitationEnergy " <<
1300 << "TResidualExcitationEnergy " <<
1301 #endif
1302 } // End-if on interactionCase
1303
1304 if ( ! Annihilation ) {
1305 if ( common.SqrtS < common.SumMasses ) {
1306 #ifdef debugAdjust
1307 G4cout << "SqrtS < SumMasses " << commo
1308 #endif
1309 return returnCode; // Unsuccessfully e
1310 }
1311 if ( interactionCase == 1 || interactio
1312 if ( common.SqrtS < common.SumMasses +
1313 #ifdef debugAdjust
1314 G4cout << "TResidualExcitationEnergy
1315 #endif
1316 common.TResidualExcitationEnergy = co
1317 #ifdef debugAdjust
1318 G4cout << "TResidualExcitationEnergy
1319 #endif
1320 Stopping = true;
1321 return returnCode; // unsuccessfully
1322 }
1323 } else if ( interactionCase == 3 ) {
1324 #ifdef debugAdjust
1325 G4cout << "SqrtS < SumMasses + PResidua
1326 << common.SqrtS << " " << common
1327 << G4endl;
1328 #endif
1329 if ( common.SqrtS < common.SumMasses
1330 + common.TResidualE
1331 Stopping = true;
1332 if ( common.PResidualExcitationEnergy
1333 common.TResidualExcitationEnergy =
1334 } else if ( common.TResidualExcitatio
1335 common.PResidualExcitationEnergy =
1336 } else {
1337 G4double Fraction = ( common.SqrtS
1338 / ( common.PResidualExcitationEn
1339 common.PResidualExcitationEnergy *=
1340 common.TResidualExcitationEnergy *=
1341 }
1342 }
1343 }
1344 } // End-if on ! Annihilation
1345
1346 if ( Annihilation ) {
1347 #ifdef debugAdjust
1348 G4cout << "SqrtS < SumMasses - TNucleonMa
1349 << common.SumMasses - common.TNucl
1350 #endif
1351 if ( common.SqrtS < common.SumMasses - co
1352 return returnCode; // unsuccessfully e
1353 }
1354 #ifdef debugAdjust
1355 G4cout << "SqrtS < SumMasses " << common.
1356 #endif
1357 if ( common.SqrtS < common.SumMasses ) {
1358 if ( interactionCase == 2 || interact
1359 common.TResidualExcitationEnergy = 0.
1360 }
1361 common.TNucleonMass = common.SqrtS -
1362 - common.TResidua
1363 #ifdef debugAdjust
1364 G4cout << "TNucleonMass " << common.TNu
1365 #endif
1366 common.SumMasses = common.SqrtS - commo
1367 Stopping = true;
1368 #ifdef debugAdjust
1369 G4cout << "SqrtS < SumMasses " << commo
1370 #endif
1371 }
1372 if ( interactionCase == 1 || interactio
1373 if ( common.SqrtS < common.SumMasses +
1374 common.TResidualExcitationEnergy = co
1375 Stopping = true;
1376 }
1377 } else if ( interactionCase == 3 ) {
1378 if ( common.SqrtS < common.SumMasses
1379 + common.TResidualE
1380 Stopping = true;
1381 if ( common.PResidualExcitationEnergy
1382 common.TResidualExcitationEnergy =
1383 } else if ( common.TResidualExcitatio
1384 common.PResidualExcitationEnergy =
1385 } else {
1386 G4double Fraction = ( common.SqrtS
1387 ( common.PResidualExcitationEnerg
1388 common.PResidualExcitationEnergy *=
1389 common.TResidualExcitationEnergy *=
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.set
1402 // New projectile
1403 if ( interactionCase == 1 ) {
1404 common.Ptmp.setE( SelectedAntiBaryon->G
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 <<
1412 #endif
1413 common.Pprojectile = common.Ptmp;
1414 common.Pprojectile.transform( common.toLa
1415 //---AR-Jul2019 : To avoid unphysical pro
1416 // original momentum of th
1417 G4LorentzVector saveSelectedAntiBaryon4Mo
1418 saveSelectedAntiBaryon4Momentum.transform
1419 //---
1420 SelectedAntiBaryon->Set4Momentum( common.
1421 // New target nucleon
1422 if ( interactionCase == 1 || interactio
1423 common.Ptmp.setE( common.TNucleonMass )
1424 } else if ( interactionCase == 2 ) {
1425 common.Ptmp.setE( SelectedTargetNucleon
1426 }
1427 #ifdef debugAdjust
1428 G4cout << "Targ stop " << common.Ptmp <<
1429 #endif
1430 common.Ptarget = common.Ptmp;
1431 common.Ptarget.transform( common.toLab );
1432 //---AR-Jul2019 : To avoid unphysical tar
1433 // momentum of the target
1434 G4LorentzVector saveSelectedTargetNucleon
1435 saveSelectedTargetNucleon4Momentum.transf
1436 //---
1437 SelectedTargetNucleon->Set4Momentum( comm
1438 // New target residual
1439 if ( interactionCase == 1 || interactio
1440 common.Ptmp.setPx( 0.0 ); common.Ptmp.s
1441 TargetResidualMassNumber = common
1442 TargetResidualCharge = common
1443 TargetResidualExcitationEnergy = common
1444 //---AR-Jul2019 : To avoid unphysical t
1445 // original momentum of
1446 // This is a rough and s
1447 //common.Ptmp.setE( common.TResidualMas
1448 common.Ptmp.setPx( -saveSelectedTargetN
1449 common.Ptmp.setPy( -saveSelectedTargetN
1450 common.Ptmp.setPz( -saveSelectedTargetN
1451 common.Ptmp.setE( std::sqrt( sqr( commo
1452 //---
1453 #ifdef debugAdjust
1454 G4cout << "Targ Resi stop " << common.P
1455 #endif
1456 common.Ptmp.transform( common.toLab );
1457 TargetResidual4Momentum = common.Ptmp;
1458 }
1459 // New projectile residual
1460 if ( interactionCase == 2 || interactio
1461 common.Ptmp.setPx( 0.0 ); common.Ptmp.s
1462 if ( interactionCase == 2 ) {
1463 ProjectileResidualMassNumber =
1464 ProjectileResidualCharge =
1465 ProjectileResidualLambdaNumber = 0; //
1466 ProjectileResidualExcitationEnergy =
1467 common.Ptmp.setE( common.TResidualMas
1468 } else {
1469 ProjectileResidualMassNumber =
1470 ProjectileResidualCharge =
1471 ProjectileResidualLambdaNumber = common
1472 ProjectileResidualExcitationEnergy =
1473 //---AR-Jul2019 : To avoid unphysical
1474 // saved original mome
1475 // This is a rough and
1476 //common.Ptmp.setE( common.PResidualM
1477 common.Ptmp.setPx( -saveSelectedAntiB
1478 common.Ptmp.setPy( -saveSelectedAntiB
1479 common.Ptmp.setPz( -saveSelectedAntiB
1480 common.Ptmp.setE( std::sqrt( sqr( com
1481 //---
1482 }
1483 #ifdef debugAdjust
1484 G4cout << "Proj Resi stop " << common.P
1485 #endif
1486 common.Ptmp.transform( common.toLab );
1487 ProjectileResidual4Momentum = common.Pt
1488 }
1489 return returnCode = 0; // successfully e
1490 } // End-if on Stopping
1491
1492 // Initializations before sampling
1493 if ( interactionCase == 1 ) {
1494 common.Mprojectile = common.Pprojectile.
1495 common.M2projectile = common.Pprojectile.
1496 common.TResidual4Momentum = common.toCms
1497 common.YtargetNucleus = common.TResidual4
1498 common.TResidualMass += common.TResidualE
1499 } else if ( interactionCase == 2 ) {
1500 common.Mtarget = common.Ptarget.mag();
1501 common.M2target = common.Ptarget.mag2();
1502 common.TResidual4Momentum = common.toCms
1503 common.YprojectileNucleus = common.TResid
1504 common.TResidualMass += common.TResidualE
1505 } else if ( interactionCase == 3 ) {
1506 common.PResidual4Momentum = common.toCms
1507 common.YprojectileNucleus = common.PResid
1508 common.TResidual4Momentum = common.toCms*
1509 common.YtargetNucleus = common.TResidual4
1510 common.PResidualMass += common.PResidualE
1511 common.TResidualMass += common.TResidualE
1512 }
1513 #ifdef debugAdjust
1514 G4cout << "YprojectileNucleus " << common.Y
1515 #endif
1516
1517 return returnCode = 1; // successfully com
1518 }
1519
1520
1521 //-------------------------------------------
1522
1523 G4bool G4FTFModel::AdjustNucleonsAlgorithm_Sa
1524
1525 // Second of the three utility methods used
1526 // This method returns "false" if it fails
1527
1528 // Ascribing of the involved nucleons Pt an
1529 G4double Dcor = theParameters->GetDofNuclea
1530 G4double DcorP = 0.0, DcorT = 0.0;
1531 if ( ProjectileResidualMassNumber != 0 ) Dc
1532 if ( TargetResidualMassNumber != 0 ) Dc
1533 G4double AveragePt2 = theParameters->GetPt2
1534 G4double maxPtSquare = theParameters->GetMa
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*(NumberOfTrie
1547 // At large number of tries it would
1548 ScaleFactor /= 2.0;
1549 DcorP *= ScaleFactor;
1550 DcorT *= ScaleFactor;
1551 AveragePt2 *= ScaleFactor;
1552 #ifdef debugAdjust
1553 //G4cout << "NumberOfTries ScaleFacto
1554 #endif
1555 }
1556
1557 // Some kinematics
1558 if ( interactionCase == 1 ) {
1559 } else if ( interactionCase == 2 ) {
1560 #ifdef debugAdjust
1561 G4cout << "ProjectileResidualMassNumb
1562 #endif
1563 if ( ProjectileResidualMassNumber > 1
1564 common.PtNucleon = GaussianPt( Aver
1565 } else {
1566 common.PtNucleon = G4ThreeVector( 0
1567 }
1568 common.PtResidual = - common.PtNucleo
1569 common.Mprojectile = std::sqrt( sqr
1570 + std::sqrt( sqr
1571 #ifdef debugAdjust
1572 G4cout << "SqrtS < Mtarget + Mproject
1573 << " " << common.Mprojectile <
1574 #endif
1575 common.M2projectile = sqr( common.Mpr
1576 if ( common.SqrtS < common.Mtarget +
1577 OuterSuccess = false;
1578 loopCondition = true;
1579 continue;
1580 }
1581 } else if ( interactionCase == 3 ) {
1582 if ( ProjectileResidualMassNumber > 1
1583 common.PtNucleonP = GaussianPt( Ave
1584 } else {
1585 common.PtNucleonP = G4ThreeVector(
1586 }
1587 common.PtResidualP = - common.PtNucle
1588 if ( TargetResidualMassNumber > 1 ) {
1589 common.PtNucleonT = GaussianPt( Ave
1590 } else {
1591 common.PtNucleonT = G4ThreeVector(
1592 }
1593 common.PtResidualT = - common.PtNucle
1594 common.Mprojectile = std::sqrt( sqr
1595 + std::sqrt( sqr
1596 common.M2projectile = sqr( common.Mpr
1597 common.Mtarget = std::sqrt( sqr( co
1598 + std::sqrt( sqr( co
1599 common.M2target = sqr( common.Mtarget
1600 if ( common.SqrtS < common.Mprojectil
1601 OuterSuccess = false;
1602 loopCondition = true;
1603 continue;
1604 }
1605 } // End-if on interactionCase
1606
1607 G4int numberOfTimesExecuteInnerLoop = 1
1608 if ( interactionCase == 3 ) numberOfTim
1609 for ( G4int iExecute = 0; iExecute < nu
1610
1611 G4bool InnerSuccess = true;
1612 G4bool isTargetToBeHandled = ( intera
1613 ( inte
1614 G4bool condition = false;
1615 if ( isTargetToBeHandled ) {
1616 condition = ( TargetResidualMassNum
1617 } else { // Projectile to be handled
1618 condition = ( ProjectileResidualMas
1619 }
1620 if ( condition ) {
1621 const G4int maxNumberOfInnerLoops =
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
1629 common.PtResidual = - common.
1630 common.Mtarget = std::sqrt(
1631 + std::sqrt(
1632 if ( common.SqrtS < common.Mp
1633 InnerSuccess = false;
1634 continue;
1635 }
1636 Xcenter = std::sqrt( sqr( com
1637 / common.Mtarget;
1638 } else {
1639 Xcenter = std::sqrt( sqr( com
1640 / common.Mtarget;
1641 }
1642 G4ThreeVector tmpX = GaussianPt
1643 common.XminusNucleon = Xcenter
1644 if ( common.XminusNucleon <= 0.
1645 InnerSuccess = false;
1646 continue;
1647 }
1648 common.XminusResidual = 1.0 - c
1649 } else { // Projectile to be han
1650 G4ThreeVector tmpX = GaussianPt
1651 G4double Xcenter = 0.0;
1652 if ( interactionCase == 2 ) {
1653 Xcenter = std::sqrt( sqr( com
1654 / common.Mprojectil
1655 } else {
1656 Xcenter = std::sqrt( sqr( com
1657 / common.Mprojectil
1658 }
1659 common.XplusNucleon = Xcenter +
1660 if ( common.XplusNucleon <= 0.0
1661 InnerSuccess = false;
1662 continue;
1663 }
1664 common.XplusResidual = 1.0 - co
1665 } // End-if on isTargetToBeHandl
1666 } while ( ( ! InnerSuccess ) &&
1667 ++innerLoopCounter < ma
1668 if ( innerLoopCounter >= maxNumberO
1669 #ifdef debugAdjust
1670 G4cout << "BAD situation: forced
1671 #endif
1672 return false;
1673 }
1674 } else { // condition is false
1675 if ( isTargetToBeHandled ) {
1676 common.XminusNucleon = 1.0;
1677 common.XminusResidual = 1.0; //
1678 } else { // Projectile to be handl
1679 common.XplusNucleon = 1.0;
1680 common.XplusResidual = 1.0; //
1681 }
1682 } // End-if on condition
1683
1684 } // End of for loop on iExecute
1685
1686 if ( interactionCase == 1 ) {
1687 common.M2target = ( sqr( common.TN
1688 / common.XminusN
1689 + ( sqr( common.TR
1690 / common.XminusR
1691 loopCondition = ( common.SqrtS < comm
1692 } else if ( interactionCase == 2 ) {
1693 #ifdef debugAdjust
1694 G4cout << "TNucleonMass PtNucleon Xpl
1695 << common.PtNucleon << " " <<
1696 << "TResidualMass PtResidual X
1697 << common.PtResidual << " " <
1698 #endif
1699 common.M2projectile = ( sqr( commo
1700 / common.Xpl
1701 + ( sqr( commo
1702 / common.Xpl
1703 #ifdef debugAdjust
1704 G4cout << "SqrtS < Mtarget + std::sqr
1705 << common.Mtarget << " " << st
1706 << common.Mtarget + std::sqrt(
1707 #endif
1708 loopCondition = ( common.SqrtS < comm
1709 } else if ( interactionCase == 3 ) {
1710 #ifdef debugAdjust
1711 G4cout << "PtNucleonP " << common.PtN
1712 << "XplusNucleon XplusResidual
1713 << " " << common.XplusResidual
1714 << "PtNucleonT " << common.PtN
1715 << "XminusNucleon XminusResidu
1716 << " " << common.XminusResidua
1717 #endif
1718 common.M2projectile = ( sqr( common
1719 / common.Xplu
1720 + ( sqr( common
1721 / common.Xplu
1722 common.M2target = ( sqr( common.TN
1723 / common.XminusN
1724 + ( sqr( common.TR
1725 / common.XminusR
1726 loopCondition = ( common.SqrtS < (
1727 + std::sqrt( common.M2target ) )
1728 } // End-if on interactionCase
1729
1730 } while ( loopCondition &&
1731 ++NumberOfTries < maxNumberOfTr
1732 if ( NumberOfTries >= maxNumberOfTries )
1733 #ifdef debugAdjust
1734 G4cout << "BAD situation: forced exit o
1735 #endif
1736 return false;
1737 }
1738
1739 // kinematics
1740 G4double Yprojectile = 0.0, YprojectileNu
1741 G4double DecayMomentum2 = sqr( common.S )
1742 - 2.0 * ( commo
1743 + com
1744 if ( interactionCase == 1 ) {
1745 common.WminusTarget = ( common.S - co
1746 + std::sqrt( De
1747 common.WplusProjectile = common.SqrtS -
1748 common.Pzprojectile = common.WplusPro
1749 - common.M2projec
1750 common.Eprojectile = common.WplusPro
1751 + common.M2projec
1752 Yprojectile = 0.5 * G4Log( ( common.
1753 / ( common.
1754 #ifdef debugAdjust
1755 G4cout << "DecayMomentum2 " << DecayMom
1756 << "WminusTarget WplusProjectile
1757 << " " << common.WplusProjectile
1758 << "Yprojectile " << Yprojectile
1759 #endif
1760 common.Mt2targetNucleon = sqr( common.T
1761 common.PztargetNucleon = - common.Wminu
1762 + common.Mt2ta
1763 / ( 2.0 * co
1764 common.EtargetNucleon = common.Wminus
1765 + common.Mt2tar
1766 / ( 2.0 * com
1767 YtargetNucleon = 0.5 * G4Log( ( commo
1768 / ( commo
1769 #ifdef debugAdjust
1770 G4cout << "YtN Ytr YtN-Ytr " << " " <<
1771 << " " << YtargetNucleon - commo
1772 << "YtN Ypr YtN-Ypr " << " " <<
1773 << " " << YtargetNucleon - Yproj
1774 #endif
1775 if ( std::abs( YtargetNucleon - common.
1776 Yprojectile < YtargetNucleon ) {
1777 OuterSuccess = false;
1778 continue;
1779 }
1780 } else if ( interactionCase == 2 ) {
1781 common.WplusProjectile = ( common.S +
1782 + std::sqrt(
1783 common.WminusTarget = common.SqrtS - co
1784 common.Pztarget = - common.WminusTarget
1785 common.Etarget = common.WminusTarget
1786 Ytarget = 0.5 * G4Log( ( common.Etarg
1787 / ( common.Etarg
1788 #ifdef debugAdjust
1789 G4cout << "DecayMomentum2 " << DecayMom
1790 << "WminusTarget WplusProjectile
1791 << " " << common.WplusProjectile
1792 << "Ytarget " << Ytarget << G4en
1793 #endif
1794 common.Mt2projectileNucleon = sqr( comm
1795 common.PzprojectileNucleon = common.W
1796 - common.M
1797 / ( 2.0
1798 common.EprojectileNucleon = common.W
1799 + common.M
1800 / ( 2.0
1801 YprojectileNucleon = 0.5 * G4Log( ( c
1802 / ( c
1803 #ifdef debugAdjust
1804 G4cout << "YpN Ypr YpN-Ypr " << " " <<
1805 << " " << YprojectileNucleon - c
1806 << "YpN Ytr YpN-Ytr " << " " <<
1807 << " " << YprojectileNucleon - Y
1808 #endif
1809 if ( std::abs( YprojectileNucleon - com
1810 Ytarget > YprojectileNucleon ) {
1811 OuterSuccess = false;
1812 continue;
1813 }
1814 } else if ( interactionCase == 3 ) {
1815 common.WplusProjectile = ( common.S +
1816 + std::sqrt(
1817 common.WminusTarget = common.SqrtS - co
1818 common.Mt2projectileNucleon = sqr( comm
1819 common.PzprojectileNucleon = common.W
1820 - common.M
1821 / ( 2.0
1822 common.EprojectileNucleon = common.W
1823 + common.M
1824 / ( 2.0
1825 YprojectileNucleon = 0.5 * G4Log( ( c
1826 / ( c
1827 common.Mt2targetNucleon = sqr( common.T
1828 common.PztargetNucleon = - common.Wminu
1829 + common.Mt2ta
1830 / ( 2.0 * co
1831 common.EtargetNucleon = common.Wminu
1832 + common.Mt2ta
1833 / ( 2.0 * co
1834 YtargetNucleon = 0.5 * G4Log( ( commo
1835 / ( commo
1836 if ( std::abs( YtargetNucleon - common.
1837 std::abs( YprojectileNucleon - com
1838 YprojectileNucleon < YtargetNucleo
1839 OuterSuccess = false;
1840 continue;
1841 }
1842 } // End-if on interactionCase
1843
1844 } while ( ( ! OuterSuccess ) &&
1845 ++loopCounter < maxNumberOfLoops
1846 if ( loopCounter >= maxNumberOfLoops ) {
1847 #ifdef debugAdjust
1848 G4cout << "BAD situation: forced exit of
1849 #endif
1850 return false;
1851 }
1852
1853 return true;
1854 }
1855
1856 //-------------------------------------------
1857
1858 void G4FTFModel::AdjustNucleonsAlgorithm_afte
1859
1860
1861
1862 // Third of the three utility methods used
1863 // and transform back.
1864
1865 // New projectile
1866 if ( interactionCase == 1 ) {
1867 common.Pprojectile.setPz( common.Pzprojec
1868 common.Pprojectile.setE( common.Eprojecti
1869 } else if ( interactionCase == 2 ) {
1870 common.Pprojectile.setPx( common.PtNucleo
1871 common.Pprojectile.setPy( common.PtNucleo
1872 common.Pprojectile.setPz( common.Pzprojec
1873 common.Pprojectile.setE( common.Eprojecti
1874 } else if ( interactionCase == 3 ) {
1875 common.Pprojectile.setPx( common.PtNucleo
1876 common.Pprojectile.setPy( common.PtNucleo
1877 common.Pprojectile.setPz( common.Pzprojec
1878 common.Pprojectile.setE( common.Eprojecti
1879 }
1880 #ifdef debugAdjust
1881 G4cout << "Proj after in CMS " << common.Pp
1882 #endif
1883 common.Pprojectile.transform( common.toLab
1884 SelectedAntiBaryon->Set4Momentum( common.Pp
1885 #ifdef debugAdjust
1886 G4cout << "Proj after in Lab " << common.Pp
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.PztargetNucl
1894 common.Ptarget.setE( common.EtargetNucleo
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.PztargetNucl
1902 common.Ptarget.setE( common.EtargetNucleo
1903 }
1904 #ifdef debugAdjust
1905 G4cout << "Targ after in CMS " << common.Pt
1906 #endif
1907 common.Ptarget.transform( common.toLab );
1908 SelectedTargetNucleon->Set4Momentum( common
1909 #ifdef debugAdjust
1910 G4cout << "Targ after in Lab " << common.Pt
1911 #endif
1912
1913 // New target residual
1914 if ( interactionCase == 1 || interactionC
1915 TargetResidualMassNumber = common.T
1916 TargetResidualCharge = common.T
1917 TargetResidualExcitationEnergy = common.T
1918 #ifdef debugAdjust
1919 G4cout << "TargetResidualMassNumber Targe
1920 << TargetResidualMassNumber << " "
1921 << TargetResidualExcitationEnergy
1922 #endif
1923 if ( TargetResidualMassNumber != 0 ) {
1924 G4double Mt2 = 0.0;
1925 if ( interactionCase == 1 ) {
1926 Mt2 = sqr( common.TResidualMass ) + c
1927 TargetResidual4Momentum.setPx( common
1928 TargetResidual4Momentum.setPy( common
1929 } else { // interactionCase == 3
1930 Mt2 = sqr( common.TResidualMass ) + c
1931 TargetResidual4Momentum.setPx( common
1932 TargetResidual4Momentum.setPy( common
1933 }
1934 G4double Pz = - common.WminusTarget * c
1935 + Mt2 / ( 2.0 * common.Wm
1936 G4double E = common.WminusTarget * c
1937 + Mt2 / ( 2.0 * common.Wm
1938 TargetResidual4Momentum.setPz( Pz );
1939 TargetResidual4Momentum.setE( E ) ;
1940 TargetResidual4Momentum.transform( comm
1941 } else {
1942 TargetResidual4Momentum = G4LorentzVect
1943 }
1944 #ifdef debugAdjust
1945 G4cout << "Tr N R " << common.Ptarget <<
1946 #endif
1947 }
1948
1949 // New projectile residual
1950 if ( interactionCase == 2 || interactionC
1951 if ( interactionCase == 2 ) {
1952 ProjectileResidualMassNumber = co
1953 ProjectileResidualCharge = co
1954 ProjectileResidualExcitationEnergy = co
1955 ProjectileResidualLambdaNumber = co
1956 } else { // interactionCase == 3
1957 ProjectileResidualMassNumber = co
1958 ProjectileResidualCharge = co
1959 ProjectileResidualExcitationEnergy = co
1960 ProjectileResidualLambdaNumber = co
1961 }
1962 #ifdef debugAdjust
1963 G4cout << "ProjectileResidualMassNumber P
1964 << ProjectileResidualMassNumber <<
1965 << ProjectileResidualLambdaNumber
1966 << ProjectileResidualExcitationEne
1967 #endif
1968 if ( ProjectileResidualMassNumber != 0 )
1969 G4double Mt2 = 0.0;
1970 if ( interactionCase == 2 ) {
1971 Mt2 = sqr( common.TResidualMass ) + c
1972 ProjectileResidual4Momentum.setPx( co
1973 ProjectileResidual4Momentum.setPy( co
1974 } else { // interactionCase == 3
1975 Mt2 = sqr( common.PResidualMass ) + c
1976 ProjectileResidual4Momentum.setPx( co
1977 ProjectileResidual4Momentum.setPy( co
1978 }
1979 G4double Pz = common.WplusProjectile
1980 - Mt2 / ( 2.0 * common.Wp
1981 G4double E = common.WplusProjectile
1982 + Mt2 / ( 2.0 * common.Wp
1983 ProjectileResidual4Momentum.setPz( Pz )
1984 ProjectileResidual4Momentum.setE( E );
1985 ProjectileResidual4Momentum.transform(
1986 } else {
1987 ProjectileResidual4Momentum = G4Lorentz
1988 }
1989 #ifdef debugAdjust
1990 G4cout << "Pr N R " << common.Pprojectile
1991 << " " << ProjectileResidual
1992 #endif
1993 }
1994
1995 }
1996
1997
1998 //===========================================
1999
2000 void G4FTFModel::BuildStrings( G4ExcitedStrin
2001 // Loop over all collisions; find all prima
2002 // (targets may be duplicate in the List (t
2003
2004 G4ExcitedString* FirstString( 0 ); // I
2005 G4ExcitedString* SecondString( 0 ); // t
2006
2007 if ( ! GetProjectileNucleus() ) {
2008
2009 std::vector< G4VSplitableHadron* > primar
2010 theParticipants.StartLoop();
2011 while ( theParticipants.Next() ) { /* Lo
2012 const G4InteractionContent& interaction
2013 // do not allow for duplicates ...
2014 if ( interaction.GetStatus() ) {
2015 if ( primaries.end() == std::find( pr
2016 in
2017 primaries.push_back( interaction.Ge
2018 }
2019 }
2020 }
2021
2022 #ifdef debugBuildString
2023 G4cout << "G4FTFModel::BuildStrings()" <<
2024 << "Number of projectile strings "
2025 #endif
2026
2027 for ( unsigned int ahadron = 0; ahadron <
2028 G4bool isProjectile( true );
2029 //G4cout << "primaries[ ahadron ] " <<
2030 //if ( primaries[ ahadron ]->GetStatus(
2031 FirstString = 0; SecondString = 0;
2032 if ( primaries[ahadron]->GetStatus() ==
2033 theExcitation->CreateStrings( primarie
2034 FirstStr
2035 NumberOfProjectileSpectatorNucleons--;
2036 } else if ( primaries[ahadron]->GetStat
2037 && primaries[ahadron]->GetSoft
2038 theExcitation->CreateStrings( primarie
2039 FirstStr
2040 NumberOfProjectileSpectatorNucleons--;
2041 } else if ( primaries[ahadron]->GetStat
2042 && primaries[ahadron]->GetSoft
2043 G4LorentzVector ParticleMomentum=prima
2044 G4KineticTrack* aTrack = new G4Kinetic
2045
2046
2047
2048 FirstString = new G4ExcitedString( aTr
2049 } else if (primaries[ahadron]->GetStatu
2050 G4LorentzVector ParticleMomentum=prima
2051 G4KineticTrack* aTrack = new G4Kinetic
2052
2053
2054
2055 FirstString = new G4ExcitedString( aTr
2056 NumberOfProjectileSpectatorNucleons--;
2057 } else {
2058 G4cout << "Something wrong in FTF Mod
2059 }
2060
2061 if ( FirstString != 0 ) strings->push_
2062 if ( SecondString != 0 ) strings->push_
2063
2064 #ifdef debugBuildString
2065 G4cout << "FirstString & SecondString?
2066 if ( FirstString->IsExcited() ) {
2067 G4cout << "Quarks on the FirstString
2068 << " " << FirstString->GetLeft
2069 } else {
2070 G4cout << "Kinetic track is stored" <
2071 }
2072 #endif
2073
2074 }
2075
2076 #ifdef debugBuildString
2077 if ( FirstString->IsExcited() ) {
2078 G4cout << "Check 1 string " << strings-
2079 << " " << strings->operator[](0)
2080 }
2081 #endif
2082
2083 std::for_each( primaries.begin(), primari
2084 primaries.clear();
2085
2086 } else { // Projectile is a nucleus
2087
2088 #ifdef debugBuildString
2089 G4cout << "Building of projectile-like st
2090 #endif
2091
2092 G4bool isProjectile = true;
2093 for ( G4int ahadron = 0; ahadron < Number
2094
2095 #ifdef debugBuildString
2096 G4cout << "Nucleon #, status, intCount
2097 << TheInvolvedNucleonsOfProjecti
2098 << " " << TheInvolvedNucleonsOfP
2099 ->GetSoftCollisionCoun
2100 #endif
2101
2102 G4VSplitableHadron* aProjectile =
2103 TheInvolvedNucleonsOfProjectile[ ah
2104
2105 #ifdef debugBuildString
2106 G4cout << G4endl << "ahadron aProjectil
2107 << " " << aProjectile->GetStatus
2108 #endif
2109
2110 FirstString = 0; SecondString = 0;
2111 if ( aProjectile->GetStatus() == 0 ) {
2112
2113 #ifdef debugBuildString
2114 G4cout << "Case1 aProjectile->GetStat
2115 #endif
2116
2117 theExcitation->CreateStrings(
2118 TheInvolvedNucleon
2119 isProjectile, Firs
2120 NumberOfProjectileSpectatorNucleons--
2121 } else if ( aProjectile->GetStatus() ==
2122 // Nucleon took part in diffractive i
2123
2124 #ifdef debugBuildString
2125 G4cout << "Case2 aProjectile->GetStat
2126 #endif
2127
2128 theExcitation->CreateStrings(
2129 TheInvolvedNucleon
2130 isProjectile, Firs
2131 NumberOfProjectileSpectatorNucleons--
2132 } else if ( aProjectile->GetStatus() ==
2133 HighEnergyInter ) {
2134 // Nucleon was considered as a parici
2135 // but the interaction was skipped du
2136 // It is now considered as an involve
2137
2138 #ifdef debugBuildString
2139 G4cout << "Case3 aProjectile->GetStat
2140 #endif
2141
2142 G4LorentzVector ParticleMomentum = aP
2143 G4KineticTrack* aTrack = new G4Kineti
2144
2145
2146
2147 FirstString = new G4ExcitedString( aT
2148
2149 #ifdef debugBuildString
2150 G4cout << " Strings are built for nuc
2151 << " the interaction was skipp
2152 #endif
2153
2154 } else if ( aProjectile->GetStatus() ==
2155 // Nucleon which was involved in the
2156
2157 #ifdef debugBuildString
2158 G4cout << "Case4 aProjectile->GetStat
2159 #endif
2160
2161 G4LorentzVector ParticleMomentum = aP
2162 G4KineticTrack* aTrack = new G4Kineti
2163
2164
2165
2166 FirstString = new G4ExcitedString( aT
2167
2168 #ifdef debugBuildString
2169 G4cout << " Strings are build for inv
2170 #endif
2171
2172 if ( aProjectile->GetStatus() == 2 )
2173 } else {
2174
2175 #ifdef debugBuildString
2176 G4cout << "Case5 " << G4endl;
2177 #endif
2178
2179 //TheInvolvedNucleonsOfProjectile[ ah
2180 //G4cout << TheInvolvedNucleonsOfProj
2181
2182 #ifdef debugBuildString
2183 G4cout << " No string" << G4endl;
2184 #endif
2185
2186 }
2187
2188 if ( FirstString != 0 ) strings->push_
2189 if ( SecondString != 0 ) strings->push_
2190 }
2191 }
2192
2193 #ifdef debugBuildString
2194 G4cout << "Building of target-like strings"
2195 #endif
2196
2197 G4bool isProjectile = false;
2198 for ( G4int ahadron = 0; ahadron < NumberOf
2199 G4VSplitableHadron* aNucleon = TheInvolve
2200
2201 #ifdef debugBuildString
2202 G4cout << "Nucleon #, status, intCount "
2203 << aNucleon->GetStatus() << " " <<
2204 #endif
2205
2206 FirstString = 0 ; SecondString = 0;
2207
2208 if ( aNucleon->GetStatus() == 0 ) { // A
2209 theExcitation->CreateStrings( aNucleon,
2210 FirstStri
2211 NumberOfTargetSpectatorNucleons--;
2212
2213 #ifdef debugBuildString
2214 G4cout << " 1 case A string is build" <
2215 #endif
2216
2217 } else if ( aNucleon->GetStatus() == 1 &
2218 // A nucleon took part in diffractive i
2219 theExcitation->CreateStrings( aNucleon,
2220 FirstStri
2221
2222 #ifdef debugBuildString
2223 G4cout << " 2 case A string is build, n
2224 #endif
2225
2226 NumberOfTargetSpectatorNucleons--;
2227
2228 } else if ( aNucleon->GetStatus() == 1 &
2229 HighEnergyInter ) {
2230 // A nucleon was considered as a partic
2231 // its interactions were skipped. It wi
2232 // at high energies.
2233
2234 G4LorentzVector ParticleMomentum = aNuc
2235 G4KineticTrack* aTrack = new G4KineticT
2236
2237
2238
2239
2240 FirstString = new G4ExcitedString( aTra
2241
2242 #ifdef debugBuildString
2243 G4cout << "3 case A string is build" <<
2244 #endif
2245
2246 } else if ( aNucleon->GetStatus() == 1 &
2247 ! HighEnergyInter ) {
2248 // A nucleon was considered as a partic
2249 // its interactions were skipped. It wi
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
2256 #endif
2257
2258 } else if ( aNucleon->GetStatus() == 2 |
2259 aNucleon->GetStatus() == 3 )
2260 G4LorentzVector ParticleMomentum = aNuc
2261 G4KineticTrack* aTrack = new G4KineticT
2262
2263
2264
2265 FirstString = new G4ExcitedString( aTra
2266
2267 #ifdef debugBuildString
2268 G4cout << "5 case A string is build" <<
2269 #endif
2270
2271 if ( aNucleon->GetStatus() == 2 ) Numbe
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_ba
2282 if ( SecondString != 0 ) strings->push_ba
2283
2284 }
2285
2286 #ifdef debugBuildString
2287 G4cout << G4endl << "theAdditionalString.si
2288 << G4endl << G4endl;
2289 #endif
2290
2291 isProjectile = true;
2292 if ( theAdditionalString.size() != 0 ) {
2293 for ( unsigned int ahadron = 0; ahadron
2294 //if ( theAdditionalString[ ahadron ]->
2295 FirstString = 0; SecondString = 0;
2296 theExcitation->CreateStrings( theAdditi
2297 FirstStri
2298 if ( FirstString != 0 ) strings->push_
2299 if ( SecondString != 0 ) strings->push_
2300 }
2301 }
2302
2303 //for ( unsigned int ahadron = 0; ahadron <
2304 // G4cout << ahadron << " " << strings->op
2305 // << " " << strings->operator[]( a
2306 //}
2307 //G4cout << "------------------------" << G
2308
2309 return;
2310 }
2311
2312
2313 //===========================================
2314
2315 void G4FTFModel::GetResiduals() {
2316 // This method is needed for the correct ap
2317
2318 #ifdef debugFTFmodel
2319 G4cout << "GetResiduals(): HighEnergyInter?
2320 << HighEnergyInter << " " << GetProj
2321 #endif
2322
2323 if ( HighEnergyInter ) {
2324
2325 #ifdef debugFTFmodel
2326 G4cout << "NumberOfInvolvedNucleonsOfTarg
2327 #endif
2328
2329 G4double DeltaExcitationE = TargetResidua
2330 G4double( Num
2331 G4LorentzVector DeltaPResidualNucleus = T
2332 G
2333
2334 for ( G4int i = 0; i < NumberOfInvolvedNu
2335 G4Nucleon* aNucleon = TheInvolvedNucleo
2336
2337 #ifdef debugFTFmodel
2338 G4VSplitableHadron* targetSplitable = a
2339 G4cout << i << " Hit? " << aNucleon->Ar
2340 if ( targetSplitable ) G4cout << i << "
2341 #endif
2342
2343 G4LorentzVector tmp = -DeltaPResidualNu
2344 aNucleon->SetMomentum( tmp );
2345 aNucleon->SetBindingEnergy( DeltaExcita
2346 }
2347
2348 if ( TargetResidualMassNumber != 0 ) {
2349 G4ThreeVector bstToCM = TargetResidual4
2350
2351 G4V3DNucleus* theTargetNucleus = GetTar
2352 G4LorentzVector residualMomentum( 0.0,
2353 G4Nucleon* aNucleon = 0;
2354 theTargetNucleus->StartLoop();
2355 while ( ( aNucleon = theTargetNucleus->
2356 if ( ! aNucleon->AreYouHit() ) {
2357 G4LorentzVector tmp = aNucleon->Get
2358 aNucleon->SetMomentum( tmp );
2359 residualMomentum += tmp;
2360 }
2361 }
2362
2363 residualMomentum /= TargetResidualMassN
2364
2365 G4double Mass = TargetResidual4Momentum
2366 G4double SumMasses = 0.0;
2367
2368 aNucleon = 0;
2369 theTargetNucleus->StartLoop();
2370 while ( ( aNucleon = theTargetNucleus->
2371 if ( ! aNucleon->AreYouHit() ) {
2372 G4LorentzVector tmp = aNucleon->Get
2373 G4double E = std::sqrt( tmp.vect().
2374 sqr( aNucle
2375 tmp.setE( E ); aNucleon->SetMoment
2376 SumMasses += E;
2377 }
2378 }
2379
2380 G4double Chigh = Mass / SumMasses; G4do
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
2389 if ( ! aNucleon->AreYouHit() ) {
2390 G4LorentzVector tmp = aNucleon->G
2391 G4double E = std::sqrt( tmp.vect(
2392 sqr( aNuc
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 < maxNumberOfLo
2402 if ( loopCounter >= maxNumberOfLoops )
2403 #ifdef debugFTFmodel
2404 G4cout << "BAD situation: forced exit
2405 << "\t return immediately from
2406 #endif
2407 return;
2408 }
2409
2410 aNucleon = 0;
2411 theTargetNucleus->StartLoop();
2412 while ( ( aNucleon = theTargetNucleus->
2413 if ( !aNucleon->AreYouHit() ) {
2414 G4LorentzVector tmp = aNucleon->Get
2415 G4double E = std::sqrt( tmp.vect().
2416 sqr( aNucle
2417 tmp.setE( E ); tmp.boost( -bstToCM
2418 aNucleon->SetMomentum( tmp );
2419 }
2420 }
2421 }
2422
2423 if ( ! GetProjectileNucleus() ) return;
2424
2425 #ifdef debugFTFmodel
2426 G4cout << "NumberOfInvolvedNucleonsOfProj
2427 << G4endl << "ProjectileResidualEx
2428 << ProjectileResidualExcitationEne
2429 #endif
2430
2431 DeltaExcitationE = ProjectileResidualExci
2432 G4double( NumberOfInvo
2433 DeltaPResidualNucleus = ProjectileResidua
2434 G4double( NumberO
2435
2436 for ( G4int i = 0; i < NumberOfInvolvedNu
2437 G4Nucleon* aNucleon = TheInvolvedNucleo
2438
2439 #ifdef debugFTFmodel
2440 G4VSplitableHadron* projSplitable = aNu
2441 G4cout << i << " Hit? " << aNucleon->Ar
2442 if ( projSplitable ) G4cout << i << "St
2443 #endif
2444
2445 G4LorentzVector tmp = -DeltaPResidualNu
2446 aNucleon->SetMomentum( tmp );
2447 aNucleon->SetBindingEnergy( DeltaExcita
2448 }
2449
2450 if ( ProjectileResidualMassNumber != 0 )
2451 G4ThreeVector bstToCM = ProjectileResid
2452
2453 G4V3DNucleus* theProjectileNucleus = Ge
2454 G4LorentzVector residualMomentum( 0.0,
2455 G4Nucleon* aNucleon = 0;
2456 theProjectileNucleus->StartLoop();
2457 while ( ( aNucleon = theProjectileNucle
2458 if ( ! aNucleon->AreYouHit() ) {
2459 G4LorentzVector tmp = aNucleon->Get
2460 aNucleon->SetMomentum( tmp );
2461 residualMomentum += tmp;
2462 }
2463 }
2464
2465 residualMomentum /= ProjectileResidualM
2466
2467 G4double Mass = ProjectileResidual4Mome
2468 G4double SumMasses= 0.0;
2469
2470 aNucleon = 0;
2471 theProjectileNucleus->StartLoop();
2472 while ( ( aNucleon = theProjectileNucle
2473 if ( ! aNucleon->AreYouHit() ) {
2474 G4LorentzVector tmp = aNucleon->Get
2475 G4double E=std::sqrt( tmp.vect().ma
2476 sqr(aNucleon-
2477 tmp.setE( E ); aNucleon->SetMoment
2478 SumMasses += E;
2479 }
2480 }
2481
2482 G4double Chigh = Mass / SumMasses; G4do
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 = theProjectileNuc
2492 if ( ! aNucleon->AreYouHit() ) {
2493 G4LorentzVector tmp = aNucleon->G
2494 G4double E = std::sqrt( tmp.vect(
2495 sqr( aNuc
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 < maxNumberOfLo
2505 if ( loopCounter >= maxNumberOfLoops )
2506 #ifdef debugFTFmodel
2507 G4cout << "BAD situation: forced exit
2508 << "\t return immediately from
2509 #endif
2510 return;
2511 }
2512
2513 aNucleon = 0;
2514 theProjectileNucleus->StartLoop();
2515 while ( ( aNucleon = theProjectileNucle
2516 if ( ! aNucleon->AreYouHit() ) {
2517 G4LorentzVector tmp = aNucleon->Get
2518 G4double E = std::sqrt( tmp.vect().
2519 sqr( aNucle
2520 tmp.setE( E ); tmp.boost( -bstToCM
2521 aNucleon->SetMomentum( tmp );
2522 }
2523 }
2524 } // End of if ( ProjectileResidualMass
2525
2526 #ifdef debugFTFmodel
2527 G4cout << "End projectile" << G4endl;
2528 #endif
2529
2530 } else { // Related to the condition: if (
2531
2532 #ifdef debugFTFmodel
2533 G4cout << "Low energy interaction: Target
2534 << "Tr ResidualMassNumber Tr Resid
2535 << TargetResidualMassNumber << " "
2536 << TargetResidualExcitationEnergy
2537 #endif
2538
2539 G4int NumberOfTargetParticipant( 0 );
2540 for ( G4int i = 0; i < NumberOfInvolvedNu
2541 G4Nucleon* aNucleon = TheInvolvedNucleo
2542 G4VSplitableHadron* targetSplitable = a
2543 if ( targetSplitable->GetSoftCollisionC
2544 }
2545
2546 G4double DeltaExcitationE( 0.0 );
2547 G4LorentzVector DeltaPResidualNucleus = G
2548
2549 if ( NumberOfTargetParticipant != 0 ) {
2550 DeltaExcitationE = TargetResidualExcita
2551 DeltaPResidualNucleus = TargetResidual4
2552 }
2553
2554 for ( G4int i = 0; i < NumberOfInvolvedNu
2555 G4Nucleon* aNucleon = TheInvolvedNucleo
2556 G4VSplitableHadron* targetSplitable = a
2557 if ( targetSplitable->GetSoftCollisionC
2558 G4LorentzVector tmp = -DeltaPResidual
2559 aNucleon->SetMomentum( tmp );
2560 aNucleon->SetBindingEnergy( DeltaExci
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 " <<
2571 << "TargetResidual4Momentum " <<
2572 #endif
2573
2574 if ( ! GetProjectileNucleus() ) return;
2575
2576 #ifdef debugFTFmodel
2577 G4cout << "Low energy interaction: Projec
2578 << "Pr ResidualMassNumber Pr Resid
2579 << ProjectileResidualMassNumber <<
2580 << ProjectileResidualExcitationEne
2581 #endif
2582
2583 G4int NumberOfProjectileParticipant( 0 );
2584 for ( G4int i = 0; i < NumberOfInvolvedNu
2585 G4Nucleon* aNucleon = TheInvolvedNucleo
2586 G4VSplitableHadron* projectileSplitable
2587 if ( projectileSplitable->GetSoftCollis
2588 }
2589
2590 #ifdef debugFTFmodel
2591 G4cout << "NumberOfProjectileParticipant"
2592 #endif
2593
2594 DeltaExcitationE = 0.0;
2595 DeltaPResidualNucleus = G4LorentzVector(
2596
2597 if ( NumberOfProjectileParticipant != 0 )
2598 DeltaExcitationE = ProjectileResidualEx
2599 DeltaPResidualNucleus = ProjectileResid
2600 }
2601 //G4cout << "DeltaExcitationE DeltaPResid
2602 // << " " << DeltaPResidualNucleus
2603 for ( G4int i = 0; i < NumberOfInvolvedNu
2604 G4Nucleon* aNucleon = TheInvolvedNucleo
2605 G4VSplitableHadron* projectileSplitable
2606 if ( projectileSplitable->GetSoftCollis
2607 G4LorentzVector tmp = -DeltaPResidual
2608 aNucleon->SetMomentum( tmp );
2609 aNucleon->SetBindingEnergy( DeltaExci
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
2620 << "ProjectileResidual4Momentum "
2621 #endif
2622
2623 } // End of the condition: if ( HighEnergy
2624
2625 #ifdef debugFTFmodel
2626 G4cout << "End GetResiduals ---------------
2627 #endif
2628
2629 }
2630
2631
2632 //===========================================
2633
2634 G4ThreeVector G4FTFModel::GaussianPt( G4doubl
2635
2636 G4double Pt2( 0.0 ), Pt( 0.0 );
2637
2638 if (AveragePt2 > 0.0) {
2639 const G4double ymax = maxPtSquare/Average
2640 if ( ymax < 200. ) {
2641 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unif
2642 } else {
2643 Pt2 = -AveragePt2 * G4Log( 1.0 - G4Unif
2644 }
2645 Pt = std::sqrt( Pt2 );
2646 }
2647
2648 G4double phi = G4UniformRand() * twopi;
2649
2650 return G4ThreeVector( Pt*std::cos(phi), Pt*
2651 }
2652
2653 //===========================================
2654
2655 G4bool G4FTFModel::
2656 ComputeNucleusProperties( G4V3DNucleus* nucle
2657 G4LorentzVector& nu
2658 G4LorentzVector& re
2659 G4double& sumMasses
2660 G4double& residualE
2661 G4double& residualM
2662 G4int& residualMass
2663 G4int& residualChar
2664
2665 // This method, which is called only by Put
2666 // - either the target nucleus (which is n
2667 // of hadronic interaction (hadron-nucle
2668 // - or the projectile nucleus or antinucl
2669 // or antinucleus-nucleus interaction.
2670 // This method assumes that the all the par
2671 // the action of this method consists in mo
2672 // first one. The return value is "false" o
2673 // is null.
2674
2675 if ( ! nucleus ) return false;
2676
2677 G4double ExcitationEnergyPerWoundedNucleon
2678 theParameters->GetExcitationEnergyPerWoun
2679
2680 // Loop over the nucleons of the nucleus.
2681 // The nucleons that have been involved in
2682 // Reggeon Cascading) will be candidate to
2683 // All the remaining nucleons will be the n
2684 // The variable sumMasses is the amount of
2685 // 1. transverse mass of each involved
2686 // 2. 20.0*MeV separation energy for ea
2687 // 3. transverse mass of the residual n
2688 // In this first evaluation of sumMasses, t
2689 // (residualExcitationEnergy, estimated by
2690 // nucleon) is not taken into account.
2691 G4int residualNumberOfLambdas = 0; // Proj
2692 G4Nucleon* aNucleon = 0;
2693 nucleus->StartLoop();
2694 while ( ( aNucleon = nucleus->GetNextNucleo
2695 nucleusMomentum += aNucleon->Get4Momentum
2696 if ( aNucleon->AreYouHit() ) { // Involv
2697 // Consider in sumMasses the nominal, i
2698 // (not the current masses, which could
2699 sumMasses += std::sqrt( sqr( aNucleon->
2700 + aNucleon->Ge
2701 sumMasses += 20.0*MeV; // Separation e
2702
2703 //residualExcitationEnergy += Excitatio
2704 residualExcitationEnergy += -Excitation
2705
2706 residualMassNumber--;
2707 // The absolute value below is needed o
2708 residualCharge -= std::abs( G4int( aNuc
2709 } else { // Spectator nucleons
2710 residualMomentum += aNucleon->Get4Momen
2711 if ( aNucleon->GetDefinition() == G4Lam
2712 aNucleon->GetDefinition() == G4AntiLambd
2713 ++residualNumberOfLambdas;
2714 }
2715 }
2716 }
2717 #ifdef debugPutOnMassShell
2718 G4cout << "ExcitationEnergyPerWoundedNucleo
2719 << "\t Residual Charge, MassNumber (
2720 << residualMassNumber << " (" << residualN
2721 << G4endl << "\t Initial Momentum "
2722 << G4endl << "\t Residual Momentum
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()
2733 } else if ( residualNumberOfLambdas ==
2734 residualMass = G4Lambda::Definition()
2735 } else {
2736 residualMass = G4Neutron::Definition(
2737 }
2738 residualExcitationEnergy = 0.0;
2739 } else {
2740 if ( residualNumberOfLambdas > 0 ) {
2741 if ( residualMassNumber == 2 ) {
2742 residualMass = G4Lambda::Definition()->Ge
2743 if ( std::abs( residualCharge ) ==
2744 residualMass += G4Proton::Definit
2745 } else if ( residualNumberOfLambdas == 1
2746 residualMass += G4Neutron::Definition()
2747 } else {
2748 residualMass += G4Lambda::Definition()-
2749 }
2750 } else {
2751 residualMass = G4HyperNucleiProperties::G
2752 residualNumberOfLambdas )
2753 }
2754 } else {
2755 residualMass = G4ParticleTable::GetPa
2756 GetIonMass( std::abs( residu
2757 }
2758 }
2759 residualMass += residualExcitationEnergy;
2760 }
2761 sumMasses += std::sqrt( sqr( residualMass )
2762 return true;
2763 }
2764
2765
2766 //===========================================
2767
2768 G4bool G4FTFModel::
2769 GenerateDeltaIsobar( const G4double sqrtS,
2770 const G4int numberOfInvo
2771 G4Nucleon* involvedNucle
2772 G4double& sumMasses ) {
2773
2774 // This method, which is called only by Put
2775 // re-interpret some of the involved nucleo
2776 // - either by replacing a proton (2212) wi
2777 // - or by replacing a neutron (2112) with
2778 // The on-shell mass of these delta-isobars
2779 // the corresponding nucleon on-shell mass.
2780 // the max number of deltas compatible with
2781 // The delta-isobars are considered with th
2782 // corresponding nucleons.
2783 // This method assumes that all the paramet
2784 // the action of this method consists in mo
2785 // sumMasses. The return value is "false" o
2786 // have unphysical values.
2787
2788 if ( sqrtS < 0.0 || numberOfInvolvedNucle
2789
2790 const G4double probDeltaIsobar = 0.05;
2791
2792 G4int maxNumberOfDeltas = G4int( (sqrtS - s
2793 G4int numberOfDeltas = 0;
2794
2795 for ( G4int i = 0; i < numberOfInvolvedNucl
2796
2797 if ( G4UniformRand() < probDeltaIsobar &
2798 numberOfDeltas++;
2799 if ( ! involvedNucleons[i] ) continue;
2800 // Skip any eventual lambda (that can b
2801 if ( involvedNucleons[i]->GetDefinition
2802 involvedNucleons[i]->GetDefinition() ==
2803 G4VSplitableHadron* splitableHadron = i
2804 G4double massNuc = std::sqrt( sqr( spli
2805 + splitab
2806 // The absolute value below is needed i
2807 G4int pdgCode = std::abs( splitableHadr
2808 const G4ParticleDefinition* old_def = s
2809 G4int newPdgCode = pdgCode/10; newPdgCo
2810 if ( splitableHadron->GetDefinition()->
2811 const G4ParticleDefinition* ptr =
2812 G4ParticleTable::GetParticleTable()->
2813 splitableHadron->SetDefinition( ptr );
2814 G4double massDelta = std::sqrt( sqr( sp
2815 + split
2816 //G4cout << i << " " << sqrtS/GeV << "
2817 // << " " << massNuc << G4endl;
2818 if ( sqrtS < sumMasses + massDelta - ma
2819 splitableHadron->SetDefinition( old_d
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 averagePt
2835 const G4double max
2836 G4double dCor,
2837 G4V3DNucleus* nucl
2838 const G4LorentzVec
2839 const G4double res
2840 const G4int residu
2841 const G4int number
2842 G4Nucleon* involve
2843 G4double& mass2 )
2844
2845 // This method, which is called only by Put
2846 // - either the target nucleons: this for
2847 // (hadron-nucleus, nucleus-nucleus, ant
2848 // - or the projectile nucleons or antinuc
2849 // nucleus-nucleus or antinucleus-nucleu
2850 // This method assumes that all the paramet
2851 // the action of this method consists in ch
2852 // whose pointers are in the vector involve
2853 // variable mass2.
2854 #ifdef debugPutOnMassShell
2855 G4cout << "G4FTFModel::SamplingNucleonKinem
2856 G4cout << " averagePt2= " << averagePt2 <<
2857 << " dCor= " << dCor << " resMass(GeV)= "
2858 << " resMassN= " << residualMassNumber
2859 << " nNuc= " << numberOfInvolvedNucl
2860 << " lv= " << pResidual << G4endl;
2861 #endif
2862
2863 if ( ! nucleus || numberOfInvolvedNucleon
2864
2865 if ( residualMassNumber == 0 && numberOfI
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)numberOfInv
2874
2875 // to avoid problems due to precision lost
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 < numberOfInvolved
2887 G4Nucleon* aNucleon = involvedNucleons[i];
2888 if ( ! aNucleon ) continue;
2889 G4ThreeVector tmpPt = GaussianPt( averagePt
2890 ptSum += tmpPt;
2891 G4LorentzVector tmp( tmpPt.x(), tmpPt.y(),
2892 aNucleon->SetMomentum( tmp );
2893 }
2894 }
2895
2896 G4double deltaPx = ( ptSum.x() - pResidua
2897 G4double deltaPy = ( ptSum.y() - pResidua
2898
2899 SumMasses = residualMass;
2900 for ( G4int i = 0; i < numberOfInvolvedNu
2901 G4Nucleon* aNucleon = involvedNucleons[
2902 if ( ! aNucleon ) continue;
2903 G4double px = aNucleon->Get4Momentum().
2904 G4double py = aNucleon->Get4Momentum().
2905 G4double MtN = std::sqrt( sqr( aNucleon
2906 + sqr( px ) +
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 < numberOfInvolvedNu
2916 G4Nucleon* aNucleon = involvedNucleons[
2917 if ( ! aNucleon ) continue;
2918
2919 G4double x = 0.0;
2920 if( 0.0 != dCor ) {
2921 G4ThreeVector tmpX = GaussianPt( dCor*dCor,
2922 x = tmpX.x();
2923 }
2924 x += aNucleon->Get4Momentum().e()/SumMa
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
2932
2933 G4LorentzVector tmp( aNucleon->Get4Mome
2934 aNucleon->Get4Mome
2935 x, aNucleon->Get4M
2936 aNucleon->SetMomentum( tmp );
2937 }
2938
2939 if ( xSum < -eps || xSum > 1.0 + eps ) su
2940 if ( ! success ) continue;
2941
2942 G4double delta = ( residualMassNumber ==
2943
2944 xSum = 1.0;
2945 mass2 = 0.0;
2946 for ( G4int i = 0; i < numberOfInvolvedNu
2947 G4Nucleon* aNucleon = involvedNucleons[
2948 if ( ! aNucleon ) continue;
2949 G4double x = aNucleon->Get4Momentum().p
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 || xS
2959 success = false;
2960 break;
2961 }
2962 }
2963 x = std::min( 1.0, std::max(x, eps) );
2964
2965 mass2 += sqr( aNucleon->Get4Momentum().
2966
2967 G4LorentzVector tmp( aNucleon->Get4Mome
2968 x, aNucleon->Get4M
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 += (
2975
2976 #ifdef debugPutOnMassShell
2977 G4cout << "success: " << success << " Mt(
2978 << std::sqrt( mass2 )/GeV << G4endl;
2979 #endif
2980
2981 } while ( ( ! success ) &&
2982 ++loopCounter < maxNumberOfLoops
2983 return ( loopCounter < maxNumberOfLoops );
2984 }
2985
2986
2987 //===========================================
2988
2989 G4bool G4FTFModel::
2990 CheckKinematics( const G4double sValue,
2991 const G4double sqrtS,
2992 const G4double projectileMas
2993 const G4double targetMass2,
2994 const G4double nucleusY,
2995 const G4bool isProjectileNuc
2996 const G4int numberOfInvolved
2997 G4Nucleon* involvedNucleons[
2998 G4double& targetWminus,
2999 G4double& projectileWplus,
3000 G4bool& success ) {
3001
3002 // This method, which is called only by Put
3003 // kinematics is acceptable or not.
3004 // This method assumes that all the paramet
3005 // notice that the input boolean parameter
3006 // only in the case of nucleus or antinucle
3007 // The action of this method consists in co
3008 // and setting the parameter success to fal
3009 // be rejeted.
3010
3011 G4double decayMomentum2 = sqr( sValue ) + s
3012 - 2.0*( sValue*(
3013 + project
3014 targetWminus = ( sValue - projectileMass2 +
3015 / 2.0 / sqrtS;
3016 projectileWplus = sqrtS - targetMass2/targe
3017 G4double projectilePz = projectileWplus/2.0
3018 G4double projectileE = projectileWplus/2.0
3019 G4double projectileY = 0.5 * G4Log( (proje
3020 (proje
3021 G4double targetPz = -targetWminus/2.0 + tar
3022 G4double targetE = targetWminus/2.0 + tar
3023 G4double targetY = 0.5 * G4Log( (targetE +
3024
3025 #ifdef debugPutOnMassShell
3026 G4cout << "decayMomentum2 " << decayMomentu
3027 << "\t targetWminus projectileWplus
3028 << "\t projectileY targetY " << proj
3029 if ( isProjectileNucleus ) {
3030 G4cout << "Order# of Wounded nucleon i, n
3031 } else {
3032 G4cout << "Order# of Wounded nucleon i, n
3033 }
3034 G4cout << G4endl;
3035 #endif
3036
3037 for ( G4int i = 0; i < numberOfInvolvedNucl
3038 G4Nucleon* aNucleon = involvedNucleons[i]
3039 if ( ! aNucleon ) continue;
3040 G4LorentzVector tmp = aNucleon->Get4Momen
3041 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.
3042 sqr( aNucleon->GetSplitabl
3043 G4double x = tmp.z();
3044 G4double pz = -targetWminus*x/2.0 + mt2/(
3045 G4double e = targetWminus*x/2.0 + mt2/(
3046 if ( isProjectileNucleus ) {
3047 pz = projectileWplus*x/2.0 - mt2/(2.0*p
3048 e = projectileWplus*x/2.0 + mt2/(2.0*p
3049 }
3050 G4double nucleonY = 0.5 * G4Log( (e + pz)
3051
3052 #ifdef debugPutOnMassShell
3053 if( isProjectileNucleus ) {
3054 G4cout << " " << i << " " << nucleonY <
3055 } else {
3056 G4cout << " " << i << " " << nucleonY <
3057 }
3058 G4cout << G4endl;
3059 #endif
3060
3061 if ( std::abs( nucleonY - nucleusY ) > 2
3062 ( isProjectileNucleus && targetY >
3063 ( ! isProjectileNucleus && project
3064 success = false;
3065 break;
3066 }
3067 }
3068 return true;
3069 }
3070
3071
3072 //===========================================
3073
3074 G4bool G4FTFModel::
3075 FinalizeKinematics( const G4double w,
3076 const G4bool isProjectile
3077 const G4LorentzRotation&
3078 const G4double residualMa
3079 const G4int residualMassN
3080 const G4int numberOfInvol
3081 G4Nucleon* involvedNucleo
3082 G4LorentzVector& residual4Momen
3083
3084 // This method, which is called only by Put
3085 // this method is called when we are sure t
3086 // acceptable.
3087 // This method assumes that all the paramet
3088 // notice that the input boolean parameter
3089 // only in the case of nucleus or antinucle
3090 // because the sign of pz (in the center-of
3091 // with respect to the case of a normal had
3092 // The action of this method consists in mo
3093 // (in the lab frame) and computing the res
3094 // frame).
3095
3096 G4ThreeVector residual3Momentum( 0.0, 0.0,
3097
3098 for ( G4int i = 0; i < numberOfInvolvedNucl
3099 G4Nucleon* aNucleon = involvedNucleons[i]
3100 if ( ! aNucleon ) continue;
3101 G4LorentzVector tmp = aNucleon->Get4Momen
3102 residual3Momentum -= tmp.vect();
3103 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.
3104 sqr( aNucleon->GetSplitabl
3105 G4double x = tmp.z();
3106 G4double pz = -w * x / 2.0 + mt2 / ( 2.
3107 G4double e = w * x / 2.0 + mt2 / ( 2.
3108 // Reverse the sign of pz in the case of
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 = aNu
3115 splitableHadron->Set4Momentum( tmp );
3116 }
3117
3118 G4double residualMt2 = sqr( residualMass )
3119 + sqr( residual3Moment
3120
3121 #ifdef debugPutOnMassShell
3122 if ( isProjectileNucleus ) {
3123 G4cout << "Wminus Proj and residual3Momen
3124 } else {
3125 G4cout << "Wplus Targ and residual3Momen
3126 }
3127 #endif
3128
3129 G4double residualPz = 0.0;
3130 G4double residualE = 0.0;
3131 if ( residualMassNumber != 0 ) {
3132 residualPz = -w * residual3Momentum.z() /
3133 residualMt2 / ( 2.0 * w * r
3134 residualE = w * residual3Momentum.z() /
3135 residualMt2 / ( 2.0 * w * r
3136 // Reverse the sign of residualPz in the
3137 if ( isProjectileNucleus ) residualPz *=
3138 }
3139
3140 residual4Momentum.setPx( residual3Momentum.
3141 residual4Momentum.setPy( residual3Momentum.
3142 residual4Momentum.setPz( residualPz );
3143 residual4Momentum.setE( residualE );
3144
3145 return true;
3146 }
3147
3148
3149 //===========================================
3150
3151 void G4FTFModel::ModelDescription( std::ostre
3152 desc << " FTF (Fritiof) Mod
3153 << "The FTF model is based on the well
3154 << "model (B. Andersson et al., Nucl.
3155 << "(1987)). Its first program impleme
3156 << "by B. Nilsson-Almquist and E. Sten
3157 << "Comm. 43, 387 (1987)). The Fritiof
3158 << "that all hadron-hadron interaction
3159 << "reactions, h_1+h_2->h_1'+h_2' wher
3160 << "are excited states of the hadrons
3161 << "mass spectra. The excited hadrons
3162 << "QCD-strings, and the corresponding
3163 << "fragmentation model is applied for
3164 << "their decays.
3165 << " The Fritiof model assumes that
3166 << "a hadron-nucleus interaction a str
3167 << "from the projectile can interact w
3168 << "nuclear nucleons and becomes into
3169 << "states. The probability of multipl
3170 << "calculated in the Glauber approxim
3171 << "of secondary particles was neglect
3172 << "to these, the original Fritiof mod
3173 << "cribe a nuclear destruction and sl
3174 << " In order to overcome the diffic
3175 << "the model by the reggeon theory in
3176 << "nuclear desctruction (Kh. Abdel-Wa
3177 << "nsky, Phys. Atom. Nucl. 60, 828 (1
3178 << "(1997)). Momenta of the nucleons e
3179 << "leus in the reggeon cascading are
3180 << "to a Fermi motion algorithm presen
3181 << "Collaboration (M.I. Adamovich et a
3182 << "A358, 337 (1997)).
3183 << " New features were also added to
3184 << "implemented in Geant4: a simulatio
3185 << "ron-nucleon scatterings, a simulat
3186 << "reactions like NN>NN* in hadron-nu
3187 << "a separate simulation of single di
3188 << " diffractive events. These allowed
3189 << "model parameter tuning a wide set
3190 << "data.
3191 }
3192
3193