| 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 #include <utility> 26 #include <utility>
27 27
28 #include "G4QGSParticipants.hh" 28 #include "G4QGSParticipants.hh"
29 #include "globals.hh" 29 #include "globals.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4LorentzVector.hh" 31 #include "G4LorentzVector.hh"
32 #include "G4V3DNucleus.hh" 32 #include "G4V3DNucleus.hh"
33 #include "G4ParticleTable.hh" 33 #include "G4ParticleTable.hh"
34 #include "G4IonTable.hh" 34 #include "G4IonTable.hh"
35 #include "G4PhysicalConstants.hh" 35 #include "G4PhysicalConstants.hh"
36 36
37 #include "G4Exp.hh" 37 #include "G4Exp.hh"
38 #include "G4Log.hh" 38 #include "G4Log.hh"
39 #include "G4Pow.hh" 39 #include "G4Pow.hh"
40 40
41 //#define debugQGSParticipants 41 //#define debugQGSParticipants
42 //#define debugPutOnMassShell 42 //#define debugPutOnMassShell
43 43
44 // Class G4QGSParticipants 44 // Class G4QGSParticipants
45 45
46 // Promoting model parameters from local varia 46 // Promoting model parameters from local variables class properties
47 G4ThreadLocal G4int G4QGSParticipants_NPart = 47 G4ThreadLocal G4int G4QGSParticipants_NPart = 0;
48 48
49 G4QGSParticipants::G4QGSParticipants() << 49 G4QGSParticipants::G4QGSParticipants() : theDiffExcitaton(),
50 : theDiffExcitaton(), ModelMode(SOFT), nCutM << 50 ModelMode(SOFT),
51 ThresholdParameter(0.0*GeV), QGSMThreshold << 51 nCutMax(7),ThresholdParameter(0.000*GeV),
52 theNucleonRadius(1.5*fermi), theCurrentVel << 52 QGSMThreshold(3*GeV),theNucleonRadius(1.5*fermi),alpha(-0.5),beta(2.5)
53 theProjectileSplitable(nullptr), Regge(nul <<
54 InteractionMode(ALL), alpha(-0.5), beta(2. <<
55 NumberOfInvolvedNucleonsOfTarget(0), Numbe <<
56 ProjectileResidual4Momentum(G4LorentzVecto <<
57 ProjectileResidualCharge(0), ProjectileRes <<
58 TargetResidual4Momentum(G4LorentzVector()) <<
59 TargetResidualCharge(0), TargetResidualExc <<
60 CofNuclearDestruction(0.0), R2ofNuclearDes <<
61 ExcitationEnergyPerWoundedNucleon(0.0), Do <<
62 Pt2ofNuclearDestruction(0.0), MaxPt2ofNucl <<
63 { 53 {
64 for (size_t i=0; i < 250; ++i) { <<
65 TheInvolvedNucleonsOfTarget[i] = nullptr; <<
66 TheInvolvedNucleonsOfProjectile[i] = nullp <<
67 } <<
68 // Parameters setting 54 // Parameters setting
69 SetCofNuclearDestruction( 1.0 ); << 55 SetCofNuclearDestruction(1.);
70 SetR2ofNuclearDestruction( 1.5*fermi*fermi ) 56 SetR2ofNuclearDestruction( 1.5*fermi*fermi );
71 SetDofNuclearDestruction( 0.3 ); 57 SetDofNuclearDestruction( 0.3 );
72 SetPt2ofNuclearDestruction( 0.075*GeV*GeV ); 58 SetPt2ofNuclearDestruction( 0.075*GeV*GeV );
73 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV ) 59 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
74 SetExcitationEnergyPerWoundedNucleon( 40.0*M 60 SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
75 61
76 sigmaPt=0.25*sqr(GeV); 62 sigmaPt=0.25*sqr(GeV);
77 } 63 }
78 64
79 G4QGSParticipants::G4QGSParticipants(const G4Q 65 G4QGSParticipants::G4QGSParticipants(const G4QGSParticipants &right)
80 : G4VParticipants(), ModelMode(right.ModelMo << 66 : G4VParticipants(),ModelMode(right.ModelMode), nCutMax(right.nCutMax),
81 ThresholdParameter(right.ThresholdParamete << 67 ThresholdParameter(right.ThresholdParameter), QGSMThreshold(right.QGSMThreshold),
82 QGSMThreshold(right.QGSMThreshold), << 68 theNucleonRadius(right.theNucleonRadius)
83 theNucleonRadius(right.theNucleonRadius), << 69 {}
84 theCurrentVelocity(right.theCurrentVelocit <<
85 theProjectileSplitable(right.theProjectile <<
86 Regge(right.Regge), InteractionMode(right. <<
87 alpha(right.alpha), beta(right.beta), sigm <<
88 NumberOfInvolvedNucleonsOfTarget(right.Num <<
89 NumberOfInvolvedNucleonsOfProjectile(right <<
90 ProjectileResidual4Momentum(right.Projecti <<
91 ProjectileResidualMassNumber(right.Project <<
92 ProjectileResidualCharge(right.ProjectileR <<
93 ProjectileResidualExcitationEnergy(right.P <<
94 TargetResidual4Momentum(right.TargetResidu <<
95 TargetResidualMassNumber(right.TargetResid <<
96 TargetResidualCharge(right.TargetResidualC <<
97 TargetResidualExcitationEnergy(right.Targe <<
98 CofNuclearDestruction(right.CofNuclearDest <<
99 R2ofNuclearDestruction(right.R2ofNuclearDe <<
100 ExcitationEnergyPerWoundedNucleon(right.Ex <<
101 DofNuclearDestruction(right.DofNuclearDest <<
102 Pt2ofNuclearDestruction(right.Pt2ofNuclear <<
103 MaxPt2ofNuclearDestruction(right.MaxPt2ofN <<
104 { <<
105 for (size_t i=0; i < 250; ++i) { <<
106 TheInvolvedNucleonsOfTarget[i] = right.The <<
107 TheInvolvedNucleonsOfProjectile[i] = right <<
108 } <<
109 } <<
110 70
111 G4QGSParticipants::~G4QGSParticipants() {} 71 G4QGSParticipants::~G4QGSParticipants() {}
112 72
113 void G4QGSParticipants::BuildInteractions(cons 73 void G4QGSParticipants::BuildInteractions(const G4ReactionProduct &thePrimary)
114 { 74 {
115 theProjectile = thePrimary; 75 theProjectile = thePrimary;
116 76
117 Regge = new G4Reggeons(theProjectile.GetDefi 77 Regge = new G4Reggeons(theProjectile.GetDefinition());
118 78
119 SetProjectileNucleus( 0 ); 79 SetProjectileNucleus( 0 );
120 80
121 NumberOfInvolvedNucleonsOfProjectile= 0; 81 NumberOfInvolvedNucleonsOfProjectile= 0;
122 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 ); 82 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 );
123 ProjectileResidualMassNumber = 0; 83 ProjectileResidualMassNumber = 0;
124 ProjectileResidualCharge = 0; 84 ProjectileResidualCharge = 0;
125 ProjectileResidualExcitationEnergy = 0.0; 85 ProjectileResidualExcitationEnergy = 0.0;
126 ProjectileResidual4Momentum = tmp; 86 ProjectileResidual4Momentum = tmp;
127 87
128 NumberOfInvolvedNucleonsOfTarget= 0; 88 NumberOfInvolvedNucleonsOfTarget= 0;
129 TargetResidualMassNumber = theNucleus- 89 TargetResidualMassNumber = theNucleus->GetMassNumber();
130 TargetResidualCharge = theNucleus- 90 TargetResidualCharge = theNucleus->GetCharge();
131 TargetResidualExcitationEnergy = 0.0; 91 TargetResidualExcitationEnergy = 0.0;
132 92
133 theNucleus->StartLoop(); 93 theNucleus->StartLoop();
134 G4Nucleon* NuclearNucleon; 94 G4Nucleon* NuclearNucleon;
135 while ( ( NuclearNucleon = theNucleus->GetNe 95 while ( ( NuclearNucleon = theNucleus->GetNextNucleon() ) ) {
136 tmp+=NuclearNucleon->Get4Momentum(); 96 tmp+=NuclearNucleon->Get4Momentum();
137 } 97 }
138 TargetResidual4Momentum = tmp; 98 TargetResidual4Momentum = tmp;
139 99
140 if ( std::abs( theProjectile.GetDefinition() 100 if ( std::abs( theProjectile.GetDefinition()->GetBaryonNumber() ) <= 1 ) {
141 // Projectile is a hadron : meson or baryo 101 // Projectile is a hadron : meson or baryon
142 ProjectileResidualMassNumber = std::abs( t 102 ProjectileResidualMassNumber = std::abs( theProjectile.GetDefinition()->GetBaryonNumber() );
143 ProjectileResidualCharge = G4int( theProje 103 ProjectileResidualCharge = G4int( theProjectile.GetDefinition()->GetPDGCharge() );
144 ProjectileResidualExcitationEnergy = 0.0; 104 ProjectileResidualExcitationEnergy = 0.0;
145 ProjectileResidual4Momentum.setVect( thePr 105 ProjectileResidual4Momentum.setVect( theProjectile.GetMomentum() );
146 ProjectileResidual4Momentum.setE( theProje 106 ProjectileResidual4Momentum.setE( theProjectile.GetTotalEnergy() );
147 } 107 }
148 108
149 #ifdef debugQGSParticipants 109 #ifdef debugQGSParticipants
150 G4cout <<G4endl<< "G4QGSParticipants::Buil 110 G4cout <<G4endl<< "G4QGSParticipants::BuildInteractions---------" << G4endl
151 << "thePrimary " << thePrimary.GetD 111 << "thePrimary " << thePrimary.GetDefinition()->GetParticleName()<<" "
152 <<ProjectileResidual4Momentum<<G4en 112 <<ProjectileResidual4Momentum<<G4endl;
153 G4cout << "Target A and Z " << theNucleus 113 G4cout << "Target A and Z " << theNucleus->GetMassNumber() <<" "<<theNucleus->GetCharge()<<" "
154 << TargetResidual4Momentum<<G4endl; 114 << TargetResidual4Momentum<<G4endl;
155 #endif 115 #endif
156 116
157 G4bool Success( true ); 117 G4bool Success( true );
158 118
159 const G4int maxNumberOfLoops = 1000; 119 const G4int maxNumberOfLoops = 1000;
160 G4int loopCounter = 0; 120 G4int loopCounter = 0;
161 do 121 do
162 { 122 {
163 const G4int maxNumberOfInternalLoops = 100 123 const G4int maxNumberOfInternalLoops = 1000;
164 G4int internalLoopCounter = 0; 124 G4int internalLoopCounter = 0;
165 do 125 do
166 { 126 {
167 if(std::abs(theProjectile.GetDefinition( 127 if(std::abs(theProjectile.GetDefinition()->GetPDGEncoding()) < 100.0)
168 { 128 {
169 SelectInteractions(theProjectile); // 129 SelectInteractions(theProjectile); // for lepton projectile
170 } 130 }
171 else 131 else
172 { 132 {
173 GetList( theProjectile ); // Get list 133 GetList( theProjectile ); // Get list of participating nucleons for hadron projectile
174 } 134 }
175 135
176 if ( theInteractions.size() == 0 ) retur 136 if ( theInteractions.size() == 0 ) return;
177 137
178 StoreInvolvedNucleon(); // Store p 138 StoreInvolvedNucleon(); // Store participating nucleon
179 139
180 ReggeonCascade(); // Make re 140 ReggeonCascade(); // Make reggeon cascading. Involve nucleons in the cascading.
181 141
182 Success = PutOnMassShell(); // Ascribe 142 Success = PutOnMassShell(); // Ascribe momenta to the involved and participating nucleons
183 143
184 if(!Success) PrepareInitialState( thePri 144 if(!Success) PrepareInitialState( thePrimary );
185 145
186 } while( (!Success) && ++internalLoopCount 146 } while( (!Success) && ++internalLoopCounter < maxNumberOfInternalLoops );
187 147
188 if ( internalLoopCounter >= maxNumberOfInt 148 if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
189 Success = false; 149 Success = false;
190 } 150 }
191 151
192 if ( Success ) { 152 if ( Success ) {
193 #ifdef debugQGSParticipants 153 #ifdef debugQGSParticipants
194 G4cout<<G4endl<<"PerformDiffractiveCol 154 G4cout<<G4endl<<"PerformDiffractiveCollisions(); if they happend." <<G4endl;
195 #endif 155 #endif
196 156
197 PerformDiffractiveCollisions(); 157 PerformDiffractiveCollisions();
198 158
199 #ifdef debugQGSParticipants 159 #ifdef debugQGSParticipants
200 G4cout<<G4endl<<"SplitHadrons();" <<G4 160 G4cout<<G4endl<<"SplitHadrons();" <<G4endl;
201 #endif 161 #endif
202 162
203 SplitHadrons(); 163 SplitHadrons();
204 164
205 if( theProjectileSplitable && theProject 165 if( theProjectileSplitable && theProjectileSplitable->GetStatus() == 0) {
206 #ifdef debugQGSParticipants 166 #ifdef debugQGSParticipants
207 G4cout<<"Perform non-Diffractive Col 167 G4cout<<"Perform non-Diffractive Collisions if they happend. Determine Parton Momenta and so on." <<G4endl;
208 #endif 168 #endif
209 Success = DeterminePartonMomenta(); 169 Success = DeterminePartonMomenta();
210 } 170 }
211 171
212 if(!Success) PrepareInitialState( thePri 172 if(!Success) PrepareInitialState( thePrimary );
213 } 173 }
214 } while( (!Success) && ++loopCounter < maxNu 174 } while( (!Success) && ++loopCounter < maxNumberOfLoops );
215 175
216 if ( loopCounter >= maxNumberOfLoops ) { 176 if ( loopCounter >= maxNumberOfLoops ) {
217 Success = false; 177 Success = false;
218 #ifdef debugQGSParticipants 178 #ifdef debugQGSParticipants
219 G4cout<<"NOT Successful ======" <<G4endl 179 G4cout<<"NOT Successful ======" <<G4endl;
220 #endif 180 #endif
221 } 181 }
222 182
223 if ( Success ) { 183 if ( Success ) {
224 CreateStrings(); // To create strings 184 CreateStrings(); // To create strings
225 185
226 GetResiduals(); // To calculate residual 186 GetResiduals(); // To calculate residual nucleus properties
227 187
228 #ifdef debugQGSParticipants 188 #ifdef debugQGSParticipants
229 G4cout<<"All O.K. ======" <<G4endl; 189 G4cout<<"All O.K. ======" <<G4endl;
230 #endif 190 #endif
231 } 191 }
232 192
233 // clean-up, if necessary 193 // clean-up, if necessary
234 #ifdef debugQGSParticipants 194 #ifdef debugQGSParticipants
235 G4cout<<"Clearing "<<G4endl; 195 G4cout<<"Clearing "<<G4endl;
236 G4cout<<"theInteractions.size() "<<theInte 196 G4cout<<"theInteractions.size() "<<theInteractions.size()<<G4endl;
237 #endif 197 #endif
238 198
239 if( Regge ) delete Regge; 199 if( Regge ) delete Regge;
240 200
241 std::for_each( theInteractions.begin(), theI 201 std::for_each( theInteractions.begin(), theInteractions.end(), DeleteInteractionContent() );
242 theInteractions.clear(); 202 theInteractions.clear();
243 203
244 // Erasing of target involved nucleons. 204 // Erasing of target involved nucleons.
245 #ifdef debugQGSParticipants 205 #ifdef debugQGSParticipants
246 G4cout<<"Erasing of involved target nucleo 206 G4cout<<"Erasing of involved target nucleons "<<NumberOfInvolvedNucleonsOfTarget<<G4endl;
247 #endif 207 #endif
248 208
249 if ( NumberOfInvolvedNucleonsOfTarget != 0 ) 209 if ( NumberOfInvolvedNucleonsOfTarget != 0 ) {
250 for ( G4int i = 0; i < NumberOfInvolvedNu 210 for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
251 G4VSplitableHadron* aNucleon = TheInvolv 211 G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfTarget[i]->GetSplitableHadron();
252 if ( (aNucleon != 0 ) && (aNucleon->GetS 212 if ( (aNucleon != 0 ) && (aNucleon->GetStatus() >= 1) ) delete aNucleon;
253 } 213 }
254 } 214 }
255 215
256 // Erasing of projectile involved nucleons. 216 // Erasing of projectile involved nucleons.
257 if ( NumberOfInvolvedNucleonsOfProjectile != 217 if ( NumberOfInvolvedNucleonsOfProjectile != 0 ) {
258 for ( G4int i = 0; i < NumberOfInvolvedNu 218 for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfProjectile; i++ ) {
259 G4VSplitableHadron* aNucleon = TheInvol 219 G4VSplitableHadron* aNucleon = TheInvolvedNucleonsOfProjectile[i]->GetSplitableHadron();
260 if ( aNucleon ) delete aNucleon; 220 if ( aNucleon ) delete aNucleon;
261 } 221 }
262 } 222 }
263 223
264 #ifdef debugQGSParticipants 224 #ifdef debugQGSParticipants
265 G4cout<<"Delition of target nucleons from 225 G4cout<<"Delition of target nucleons from soft interactions "<<theTargets.size()
266 <<G4endl<<G4endl; 226 <<G4endl<<G4endl;
267 #endif 227 #endif
268 std::for_each(theTargets.begin(), theTargets 228 std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
269 theTargets.clear(); 229 theTargets.clear();
270 230
271 if ( theProjectileSplitable ) { 231 if ( theProjectileSplitable ) {
272 delete theProjectileSplitable; 232 delete theProjectileSplitable;
273 theProjectileSplitable = 0; 233 theProjectileSplitable = 0;
274 } 234 }
275 } 235 }
276 236
277 //============================================ 237 //===========================================================
278 void G4QGSParticipants::PrepareInitialState( c 238 void G4QGSParticipants::PrepareInitialState( const G4ReactionProduct& thePrimary )
279 { 239 {
280 // Clearing of the arrays 240 // Clearing of the arrays
281 // Erasing of the projectile 241 // Erasing of the projectile
282 G4InteractionContent* anIniteraction = theIn 242 G4InteractionContent* anIniteraction = theInteractions[0];
283 G4VSplitableHadron* pProjectile = anIniterac 243 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
284 if( pProjectile ) delete pProjectile; 244 if( pProjectile ) delete pProjectile;
285 245
286 std::for_each(theInteractions.begin(), theIn 246 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
287 theInteractions.clear(); 247 theInteractions.clear();
288 248
289 // Erasing of the envolved nucleons and targ 249 // Erasing of the envolved nucleons and target nucleons from diffraction dissociations
290 theNucleus->StartLoop(); 250 theNucleus->StartLoop();
291 G4Nucleon* aNucleon; 251 G4Nucleon* aNucleon;
292 while ( ( aNucleon = theNucleus->GetNextNucl 252 while ( ( aNucleon = theNucleus->GetNextNucleon() ) )
293 { 253 {
294 if ( aNucleon->AreYouHit() ) { 254 if ( aNucleon->AreYouHit() ) {
295 G4VSplitableHadron* splaNucleon = aNucle 255 G4VSplitableHadron* splaNucleon = aNucleon->GetSplitableHadron();
296 if ( (splaNucleon != 0) && (splaNucleon- 256 if ( (splaNucleon != 0) && (splaNucleon->GetStatus() >=1) ) delete splaNucleon;
297 aNucleon->Hit(nullptr); 257 aNucleon->Hit(nullptr);
298 NumberOfInvolvedNucleonsOfTarget--; 258 NumberOfInvolvedNucleonsOfTarget--;
299 } 259 }
300 } 260 }
301 261
302 // Erasing of nuclear nucleons participated 262 // Erasing of nuclear nucleons participated in soft interactions
303 std::for_each(theTargets.begin(), theTargets 263 std::for_each(theTargets.begin(), theTargets.end(), DeleteSplitableHadron());
304 theTargets.clear(); 264 theTargets.clear();
305 265
306 // Preparation to a new attempt 266 // Preparation to a new attempt
307 theProjectile = thePrimary; 267 theProjectile = thePrimary;
308 268
309 theNucleus->Init(theNucleus->GetMassNumber() 269 theNucleus->Init(theNucleus->GetMassNumber(), theNucleus->GetCharge());
310 theNucleus->SortNucleonsIncZ(); 270 theNucleus->SortNucleonsIncZ();
311 DoLorentzBoost(-theCurrentVelocity); // Lor 271 DoLorentzBoost(-theCurrentVelocity); // Lorentz boost of the target nucleus
312 272
313 if (theNucleus->GetMassNumber() == 1) 273 if (theNucleus->GetMassNumber() == 1)
314 { 274 {
315 G4ThreeVector aPos = G4ThreeVector(0.,0.,0 275 G4ThreeVector aPos = G4ThreeVector(0.,0.,0.);
316 theNucleus->StartLoop(); 276 theNucleus->StartLoop();
317 G4Nucleon* tNucleon=theNucleus->GetNextNuc 277 G4Nucleon* tNucleon=theNucleus->GetNextNucleon();
318 tNucleon->SetPosition(aPos); 278 tNucleon->SetPosition(aPos);
319 } 279 }
320 280
321 G4LorentzVector Tmp( 0.0, 0.0, 0.0, 0.0 ); 281 G4LorentzVector Tmp( 0.0, 0.0, 0.0, 0.0 );
322 NumberOfInvolvedNucleonsOfTarget= 0; 282 NumberOfInvolvedNucleonsOfTarget= 0;
323 TargetResidualMassNumber = theNucleus- 283 TargetResidualMassNumber = theNucleus->GetMassNumber();
324 TargetResidualCharge = theNucleus- 284 TargetResidualCharge = theNucleus->GetCharge();
325 TargetResidualExcitationEnergy = 0.0; 285 TargetResidualExcitationEnergy = 0.0;
326 286
327 G4Nucleon* NuclearNucleon; 287 G4Nucleon* NuclearNucleon;
328 while ( ( NuclearNucleon = theNucleus->GetNe 288 while ( ( NuclearNucleon = theNucleus->GetNextNucleon() ) )
329 {Tmp+=NuclearNucleon->Get4Moment 289 {Tmp+=NuclearNucleon->Get4Momentum();}
330 290
331 TargetResidual4Momentum = Tmp; 291 TargetResidual4Momentum = Tmp;
332 } 292 }
333 293
334 //============================================ 294 //===========================================================
335 void G4QGSParticipants::GetList( const G4React 295 void G4QGSParticipants::GetList( const G4ReactionProduct& thePrimary ) {
336 #ifdef debugQGSParticipants 296 #ifdef debugQGSParticipants
337 G4cout<<G4endl<<"G4QGSParticipants::GetLis 297 G4cout<<G4endl<<"G4QGSParticipants::GetList +++++++++++++"<<G4endl;
338 #endif 298 #endif
339 299
340 // Direction: True - Proj, False - Target 300 // Direction: True - Proj, False - Target
341 theProjectileSplitable = new G4QGSMSplitable 301 theProjectileSplitable = new G4QGSMSplitableHadron(thePrimary, TRUE);
342 theProjectileSplitable->SetStatus(1); 302 theProjectileSplitable->SetStatus(1);
343 303
344 G4LorentzVector aPrimaryMomentum(thePrimary. 304 G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
345 G4LorentzVector aNucleonMomentum(0.,0.,0., 9 305 G4LorentzVector aNucleonMomentum(0.,0.,0., 938.0*MeV);
346 306
347 G4double SS=(aPrimaryMomentum + aNucleonMome 307 G4double SS=(aPrimaryMomentum + aNucleonMomentum).mag2();
348 308
349 Regge->SetS(SS); 309 Regge->SetS(SS);
350 310
351 //-------------------------------------- 311 //--------------------------------------
352 theNucleus->StartLoop(); 312 theNucleus->StartLoop();
353 G4Nucleon * tNucleon = theNucleus->GetNextNu 313 G4Nucleon * tNucleon = theNucleus->GetNextNucleon();
354 314
355 if ( ! tNucleon ) { 315 if ( ! tNucleon ) {
356 #ifdef debugQGSParticipants 316 #ifdef debugQGSParticipants
357 G4cout << "QGSM - BAD situation: pNucleo 317 G4cout << "QGSM - BAD situation: pNucleon is NULL ! Leaving immediately!" << G4endl;
358 #endif 318 #endif
359 return; 319 return;
360 } 320 }
361 321
362 G4double theNucleusOuterR = theNucleus->GetO 322 G4double theNucleusOuterR = theNucleus->GetOuterRadius();
363 323
364 if (theNucleus->GetMassNumber() == 1) 324 if (theNucleus->GetMassNumber() == 1)
365 { 325 {
366 G4ThreeVector aPos = G4ThreeVector(0.,0.,0 326 G4ThreeVector aPos = G4ThreeVector(0.,0.,0.);
367 tNucleon->SetPosition(aPos); 327 tNucleon->SetPosition(aPos);
368 theNucleusOuterR = 0.; 328 theNucleusOuterR = 0.;
369 } 329 }
370 330
371 // Determination of participating nucleons o 331 // Determination of participating nucleons of nucleus ------------------------------------
372 332
373 std::for_each(theInteractions.begin(), theIn 333 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
374 theInteractions.clear(); 334 theInteractions.clear();
375 335
>> 336 G4int totalCuts = 0;
376 G4int MaxPower=thePrimary.GetMomentum().mag( 337 G4int MaxPower=thePrimary.GetMomentum().mag()/(3.3*GeV); if(MaxPower < 1) MaxPower=1;
377 338
378 const G4int maxNumberOfLoops = 1000; 339 const G4int maxNumberOfLoops = 1000;
379 340
380 G4int NumberOfTries = 0; 341 G4int NumberOfTries = 0;
381 G4double Scale = 1.0; 342 G4double Scale = 1.0;
382 343
383 G4int loopCounter = -1; 344 G4int loopCounter = -1;
384 while( (theInteractions.size() == 0) && ++lo 345 while( (theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops )
385 { 346 {
386 InteractionMode = ALL; // Mode = ALL, WIT 347 InteractionMode = ALL; // Mode = ALL, WITHOUT_R, NON_DIFF
387 348
388 // choose random impact parameter of a col 349 // choose random impact parameter of a collision
389 std::pair<G4double, G4double> theImpactPar 350 std::pair<G4double, G4double> theImpactParameter;
390 351
391 NumberOfTries++; 352 NumberOfTries++;
392 if( NumberOfTries == 100*(NumberOfTries/10 353 if( NumberOfTries == 100*(NumberOfTries/100) ) Scale /=2.0;
393 354
394 theImpactParameter = theNucleus->ChooseImp 355 theImpactParameter = theNucleus->ChooseImpactXandY(theNucleusOuterR/Scale + theNucleonRadius);
395 G4double impactX = theImpactParameter.firs 356 G4double impactX = theImpactParameter.first;
396 G4double impactY = theImpactParameter.seco 357 G4double impactY = theImpactParameter.second;
397 358
398 #ifdef debugQGSParticipants 359 #ifdef debugQGSParticipants
399 G4cout<<"InteractionMode "<<InteractionM 360 G4cout<<"InteractionMode "<<InteractionMode<<G4endl;
400 G4cout<<"Impact parameter (fm ) "<<std:: 361 G4cout<<"Impact parameter (fm ) "<<std::sqrt(sqr(impactX)+sqr(impactY))/fermi<<" "<<G4endl;
401 G4int nucleonCount = -1; <<
402 #endif 362 #endif
403 363
404 // loop over nucleons to find collisions 364 // loop over nucleons to find collisions
405 theNucleus->StartLoop(); 365 theNucleus->StartLoop();
>> 366 G4int nucleonCount = -1;
406 G4QGSParticipants_NPart = 0; 367 G4QGSParticipants_NPart = 0;
407 368
408 G4double Power=MaxPower; 369 G4double Power=MaxPower;
409 370
410 while( (tNucleon = theNucleus->GetNextNucl 371 while( (tNucleon = theNucleus->GetNextNucleon()) )
411 { 372 {
412 if(Power <= 0.) break; 373 if(Power <= 0.) break;
>> 374 nucleonCount++;
413 375
414 G4LorentzVector nucleonMomentum=tNucleon 376 G4LorentzVector nucleonMomentum=tNucleon->Get4Momentum();
415 377
416 G4double Distance2 = sqr(impactX - tNucl 378 G4double Distance2 = sqr(impactX - tNucleon->GetPosition().x()) +
417 sqr(impactY - tNucl 379 sqr(impactY - tNucleon->GetPosition().y());
418 380
419 G4double Pint(0.); // 381 G4double Pint(0.); // A probability of interaction at given impact parameter
420 G4double Pprd(0.), Ptrd(0.), Pdd(0.); // 382 G4double Pprd(0.), Ptrd(0.), Pdd(0.); // Probabilities of Proj. diffr., Target diffr., Double diffr.
421 G4double Pnd (0.), Pnvr(0.); // 383 G4double Pnd (0.), Pnvr(0.); // Probabilities of non-diffr. and quark exchange
422 G4int NcutPomerons(0); // 384 G4int NcutPomerons(0); // Number of cutted pomerons
423 385
424 Regge->GetProbabilities(std::sqrt(Distan 386 Regge->GetProbabilities(std::sqrt(Distance2), InteractionMode,
425 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr); 387 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr);
426 #ifdef debugQGSParticipants 388 #ifdef debugQGSParticipants
427 nucleonCount++; <<
428 G4cout<<"Nucleon & its impact paramete 389 G4cout<<"Nucleon & its impact parameter: "<<nucleonCount<<" "<<std::sqrt(Distance2)/fermi<<" (fm)"<<G4endl;
429 G4cout<<"Probability of interaction: 390 G4cout<<"Probability of interaction: "<<Pint<<G4endl;
430 G4cout<<"Probability of PrD, TrD, DD: "<< 391 G4cout<<"Probability of PrD, TrD, DD: "<<Pprd<<" "<<Ptrd<<" "<<Pdd<<G4endl;
431 G4cout<<"Probability of NonDiff, QuarkExc.: 392 G4cout<<"Probability of NonDiff, QuarkExc.: "<<Pnd<<" "<<Pnvr<<" in inel. inter."<<G4endl;
432 #endif 393 #endif
433 394
434 if (Pint > G4UniformRand()) 395 if (Pint > G4UniformRand())
435 { // An inte 396 { // An interaction is happend.
436 397
437 G4double rndNumber = G4UniformRand(); 398 G4double rndNumber = G4UniformRand();
438 G4int InteractionType(0); 399 G4int InteractionType(0);
439 400
440 if((InteractionMode==ALL)||(Interactio 401 if((InteractionMode==ALL)||(InteractionMode==WITHOUT_R)) // Mode = ALL, WITHOUT_R, NON_DIFF
441 { 402 {
442 if( rndNumber < Pprd ) { 403 if( rndNumber < Pprd ) {InteractionType = PrD; InteractionMode = WITHOUT_R;}
443 else if( rndNumber < Pprd+Ptrd) { 404 else if( rndNumber < Pprd+Ptrd) {InteractionType = TrD; InteractionMode = WITHOUT_R;}
444 else if( rndNumber < Pprd+Ptrd+Pdd) { 405 else if( rndNumber < Pprd+Ptrd+Pdd) {InteractionType = DD; InteractionMode = WITHOUT_R;}
445 else if( rndNumber < Pprd+Ptrd+Pdd+Pnd ) { 406 else if( rndNumber < Pprd+Ptrd+Pdd+Pnd ) {InteractionType = NonD; InteractionMode = NON_DIFF;
446 NcutPomerons = Regge->ncPom 407 NcutPomerons = Regge->ncPomerons(); }
447 else {Intera 408 else {InteractionType = Qexc; InteractionMode = ALL; }
448 } 409 }
449 else // InteractionMode == NON_DIFF 410 else // InteractionMode == NON_DIFF
450 { 411 {
451 InteractionMode = NON_DIFF; 412 InteractionMode = NON_DIFF;
452 if( rndNumber < Ptrd ) {Interact 413 if( rndNumber < Ptrd ) {InteractionType = TrD; }
453 else if( rndNumber < Ptrd + Pnd) {Interact 414 else if( rndNumber < Ptrd + Pnd) {InteractionType = NonD; NcutPomerons = Regge->ncPomerons();}
454 } 415 }
455 416
456 if( (InteractionType == NonD) && (Ncut 417 if( (InteractionType == NonD) && (NcutPomerons == 0)) continue;
457 418
458 G4QGSParticipants_NPart ++; 419 G4QGSParticipants_NPart ++;
459 G4QGSMSplitableHadron* aTargetSPB = ne 420 G4QGSMSplitableHadron* aTargetSPB = new G4QGSMSplitableHadron(*tNucleon);
460 tNucleon->Hit(aTargetSPB); 421 tNucleon->Hit(aTargetSPB);
461 422
462 #ifdef debugQGSParticipants 423 #ifdef debugQGSParticipants
463 G4cout<<"An interaction is happend." 424 G4cout<<"An interaction is happend."<<G4endl;
464 G4cout<<"Target nucleon - "<<nucleon 425 G4cout<<"Target nucleon - "<<nucleonCount<<" "
465 <<tNucleon->GetDefinition()->G 426 <<tNucleon->GetDefinition()->GetParticleName()<<G4endl;
466 G4cout<<"Interaction type:"<<Interac 427 G4cout<<"Interaction type:"<<InteractionType
467 <<" (0 -PrD, 1 - TrD, 2 - DD, 428 <<" (0 -PrD, 1 - TrD, 2 - DD, 3 - NonD, 4 - Qexc)"<<G4endl;
468 G4cout<<"New Inter. mode:"<<Interac 429 G4cout<<"New Inter. mode:"<<InteractionMode
469 <<" (0 -ALL, 1 - WITHOUT_R, 2 430 <<" (0 -ALL, 1 - WITHOUT_R, 2 - NON_DIFF)"<<G4endl;
470 if( InteractionType == NonD ) 431 if( InteractionType == NonD )
471 G4cout<<"Number of cutted pomerons 432 G4cout<<"Number of cutted pomerons: "<<NcutPomerons<<G4endl;
472 #endif 433 #endif
473 434
474 if((InteractionType == PrD) || (Intera 435 if((InteractionType == PrD) || (InteractionType == TrD) || (InteractionType == DD) ||
475 (InteractionType == Qexc)) 436 (InteractionType == Qexc))
476 { // 437 { // diffractive-like interaction occurs
477 #ifdef debugQGSParticipants 438 #ifdef debugQGSParticipants
478 G4cout<<"Diffractive-like interact 439 G4cout<<"Diffractive-like interaction occurs"<<G4endl;
479 #endif 440 #endif
480 441
481 G4InteractionContent * aInteraction 442 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
482 theProjectileSplitable->SetStatus(1* 443 theProjectileSplitable->SetStatus(1*theProjectileSplitable->GetStatus());
483 444
484 aInteraction->SetTarget(aTargetSPB); 445 aInteraction->SetTarget(aTargetSPB);
485 aInteraction->SetTargetNucleon(tNucl 446 aInteraction->SetTargetNucleon(tNucleon);
486 aTargetSPB->SetCollisionCount(0); 447 aTargetSPB->SetCollisionCount(0);
487 aTargetSPB->SetStatus(1); 448 aTargetSPB->SetStatus(1);
488 449
489 aInteraction->SetNumberOfDiffractive 450 aInteraction->SetNumberOfDiffractiveCollisions(1);
490 aInteraction->SetNumberOfSoftCollisi 451 aInteraction->SetNumberOfSoftCollisions(0);
491 aInteraction->SetStatus(InteractionT 452 aInteraction->SetStatus(InteractionType);
492 theInteractions.push_back(aInteracti 453 theInteractions.push_back(aInteraction);
493 } 454 }
494 else 455 else
495 { // non 456 { // nondiffractive interaction occurs
496 #ifdef debugQGSParticipants 457 #ifdef debugQGSParticipants
497 G4cout<<"Non-diffractive interacti 458 G4cout<<"Non-diffractive interaction occurs, max NcutPomerons "<<NcutPomerons<<G4endl;
498 #endif 459 #endif
499 460
500 G4int nCuts; 461 G4int nCuts;
501 462
502 G4int Vncut=0; 463 G4int Vncut=0;
503 for(nCuts = 0; nCuts < NcutPomerons; nCuts 464 for(nCuts = 0; nCuts < NcutPomerons; nCuts++)
504 { 465 {
505 if( G4UniformRand() < Power/MaxPower ){V 466 if( G4UniformRand() < Power/MaxPower ){Vncut++; Power--; if(Power <= 0.) break;}
506 } 467 }
507 nCuts=Vncut; 468 nCuts=Vncut;
508 469
509 if( nCuts == 0 ) {delete aTargetSPB; tNucl 470 if( nCuts == 0 ) {delete aTargetSPB; tNucleon->Hit(nullptr); continue;}
510 471
>> 472 totalCuts += nCuts;
511 #ifdef debugQGSParticipants 473 #ifdef debugQGSParticipants
512 G4cout<<"Number of cuts in the int 474 G4cout<<"Number of cuts in the interaction "<<nCuts<<G4endl;
513 #endif 475 #endif
514 476
515 aTargetSPB->IncrementCollisionCount(nCuts) 477 aTargetSPB->IncrementCollisionCount(nCuts);
516 aTargetSPB->SetStatus(0); 478 aTargetSPB->SetStatus(0);
517 theTargets.push_back(aTargetSPB); 479 theTargets.push_back(aTargetSPB);
518 480
519 theProjectileSplitable->IncrementCollision 481 theProjectileSplitable->IncrementCollisionCount(nCuts);
520 theProjectileSplitable->SetStatus(0* 482 theProjectileSplitable->SetStatus(0*theProjectileSplitable->GetStatus());
521 483
522 G4InteractionContent * aInteraction = new 484 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
523 aInteraction->SetTarget(aTargetSPB); 485 aInteraction->SetTarget(aTargetSPB);
524 aInteraction->SetTargetNucleon(tNucl 486 aInteraction->SetTargetNucleon(tNucleon);
525 aInteraction->SetNumberOfSoftCollisions(nC 487 aInteraction->SetNumberOfSoftCollisions(nCuts);
526 aInteraction->SetStatus(InteractionT 488 aInteraction->SetStatus(InteractionType);
527 theInteractions.push_back(aInteraction); 489 theInteractions.push_back(aInteraction);
528 } 490 }
529 } // End of if (Pint > G4UniformRand( 491 } // End of if (Pint > G4UniformRand())
530 } // End of while( (tNucleon = theNucl 492 } // End of while( (tNucleon = theNucleus->GetNextNucleon()) )
531 493
532 #ifdef debugQGSParticipants 494 #ifdef debugQGSParticipants
533 G4cout << G4endl<<"Number of wounded nuc 495 G4cout << G4endl<<"Number of wounded nucleons "<<G4QGSParticipants_NPart<<G4endl;
534 #endif 496 #endif
535 497
536 } // End of while( (theInteractions.size() 498 } // End of while( (theInteractions.size() == 0) && ++loopCounter < maxNumberOfLoops )
537 499
538 if ( loopCounter >= maxNumberOfLoops ) { 500 if ( loopCounter >= maxNumberOfLoops ) {
539 #ifdef debugQGSParticipants 501 #ifdef debugQGSParticipants
540 G4cout <<"BAD situation: forced loop exi 502 G4cout <<"BAD situation: forced loop exit!" << G4endl;
541 #endif 503 #endif
542 // Perhaps there is something to set here. 504 // Perhaps there is something to set here...
543 // Decrease impact parameter ?? 505 // Decrease impact parameter ??
544 // Select collisions with only diffraction 506 // Select collisions with only diffraction ??
545 // Selecy only non-diffractive interaction 507 // Selecy only non-diffractive interactions ??
546 } 508 }
547 //------------------------------------------ 509 //------------------------------------------------------------
548 std::vector<G4InteractionContent*>::iterator 510 std::vector<G4InteractionContent*>::iterator i;
549 511
550 if( theInteractions.size() != 0) 512 if( theInteractions.size() != 0)
551 { 513 {
552 if( InteractionMode == ALL ) // It can be 514 if( InteractionMode == ALL ) // It can be if all interactions were quark-exchange.
553 { // Only the 515 { // Only the first one will be saved, all other will be erased.
554 i = theInteractions.end()-1; 516 i = theInteractions.end()-1;
555 517
556 while ( theInteractions.size() != 1 ) 518 while ( theInteractions.size() != 1 )
557 { 519 {
558 G4InteractionContent* anInteraction = 520 G4InteractionContent* anInteraction = *i;
559 G4Nucleon * pNucleon = anInteraction-> 521 G4Nucleon * pNucleon = anInteraction->GetTargetNucleon(); pNucleon->Hit(nullptr);
560 delete anInteraction->GetTarget(); 522 delete anInteraction->GetTarget();
561 delete *i; 523 delete *i;
562 i=theInteractions.erase(i); 524 i=theInteractions.erase(i);
563 i--; 525 i--;
564 } 526 }
565 } 527 }
566 else 528 else
567 { // All quark 529 { // All quark exchanges will be erased
568 i = theInteractions.begin(); 530 i = theInteractions.begin();
569 while ( i != theInteractions.end() ) 531 while ( i != theInteractions.end() )
570 { 532 {
571 G4InteractionContent* anInteraction = 533 G4InteractionContent* anInteraction = *i;
572 534
573 if( anInteraction->GetStatus() == Qexc 535 if( anInteraction->GetStatus() == Qexc )
574 { 536 {
575 G4Nucleon* aTargetNucleon = a 537 G4Nucleon* aTargetNucleon = anInteraction->GetTargetNucleon();
576 aTargetNucleon->Hit(nullptr); 538 aTargetNucleon->Hit(nullptr);
577 539
578 delete anInteraction->GetTarget(); 540 delete anInteraction->GetTarget();
579 delete *i; 541 delete *i;
580 i=theInteractions.erase(i); 542 i=theInteractions.erase(i);
581 } 543 }
582 else 544 else
583 { 545 {
584 i++; 546 i++;
585 } 547 }
586 } 548 }
587 } 549 }
>> 550
>> 551 #ifdef debugQGSParticipants
>> 552 G4cout <<"Total number of cuts "<< totalCuts <<G4endl;
>> 553 #endif
588 } 554 }
589 } 555 }
590 556
591 //============================================ 557 //=============================================================
592 void G4QGSParticipants::StoreInvolvedNucleon() 558 void G4QGSParticipants::StoreInvolvedNucleon()
593 { //To store nucleons involved in the interact 559 { //To store nucleons involved in the interaction
594 560
595 NumberOfInvolvedNucleonsOfTarget = 0; 561 NumberOfInvolvedNucleonsOfTarget = 0;
596 562
597 theNucleus->StartLoop(); 563 theNucleus->StartLoop();
598 564
599 G4Nucleon* aNucleon; 565 G4Nucleon* aNucleon;
600 while ( ( aNucleon = theNucleus->GetNextNucl 566 while ( ( aNucleon = theNucleus->GetNextNucleon() ) ) {
601 if ( aNucleon->AreYouHit() ) { 567 if ( aNucleon->AreYouHit() ) {
602 TheInvolvedNucleonsOfTarget[NumberOfInvo 568 TheInvolvedNucleonsOfTarget[NumberOfInvolvedNucleonsOfTarget] = aNucleon;
603 NumberOfInvolvedNucleonsOfTarget++; 569 NumberOfInvolvedNucleonsOfTarget++;
604 } 570 }
605 } 571 }
606 572
607 #ifdef debugQGSParticipants 573 #ifdef debugQGSParticipants
608 G4cout << G4endl<<"G4QGSParticipants::Stor 574 G4cout << G4endl<<"G4QGSParticipants::StoreInvolvedNucleon() if they were "<<G4endl
609 <<"Stored # of wounded nucleons of 575 <<"Stored # of wounded nucleons of target "
610 << NumberOfInvolvedNucleonsOfTarget 576 << NumberOfInvolvedNucleonsOfTarget <<G4endl;
611 #endif 577 #endif
612 return; 578 return;
613 } 579 }
614 580
615 //============================================ 581 //=============================================================
616 582
617 void G4QGSParticipants::ReggeonCascade() 583 void G4QGSParticipants::ReggeonCascade()
618 { // Implementation of the reggeon theory insp 584 { // Implementation of the reggeon theory inspired model of nuclear destruction
619 #ifdef debugQGSParticipants 585 #ifdef debugQGSParticipants
620 G4cout << G4endl<<"Reggeon cascading ..... 586 G4cout << G4endl<<"Reggeon cascading ........."<<G4endl;
621 G4cout<<"C of nucl. desctruction "<<GetCof 587 G4cout<<"C of nucl. desctruction "<<GetCofNuclearDestruction()
622 <<" R2 "<<GetR2ofNuclearDestruction( 588 <<" R2 "<<GetR2ofNuclearDestruction()/fermi/fermi<<" fermi^2"<<G4endl;
623 #endif 589 #endif
624 590
625 G4int InitNINt = NumberOfInvolvedNucleonsOfT 591 G4int InitNINt = NumberOfInvolvedNucleonsOfTarget;
626 592
627 // Reggeon cascading in target nucleus 593 // Reggeon cascading in target nucleus
628 for ( G4int InvTN = 0; InvTN < InitNINt; Inv 594 for ( G4int InvTN = 0; InvTN < InitNINt; InvTN++ ) {
629 G4Nucleon* aTargetNucleon = TheInvolvedNuc 595 G4Nucleon* aTargetNucleon = TheInvolvedNucleonsOfTarget[ InvTN ];
630 596
631 G4double CreationTime = aTargetNucleon->Ge 597 G4double CreationTime = aTargetNucleon->GetSplitableHadron()->GetTimeOfCreation();
632 598
633 G4double XofWoundedNucleon = aTargetNucleo 599 G4double XofWoundedNucleon = aTargetNucleon->GetPosition().x();
634 G4double YofWoundedNucleon = aTargetNucleo 600 G4double YofWoundedNucleon = aTargetNucleon->GetPosition().y();
635 601
636 G4V3DNucleus* theTargetNucleus = theNucleu 602 G4V3DNucleus* theTargetNucleus = theNucleus;
637 theTargetNucleus->StartLoop(); 603 theTargetNucleus->StartLoop();
638 604
639 #ifdef debugQGSParticipants << 605 G4int TrgNuc=0;
640 G4int TrgNuc=0; <<
641 #endif <<
642 G4Nucleon* Neighbour(0); 606 G4Nucleon* Neighbour(0);
643 while ( ( Neighbour = theTargetNucleus->Ge 607 while ( ( Neighbour = theTargetNucleus->GetNextNucleon() ) ) {
644 #ifdef debugQGSParticipants << 608 TrgNuc++;
645 TrgNuc++; <<
646 #endif <<
647 if ( ! Neighbour->AreYouHit() ) { 609 if ( ! Neighbour->AreYouHit() ) {
648 G4double impact2 = sqr( XofWoundedNucl 610 G4double impact2 = sqr( XofWoundedNucleon - Neighbour->GetPosition().x() ) +
649 sqr( YofWoundedNucl 611 sqr( YofWoundedNucleon - Neighbour->GetPosition().y() );
650 612
651 if ( G4UniformRand() < GetCofNuclearDe 613 if ( G4UniformRand() < GetCofNuclearDestruction() *
652 G4Exp( -impact2 614 G4Exp( -impact2 / GetR2ofNuclearDestruction() )
653 ) { 615 ) {
654 // The neighbour nucleon is involved 616 // The neighbour nucleon is involved in the reggeon cascade
655 #ifdef debugQGSParticipants 617 #ifdef debugQGSParticipants
656 G4cout<<"Target nucleon involved i 618 G4cout<<"Target nucleon involved in reggeon cascading No "<<TrgNuc<<" "
657 <<Neighbour->GetDefinition() 619 <<Neighbour->GetDefinition()->GetParticleName()<<G4endl;
658 #endif 620 #endif
659 TheInvolvedNucleonsOfTarget[ NumberO 621 TheInvolvedNucleonsOfTarget[ NumberOfInvolvedNucleonsOfTarget ] = Neighbour;
660 NumberOfInvolvedNucleonsOfTarget++; 622 NumberOfInvolvedNucleonsOfTarget++;
661 623
662 G4QGSMSplitableHadron* targetSplitab 624 G4QGSMSplitableHadron* targetSplitable = new G4QGSMSplitableHadron( *Neighbour );
663 625
664 Neighbour->Hit( targetSplitable ); 626 Neighbour->Hit( targetSplitable );
665 targetSplitable->SetTimeOfCreation( 627 targetSplitable->SetTimeOfCreation( CreationTime );
666 targetSplitable->SetStatus( 2 ); 628 targetSplitable->SetStatus( 2 );
667 targetSplitable->SetCollisionCount(0 629 targetSplitable->SetCollisionCount(0);
668 630
669 G4InteractionContent * anInteraction 631 G4InteractionContent * anInteraction = new G4InteractionContent(theProjectileSplitable);
670 anInteraction->SetTarget(targetSplit 632 anInteraction->SetTarget(targetSplitable);
671 anInteraction->SetTargetNucleon(Neig 633 anInteraction->SetTargetNucleon(Neighbour);
672 634
673 anInteraction->SetNumberOfDiffractiv 635 anInteraction->SetNumberOfDiffractiveCollisions(1);
674 anInteraction->SetNumberOfSoftCollis 636 anInteraction->SetNumberOfSoftCollisions(0);
675 anInteraction->SetStatus(3); 637 anInteraction->SetStatus(3);
676 theInteractions.push_back(anInteract 638 theInteractions.push_back(anInteraction);
677 } 639 }
678 } 640 }
679 } 641 }
680 } 642 }
681 643
682 #ifdef debugQGSParticipants 644 #ifdef debugQGSParticipants
683 G4cout <<"Number of new involved nucleons 645 G4cout <<"Number of new involved nucleons "<<NumberOfInvolvedNucleonsOfTarget - InitNINt<<G4endl;
684 #endif 646 #endif
685 return; 647 return;
686 } 648 }
687 649
688 //============================================ 650 //============================================================================
689 651
690 G4bool G4QGSParticipants::PutOnMassShell() { 652 G4bool G4QGSParticipants::PutOnMassShell() {
691 653
692 G4bool isProjectileNucleus = false; 654 G4bool isProjectileNucleus = false;
693 if ( GetProjectileNucleus() ) { 655 if ( GetProjectileNucleus() ) {
694 isProjectileNucleus = true; 656 isProjectileNucleus = true;
695 } 657 }
696 658
697 #ifdef debugPutOnMassShell 659 #ifdef debugPutOnMassShell
698 G4cout <<G4endl<< "PutOnMassShell start .. 660 G4cout <<G4endl<< "PutOnMassShell start ..............." << G4endl;
699 if ( isProjectileNucleus ) {G4cout << "Put 661 if ( isProjectileNucleus ) {G4cout << "PutOnMassShell for Nucleus_Nucleus " << G4endl;}
700 #endif 662 #endif
701 663
702 G4LorentzVector Pprojectile( theProjectile.G 664 G4LorentzVector Pprojectile( theProjectile.GetMomentum(), theProjectile.GetTotalEnergy() );
703 if ( Pprojectile.z() < 0.0 ) { 665 if ( Pprojectile.z() < 0.0 ) {
704 return false; 666 return false;
705 } 667 }
706 668
707 G4bool isOk = true; 669 G4bool isOk = true;
708 670
709 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 671 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 );
710 G4LorentzVector PtargetResidual( 0.0, 0.0, 0 672 G4LorentzVector PtargetResidual( 0.0, 0.0, 0.0, 0.0 );
711 G4double SumMasses = 0.0; 673 G4double SumMasses = 0.0;
712 G4V3DNucleus* theTargetNucleus = GetTargetNu 674 G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
713 G4double TargetResidualMass = 0.0; 675 G4double TargetResidualMass = 0.0;
714 676
715 #ifdef debugPutOnMassShell 677 #ifdef debugPutOnMassShell
716 G4cout << "Target : "; 678 G4cout << "Target : ";
717 #endif 679 #endif
718 680
719 isOk = ComputeNucleusProperties( theTargetNu 681 isOk = ComputeNucleusProperties( theTargetNucleus, Ptarget, PtargetResidual, SumMasses,
720 TargetResid 682 TargetResidualExcitationEnergy, TargetResidualMass,
721 TargetResid 683 TargetResidualMassNumber, TargetResidualCharge );
722 684
723 if ( ! isOk ) return false; 685 if ( ! isOk ) return false;
724 686
725 G4double Mprojectile = 0.0; 687 G4double Mprojectile = 0.0;
726 G4double M2projectile = 0.0; 688 G4double M2projectile = 0.0;
727 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 ); 689 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 );
728 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0 690 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0, 0.0 );
729 G4V3DNucleus* thePrNucleus = GetProjectileNu 691 G4V3DNucleus* thePrNucleus = GetProjectileNucleus();
730 G4double PrResidualMass = 0.0; 692 G4double PrResidualMass = 0.0;
731 693
732 if ( ! isProjectileNucleus ) { // hadron-nu 694 if ( ! isProjectileNucleus ) { // hadron-nucleus collision
733 Mprojectile = Pprojectile.mag(); 695 Mprojectile = Pprojectile.mag();
734 M2projectile = Pprojectile.mag2(); 696 M2projectile = Pprojectile.mag2();
735 SumMasses += Mprojectile + 20.0*MeV; 697 SumMasses += Mprojectile + 20.0*MeV; // Maybe DM must be larger?
736 } else { // nucleus-nucleus or antinucleus- 698 } else { // nucleus-nucleus or antinucleus-nucleus collision
737 699
738 #ifdef debugPutOnMassShell 700 #ifdef debugPutOnMassShell
739 G4cout << "Projectile : "; 701 G4cout << "Projectile : ";
740 #endif 702 #endif
741 703
742 isOk = ComputeNucleusProperties( thePrNucl 704 isOk = ComputeNucleusProperties( thePrNucleus, Pproj, PprojResidual, SumMasses,
743 Projectil 705 ProjectileResidualExcitationEnergy, PrResidualMass,
744 Projectil 706 ProjectileResidualMassNumber, ProjectileResidualCharge );
745 if ( ! isOk ) return false; 707 if ( ! isOk ) return false;
746 } 708 }
747 709
748 G4LorentzVector Psum = Pprojectile + Ptarget 710 G4LorentzVector Psum = Pprojectile + Ptarget;
749 G4double SqrtS = Psum.mag(); 711 G4double SqrtS = Psum.mag();
750 G4double S = Psum.mag2(); 712 G4double S = Psum.mag2();
751 713
752 #ifdef debugPutOnMassShell 714 #ifdef debugPutOnMassShell
753 G4cout << "Pproj "<<Pprojectile<<G4endl; 715 G4cout << "Pproj "<<Pprojectile<<G4endl;
754 G4cout << "Ptarg "<<Ptarget<<G4endl; 716 G4cout << "Ptarg "<<Ptarget<<G4endl;
755 G4cout << "Psum " << Psum/GeV << " GeV" << 717 G4cout << "Psum " << Psum/GeV << " GeV" << G4endl << "SqrtS " << SqrtS/GeV << " GeV" << G4endl
756 << "SumMasses, PrResidualMass and T 718 << "SumMasses, PrResidualMass and TargetResidualMass " << SumMasses/GeV << " "
757 << PrResidualMass/GeV << " " << Tar 719 << PrResidualMass/GeV << " " << TargetResidualMass/GeV << " GeV" << G4endl;
758 G4cout << "Ptar res. "<<PtargetResidual<<G 720 G4cout << "Ptar res. "<<PtargetResidual<<G4endl;
759 #endif 721 #endif
760 722
761 if ( SqrtS < SumMasses ) { 723 if ( SqrtS < SumMasses ) {
762 return false; // It is impossible to simu 724 return false; // It is impossible to simulate after putting nuclear nucleons on mass-shell.
763 } 725 }
764 726
765 // Try to consider also the excitation energ 727 // Try to consider also the excitation energy of the residual nucleus, if this is
766 // possible, with the available energy; othe 728 // possible, with the available energy; otherwise, set the excitation energy to zero.
767 729
768 G4double savedSumMasses = SumMasses; 730 G4double savedSumMasses = SumMasses;
769 if ( isProjectileNucleus ) { 731 if ( isProjectileNucleus ) {
770 SumMasses -= std::sqrt( sqr( PrResidualMas 732 SumMasses -= std::sqrt( sqr( PrResidualMass ) + PprojResidual.perp2() );
771 SumMasses += std::sqrt( sqr( PrResidualMas 733 SumMasses += std::sqrt( sqr( PrResidualMass + ProjectileResidualExcitationEnergy )
772 + PprojResidual.pe 734 + PprojResidual.perp2() );
773 } 735 }
774 SumMasses -= std::sqrt( sqr( TargetResidualM 736 SumMasses -= std::sqrt( sqr( TargetResidualMass ) + PtargetResidual.perp2() );
775 SumMasses += std::sqrt( sqr( TargetResidualM 737 SumMasses += std::sqrt( sqr( TargetResidualMass + TargetResidualExcitationEnergy )
776 + PtargetResidual.pe 738 + PtargetResidual.perp2() );
777 739
778 if ( SqrtS < SumMasses ) { 740 if ( SqrtS < SumMasses ) {
779 SumMasses = savedSumMasses; 741 SumMasses = savedSumMasses;
780 if ( isProjectileNucleus ) { 742 if ( isProjectileNucleus ) {
781 ProjectileResidualExcitationEnergy = 0.0 743 ProjectileResidualExcitationEnergy = 0.0;
782 } 744 }
783 TargetResidualExcitationEnergy = 0.0; 745 TargetResidualExcitationEnergy = 0.0;
784 } 746 }
785 747
786 TargetResidualMass += TargetResidualExcitati 748 TargetResidualMass += TargetResidualExcitationEnergy;
787 if ( isProjectileNucleus ) { 749 if ( isProjectileNucleus ) {
788 PrResidualMass += ProjectileResidualExcita 750 PrResidualMass += ProjectileResidualExcitationEnergy;
789 } 751 }
790 752
791 #ifdef debugPutOnMassShell 753 #ifdef debugPutOnMassShell
792 if ( isProjectileNucleus ) { 754 if ( isProjectileNucleus ) {
793 G4cout << "PrResidualMass ProjResidualEx 755 G4cout << "PrResidualMass ProjResidualExcitationEnergy " << PrResidualMass/GeV << " "
794 << ProjectileResidualExcitationEnergy < 756 << ProjectileResidualExcitationEnergy << " MeV" << G4endl;
795 } 757 }
796 G4cout << "TargetResidualMass TargetResidu 758 G4cout << "TargetResidualMass TargetResidualExcitationEnergy " << TargetResidualMass/GeV << " GeV "
797 << TargetResidualExcitationEnergy < 759 << TargetResidualExcitationEnergy << " MeV" << G4endl
798 << "Sum masses " << SumMasses/GeV < 760 << "Sum masses " << SumMasses/GeV << G4endl;
799 #endif 761 #endif
800 762
801 // Sampling of nucleons what can transfer to 763 // Sampling of nucleons what can transfer to delta-isobars
802 if ( isProjectileNucleus && thePrNucleus-> 764 if ( isProjectileNucleus && thePrNucleus->GetMassNumber() != 1 ) {
803 isOk = GenerateDeltaIsobar( SqrtS, NumberO 765 isOk = GenerateDeltaIsobar( SqrtS, NumberOfInvolvedNucleonsOfProjectile,
804 TheInvolvedNuc 766 TheInvolvedNucleonsOfProjectile, SumMasses );
805 } 767 }
806 if ( theTargetNucleus->GetMassNumber() != 1 768 if ( theTargetNucleus->GetMassNumber() != 1 ) {
807 isOk = isOk && 769 isOk = isOk &&
808 GenerateDeltaIsobar( SqrtS, NumberO 770 GenerateDeltaIsobar( SqrtS, NumberOfInvolvedNucleonsOfTarget,
809 TheInvolvedNuc 771 TheInvolvedNucleonsOfTarget, SumMasses );
810 } 772 }
811 if ( ! isOk ) return false; 773 if ( ! isOk ) return false;
812 774
813 // Now we know that it is kinematically poss 775 // Now we know that it is kinematically possible to produce a final state made
814 // of the involved nucleons (or correspondin 776 // of the involved nucleons (or corresponding delta-isobars) and a residual nucleus.
815 // We have to sample the kinematical variabl 777 // We have to sample the kinematical variables which will allow to define the 4-momenta
816 // of the final state. The sampled kinematic 778 // of the final state. The sampled kinematical variables refer to the center-of-mass frame.
817 // Notice that the sampling of the transvers 779 // Notice that the sampling of the transverse momentum corresponds to take into account
818 // Fermi motion. 780 // Fermi motion.
819 781
820 // If target is nucleon - return ? 782 // If target is nucleon - return ?
821 783
822 G4LorentzRotation toCms( -1*Psum.boostVector 784 G4LorentzRotation toCms( -1*Psum.boostVector() );
823 G4LorentzVector Ptmp = toCms*Pprojectile; 785 G4LorentzVector Ptmp = toCms*Pprojectile;
824 if ( Ptmp.pz() <= 0.0 ) { // "String" movin 786 if ( Ptmp.pz() <= 0.0 ) { // "String" moving backwards in c.m.s., abort collision!
825 return false; 787 return false;
826 } 788 }
827 789
828 G4LorentzRotation toLab( toCms.inverse() ); 790 G4LorentzRotation toLab( toCms.inverse() );
829 791
830 G4double YprojectileNucleus = 0.0; 792 G4double YprojectileNucleus = 0.0;
831 if ( isProjectileNucleus ) { 793 if ( isProjectileNucleus ) {
832 Ptmp = toCms*Pproj; 794 Ptmp = toCms*Pproj;
833 YprojectileNucleus = Ptmp.rapidity(); 795 YprojectileNucleus = Ptmp.rapidity();
834 } 796 }
835 Ptmp = toCms*Ptarget; 797 Ptmp = toCms*Ptarget;
836 G4double YtargetNucleus = Ptmp.rapidity(); 798 G4double YtargetNucleus = Ptmp.rapidity();
837 799
838 // Ascribing of the involved nucleons Pt and 800 // Ascribing of the involved nucleons Pt and Xminus
839 G4double DcorP = 0.0; 801 G4double DcorP = 0.0;
840 if ( isProjectileNucleus ) { 802 if ( isProjectileNucleus ) {
841 DcorP = GetDofNuclearDestruction() / thePr 803 DcorP = GetDofNuclearDestruction() / thePrNucleus->GetMassNumber();
842 } 804 }
843 G4double DcorT = GetDofNuclearDestruct 805 G4double DcorT = GetDofNuclearDestruction() / theTargetNucleus->GetMassNumber();
844 G4double AveragePt2 = GetPt2ofNuclearDestru 806 G4double AveragePt2 = GetPt2ofNuclearDestruction();
845 G4double maxPtSquare = GetMaxPt2ofNuclearDes 807 G4double maxPtSquare = GetMaxPt2ofNuclearDestruction();
846 808
847 #ifdef debugPutOnMassShell 809 #ifdef debugPutOnMassShell
848 if ( isProjectileNucleus ) { 810 if ( isProjectileNucleus ) {
849 G4cout << "Y projectileNucleus " << Ypro 811 G4cout << "Y projectileNucleus " << YprojectileNucleus << G4endl;
850 } 812 }
851 G4cout << "Y targetNucleus " << Ytarge 813 G4cout << "Y targetNucleus " << YtargetNucleus << G4endl
852 << "Dcor " << GetDofNuclearDestruct 814 << "Dcor " << GetDofNuclearDestruction()
853 << " DcorP DcorT " << DcorP << " " 815 << " DcorP DcorT " << DcorP << " " << DcorT << " AveragePt2 " << AveragePt2 << G4endl;
854 #endif 816 #endif
855 817
856 G4double M2proj = M2projectile; // Initiali 818 G4double M2proj = M2projectile; // Initialization needed only for hadron-nucleus collisions
857 G4double WplusProjectile = 0.0; 819 G4double WplusProjectile = 0.0;
858 G4double M2target = 0.0; 820 G4double M2target = 0.0;
859 G4double WminusTarget = 0.0; 821 G4double WminusTarget = 0.0;
860 G4int NumberOfTries = 0; 822 G4int NumberOfTries = 0;
861 G4double ScaleFactor = 1.0; 823 G4double ScaleFactor = 1.0;
862 G4bool OuterSuccess = true; 824 G4bool OuterSuccess = true;
863 825
864 const G4int maxNumberOfLoops = 1000; 826 const G4int maxNumberOfLoops = 1000;
865 G4int loopCounter = 0; 827 G4int loopCounter = 0;
866 do { 828 do {
867 G4double sqrtM2proj = 0.0, sqrtM2target = 829 G4double sqrtM2proj = 0.0, sqrtM2target = 0.0;
868 OuterSuccess = true; 830 OuterSuccess = true;
869 const G4int maxNumberOfTries = 1000; 831 const G4int maxNumberOfTries = 1000;
870 do { 832 do {
871 NumberOfTries++; 833 NumberOfTries++;
872 if ( NumberOfTries == 100*(NumberOfTries 834 if ( NumberOfTries == 100*(NumberOfTries/100) ) {
873 // After many tries, it is convenient 835 // After many tries, it is convenient to reduce the values of DcorP, DcorT and
874 // AveragePt2, so that the sampled mom 836 // AveragePt2, so that the sampled momenta (respectively, pz, and pt) of the
875 // involved nucleons (or corresponding delta 837 // involved nucleons (or corresponding delta-isomers) are smaller, and therefore
876 // it is more likely to satisfy the mo 838 // it is more likely to satisfy the momentum conservation.
877 ScaleFactor /= 2.0; 839 ScaleFactor /= 2.0;
878 DcorP *= ScaleFactor; 840 DcorP *= ScaleFactor;
879 DcorT *= ScaleFactor; 841 DcorT *= ScaleFactor;
880 AveragePt2 *= ScaleFactor; 842 AveragePt2 *= ScaleFactor;
881 } 843 }
882 if ( isProjectileNucleus ) { 844 if ( isProjectileNucleus ) {
883 // Sampling of kinematical properties 845 // Sampling of kinematical properties of projectile nucleons
884 isOk = SamplingNucleonKinematics( Aver 846 isOk = SamplingNucleonKinematics( AveragePt2, maxPtSquare, DcorP,
885 theP 847 thePrNucleus, PprojResidual,
886 PrRe 848 PrResidualMass, ProjectileResidualMassNumber,
887 Numb 849 NumberOfInvolvedNucleonsOfProjectile,
888 TheI 850 TheInvolvedNucleonsOfProjectile, M2proj );
889 } 851 }
890 // Sampling of kinematical properties of 852 // Sampling of kinematical properties of target nucleons
891 isOk = isOk && 853 isOk = isOk &&
892 SamplingNucleonKinematics( Averag 854 SamplingNucleonKinematics( AveragePt2, maxPtSquare, DcorT,
893 theTar 855 theTargetNucleus, PtargetResidual,
894 Target 856 TargetResidualMass, TargetResidualMassNumber,
895 Number 857 NumberOfInvolvedNucleonsOfTarget,
896 TheInv 858 TheInvolvedNucleonsOfTarget, M2target );
897 859
898 if ( M2proj < 0.0 ) { 860 if ( M2proj < 0.0 ) {
899 if( M2proj < -0.000001 ) { << 861 G4ExceptionDescription ed;
900 G4ExceptionDescription ed; << 862 ed << "Projectile " << theProjectile.GetDefinition()->GetParticleName()
901 ed << "Projectile " << theProjectile.GetDe << 863 << " Target (Z,A)=(" << theTargetNucleus->GetCharge() << "," << theTargetNucleus->GetMassNumber()
902 << " Target (Z,A)=(" << theTargetNucle <<
903 << ") M2proj=" << M2proj << " -> 864 << ") M2proj=" << M2proj << " -> sets it to 0.0 !" << G4endl;
904 G4Exception( "G4QGSParticipants::PutOnMass << 865 G4Exception( "G4QGSParticipants::PutOnMassShell(): negative projectile squared mass!",
905 "HAD_QGSPARTICIPANTS_002", JustWarn << 866 "HAD_QGSPARTICIPANTS_002", JustWarning, ed );
906 } <<
907 M2proj = 0.0; 867 M2proj = 0.0;
908 } 868 }
909 sqrtM2proj = std::sqrt( M2proj ); 869 sqrtM2proj = std::sqrt( M2proj );
910 if ( M2target < 0.0 ) { 870 if ( M2target < 0.0 ) {
911 G4ExceptionDescription ed; 871 G4ExceptionDescription ed;
912 ed << "Projectile " << theProjectile.G 872 ed << "Projectile " << theProjectile.GetDefinition()->GetParticleName()
913 << " Target (Z,A)=(" << theTargetN 873 << " Target (Z,A)=(" << theTargetNucleus->GetCharge() << "," << theTargetNucleus->GetMassNumber()
914 << ") M2target=" << M2target << " 874 << ") M2target=" << M2target << " -> sets it to 0.0 !" << G4endl;
915 G4Exception( "G4QGSParticipants::PutOn 875 G4Exception( "G4QGSParticipants::PutOnMassShell(): negative target squared mass!",
916 "HAD_QGSPARTICIPANTS_003", 876 "HAD_QGSPARTICIPANTS_003", JustWarning, ed );
917 M2target = 0.0; 877 M2target = 0.0;
918 }; 878 };
919 sqrtM2target = std::sqrt( M2target ); 879 sqrtM2target = std::sqrt( M2target );
920 880
921 #ifdef debugPutOnMassShell 881 #ifdef debugPutOnMassShell
922 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << S 882 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << SqrtS/GeV << " "
923 << ( sqrtM2proj + sqrtM2target 883 << ( sqrtM2proj + sqrtM2target )/GeV << " "
924 << sqrtM2proj/GeV << " " << sqr 884 << sqrtM2proj/GeV << " " << sqrtM2target/GeV << G4endl;
925 #endif 885 #endif
926 886
927 if ( ! isOk ) return false; 887 if ( ! isOk ) return false;
928 } while ( ( SqrtS < ( sqrtM2proj + sqrtM2t 888 } while ( ( SqrtS < ( sqrtM2proj + sqrtM2target ) ) &&
929 ++NumberOfTries < maxNumberOfTri 889 ++NumberOfTries < maxNumberOfTries ); /* Loop checking, 07.08.2015, A.Ribon */
930 if ( NumberOfTries >= maxNumberOfTries ) { 890 if ( NumberOfTries >= maxNumberOfTries ) {
931 return false; 891 return false;
932 } 892 }
933 if ( isProjectileNucleus ) { 893 if ( isProjectileNucleus ) {
934 isOk = CheckKinematics( S, SqrtS, M2proj 894 isOk = CheckKinematics( S, SqrtS, M2proj, M2target, YprojectileNucleus, true,
935 NumberOfInvolved 895 NumberOfInvolvedNucleonsOfProjectile,
936 TheInvolvedNucle 896 TheInvolvedNucleonsOfProjectile,
937 WminusTarget, Wp 897 WminusTarget, WplusProjectile, OuterSuccess );
938 } 898 }
939 isOk = isOk && 899 isOk = isOk &&
940 CheckKinematics( S, SqrtS, M2proj, 900 CheckKinematics( S, SqrtS, M2proj, M2target, YtargetNucleus, false,
941 NumberOfInvolvedNu 901 NumberOfInvolvedNucleonsOfTarget, TheInvolvedNucleonsOfTarget,
942 WminusTarget, Wplu 902 WminusTarget, WplusProjectile, OuterSuccess );
943 if ( ! isOk ) return false; 903 if ( ! isOk ) return false;
944 } while ( ( ! OuterSuccess ) && 904 } while ( ( ! OuterSuccess ) &&
945 ++loopCounter < maxNumberOfLoops ) 905 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
946 if ( loopCounter >= maxNumberOfLoops ) { 906 if ( loopCounter >= maxNumberOfLoops ) {
947 return false; 907 return false;
948 } 908 }
949 909
950 // Now the sampling is completed, and we can 910 // Now the sampling is completed, and we can determine the kinematics of the
951 // whole system. This is done first in the c 911 // whole system. This is done first in the center-of-mass frame, and then it is boosted
952 // to the lab frame. The transverse momentum 912 // to the lab frame. The transverse momentum of the residual nucleus is determined as
953 // the recoil of each hadron (nucleon or del 913 // the recoil of each hadron (nucleon or delta) which is emitted, i.e. in such a way
954 // to conserve (by construction) the transve 914 // to conserve (by construction) the transverse momentum.
955 915
956 if ( ! isProjectileNucleus ) { // hadron-nu 916 if ( ! isProjectileNucleus ) { // hadron-nucleus collision
957 917
958 G4double Pzprojectile = WplusProjectile/2. 918 G4double Pzprojectile = WplusProjectile/2.0 - M2projectile/2.0/WplusProjectile;
959 G4double Eprojectile = WplusProjectile/2. 919 G4double Eprojectile = WplusProjectile/2.0 + M2projectile/2.0/WplusProjectile;
960 Pprojectile.setPz( Pzprojectile ); 920 Pprojectile.setPz( Pzprojectile );
961 Pprojectile.setE( Eprojectile ); 921 Pprojectile.setE( Eprojectile );
962 922
963 #ifdef debugPutOnMassShell 923 #ifdef debugPutOnMassShell
964 G4cout << "Proj after in CMS " << Pproje 924 G4cout << "Proj after in CMS " << Pprojectile/GeV <<" GeV"<< G4endl;
965 #endif 925 #endif
966 926
967 Pprojectile.transform( toLab ); 927 Pprojectile.transform( toLab );
968 theProjectile.SetMomentum( Pprojectile.vec 928 theProjectile.SetMomentum( Pprojectile.vect() );
969 theProjectile.SetTotalEnergy( Pprojectile. 929 theProjectile.SetTotalEnergy( Pprojectile.e() );
970 930
971 if ( theProjectileSplitable ) theProjectil 931 if ( theProjectileSplitable ) theProjectileSplitable->Set4Momentum(Pprojectile);
972 932
973 #ifdef debugPutOnMassShell 933 #ifdef debugPutOnMassShell
974 G4cout << "Final proj. mom in Lab. " <<t 934 G4cout << "Final proj. mom in Lab. " <<theProjectile.GetMomentum()/GeV<<" "
975 <<t 935 <<theProjectile.GetTotalEnergy()/GeV<<" GeV"<<G4endl;
976 #endif 936 #endif
977 937
978 } else { // nucleus-nucleus or antinucleus- 938 } else { // nucleus-nucleus or antinucleus-nucleus collision
979 939
980 isOk = FinalizeKinematics( WplusProjectile 940 isOk = FinalizeKinematics( WplusProjectile, true, toLab, PrResidualMass,
981 ProjectileResid 941 ProjectileResidualMassNumber, NumberOfInvolvedNucleonsOfProjectile,
982 TheInvolvedNucl 942 TheInvolvedNucleonsOfProjectile, ProjectileResidual4Momentum );
983 943
984 #ifdef debugPutOnMassShell 944 #ifdef debugPutOnMassShell
985 G4cout << "Projectile Residual4Momentum 945 G4cout << "Projectile Residual4Momentum in CMS " << ProjectileResidual4Momentum/GeV <<" GeV"<< G4endl;
986 #endif 946 #endif
987 947
988 if ( ! isOk ) return false; 948 if ( ! isOk ) return false;
989 949
990 ProjectileResidual4Momentum.transform( toL 950 ProjectileResidual4Momentum.transform( toLab );
991 951
992 #ifdef debugPutOnMassShell 952 #ifdef debugPutOnMassShell
993 G4cout << "Projectile Residual4Momentum 953 G4cout << "Projectile Residual4Momentum in Lab " << ProjectileResidual4Momentum/GeV <<" GeV"<< G4endl;
994 #endif 954 #endif
995 955
996 } 956 }
997 957
998 isOk = FinalizeKinematics( WminusTarget, fal 958 isOk = FinalizeKinematics( WminusTarget, false, toLab, TargetResidualMass,
999 TargetResidualMas 959 TargetResidualMassNumber, NumberOfInvolvedNucleonsOfTarget,
1000 TheInvolvedNucle 960 TheInvolvedNucleonsOfTarget, TargetResidual4Momentum );
1001 961
1002 #ifdef debugPutOnMassShell 962 #ifdef debugPutOnMassShell
1003 G4cout << "Target Residual4Momentum in CM 963 G4cout << "Target Residual4Momentum in CMS " << TargetResidual4Momentum/GeV << " GeV "<< G4endl;
1004 #endif 964 #endif
1005 965
1006 if ( ! isOk ) return false; 966 if ( ! isOk ) return false;
1007 967
1008 TargetResidual4Momentum.transform( toLab ); 968 TargetResidual4Momentum.transform( toLab );
1009 969
1010 #ifdef debugPutOnMassShell 970 #ifdef debugPutOnMassShell
1011 G4cout << "Target Residual4Momentum in La 971 G4cout << "Target Residual4Momentum in Lab " << TargetResidual4Momentum/GeV << " GeV "<< G4endl;
1012 #endif 972 #endif
1013 973
1014 return true; 974 return true;
1015 975
1016 } 976 }
1017 977
1018 //=========================================== 978 //============================================================================
1019 979
1020 G4ThreeVector G4QGSParticipants::GaussianPt( 980 G4ThreeVector G4QGSParticipants::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const {
1021 // @@ this method is used in FTFModel as w 981 // @@ this method is used in FTFModel as well. Should go somewhere common!
1022 982
1023 G4double Pt2( 0.0 ), Pt(0.0); << 983 G4double Pt2( 0.0 );
1024 if ( AveragePt2 > 0.0 ) { << 984 if ( AveragePt2 <= 0.0 ) {
1025 G4double x = maxPtSquare/AveragePt2; << 985 Pt2 = 0.0;
1026 Pt2 = (x < 200) ? << 986 } else {
1027 -AveragePt2 * G4Log( 1.0 + G4UniformRan << 987 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() * ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) );
1028 : -AveragePt2 * G4Log( 1.0 - G4UniformR <<
1029 Pt = std::sqrt( Pt2 ); <<
1030 } 988 }
>> 989 G4double Pt = std::sqrt( Pt2 );
1031 G4double phi = G4UniformRand() * twopi; 990 G4double phi = G4UniformRand() * twopi;
1032 991
1033 return G4ThreeVector( Pt*std::cos(phi), Pt* 992 return G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0.0 );
1034 } 993 }
1035 //=========================================== 994 //============================================================================
1036 995
1037 G4bool G4QGSParticipants:: 996 G4bool G4QGSParticipants::
1038 ComputeNucleusProperties( G4V3DNucleus* nucle 997 ComputeNucleusProperties( G4V3DNucleus* nucleus, // input parameter
1039 G4LorentzVector& nu 998 G4LorentzVector& nucleusMomentum, // input & output parameter
1040 G4LorentzVector& re 999 G4LorentzVector& residualMomentum, // input & output parameter
1041 G4double& sumMasses 1000 G4double& sumMasses, // input & output parameter
1042 G4double& residualE 1001 G4double& residualExcitationEnergy, // input & output parameter
1043 G4double& residualM 1002 G4double& residualMass, // input & output parameter
1044 G4int& residualMass 1003 G4int& residualMassNumber, // input & output parameter
1045 G4int& residualChar 1004 G4int& residualCharge ) { // input & output parameter
1046 1005
1047 // This method, which is called only by Put 1006 // This method, which is called only by PutOnMassShell, computes some nucleus properties for:
1048 // - either the target nucleus (which is n 1007 // - either the target nucleus (which is never an antinucleus): this for any kind
1049 // of hadronic interaction (hadron-nucle 1008 // of hadronic interaction (hadron-nucleus, nucleus-nucleus, antinucleus-nucleus);
1050 // - or the projectile nucleus or antinucl 1009 // - or the projectile nucleus or antinucleus: this only in the case of nucleus-nucleus
1051 // or antinucleus-nucleus interaction. 1010 // or antinucleus-nucleus interaction.
1052 // This method assumes that the all the par 1011 // This method assumes that the all the parameters have been initialized by the caller;
1053 // the action of this method consists in mo 1012 // the action of this method consists in modifying all these parameters, except the
1054 // first one. The return value is "false" o 1013 // first one. The return value is "false" only in the case the pointer to the nucleus
1055 // is null. 1014 // is null.
1056 1015
1057 if ( ! nucleus ) return false; 1016 if ( ! nucleus ) return false;
1058 1017
1059 G4double ExcitationEPerWoundedNucleon = Get 1018 G4double ExcitationEPerWoundedNucleon = GetExcitationEnergyPerWoundedNucleon();
1060 1019
1061 // Loop over the nucleons of the nucleus. 1020 // Loop over the nucleons of the nucleus.
1062 // The nucleons that have been involved in 1021 // The nucleons that have been involved in the interaction (either from Glauber or
1063 // Reggeon Cascading) will be candidate to 1022 // Reggeon Cascading) will be candidate to be emitted.
1064 // All the remaining nucleons will be the n 1023 // All the remaining nucleons will be the nucleons of the candidate residual nucleus.
1065 // The variable sumMasses is the amount of 1024 // The variable sumMasses is the amount of energy corresponding to:
1066 // 1. transverse mass of each involved 1025 // 1. transverse mass of each involved nucleon
1067 // 2. 20.0*MeV separation energy for ea 1026 // 2. 20.0*MeV separation energy for each involved nucleon
1068 // 3. transverse mass of the residual n 1027 // 3. transverse mass of the residual nucleus
1069 // In this first evaluation of sumMasses, t 1028 // In this first evaluation of sumMasses, the excitation energy of the residual nucleus
1070 // (residualExcitationEnergy, estimated by 1029 // (residualExcitationEnergy, estimated by adding a constant value to each involved
1071 // nucleon) is not taken into account. 1030 // nucleon) is not taken into account.
1072 G4Nucleon* aNucleon = 0; 1031 G4Nucleon* aNucleon = 0;
1073 nucleus->StartLoop(); 1032 nucleus->StartLoop();
1074 while ( ( aNucleon = nucleus->GetNextNucleo 1033 while ( ( aNucleon = nucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
1075 nucleusMomentum += aNucleon->Get4Momentum 1034 nucleusMomentum += aNucleon->Get4Momentum();
1076 if ( aNucleon->AreYouHit() ) { // Involv 1035 if ( aNucleon->AreYouHit() ) { // Involved nucleons
1077 // Consider in sumMasses the nominal, i 1036 // Consider in sumMasses the nominal, i.e. on-shell, masses of the nucleons
1078 // (not the current masses, which could 1037 // (not the current masses, which could be different because the nucleons are off-shell).
1079 1038
1080 sumMasses += std::sqrt( sqr( aNucleon-> 1039 sumMasses += std::sqrt( sqr( aNucleon->GetDefinition()->GetPDGMass() )
1081 + aNucleon->Ge 1040 + aNucleon->Get4Momentum().perp2() );
1082 sumMasses += 20.0*MeV; // Separation e 1041 sumMasses += 20.0*MeV; // Separation energy for a nucleon
1083 1042
1084 //residualExcitationEnergy += Excitatio 1043 //residualExcitationEnergy += ExcitationEPerWoundedNucleon;
1085 residualExcitationEnergy += -Excitation 1044 residualExcitationEnergy += -ExcitationEPerWoundedNucleon*G4Log( G4UniformRand());
1086 residualMassNumber--; 1045 residualMassNumber--;
1087 // The absolute value below is needed o 1046 // The absolute value below is needed only in the case of anti-nucleus.
1088 residualCharge -= std::abs( G4int( aNuc 1047 residualCharge -= std::abs( G4int( aNucleon->GetDefinition()->GetPDGCharge() ) );
1089 } else { // Spectator nucleons 1048 } else { // Spectator nucleons
1090 residualMomentum += aNucleon->Get4Momen 1049 residualMomentum += aNucleon->Get4Momentum();
1091 } 1050 }
1092 } 1051 }
1093 #ifdef debugPutOnMassShell 1052 #ifdef debugPutOnMassShell
1094 G4cout << "ExcitationEnergyPerWoundedNucl 1053 G4cout << "ExcitationEnergyPerWoundedNucleon " << ExcitationEPerWoundedNucleon <<" MeV"<<G4endl
1095 << "\t Residual Charge, MassNumber 1054 << "\t Residual Charge, MassNumber " << residualCharge << " " << residualMassNumber
1096 << G4endl << "\t Initial Momentum 1055 << G4endl << "\t Initial Momentum " << nucleusMomentum/GeV<<" GeV"
1097 << G4endl << "\t Residual Momentum 1056 << G4endl << "\t Residual Momentum " << residualMomentum/GeV<<" GeV"<<G4endl;
1098 #endif 1057 #endif
1099 1058
1100 residualMomentum.setPz( 0.0 ); 1059 residualMomentum.setPz( 0.0 );
1101 residualMomentum.setE( 0.0 ); 1060 residualMomentum.setE( 0.0 );
1102 if ( residualMassNumber == 0 ) { 1061 if ( residualMassNumber == 0 ) {
1103 residualMass = 0.0; 1062 residualMass = 0.0;
1104 residualExcitationEnergy = 0.0; 1063 residualExcitationEnergy = 0.0;
1105 } else { 1064 } else {
1106 residualMass = G4ParticleTable::GetPartic 1065 residualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
1107 GetIonMass( residualChar 1066 GetIonMass( residualCharge, residualMassNumber );
1108 if ( residualMassNumber == 1 ) { 1067 if ( residualMassNumber == 1 ) {
1109 residualExcitationEnergy = 0.0; 1068 residualExcitationEnergy = 0.0;
1110 } 1069 }
1111 } 1070 }
1112 sumMasses += std::sqrt( sqr( residualMass ) 1071 sumMasses += std::sqrt( sqr( residualMass ) + residualMomentum.perp2() );
1113 return true; 1072 return true;
1114 } 1073 }
1115 1074
1116 1075
1117 //=========================================== 1076 //============================================================================
1118 1077
1119 G4bool G4QGSParticipants:: 1078 G4bool G4QGSParticipants::
1120 GenerateDeltaIsobar( const G4double sqrtS, 1079 GenerateDeltaIsobar( const G4double sqrtS, // input parameter
1121 const G4int numberOfInvo 1080 const G4int numberOfInvolvedNucleons, // input parameter
1122 G4Nucleon* involvedNucle 1081 G4Nucleon* involvedNucleons[], // input & output parameter
1123 G4double& sumMasses ) { 1082 G4double& sumMasses ) { // input & output parameter
1124 1083
1125 // This method, which is called only by Put 1084 // This method, which is called only by PutOnMassShell, check whether is possible to
1126 // re-interpret some of the involved nucleo 1085 // re-interpret some of the involved nucleons as delta-isobars:
1127 // - either by replacing a proton (2212) wi 1086 // - either by replacing a proton (2212) with a Delta+ (2214),
1128 // - or by replacing a neutron (2112) with 1087 // - or by replacing a neutron (2112) with a Delta0 (2114).
1129 // The on-shell mass of these delta-isobars 1088 // The on-shell mass of these delta-isobars is ~1232 MeV, so ~292-294 MeV heavier than
1130 // the corresponding nucleon on-shell mass. 1089 // the corresponding nucleon on-shell mass. However 400.0*MeV is considered to estimate
1131 // the max number of deltas compatible with 1090 // the max number of deltas compatible with the available energy.
1132 // The delta-isobars are considered with th 1091 // The delta-isobars are considered with the same transverse momentum as their
1133 // corresponding nucleons. 1092 // corresponding nucleons.
1134 // This method assumes that all the paramet 1093 // This method assumes that all the parameters have been initialized by the caller;
1135 // the action of this method consists in mo 1094 // the action of this method consists in modifying (eventually) involveNucleons and
1136 // sumMasses. The return value is "false" o 1095 // sumMasses. The return value is "false" only in the case that the input parameters
1137 // have unphysical values. 1096 // have unphysical values.
1138 1097
1139 if ( sqrtS < 0.0 || numberOfInvolvedNucle 1098 if ( sqrtS < 0.0 || numberOfInvolvedNucleons <= 0 || sumMasses < 0.0 ) return false;
1140 1099
1141 //const G4double ProbDeltaIsobar = 0.05; / 1100 //const G4double ProbDeltaIsobar = 0.05; // Uzhi 6.07.2012
1142 //const G4double ProbDeltaIsobar = 0.25; / 1101 //const G4double ProbDeltaIsobar = 0.25; // Uzhi 13.06.2013
1143 const G4double probDeltaIsobar = 0.10; // 1102 const G4double probDeltaIsobar = 0.10; // A.R. 07.08.2013
1144 1103
1145 G4int maxNumberOfDeltas = G4int( (sqrtS - s 1104 G4int maxNumberOfDeltas = G4int( (sqrtS - sumMasses)/(400.0*MeV) );
1146 G4int numberOfDeltas = 0; 1105 G4int numberOfDeltas = 0;
1147 1106
1148 for ( G4int i = 0; i < numberOfInvolvedNucl 1107 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1149 //G4cout << "i maxNumberOfDeltas probDelt 1108 //G4cout << "i maxNumberOfDeltas probDeltaIsobar " << i << " " << maxNumberOfDeltas
1150 // << " " << probDeltaIsobar << G4e 1109 // << " " << probDeltaIsobar << G4endl;
1151 if ( G4UniformRand() < probDeltaIsobar & 1110 if ( G4UniformRand() < probDeltaIsobar && numberOfDeltas < maxNumberOfDeltas ) {
1152 numberOfDeltas++; 1111 numberOfDeltas++;
1153 if ( ! involvedNucleons[i] ) continue; 1112 if ( ! involvedNucleons[i] ) continue;
1154 G4VSplitableHadron* splitableHadron = i 1113 G4VSplitableHadron* splitableHadron = involvedNucleons[i]->GetSplitableHadron();
1155 G4double massNuc = std::sqrt( sqr( spli 1114 G4double massNuc = std::sqrt( sqr( splitableHadron->GetDefinition()->GetPDGMass() )
1156 + splitab 1115 + splitableHadron->Get4Momentum().perp2() );
1157 //AR The absolute value below is needed 1116 //AR The absolute value below is needed in the case of an antinucleus.
1158 G4int pdgCode = std::abs( splitableHadr 1117 G4int pdgCode = std::abs( splitableHadron->GetDefinition()->GetPDGEncoding() );
1159 const G4ParticleDefinition* old_def = s 1118 const G4ParticleDefinition* old_def = splitableHadron->GetDefinition();
1160 G4int newPdgCode = pdgCode/10; newPdgCo 1119 G4int newPdgCode = pdgCode/10; newPdgCode = newPdgCode*10 + 4; // Delta
1161 if ( splitableHadron->GetDefinition()-> 1120 if ( splitableHadron->GetDefinition()->GetPDGEncoding() < 0 ) newPdgCode *= -1;
1162 const G4ParticleDefinition* ptr = 1121 const G4ParticleDefinition* ptr =
1163 G4ParticleTable::GetParticleTable()-> 1122 G4ParticleTable::GetParticleTable()->FindParticle( newPdgCode );
1164 splitableHadron->SetDefinition( ptr ); 1123 splitableHadron->SetDefinition( ptr );
1165 G4double massDelta = std::sqrt( sqr( sp 1124 G4double massDelta = std::sqrt( sqr( splitableHadron->GetDefinition()->GetPDGMass() )
1166 + split 1125 + splitableHadron->Get4Momentum().perp2() );
1167 //G4cout << i << " " << sqrtS/GeV << " 1126 //G4cout << i << " " << sqrtS/GeV << " " << sumMasses/GeV << " " << massDelta/GeV
1168 // << " " << massNuc << G4endl; 1127 // << " " << massNuc << G4endl;
1169 if ( sqrtS < sumMasses + massDelta - ma 1128 if ( sqrtS < sumMasses + massDelta - massNuc ) { // Change cannot be accepted!
1170 splitableHadron->SetDefinition( old_d 1129 splitableHadron->SetDefinition( old_def );
1171 break; 1130 break;
1172 } else { // Change is accepted 1131 } else { // Change is accepted
1173 sumMasses += ( massDelta - massNuc ); 1132 sumMasses += ( massDelta - massNuc );
1174 } 1133 }
1175 } 1134 }
1176 } 1135 }
1177 //G4cout << "maxNumberOfDeltas numberOfDelt 1136 //G4cout << "maxNumberOfDeltas numberOfDeltas " << maxNumberOfDeltas << " "
1178 // << numberOfDeltas << G4endl; 1137 // << numberOfDeltas << G4endl;
1179 return true; 1138 return true;
1180 } 1139 }
1181 1140
1182 1141
1183 //=========================================== 1142 //============================================================================
1184 1143
1185 G4bool G4QGSParticipants:: 1144 G4bool G4QGSParticipants::
1186 SamplingNucleonKinematics( G4double averagePt 1145 SamplingNucleonKinematics( G4double averagePt2, // input parameter
1187 const G4double max 1146 const G4double maxPt2, // input parameter
1188 G4double dCor, 1147 G4double dCor, // input parameter
1189 G4V3DNucleus* nucl 1148 G4V3DNucleus* nucleus, // input parameter
1190 const G4LorentzVec 1149 const G4LorentzVector& pResidual, // input parameter
1191 const G4double res 1150 const G4double residualMass, // input parameter
1192 const G4int residu 1151 const G4int residualMassNumber, // input parameter
1193 const G4int number 1152 const G4int numberOfInvolvedNucleons, // input parameter
1194 G4Nucleon* involve 1153 G4Nucleon* involvedNucleons[], // input & output parameter
1195 G4double& mass2 ) 1154 G4double& mass2 ) { // output parameter
1196 1155
1197 // This method, which is called only by Put 1156 // This method, which is called only by PutOnMassShell, does the sampling of:
1198 // - either the target nucleons: this for 1157 // - either the target nucleons: this for any kind of hadronic interactions
1199 // (hadron-nucleus, nucleus-nucleus, ant 1158 // (hadron-nucleus, nucleus-nucleus, antinucleus-nucleus);
1200 // - or the projectile nucleons or antinuc 1159 // - or the projectile nucleons or antinucleons: this only in the case of
1201 // nucleus-nucleus or antinucleus-nucleu 1160 // nucleus-nucleus or antinucleus-nucleus interactions, respectively.
1202 // This method assumes that all the paramet 1161 // This method assumes that all the parameters have been initialized by the caller;
1203 // the action of this method consists in ch 1162 // the action of this method consists in changing the properties of the nucleons
1204 // whose pointers are in the vector involve 1163 // whose pointers are in the vector involvedNucleons, as well as changing the
1205 // variable mass2. 1164 // variable mass2.
1206 1165
1207 if ( ! nucleus ) return false; 1166 if ( ! nucleus ) return false;
1208 1167
1209 if ( residualMassNumber == 0 && numberOfI 1168 if ( residualMassNumber == 0 && numberOfInvolvedNucleons == 1 ) {
1210 dCor = 0.0; 1169 dCor = 0.0;
1211 averagePt2 = 0.0; 1170 averagePt2 = 0.0;
1212 } 1171 }
1213 1172
1214 G4bool success = true; 1173 G4bool success = true;
1215 1174
1216 G4double SumMasses = residualMass; 1175 G4double SumMasses = residualMass;
1217 for ( G4int i = 0; i < numberOfInvolvedNucl 1176 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1218 G4Nucleon* aNucleon = involvedNucleons[i] 1177 G4Nucleon* aNucleon = involvedNucleons[i];
1219 if ( ! aNucleon ) continue; 1178 if ( ! aNucleon ) continue;
1220 SumMasses += aNucleon->GetSplitableHadron 1179 SumMasses += aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
1221 } 1180 }
1222 1181
1223 const G4int maxNumberOfLoops = 1000; 1182 const G4int maxNumberOfLoops = 1000;
1224 G4int loopCounter = 0; 1183 G4int loopCounter = 0;
1225 do { 1184 do {
1226 1185
1227 success = true; 1186 success = true;
1228 G4ThreeVector ptSum( 0.0, 0.0, 0.0 ); 1187 G4ThreeVector ptSum( 0.0, 0.0, 0.0 );
1229 G4double xSum = 0.0; 1188 G4double xSum = 0.0;
1230 1189
1231 for ( G4int i = 0; i < numberOfInvolvedNu 1190 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1232 G4Nucleon* aNucleon = involvedNucleons[ 1191 G4Nucleon* aNucleon = involvedNucleons[i];
1233 if ( ! aNucleon ) continue; 1192 if ( ! aNucleon ) continue;
1234 G4ThreeVector tmpPt = GaussianPt( avera 1193 G4ThreeVector tmpPt = GaussianPt( averagePt2, maxPt2 );
1235 ptSum += tmpPt; 1194 ptSum += tmpPt;
1236 G4ThreeVector tmpX = GaussianPt( dCor*d 1195 G4ThreeVector tmpX = GaussianPt( dCor*dCor, 1.0 );
1237 G4double x = tmpX.x() + 1196 G4double x = tmpX.x() +
1238 aNucleon->GetSplitableHadr 1197 aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()/SumMasses;
1239 if ( x < 0.0 || x > 1.0 ) { 1198 if ( x < 0.0 || x > 1.0 ) {
1240 success = false; 1199 success = false;
1241 break; 1200 break;
1242 } 1201 }
1243 xSum += x; 1202 xSum += x;
1244 //AR The energy is in the lab (instead 1203 //AR The energy is in the lab (instead of cms) frame but it will not be used.
1245 G4LorentzVector tmp( tmpPt.x(), tmpPt.y 1204 G4LorentzVector tmp( tmpPt.x(), tmpPt.y(), x, aNucleon->Get4Momentum().e() );
1246 aNucleon->SetMomentum( tmp ); 1205 aNucleon->SetMomentum( tmp );
1247 } 1206 }
1248 1207
1249 if ( xSum < 0.0 || xSum > 1.0 ) success 1208 if ( xSum < 0.0 || xSum > 1.0 ) success = false;
1250 1209
1251 if ( ! success ) continue; 1210 if ( ! success ) continue;
1252 1211
1253 G4double deltaPx = ( ptSum.x() - pResidua 1212 G4double deltaPx = ( ptSum.x() - pResidual.x() ) / numberOfInvolvedNucleons;
1254 G4double deltaPy = ( ptSum.y() - pResidua 1213 G4double deltaPy = ( ptSum.y() - pResidual.y() ) / numberOfInvolvedNucleons;
1255 G4double delta = 0.0; 1214 G4double delta = 0.0;
1256 if ( residualMassNumber == 0 ) { 1215 if ( residualMassNumber == 0 ) {
1257 delta = ( xSum - 1.0 ) / numberOfInvolv 1216 delta = ( xSum - 1.0 ) / numberOfInvolvedNucleons;
1258 } else { 1217 } else {
1259 delta = 0.0; 1218 delta = 0.0;
1260 } 1219 }
1261 1220
1262 xSum = 1.0; 1221 xSum = 1.0;
1263 mass2 = 0.0; 1222 mass2 = 0.0;
1264 for ( G4int i = 0; i < numberOfInvolvedNu 1223 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1265 G4Nucleon* aNucleon = involvedNucleons[ 1224 G4Nucleon* aNucleon = involvedNucleons[i];
1266 if ( ! aNucleon ) continue; 1225 if ( ! aNucleon ) continue;
1267 G4double x = aNucleon->Get4Momentum().p 1226 G4double x = aNucleon->Get4Momentum().pz() - delta;
1268 xSum -= x; 1227 xSum -= x;
1269 if ( residualMassNumber == 0 ) { 1228 if ( residualMassNumber == 0 ) {
1270 if ( x <= 0.0 || x > 1.0 ) { 1229 if ( x <= 0.0 || x > 1.0 ) {
1271 success = false; 1230 success = false;
1272 break; 1231 break;
1273 } 1232 }
1274 } else { 1233 } else {
1275 if ( x <= 0.0 || x > 1.0 || xSum 1234 if ( x <= 0.0 || x > 1.0 || xSum <= 0.0 || xSum > 1.0 ) {
1276 success = false; 1235 success = false;
1277 break; 1236 break;
1278 } 1237 }
1279 } 1238 }
1280 G4double px = aNucleon->Get4Momentum(). 1239 G4double px = aNucleon->Get4Momentum().px() - deltaPx;
1281 G4double py = aNucleon->Get4Momentum(). 1240 G4double py = aNucleon->Get4Momentum().py() - deltaPy;
1282 mass2 += ( sqr( aNucleon->GetSplitableH 1241 mass2 += ( sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() )
1283 + sqr( px ) + sqr( py ) ) 1242 + sqr( px ) + sqr( py ) ) / x;
1284 G4LorentzVector tmp( px, py, x, aNucleo 1243 G4LorentzVector tmp( px, py, x, aNucleon->Get4Momentum().e() );
1285 aNucleon->SetMomentum( tmp ); 1244 aNucleon->SetMomentum( tmp );
1286 } 1245 }
1287 1246
1288 if ( success && residualMassNumber != 0 1247 if ( success && residualMassNumber != 0 ) {
1289 mass2 += ( sqr( residualMass ) + pResid 1248 mass2 += ( sqr( residualMass ) + pResidual.perp2() ) / xSum;
1290 } 1249 }
1291 1250
1292 #ifdef debugPutOnMassShell 1251 #ifdef debugPutOnMassShell
1293 G4cout << "success " << success << G4en 1252 G4cout << "success " << success << G4endl << " Mt " << std::sqrt( mass2 )/GeV << G4endl;
1294 #endif 1253 #endif
1295 1254
1296 } while ( ( ! success ) && 1255 } while ( ( ! success ) &&
1297 ++loopCounter < maxNumberOfLoops 1256 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
1298 if ( loopCounter >= maxNumberOfLoops ) { 1257 if ( loopCounter >= maxNumberOfLoops ) {
1299 success = false; 1258 success = false;
1300 } 1259 }
1301 1260
1302 return success; 1261 return success;
1303 } 1262 }
1304 1263
1305 1264
1306 //=========================================== 1265 //============================================================================
1307 1266
1308 G4bool G4QGSParticipants:: 1267 G4bool G4QGSParticipants::
1309 CheckKinematics( const G4double sValue, 1268 CheckKinematics( const G4double sValue, // input parameter
1310 const G4double sqrtS, 1269 const G4double sqrtS, // input parameter
1311 const G4double projectileMas 1270 const G4double projectileMass2, // input parameter
1312 const G4double targetMass2, 1271 const G4double targetMass2, // input parameter
1313 const G4double nucleusY, 1272 const G4double nucleusY, // input parameter
1314 const G4bool isProjectileNuc 1273 const G4bool isProjectileNucleus, // input parameter
1315 const G4int numberOfInvolved 1274 const G4int numberOfInvolvedNucleons, // input parameter
1316 G4Nucleon* involvedNucleons[ 1275 G4Nucleon* involvedNucleons[], // input parameter
1317 G4double& targetWminus, 1276 G4double& targetWminus, // output parameter
1318 G4double& projectileWplus, 1277 G4double& projectileWplus, // output parameter
1319 G4bool& success ) { 1278 G4bool& success ) { // input & output parameter
1320 1279
1321 // This method, which is called only by Put 1280 // This method, which is called only by PutOnMassShell, checks whether the
1322 // kinematics is acceptable or not. 1281 // kinematics is acceptable or not.
1323 // This method assumes that all the paramet 1282 // This method assumes that all the parameters have been initialized by the caller;
1324 // notice that the input boolean parameter 1283 // notice that the input boolean parameter isProjectileNucleus is meant to be true
1325 // only in the case of nucleus or antinucle 1284 // only in the case of nucleus or antinucleus projectile.
1326 // The action of this method consists in co 1285 // The action of this method consists in computing targetWminus and projectileWplus
1327 // and setting the parameter success to fal 1286 // and setting the parameter success to false in the case that the kinematics should
1328 // be rejeted. 1287 // be rejeted.
1329 1288
1330 G4double decayMomentum2 = sqr( sValue ) + s 1289 G4double decayMomentum2 = sqr( sValue ) + sqr( projectileMass2 ) + sqr( targetMass2 )
1331 - 2.0*sValue*proj 1290 - 2.0*sValue*projectileMass2 - 2.0*sValue*targetMass2
1332 - 2.0*projectileM 1291 - 2.0*projectileMass2*targetMass2;
1333 targetWminus = ( sValue - projectileMass2 + 1292 targetWminus = ( sValue - projectileMass2 + targetMass2 + std::sqrt( decayMomentum2 ) )
1334 / 2.0 / sqrtS; 1293 / 2.0 / sqrtS;
1335 projectileWplus = sqrtS - targetMass2/targe 1294 projectileWplus = sqrtS - targetMass2/targetWminus;
1336 G4double projectilePz = projectileWplus/2.0 1295 G4double projectilePz = projectileWplus/2.0 - projectileMass2/2.0/projectileWplus;
1337 G4double projectileE = projectileWplus/2.0 1296 G4double projectileE = projectileWplus/2.0 + projectileMass2/2.0/projectileWplus;
1338 G4double projectileY(1.0e5); 1297 G4double projectileY(1.0e5);
1339 if (projectileE - projectilePz > 0.) { 1298 if (projectileE - projectilePz > 0.) {
1340 projectileY = 0.5 * G4Log( (proje 1299 projectileY = 0.5 * G4Log( (projectileE + projectilePz)/
1341 (proje 1300 (projectileE - projectilePz) );
1342 } 1301 }
1343 G4double targetPz = -targetWminus/2.0 + tar 1302 G4double targetPz = -targetWminus/2.0 + targetMass2/2.0/targetWminus;
1344 G4double targetE = targetWminus/2.0 + tar 1303 G4double targetE = targetWminus/2.0 + targetMass2/2.0/targetWminus;
1345 G4double targetY = 0.5 * G4Log( (targetE + 1304 G4double targetY = 0.5 * G4Log( (targetE + targetPz)/(targetE - targetPz) );
1346 1305
1347 #ifdef debugPutOnMassShell 1306 #ifdef debugPutOnMassShell
1348 G4cout << "decayMomentum2 " << decayMomen 1307 G4cout << "decayMomentum2 " << decayMomentum2 << G4endl
1349 << "\t targetWminus projectileWplu 1308 << "\t targetWminus projectileWplus " << targetWminus << " " << projectileWplus << G4endl
1350 << "\t projectileY targetY " << pr 1309 << "\t projectileY targetY " << projectileY << " " << targetY << G4endl;
1351 #endif 1310 #endif
1352 1311
1353 for ( G4int i = 0; i < numberOfInvolvedNucl 1312 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1354 G4Nucleon* aNucleon = involvedNucleons[i] 1313 G4Nucleon* aNucleon = involvedNucleons[i];
1355 if ( ! aNucleon ) continue; 1314 if ( ! aNucleon ) continue;
1356 G4LorentzVector tmp = aNucleon->Get4Momen 1315 G4LorentzVector tmp = aNucleon->Get4Momentum();
1357 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. 1316 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.y() ) +
1358 sqr( aNucleon->GetSplitabl 1317 sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() );
1359 G4double x = tmp.z(); 1318 G4double x = tmp.z();
1360 G4double pz = -targetWminus*x/2.0 + mt2/( 1319 G4double pz = -targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
1361 G4double e = targetWminus*x/2.0 + mt2/( 1320 G4double e = targetWminus*x/2.0 + mt2/(2.0*targetWminus*x);
1362 if ( isProjectileNucleus ) { 1321 if ( isProjectileNucleus ) {
1363 pz = projectileWplus*x/2.0 - mt2/(2.0*p 1322 pz = projectileWplus*x/2.0 - mt2/(2.0*projectileWplus*x);
1364 e = projectileWplus*x/2.0 + mt2/(2.0*p 1323 e = projectileWplus*x/2.0 + mt2/(2.0*projectileWplus*x);
1365 } 1324 }
1366 G4double nucleonY = 0.5 * G4Log( (e + pz) 1325 G4double nucleonY = 0.5 * G4Log( (e + pz)/(e - pz) );
1367 1326
1368 #ifdef debugPutOnMassShell 1327 #ifdef debugPutOnMassShell
1369 G4cout << "i nY pY nY-AY AY " << i << " 1328 G4cout << "i nY pY nY-AY AY " << i << " " << nucleonY << " " << projectileY <<G4endl;
1370 #endif 1329 #endif
1371 1330
1372 if ( std::abs( nucleonY - nucleusY ) > 2 1331 if ( std::abs( nucleonY - nucleusY ) > 2 ||
1373 ( isProjectileNucleus && targetY > 1332 ( isProjectileNucleus && targetY > nucleonY ) ||
1374 ( ! isProjectileNucleus && project 1333 ( ! isProjectileNucleus && projectileY < nucleonY ) ) {
1375 success = false; 1334 success = false;
1376 break; 1335 break;
1377 } 1336 }
1378 } 1337 }
1379 return true; 1338 return true;
1380 } 1339 }
1381 1340
1382 1341
1383 //=========================================== 1342 //============================================================================
1384 1343
1385 G4bool G4QGSParticipants:: 1344 G4bool G4QGSParticipants::
1386 FinalizeKinematics( const G4double w, 1345 FinalizeKinematics( const G4double w, // input parameter
1387 const G4bool isProjectile 1346 const G4bool isProjectileNucleus, // input parameter
1388 const G4LorentzRotation& 1347 const G4LorentzRotation& boostFromCmsToLab, // input parameter
1389 const G4double residualMa 1348 const G4double residualMass, // input parameter
1390 const G4int residualMassN 1349 const G4int residualMassNumber, // input parameter
1391 const G4int numberOfInvol 1350 const G4int numberOfInvolvedNucleons, // input parameter
1392 G4Nucleon* involvedNucleo 1351 G4Nucleon* involvedNucleons[], // input & output parameter
1393 G4LorentzVector& residual4Momen 1352 G4LorentzVector& residual4Momentum ) { // output parameter
1394 1353
1395 // This method, which is called only by Put 1354 // This method, which is called only by PutOnMassShell, finalizes the kinematics:
1396 // this method is called when we are sure t 1355 // this method is called when we are sure that the sampling of the kinematics is
1397 // acceptable. 1356 // acceptable.
1398 // This method assumes that all the paramet 1357 // This method assumes that all the parameters have been initialized by the caller;
1399 // notice that the input boolean parameter 1358 // notice that the input boolean parameter isProjectileNucleus is meant to be true
1400 // only in the case of nucleus or antinucle 1359 // only in the case of nucleus or antinucleus projectile: this information is needed
1401 // because the sign of pz (in the center-of 1360 // because the sign of pz (in the center-of-mass frame) in this case is opposite
1402 // with respect to the case of a normal had 1361 // with respect to the case of a normal hadron projectile.
1403 // The action of this method consists in mo 1362 // The action of this method consists in modifying the momenta of the nucleons
1404 // (in the lab frame) and computing the res 1363 // (in the lab frame) and computing the residual 4-momentum (in the center-of-mass
1405 // frame). 1364 // frame).
1406 1365
1407 G4ThreeVector residual3Momentum( 0.0, 0.0, 1366 G4ThreeVector residual3Momentum( 0.0, 0.0, 1.0 );
1408 1367
1409 for ( G4int i = 0; i < numberOfInvolvedNucl 1368 for ( G4int i = 0; i < numberOfInvolvedNucleons; i++ ) {
1410 G4Nucleon* aNucleon = involvedNucleons[i] 1369 G4Nucleon* aNucleon = involvedNucleons[i];
1411 if ( ! aNucleon ) continue; 1370 if ( ! aNucleon ) continue;
1412 G4LorentzVector tmp = aNucleon->Get4Momen 1371 G4LorentzVector tmp = aNucleon->Get4Momentum();
1413 residual3Momentum -= tmp.vect(); 1372 residual3Momentum -= tmp.vect();
1414 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. 1373 G4double mt2 = sqr( tmp.x() ) + sqr( tmp.y() ) +
1415 sqr( aNucleon->GetSplitabl 1374 sqr( aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass() );
1416 G4double x = tmp.z(); 1375 G4double x = tmp.z();
1417 G4double pz = -w * x / 2.0 + mt2 / ( 2. 1376 G4double pz = -w * x / 2.0 + mt2 / ( 2.0 * w * x );
1418 G4double e = w * x / 2.0 + mt2 / ( 2. 1377 G4double e = w * x / 2.0 + mt2 / ( 2.0 * w * x );
1419 // Reverse the sign of pz in the case of 1378 // Reverse the sign of pz in the case of nucleus or antinucleus projectile
1420 if ( isProjectileNucleus ) pz *= -1.0; 1379 if ( isProjectileNucleus ) pz *= -1.0;
1421 tmp.setPz( pz ); 1380 tmp.setPz( pz );
1422 tmp.setE( e ); 1381 tmp.setE( e );
1423 tmp.transform( boostFromCmsToLab ); 1382 tmp.transform( boostFromCmsToLab );
1424 aNucleon->SetMomentum( tmp ); 1383 aNucleon->SetMomentum( tmp );
1425 G4VSplitableHadron* splitableHadron = aNu 1384 G4VSplitableHadron* splitableHadron = aNucleon->GetSplitableHadron();
1426 splitableHadron->Set4Momentum( tmp ); 1385 splitableHadron->Set4Momentum( tmp );
1427 #ifdef debugPutOnMassShell 1386 #ifdef debugPutOnMassShell
1428 G4cout << "Target involved nucleon No, 1387 G4cout << "Target involved nucleon No, name, 4Mom "
1429 << i<<" "<<aNucleon->GetDefinitio 1388 << i<<" "<<aNucleon->GetDefinition()->GetParticleName()<<" "<<tmp<< G4endl;
1430 #endif 1389 #endif
1431 } 1390 }
1432 1391
1433 G4double residualMt2 = sqr( residualMass ) 1392 G4double residualMt2 = sqr( residualMass ) + sqr( residual3Momentum.x() )
1434 + sqr( residual3Moment 1393 + sqr( residual3Momentum.y() );
1435 1394
1436 #ifdef debugPutOnMassShell 1395 #ifdef debugPutOnMassShell
1437 G4cout <<G4endl<< "w residual3Momentum.z( 1396 G4cout <<G4endl<< "w residual3Momentum.z() " << w << " " << residual3Momentum.z() << G4endl;
1438 #endif 1397 #endif
1439 1398
1440 G4double residualPz = 0.0; 1399 G4double residualPz = 0.0;
1441 G4double residualE = 0.0; 1400 G4double residualE = 0.0;
1442 if ( residualMassNumber != 0 ) { 1401 if ( residualMassNumber != 0 ) {
1443 residualPz = -w * residual3Momentum.z() / 1402 residualPz = -w * residual3Momentum.z() / 2.0 +
1444 residualMt2 / ( 2.0 * w * r 1403 residualMt2 / ( 2.0 * w * residual3Momentum.z() );
1445 residualE = w * residual3Momentum.z() / 1404 residualE = w * residual3Momentum.z() / 2.0 +
1446 residualMt2 / ( 2.0 * w * r 1405 residualMt2 / ( 2.0 * w * residual3Momentum.z() );
1447 // Reverse the sign of residualPz in the 1406 // Reverse the sign of residualPz in the case of nucleus or antinucleus projectile
1448 if ( isProjectileNucleus ) residualPz *= 1407 if ( isProjectileNucleus ) residualPz *= -1.0;
1449 } 1408 }
1450 1409
1451 residual4Momentum.setPx( residual3Momentum. 1410 residual4Momentum.setPx( residual3Momentum.x() );
1452 residual4Momentum.setPy( residual3Momentum. 1411 residual4Momentum.setPy( residual3Momentum.y() );
1453 residual4Momentum.setPz( residualPz ); 1412 residual4Momentum.setPz( residualPz );
1454 residual4Momentum.setE( residualE ); 1413 residual4Momentum.setE( residualE );
1455 1414
1456 return true; 1415 return true;
1457 } 1416 }
1458 1417
1459 //=========================================== 1418 //======================================================
1460 void G4QGSParticipants::PerformDiffractiveCol 1419 void G4QGSParticipants::PerformDiffractiveCollisions()
1461 { 1420 {
1462 #ifdef debugQGSParticipants 1421 #ifdef debugQGSParticipants
1463 G4cout<<G4endl<<"PerformDiffractiveCollisi 1422 G4cout<<G4endl<<"PerformDiffractiveCollisions()......"<<G4endl
1464 <<"theInteractions.size() "<<theInte 1423 <<"theInteractions.size() "<<theInteractions.size()<<G4endl;
1465 #endif 1424 #endif
1466 1425
1467 unsigned int i; 1426 unsigned int i;
1468 for (i = 0; i < theInteractions.size(); i++ 1427 for (i = 0; i < theInteractions.size(); i++)
1469 { 1428 {
1470 G4InteractionContent* anIniteraction = th 1429 G4InteractionContent* anIniteraction = theInteractions[i];
1471 #ifdef debugQGSParticipants 1430 #ifdef debugQGSParticipants
1472 G4cout<<"Interaction # and its status " 1431 G4cout<<"Interaction # and its status "
1473 <<i<<" "<<theInteractions[i]->Get 1432 <<i<<" "<<theInteractions[i]->GetStatus()<<G4endl;
1474 #endif 1433 #endif
1475 1434
1476 G4int InterStatus = theInteractions[i]->G 1435 G4int InterStatus = theInteractions[i]->GetStatus();
1477 if ( (InterStatus == PrD) || (InterStatus 1436 if ( (InterStatus == PrD) || (InterStatus == TrD) || (InterStatus == DD))
1478 { // Selection of diffractive interactio 1437 { // Selection of diffractive interactions
1479 #ifdef debugQGSParticipants 1438 #ifdef debugQGSParticipants
1480 G4cout<<"Simulation of diffractive in 1439 G4cout<<"Simulation of diffractive interaction. "<<InterStatus <<" PrD/TrD/DD/ND/Qech - 0,1,2,3,4"<<G4endl;
1481 #endif 1440 #endif
1482 1441
1483 G4VSplitableHadron* aTarget = anInitera 1442 G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
1484 1443
1485 #ifdef debugQGSParticipants 1444 #ifdef debugQGSParticipants
1486 G4cout<<"The proj. before inter " 1445 G4cout<<"The proj. before inter "
1487 <<theProjectileSplitable->Get4M 1446 <<theProjectileSplitable->Get4Momentum()<<" "
1488 <<theProjectileSplitable->Get4M 1447 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1489 G4cout<<"The targ. before inter " <<a 1448 G4cout<<"The targ. before inter " <<aTarget->Get4Momentum()<<" "
1490 <<aTarget->Get4Momentum().mag() 1449 <<aTarget->Get4Momentum().mag()<<G4endl;
1491 #endif 1450 #endif
1492 1451
1493 if ( InterStatus == PrD ) 1452 if ( InterStatus == PrD )
1494 theSingleDiffExcitation.ExcitePartici 1453 theSingleDiffExcitation.ExciteParticipants(theProjectileSplitable, aTarget, TRUE);
1495 1454
1496 if ( InterStatus == TrD ) 1455 if ( InterStatus == TrD )
1497 theSingleDiffExcitation.ExcitePartici 1456 theSingleDiffExcitation.ExciteParticipants(theProjectileSplitable, aTarget, FALSE);
1498 1457
1499 if ( InterStatus == DD ) 1458 if ( InterStatus == DD )
1500 theDiffExcitaton.ExciteParticipants(t 1459 theDiffExcitaton.ExciteParticipants(theProjectileSplitable, aTarget);
1501 1460
1502 #ifdef debugQGSParticipants 1461 #ifdef debugQGSParticipants
1503 G4cout<<"The proj. after inter " <<t 1462 G4cout<<"The proj. after inter " <<theProjectileSplitable->Get4Momentum()<<" "
1504 <<theProjectileSplitable->Get4M 1463 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1505 G4cout<<"The targ. after inter " <<aTarget 1464 G4cout<<"The targ. after inter " <<aTarget->Get4Momentum()<<" "
1506 <<aTarget->Get4Momentum().mag() 1465 <<aTarget->Get4Momentum().mag()<<G4endl;
1507 #endif 1466 #endif
1508 } 1467 }
1509 1468
1510 if ( InterStatus == Qexc ) 1469 if ( InterStatus == Qexc )
1511 { // Quark exchange process 1470 { // Quark exchange process
1512 #ifdef debugQGSParticipants 1471 #ifdef debugQGSParticipants
1513 G4cout<<"Simulation of interaction wi 1472 G4cout<<"Simulation of interaction with quark exchange."<<G4endl;
1514 #endif 1473 #endif
1515 G4VSplitableHadron* aTarget = anInitera 1474 G4VSplitableHadron* aTarget = anIniteraction->GetTarget();
1516 1475
1517 #ifdef debugQGSParticipants 1476 #ifdef debugQGSParticipants
1518 G4cout<<"The proj. before inter " <<t 1477 G4cout<<"The proj. before inter " <<theProjectileSplitable->Get4Momentum()<<" "
1519 <<theProjectileSplitable->Get4M 1478 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1520 G4cout<<"The targ. before inter "<<aT 1479 G4cout<<"The targ. before inter "<<aTarget->Get4Momentum()<<" "
1521 <<aTarget->Get4Momentum().mag() 1480 <<aTarget->Get4Momentum().mag()<<G4endl;
1522 #endif 1481 #endif
1523 1482
1524 theQuarkExchange.ExciteParticipants(the 1483 theQuarkExchange.ExciteParticipants(theProjectileSplitable, aTarget);
1525 1484
1526 #ifdef debugQGSParticipants 1485 #ifdef debugQGSParticipants
1527 G4cout<<"The proj. after inter " <<t 1486 G4cout<<"The proj. after inter " <<theProjectileSplitable->Get4Momentum()<<" "
1528 <<theProjectileSplitable->Get4M 1487 <<theProjectileSplitable->Get4Momentum().mag()<<G4endl;
1529 G4cout<<"The targ. after inter " <<aTarget 1488 G4cout<<"The targ. after inter " <<aTarget->Get4Momentum()<<" "
1530 <<aTarget->Get4Momentum().mag() 1489 <<aTarget->Get4Momentum().mag()<<G4endl;
1531 #endif 1490 #endif
1532 } 1491 }
1533 } 1492 }
1534 } 1493 }
1535 1494
1536 //=========================================== 1495 //======================================================
1537 G4bool G4QGSParticipants::DeterminePartonMome 1496 G4bool G4QGSParticipants::DeterminePartonMomenta()
1538 { 1497 {
1539 if ( ! theProjectileSplitable ) return fals 1498 if ( ! theProjectileSplitable ) return false;
1540 1499
1541 const G4double aHugeValue = 1.0e+10; 1500 const G4double aHugeValue = 1.0e+10;
1542 1501
1543 #ifdef debugQGSParticipants 1502 #ifdef debugQGSParticipants
1544 G4cout<<G4endl<<"DeterminePartonMomenta() 1503 G4cout<<G4endl<<"DeterminePartonMomenta()......"<<G4endl;
1545 G4cout<<"theProjectile status (0 -nondiff 1504 G4cout<<"theProjectile status (0 -nondiffr, #0 diffr./reggeon): "<<theProjectileSplitable->GetStatus()<<G4endl;
1546 #endif 1505 #endif
1547 1506
1548 if (theProjectileSplitable->GetStatus() != 1507 if (theProjectileSplitable->GetStatus() != 0) {return false;} // There were only diffractive interactions.
1549 1508
1550 G4LorentzVector Projectile4Momentum = theP 1509 G4LorentzVector Projectile4Momentum = theProjectileSplitable->Get4Momentum();
1551 G4LorentzVector Psum = Projectile4Momentum; 1510 G4LorentzVector Psum = Projectile4Momentum;
1552 1511
1553 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(35 1512 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(350.), VqqM_tr(700);
1554 if (std::abs(theProjectile.GetDefinition()- 1513 if (std::abs(theProjectile.GetDefinition()->GetBaryonNumber()) != 0) {VqM_pr=350*MeV; VaqM_pr=700*MeV;}
1555 1514
1556 #ifdef debugQGSParticipants 1515 #ifdef debugQGSParticipants
1557 G4cout<<"Projectile 4 momentum "<<Psum<<G 1516 G4cout<<"Projectile 4 momentum "<<Psum<<G4endl
1558 <<"Target nucleon momenta at start" 1517 <<"Target nucleon momenta at start"<<G4endl;
1559 G4int NuclNo=0; <<
1560 #endif 1518 #endif
1561 1519
1562 std::vector<G4VSplitableHadron*>::iterator 1520 std::vector<G4VSplitableHadron*>::iterator i;
>> 1521 G4int NuclNo=0;
1563 1522
1564 for (i = theTargets.begin(); i != theTarget 1523 for (i = theTargets.begin(); i != theTargets.end(); i++ )
1565 { 1524 {
1566 Psum += (*i)->Get4Momentum(); 1525 Psum += (*i)->Get4Momentum();
1567 #ifdef debugQGSParticipants 1526 #ifdef debugQGSParticipants
1568 G4cout<<"Nusleus nucleon # and its 4Mom 1527 G4cout<<"Nusleus nucleon # and its 4Mom. "<<NuclNo<<" "<<(*i)->Get4Momentum()<<G4endl;
1569 NuclNo++; <<
1570 #endif 1528 #endif
>> 1529 NuclNo++;
1571 } 1530 }
1572 1531
1573 G4LorentzRotation toCms( -1*Psum.boostVecto 1532 G4LorentzRotation toCms( -1*Psum.boostVector() );
1574 1533
1575 G4LorentzVector Ptmp = toCms*Projectile4Mom 1534 G4LorentzVector Ptmp = toCms*Projectile4Momentum;
1576 1535
1577 toCms.rotateZ( -1*Ptmp.phi() ); 1536 toCms.rotateZ( -1*Ptmp.phi() );
1578 toCms.rotateY( -1*Ptmp.theta() ); 1537 toCms.rotateY( -1*Ptmp.theta() );
1579 G4LorentzRotation toLab(toCms.inverse()); 1538 G4LorentzRotation toLab(toCms.inverse());
1580 Projectile4Momentum.transform( toCms ); 1539 Projectile4Momentum.transform( toCms );
1581 // Ptarget.transform( toCms ); 1540 // Ptarget.transform( toCms );
1582 1541
1583 #ifdef debugQGSParticipants 1542 #ifdef debugQGSParticipants
1584 G4cout<<G4endl<<"In CMS---------------"<< 1543 G4cout<<G4endl<<"In CMS---------------"<<G4endl;
1585 G4cout<<"Projectile 4 Mom "<<Projectile4M 1544 G4cout<<"Projectile 4 Mom "<<Projectile4Momentum<<G4endl;
1586 NuclNo=0; <<
1587 #endif 1545 #endif
1588 1546
>> 1547 NuclNo=0;
1589 G4LorentzVector Target4Momentum(0.,0.,0.,0. 1548 G4LorentzVector Target4Momentum(0.,0.,0.,0.);
1590 for(i = theTargets.begin(); i != theTargets 1549 for(i = theTargets.begin(); i != theTargets.end(); i++ )
1591 { 1550 {
1592 G4LorentzVector tmp= (*i)->Get4Momentum() 1551 G4LorentzVector tmp= (*i)->Get4Momentum(); tmp.transform( toCms );
1593 (*i)->Set4Momentum( tmp ); 1552 (*i)->Set4Momentum( tmp );
1594 #ifdef debugQGSParticipants 1553 #ifdef debugQGSParticipants
1595 G4cout<<"Target nucleon # and 4Mom "<<" 1554 G4cout<<"Target nucleon # and 4Mom "<<" "<<NuclNo<<" "<<(*i)->Get4Momentum()<<G4endl;
1596 NuclNo++; <<
1597 #endif 1555 #endif
1598 Target4Momentum += tmp; 1556 Target4Momentum += tmp;
>> 1557 NuclNo++;
1599 } 1558 }
1600 1559
1601 G4double S = Psum.mag2(); 1560 G4double S = Psum.mag2();
1602 G4double SqrtS = std::sqrt(S); 1561 G4double SqrtS = std::sqrt(S);
1603 1562
1604 #ifdef debugQGSParticipants 1563 #ifdef debugQGSParticipants
1605 G4cout<<"Sum of target nucleons 4 momentu 1564 G4cout<<"Sum of target nucleons 4 momentum "<<Target4Momentum<<G4endl<<G4endl;
1606 G4cout<<"Target nucleons mom: px, py, z_1 1565 G4cout<<"Target nucleons mom: px, py, z_1, m_i"<<G4endl;
1607 NuclNo=0; <<
1608 #endif 1566 #endif
1609 1567
1610 //G4double PplusProjectile = Projectile4Mom 1568 //G4double PplusProjectile = Projectile4Momentum.plus();
1611 G4double PminusTarget = Target4Momentum. 1569 G4double PminusTarget = Target4Momentum.minus();
>> 1570 NuclNo=0;
1612 1571
1613 for(i = theTargets.begin(); i != theTargets 1572 for(i = theTargets.begin(); i != theTargets.end(); i++ )
1614 { 1573 {
1615 G4LorentzVector tmp = (*i)->Get4Momentum( 1574 G4LorentzVector tmp = (*i)->Get4Momentum(); // tmp.boost(bstToCM);
1616 1575
1617 //AR-19Jan2017 : the following line is ca 1576 //AR-19Jan2017 : the following line is causing a strange crash when Geant4
1618 // is built in optimized mo 1577 // is built in optimized mode.
1619 // To fix it, I get the mas 1578 // To fix it, I get the mass square instead of directly the
1620 // mass from the Lorentz ve 1579 // mass from the Lorentz vector, and then I take care of the
1621 // square root. If the mass 1580 // square root. If the mass square is negative, a JustWarning
1622 // exception is thrown, and 1581 // exception is thrown, and the mass is set to 0.
1623 //G4double Mass = tmp.mag(); 1582 //G4double Mass = tmp.mag();
1624 G4double Mass2 = tmp.mag2(); 1583 G4double Mass2 = tmp.mag2();
1625 G4double Mass = 0.0; 1584 G4double Mass = 0.0;
1626 if ( Mass2 < 0.0 ) { 1585 if ( Mass2 < 0.0 ) {
1627 G4ExceptionDescription ed; 1586 G4ExceptionDescription ed;
1628 ed << "Projectile " << theProjectile.Ge 1587 ed << "Projectile " << theProjectile.GetDefinition()->GetParticleName()
1629 << "Â 4-momentum " << Psum << G4end 1588 << "Â 4-momentum " << Psum << G4endl;
1630 ed << "LorentzVector tmp " << tmp << " 1589 ed << "LorentzVector tmp " << tmp << " with Mass2 " << Mass2 << G4endl;
1631 G4Exception( "G4QGSParticipants::Determ 1590 G4Exception( "G4QGSParticipants::DeterminePartonMomenta(): 4-momentum with negative mass!",
1632 "HAD_QGSPARTICIPANTS_001", 1591 "HAD_QGSPARTICIPANTS_001", JustWarning, ed );
1633 } else { 1592 } else {
1634 Mass = std::sqrt( Mass2 ); 1593 Mass = std::sqrt( Mass2 );
1635 } 1594 }
1636 1595
1637 tmp.setPz(tmp.minus()/PminusTarget); tm 1596 tmp.setPz(tmp.minus()/PminusTarget); tmp.setE(Mass);
1638 (*i)->Set4Momentum(tmp); 1597 (*i)->Set4Momentum(tmp);
1639 #ifdef debugQGSParticipants 1598 #ifdef debugQGSParticipants
1640 G4cout<<"Target nucleons # and mom: "<< 1599 G4cout<<"Target nucleons # and mom: "<<NuclNo<<" "<<(*i)->Get4Momentum()<<G4endl;
1641 NuclNo++; <<
1642 #endif 1600 #endif
>> 1601 NuclNo++;
1643 } 1602 }
1644 1603
1645 //+++++++++++++++++++++++++++++++++++++++++ 1604 //+++++++++++++++++++++++++++++++++++++++++++
1646 1605
1647 G4double SigPt = sigmaPt; 1606 G4double SigPt = sigmaPt;
1648 G4Parton* aParton(0); 1607 G4Parton* aParton(0);
1649 G4ThreeVector aPtVector(0.,0.,0.); 1608 G4ThreeVector aPtVector(0.,0.,0.);
1650 G4LorentzVector tmp(0.,0.,0.,0.); 1609 G4LorentzVector tmp(0.,0.,0.,0.);
1651 1610
1652 G4double Mt(0.); 1611 G4double Mt(0.);
1653 G4double ProjSumMt(0.), ProjSumMt2perX(0.); 1612 G4double ProjSumMt(0.), ProjSumMt2perX(0.);
1654 G4double TargSumMt(0.), TargSumMt2perX(0.); 1613 G4double TargSumMt(0.), TargSumMt2perX(0.);
1655 1614
1656 1615
1657 G4double aBeta = beta; // Member of the c 1616 G4double aBeta = beta; // Member of the class
1658 1617
1659 const G4ParticleDefinition* theProjectileDe 1618 const G4ParticleDefinition* theProjectileDefinition = theProjectileSplitable->GetDefinition();
1660 if (theProjectileDefinition == G4PionMinus: 1619 if (theProjectileDefinition == G4PionMinus::PionMinusDefinition()) aBeta = 1.;
1661 if (theProjectileDefinition == G4Gamma::Gam 1620 if (theProjectileDefinition == G4Gamma::GammaDefinition()) aBeta = 1.;
1662 if (theProjectileDefinition == G4PionPlus:: 1621 if (theProjectileDefinition == G4PionPlus::PionPlusDefinition()) aBeta = 1.;
1663 if (theProjectileDefinition == G4PionZero:: 1622 if (theProjectileDefinition == G4PionZero::PionZeroDefinition()) aBeta = 1.;
1664 if (theProjectileDefinition == G4KaonPlus:: 1623 if (theProjectileDefinition == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.;
1665 if (theProjectileDefinition == G4KaonMinus: 1624 if (theProjectileDefinition == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.;
1666 1625
1667 G4double Xmin = 0.; 1626 G4double Xmin = 0.;
1668 1627
1669 G4bool Success = true; G4int attempt = 0; 1628 G4bool Success = true; G4int attempt = 0;
1670 const G4int maxNumberOfAttempts = 1000; 1629 const G4int maxNumberOfAttempts = 1000;
1671 do 1630 do
1672 { 1631 {
1673 attempt++; if( attempt == 100*(attempt/ 1632 attempt++; if( attempt == 100*(attempt/100) ) {SigPt/=2.;}
1674 1633
1675 ProjSumMt=0.; ProjSumMt2perX=0.; 1634 ProjSumMt=0.; ProjSumMt2perX=0.;
1676 TargSumMt=0.; TargSumMt2perX=0.; 1635 TargSumMt=0.; TargSumMt2perX=0.;
1677 1636
1678 Success = true; 1637 Success = true;
1679 G4int nSeaPair = theProjectileSplitable-> 1638 G4int nSeaPair = theProjectileSplitable->GetSoftCollisionCount()-1;
1680 #ifdef debugQGSParticipants 1639 #ifdef debugQGSParticipants
1681 G4cout<<"attempt ---------------------- 1640 G4cout<<"attempt ------------------------ "<<attempt<<G4endl;
1682 G4cout<<"nSeaPair of proj "<<nSeaPair<< 1641 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4endl;
1683 #endif 1642 #endif
1684 1643
1685 G4double SumPx = 0.; 1644 G4double SumPx = 0.;
1686 G4double SumPy = 0.; 1645 G4double SumPy = 0.;
1687 G4double SumZ = 0.; 1646 G4double SumZ = 0.;
1688 G4int NumberOfUnsampledSeaQuarks = 2*nSea 1647 G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1689 1648
1690 G4double Qmass=0.; 1649 G4double Qmass=0.;
1691 for (G4int aSeaPair = 0; aSeaPair < nSeaP 1650 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1692 { 1651 {
1693 aParton = theProjectileSplitable->GetNe 1652 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1694 #ifdef debugQGSParticipants 1653 #ifdef debugQGSParticipants
1695 G4cout<<"Sea quarks: "<<aSeaPair<<" " 1654 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1696 #endif 1655 #endif
1697 aPtVector = GaussianPt(SigPt, aHugeValu 1656 aPtVector = GaussianPt(SigPt, aHugeValue);
1698 tmp.setPx(aPtVector.x()); tmp.setPy(aPt 1657 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1699 SumPx += aPtVector.x(); SumPy += aPtV 1658 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1700 Mt = std::sqrt(aPtVector.mag2()+sqr(Qma 1659 Mt = std::sqrt(aPtVector.mag2()+sqr(Qmass));
1701 ProjSumMt += Mt; 1660 ProjSumMt += Mt;
1702 1661
1703 // Sampling of Z fraction 1662 // Sampling of Z fraction
1704 tmp.setPz(SampleX(Xmin, NumberOfUnsampl 1663 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1705 SumZ += tmp.z(); 1664 SumZ += tmp.z();
1706 1665
1707 NumberOfUnsampledSeaQuarks--; 1666 NumberOfUnsampledSeaQuarks--;
1708 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); 1667 ProjSumMt2perX +=sqr(Mt)/tmp.pz();
1709 tmp.setE(sqr(Mt)); 1668 tmp.setE(sqr(Mt));
1710 aParton->Set4Momentum(tmp); 1669 aParton->Set4Momentum(tmp);
1711 1670
1712 aParton = theProjectileSplitable->GetNe 1671 aParton = theProjectileSplitable->GetNextAntiParton(); // for anti-quarks
1713 #ifdef debugQGSParticipants 1672 #ifdef debugQGSParticipants
1714 G4cout<<" "<<aParton->GetDefinition() 1673 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1715 G4cout<<" "<<tmp<<" "<<S 1674 G4cout<<" "<<tmp<<" "<<SumZ<<G4endl;
1716 #endif 1675 #endif
1717 aPtVector = GaussianPt(SigPt, aHugeValu 1676 aPtVector = GaussianPt(SigPt, aHugeValue);
1718 tmp.setPx(aPtVector.x()); tmp.setPy(aPt 1677 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1719 SumPx += aPtVector.x(); SumPy += aPtV 1678 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1720 Mt = std::sqrt(aPtVector.mag2()+sqr(Qma 1679 Mt = std::sqrt(aPtVector.mag2()+sqr(Qmass));
1721 ProjSumMt += Mt; 1680 ProjSumMt += Mt;
1722 1681
1723 // Sampling of Z fraction 1682 // Sampling of Z fraction
1724 tmp.setPz(SampleX(Xmin, NumberOfUnsampl 1683 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1725 SumZ += tmp.z(); 1684 SumZ += tmp.z();
1726 1685
1727 NumberOfUnsampledSeaQuarks--; 1686 NumberOfUnsampledSeaQuarks--;
1728 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); 1687 ProjSumMt2perX +=sqr(Mt)/tmp.pz();
1729 tmp.setE(sqr(Mt)); 1688 tmp.setE(sqr(Mt));
1730 aParton->Set4Momentum(tmp); 1689 aParton->Set4Momentum(tmp);
1731 #ifdef debugQGSParticipants 1690 #ifdef debugQGSParticipants
1732 G4cout<<" "<<tmp<<" "<<S 1691 G4cout<<" "<<tmp<<" "<<SumZ<<G4endl;
1733 #endif 1692 #endif
1734 } 1693 }
1735 1694
1736 // For valence quark 1695 // For valence quark
1737 aParton = theProjectileSplitable->GetNext 1696 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1738 #ifdef debugQGSParticipants 1697 #ifdef debugQGSParticipants
1739 G4cout<<"Val quark of Pr"<<" "<<aParton 1698 G4cout<<"Val quark of Pr"<<" "<<aParton->GetDefinition()->GetParticleName();
1740 #endif 1699 #endif
1741 aPtVector = GaussianPt(SigPt, aHugeValue) 1700 aPtVector = GaussianPt(SigPt, aHugeValue);
1742 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVe 1701 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1743 SumPx += aPtVector.x(); SumPy += aPtVec 1702 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1744 Mt = std::sqrt(aPtVector.mag2()+sqr(VqM_p 1703 Mt = std::sqrt(aPtVector.mag2()+sqr(VqM_pr));
1745 ProjSumMt += Mt; 1704 ProjSumMt += Mt;
1746 1705
1747 // Sampling of Z fraction 1706 // Sampling of Z fraction
1748 tmp.setPz(SampleX(Xmin, NumberOfUnsampled 1707 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1749 SumZ += tmp.z(); 1708 SumZ += tmp.z();
1750 1709
1751 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); 1710 ProjSumMt2perX +=sqr(Mt)/tmp.pz();
1752 tmp.setE(sqr(Mt)); 1711 tmp.setE(sqr(Mt));
1753 aParton->Set4Momentum(tmp); 1712 aParton->Set4Momentum(tmp);
1754 1713
1755 // For valence di-quark 1714 // For valence di-quark
1756 aParton = theProjectileSplitable->GetNext 1715 aParton = theProjectileSplitable->GetNextAntiParton();
1757 #ifdef debugQGSParticipants 1716 #ifdef debugQGSParticipants
1758 G4cout<<" "<<aParton->GetDefinition()-> 1717 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1759 G4cout<<" "<<tmp<<" "<<Sum 1718 G4cout<<" "<<tmp<<" "<<SumZ<<" (z-fraction)"<<G4endl;
1760 #endif 1719 #endif
1761 tmp.setPx(-SumPx); tmp.setPy(-SumPy); 1720 tmp.setPx(-SumPx); tmp.setPy(-SumPy);
1762 //Uzhi 2019 Mt = std::sqrt(aPtVector.mag 1721 //Uzhi 2019 Mt = std::sqrt(aPtVector.mag2()+sqr(VaqM_pr));
1763 Mt = std::sqrt( sqr(SumPx) + sqr(SumPy) + 1722 Mt = std::sqrt( sqr(SumPx) + sqr(SumPy) + sqr(VaqM_pr) ); //Uzhi 2019
1764 ProjSumMt += Mt; 1723 ProjSumMt += Mt;
1765 tmp.setPz(1.-SumZ); 1724 tmp.setPz(1.-SumZ);
1766 1725
1767 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); // QQ 1726 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); // QQmass=750 MeV
1768 tmp.setE(sqr(Mt)); 1727 tmp.setE(sqr(Mt));
1769 aParton->Set4Momentum(tmp); 1728 aParton->Set4Momentum(tmp);
1770 #ifdef debugQGSParticipants 1729 #ifdef debugQGSParticipants
1771 G4cout<<" "<<tmp<<" "<<Sum 1730 G4cout<<" "<<tmp<<" "<<SumZ+(1.-SumZ)<<" (z-fraction)"<<G4endl;
1772 NuclNo=0; <<
1773 #endif 1731 #endif
1774 1732
1775 // End of work with the projectile 1733 // End of work with the projectile
1776 1734
1777 // Work with target nucleons 1735 // Work with target nucleons
1778 1736
>> 1737 NuclNo=0;
1779 for(i = theTargets.begin(); i != theTarge 1738 for(i = theTargets.begin(); i != theTargets.end(); i++ )
1780 { 1739 {
1781 nSeaPair = (*i)->GetSoftCollisionCount( 1740 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1782 #ifdef debugQGSParticipants 1741 #ifdef debugQGSParticipants
1783 G4cout<<"nSeaPair of target N "<<nSea 1742 G4cout<<"nSeaPair of target N "<<nSeaPair<<G4endl
1784 <<"Target nucleon 4Mom "<<(*i)- 1743 <<"Target nucleon 4Mom "<<(*i)->Get4Momentum()<<G4endl;
1785 #endif 1744 #endif
1786 1745
1787 SumPx = (*i)->Get4Momentum().px() * (-1 1746 SumPx = (*i)->Get4Momentum().px() * (-1.);
1788 SumPy = (*i)->Get4Momentum().py() * (-1 1747 SumPy = (*i)->Get4Momentum().py() * (-1.);
1789 SumZ = 0.; 1748 SumZ = 0.;
1790 1749
>> 1750 G4double SumZw=0.;
1791 NumberOfUnsampledSeaQuarks = 2*nSeaPair 1751 NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1792 1752
1793 Qmass=0; 1753 Qmass=0;
1794 for (G4int aSeaPair = 0; aSeaPair < nSe 1754 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1795 { 1755 {
1796 aParton = (*i)->GetNextParton(); // 1756 aParton = (*i)->GetNextParton(); // for quarks
1797 #ifdef debugQGSParticipants 1757 #ifdef debugQGSParticipants
1798 G4cout<<"Sea quarks: "<<aSeaPair<<" 1758 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1799 #endif 1759 #endif
1800 aPtVector = GaussianPt(SigPt, aHugeVa 1760 aPtVector = GaussianPt(SigPt, aHugeValue);
1801 tmp.setPx(aPtVector.x()); tmp.setPy(a 1761 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1802 SumPx += aPtVector.x(); SumPy += aP 1762 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1803 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma 1763 Mt=std::sqrt(aPtVector.mag2()+sqr(Qmass));
1804 TargSumMt += Mt; 1764 TargSumMt += Mt;
1805 1765
1806 // Sampling of Z fraction 1766 // Sampling of Z fraction
1807 tmp.setPz(SampleX(Xmin, NumberOfUnsam 1767 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1808 SumZ += tmp.z(); 1768 SumZ += tmp.z();
1809 tmp.setPz((*i)->Get4Momentum().pz()*t 1769 tmp.setPz((*i)->Get4Momentum().pz()*tmp.pz());
>> 1770 SumZw+=tmp.pz();
1810 NumberOfUnsampledSeaQuarks--; 1771 NumberOfUnsampledSeaQuarks--;
1811 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1772 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1812 tmp.setE(sqr(Mt)); 1773 tmp.setE(sqr(Mt));
1813 aParton->Set4Momentum(tmp); 1774 aParton->Set4Momentum(tmp);
1814 1775
1815 aParton = (*i)->GetNextAntiParton(); 1776 aParton = (*i)->GetNextAntiParton(); // for anti-quarks
1816 #ifdef debugQGSParticipants 1777 #ifdef debugQGSParticipants
1817 G4cout<<" "<<aParton->GetDefinition 1778 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1818 G4cout<<" "<<tmp<<" "< << 1779 G4cout<<" "<<tmp<<" "<<SumZw<<" "<<SumZ<<G4endl;
1819 #endif 1780 #endif
1820 aPtVector = GaussianPt(SigPt, aHugeVa 1781 aPtVector = GaussianPt(SigPt, aHugeValue);
1821 tmp.setPx(aPtVector.x()); tmp.setPy(a 1782 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1822 SumPx += aPtVector.x(); SumPy += aP 1783 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1823 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma 1784 Mt=std::sqrt(aPtVector.mag2()+sqr(Qmass));
1824 TargSumMt += Mt; 1785 TargSumMt += Mt;
1825 1786
1826 // Sampling of Z fraction 1787 // Sampling of Z fraction
1827 tmp.setPz(SampleX(Xmin, NumberOfUnsam 1788 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1828 SumZ += tmp.z(); 1789 SumZ += tmp.z();
1829 tmp.setPz((*i)->Get4Momentum().pz()*t 1790 tmp.setPz((*i)->Get4Momentum().pz()*tmp.pz());
>> 1791 SumZw+=tmp.pz();
1830 NumberOfUnsampledSeaQuarks--; 1792 NumberOfUnsampledSeaQuarks--;
1831 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1793 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1832 tmp.setE(sqr(Mt)); 1794 tmp.setE(sqr(Mt));
1833 aParton->Set4Momentum(tmp); 1795 aParton->Set4Momentum(tmp);
1834 #ifdef debugQGSParticipants 1796 #ifdef debugQGSParticipants
1835 G4cout<<" "<<tmp<<" "< << 1797 G4cout<<" "<<tmp<<" "<<SumZw<<" "<<SumZ<<G4endl;
1836 #endif 1798 #endif
1837 } 1799 }
1838 1800
1839 // Valence quark 1801 // Valence quark
1840 aParton = (*i)->GetNextParton(); // f 1802 aParton = (*i)->GetNextParton(); // for quarks
1841 #ifdef debugQGSParticipants 1803 #ifdef debugQGSParticipants
1842 G4cout<<"Val quark of Tr"<<" "<<aPart 1804 G4cout<<"Val quark of Tr"<<" "<<aParton->GetDefinition()->GetParticleName();
1843 #endif 1805 #endif
1844 aPtVector = GaussianPt(SigPt, aHugeValu 1806 aPtVector = GaussianPt(SigPt, aHugeValue);
1845 tmp.setPx(aPtVector.x()); tmp.setPy(aPt 1807 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVector.y());
1846 SumPx += aPtVector.x(); SumPy += aPtV 1808 SumPx += aPtVector.x(); SumPy += aPtVector.y();
1847 Mt=std::sqrt(aPtVector.mag2()+sqr(VqM_t 1809 Mt=std::sqrt(aPtVector.mag2()+sqr(VqM_tr));
1848 TargSumMt += Mt; 1810 TargSumMt += Mt;
1849 1811
1850 // Sampling of Z fraction 1812 // Sampling of Z fraction
1851 tmp.setPz(SampleX(Xmin, NumberOfUnsampl 1813 tmp.setPz(SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta)*(1.0-SumZ));
1852 SumZ += tmp.z(); 1814 SumZ += tmp.z();
1853 tmp.setPz((*i)->Get4Momentum().pz()*tmp 1815 tmp.setPz((*i)->Get4Momentum().pz()*tmp.pz());
>> 1816 SumZw+=tmp.pz();
1854 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1817 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1855 tmp.setE(sqr(Mt)); 1818 tmp.setE(sqr(Mt));
1856 aParton->Set4Momentum(tmp); 1819 aParton->Set4Momentum(tmp);
1857 1820
1858 // Valence di-quark 1821 // Valence di-quark
1859 aParton = (*i)->GetNextAntiParton(); 1822 aParton = (*i)->GetNextAntiParton(); // for quarks
1860 #ifdef debugQGSParticipants 1823 #ifdef debugQGSParticipants
1861 G4cout<<" "<<aParton->GetDefinition() 1824 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1862 G4cout<<" "<<tmp<<" "<<S << 1825 G4cout<<" "<<tmp<<" "<<SumZw<<" (sum z-fracs) "<<SumZ<<" (total z-sum) "<<G4endl;
1863 #endif 1826 #endif
1864 tmp.setPx(-SumPx); tmp 1827 tmp.setPx(-SumPx); tmp.setPy(-SumPy);
1865 //Uzhi 2019 Mt=std::sqrt(aPtVector.mag 1828 //Uzhi 2019 Mt=std::sqrt(aPtVector.mag2()+sqr(VqqM_tr));
1866 Mt=std::sqrt( sqr(SumPx) + sqr(SumPy) + 1829 Mt=std::sqrt( sqr(SumPx) + sqr(SumPy) + sqr(VqqM_tr) ); //Uzhi 2019
1867 TargSumMt += Mt; 1830 TargSumMt += Mt;
1868 1831
1869 tmp.setPz((*i)->Get4Momentum().pz()*(1. 1832 tmp.setPz((*i)->Get4Momentum().pz()*(1.0 - SumZ));
>> 1833 SumZw+=tmp.pz();
1870 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1834 TargSumMt2perX +=sqr(Mt)/tmp.pz();
1871 tmp.setE(sqr(Mt)); 1835 tmp.setE(sqr(Mt));
1872 aParton->Set4Momentum(tmp); 1836 aParton->Set4Momentum(tmp);
1873 #ifdef debugQGSParticipants 1837 #ifdef debugQGSParticipants
1874 G4cout<<" "<<tmp<<" "<<1 << 1838 G4cout<<" "<<tmp<<" "<<SumZw<<" "<<1.0<<" "<<(*i)->Get4Momentum().pz()<<G4endl;
1875 #endif 1839 #endif
1876 1840
1877 } // End of for(i = theTargets.begin(); 1841 } // End of for(i = theTargets.begin(); i != theTargets.end(); i++ )
1878 1842
1879 if( ProjSumMt + TargSumMt > Sqr 1843 if( ProjSumMt + TargSumMt > SqrtS ) {
1880 Success = false; continue;} 1844 Success = false; continue;}
1881 if( std::sqrt(ProjSumMt2perX) + std::sqrt 1845 if( std::sqrt(ProjSumMt2perX) + std::sqrt(TargSumMt2perX) > SqrtS ) {
1882 Success = false; continue;} 1846 Success = false; continue;}
1883 1847
1884 } while( (!Success) && 1848 } while( (!Success) &&
1885 attempt < maxNumberOfAttempts ); 1849 attempt < maxNumberOfAttempts ); /* Loop checking, 07.08.2015, A.Ribon */
1886 1850
1887 if ( attempt >= maxNumberOfAttempts ) { 1851 if ( attempt >= maxNumberOfAttempts ) {
1888 return false; 1852 return false;
1889 } 1853 }
1890 1854
1891 //+++++++++++++++++++++++++++++++++++++++++ 1855 //+++++++++++++++++++++++++++++++++++++++++++
1892 1856
1893 G4double DecayMomentum2 = sqr(S) + sqr(Proj 1857 G4double DecayMomentum2 = sqr(S) + sqr(ProjSumMt2perX) + sqr(TargSumMt2perX)
1894 - 2.0*S*ProjSumMt2perX - 2.0*S 1858 - 2.0*S*ProjSumMt2perX - 2.0*S*TargSumMt2perX - 2.0*ProjSumMt2perX*TargSumMt2perX;
1895 1859
1896 G4double targetWminus=( S - ProjSumMt2perX 1860 G4double targetWminus=( S - ProjSumMt2perX + TargSumMt2perX + std::sqrt( DecayMomentum2 ))/2.0/SqrtS;
1897 G4double projectileWplus = SqrtS - TargSumM 1861 G4double projectileWplus = SqrtS - TargSumMt2perX/targetWminus;
1898 1862
1899 G4LorentzVector Tmp(0.,0.,0.,0.); 1863 G4LorentzVector Tmp(0.,0.,0.,0.);
1900 G4double z(0.); 1864 G4double z(0.);
1901 1865
1902 G4int nSeaPair = theProjectileSplitable->Ge 1866 G4int nSeaPair = theProjectileSplitable->GetSoftCollisionCount()-1;
1903 #ifdef debugQGSParticipants 1867 #ifdef debugQGSParticipants
1904 G4cout<<"Backward transformation ======== 1868 G4cout<<"Backward transformation ===================="<<G4endl;
1905 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4 1869 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4endl;
1906 #endif 1870 #endif
1907 1871
1908 for (G4int aSeaPair = 0; aSeaPair < nSeaPai 1872 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1909 { 1873 {
1910 aParton = theProjectileSplitable->GetNext 1874 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1911 #ifdef debugQGSParticipants 1875 #ifdef debugQGSParticipants
1912 G4cout<<"Sea quarks: "<<aSeaPair<<" "<< 1876 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1913 #endif 1877 #endif
1914 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1878 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1915 1879
1916 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e() 1880 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1917 Tmp.setE( projectileWplus*z/2.0 + Tmp.e() 1881 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1918 Tmp.transform( toLab ); 1882 Tmp.transform( toLab );
1919 1883
1920 aParton->Set4Momentum(Tmp); 1884 aParton->Set4Momentum(Tmp);
1921 1885
1922 aParton = theProjectileSplitable->GetNext 1886 aParton = theProjectileSplitable->GetNextAntiParton(); // for anti-quarks
1923 #ifdef debugQGSParticipants 1887 #ifdef debugQGSParticipants
1924 G4cout<<" "<<aParton->GetDefinition()-> 1888 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1925 G4cout<<" "<<Tmp<<" "<<Tmp 1889 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
1926 #endif 1890 #endif
1927 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1891 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1928 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e() 1892 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1929 Tmp.setE( projectileWplus*z/2.0 + Tmp.e() 1893 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1930 Tmp.transform( toLab ); 1894 Tmp.transform( toLab );
1931 1895
1932 aParton->Set4Momentum(Tmp); 1896 aParton->Set4Momentum(Tmp);
1933 #ifdef debugQGSParticipants 1897 #ifdef debugQGSParticipants
1934 G4cout<<" "<<Tmp<<" "<<Tmp 1898 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
1935 #endif 1899 #endif
1936 } 1900 }
1937 1901
1938 // For valence quark 1902 // For valence quark
1939 aParton = theProjectileSplitable->GetNextPa 1903 aParton = theProjectileSplitable->GetNextParton(); // for quarks
1940 #ifdef debugQGSParticipants 1904 #ifdef debugQGSParticipants
1941 G4cout<<"Val quark of Pr"<<" "<<aParton-> 1905 G4cout<<"Val quark of Pr"<<" "<<aParton->GetDefinition()->GetParticleName();
1942 #endif 1906 #endif
1943 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1907 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1944 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/( 1908 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1945 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/( 1909 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1946 Tmp.transform( toLab ); 1910 Tmp.transform( toLab );
1947 1911
1948 aParton->Set4Momentum(Tmp); 1912 aParton->Set4Momentum(Tmp);
1949 1913
1950 // For valence di-quark 1914 // For valence di-quark
1951 aParton = theProjectileSplitable->GetNextAn 1915 aParton = theProjectileSplitable->GetNextAntiParton();
1952 #ifdef debugQGSParticipants 1916 #ifdef debugQGSParticipants
1953 G4cout<<" "<<aParton->GetDefinition()->Ge 1917 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1954 G4cout<<" "<<Tmp<<" "<<Tmp.m 1918 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
1955 #endif 1919 #endif
1956 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1920 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1957 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/( 1921 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/(2.0*z*projectileWplus));
1958 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/( 1922 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/(2.0*z*projectileWplus));
1959 Tmp.transform( toLab ); 1923 Tmp.transform( toLab );
1960 1924
1961 aParton->Set4Momentum(Tmp); 1925 aParton->Set4Momentum(Tmp);
1962 1926
1963 #ifdef debugQGSParticipants 1927 #ifdef debugQGSParticipants
1964 G4cout<<" "<<Tmp<<" "<<Tmp.m 1928 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
1965 NuclNo=0; <<
1966 #endif 1929 #endif
1967 1930
1968 // End of work with the projectile 1931 // End of work with the projectile
1969 1932
1970 // Work with target nucleons 1933 // Work with target nucleons
>> 1934 NuclNo=0;
1971 for(i = theTargets.begin(); i != theTargets 1935 for(i = theTargets.begin(); i != theTargets.end(); i++ )
1972 { 1936 {
1973 nSeaPair = (*i)->GetSoftCollisionCount()- 1937 nSeaPair = (*i)->GetSoftCollisionCount()-1;
1974 #ifdef debugQGSParticipants 1938 #ifdef debugQGSParticipants
1975 G4cout<<"nSeaPair of target and N# "<<n 1939 G4cout<<"nSeaPair of target and N# "<<nSeaPair<<" "<<NuclNo<<G4endl;
1976 NuclNo++; <<
1977 #endif 1940 #endif
>> 1941 NuclNo++;
1978 for (G4int aSeaPair = 0; aSeaPair < nSeaP 1942 for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1979 { 1943 {
1980 aParton = (*i)->GetNextParton(); // f 1944 aParton = (*i)->GetNextParton(); // for quarks
1981 #ifdef debugQGSParticipants 1945 #ifdef debugQGSParticipants
1982 G4cout<<"Sea quarks: "<<aSeaPair<<" " 1946 G4cout<<"Sea quarks: "<<aSeaPair<<" "<<aParton->GetDefinition()->GetParticleName();
1983 #endif 1947 #endif
1984 Tmp =aParton->Get4Momentum(); z=Tmp.z() 1948 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1985 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e() 1949 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1986 Tmp.setE( targetWminus*z/2.0 + Tmp.e() 1950 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1987 Tmp.transform( toLab ); 1951 Tmp.transform( toLab );
1988 1952
1989 aParton->Set4Momentum(Tmp); 1953 aParton->Set4Momentum(Tmp);
1990 1954
1991 aParton = (*i)->GetNextAntiParton(); 1955 aParton = (*i)->GetNextAntiParton(); // for quarks
1992 #ifdef debugQGSParticipants 1956 #ifdef debugQGSParticipants
1993 G4cout<<" "<<aParton->GetDefinition() 1957 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
1994 G4cout<<" "<<Tmp<<" "<<T 1958 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
1995 #endif 1959 #endif
1996 Tmp =aParton->Get4Momentum(); z=Tmp.z() 1960 Tmp =aParton->Get4Momentum(); z=Tmp.z();
1997 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e() 1961 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1998 Tmp.setE( targetWminus*z/2.0 + Tmp.e() 1962 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
1999 Tmp.transform( toLab ); 1963 Tmp.transform( toLab );
2000 1964
2001 aParton->Set4Momentum(Tmp); 1965 aParton->Set4Momentum(Tmp);
2002 #ifdef debugQGSParticipants 1966 #ifdef debugQGSParticipants
2003 G4cout<<" "<<Tmp<<" "<<T 1967 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<G4endl;
2004 #endif 1968 #endif
2005 } 1969 }
2006 1970
2007 // Valence quark 1971 // Valence quark
2008 1972
2009 aParton = (*i)->GetNextParton(); // for 1973 aParton = (*i)->GetNextParton(); // for quarks
2010 #ifdef debugQGSParticipants 1974 #ifdef debugQGSParticipants
2011 G4cout<<"Val quark of Tr"<<" "<<aParton 1975 G4cout<<"Val quark of Tr"<<" "<<aParton->GetDefinition()->GetParticleName();
2012 #endif 1976 #endif
2013 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1977 Tmp =aParton->Get4Momentum(); z=Tmp.z();
2014 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/( 1978 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2015 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/( 1979 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2016 Tmp.transform( toLab ); 1980 Tmp.transform( toLab );
2017 1981
2018 aParton->Set4Momentum(Tmp); 1982 aParton->Set4Momentum(Tmp);
2019 1983
2020 // Valence di-quark 1984 // Valence di-quark
2021 aParton = (*i)->GetNextAntiParton(); // 1985 aParton = (*i)->GetNextAntiParton(); // for quarks
2022 #ifdef debugQGSParticipants 1986 #ifdef debugQGSParticipants
2023 G4cout<<" "<<aParton->GetDefinition()-> 1987 G4cout<<" "<<aParton->GetDefinition()->GetParticleName()<<G4endl;
2024 G4cout<<" "<<Tmp<<" "<<Tmp 1988 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
2025 #endif 1989 #endif
2026 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1990 Tmp =aParton->Get4Momentum(); z=Tmp.z();
2027 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/( 1991 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2028 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/( 1992 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/(2.0*z*targetWminus));
2029 Tmp.transform( toLab ); 1993 Tmp.transform( toLab );
2030 1994
2031 aParton->Set4Momentum(Tmp); 1995 aParton->Set4Momentum(Tmp);
2032 #ifdef debugQGSParticipants 1996 #ifdef debugQGSParticipants
2033 G4cout<<" "<<Tmp<<" "<<Tmp 1997 G4cout<<" "<<Tmp<<" "<<Tmp.mag()<<" (mass)"<<G4endl;
2034 NuclNo++; <<
2035 #endif 1998 #endif
>> 1999 NuclNo++;
2036 } // End of for(i = theTargets.begin(); i 2000 } // End of for(i = theTargets.begin(); i != theTargets.end(); i++ )
2037 2001
2038 return true; 2002 return true;
2039 } 2003 }
2040 2004
2041 //=========================================== 2005 //======================================================
2042 G4double G4QGSParticipants:: 2006 G4double G4QGSParticipants::
2043 SampleX(G4double, G4int nSea, G4int, G4double << 2007 SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
2044 { 2008 {
>> 2009 G4double Xmin=anXmin; G4int Nsea=totalSea; Xmin*=1.; Nsea++; // Must be erased
2045 G4double Oalfa = 1./(alpha + 1.); 2010 G4double Oalfa = 1./(alpha + 1.);
2046 G4double Obeta = 1./(aBeta + (alpha + 1.)*n 2011 G4double Obeta = 1./(aBeta + (alpha + 1.)*nSea + 1.); // ?
2047 2012
2048 G4double Ksi1, Ksi2, r1, r2, r12; 2013 G4double Ksi1, Ksi2, r1, r2, r12;
2049 const G4int maxNumberOfLoops = 1000; 2014 const G4int maxNumberOfLoops = 1000;
2050 G4int loopCounter = 0; 2015 G4int loopCounter = 0;
2051 do 2016 do
2052 { 2017 {
2053 Ksi1 = G4UniformRand(); r1 = G4Pow::GetIn 2018 Ksi1 = G4UniformRand(); r1 = G4Pow::GetInstance()->powA(Ksi1,Oalfa);
2054 Ksi2 = G4UniformRand(); r2 = G4Pow::GetIn 2019 Ksi2 = G4UniformRand(); r2 = G4Pow::GetInstance()->powA(Ksi2,Obeta);
2055 r12=r1+r2; 2020 r12=r1+r2;
2056 } while( ( r12 > 1.) && 2021 } while( ( r12 > 1.) &&
2057 ++loopCounter < maxNumberOfLoops ) 2022 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
2058 if ( loopCounter >= maxNumberOfLoops ) { 2023 if ( loopCounter >= maxNumberOfLoops ) {
2059 return 0.5; // Just an acceptable value, 2024 return 0.5; // Just an acceptable value, without any physics consideration.
2060 } 2025 }
2061 2026
2062 G4double result = r1/r12; 2027 G4double result = r1/r12;
2063 return result; 2028 return result;
2064 } 2029 }
2065 2030
2066 //=========================================== 2031 //======================================================
2067 void G4QGSParticipants::CreateStrings() 2032 void G4QGSParticipants::CreateStrings()
2068 { 2033 {
2069 2034
2070 #ifdef debugQGSParticipants 2035 #ifdef debugQGSParticipants
2071 G4cout<<"CreateStrings() ................ 2036 G4cout<<"CreateStrings() ..................."<<G4endl;
2072 #endif 2037 #endif
2073 2038
2074 if ( ! theProjectileSplitable ) { 2039 if ( ! theProjectileSplitable ) {
2075 #ifdef debugQGSParticipants 2040 #ifdef debugQGSParticipants
2076 G4cout<<"BAD situation: theProjectileSp 2041 G4cout<<"BAD situation: theProjectileSplitable is NULL ! Returning immediately"<<G4endl;
2077 #endif 2042 #endif
2078 return; 2043 return;
2079 } 2044 }
2080 2045
2081 #ifdef debugQGSParticipants 2046 #ifdef debugQGSParticipants
2082 G4cout<<"theProjectileSplitable->GetStatu 2047 G4cout<<"theProjectileSplitable->GetStatus() "<<theProjectileSplitable->GetStatus()<<G4endl;
2083 G4LorentzVector str4Mom; 2048 G4LorentzVector str4Mom;
2084 #endif 2049 #endif
2085 2050
2086 if( theProjectileSplitable->GetStatus() == 2051 if( theProjectileSplitable->GetStatus() == 1 ) // The projectile has participated only in diffr. inter.
2087 { 2052 {
2088 G4ThreeVector Position = theProjectileSpl 2053 G4ThreeVector Position = theProjectileSplitable->GetPosition();
2089 2054
2090 G4PartonPair * aPair = new G4PartonPair(t 2055 G4PartonPair * aPair = new G4PartonPair(theProjectileSplitable->GetNextParton(),
2091 t 2056 theProjectileSplitable->GetNextAntiParton(),
2092 G4PartonPair::DIFFRACTIVE 2057 G4PartonPair::DIFFRACTIVE, G4PartonPair::TARGET);
2093 #ifdef debugQGSParticipants 2058 #ifdef debugQGSParticipants
2094 G4cout << "Pr. Diffr. String: Qs 4mom X 2059 G4cout << "Pr. Diffr. String: Qs 4mom X " <<G4endl;
2095 G4cout << " " << aPair->Ge 2060 G4cout << " " << aPair->GetParton1()->GetPDGcode() << " "
2096 << aPair->GetParton1()->Get4Momentum 2061 << aPair->GetParton1()->Get4Momentum() << " "
2097 << aPair->GetParton1()->GetX() 2062 << aPair->GetParton1()->GetX() << " " << G4endl;
2098 G4cout << " " << aPair->Ge 2063 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2099 << aPair->GetParton2()->Get4Momentum 2064 << aPair->GetParton2()->Get4Momentum() << " "
2100 << aPair->GetParton2()->GetX() 2065 << aPair->GetParton2()->GetX() << " " << G4endl;
2101 str4Mom += aPair->GetParton1()->Get4Mom 2066 str4Mom += aPair->GetParton1()->Get4Momentum();
2102 str4Mom += aPair->GetParton2()->Get4Mom 2067 str4Mom += aPair->GetParton2()->Get4Momentum();
2103 #endif 2068 #endif
2104 2069
2105 thePartonPairs.push_back(aPair); 2070 thePartonPairs.push_back(aPair);
2106 } 2071 }
2107 2072
2108 G4int N_EnvTarg = NumberOfInvolvedNucleonsO 2073 G4int N_EnvTarg = NumberOfInvolvedNucleonsOfTarget;
2109 2074
2110 for ( G4int i = 0; i < N_EnvTarg; i++ ) { 2075 for ( G4int i = 0; i < N_EnvTarg; i++ ) {
2111 G4Nucleon* aTargetNucleon = TheInvolvedNu 2076 G4Nucleon* aTargetNucleon = TheInvolvedNucleonsOfTarget[ i ];
2112 2077
2113 G4VSplitableHadron* HitTargetNucleon = aT 2078 G4VSplitableHadron* HitTargetNucleon = aTargetNucleon->GetSplitableHadron();
2114 2079
2115 #ifdef debugQGSParticipants 2080 #ifdef debugQGSParticipants
2116 G4cout<<"Involved Nucleon # and its sta 2081 G4cout<<"Involved Nucleon # and its status "<<i<<" "<<HitTargetNucleon->GetStatus()<<G4endl;
2117 #endif 2082 #endif
2118 if( HitTargetNucleon->GetStatus() >= 1) 2083 if( HitTargetNucleon->GetStatus() >= 1) // Create diffractive string
2119 { 2084 {
2120 G4ThreeVector Position = HitTargetN 2085 G4ThreeVector Position = HitTargetNucleon->GetPosition();
2121 2086
2122 G4PartonPair * aPair = new G4PartonPair 2087 G4PartonPair * aPair = new G4PartonPair(HitTargetNucleon->GetNextParton(),
2123 2088 HitTargetNucleon->GetNextAntiParton(),
2124 G4PartonPair::DIFFRACTI 2089 G4PartonPair::DIFFRACTIVE, G4PartonPair::TARGET);
2125 #ifdef debugQGSParticipants 2090 #ifdef debugQGSParticipants
2126 G4cout << "Tr. Diffr. String: Qs 4mom 2091 G4cout << "Tr. Diffr. String: Qs 4mom X " <<G4endl;
2127 G4cout << "Diffr. String " << aPair->GetPar 2092 G4cout << "Diffr. String " << aPair->GetParton1()->GetPDGcode() << " "
2128 << aPair->GetParton1()->Get4Moment 2093 << aPair->GetParton1()->Get4Momentum() << " "
2129 << aPair->GetParton1()->GetX() 2094 << aPair->GetParton1()->GetX() << " " << G4endl;
2130 G4cout << " " << aPair->GetPar 2095 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2131 << aPair->GetParton2()->Get4Moment 2096 << aPair->GetParton2()->Get4Momentum() << " "
2132 << aPair->GetParton2()->GetX() 2097 << aPair->GetParton2()->GetX() << " " << G4endl;
2133 2098
2134 str4Mom += aPair->GetParton1()->Get4Momentu 2099 str4Mom += aPair->GetParton1()->Get4Momentum();
2135 str4Mom += aPair->GetParton2()->Get4Momentu 2100 str4Mom += aPair->GetParton2()->Get4Momentum();
2136 #endif 2101 #endif
2137 2102
2138 thePartonPairs.push_back(aPair); 2103 thePartonPairs.push_back(aPair);
2139 } 2104 }
2140 } 2105 }
2141 2106
2142 //----------------------------------------- 2107 //-----------------------------------------
2143 #ifdef debugQGSParticipants 2108 #ifdef debugQGSParticipants
2144 G4int IntNo=0; <<
2145 G4cout<<"Strings created in soft interact 2109 G4cout<<"Strings created in soft interactions"<<G4endl;
2146 #endif 2110 #endif
2147 std::vector<G4InteractionContent*>::iterato 2111 std::vector<G4InteractionContent*>::iterator i;
>> 2112 G4int IntNo=0;
2148 i = theInteractions.begin(); 2113 i = theInteractions.begin();
2149 while ( i != theInteractions.end() ) /* Lo 2114 while ( i != theInteractions.end() ) /* Loop checking, 07.08.2015, A.Ribon */
2150 { 2115 {
2151 G4InteractionContent* anIniteraction = *i 2116 G4InteractionContent* anIniteraction = *i;
2152 G4PartonPair * aPair = NULL; 2117 G4PartonPair * aPair = NULL;
2153 2118
2154 #ifdef debugQGSParticipants 2119 #ifdef debugQGSParticipants
2155 G4cout<<"An interaction # and soft coll 2120 G4cout<<"An interaction # and soft coll. # "<<IntNo<<" "
2156 <<anIniteraction->GetNumberOfSoft 2121 <<anIniteraction->GetNumberOfSoftCollisions()<<G4endl;
2157 IntNo++; <<
2158 #endif 2122 #endif
>> 2123 IntNo++;
2159 if (anIniteraction->GetNumberOfSoftCollis 2124 if (anIniteraction->GetNumberOfSoftCollisions())
2160 { 2125 {
2161 G4VSplitableHadron* pProjectile = anIni 2126 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
2162 G4VSplitableHadron* pTarget = anIni 2127 G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
2163 2128
2164 for (G4int j = 0; j < anIniteraction->G 2129 for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
2165 { 2130 {
2166 aPair = new G4PartonPair(pTarget->Get 2131 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
2167 G4PartonPair::SOFT, G4PartonPair::TA 2132 G4PartonPair::SOFT, G4PartonPair::TARGET);
2168 #ifdef debugQGSParticipants 2133 #ifdef debugQGSParticipants
2169 G4cout << "SoftPair " << aPair->GetParton 2134 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2170 << aPair->GetParton1()->Get4Momentum() < 2135 << aPair->GetParton1()->Get4Momentum() << " "
2171 << aPair->GetParton1()->Get4Momentum().m 2136 << aPair->GetParton1()->Get4Momentum().mag()<<G4endl;
2172 G4cout << " " << aPair->GetParton 2137 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2173 << aPair->GetParton2()->Get4Momentum() < 2138 << aPair->GetParton2()->Get4Momentum() << " "
2174 <<aPair->GetParton2()->Get4Momentum().m 2139 <<aPair->GetParton2()->Get4Momentum().mag()<<G4endl;
2175 str4Mom += aPair->GetParton1()->Get4Momen 2140 str4Mom += aPair->GetParton1()->Get4Momentum();
2176 str4Mom += aPair->GetParton2()->Get4Momen 2141 str4Mom += aPair->GetParton2()->Get4Momentum();
2177 #endif 2142 #endif
2178 2143
2179 thePartonPairs.push_back(aPair); 2144 thePartonPairs.push_back(aPair);
2180 2145
2181 aPair = new G4PartonPair(pProjectile->GetNe 2146 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
2182 G4PartonPair::SOFT, G4PartonPair::PR 2147 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
2183 #ifdef debugQGSParticipants 2148 #ifdef debugQGSParticipants
2184 G4cout << "SoftPair " << aPair->GetParton 2149 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2185 << aPair->GetParton1()->Get4Momentum() < 2150 << aPair->GetParton1()->Get4Momentum() << " "
2186 << aPair->GetParton1()->Get4Moment 2151 << aPair->GetParton1()->Get4Momentum().mag()<<G4endl;
2187 G4cout << " " << aPair->GetParton 2152 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2188 << aPair->GetParton2()->Get4Momentum() < 2153 << aPair->GetParton2()->Get4Momentum() << " "
2189 << aPair->GetParton2()->Get4Momentum().m 2154 << aPair->GetParton2()->Get4Momentum().mag()<<G4endl;
2190 #endif 2155 #endif
2191 #ifdef debugQGSParticipants 2156 #ifdef debugQGSParticipants
2192 str4Mom += aPair->GetParton1()->Get4Momen 2157 str4Mom += aPair->GetParton1()->Get4Momentum();
2193 str4Mom += aPair->GetParton2()->Get4Momen 2158 str4Mom += aPair->GetParton2()->Get4Momentum();
2194 #endif 2159 #endif
2195 2160
2196 thePartonPairs.push_back(aPair); 2161 thePartonPairs.push_back(aPair);
2197 } 2162 }
2198 2163
2199 delete *i; 2164 delete *i;
2200 i=theInteractions.erase(i); // i now 2165 i=theInteractions.erase(i); // i now points to the next interaction
2201 } 2166 }
2202 else 2167 else
2203 { 2168 {
2204 i++; 2169 i++;
2205 } 2170 }
2206 } // End of while ( i != theInteractions.e 2171 } // End of while ( i != theInteractions.end() )
2207 #ifdef debugQGSParticipants 2172 #ifdef debugQGSParticipants
2208 G4cout << "Sum of strings 4 momenta " << 2173 G4cout << "Sum of strings 4 momenta " << str4Mom << G4endl<<G4endl;
2209 #endif 2174 #endif
2210 } 2175 }
2211 2176
2212 //=========================================== 2177 //============================================================================
2213 2178
2214 void G4QGSParticipants::GetResiduals() { 2179 void G4QGSParticipants::GetResiduals() {
2215 // This method is needed for the correct ap 2180 // This method is needed for the correct application of G4PrecompoundModelInterface
2216 2181
2217 #ifdef debugQGSParticipants 2182 #ifdef debugQGSParticipants
2218 G4cout << "GetResiduals(): GetProjectileN 2183 G4cout << "GetResiduals(): GetProjectileNucleus()? " << GetProjectileNucleus() << G4endl;
2219 #endif 2184 #endif
2220 2185
2221 #ifdef debugQGSParticipants 2186 #ifdef debugQGSParticipants
2222 G4cout << "NumberOfInvolvedNucleonsOfTarg 2187 G4cout << "NumberOfInvolvedNucleonsOfTarget "<< NumberOfInvolvedNucleonsOfTarget << G4endl;
2223 #endif 2188 #endif
2224 2189
2225 G4double DeltaExcitationE = TargetResidualE 2190 G4double DeltaExcitationE = TargetResidualExcitationEnergy / G4double( NumberOfInvolvedNucleonsOfTarget );
2226 G4LorentzVector DeltaPResidualNucleus = Tar 2191 G4LorentzVector DeltaPResidualNucleus = TargetResidual4Momentum /
2227 G4d 2192 G4double( NumberOfInvolvedNucleonsOfTarget );
2228 2193
2229 for ( G4int i = 0; i < NumberOfInvolvedNucl 2194 for ( G4int i = 0; i < NumberOfInvolvedNucleonsOfTarget; i++ ) {
2230 G4Nucleon* aNucleon = TheInvolvedNucleons 2195 G4Nucleon* aNucleon = TheInvolvedNucleonsOfTarget[i];
2231 2196
2232 #ifdef debugQGSParticipants 2197 #ifdef debugQGSParticipants
2233 G4VSplitableHadron* targetSplitable = a 2198 G4VSplitableHadron* targetSplitable = aNucleon->GetSplitableHadron();
2234 G4cout << i << " Hit? " << aNucleon->Ar 2199 G4cout << i << " Hit? " << aNucleon->AreYouHit() << " " << targetSplitable << G4endl;
2235 if ( targetSplitable ) G4cout << i << " 2200 if ( targetSplitable ) G4cout << i << "Status " << targetSplitable->GetStatus() << G4endl;
2236 #endif 2201 #endif
2237 2202
2238 G4LorentzVector tmp = -DeltaPResidualNucl 2203 G4LorentzVector tmp = -DeltaPResidualNucleus;
2239 aNucleon->SetMomentum( tmp ); 2204 aNucleon->SetMomentum( tmp );
2240 aNucleon->SetBindingEnergy( DeltaExcitati 2205 aNucleon->SetBindingEnergy( DeltaExcitationE );
2241 } 2206 }
2242 2207
2243 //------------------------------------- 2208 //-------------------------------------
2244 if( TargetResidualMassNumber != 0 ) 2209 if( TargetResidualMassNumber != 0 )
2245 { 2210 {
2246 G4ThreeVector bstToCM =TargetResidual4Mom 2211 G4ThreeVector bstToCM =TargetResidual4Momentum.findBoostToCM();
2247 2212
2248 G4V3DNucleus* theTargetNucleus = GetTarge 2213 G4V3DNucleus* theTargetNucleus = GetTargetNucleus();
2249 G4LorentzVector residualMomentum(0.,0.,0. 2214 G4LorentzVector residualMomentum(0.,0.,0.,0.);
2250 G4Nucleon* aNucleon = 0; 2215 G4Nucleon* aNucleon = 0;
2251 theTargetNucleus->StartLoop(); 2216 theTargetNucleus->StartLoop();
2252 while ( ( aNucleon = theTargetNucleus->Ge 2217 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2253 if ( !aNucleon->AreYouHit() ) { 2218 if ( !aNucleon->AreYouHit() ) {
2254 G4LorentzVector tmp=aNucleon->Get4Mom 2219 G4LorentzVector tmp=aNucleon->Get4Momentum(); tmp.boost(bstToCM);
2255 aNucleon->SetMomentum(tmp); 2220 aNucleon->SetMomentum(tmp);
2256 residualMomentum +=tmp; 2221 residualMomentum +=tmp;
2257 } 2222 }
2258 } 2223 }
2259 2224
2260 residualMomentum/=TargetResidualMassNumbe 2225 residualMomentum/=TargetResidualMassNumber;
2261 2226
2262 G4double Mass = TargetResidual4Momentum.m 2227 G4double Mass = TargetResidual4Momentum.mag();
2263 G4double SumMasses=0.; 2228 G4double SumMasses=0.;
2264 2229
2265 aNucleon = 0; 2230 aNucleon = 0;
2266 theTargetNucleus->StartLoop(); 2231 theTargetNucleus->StartLoop();
2267 while ( ( aNucleon = theTargetNucleus->Ge 2232 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2268 if ( !aNucleon->AreYouHit() ) { 2233 if ( !aNucleon->AreYouHit() ) {
2269 G4LorentzVector tmp=aNucleon->Get4Mom 2234 G4LorentzVector tmp=aNucleon->Get4Momentum() - residualMomentum;
2270 G4double E=std::sqrt(tmp.vect().mag2( 2235 G4double E=std::sqrt(tmp.vect().mag2()+
2271 sqr(aNucleon->Ge 2236 sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy()));
2272 tmp.setE(E); aNucleon->SetMomentum(t 2237 tmp.setE(E); aNucleon->SetMomentum(tmp);
2273 SumMasses+=E; 2238 SumMasses+=E;
2274 } 2239 }
2275 } 2240 }
2276 2241
2277 G4double Chigh=Mass/SumMasses; G4double C 2242 G4double Chigh=Mass/SumMasses; G4double Clow=0; G4double C;
2278 const G4int maxNumberOfLoops = 1000; 2243 const G4int maxNumberOfLoops = 1000;
2279 G4int loopCounter = 0; 2244 G4int loopCounter = 0;
2280 do 2245 do
2281 { 2246 {
2282 C=(Chigh+Clow)/2.; 2247 C=(Chigh+Clow)/2.;
2283 2248
2284 SumMasses=0.; 2249 SumMasses=0.;
2285 aNucleon = 0; 2250 aNucleon = 0;
2286 theTargetNucleus->StartLoop(); 2251 theTargetNucleus->StartLoop();
2287 while ( ( aNucleon = theTargetNucleus-> 2252 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2288 if ( !aNucleon->AreYouHit() ) { 2253 if ( !aNucleon->AreYouHit() ) {
2289 G4LorentzVector tmp=aNucleon->Get4M 2254 G4LorentzVector tmp=aNucleon->Get4Momentum();
2290 G4double E=std::sqrt(tmp.vect().mag 2255 G4double E=std::sqrt(tmp.vect().mag2()*sqr(C)+
2291 sqr(aNucleon-> 2256 sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy()));
2292 SumMasses+=E; 2257 SumMasses+=E;
2293 } 2258 }
2294 } 2259 }
2295 2260
2296 if(SumMasses > Mass) {Chigh=C;} 2261 if(SumMasses > Mass) {Chigh=C;}
2297 else {Clow =C;} 2262 else {Clow =C;}
2298 2263
2299 } while( (Chigh-Clow > 0.01) && 2264 } while( (Chigh-Clow > 0.01) &&
2300 ++loopCounter < maxNumberOfLoops 2265 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
2301 if ( loopCounter >= maxNumberOfLoops ) { 2266 if ( loopCounter >= maxNumberOfLoops ) {
2302 #ifdef debugQGSParticipants 2267 #ifdef debugQGSParticipants
2303 G4cout <<"BAD situation: forced loop 2268 G4cout <<"BAD situation: forced loop exit!" << G4endl;
2304 #endif 2269 #endif
2305 // Perhaps there is something to set he 2270 // Perhaps there is something to set here...
2306 } else { 2271 } else {
2307 aNucleon = 0; 2272 aNucleon = 0;
2308 theTargetNucleus->StartLoop(); 2273 theTargetNucleus->StartLoop();
2309 while ( ( aNucleon = theTargetNucleus-> 2274 while ( ( aNucleon = theTargetNucleus->GetNextNucleon() ) ) { /* Loop checking, 07.08.2015, A.Ribon */
2310 if ( !aNucleon->AreYouHit() ) { 2275 if ( !aNucleon->AreYouHit() ) {
2311 G4LorentzVector tmp=aNucleon->Get4M 2276 G4LorentzVector tmp=aNucleon->Get4Momentum()*C;
2312 G4double E=std::sqrt(tmp.vect().mag 2277 G4double E=std::sqrt(tmp.vect().mag2()+
2313 sqr(aNucleon-> 2278 sqr(aNucleon->GetDefinition()->GetPDGMass()-aNucleon->GetBindingEnergy()));
2314 tmp.setE(E); tmp.boost(-bstToCM); 2279 tmp.setE(E); tmp.boost(-bstToCM);
2315 aNucleon->SetMomentum(tmp); 2280 aNucleon->SetMomentum(tmp);
2316 } 2281 }
2317 } 2282 }
2318 } 2283 }
2319 2284
2320 } // End of if( TargetResidualMassNumber ! 2285 } // End of if( TargetResidualMassNumber != 0 )
2321 //------------------------------------- 2286 //-------------------------------------
2322 2287
2323 #ifdef debugQGSParticipants 2288 #ifdef debugQGSParticipants
2324 G4cout << "End GetResiduals ------------- 2289 G4cout << "End GetResiduals -----------------" << G4endl;
2325 #endif 2290 #endif
2326 2291
2327 } 2292 }
2328 2293
2329 2294
2330 //=========================================== 2295 //======================================================
2331 void G4QGSParticipants::PerformSoftCollisions 2296 void G4QGSParticipants::PerformSoftCollisions() // It is not used
2332 { 2297 {
2333 std::vector<G4InteractionContent*>::iterato 2298 std::vector<G4InteractionContent*>::iterator i;
2334 G4LorentzVector str4Mom; 2299 G4LorentzVector str4Mom;
2335 i = theInteractions.begin(); 2300 i = theInteractions.begin();
2336 while ( i != theInteractions.end() ) /* Lo 2301 while ( i != theInteractions.end() ) /* Loop checking, 07.08.2015, A.Ribon */
2337 { 2302 {
2338 G4InteractionContent* anIniteraction = *i 2303 G4InteractionContent* anIniteraction = *i;
2339 G4PartonPair * aPair = NULL; 2304 G4PartonPair * aPair = NULL;
2340 if (anIniteraction->GetNumberOfSoftCollis 2305 if (anIniteraction->GetNumberOfSoftCollisions())
2341 { 2306 {
2342 G4VSplitableHadron* pProjectile = anIni 2307 G4VSplitableHadron* pProjectile = anIniteraction->GetProjectile();
2343 G4VSplitableHadron* pTarget = anIni 2308 G4VSplitableHadron* pTarget = anIniteraction->GetTarget();
2344 for (G4int j = 0; j < anIniteraction->G 2309 for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
2345 { 2310 {
2346 aPair = new G4PartonPair(pTarget->Get 2311 aPair = new G4PartonPair(pTarget->GetNextParton(), pProjectile->GetNextAntiParton(),
2347 G4PartonPair::SOFT, G4PartonPair::TA 2312 G4PartonPair::SOFT, G4PartonPair::TARGET);
2348 #ifdef debugQGSParticipants 2313 #ifdef debugQGSParticipants
2349 G4cout << "SoftPair " << aPair->GetParton 2314 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2350 << aPair->GetParton1()->Get4Momentum() < 2315 << aPair->GetParton1()->Get4Momentum() << " "
2351 << aPair->GetParton1()->GetX() << " " << 2316 << aPair->GetParton1()->GetX() << " " << G4endl;
2352 G4cout << " " << aPair->GetParton 2317 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2353 << aPair->GetParton2()->Get4Momentum() < 2318 << aPair->GetParton2()->Get4Momentum() << " "
2354 << aPair->GetParton2()->GetX() << " " << 2319 << aPair->GetParton2()->GetX() << " " << G4endl;
2355 #endif 2320 #endif
2356 #ifdef debugQGSParticipants 2321 #ifdef debugQGSParticipants
2357 str4Mom += aPair->GetParton1()->Get4Momen 2322 str4Mom += aPair->GetParton1()->Get4Momentum();
2358 str4Mom += aPair->GetParton2()->Get4Momen 2323 str4Mom += aPair->GetParton2()->Get4Momentum();
2359 #endif 2324 #endif
2360 thePartonPairs.push_back(aPair); 2325 thePartonPairs.push_back(aPair);
2361 aPair = new G4PartonPair(pProjectile->GetNe 2326 aPair = new G4PartonPair(pProjectile->GetNextParton(), pTarget->GetNextAntiParton(),
2362 G4PartonPair::SOFT, G4PartonPair::PR 2327 G4PartonPair::SOFT, G4PartonPair::PROJECTILE);
2363 #ifdef debugQGSParticipants 2328 #ifdef debugQGSParticipants
2364 G4cout << "SoftPair " << aPair->GetParton 2329 G4cout << "SoftPair " << aPair->GetParton1()->GetPDGcode() << " "
2365 << aPair->GetParton1()->Get4Momentum() < 2330 << aPair->GetParton1()->Get4Momentum() << " "
2366 << aPair->GetParton1()->GetX() << " " << 2331 << aPair->GetParton1()->GetX() << " " << G4endl;
2367 G4cout << " " << aPair->GetParton 2332 G4cout << " " << aPair->GetParton2()->GetPDGcode() << " "
2368 << aPair->GetParton2()->Get4Momentum() < 2333 << aPair->GetParton2()->Get4Momentum() << " "
2369 << aPair->GetParton2()->GetX() << " " << 2334 << aPair->GetParton2()->GetX() << " " << G4endl;
2370 #endif 2335 #endif
2371 #ifdef debugQGSParticipants 2336 #ifdef debugQGSParticipants
2372 str4Mom += aPair->GetParton1()->Get4Momen 2337 str4Mom += aPair->GetParton1()->Get4Momentum();
2373 str4Mom += aPair->GetParton2()->Get4Momen 2338 str4Mom += aPair->GetParton2()->Get4Momentum();
2374 #endif 2339 #endif
2375 thePartonPairs.push_back(aPair); 2340 thePartonPairs.push_back(aPair);
2376 } 2341 }
2377 delete *i; 2342 delete *i;
2378 i=theInteractions.erase(i); // i now 2343 i=theInteractions.erase(i); // i now points to the next interaction
2379 } else { 2344 } else {
2380 i++; 2345 i++;
2381 } 2346 }
2382 } 2347 }
2383 #ifdef debugQGSParticipants 2348 #ifdef debugQGSParticipants
2384 G4cout << " string 4 mom " << str4Mom << 2349 G4cout << " string 4 mom " << str4Mom << G4endl;
2385 #endif 2350 #endif
2386 } 2351 }
2387 2352
2388 2353
2389 //******************************************* 2354 //************************************************
2390 G4VSplitableHadron* G4QGSParticipants::Select 2355 G4VSplitableHadron* G4QGSParticipants::SelectInteractions(const G4ReactionProduct &thePrimary)
2391 { 2356 {
2392 // Check reaction threshold - goes to Chec 2357 // Check reaction threshold - goes to CheckThreshold
2393 2358
2394 theProjectileSplitable = new G4QGSMSplitabl 2359 theProjectileSplitable = new G4QGSMSplitableHadron(thePrimary, TRUE);
2395 theProjectileSplitable->SetStatus(1); 2360 theProjectileSplitable->SetStatus(1);
2396 2361
2397 G4LorentzVector aPrimaryMomentum(thePrimary 2362 G4LorentzVector aPrimaryMomentum(thePrimary.GetMomentum(), thePrimary.GetTotalEnergy());
2398 G4LorentzVector aTargetNMomentum(0.,0.,0.,9 2363 G4LorentzVector aTargetNMomentum(0.,0.,0.,938.);
2399 2364
2400 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPri 2365 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPrimaryMomentum.e()<1)) )
2401 { 2366 {
2402 throw G4HadronicException(__FILE__, __LIN 2367 throw G4HadronicException(__FILE__, __LINE__,
2403 "G4GammaParticipants::SelectInter 2368 "G4GammaParticipants::SelectInteractions: primary nan energy.");
2404 } 2369 }
2405 G4double S = (aPrimaryMomentum + aTargetNMo 2370 G4double S = (aPrimaryMomentum + aTargetNMomentum).mag2();
2406 G4double ThresholdMass = thePrimary.GetMass 2371 G4double ThresholdMass = thePrimary.GetMass() + 938.;
2407 ModelMode = SOFT; 2372 ModelMode = SOFT;
2408 2373
2409 #ifdef debugQGSParticipants 2374 #ifdef debugQGSParticipants
2410 G4cout << "Gamma projectile - SelectInter 2375 G4cout << "Gamma projectile - SelectInteractions " << G4endl;
2411 G4cout<<"Energy and Nucleus "<<thePrimary 2376 G4cout<<"Energy and Nucleus "<<thePrimary.GetTotalEnergy()<<" "<<theNucleus->GetMassNumber()<<G4endl;
2412 G4cout << "SqrtS ThresholdMass ModelMode 2377 G4cout << "SqrtS ThresholdMass ModelMode " <<std::sqrt(S)<<" "<<ThresholdMass<<" "<<ModelMode<< G4endl;
2413 #endif 2378 #endif
2414 2379
2415 if (sqr(ThresholdMass + ThresholdParameter) 2380 if (sqr(ThresholdMass + ThresholdParameter) > S)
2416 { 2381 {
2417 ModelMode = DIFFRACTIVE; 2382 ModelMode = DIFFRACTIVE;
2418 } 2383 }
2419 if (sqr(ThresholdMass + QGSMThreshold) > S) 2384 if (sqr(ThresholdMass + QGSMThreshold) > S) // thus only diffractive in cascade!
2420 { 2385 {
2421 ModelMode = DIFFRACTIVE; 2386 ModelMode = DIFFRACTIVE;
2422 } 2387 }
2423 2388
2424 // first find the collisions HPW 2389 // first find the collisions HPW
2425 std::for_each(theInteractions.begin(), theI 2390 std::for_each(theInteractions.begin(), theInteractions.end(), DeleteInteractionContent());
2426 theInteractions.clear(); 2391 theInteractions.clear();
>> 2392 G4int totalCuts = 0;
2427 2393
2428 G4int theCurrent = G4int(theNucleus->GetMas 2394 G4int theCurrent = G4int(theNucleus->GetMassNumber()*G4UniformRand());
2429 G4int NucleonNo=0; 2395 G4int NucleonNo=0;
2430 2396
2431 theNucleus->StartLoop(); 2397 theNucleus->StartLoop();
2432 G4Nucleon * pNucleon = 0; 2398 G4Nucleon * pNucleon = 0;
2433 2399
2434 while( (pNucleon = theNucleus->GetNextNucle 2400 while( (pNucleon = theNucleus->GetNextNucleon()) ) /* Loop checking, 07.08.2015, A.Ribon */
2435 { if(NucleonNo == theCurrent) break; Nucleo 2401 { if(NucleonNo == theCurrent) break; NucleonNo++; }
2436 2402
2437 if ( pNucleon ) { 2403 if ( pNucleon ) {
2438 G4QGSMSplitableHadron* aTarget = new G4QG 2404 G4QGSMSplitableHadron* aTarget = new G4QGSMSplitableHadron(*pNucleon);
2439 2405
2440 pNucleon->Hit(aTarget); 2406 pNucleon->Hit(aTarget);
2441 2407
2442 if ( (0.06 > G4UniformRand() &&(ModelMode 2408 if ( (0.06 > G4UniformRand() &&(ModelMode==SOFT)) || (ModelMode==DIFFRACTIVE ) )
2443 { 2409 {
2444 G4InteractionContent * aInteraction = n 2410 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
2445 theProjectileSplitable->SetStatus(1*the 2411 theProjectileSplitable->SetStatus(1*theProjectileSplitable->GetStatus());
2446 2412
2447 aInteraction->SetTarget(aTarget); 2413 aInteraction->SetTarget(aTarget);
2448 aInteraction->SetTargetNucleon(pNucleon 2414 aInteraction->SetTargetNucleon(pNucleon);
2449 aTarget->SetCollisionCount(0); 2415 aTarget->SetCollisionCount(0);
2450 aTarget->SetStatus(1); 2416 aTarget->SetStatus(1);
2451 2417
2452 aInteraction->SetNumberOfDiffractiveCol 2418 aInteraction->SetNumberOfDiffractiveCollisions(1);
2453 aInteraction->SetNumberOfSoftCollisions 2419 aInteraction->SetNumberOfSoftCollisions(0);
2454 aInteraction->SetStatus(1); 2420 aInteraction->SetStatus(1);
2455 2421
2456 theInteractions.push_back(aInteraction) 2422 theInteractions.push_back(aInteraction);
>> 2423 totalCuts += 1;
2457 } 2424 }
2458 else 2425 else
2459 { 2426 {
2460 // nondiffractive soft interaction occu 2427 // nondiffractive soft interaction occurs
2461 aTarget->IncrementCollisionCount(1); 2428 aTarget->IncrementCollisionCount(1);
2462 aTarget->SetStatus(0); 2429 aTarget->SetStatus(0);
2463 theTargets.push_back(aTarget); 2430 theTargets.push_back(aTarget);
2464 2431
2465 theProjectileSplitable->IncrementCollis 2432 theProjectileSplitable->IncrementCollisionCount(1);
2466 theProjectileSplitable->SetStatus(0*the 2433 theProjectileSplitable->SetStatus(0*theProjectileSplitable->GetStatus());
2467 2434
2468 G4InteractionContent * aInteraction = n 2435 G4InteractionContent * aInteraction = new G4InteractionContent(theProjectileSplitable);
2469 aInteraction->SetTarget(aTarget); 2436 aInteraction->SetTarget(aTarget);
2470 aInteraction->SetTargetNucleon(pNucleon 2437 aInteraction->SetTargetNucleon(pNucleon);
2471 aInteraction->SetNumberOfSoftCollisions 2438 aInteraction->SetNumberOfSoftCollisions(1);
2472 aInteraction->SetStatus(3); 2439 aInteraction->SetStatus(3);
2473 theInteractions.push_back(aInteraction) 2440 theInteractions.push_back(aInteraction);
>> 2441 totalCuts += 1;
2474 } 2442 }
2475 } 2443 }
2476 return theProjectileSplitable; 2444 return theProjectileSplitable;
2477 } 2445 }
2478 2446