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