| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4FTFModel.cc,v 1.7 2007/04/24 10:32:59 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
27 // 29 //
28 30
29 // ------------------------------------------- 31 // ------------------------------------------------------------
30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file
31 // 33 //
32 // ---------------- G4FTFModel ---------- 34 // ---------------- G4FTFModel ----------------
33 // by Gunter Folger, May 1998. 35 // by Gunter Folger, May 1998.
34 // class implementing the excitation in 36 // class implementing the excitation in the FTF Parton String Model
35 // <<
36 // Vladimir Uzhinsky, November <<
37 // simulation of nucleus-nucleus interac <<
38 // ------------------------------------------- 37 // ------------------------------------------------------------
39 38
40 #include <utility> <<
41 <<
42 #include "G4FTFModel.hh" 39 #include "G4FTFModel.hh"
43 #include "G4ios.hh" <<
44 #include "G4PhysicalConstants.hh" <<
45 #include "G4SystemOfUnits.hh" <<
46 #include "G4FTFParameters.hh" <<
47 #include "G4FTFParticipants.hh" 40 #include "G4FTFParticipants.hh"
48 #include "G4DiffractiveSplitableHadron.hh" <<
49 #include "G4InteractionContent.hh" 41 #include "G4InteractionContent.hh"
50 #include "G4LorentzRotation.hh" 42 #include "G4LorentzRotation.hh"
51 #include "G4ParticleDefinition.hh" 43 #include "G4ParticleDefinition.hh"
52 #include "G4ParticleTable.hh" << 44 #include "G4ios.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 45
840 return true; << 46 // Class G4FTFModel
841 47
>> 48 G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()) // Uzhi
>> 49 {
>> 50 G4VPartonStringModel::SetThisPointer(this);
842 } 51 }
843 52
844 << 53 G4FTFModel::G4FTFModel(G4double a, G4double b, G4double c):theExcitation(new G4DiffractiveExcitation())
845 //============================================ << 54 {
846 << 55 G4VPartonStringModel::SetThisPointer(this);
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 } 56 }
1029 57
1030 << 58 G4FTFModel::G4FTFModel(G4DiffractiveExcitation * anExcitation)
1031 //=========================================== << 59 :
1032 << 60 theExcitation(anExcitation)
1033 G4bool G4FTFModel::AdjustNucleons( G4VSplitab << 61 {
1034 G4Nucleon* << 62 G4VPartonStringModel::SetThisPointer(this);
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 } 63 }
1146 64
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 65
1520 66
1521 //------------------------------------------- << 67 G4FTFModel::~G4FTFModel()
>> 68 {}
1522 69
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 70
1853 return true; << 71 const G4FTFModel & G4FTFModel::operator=(const G4FTFModel &)
>> 72 {
>> 73 throw G4HadronicException(__FILE__, __LINE__, "G4FTFModel::operator= is not meant to be accessed ");
>> 74 return *this;
1854 } 75 }
1855 76
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 77
>> 78 int G4FTFModel::operator==(const G4FTFModel &right) const
>> 79 {
>> 80 return this==&right;
1995 } 81 }
1996 82
1997 << 83 int G4FTFModel::operator!=(const G4FTFModel &right) const
1998 //=========================================== << 84 {
1999 << 85 return this!=&right;
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 } 86 }
2311 87
2312 << 88 void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
2313 //=========================================== << 89 {
2314 << 90 theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); // Uzhi N-mass number Z-charge
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 91
2590 #ifdef debugFTFmodel << 92 theProjectile = aProjectile;
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 } 93 }
2630 94
>> 95 G4ExcitedStringVector * G4FTFModel::GetStrings()
>> 96 {
>> 97 theParticipants.BuildInteractions(theProjectile);
2631 98
2632 //=========================================== << 99 if (! ExciteParticipants()) return NULL;;
2633 <<
2634 G4ThreeVector G4FTFModel::GaussianPt( G4doubl <<
2635 <<
2636 G4double Pt2( 0.0 ), Pt( 0.0 ); <<
2637 100
2638 if (AveragePt2 > 0.0) { << 101 G4ExcitedStringVector * theStrings = BuildStrings();
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 102
2648 G4double phi = G4UniformRand() * twopi; << 103 return theStrings;
2649 <<
2650 return G4ThreeVector( Pt*std::cos(phi), Pt* <<
2651 } 104 }
2652 105
2653 //=========================================== << 106 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} };
2654 << 107
2655 G4bool G4FTFModel:: << 108 G4ExcitedStringVector * G4FTFModel::BuildStrings()
2656 ComputeNucleusProperties( G4V3DNucleus* nucle << 109 {
2657 G4LorentzVector& nu << 110
2658 G4LorentzVector& re << 111 // Loop over all collisions; find all primaries, and all target ( targets may
2659 G4double& sumMasses << 112 // be duplicate in the List ( to unique G4VSplitableHadrons)
2660 G4double& residualE << 113
2661 G4double& residualM << 114 G4ExcitedStringVector * strings;
2662 G4int& residualMass << 115 strings = new G4ExcitedStringVector();
2663 G4int& residualChar << 116
2664 << 117 std::vector<G4VSplitableHadron *> primaries;
2665 // This method, which is called only by Put << 118 std::vector<G4VSplitableHadron *> targets;
2666 // - either the target nucleus (which is n << 119
2667 // of hadronic interaction (hadron-nucle << 120 theParticipants.StartLoop(); // restart a loop
2668 // - or the projectile nucleus or antinucl << 121 while ( theParticipants.Next() )
2669 // or antinucleus-nucleus interaction. << 122 {
2670 // This method assumes that the all the par << 123 const G4InteractionContent & interaction=theParticipants.GetInteraction();
2671 // the action of this method consists in mo << 124 // do not allow for duplicates ...
2672 // first one. The return value is "false" o << 125 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) )
2673 // is null. << 126 primaries.push_back(interaction.GetProjectile());
2674 << 127
2675 if ( ! nucleus ) return false; << 128 if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) )
2676 << 129 targets.push_back(interaction.GetTarget());
2677 G4double ExcitationEnergyPerWoundedNucleon << 130 }
2678 theParameters->GetExcitationEnergyPerWoun << 131
2679 << 132
2680 // Loop over the nucleons of the nucleus. << 133 // G4cout << "BuildStrings prim/targ " << primaries.entries() << " , " <<
2681 // The nucleons that have been involved in << 134 // targets.entries() << G4endl;
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 135
2830 136
2831 //=========================================== << 137 unsigned int ahadron;
>> 138 for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
>> 139 {
>> 140 G4bool isProjectile=true;
>> 141 strings->push_back(theExcitation->String(primaries[ahadron], isProjectile));
>> 142 }
>> 143 for ( ahadron=0; ahadron < targets.size() ; ahadron++)
>> 144 {
>> 145 G4bool isProjectile=false;
>> 146 strings->push_back(theExcitation->String(targets[ahadron], isProjectile));
>> 147 }
2832 148
2833 G4bool G4FTFModel:: << 149 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
2834 SamplingNucleonKinematics( G4double averagePt << 150 primaries.clear();
2835 const G4double max << 151 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
2836 G4double dCor, << 152 targets.clear();
2837 G4V3DNucleus* nucl << 153
2838 const G4LorentzVec << 154 return strings;
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 } 155 }
2985 156
>> 157 G4bool G4FTFModel::ExciteParticipants()
>> 158 {
>> 159
>> 160 while (theParticipants.Next())
>> 161 {
>> 162 const G4InteractionContent & collision=theParticipants.GetInteraction();
>> 163
>> 164 //G4cout << " soft colls : " << collision.GetNumberOfSoftCollisions() << G4endl; // Uzhi no match
>> 165 G4VSplitableHadron * projectile=collision.GetProjectile();
>> 166 G4VSplitableHadron * target=collision.GetTarget();
>> 167
>> 168 if ( ! theExcitation->ExciteParticipants(projectile, target) )
>> 169 {
>> 170 // give up, clean up
>> 171 std::vector<G4VSplitableHadron *> primaries;
>> 172 std::vector<G4VSplitableHadron *> targets;
>> 173 theParticipants.StartLoop(); // restart a loop
>> 174 while ( theParticipants.Next() )
>> 175 {
>> 176 const G4InteractionContent & interaction=theParticipants.GetInteraction();
>> 177 // do not allow for duplicates ...
>> 178 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) )
>> 179 primaries.push_back(interaction.GetProjectile());
>> 180
>> 181 if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) )
>> 182 targets.push_back(interaction.GetTarget());
>> 183 }
>> 184 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
>> 185 primaries.clear();
>> 186 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
>> 187 targets.clear();
2986 188
2987 //=========================================== << 189 return false;
2988 << 190 }
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 191
3151 void G4FTFModel::ModelDescription( std::ostre << 192 }
3152 desc << " FTF (Fritiof) Mod << 193 return true;
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 } 194 }
3192 <<
3193 195