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