| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // ------------------------------------------------------------
28 // GEANT 4 class implemetation file 28 // GEANT 4 class implemetation file
29 // 29 //
30 // ---------------- G4QuarkExchange ----- 30 // ---------------- G4QuarkExchange --------------
31 // by V. Uzhinsky, October 2016. 31 // by V. Uzhinsky, October 2016.
32 // QuarkExchange is used by strings mode 32 // QuarkExchange is used by strings models.
33 // Take a projectile and a target. 33 // Take a projectile and a target.
34 //Simulate Q exchange with excitation of proje 34 //Simulate Q exchange with excitation of projectile or target.
35 // ------------------------------------------- 35 // ------------------------------------------------------------
36 36
37 #include "G4QuarkExchange.hh" 37 #include "G4QuarkExchange.hh"
38 #include "globals.hh" 38 #include "globals.hh"
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 #include "Randomize.hh" 41 #include "Randomize.hh"
42 #include "G4LorentzRotation.hh" 42 #include "G4LorentzRotation.hh"
43 #include "G4ThreeVector.hh" 43 #include "G4ThreeVector.hh"
44 #include "G4ParticleDefinition.hh" 44 #include "G4ParticleDefinition.hh"
45 #include "G4VSplitableHadron.hh" 45 #include "G4VSplitableHadron.hh"
46 #include "G4ExcitedString.hh" 46 #include "G4ExcitedString.hh"
47 47
48 #include "G4ParticleTable.hh" // Uzhi June 20 <<
49 <<
50 #include "G4Log.hh" 48 #include "G4Log.hh"
51 #include "G4Pow.hh" 49 #include "G4Pow.hh"
52 50
53 //#define debugQuarkExchange 51 //#define debugQuarkExchange
54 52
55 G4QuarkExchange::G4QuarkExchange() << 53 G4QuarkExchange::G4QuarkExchange(){}
56 { <<
57 StrangeSuppress = (1.0-0.04)/2.0; // Uzhi J <<
58 // <<
59 } <<
60 54
61 G4QuarkExchange::~G4QuarkExchange(){} 55 G4QuarkExchange::~G4QuarkExchange(){}
62 56
63 G4bool G4QuarkExchange:: 57 G4bool G4QuarkExchange::
64 ExciteParticipants(G4VSplitableHadron *project 58 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
65 { 59 {
66 #ifdef debugQuarkExchange 60 #ifdef debugQuarkExchange
67 G4cout<<G4endl<<"G4QuarkExchange::ExcitePa 61 G4cout<<G4endl<<"G4QuarkExchange::ExciteParticipants"<<G4endl;
68 #endif 62 #endif
69 63
70 G4LorentzVector Pprojectile = projectile->Ge 64 G4LorentzVector Pprojectile = projectile->Get4Momentum();
71 G4double Mprojectile = projectile->Ge 65 G4double Mprojectile = projectile->GetDefinition()->GetPDGMass();
72 G4double Mprojectile2 = sqr(Mprojectil 66 G4double Mprojectile2 = sqr(Mprojectile);
73 67
74 G4LorentzVector Ptarget = target->Get4Moment 68 G4LorentzVector Ptarget = target->Get4Momentum();
75 G4double Mtarget = target->GetDefinit 69 G4double Mtarget = target->GetDefinition()->GetPDGMass();
76 G4double Mtarget2 = sqr(Mtarget); 70 G4double Mtarget2 = sqr(Mtarget);
77 71
78 #ifdef debugQuarkExchange 72 #ifdef debugQuarkExchange
79 G4cout<<"Proj Targ "<<projectile->GetDefin 73 G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl;
80 G4cout<<"Proj. 4-Mom "<<Pprojectile<<" "<< 74 G4cout<<"Proj. 4-Mom "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl
81 <<"Targ. 4-Mom "<<Ptarget <<" "<< 75 <<"Targ. 4-Mom "<<Ptarget <<" "<<Ptarget.mag() <<G4endl;
82 #endif 76 #endif
83 77
84 G4LorentzVector Psum=Pprojectile+Ptarget; 78 G4LorentzVector Psum=Pprojectile+Ptarget;
85 G4double SqrtS=Psum.mag(); 79 G4double SqrtS=Psum.mag();
86 G4double S =Psum.mag2(); 80 G4double S =Psum.mag2();
87 81
88 #ifdef debugQuarkExchange 82 #ifdef debugQuarkExchange
89 G4cout<<"SS Mpr Mtr SqrtS-Mprojectile-Mtar 83 G4cout<<"SS Mpr Mtr SqrtS-Mprojectile-Mtarget "<<SqrtS<<" "<<Mprojectile<<" "<<Mtarget
90 <<" "<<SqrtS-Mprojectile-Mtarget<<G4 84 <<" "<<SqrtS-Mprojectile-Mtarget<<G4endl;
91 #endif 85 #endif
92 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) 86 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) {
93 #ifdef debugQuarkExchange 87 #ifdef debugQuarkExchange
94 G4cerr<<"Energy is too small for quark e 88 G4cerr<<"Energy is too small for quark exchange!"<<G4endl;
95 G4cerr<<"Projectile: "<<projectile->GetD 89 G4cerr<<"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<" "
96 <<Pprojectile<<" "<<Pprojectile.ma 90 <<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
97 G4cerr<<"Target: "<<target->GetDefin 91 G4cerr<<"Target: "<<target->GetDefinition()->GetPDGEncoding()<<" "
98 <<Ptarget<<" "<<Ptarget.mag()<<G4e 92 <<Ptarget<<" "<<Ptarget.mag()<<G4endl;
99 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = 93 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "<<Pprojectile.mag()+Ptarget.mag()<<G4endl;
100 #endif 94 #endif
101 return true; 95 return true;
102 } 96 }
103 97
104 G4LorentzRotation toCms(-1*Psum.boostVector( 98 G4LorentzRotation toCms(-1*Psum.boostVector());
105 99
106 G4LorentzVector Ptmp=toCms*Pprojectile; 100 G4LorentzVector Ptmp=toCms*Pprojectile;
107 101
108 if ( Ptmp.pz() <= 0. ) 102 if ( Ptmp.pz() <= 0. )
109 { 103 {
110 // "String" moving backwards in CMS, abor 104 // "String" moving backwards in CMS, abort collision !!
111 // G4cout << " abort Collision!! " 105 // G4cout << " abort Collision!! " << G4endl;
112 return false; 106 return false;
113 } 107 }
114 108
115 toCms.rotateZ(-1*Ptmp.phi()); 109 toCms.rotateZ(-1*Ptmp.phi());
116 toCms.rotateY(-1*Ptmp.theta()); 110 toCms.rotateY(-1*Ptmp.theta());
117 111
118 G4LorentzRotation toLab(toCms.inverse()); 112 G4LorentzRotation toLab(toCms.inverse());
119 113
120 Pprojectile.transform(toCms); 114 Pprojectile.transform(toCms);
121 Ptarget.transform(toCms); 115 Ptarget.transform(toCms);
122 116
123 #ifdef debugQuarkExchange 117 #ifdef debugQuarkExchange
124 G4cout << "Pprojectile in CMS " << Pproje 118 G4cout << "Pprojectile in CMS " << Pprojectile << G4endl;
125 G4cout << "Ptarget in CMS " << Ptarge 119 G4cout << "Ptarget in CMS " << Ptarget << G4endl;
126 #endif 120 #endif
127 G4double maxPtSquare=sqr(Ptarget.pz()); 121 G4double maxPtSquare=sqr(Ptarget.pz());
128 122
129 G4double ProjectileMinDiffrMass = Pprojectil 123 G4double ProjectileMinDiffrMass = Pprojectile.mag()/GeV;
130 G4double TargetMinDiffrMass = Ptarget.ma 124 G4double TargetMinDiffrMass = Ptarget.mag()/GeV;
131 125
132 G4double AveragePt2(0.); 126 G4double AveragePt2(0.);
133 127
134 G4int PDGcode=projectile->GetDefinition() 128 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
135 G4int absPDGcode=std::abs(PDGcode); 129 G4int absPDGcode=std::abs(PDGcode);
136 130
137 G4bool ProjectileDiffraction = true; 131 G4bool ProjectileDiffraction = true;
138 132
139 // Also for heavy hadrons, assume 50% probab <<
140 if ( absPDGcode > 1000 ) 133 if ( absPDGcode > 1000 ) { ProjectileDiffraction = G4UniformRand() <= 0.5; }
141 if ( (absPDGcode == 211) || (absPDGcode == 1 134 if ( (absPDGcode == 211) || (absPDGcode == 111) ) { ProjectileDiffraction = G4UniformRand() <= 0.66; }
142 if ( (absPDGcode == 321) || (absPDGcode == 3 135 if ( (absPDGcode == 321) || (absPDGcode == 311) ||
143 ( PDGcode == 130) || ( PDGcode == 3 136 ( PDGcode == 130) || ( PDGcode == 310) ) { ProjectileDiffraction = G4UniformRand() <= 0.5; }
144 if ( absPDGcode > 400 && absPDGcode < 600 <<
145 <<
146 //G4cout<<"ProjectileDiffr "<<ProjectileDiff <<
147 137
148 if ( ProjectileDiffraction ) { 138 if ( ProjectileDiffraction ) {
149 if ( absPDGcode > 1000 ) 139 if ( absPDGcode > 1000 ) //------Projectile is baryon --------
150 { 140 {
151 if ( absPDGcode > 4000 && absPDGcode < 6 << 141 ProjectileMinDiffrMass = 1.16; // GeV
152 { << 142 AveragePt2 = 0.3; // GeV^2
153 ProjectileMinDiffrMass = projectile->G <<
154 AveragePt2 = 0.3; <<
155 } <<
156 else <<
157 { <<
158 ProjectileMinDiffrMass = 1.16; <<
159 AveragePt2 = 0.3; <<
160 } <<
161 } 143 }
162 else if( absPDGcode == 211 || absPDGcode = 144 else if( absPDGcode == 211 || absPDGcode == 111) //------Projectile is Pion -----------
163 { 145 {
164 ProjectileMinDiffrMass = 1.0; 146 ProjectileMinDiffrMass = 1.0; // GeV
165 AveragePt2 = 0.3; 147 AveragePt2 = 0.3; // GeV^2
166 } 148 }
167 else if( absPDGcode == 321 || absPDGcode = 149 else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon
168 { 150 {
169 ProjectileMinDiffrMass = 1.1; 151 ProjectileMinDiffrMass = 1.1; // GeV
170 AveragePt2 = 0.3; 152 AveragePt2 = 0.3; // GeV^2
171 } 153 }
172 else if( absPDGcode == 22) 154 else if( absPDGcode == 22) //------Projectile is Gamma -----------
173 { 155 {
174 ProjectileMinDiffrMass = 0.25; 156 ProjectileMinDiffrMass = 0.25; // GeV
175 AveragePt2 = 0.36; 157 AveragePt2 = 0.36; // GeV^2
176 } 158 }
177 else if( absPDGcode > 400 && absPDGcode < <<
178 { <<
179 ProjectileMinDiffrMass = projectile->Get <<
180 AveragePt2 = 0.3; <<
181 } <<
182 else 159 else //------Projectile is undefined, Nucleon assumed
183 { 160 {
184 ProjectileMinDiffrMass = 1.1; 161 ProjectileMinDiffrMass = 1.1; // GeV
185 AveragePt2 = 0.3; 162 AveragePt2 = 0.3; // GeV^2
186 }; 163 };
187 164
188 ProjectileMinDiffrMass = ProjectileMinDiff 165 ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV;
189 Mprojectile2=sqr(ProjectileMinDiffrMass); 166 Mprojectile2=sqr(ProjectileMinDiffrMass);
190 167
191 if (G4UniformRand() <= 0.5) TargetMinDiffr 168 if (G4UniformRand() <= 0.5) TargetMinDiffrMass += 0.22;
192 TargetMinDiffrMass *= GeV; 169 TargetMinDiffrMass *= GeV;
193 Mtarget2 = sqr( TargetMinDiffrMass) ; 170 Mtarget2 = sqr( TargetMinDiffrMass) ;
194 } 171 }
195 else 172 else
196 { 173 {
197 if (G4UniformRand() <= 0.5) ProjectileMinD 174 if (G4UniformRand() <= 0.5) ProjectileMinDiffrMass += 0.22;
198 ProjectileMinDiffrMass *=GeV; 175 ProjectileMinDiffrMass *=GeV;
199 Mprojectile2=sqr(ProjectileMinDiffrMass); 176 Mprojectile2=sqr(ProjectileMinDiffrMass);
200 177
201 TargetMinDiffrMass = 1.16*GeV; 178 TargetMinDiffrMass = 1.16*GeV; // For target nucleon
202 Mtarget2 = sqr( TargetMinDiffrMass) ; 179 Mtarget2 = sqr( TargetMinDiffrMass) ;
203 AveragePt2 = 0.3; 180 AveragePt2 = 0.3; // GeV^2
204 } // end of if ( ProjectileDiffraction ) 181 } // end of if ( ProjectileDiffraction )
205 182
206 AveragePt2 = AveragePt2 * GeV*GeV; 183 AveragePt2 = AveragePt2 * GeV*GeV;
207 184
208 if ( SqrtS - (ProjectileMinDiffrMass+TargetM << 185 if ( SqrtS - (ProjectileMinDiffrMass+TargetMinDiffrMass) < 220* MeV ) return false;
209 186
210 //----------------------- 187 //-----------------------
211 G4double Pt2, PZcms, PZcms2; 188 G4double Pt2, PZcms, PZcms2;
212 G4double ProjMassT2, ProjMassT; 189 G4double ProjMassT2, ProjMassT;
213 G4double TargMassT2, TargMassT; 190 G4double TargMassT2, TargMassT;
214 G4double PMinusMin, PMinusMax, sqrtPMinusMi 191 G4double PMinusMin, PMinusMax, sqrtPMinusMin, sqrtPMinusMax;
215 G4double TPlusMin, TPlusMax, sqrtTPlusMin 192 G4double TPlusMin, TPlusMax, sqrtTPlusMin, sqrtTPlusMax;
216 G4double PMinusNew, PPlusNew, TPlusNew(0.), 193 G4double PMinusNew, PPlusNew, TPlusNew(0.), TMinusNew;
217 194
218 G4LorentzVector Qmomentum; 195 G4LorentzVector Qmomentum;
219 G4double Qminus, Qplus; 196 G4double Qminus, Qplus;
220 197
221 G4double x(0.), y(0.); 198 G4double x(0.), y(0.);
222 G4int whilecount=0; 199 G4int whilecount=0;
223 do { 200 do {
224 whilecount++; 201 whilecount++;
225 202
226 if (whilecount > 1000 ) 203 if (whilecount > 1000 )
227 { 204 {
228 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 205 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
229 return false; // Ignore this intera 206 return false; // Ignore this interaction
230 } 207 }
231 208
232 // Generate pt 209 // Generate pt
233 Qmomentum=G4LorentzVector(GaussianPt(Avera 210 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
234 211
235 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma 212 Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
236 ProjMassT2 = Mprojectile2 + Pt2; 213 ProjMassT2 = Mprojectile2 + Pt2;
237 ProjMassT = std::sqrt( ProjMassT2 ); 214 ProjMassT = std::sqrt( ProjMassT2 );
238 TargMassT2 = Mtarget2 + Pt2; 215 TargMassT2 = Mtarget2 + Pt2;
239 TargMassT = std::sqrt( TargMassT2 ); 216 TargMassT = std::sqrt( TargMassT2 );
240 217
241 #ifdef debugQuarkExchange 218 #ifdef debugQuarkExchange
242 G4cout<<"whilecount Pt2 ProjMassT Tar 219 G4cout<<"whilecount Pt2 ProjMassT TargMassT SqrtS S ProjectileDiffraction"<<G4endl;
243 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjM 220 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl;
244 #endif 221 #endif
245 222
246 if ( SqrtS < ProjMassT + TargMassT + 220.0 223 if ( SqrtS < ProjMassT + TargMassT + 220.0*MeV ) continue;
247 224
248 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + T 225 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
249 - 2.0*S*ProjMassT2 - 2.0*S*Targ 226 - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
250 227
251 if ( PZcms2 < 0 ) continue; 228 if ( PZcms2 < 0 ) continue;
252 229
253 PZcms = std::sqrt( PZcms2 ); 230 PZcms = std::sqrt( PZcms2 );
254 231
255 if ( ProjectileDiffraction ) 232 if ( ProjectileDiffraction )
256 { // The projectile will fragment, the tar 233 { // The projectile will fragment, the target will saved.
257 PMinusMin = std::sqrt( ProjMassT2 + PZcm 234 PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
258 PMinusMax = SqrtS - TargMassT; 235 PMinusMax = SqrtS - TargMassT;
259 sqrtPMinusMin = std::sqrt(PMinusMin); sq 236 sqrtPMinusMin = std::sqrt(PMinusMin); sqrtPMinusMax = std::sqrt(PMinusMax);
260 237
261 if ( absPDGcode > 1000 ) 238 if ( absPDGcode > 1000 )
262 { 239 {
263 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinus 240 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinusMin/PMinusMax)
264 * G4Pow:: 241 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
265 } 242 }
266 else if ( (absPDGcode == 211) || (absPDG 243 else if ( (absPDGcode == 211) || (absPDGcode == 111) )
267 { 244 {
268 while (true) 245 while (true)
269 { 246 {
270 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtP 247 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusMin)*G4UniformRand();
271 y=G4UniformRand(); 248 y=G4UniformRand();
272 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) 249 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) break;
273 } 250 }
274 PMinusNew = sqr(x); 251 PMinusNew = sqr(x);
275 } 252 }
276 else if ( (absPDGcode == 321) || (absPDG 253 else if ( (absPDGcode == 321) || (absPDGcode == 311) ||
277 ( PDGcode == 130) || ( PDGcode = 254 ( PDGcode == 130) || ( PDGcode == 310) )
278 { // For K-mesons it must be found !!! 255 { // For K-mesons it must be found !!! Uzhi
279 while (true) 256 while (true)
280 { 257 {
281 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusM 258 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusMin)*G4UniformRand();
282 y=G4UniformRand(); 259 y=G4UniformRand();
283 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) 260 if ( y < 1.0-0.7 * x/sqrtPMinusMax ) break;
284 } 261 }
285 PMinusNew = sqr(x); 262 PMinusNew = sqr(x);
286 } 263 }
287 else 264 else
288 { 265 {
289 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinus 266 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinusMin/PMinusMax)
290 * G4Pow::GetI 267 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
291 }; 268 };
292 269
293 TMinusNew = SqrtS - PMinusNew; 270 TMinusNew = SqrtS - PMinusNew;
294 271
295 Qminus = Ptarget.minus() - TMinusNew; 272 Qminus = Ptarget.minus() - TMinusNew;
296 TPlusNew = TargMassT2 / TMinusNew; 273 TPlusNew = TargMassT2 / TMinusNew;
297 Qplus = Ptarget.plus() - TPlusNew; 274 Qplus = Ptarget.plus() - TPlusNew;
298 275
299 } 276 }
300 else 277 else
301 { // The target will fragment, the project 278 { // The target will fragment, the projectile will saved.
302 TPlusMin = std::sqrt( TargMassT2 + PZcms 279 TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
303 TPlusMax = SqrtS - ProjMassT; 280 TPlusMax = SqrtS - ProjMassT;
304 sqrtTPlusMin = std::sqrt(TPlusMin); sqrt 281 sqrtTPlusMin = std::sqrt(TPlusMin); sqrtTPlusMax = std::sqrt(TPlusMax);
305 282
306 if ( absPDGcode > 1000 ) 283 if ( absPDGcode > 1000 )
307 { 284 {
308 TPlusNew = TPlusMax * (1.0 - (1.0 - TP 285 TPlusNew = TPlusMax * (1.0 - (1.0 - TPlusMin/TPlusMax)
309 * G4Pow::G 286 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
310 } 287 }
311 else if ( (absPDGcode == 211) || (absPDG 288 else if ( (absPDGcode == 211) || (absPDGcode == 111) )
312 { 289 {
313 while (true) 290 while (true)
314 { 291 {
315 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPl 292 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPlusMin)*G4UniformRand();
316 y=G4UniformRand(); 293 y=G4UniformRand();
317 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) 294 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break;
318 } 295 }
319 TPlusNew = sqr(x); 296 TPlusNew = sqr(x);
320 } 297 }
321 else if ( (absPDGcode == 321) || (absPDG 298 else if ( (absPDGcode == 321) || (absPDGcode == 311) ||
322 ( PDGcode == 130) || ( PDGcode = 299 ( PDGcode == 130) || ( PDGcode == 310) )
323 { // For K-mesons it must be found !!! U 300 { // For K-mesons it must be found !!! Uzhi
324 while (true) 301 while (true)
325 { 302 {
326 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPl 303 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPlusMin)*G4UniformRand();
327 y=G4UniformRand(); 304 y=G4UniformRand();
328 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break; 305 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break;
329 } 306 }
330 } 307 }
331 else 308 else
332 { 309 {
333 TPlusNew = TPlusMax * (1.0 - (1.0 - TP 310 TPlusNew = TPlusMax * (1.0 - (1.0 - TPlusMin/TPlusMax)
334 * G4Pow::G 311 * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
335 }; 312 };
336 313
337 PPlusNew = SqrtS - TPlusNew; 314 PPlusNew = SqrtS - TPlusNew;
338 315
339 Qplus = PPlusNew - Pprojectile.plus(); 316 Qplus = PPlusNew - Pprojectile.plus();
340 PMinusNew = ProjMassT2 / PPlusNew; 317 PMinusNew = ProjMassT2 / PPlusNew;
341 Qminus = PMinusNew - Pprojectile.minus() 318 Qminus = PMinusNew - Pprojectile.minus();
342 } 319 }
343 320
344 Qmomentum.setPz( (Qplus - Qminus)/2 ); 321 Qmomentum.setPz( (Qplus - Qminus)/2 );
345 Qmomentum.setE( (Qplus + Qminus)/2 ); 322 Qmomentum.setE( (Qplus + Qminus)/2 );
346 323
347 #ifdef debugQuarkExchange 324 #ifdef debugQuarkExchange
348 G4cout<<"ProjectileDiffraction (Pproject 325 G4cout<<"ProjectileDiffraction (Pprojectile + Qmomentum).mag2() Mprojectile2"<<G4endl;
349 G4cout<<ProjectileDiffraction<<" "<<( Pp 326 G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl;
350 G4cout<<"TargetDiffraction (Ptarget 327 G4cout<<"TargetDiffraction (Ptarget - Qmomentum).mag2() Mtarget2"<<G4endl;
351 G4cout<<!ProjectileDiffraction<<" "<<( P 328 G4cout<<!ProjectileDiffraction<<" "<<( Ptarget - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl;
352 #endif 329 #endif
353 330
354 } while ( ( ProjectileDiffraction&&( Pprojec 331 } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) ||
355 (!ProjectileDiffraction&&( Ptarget 332 (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) );
356 // Repeat the sampling because there was n 333 // Repeat the sampling because there was not any excitation
357 334
358 Pprojectile += Qmomentum; 335 Pprojectile += Qmomentum;
359 336
360 Ptarget -= Qmomentum; 337 Ptarget -= Qmomentum;
361 338
362 // Transform back and update SplitableHadron 339 // Transform back and update SplitableHadron Participant.
363 Pprojectile.transform(toLab); 340 Pprojectile.transform(toLab);
364 Ptarget.transform(toLab); 341 Ptarget.transform(toLab);
365 342
366 #ifdef debugQuarkExchange 343 #ifdef debugQuarkExchange
367 G4cout << "Pprojectile in Lab. " << Pproj 344 G4cout << "Pprojectile in Lab. " << Pprojectile << G4endl;
368 G4cout << "Ptarget in Lab. " << Ptarg 345 G4cout << "Ptarget in Lab. " << Ptarget << G4endl;
369 G4cout << "G4QuarkExchange: Projectile mas 346 G4cout << "G4QuarkExchange: Projectile mass " << Pprojectile.mag() << G4endl;
370 G4cout << "G4QuarkExchange: Target mass 347 G4cout << "G4QuarkExchange: Target mass " << Ptarget.mag() << G4endl;
371 #endif 348 #endif
372 349
373 target->Set4Momentum(Ptarget); 350 target->Set4Momentum(Ptarget);
374 projectile->Set4Momentum(Pprojectile); 351 projectile->Set4Momentum(Pprojectile);
375 352
376 //=================================== Quark 353 //=================================== Quark exchange ================================
377 projectile->SplitUp(); 354 projectile->SplitUp();
378 target->SplitUp(); 355 target->SplitUp();
379 356
380 G4Parton* PrQuark = nullptr; << 357 G4Parton* PrQuark = projectile->GetNextParton();
381 G4Parton* TrQuark = nullptr; << 358 G4Parton* TrQuark = target->GetNextParton();
382 359
383 if( projectile->GetDefinition()->GetBaryonNu << 360 G4ParticleDefinition * Tmp = PrQuark->GetDefinition();
384 // Quark exchange ---- << 361 PrQuark->SetDefinition(TrQuark->GetDefinition());
385 PrQuark = projectile->GetNextParton(); << 362 TrQuark->SetDefinition(Tmp);
386 TrQuark = target->GetNextParton(); <<
387 G4ParticleDefinition * Tmp = PrQuark->GetD <<
388 PrQuark->SetDefinition(TrQuark->GetDefinit <<
389 TrQuark->SetDefinition(Tmp); <<
390 return true; <<
391 } <<
392 <<
393 // Quark exchage for projectile anti-baryon <<
394 // This part added by Uzhi June 2020 <<
395 PrQuark = projectile->GetNextAntiParton(); <<
396 TrQuark = target->GetNextParton(); <<
397 if( -PrQuark->GetDefinition()->GetPDGEncodin <<
398 G4int QuarkCode = 1 + (int)(G4UniformRand( <<
399 G4ParticleDefinition* AQpr = G4ParticleTab <<
400 G4ParticleDefinition* Qtr = G4ParticleTab <<
401 if( (AQpr != nullptr) && (Qtr != nullptr) <<
402 PrQuark->SetDefinition(AQpr); <<
403 TrQuark->SetDefinition( Qtr); <<
404 } <<
405 } <<
406 363
407 return true; 364 return true;
408 } 365 }
409 366
410 367
411 // --------- private methods ----------------- 368 // --------- private methods ----------------------
412 369
413 G4ThreeVector G4QuarkExchange::GaussianPt(G4do 370 G4ThreeVector G4QuarkExchange::GaussianPt(G4double widthSquare, G4double maxPtSquare) const
414 { // @@ this method is used in FTF 371 { // @@ this method is used in FTFModel as well. Should go somewhere common!
415 372
416 G4double pt2; 373 G4double pt2;
417 374
418 const G4int maxNumberOfLoops = 1000; 375 const G4int maxNumberOfLoops = 1000;
419 G4int loopCounter = 0; 376 G4int loopCounter = 0;
420 do { 377 do {
421 pt2=-widthSquare * G4Log( G4UniformRand() 378 pt2=-widthSquare * G4Log( G4UniformRand() );
422 } while ( ( pt2 > maxPtSquare) && ++loopCoun 379 } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
423 if ( loopCounter >= maxNumberOfLoops ) { 380 if ( loopCounter >= maxNumberOfLoops ) {
424 pt2 = 0.99*maxPtSquare; // Just an accept 381 pt2 = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration.
425 } 382 }
426 383
427 pt2=std::sqrt(pt2); 384 pt2=std::sqrt(pt2);
428 385
429 G4double phi=G4UniformRand() * twopi; 386 G4double phi=G4UniformRand() * twopi;
430 387
431 return G4ThreeVector (pt2*std::cos(phi), pt2 388 return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);
432 } 389 }
433 390
434 391