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