00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include <utility>
00040
00041 #include "G4FTFModel.hh"
00042 #include "G4ios.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4FTFParameters.hh"
00046 #include "G4FTFParticipants.hh"
00047 #include "G4DiffractiveSplitableHadron.hh"
00048 #include "G4InteractionContent.hh"
00049 #include "G4LorentzRotation.hh"
00050 #include "G4ParticleDefinition.hh"
00051 #include "G4ParticleTable.hh"
00052 #include "G4IonTable.hh"
00053
00054
00055
00056 G4FTFModel::G4FTFModel(const G4String& modelName):G4VPartonStringModel(modelName),
00057 theExcitation(new G4DiffractiveExcitation()),
00058 theElastic(new G4ElasticHNScattering()),
00059 theAnnihilation(new G4FTFAnnihilation())
00060 {
00061 G4VPartonStringModel::SetThisPointer(this);
00062 theParameters=0;
00063 NumberOfInvolvedNucleon=0;
00064 NumberOfInvolvedNucleonOfProjectile=0;
00065 SetEnergyMomentumCheckLevels(2*perCent, 150*MeV);
00066 }
00067
00068 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
00069
00070 G4FTFModel::~G4FTFModel()
00071 {
00072
00073
00074
00075 if( theParameters != 0 ) delete theParameters;
00076 if( theExcitation != 0 ) delete theExcitation;
00077 if( theElastic != 0 ) delete theElastic;
00078 if( theAnnihilation != 0 ) delete theAnnihilation;
00079
00080
00081 if(theAdditionalString.size() != 0)
00082 {
00083 std::for_each(theAdditionalString.begin(), theAdditionalString.end(),
00084 DeleteVSplitableHadron());
00085 }
00086 theAdditionalString.clear();
00087
00088
00089 if( NumberOfInvolvedNucleon != 0)
00090 {
00091 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
00092 {
00093 G4VSplitableHadron * aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
00094 if(aNucleon) delete aNucleon;
00095 }
00096 }
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 }
00110
00111
00112 void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
00113 {
00114 theProjectile = aProjectile;
00115
00116 G4double PlabPerParticle(0.);
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 theParticipants.SetProjectileNucleus(0);
00128 theParticipants.Init(aNucleus.GetA_asInt(),aNucleus.GetZ_asInt());
00129
00130 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
00131 {
00132 PlabPerParticle=theProjectile.GetMomentum().z();
00133
00134
00135
00136 }
00137
00138
00139 if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
00140 {
00141 theParticipants.InitProjectileNucleus(
00142 theProjectile.GetDefinition()->GetBaryonNumber(),
00143 (G4int) theProjectile.GetDefinition()->GetPDGCharge() );
00144
00145 G4ThreeVector BoostVector=theProjectile.GetMomentum()/theProjectile.GetTotalEnergy();
00146 theParticipants.theProjectileNucleus->DoLorentzBoost(BoostVector);
00147
00148 PlabPerParticle=theProjectile.GetMomentum().z()/
00149 theProjectile.GetDefinition()->GetBaryonNumber();
00150
00151
00152
00153 }
00154
00155 if(theProjectile.GetDefinition()->GetBaryonNumber() < -1)
00156 {
00157 theParticipants.InitProjectileNucleus(
00158 std::abs( theProjectile.GetDefinition()->GetBaryonNumber()),
00159 std::abs((G4int) theProjectile.GetDefinition()->GetPDGCharge()) );
00160
00161 G4ThreeVector BoostVector=theProjectile.GetMomentum()/theProjectile.GetTotalEnergy();
00162
00163 theParticipants.theProjectileNucleus->StartLoop();
00164 G4Nucleon * aNucleon;
00165 while ( (aNucleon = theParticipants.theProjectileNucleus->GetNextNucleon()) )
00166 {
00167 if(aNucleon->GetDefinition()->GetPDGEncoding() == 2212)
00168 {aNucleon->SetParticleType(G4AntiProton::AntiProton());}
00169
00170 if(aNucleon->GetDefinition()->GetPDGEncoding() == 2112)
00171 {aNucleon->SetParticleType(G4AntiNeutron::AntiNeutron());}
00172 }
00173
00174 theParticipants.theProjectileNucleus->DoLorentzBoost(BoostVector);
00175
00176 PlabPerParticle= theProjectile.GetMomentum().z()/
00177 std::abs(theProjectile.GetDefinition()->GetBaryonNumber());
00178
00179
00180
00181
00182 }
00183
00184
00185 if( theParameters != 0 ) delete theParameters;
00186 theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
00187 aNucleus.GetA_asInt(),aNucleus.GetZ_asInt(),
00188 PlabPerParticle);
00189
00190
00191
00192
00193
00194
00195
00196 if(theAdditionalString.size() != 0)
00197 {
00198 std::for_each(theAdditionalString.begin(), theAdditionalString.end(),
00199 DeleteVSplitableHadron());
00200 }
00201 theAdditionalString.clear();
00202
00203 }
00204
00205
00206 G4ExcitedStringVector * G4FTFModel::GetStrings()
00207 {
00208 G4ExcitedStringVector * theStrings(0);
00209
00210 theParticipants.GetList(theProjectile,theParameters);
00211
00212
00213 G4bool Success(true);
00214
00215 if((std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1) &&
00216 (theProjectile.GetDefinition()->GetBaryonNumber() != -1) )
00217 {
00218
00219 ReggeonCascade();
00220
00221 Success=PutOnMassShell();
00222
00223 GetResidualNucleus();
00224 }
00225
00226
00227 if(theProjectile.GetDefinition()->GetBaryonNumber() > 1)
00228 {
00229
00230 StoreInvolvedNucleon();
00231 ReggeonCascade();
00232 Success=PutOnMassShell();
00233 GetResidualNucleus();
00234 }
00235
00236
00237 if(theProjectile.GetDefinition()->GetBaryonNumber() <= -1)
00238 {
00239
00240
00241 StoreInvolvedNucleon();
00242 if(theProjectile.GetTotalMomentum()/
00243 std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 5000.*MeV)
00244 {
00245
00246
00247 ReggeonCascade();
00248
00249 Success=PutOnMassShell();
00250
00251 GetResidualNucleus();
00252 }
00253 else
00254 {
00255
00256
00257 Success=true;
00258 }
00259 }
00260
00261 Success=Success && ExciteParticipants();
00262
00263
00264
00265 if( Success )
00266 {
00267 theStrings = BuildStrings();
00268
00269 if( theParameters != 0 )
00270 {
00271 delete theParameters;
00272 theParameters=0;
00273 }
00274 }
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295 if(!Success)
00296 {
00297
00298
00299 std::vector<G4VSplitableHadron *> primaries;
00300
00301 theParticipants.StartLoop();
00302 while ( theParticipants.Next() )
00303 {
00304 const G4InteractionContent & interaction=theParticipants.GetInteraction();
00305
00306
00307 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
00308 interaction.GetProjectile()) )
00309 primaries.push_back(interaction.GetProjectile());
00310 }
00311 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
00312 primaries.clear();
00313 }
00314
00315 G4VSplitableHadron * aNucleon = 0;
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
00331 {
00332 aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
00333 if(aNucleon) delete aNucleon;
00334 }
00335
00336 NumberOfInvolvedNucleon=0;
00337
00338 return theStrings;
00339
00340 }
00341
00342
00343 void G4FTFModel::StoreInvolvedNucleon()
00344 {
00345 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) <= 1)
00346 {
00347
00348 NumberOfInvolvedNucleon=0;
00349
00350 theParticipants.StartLoop();
00351
00352 while (theParticipants.Next())
00353 {
00354
00355 const G4InteractionContent & collision=theParticipants.GetInteraction();
00356
00357 G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
00358
00359
00360 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
00361 NumberOfInvolvedNucleon++;
00362
00363 }
00364
00365 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
00366
00367
00368 theParticipants.StartLoop();
00369
00370 theParticipants.Next();
00371 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
00372
00373
00374
00375
00376 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
00377
00378 primary->SetTimeOfCreation(0.);
00379
00380 G4double ZcoordinateOfPreviousCollision(0.);
00381 G4double ZcoordinateOfCurrentInteraction(0.);
00382 G4double TimeOfPreviousCollision(0.);
00383 G4double TimeOfCurrentCollision(0);
00384
00385 theParticipants.theNucleus->StartLoop();
00386 G4Nucleon * aNucleon;
00387 G4bool theFirstInvolvedNucleon(true);
00388 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
00389 {
00390
00391 if(aNucleon->AreYouHit())
00392 {
00393 if(theFirstInvolvedNucleon)
00394 {
00395 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
00396
00397 theFirstInvolvedNucleon=false;
00398 }
00399
00400 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
00401
00402
00403
00404
00405
00406 if ( betta_z > 1.0e-10 ) {
00407 TimeOfCurrentCollision=TimeOfPreviousCollision+
00408 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
00409 } else {
00410 TimeOfCurrentCollision=TimeOfPreviousCollision;
00411 }
00412
00413
00414
00415 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
00416
00417 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
00418 TimeOfPreviousCollision=TimeOfCurrentCollision;
00419 }
00420 }
00421
00422 return;
00423 }
00424
00425
00426
00427 NumberOfInvolvedNucleonOfProjectile=0;
00428
00429 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
00430 ProjectileNucleus->StartLoop();
00431
00432 G4Nucleon * ProjectileNucleon;
00433 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
00434 {
00435 if ( ProjectileNucleon->AreYouHit() )
00436 {
00437 TheInvolvedNucleonOfProjectile[NumberOfInvolvedNucleonOfProjectile]=
00438 ProjectileNucleon;
00439 NumberOfInvolvedNucleonOfProjectile++;
00440 }
00441 }
00442
00443
00444 NumberOfInvolvedNucleon=0;
00445
00446 G4V3DNucleus * TargetNucleus =theParticipants.GetWoundedNucleus();
00447 TargetNucleus->StartLoop();
00448
00449 G4Nucleon * TargetNucleon;
00450 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
00451 {
00452 if ( TargetNucleon->AreYouHit() )
00453 {
00454 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
00455 NumberOfInvolvedNucleon++;
00456 }
00457 }
00458
00459
00460 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
00461 return;
00462 }
00463
00464
00465 void G4FTFModel::ReggeonCascade()
00466 {
00467
00468
00469 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 1) return;
00470
00471
00472
00473 NumberOfInvolvedNucleon=0;
00474
00475 theParticipants.StartLoop();
00476 while (theParticipants.Next())
00477 {
00478 const G4InteractionContent & collision=theParticipants.GetInteraction();
00479
00480 G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
00481
00482 TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
00483 NumberOfInvolvedNucleon++;
00484
00485 G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
00486 G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
00487
00488 theParticipants.theNucleus->StartLoop();
00489 G4Nucleon * Neighbour(0);
00490
00491 while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) )
00492 {
00493 if(!Neighbour->AreYouHit())
00494 {
00495 G4double impact2= sqr(XofWoundedNucleon - Neighbour->GetPosition().x()) +
00496 sqr(YofWoundedNucleon - Neighbour->GetPosition().y());
00497
00498 if(G4UniformRand() < theParameters->GetCofNuclearDestruction()*
00499 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction()))
00500 {
00501
00502 TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
00503 NumberOfInvolvedNucleon++;
00504
00505
00506 G4VSplitableHadron *targetSplitable;
00507 targetSplitable = new G4DiffractiveSplitableHadron(*Neighbour);
00508
00509 Neighbour->Hit(targetSplitable);
00510 targetSplitable->SetStatus(2);
00511 }
00512 }
00513 }
00514 }
00515
00516 NumberOfInvolvedTargetNucleon=NumberOfInvolvedNucleon;
00517
00518
00519 theParticipants.StartLoop();
00520 theParticipants.Next();
00521 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
00522 G4double betta_z=primary->Get4Momentum().pz()/primary->Get4Momentum().e();
00523 primary->SetTimeOfCreation(0.);
00524
00525 G4double ZcoordinateOfPreviousCollision(0.);
00526 G4double ZcoordinateOfCurrentInteraction(0.);
00527 G4double TimeOfPreviousCollision(0.);
00528 G4double TimeOfCurrentCollision(0);
00529
00530 theParticipants.theNucleus->StartLoop();
00531 G4Nucleon * aNucleon;
00532 G4bool theFirstInvolvedNucleon(true);
00533 while ( (aNucleon = theParticipants.theNucleus->GetNextNucleon()) )
00534 {
00535 if(aNucleon->AreYouHit())
00536 {
00537 if(theFirstInvolvedNucleon)
00538 {
00539 ZcoordinateOfPreviousCollision=aNucleon->GetPosition().z();
00540 theFirstInvolvedNucleon=false;
00541 }
00542
00543 ZcoordinateOfCurrentInteraction=aNucleon->GetPosition().z();
00544 TimeOfCurrentCollision=TimeOfPreviousCollision+
00545 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
00546
00547 aNucleon->GetSplitableHadron()->SetTimeOfCreation(TimeOfCurrentCollision);
00548
00549 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
00550 TimeOfPreviousCollision=TimeOfCurrentCollision;
00551 }
00552 }
00553
00554
00555
00556 }
00557
00558
00559 G4bool G4FTFModel::PutOnMassShell()
00560 {
00561
00562 if(std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) > 1)
00563 {
00564
00565 G4LorentzVector P_total(0.,0.,0.,0.);
00566
00567 G4LorentzVector P_participants(0.,0.,0.,0.);
00568 G4int MultiplicityOfParticipants(0);
00569
00570 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
00571 ProjectileNucleus->StartLoop();
00572
00573 G4Nucleon * ProjectileNucleon;
00574 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
00575 {
00576 if ( ProjectileNucleon->AreYouHit() )
00577 {
00578 P_total+=ProjectileNucleon->Get4Momentum();
00579 MultiplicityOfParticipants++;
00580 P_participants+=ProjectileNucleon->Get4Momentum();
00581 }
00582 }
00583
00584
00585 G4int ResidualMassNumber(0);
00586 G4int ResidualCharge(0);
00587 ResidualExcitationEnergy=0.;
00588 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
00589
00590 G4double ExcitationEnergyPerWoundedNucleon=
00591 theParameters->GetExcitationEnergyPerWoundedNucleon();
00592
00593 G4V3DNucleus * TargetNucleus =theParticipants.GetWoundedNucleus();
00594 TargetNucleus->StartLoop();
00595
00596 G4Nucleon * TargetNucleon;
00597 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
00598 {
00599 P_total+=TargetNucleon->Get4Momentum();
00600 if ( TargetNucleon->AreYouHit() )
00601 {
00602 MultiplicityOfParticipants++;
00603 P_participants+=TargetNucleon->Get4Momentum();
00604 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
00605 }
00606 else
00607 {
00608 ResidualMassNumber+=1;
00609 ResidualCharge+=(G4int) TargetNucleon->GetDefinition()->GetPDGCharge();
00610 PnuclearResidual+=TargetNucleon->Get4Momentum();
00611 }
00612 }
00613
00614 G4double ResidualMass(0.);
00615 if(ResidualMassNumber == 0)
00616 {
00617 ResidualMass=0.;
00618 ResidualExcitationEnergy=0.;
00619 }
00620 else
00621 {
00622 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
00623 GetIonMass(ResidualCharge ,ResidualMassNumber);
00624 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
00625 }
00626
00627
00628
00629
00630
00631
00632
00633
00634 G4double SqrtS=P_total.mag();
00635 G4double S=P_total.mag2();
00636
00637
00638 {
00639 G4double MassOfParticipants=P_participants.mag();
00640 G4double MassOfParticipants2=sqr(MassOfParticipants);
00641
00642
00643 if(SqrtS < ResidualMass + MassOfParticipants) {return false;}
00644
00645 if(SqrtS < ResidualMass+ResidualExcitationEnergy + MassOfParticipants)
00646 {ResidualExcitationEnergy=0.;}
00647
00648 ResidualMass +=ResidualExcitationEnergy;
00649
00650
00651 G4double ResidualMass2=sqr(ResidualMass);
00652
00653 G4LorentzRotation toCms(-1*P_total.boostVector());
00654
00655 G4LorentzVector Pcms=toCms*P_participants;
00656
00657
00658 if ( Pcms.pz() <= 0. )
00659 {
00660
00661 return false;
00662 }
00663
00664 toCms.rotateZ(-1*Pcms.phi());
00665 toCms.rotateY(-1*Pcms.theta());
00666
00667
00668 G4LorentzRotation toLab(toCms.inverse());
00669
00670 G4double DecayMomentum2=
00671 sqr(S)+sqr(MassOfParticipants2)+ sqr(ResidualMass2) -
00672 2.*S*MassOfParticipants2 - 2.*S*ResidualMass2
00673 -2.*MassOfParticipants2*ResidualMass2;
00674
00675 if(DecayMomentum2 < 0.) return false;
00676
00677 DecayMomentum2/=(4.*S);
00678 G4double DecayMomentum = std::sqrt(DecayMomentum2);
00679
00680
00681 G4LorentzVector New_P_participants(0.,0.,DecayMomentum,
00682 std::sqrt(DecayMomentum2+MassOfParticipants2));
00683 G4LorentzVector New_PnuclearResidual(0.,0.,-DecayMomentum,
00684 std::sqrt(DecayMomentum2+ResidualMass2));
00685
00686
00687
00688 New_P_participants.transform(toLab);
00689 New_PnuclearResidual.transform(toLab);
00690
00691
00692
00693
00694 G4LorentzVector DeltaP_participants=(New_P_participants - P_participants)/
00695 ((G4double) MultiplicityOfParticipants);
00696
00697
00698
00699 ProjectileNucleus->StartLoop();
00700 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
00701 {
00702 if ( ProjectileNucleon->AreYouHit() )
00703 {
00704 G4LorentzVector Ptmp=ProjectileNucleon->Get4Momentum() + DeltaP_participants;
00705 ProjectileNucleon->GetSplitableHadron()->Set4Momentum(Ptmp);
00706 ProjectileNucleon->SetMomentum(Ptmp);
00707 }
00708 }
00709
00710
00711 TargetNucleus->StartLoop();
00712 while ( (TargetNucleon=TargetNucleus->GetNextNucleon()) )
00713 {
00714 if ( TargetNucleon->AreYouHit() )
00715 {
00716 G4LorentzVector Ptmp=TargetNucleon->Get4Momentum() + DeltaP_participants;
00717 TargetNucleon->GetSplitableHadron()->Set4Momentum(Ptmp);
00718 }
00719 }
00720
00721 Residual4Momentum = New_PnuclearResidual;
00722
00723 }
00724
00725 return true;
00726 }
00727
00728
00729
00730 theParticipants.StartLoop();
00731 theParticipants.Next();
00732 G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
00733 G4LorentzVector Pprojectile=primary->Get4Momentum();
00734
00735
00736
00737
00738
00739
00740 if(Pprojectile.z() < 0.){return false;}
00741
00742 G4double Mprojectile = Pprojectile.mag();
00743 G4double M2projectile = Pprojectile.mag2();
00744
00745 G4LorentzVector Psum = Pprojectile;
00746
00747 G4double SumMasses = Mprojectile + 20.*MeV;
00748
00749
00750
00751
00752 G4V3DNucleus *theNucleus = GetWoundedNucleus();
00753 G4int ResidualMassNumber=theNucleus->GetMassNumber();
00754 G4int ResidualCharge =theNucleus->GetCharge();
00755
00756 ResidualExcitationEnergy=0.;
00757 G4LorentzVector Ptarget(0.,0.,0.,0.);
00758 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
00759
00760 G4double ExcitationEnergyPerWoundedNucleon=
00761 theParameters->GetExcitationEnergyPerWoundedNucleon();
00762
00763
00764
00765 theNucleus->StartLoop();
00766
00767 while (G4Nucleon * aNucleon = theNucleus->GetNextNucleon())
00768 {
00769 Ptarget+=aNucleon->Get4Momentum();
00770
00771 if(aNucleon->AreYouHit())
00772 {
00773
00774
00775
00776
00777 SumMasses += std::sqrt(sqr(aNucleon->GetDefinition()->GetPDGMass())
00778 + aNucleon->Get4Momentum().perp2());
00779 SumMasses += 20.*MeV;
00780 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791 ResidualMassNumber--;
00792 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
00793 }
00794 else
00795 {
00796 PnuclearResidual += aNucleon->Get4Momentum();
00797 }
00798 }
00799
00800 Psum += Ptarget;
00801 PnuclearResidual.setPz(0.); PnuclearResidual.setE(0.);
00802
00803
00804 G4double ResidualMass(0.);
00805 if(ResidualMassNumber == 0)
00806 {
00807 ResidualMass=0.;
00808 ResidualExcitationEnergy=0.;
00809 }
00810 else
00811 {
00812 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->
00813 GetIonMass(ResidualCharge ,ResidualMassNumber);
00814 if(ResidualMassNumber == 1) {ResidualExcitationEnergy=0.;}
00815 }
00816
00817
00818
00819 SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.perp2());
00820
00821
00822
00823
00824 G4double SqrtS=Psum.mag();
00825 G4double S=Psum.mag2();
00826
00827
00828
00829 if(SqrtS < SumMasses) {return false;}
00830
00831
00832
00833 SumMasses -= std::sqrt(sqr(ResidualMass)+PnuclearResidual.perp2());
00834 SumMasses += std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
00835 +PnuclearResidual.perp2());
00836 if(SqrtS < SumMasses)
00837 {
00838 SumMasses -= std::sqrt(sqr(ResidualMass+ResidualExcitationEnergy)
00839 +PnuclearResidual.perp2());
00840 SumMasses += std::sqrt(sqr(ResidualMass)+PnuclearResidual.perp2());
00841 ResidualExcitationEnergy=0.;
00842 }
00843
00844 ResidualMass +=ResidualExcitationEnergy;
00845
00846
00847
00848
00849 G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV));
00850 G4int NumberOfDeltas(0);
00851
00852 if(theNucleus->GetMassNumber() != 1)
00853 {
00854
00855 G4double ProbDeltaIsobar(0.05);
00856 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
00857 {
00858 if((G4UniformRand() < ProbDeltaIsobar)&&(NumberOfDeltas < MaxNumberOfDeltas))
00859 {
00860 NumberOfDeltas++;
00861 G4VSplitableHadron * targetSplitable=TheInvolvedNucleon[i]->GetSplitableHadron();
00862 G4double MassDec=std::sqrt(sqr(targetSplitable->GetDefinition()->GetPDGMass())
00863 + targetSplitable->Get4Momentum().perp2());
00864
00865 G4int PDGcode = targetSplitable->GetDefinition()->GetPDGEncoding();
00866 G4ParticleDefinition* Old_def = targetSplitable->GetDefinition();
00867
00868 G4int newPDGcode = PDGcode/10; newPDGcode=newPDGcode*10+4;
00869 G4ParticleDefinition* ptr =
00870 G4ParticleTable::GetParticleTable()->FindParticle(newPDGcode);
00871 targetSplitable->SetDefinition(ptr);
00872 G4double MassInc=std::sqrt(sqr(targetSplitable->GetDefinition()->GetPDGMass())
00873 + targetSplitable->Get4Momentum().perp2());
00874 if(SqrtS < SumMasses + MassInc - MassDec)
00875 {
00876 targetSplitable->SetDefinition(Old_def);
00877 ProbDeltaIsobar = 0.;
00878 } else
00879 {
00880 SumMasses += (MassInc - MassDec);
00881 }
00882 }
00883 }
00884 }
00885
00886
00887
00888 G4LorentzRotation toCms(-1*Psum.boostVector());
00889 G4LorentzVector Ptmp=toCms*Pprojectile;
00890 if ( Ptmp.pz() <= 0. )
00891 {
00892
00893 return false;
00894 }
00895
00896 G4LorentzRotation toLab(toCms.inverse());
00897
00898 Ptmp=toCms*Ptarget;
00899 G4double YtargetNucleus=Ptmp.rapidity();
00900
00901
00902
00903 G4double Dcor = theParameters->GetDofNuclearDestruction()/
00904 theNucleus->GetMassNumber();
00905
00906 G4double AveragePt2 = theParameters->GetPt2ofNuclearDestruction();
00907 G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
00908
00909 G4double M2target(0.);
00910 G4double WminusTarget(0.);
00911 G4double WplusProjectile(0.);
00912
00913 G4int NumberOfTries(0);
00914 G4double ScaleFactor(1.);
00915 G4bool OuterSuccess(true);
00916 do
00917 {
00918 OuterSuccess=true;
00919
00920 do
00921 {
00922
00923 NumberOfTries++;
00924
00925 if(NumberOfTries == 100*(NumberOfTries/100))
00926 {
00927 ScaleFactor/=2.;
00928 Dcor *=ScaleFactor;
00929 AveragePt2 *=ScaleFactor;
00930 }
00931
00932 G4ThreeVector PtSum(0.,0.,0.);
00933 G4double XminusSum(0.);
00934 G4double Xminus(0.);
00935 G4bool InerSuccess=true;
00936
00937 do
00938 {
00939 InerSuccess=true;
00940
00941 PtSum =G4ThreeVector(0.,0.,0.);
00942 XminusSum=0.;
00943
00944 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
00945 {
00946 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
00947
00948 G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare);
00949 PtSum += tmpPt;
00950
00951 G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.);
00952 Xminus=tmpX.x();
00953 XminusSum+=Xminus;
00954
00955 G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,
00956 aNucleon->Get4Momentum().e());
00957
00958 aNucleon->SetMomentum(tmp);
00959 }
00960
00961
00962 G4double DeltaX = (PtSum.x()-PnuclearResidual.x())/NumberOfInvolvedNucleon;
00963 G4double DeltaY = (PtSum.y()-PnuclearResidual.y())/NumberOfInvolvedNucleon;
00964 G4double DeltaXminus(0.);
00965
00966 if(ResidualMassNumber == 0)
00967 {
00968 DeltaXminus = (XminusSum-1.)/NumberOfInvolvedNucleon;
00969 }
00970 else
00971 {
00972 DeltaXminus = -1./theNucleus->GetMassNumber();
00973 }
00974
00975 XminusSum=1.;
00976 M2target =0.;
00977
00978 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
00979 {
00980 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
00981
00982 Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
00983 XminusSum-=Xminus;
00984
00985 if(ResidualMassNumber == 0)
00986 {
00987 if((Xminus <= 0.) || (Xminus > 1.)) {InerSuccess=false; break;}
00988 } else
00989 {
00990 if((Xminus <= 0.) || (Xminus > 1.) ||
00991 (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
00992 }
00993
00994 G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
00995 G4double Py=aNucleon->Get4Momentum().py() - DeltaY;
00996
00997
00998
00999 M2target +=(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
01000 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() +
01001 Px*Px + Py*Py)/Xminus;
01002
01003
01004
01005
01006
01007
01008
01009
01010 G4LorentzVector tmp(Px,Py,Xminus,aNucleon->Get4Momentum().e());
01011 aNucleon->SetMomentum(tmp);
01012 }
01013
01014 if(InerSuccess && (ResidualMassNumber != 0))
01015 {
01016 M2target +=(sqr(ResidualMass) + PnuclearResidual.perp2())/XminusSum;
01017 }
01018
01019
01020 } while(!InerSuccess);
01021 } while (SqrtS < Mprojectile + std::sqrt(M2target));
01022
01023 G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
01024 -2.*S*M2projectile - 2.*S*M2target
01025 -2.*M2projectile*M2target;
01026
01027 WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
01028 WplusProjectile=SqrtS - M2target/WminusTarget;
01029
01030 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
01031 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
01032 G4double Yprojectile=0.5*std::log((Eprojectile+Pzprojectile)/
01033 (Eprojectile-Pzprojectile));
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
01044 {
01045 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
01046 G4LorentzVector tmp=aNucleon->Get4Momentum();
01047
01048 G4double Mt2(0.);
01049
01050
01051
01052 Mt2 = sqr(tmp.x())+sqr(tmp.y())+
01053 sqr(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass());
01054
01055
01056
01057
01058
01059
01060
01061 G4double Xminus=tmp.z();
01062
01063 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
01064 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
01065 G4double YtargetNucleon=0.5*std::log((E+Pz)/(E-Pz));
01066
01067
01068
01069
01070 if((std::abs(YtargetNucleon-YtargetNucleus) > 2) ||
01071 (Yprojectile < YtargetNucleon)) {OuterSuccess=false; break;}
01072
01073 }
01074
01075 } while(!OuterSuccess);
01076
01077
01078 G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
01079 G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
01080 Pprojectile.setPz(Pzprojectile); Pprojectile.setE(Eprojectile);
01081
01082
01083 Pprojectile.transform(toLab);
01084 primary->Set4Momentum(Pprojectile);
01085
01086
01087
01088 G4ThreeVector Residual3Momentum(0.,0.,1.);
01089
01090 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
01091 {
01092 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
01093 G4LorentzVector tmp=aNucleon->Get4Momentum();
01094
01095 Residual3Momentum-=tmp.vect();
01096
01097 G4double Mt2(0.);
01098
01099
01100
01101 Mt2 = sqr(tmp.x())+sqr(tmp.y())+
01102 sqr(aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass());
01103
01104
01105
01106
01107
01108
01109
01110 G4double Xminus=tmp.z();
01111
01112 G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
01113 G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
01114
01115 tmp.setPz(Pz);
01116 tmp.setE(E);
01117
01118 tmp.transform(toLab);
01119
01120 aNucleon->SetMomentum(tmp);
01121
01122 G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron();
01123 targetSplitable->Set4Momentum(tmp);
01124
01125 }
01126
01127
01128 G4double Mt2Residual=sqr(ResidualMass) +
01129 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y());
01130
01131
01132 G4double PzResidual=0.;
01133 G4double EResidual =0.;
01134 if(ResidualMassNumber != 0)
01135 {
01136 PzResidual=-WminusTarget*Residual3Momentum.z()/2. +
01137 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
01138 EResidual = WminusTarget*Residual3Momentum.z()/2. +
01139 Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
01140 }
01141
01142 Residual4Momentum.setPx(Residual3Momentum.x());
01143 Residual4Momentum.setPy(Residual3Momentum.y());
01144 Residual4Momentum.setPz(PzResidual);
01145 Residual4Momentum.setE(EResidual);
01146
01147
01148 Residual4Momentum.transform(toLab);
01149
01150
01151
01152 return true;
01153 }
01154
01155
01156 G4bool G4FTFModel::ExciteParticipants()
01157 {
01158
01159 G4bool Successfull(true);
01160
01161 theParticipants.StartLoop();
01162 G4int CurrentInteraction(0);
01163
01164 G4int MaxNumOfInelCollisions=G4int(theParameters->GetMaxNumberOfCollisions());
01165
01166 G4double NumberOfInel(0.);
01167
01168 if(MaxNumOfInelCollisions > 0)
01169 {
01170 G4double ProbMaxNumber=theParameters->GetMaxNumberOfCollisions()-
01171 MaxNumOfInelCollisions;
01172 if(G4UniformRand() < ProbMaxNumber) {MaxNumOfInelCollisions++;}
01173 NumberOfInel=MaxNumOfInelCollisions;
01174
01175 } else
01176 {
01177
01178 if(theParticipants.theNucleus->GetMassNumber() > 1)
01179 {
01180 NumberOfInel = theParameters->GetProbOfInteraction();
01181 MaxNumOfInelCollisions = 1;
01182 } else
01183 {
01184 NumberOfInel = 1.;
01185 MaxNumOfInelCollisions = 1;
01186
01187 }
01188 }
01189
01190
01191
01192 while (theParticipants.Next())
01193 {
01194 CurrentInteraction++;
01195 const G4InteractionContent & collision=theParticipants.GetInteraction();
01196
01197 G4VSplitableHadron * projectile=collision.GetProjectile();
01198 G4VSplitableHadron * target=collision.GetTarget();
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210 if(collision.GetStatus())
01211 {
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
01223 {
01224
01225 Successfull = theElastic->ElasticScattering(projectile, target, theParameters)
01226 || Successfull;
01227 }
01228 else if(G4UniformRand() > theParameters->GetProbabilityOfAnnihilation())
01229 {
01230
01231
01232 if(G4UniformRand()< NumberOfInel/MaxNumOfInelCollisions)
01233 {
01234 if(theExcitation->ExciteParticipants(projectile, target,
01235 theParameters, theElastic))
01236 {
01237 NumberOfInel--;
01238
01239
01240
01241 } else
01242 {
01243
01244
01245 Successfull = theElastic->ElasticScattering(projectile, target, theParameters)
01246 || Successfull;
01247
01248
01249
01250
01251
01252
01253
01254
01255 }
01256 } else
01257 {
01258
01259 Successfull = theElastic->ElasticScattering(projectile, target, theParameters)
01260 || Successfull;
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278 }
01279 }
01280 else
01281 {
01282
01283
01284
01285
01286
01287 {
01288 while (theParticipants.Next())
01289 {
01290 G4InteractionContent & acollision=theParticipants.GetInteraction();
01291
01292 G4VSplitableHadron * NextProjectileNucleon=acollision.GetProjectile();
01293 G4VSplitableHadron * NextTargetNucleon =acollision.GetTarget();
01294 if((projectile == NextProjectileNucleon) ||
01295 (target == NextTargetNucleon )) acollision.SetStatus(0);
01296
01297 }
01298
01299 theParticipants.StartLoop();
01300 for(G4int I=0; I < CurrentInteraction; I++) theParticipants.Next();
01301
01302
01303
01304
01305
01306
01307 }
01308 G4VSplitableHadron *AdditionalString=0;
01309 if(theAnnihilation->Annihilate(projectile, target, AdditionalString, theParameters))
01310 {
01311 Successfull = Successfull || true;
01312
01313
01314
01315
01316 if(AdditionalString != 0) theAdditionalString.push_back(AdditionalString);
01317
01318
01319
01320 } else
01321 {
01322
01323
01324
01325
01326
01327 }
01328 }
01329
01330 }
01331
01332 }
01333
01334
01335
01336
01337 return Successfull;
01338 }
01339
01340
01341 void G4FTFModel::AjustTargetNucleonForAnnihilation(G4VSplitableHadron *SelectedAntiBaryon,
01342 G4VSplitableHadron *SelectedTargetNucleon)
01343 {
01344 G4LorentzVector Pparticipants=SelectedAntiBaryon->Get4Momentum()+
01345 SelectedTargetNucleon->Get4Momentum();
01346
01347 G4V3DNucleus *theNucleus = GetWoundedNucleus();
01348
01349
01350 ResidualExcitationEnergy=0.;
01351 G4int ResidualCharge =theNucleus->GetCharge();
01352 G4int ResidualMassNumber=theNucleus->GetMassNumber();
01353
01354 G4ThreeVector P3nuclearResidual(0.,0.,0.);
01355 G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
01356
01357
01358 G4double ExcitationEnergyPerWoundedNucleon=
01359 theParameters->GetExcitationEnergyPerWoundedNucleon();
01360
01361 G4Nucleon * aNucleon;
01362 theNucleus->StartLoop();
01363 G4int NumberOfHoles(0);
01364
01365 while ((aNucleon = theNucleus->GetNextNucleon()))
01366 {
01367 G4int CurrentStatus=0;
01368 if(aNucleon->AreYouHit())
01369 {
01370 if(aNucleon->GetSplitableHadron() == SelectedTargetNucleon)
01371 {CurrentStatus=1;}
01372 else
01373 {
01374 if(aNucleon->GetSplitableHadron()->GetSoftCollisionCount() == 0)
01375 {CurrentStatus=0;}
01376 else {CurrentStatus=1;}
01377 }
01378 }
01379
01380 if(CurrentStatus != 0)
01381 {
01382
01383 NumberOfHoles++;
01384 ResidualExcitationEnergy+=ExcitationEnergyPerWoundedNucleon;
01385 ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
01386 ResidualMassNumber--;
01387 }
01388 else
01389 {
01390 PnuclearResidual+=aNucleon->Get4Momentum();
01391 }
01392 }
01393
01394
01395
01396 G4LorentzVector Psum=Pparticipants + PnuclearResidual;
01397
01398
01399 G4LorentzRotation toCms(-1*Psum.boostVector());
01400
01401 G4LorentzVector Ptmp=toCms*Psum;
01402
01403 toCms.rotateZ(-1*Ptmp.phi());
01404 toCms.rotateY(-1*Ptmp.theta());
01405
01406 G4LorentzRotation toLab(toCms.inverse());
01407
01408
01409 G4double SqrtS=Psum.mag();
01410 G4double S =sqr(SqrtS);
01411
01412 G4double ResidualMass(0.);
01413 if(ResidualMassNumber != 0)
01414 {
01415 ResidualMass=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(
01416 ResidualCharge,ResidualMassNumber);
01417 } else {return;}
01418
01419
01420
01421 if(ResidualMass > SqrtS) {return;}
01422 else
01423 {
01424 if(ResidualMass+ResidualExcitationEnergy > SqrtS)
01425 ResidualExcitationEnergy = SqrtS-ResidualMass;
01426 }
01427
01428 ResidualMass+=ResidualExcitationEnergy;
01429 G4double ResidualMass2=sqr(ResidualMass);
01430
01431
01432
01433 G4double ParticipantMass=Pparticipants.mag();
01434 G4double ParticipantMass2=sqr(ParticipantMass);
01435
01436 if(ResidualMass + ParticipantMass > SqrtS) ParticipantMass=SqrtS-ResidualMass;
01437
01438
01439
01440
01441 G4double DecayMomentum2=
01442 sqr(S)+sqr(ParticipantMass2)+ sqr(ResidualMass2) -
01443 2.*S*ParticipantMass2 - 2.*S*ResidualMass2
01444 -2.*ParticipantMass2*ResidualMass2;
01445
01446 if(DecayMomentum2 < 0.) return;
01447
01448 DecayMomentum2/=(4.*S);
01449 G4double DecayMomentum = std::sqrt(DecayMomentum2);
01450
01451
01452 G4LorentzVector New_Pparticipants(0.,0.,DecayMomentum,
01453 std::sqrt(DecayMomentum2+ParticipantMass2));
01454
01455 G4LorentzVector New_PnuclearResidual(0.,0.,-DecayMomentum,
01456 std::sqrt(DecayMomentum2+ResidualMass2));
01457
01458
01459
01460
01461 New_Pparticipants.transform(toLab);
01462 New_PnuclearResidual.transform(toLab);
01463
01464
01465
01466 G4LorentzVector DeltaP_participants=(Pparticipants - New_Pparticipants)/2.;
01467 G4LorentzVector DeltaP_nuclearResidual=(PnuclearResidual - New_PnuclearResidual)/
01468 (G4double) ResidualMassNumber;
01469
01470
01471 Ptmp=SelectedAntiBaryon->Get4Momentum() - DeltaP_participants;
01472 SelectedAntiBaryon->Set4Momentum(Ptmp);
01473
01474 Ptmp=SelectedTargetNucleon->Get4Momentum() - DeltaP_participants;
01475 SelectedTargetNucleon->Set4Momentum(Ptmp);
01476
01477
01478
01479
01480 G4double DeltaExcitationEnergy = 0.0;
01481 if ( NumberOfHoles != 0 ) {
01482 DeltaExcitationEnergy = ResidualExcitationEnergy / ((G4double) NumberOfHoles);
01483 }
01484
01485
01486 theNucleus->StartLoop();
01487 while ((aNucleon = theNucleus->GetNextNucleon()))
01488 {
01489 G4int CurrentStatus=0;
01490 if(aNucleon->AreYouHit())
01491 {
01492 if(aNucleon->GetSplitableHadron() == SelectedTargetNucleon)
01493 {CurrentStatus=1;}
01494 else
01495 {
01496 if(aNucleon->GetSplitableHadron()->GetSoftCollisionCount() == 0)
01497 {CurrentStatus=0;}
01498 else {CurrentStatus=1;}
01499 }
01500 }
01501
01502 if(CurrentStatus != 0)
01503 {
01504
01505 aNucleon->SetBindingEnergy(DeltaExcitationEnergy);
01506 }
01507 else
01508 {
01509 Ptmp=aNucleon->Get4Momentum() - DeltaP_nuclearResidual;
01510 aNucleon->SetMomentum(Ptmp);
01511 }
01512 }
01513
01514
01515 return;
01516 }
01517
01518
01519 G4ExcitedStringVector * G4FTFModel::BuildStrings()
01520 {
01521
01522
01523
01524 G4ExcitedStringVector * strings;
01525 strings = new G4ExcitedStringVector();
01526
01527 std::vector<G4VSplitableHadron *> primaries;
01528
01529 G4ExcitedString * FirstString(0);
01530 G4ExcitedString * SecondString(0);
01531
01532 theParticipants.StartLoop();
01533
01534 while ( theParticipants.Next() )
01535 {
01536 const G4InteractionContent & interaction=theParticipants.GetInteraction();
01537
01538 {
01539
01540 if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
01541 interaction.GetProjectile()) )
01542 primaries.push_back(interaction.GetProjectile());
01543 }
01544 }
01545
01546
01547 for (unsigned int ahadron=0; ahadron < primaries.size() ; ahadron++)
01548 {
01549 G4bool isProjectile(0);
01550
01551 if(primaries[ahadron]->GetStatus() <= 1) {isProjectile=true; }
01552
01553
01554 FirstString=0; SecondString=0;
01555 theExcitation->CreateStrings(primaries[ahadron], isProjectile,
01556 FirstString, SecondString,
01557 theParameters);
01558
01559 if(FirstString != 0) strings->push_back(FirstString);
01560 if(SecondString != 0) strings->push_back(SecondString);
01561
01562
01563
01564 }
01565
01566
01567
01568
01569 G4V3DNucleus * ProjectileNucleus =theParticipants.GetProjectileNucleus();
01570 if(ProjectileNucleus)
01571 {
01572 ProjectileNucleus->StartLoop();
01573
01574 G4Nucleon * ProjectileNucleon;
01575 while ( (ProjectileNucleon=ProjectileNucleus->GetNextNucleon()) )
01576 {
01577 if ( !ProjectileNucleon->AreYouHit() )
01578 {
01579
01580 G4VSplitableHadron * ProjectileSplitable=0;
01581 ProjectileSplitable= new G4DiffractiveSplitableHadron(*ProjectileNucleon);
01582 ProjectileNucleon->Hit(0);
01583
01584 G4bool isProjectile=true;
01585 FirstString=0; SecondString=0;
01586 theExcitation->CreateStrings(ProjectileSplitable,
01587 isProjectile,
01588 FirstString, SecondString,
01589 theParameters);
01590 if(FirstString != 0) strings->push_back(FirstString);
01591 if(SecondString != 0) strings->push_back(SecondString);
01592
01593 delete ProjectileSplitable;
01594 }
01595 }
01596 }
01597
01598
01599 if(theAdditionalString.size() != 0)
01600 {
01601 for (unsigned int ahadron=0; ahadron < theAdditionalString.size() ; ahadron++)
01602 {
01603 G4bool isProjectile(0);
01604
01605 if(theAdditionalString[ahadron]->GetStatus() <= 1) {isProjectile=true; }
01606
01607
01608 FirstString=0; SecondString=0;
01609 theExcitation->CreateStrings(theAdditionalString[ahadron], isProjectile,
01610 FirstString, SecondString,
01611 theParameters);
01612
01613 if(FirstString != 0) strings->push_back(FirstString);
01614 if(SecondString != 0) strings->push_back(SecondString);
01615
01616
01617 }
01618 }
01619
01620
01621
01622
01623 for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
01624 {
01625
01626 if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() ==4)
01627 {
01628
01629 delete TheInvolvedNucleon[ahadron]->GetSplitableHadron();
01630 G4VSplitableHadron *aHit=0;
01631 TheInvolvedNucleon[ahadron]->Hit(aHit);
01632 }
01633 else if((TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() <=1) &&
01634 (TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetSoftCollisionCount() ==0))
01635 {
01636
01637
01638 delete TheInvolvedNucleon[ahadron]->GetSplitableHadron();
01639 G4VSplitableHadron *aHit=0;
01640 TheInvolvedNucleon[ahadron]->Hit(aHit);
01641 }
01642 else if((TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() <=1) &&
01643 (TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetSoftCollisionCount() !=0))
01644 {
01645
01646 G4bool isProjectile=false;
01647 FirstString=0; SecondString=0;
01648 theExcitation->CreateStrings(
01649 TheInvolvedNucleon[ahadron]->GetSplitableHadron(),
01650 isProjectile,
01651 FirstString, SecondString,
01652 theParameters);
01653 if(FirstString != 0) strings->push_back(FirstString);
01654 if(SecondString != 0) strings->push_back(SecondString);
01655 }
01656 else if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() ==2)
01657 {
01658
01659 G4bool isProjectile=false;
01660 FirstString=0; SecondString=0;
01661 theExcitation->CreateStrings(
01662 TheInvolvedNucleon[ahadron]->GetSplitableHadron(),
01663 isProjectile,
01664 FirstString, SecondString,
01665 theParameters);
01666 if(FirstString != 0) strings->push_back(FirstString);
01667 if(SecondString != 0) strings->push_back(SecondString);
01668
01669 }
01670 else if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() ==3)
01671 {
01672
01673 TheInvolvedNucleon[ahadron]->SetBindingEnergy(theParameters->GetExcitationEnergyPerWoundedNucleon());
01674 }
01675 else {}
01676
01677 }
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
01688 primaries.clear();
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707 return strings;
01708 }
01709
01710 void G4FTFModel::GetResidualNucleus()
01711 {
01712 G4double DeltaExcitationE=ResidualExcitationEnergy/
01713 (G4double) NumberOfInvolvedNucleon;
01714 G4LorentzVector DeltaPResidualNucleus = Residual4Momentum/
01715 (G4double) NumberOfInvolvedNucleon;
01716
01717 for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
01718 {
01719 G4Nucleon * aNucleon = TheInvolvedNucleon[i];
01720
01721 G4LorentzVector tmp=-DeltaPResidualNucleus;
01722 aNucleon->SetMomentum(tmp);
01723 aNucleon->SetBindingEnergy(DeltaExcitationE);
01724 }
01725
01726 }
01727
01728
01729 G4ThreeVector G4FTFModel::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
01730 {
01731
01732 G4double Pt2(0.);
01733 if(AveragePt2 <= 0.) {Pt2=0.;}
01734 else
01735 {
01736 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
01737 (std::exp(-maxPtSquare/AveragePt2)-1.));
01738 }
01739 G4double Pt=std::sqrt(Pt2);
01740 G4double phi=G4UniformRand() * twopi;
01741
01742 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
01743 }
01744
01745 void G4FTFModel::ModelDescription(std::ostream& desc) const
01746 {
01747 desc << "please add description here" << G4endl;
01748 }