| 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 "G4QGSMSplitableHadron.hh" 26 #include "G4QGSMSplitableHadron.hh"
27 #include "G4PhysicalConstants.hh" 27 #include "G4PhysicalConstants.hh"
28 #include "G4SystemOfUnits.hh" 28 #include "G4SystemOfUnits.hh"
29 #include "G4ParticleTable.hh" 29 #include "G4ParticleTable.hh"
30 #include "G4PionPlus.hh" 30 #include "G4PionPlus.hh"
31 #include "G4PionMinus.hh" 31 #include "G4PionMinus.hh"
32 #include "G4Gamma.hh" 32 #include "G4Gamma.hh"
33 #include "G4PionZero.hh" 33 #include "G4PionZero.hh"
34 #include "G4KaonPlus.hh" 34 #include "G4KaonPlus.hh"
35 #include "G4KaonMinus.hh" 35 #include "G4KaonMinus.hh"
36 36
37 #include "G4Log.hh" 37 #include "G4Log.hh"
38 #include "G4Pow.hh" 38 #include "G4Pow.hh"
39 39
40 // based on prototype by Maxim Komogorov 40 // based on prototype by Maxim Komogorov
41 // Splitting into methods, and centralizing of 41 // Splitting into methods, and centralizing of model parameters HPW Feb 1999
42 // restructuring HPW Feb 1999 42 // restructuring HPW Feb 1999
43 // fixing bug in the sampling of 'x', HPW Feb 43 // fixing bug in the sampling of 'x', HPW Feb 1999
44 // fixing bug in sampling pz, HPW Feb 1999. 44 // fixing bug in sampling pz, HPW Feb 1999.
45 // Code now also good for p-nucleus scattering 45 // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999.
46 // Using Parton more directly, HPW Feb 1999. 46 // Using Parton more directly, HPW Feb 1999.
47 // Shortening the algorithm for sampling x, HP 47 // Shortening the algorithm for sampling x, HPW Feb 1999.
48 // sampling of x replaced by formula, taking X 48 // sampling of x replaced by formula, taking X_min into account in the correlated sampling. HPW, Feb 1999.
49 // logic much clearer now. HPW Feb 1999 49 // logic much clearer now. HPW Feb 1999
50 // Removed the ordering problem. No Direction 50 // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99.
51 // Fixing p-t distributions for scattering of 51 // Fixing p-t distributions for scattering of nuclei.
52 // Separating out parameters. 52 // Separating out parameters.
53 53
54 void G4QGSMSplitableHadron::InitParameters() 54 void G4QGSMSplitableHadron::InitParameters()
55 { 55 {
56 // changing rapidity distribution for all 56 // changing rapidity distribution for all
57 alpha = -0.5; // Note that this number is st 57 alpha = -0.5; // Note that this number is still assumed in the algorithm
58 // needs to be generalized. 58 // needs to be generalized.
59 // changing rapidity distribution for projec 59 // changing rapidity distribution for projectile like
60 beta = 2.5;// Note that this number is still 60 beta = 2.5;// Note that this number is still assumed in the algorithm
61 // needs to be generalized. 61 // needs to be generalized.
62 theMinPz = 0.5*G4PionMinus::PionMinus()->Get 62 theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass();
63 // theMinPz = 0.1*G4PionMinus::PionMinus()- 63 // theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass();
64 // theMinPz = G4PionMinus::PionMinus()->Get 64 // theMinPz = G4PionMinus::PionMinus()->GetPDGMass();
65 // as low as possible, otherwise, we have un 65 // as low as possible, otherwise, we have unphysical boundary conditions in the sampling.
66 StrangeSuppress = 0.48; 66 StrangeSuppress = 0.48;
67 sigmaPt = 0.*GeV; // widens eta slightly, if 67 sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7,
68 // but Maxim's algorithm breaks energy conse 68 // but Maxim's algorithm breaks energy conservation to be revised.
69 widthOfPtSquare = 0.5*sqr(GeV); 69 widthOfPtSquare = 0.5*sqr(GeV);
70 Direction = FALSE; 70 Direction = FALSE;
71 minTransverseMass = 1*keV; 71 minTransverseMass = 1*keV;
72 iP =0;// Color.begin(); 72 iP =0;// Color.begin();
73 iAP =0;// AntiColor.begin(); 73 iAP =0;// AntiColor.begin();
74 } 74 }
75 75
76 G4QGSMSplitableHadron::G4QGSMSplitableHadron() 76 G4QGSMSplitableHadron::G4QGSMSplitableHadron()
77 { 77 {
78 InitParameters(); 78 InitParameters();
79 } 79 }
80 80
81 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 81 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection)
82 :G4VSplitableHadron(aPrimary) 82 :G4VSplitableHadron(aPrimary)
83 { 83 {
84 InitParameters(); 84 InitParameters();
85 Direction = aDirection; 85 Direction = aDirection;
86 } 86 }
87 87
88 88
89 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 89 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary)
90 : G4VSplitableHadron(aPrimary) 90 : G4VSplitableHadron(aPrimary)
91 { 91 {
92 InitParameters(); 92 InitParameters();
93 } 93 }
94 94
95 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 95 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon)
96 : G4VSplitableHadron(aNucleon) 96 : G4VSplitableHadron(aNucleon)
97 { 97 {
98 InitParameters(); 98 InitParameters();
99 } 99 }
100 100
101 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 101 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection)
102 : G4VSplitableHadron(aNucleon) 102 : G4VSplitableHadron(aNucleon)
103 { 103 {
104 InitParameters(); 104 InitParameters();
105 Direction = aDirection; 105 Direction = aDirection;
106 } 106 }
107 107
108 G4QGSMSplitableHadron::~G4QGSMSplitableHadron( 108 G4QGSMSplitableHadron::~G4QGSMSplitableHadron() {}
109 109
110 110
111 //******************************************** 111 //**************************************************************************************************************************
112 112
113 void G4QGSMSplitableHadron::SplitUp() 113 void G4QGSMSplitableHadron::SplitUp()
114 { 114 {
115 if (IsSplit()) return; 115 if (IsSplit()) return;
116 Splitting(); // To mark that a hadron is sp 116 Splitting(); // To mark that a hadron is split
117 if (Color.size()!=0) return; 117 if (Color.size()!=0) return;
118 if (GetSoftCollisionCount() == 0) // GetSof 118 if (GetSoftCollisionCount() == 0) // GetSoftCollisionCount() from G4VSplitableHadron
119 { 119 {
120 DiffractiveSplitUp(); 120 DiffractiveSplitUp();
121 } 121 }
122 else 122 else
123 { 123 {
124 SoftSplitUp(); 124 SoftSplitUp();
125 } 125 }
126 } 126 }
127 127
128 void G4QGSMSplitableHadron::DiffractiveSplitUp 128 void G4QGSMSplitableHadron::DiffractiveSplitUp()
129 { 129 {
130 // take the particle definitions and get the 130 // take the particle definitions and get the partons HPW
131 G4Parton * Left = NULL; 131 G4Parton * Left = NULL;
132 G4Parton * Right = NULL; 132 G4Parton * Right = NULL;
133 GetValenceQuarkFlavors(GetDefinition(), Left 133 GetValenceQuarkFlavors(GetDefinition(), Left, Right);
134 Left->SetPosition(GetPosition()); 134 Left->SetPosition(GetPosition());
135 Right->SetPosition(GetPosition()); 135 Right->SetPosition(GetPosition());
136 136
137 G4LorentzVector HadronMom = Get4Momentum(); 137 G4LorentzVector HadronMom = Get4Momentum();
138 138
139 G4double maxAvailMomentum2 = sqr(HadronMom.m 139 G4double maxAvailMomentum2 = sqr(HadronMom.mag()/2.);
140 140
141 G4ThreeVector pt(minTransverseMass, minTrans 141 G4ThreeVector pt(minTransverseMass, minTransverseMass, 0);
142 if (maxAvailMomentum2/widthOfPtSquare>0.01) 142 if (maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2);
143 143
144 G4LorentzVector LeftMom(pt, 0.); 144 G4LorentzVector LeftMom(pt, 0.);
145 G4LorentzVector RightMom; 145 G4LorentzVector RightMom;
146 RightMom.setPx(HadronMom.px() - pt.x()); 146 RightMom.setPx(HadronMom.px() - pt.x());
147 RightMom.setPy(HadronMom.py() - pt.y()); 147 RightMom.setPy(HadronMom.py() - pt.y());
148 148
149 G4double Local1 = HadronMom.minus() + (Right 149 G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus();
150 G4double Local2 = std::sqrt(std::max(0., sqr 150 G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus()));
151 151
152 if (Direction) Local2 = -Local2; 152 if (Direction) Local2 = -Local2;
153 G4double RightMinus = 0.5*(Local1 + Local2 153 G4double RightMinus = 0.5*(Local1 + Local2);
154 G4double LeftMinus = HadronMom.minus() - Rig 154 G4double LeftMinus = HadronMom.minus() - RightMinus;
155 155
156 if (LeftMinus <= 0.) { 156 if (LeftMinus <= 0.) {
157 RightMinus = 0.5*(Local1 - Local2); 157 RightMinus = 0.5*(Local1 - Local2);
158 LeftMinus = HadronMom.minus() - RightMinus 158 LeftMinus = HadronMom.minus() - RightMinus;
159 } 159 }
160 160
161 G4double LeftPlus = LeftMom.perp2()/LeftMin 161 G4double LeftPlus = LeftMom.perp2()/LeftMinus;
162 G4double RightPlus = HadronMom.plus() - Left 162 G4double RightPlus = HadronMom.plus() - LeftPlus;
163 163
164 LeftMom.setPz(0.5*(LeftPlus - LeftMinus)); 164 LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
165 LeftMom.setE (0.5*(LeftPlus + LeftMinus)); 165 LeftMom.setE (0.5*(LeftPlus + LeftMinus));
166 RightMom.setPz(0.5*(RightPlus - RightMinus)) 166 RightMom.setPz(0.5*(RightPlus - RightMinus));
167 RightMom.setE (0.5*(RightPlus + RightMinus)) 167 RightMom.setE (0.5*(RightPlus + RightMinus));
168 168
169 Left->Set4Momentum(LeftMom); 169 Left->Set4Momentum(LeftMom);
170 Right->Set4Momentum(RightMom); 170 Right->Set4Momentum(RightMom);
171 171
172 Color.push_back(Left); 172 Color.push_back(Left);
173 AntiColor.push_back(Right); 173 AntiColor.push_back(Right);
174 iP=0; iAP=0; 174 iP=0; iAP=0;
175 } 175 }
176 176
177 177
178 void G4QGSMSplitableHadron::SoftSplitUp() 178 void G4QGSMSplitableHadron::SoftSplitUp()
179 { 179 {
180 G4int nSeaPair = GetSoftCollisionCount()-1; 180 G4int nSeaPair = GetSoftCollisionCount()-1;
181 181
182 G4LorentzVector tmp(0., 0., 0., 0.); 182 G4LorentzVector tmp(0., 0., 0., 0.);
183 183
184 G4int aSeaPair; 184 G4int aSeaPair;
185 for (aSeaPair = 0; aSeaPair < nSeaPair; aSea 185 for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
186 { 186 {
187 // choose quark flavour, d:u:s = 1:1:(1/S 187 // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
188 G4int aPDGCode = 1 + (G4int)(G4UniformRand 188 G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
189 189
190 // BuildSeaQuark() determines quark spin, 190 // BuildSeaQuark() determines quark spin, isospin and colour
191 // via parton-constructor G4Parton(aPDGCo 191 // via parton-constructor G4Parton(aPDGCode)
192 G4Parton * aParton = BuildSeaQuark(false, 192 G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
193 193
194 G4int firstPartonColour = aParton->GetColo 194 G4int firstPartonColour = aParton->GetColour();
195 G4double firstPartonSpinZ = aParton->GetSp 195 G4double firstPartonSpinZ = aParton->GetSpinZ();
196 196
197 aParton->Set4Momentum(tmp); 197 aParton->Set4Momentum(tmp);
198 Color.push_back(aParton); 198 Color.push_back(aParton);
199 199
200 // create anti-quark 200 // create anti-quark
201 aParton = BuildSeaQuark(true, aPDGCode, nS 201 aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
202 aParton->SetSpinZ(-firstPartonSpinZ); 202 aParton->SetSpinZ(-firstPartonSpinZ);
203 aParton->SetColour(-firstPartonColour); 203 aParton->SetColour(-firstPartonColour);
204 AntiColor.push_back(aParton); 204 AntiColor.push_back(aParton);
205 } 205 }
206 206
207 // Valence quark 207 // Valence quark
208 G4Parton* pColorParton = NULL; 208 G4Parton* pColorParton = NULL;
209 G4Parton* pAntiColorParton = NULL; 209 G4Parton* pAntiColorParton = NULL;
210 GetValenceQuarkFlavors(GetDefinition(), pCol 210 GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton);
211 211
212 pColorParton->Set4Momentum(tmp); 212 pColorParton->Set4Momentum(tmp);
213 pAntiColorParton->Set4Momentum(tmp); 213 pAntiColorParton->Set4Momentum(tmp);
214 214
215 Color.push_back(pColorParton); 215 Color.push_back(pColorParton);
216 AntiColor.push_back(pAntiColorParton); 216 AntiColor.push_back(pAntiColorParton);
217 217
218 iP=0; iAP=0; 218 iP=0; iAP=0;
219 219
220 return; 220 return;
221 } 221 }
222 222
223 223
224 void G4QGSMSplitableHadron::GetValenceQuarkFla 224 void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart,
225 225 G4Parton *& Parton1, G4Parton *& Parton2)
226 { 226 {
227 // Note! convention aEnd = q or (qq)bar and 227 // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
228 G4int aEnd=0; << 228 G4int aEnd;
229 G4int bEnd=0; << 229 G4int bEnd;
230 G4int HadronEncoding = aPart->GetPDGEncoding 230 G4int HadronEncoding = aPart->GetPDGEncoding();
231 if (aPart->GetBaryonNumber() == 0) 231 if (aPart->GetBaryonNumber() == 0)
232 { 232 {
233 theMesonSplitter.SplitMeson(HadronEncoding << 233 theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd);
234 } 234 }
235 else 235 else
236 { 236 {
237 theBaryonSplitter.SplitBarion(HadronEncodi << 237 theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd);
238 } 238 }
239 239
240 Parton1 = new G4Parton(aEnd); 240 Parton1 = new G4Parton(aEnd);
241 Parton1->SetPosition(GetPosition()); 241 Parton1->SetPosition(GetPosition());
242 242
243 Parton2 = new G4Parton(bEnd); 243 Parton2 = new G4Parton(bEnd);
244 Parton2->SetPosition(GetPosition()); 244 Parton2->SetPosition(GetPosition());
245 245
246 // colour of parton 1 choosen at random by G 246 // colour of parton 1 choosen at random by G4Parton(aEnd)
247 // colour of parton 2 is the opposite: 247 // colour of parton 2 is the opposite:
248 248
249 Parton2->SetColour(-(Parton1->GetColour())); 249 Parton2->SetColour(-(Parton1->GetColour()));
250 250
251 // isospin-3 of both partons is handled by G 251 // isospin-3 of both partons is handled by G4Parton(PDGCode)
252 252
253 // spin-3 of parton 1 and 2 choosen at rando 253 // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd)
254 // spin-3 of parton 2 may be constrained by 254 // spin-3 of parton 2 may be constrained by spin of original particle:
255 255
256 if ( std::abs(Parton1->GetSpinZ() + Parton2- 256 if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin())
257 { 257 {
258 Parton2->SetSpinZ(-(Parton2->GetSpinZ())); 258 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
259 } 259 }
260 } 260 }
261 261
262 262
263 G4ThreeVector G4QGSMSplitableHadron::GaussianP 263 G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
264 { 264 {
265 G4double R; 265 G4double R;
266 const G4int maxNumberOfLoops = 1000; 266 const G4int maxNumberOfLoops = 1000;
267 G4int loopCounter = -1; 267 G4int loopCounter = -1;
268 while( ((R = -widthSquare*G4Log(G4UniformRan 268 while( ((R = -widthSquare*G4Log(G4UniformRand())) > maxPtSquare) &&
269 ++loopCounter < maxNumberOfLoops ) {; 269 ++loopCounter < maxNumberOfLoops ) {;} /* Loop checking, 07.08.2015, A.Ribon */
270 if ( loopCounter >= maxNumberOfLoops ) { 270 if ( loopCounter >= maxNumberOfLoops ) {
271 R = 0.99*maxPtSquare; // Just an acceptab 271 R = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration.
272 } 272 }
273 R = std::sqrt(R); 273 R = std::sqrt(R);
274 G4double phi = twopi*G4UniformRand(); 274 G4double phi = twopi*G4UniformRand();
275 return G4ThreeVector (R*std::cos(phi), R*std 275 return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.);
276 } 276 }
277 277
278 278
279 G4Parton * G4QGSMSplitableHadron:: 279 G4Parton * G4QGSMSplitableHadron::
280 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCo 280 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/)
281 { 281 {
282 if (isAntiQuark) aPDGCode*=-1; 282 if (isAntiQuark) aPDGCode*=-1;
283 G4Parton* result = new G4Parton(aPDGCode); 283 G4Parton* result = new G4Parton(aPDGCode);
284 result->SetPosition(GetPosition()); 284 result->SetPosition(GetPosition());
285 G4ThreeVector aPtVector = GaussianPt(sigmaPt 285 G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
286 G4LorentzVector a4Momentum(aPtVector, 0); 286 G4LorentzVector a4Momentum(aPtVector, 0);
287 result->Set4Momentum(a4Momentum); 287 result->Set4Momentum(a4Momentum);
288 return result; 288 return result;
289 } 289 }
290 290
291 291
292 G4double G4QGSMSplitableHadron:: 292 G4double G4QGSMSplitableHadron::
293 SampleX(G4double anXmin, G4int nSea, G4int tot 293 SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
294 { 294 {
295 G4double result; 295 G4double result;
296 G4double x1, x2; 296 G4double x1, x2;
297 G4double ymax = 0; 297 G4double ymax = 0;
298 for(G4int ii=1; ii<100; ii++) 298 for(G4int ii=1; ii<100; ii++)
299 { 299 {
300 G4double y = G4Pow::GetInstance()->powA(1. 300 G4double y = G4Pow::GetInstance()->powA(1./G4double(ii), alpha);
301 y *= G4Pow::GetInstance()->powN( G4Pow::Ge 301 y *= G4Pow::GetInstance()->powN( G4Pow::GetInstance()->powA(1-anXmin-totalSea*anXmin, alpha+1) -
302 G4Pow::Ge 302 G4Pow::GetInstance()->powA(anXmin, alpha+1), nSea );
303 y *= G4Pow::GetInstance()->powA(1-anXmin-t 303 y *= G4Pow::GetInstance()->powA(1-anXmin-totalSea*anXmin, aBeta+1) -
304 G4Pow::GetInstance()->powA(anXmin, aB 304 G4Pow::GetInstance()->powA(anXmin, aBeta+1);
305 if (y>ymax) ymax = y; 305 if (y>ymax) ymax = y;
306 } 306 }
307 G4double y; 307 G4double y;
308 G4double xMax=1-(totalSea+1)*anXmin; 308 G4double xMax=1-(totalSea+1)*anXmin;
309 if (anXmin > xMax) 309 if (anXmin > xMax)
310 { 310 {
311 throw G4HadronicException(__FILE__, __LINE 311 throw G4HadronicException(__FILE__, __LINE__,
312 "G4QGSMSplitableHadron - Fatal: Ca 312 "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints.");
313 } 313 }
314 const G4int maxNumberOfLoops = 1000; 314 const G4int maxNumberOfLoops = 1000;
315 G4int loopCounter = 0; 315 G4int loopCounter = 0;
316 do 316 do
317 { 317 {
318 x1 = G4RandFlat::shoot(anXmin, xMax); 318 x1 = G4RandFlat::shoot(anXmin, xMax);
319 y = G4Pow::GetInstance()->powA(x1, alpha); 319 y = G4Pow::GetInstance()->powA(x1, alpha);
320 y *= G4Pow::GetInstance()->powN( G4Pow::Ge 320 y *= G4Pow::GetInstance()->powN( G4Pow::GetInstance()->powA(1-x1-totalSea*anXmin, alpha+1) -
321 G4Pow::Ge 321 G4Pow::GetInstance()->powA(anXmin, alpha+1), nSea );
322 y *= G4Pow::GetInstance()->powA(1-x1-total 322 y *= G4Pow::GetInstance()->powA(1-x1-totalSea*anXmin, aBeta+1) -
323 G4Pow::GetInstance()->powA(anXmin, aB 323 G4Pow::GetInstance()->powA(anXmin, aBeta+1);
324 x2 = ymax*G4UniformRand(); 324 x2 = ymax*G4UniformRand();
325 } while( (x2>y) && ++loopCounter < maxNumber 325 } while( (x2>y) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
326 if ( loopCounter >= maxNumberOfLoops ) { 326 if ( loopCounter >= maxNumberOfLoops ) {
327 x1 = 0.5*( anXmin + xMax ); // Just an ac 327 x1 = 0.5*( anXmin + xMax ); // Just an acceptable value, without any physics consideration.
328 } 328 }
329 result = x1; 329 result = x1;
330 return result; 330 return result;
331 } 331 }
332 332