1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | #include "G4QGSMSplitableHadron.hh" |
---|
27 | #include "G4ParticleTable.hh" |
---|
28 | #include "G4PionPlus.hh" |
---|
29 | #include "G4PionMinus.hh" |
---|
30 | #include "G4Gamma.hh" |
---|
31 | #include "G4PionZero.hh" |
---|
32 | #include "G4KaonPlus.hh" |
---|
33 | #include "G4KaonMinus.hh" |
---|
34 | |
---|
35 | // based on prototype by Maxim Komogorov |
---|
36 | // Splitting into methods, and centralizing of model parameters HPW Feb 1999 |
---|
37 | // restructuring HPW Feb 1999 |
---|
38 | // fixing bug in the sampling of 'x', HPW Feb 1999 |
---|
39 | // fixing bug in sampling pz, HPW Feb 1999. |
---|
40 | // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999. |
---|
41 | // Using Parton more directly, HPW Feb 1999. |
---|
42 | // Shortening the algorithm for sampling x, HPW Feb 1999. |
---|
43 | // sampling of x replaced by formula, taking X_min into account in the correlated sampling. HPW, Feb 1999. |
---|
44 | // logic much clearer now. HPW Feb 1999 |
---|
45 | // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99. |
---|
46 | // Fixing p-t distributions for scattering of nuclei. |
---|
47 | // Separating out parameters. |
---|
48 | |
---|
49 | void G4QGSMSplitableHadron::InitParameters() |
---|
50 | { |
---|
51 | // changing rapidity distribution for all |
---|
52 | alpha = -0.5; // Note that this number is still assumed in the algorithm |
---|
53 | // needs to be generalized. |
---|
54 | // changing rapidity distribution for projectile like |
---|
55 | beta = 2.5;// Note that this number is still assumed in the algorithm |
---|
56 | // needs to be generalized. |
---|
57 | theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass(); |
---|
58 | // theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass(); |
---|
59 | // theMinPz = G4PionMinus::PionMinus()->GetPDGMass(); |
---|
60 | // as low as possible, otherwise, we have unphysical boundary conditions in the sampling. |
---|
61 | StrangeSuppress = 0.48; |
---|
62 | sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7, |
---|
63 | // but Maxim's algorithm breaks energy conservation |
---|
64 | // to be revised. |
---|
65 | widthOfPtSquare = 0.01*GeV*GeV; |
---|
66 | Direction = FALSE; |
---|
67 | minTransverseMass = 1*keV; |
---|
68 | } |
---|
69 | |
---|
70 | G4QGSMSplitableHadron::G4QGSMSplitableHadron() |
---|
71 | { |
---|
72 | InitParameters(); |
---|
73 | } |
---|
74 | |
---|
75 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection) |
---|
76 | :G4VSplitableHadron(aPrimary) |
---|
77 | { |
---|
78 | InitParameters(); |
---|
79 | Direction = aDirection; |
---|
80 | } |
---|
81 | |
---|
82 | |
---|
83 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary) |
---|
84 | : G4VSplitableHadron(aPrimary) |
---|
85 | { |
---|
86 | InitParameters(); |
---|
87 | } |
---|
88 | |
---|
89 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon) |
---|
90 | : G4VSplitableHadron(aNucleon) |
---|
91 | { |
---|
92 | InitParameters(); |
---|
93 | } |
---|
94 | |
---|
95 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection) |
---|
96 | : G4VSplitableHadron(aNucleon) |
---|
97 | { |
---|
98 | InitParameters(); |
---|
99 | Direction = aDirection; |
---|
100 | } |
---|
101 | |
---|
102 | G4QGSMSplitableHadron::~G4QGSMSplitableHadron(){} |
---|
103 | |
---|
104 | const G4QGSMSplitableHadron & G4QGSMSplitableHadron::operator=(const G4QGSMSplitableHadron &) |
---|
105 | { |
---|
106 | throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron::operator= meant to not be accessable"); |
---|
107 | return *this; |
---|
108 | } |
---|
109 | |
---|
110 | |
---|
111 | //************************************************************************************************************************** |
---|
112 | |
---|
113 | void G4QGSMSplitableHadron::SplitUp() |
---|
114 | { |
---|
115 | if (IsSplit()) return; |
---|
116 | Splitting(); |
---|
117 | if (Color.size()!=0) return; |
---|
118 | if (GetSoftCollisionCount() == 0) |
---|
119 | { |
---|
120 | DiffractiveSplitUp(); |
---|
121 | } |
---|
122 | else |
---|
123 | { |
---|
124 | SoftSplitUp(); |
---|
125 | } |
---|
126 | } |
---|
127 | |
---|
128 | void G4QGSMSplitableHadron::DiffractiveSplitUp() |
---|
129 | { |
---|
130 | // take the particle definitions and get the partons HPW |
---|
131 | G4Parton * Left = NULL; |
---|
132 | G4Parton * Right = NULL; |
---|
133 | GetValenceQuarkFlavors(GetDefinition(), Left, Right); |
---|
134 | Left->SetPosition(GetPosition()); |
---|
135 | Right->SetPosition(GetPosition()); |
---|
136 | |
---|
137 | G4LorentzVector HadronMom = Get4Momentum(); |
---|
138 | //std::cout << "DSU 1 - "<<HadronMom<<std::endl; |
---|
139 | |
---|
140 | // momenta of string ends |
---|
141 | G4double pt2 = HadronMom.perp2(); |
---|
142 | G4double transverseMass2 = HadronMom.plus()*HadronMom.minus(); |
---|
143 | G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2)); |
---|
144 | G4ThreeVector pt(minTransverseMass, minTransverseMass, 0); |
---|
145 | if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2); |
---|
146 | //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl; |
---|
147 | |
---|
148 | G4LorentzVector LeftMom(pt, 0.); |
---|
149 | G4LorentzVector RightMom; |
---|
150 | RightMom.setPx(HadronMom.px() - pt.x()); |
---|
151 | RightMom.setPy(HadronMom.py() - pt.y()); |
---|
152 | //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl; |
---|
153 | |
---|
154 | G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus(); |
---|
155 | G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus())); |
---|
156 | //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl; |
---|
157 | if (Direction) Local2 = -Local2; |
---|
158 | G4double RightMinus = 0.5*(Local1 + Local2); |
---|
159 | G4double LeftMinus = HadronMom.minus() - RightMinus; |
---|
160 | //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl; |
---|
161 | |
---|
162 | G4double LeftPlus = LeftMom.perp2()/LeftMinus; |
---|
163 | G4double RightPlus = HadronMom.plus() - LeftPlus; |
---|
164 | //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl; |
---|
165 | LeftMom.setPz(0.5*(LeftPlus - LeftMinus)); |
---|
166 | LeftMom.setE (0.5*(LeftPlus + LeftMinus)); |
---|
167 | RightMom.setPz(0.5*(RightPlus - RightMinus)); |
---|
168 | RightMom.setE (0.5*(RightPlus + RightMinus)); |
---|
169 | //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl; |
---|
170 | Left->Set4Momentum(LeftMom); |
---|
171 | Right->Set4Momentum(RightMom); |
---|
172 | Color.push_back(Left); |
---|
173 | AntiColor.push_back(Right); |
---|
174 | } |
---|
175 | |
---|
176 | |
---|
177 | void G4QGSMSplitableHadron::SoftSplitUp() |
---|
178 | { |
---|
179 | //... sample transversal momenta for sea and valence quarks |
---|
180 | G4double phi, pts; |
---|
181 | G4double SumPy = 0.; |
---|
182 | G4double SumPx = 0.; |
---|
183 | G4ThreeVector Pos = GetPosition(); |
---|
184 | G4int nSeaPair = GetSoftCollisionCount()-1; |
---|
185 | |
---|
186 | // here the condition,to ensure viability of splitting, also in cases |
---|
187 | // where difractive excitation occured together with soft scattering. |
---|
188 | // G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus(); |
---|
189 | // G4double Xmin = theMinPz/LightConeMomentum; |
---|
190 | G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() ); |
---|
191 | while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95; |
---|
192 | |
---|
193 | G4int aSeaPair; |
---|
194 | for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++) |
---|
195 | { |
---|
196 | // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2) |
---|
197 | |
---|
198 | G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); |
---|
199 | |
---|
200 | // BuildSeaQuark() determines quark spin, isospin and colour |
---|
201 | // via parton-constructor G4Parton(aPDGCode) |
---|
202 | |
---|
203 | G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair); |
---|
204 | |
---|
205 | // G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl; |
---|
206 | |
---|
207 | // G4cerr << "Parton 1: " |
---|
208 | // << " PDGcode: " << aPDGCode |
---|
209 | // << " - Name: " << aParton->GetDefinition()->GetParticleName() |
---|
210 | // << " - Type: " << aParton->GetDefinition()->GetParticleType() |
---|
211 | // << " - Spin-3: " << aParton->GetSpinZ() |
---|
212 | // << " - Colour: " << aParton->GetColour() << G4endl; |
---|
213 | |
---|
214 | // save colour a spin-3 for anti-quark |
---|
215 | |
---|
216 | G4int firstPartonColour = aParton->GetColour(); |
---|
217 | G4double firstPartonSpinZ = aParton->GetSpinZ(); |
---|
218 | |
---|
219 | SumPx += aParton->Get4Momentum().px(); |
---|
220 | SumPy += aParton->Get4Momentum().py(); |
---|
221 | Color.push_back(aParton); |
---|
222 | |
---|
223 | // create anti-quark |
---|
224 | |
---|
225 | aParton = BuildSeaQuark(true, aPDGCode, nSeaPair); |
---|
226 | aParton->SetSpinZ(-firstPartonSpinZ); |
---|
227 | aParton->SetColour(-firstPartonColour); |
---|
228 | |
---|
229 | // G4cerr << "Parton 2: " |
---|
230 | // << " PDGcode: " << -aPDGCode |
---|
231 | // << " - Name: " << aParton->GetDefinition()->GetParticleName() |
---|
232 | // << " - Type: " << aParton->GetDefinition()->GetParticleType() |
---|
233 | // << " - Spin-3: " << aParton->GetSpinZ() |
---|
234 | // << " - Colour: " << aParton->GetColour() << G4endl; |
---|
235 | // G4cerr << "------------" << G4endl; |
---|
236 | |
---|
237 | SumPx += aParton->Get4Momentum().px(); |
---|
238 | SumPy += aParton->Get4Momentum().py(); |
---|
239 | AntiColor.push_back(aParton); |
---|
240 | } |
---|
241 | // Valence quark |
---|
242 | G4Parton* pColorParton = NULL; |
---|
243 | G4Parton* pAntiColorParton = NULL; |
---|
244 | GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton); |
---|
245 | G4int ColorEncoding = pColorParton->GetPDGcode(); |
---|
246 | G4int AntiColorEncoding = pAntiColorParton->GetPDGcode(); |
---|
247 | |
---|
248 | pts = sigmaPt*std::sqrt(-std::log(G4UniformRand())); |
---|
249 | phi = 2.*pi*G4UniformRand(); |
---|
250 | G4double Px = pts*std::cos(phi); |
---|
251 | G4double Py = pts*std::sin(phi); |
---|
252 | SumPx += Px; |
---|
253 | SumPy += Py; |
---|
254 | |
---|
255 | if (ColorEncoding < 0) // use particle definition |
---|
256 | { |
---|
257 | G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0); |
---|
258 | pColorParton->Set4Momentum(ColorMom); |
---|
259 | G4LorentzVector AntiColorMom(Px, Py, 0, 0); |
---|
260 | pAntiColorParton->Set4Momentum(AntiColorMom); |
---|
261 | } |
---|
262 | else |
---|
263 | { |
---|
264 | G4LorentzVector ColorMom(Px, Py, 0, 0); |
---|
265 | pColorParton->Set4Momentum(ColorMom); |
---|
266 | G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0); |
---|
267 | pAntiColorParton->Set4Momentum(AntiColorMom); |
---|
268 | } |
---|
269 | Color.push_back(pColorParton); |
---|
270 | AntiColor.push_back(pAntiColorParton); |
---|
271 | |
---|
272 | // Sample X |
---|
273 | G4int nAttempt = 0; |
---|
274 | G4double SumX = 0; |
---|
275 | G4double aBeta = beta; |
---|
276 | G4double ColorX, AntiColorX; |
---|
277 | G4double HPWtest = 0; |
---|
278 | if (GetDefinition() == G4PionMinus::PionMinusDefinition()) aBeta = 1.; |
---|
279 | if (GetDefinition() == G4Gamma::GammaDefinition()) aBeta = 1.; |
---|
280 | if (GetDefinition() == G4PionPlus::PionPlusDefinition()) aBeta = 1.; |
---|
281 | if (GetDefinition() == G4PionZero::PionZeroDefinition()) aBeta = 1.; |
---|
282 | if (GetDefinition() == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.; |
---|
283 | if (GetDefinition() == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.; |
---|
284 | do |
---|
285 | { |
---|
286 | SumX = 0; |
---|
287 | nAttempt++; |
---|
288 | G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair; |
---|
289 | G4double beta1 = beta; |
---|
290 | if (std::abs(ColorEncoding) <= 1000 && std::abs(AntiColorEncoding) <= 1000) beta1 = 1.; //... in a meson |
---|
291 | ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
292 | HPWtest = ColorX; |
---|
293 | while (ColorX < Xmin || ColorX > 1.|| 1. - ColorX <= Xmin) {;} |
---|
294 | Color.back()->SetX(SumX = ColorX);// this is the valenz quark. |
---|
295 | for(G4int aPair = 0; aPair < nSeaPair; aPair++) |
---|
296 | { |
---|
297 | NumberOfUnsampledSeaQuarks--; |
---|
298 | ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
299 | Color[aPair]->SetX(ColorX); |
---|
300 | SumX += ColorX; |
---|
301 | NumberOfUnsampledSeaQuarks--; |
---|
302 | AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
303 | AntiColor[aPair]->SetX(AntiColorX); // the 'sea' partons |
---|
304 | SumX += AntiColorX; |
---|
305 | if (1. - SumX <= Xmin) break; |
---|
306 | } |
---|
307 | } |
---|
308 | while (1. - SumX <= Xmin); |
---|
309 | (*(AntiColor.end()-1))->SetX(1. - SumX); // the di-quark takes the rest, then go to momentum |
---|
310 | /// and here is the bug ;-) @@@@@@@@@@@@@ |
---|
311 | if(getenv("debug_QGSMSplitableHadron") )G4cout << "particle energy at split = "<<Get4Momentum().t()<<G4endl; |
---|
312 | G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus()); |
---|
313 | G4double lightCone2 = ((!Direction) ? Get4Momentum().plus() : Get4Momentum().minus()); |
---|
314 | // lightCone -= 0.5*Get4Momentum().m(); |
---|
315 | // hpw testing @@@@@ lightCone = 2.*Get4Momentum().t(); |
---|
316 | if(getenv("debug_QGSMSplitableHadron") )G4cout << "Light cone = "<<lightCone<<G4endl; |
---|
317 | for(aSeaPair = 0; aSeaPair < nSeaPair+1; aSeaPair++) |
---|
318 | { |
---|
319 | G4Parton* aParton = Color[aSeaPair]; |
---|
320 | aParton->DefineMomentumInZ(lightCone, lightCone2, Direction); |
---|
321 | |
---|
322 | aParton = AntiColor[aSeaPair]; |
---|
323 | aParton->DefineMomentumInZ(lightCone, lightCone2, Direction); |
---|
324 | } |
---|
325 | //--DEBUG-- cout <<G4endl<<"XSAMPLE "<<HPWtest<<G4endl; |
---|
326 | return; |
---|
327 | } |
---|
328 | |
---|
329 | void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart, G4Parton *& Parton1, G4Parton *& Parton2) |
---|
330 | { |
---|
331 | // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq. |
---|
332 | G4int aEnd; |
---|
333 | G4int bEnd; |
---|
334 | G4int HadronEncoding = aPart->GetPDGEncoding(); |
---|
335 | if (aPart->GetBaryonNumber() == 0) |
---|
336 | { |
---|
337 | theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd); |
---|
338 | } |
---|
339 | else |
---|
340 | { |
---|
341 | theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd); |
---|
342 | } |
---|
343 | |
---|
344 | Parton1 = new G4Parton(aEnd); |
---|
345 | Parton1->SetPosition(GetPosition()); |
---|
346 | |
---|
347 | // G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl; |
---|
348 | // G4cerr << "Parton 1: " |
---|
349 | // << " PDGcode: " << aEnd |
---|
350 | // << " - Name: " << Parton1->GetDefinition()->GetParticleName() |
---|
351 | // << " - Type: " << Parton1->GetDefinition()->GetParticleType() |
---|
352 | // << " - Spin-3: " << Parton1->GetSpinZ() |
---|
353 | // << " - Colour: " << Parton1->GetColour() << G4endl; |
---|
354 | |
---|
355 | Parton2 = new G4Parton(bEnd); |
---|
356 | Parton2->SetPosition(GetPosition()); |
---|
357 | |
---|
358 | // G4cerr << "Parton 2: " |
---|
359 | // << " PDGcode: " << bEnd |
---|
360 | // << " - Name: " << Parton2->GetDefinition()->GetParticleName() |
---|
361 | // << " - Type: " << Parton2->GetDefinition()->GetParticleType() |
---|
362 | // << " - Spin-3: " << Parton2->GetSpinZ() |
---|
363 | // << " - Colour: " << Parton2->GetColour() << G4endl; |
---|
364 | // G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl; |
---|
365 | |
---|
366 | // colour of parton 1 choosen at random by G4Parton(aEnd) |
---|
367 | // colour of parton 2 is the opposite: |
---|
368 | |
---|
369 | Parton2->SetColour(-(Parton1->GetColour())); |
---|
370 | |
---|
371 | // isospin-3 of both partons is handled by G4Parton(PDGCode) |
---|
372 | |
---|
373 | // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd) |
---|
374 | // spin-3 of parton 2 may be constrained by spin of original particle: |
---|
375 | |
---|
376 | if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin()) |
---|
377 | { |
---|
378 | Parton2->SetSpinZ(-(Parton2->GetSpinZ())); |
---|
379 | } |
---|
380 | |
---|
381 | // G4cerr << "Parton 2: " |
---|
382 | // << " PDGcode: " << bEnd |
---|
383 | // << " - Name: " << Parton2->GetDefinition()->GetParticleName() |
---|
384 | // << " - Type: " << Parton2->GetDefinition()->GetParticleType() |
---|
385 | // << " - Spin-3: " << Parton2->GetSpinZ() |
---|
386 | // << " - Colour: " << Parton2->GetColour() << G4endl; |
---|
387 | // G4cerr << "------------" << G4endl; |
---|
388 | |
---|
389 | } |
---|
390 | |
---|
391 | |
---|
392 | G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare) |
---|
393 | { |
---|
394 | G4double R; |
---|
395 | while((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare) {;} |
---|
396 | R = std::sqrt(R); |
---|
397 | G4double phi = twopi*G4UniformRand(); |
---|
398 | return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.); |
---|
399 | } |
---|
400 | |
---|
401 | G4Parton * G4QGSMSplitableHadron:: |
---|
402 | BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/) |
---|
403 | { |
---|
404 | if (isAntiQuark) aPDGCode*=-1; |
---|
405 | G4Parton* result = new G4Parton(aPDGCode); |
---|
406 | result->SetPosition(GetPosition()); |
---|
407 | G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX); |
---|
408 | G4LorentzVector a4Momentum(aPtVector, 0); |
---|
409 | result->Set4Momentum(a4Momentum); |
---|
410 | return result; |
---|
411 | } |
---|
412 | |
---|
413 | G4double G4QGSMSplitableHadron:: |
---|
414 | SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta) |
---|
415 | { |
---|
416 | G4double result; |
---|
417 | G4double x1, x2; |
---|
418 | G4double ymax = 0; |
---|
419 | for(G4int ii=1; ii<100; ii++) |
---|
420 | { |
---|
421 | G4double y = std::pow(1./G4double(ii), alpha); |
---|
422 | y *= std::pow( std::pow(1-anXmin-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea); |
---|
423 | y *= std::pow(1-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); |
---|
424 | if(y>ymax) ymax = y; |
---|
425 | } |
---|
426 | G4double y; |
---|
427 | G4double xMax=1-(totalSea+1)*anXmin; |
---|
428 | if(anXmin > xMax) |
---|
429 | { |
---|
430 | G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl; |
---|
431 | throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints."); |
---|
432 | } |
---|
433 | do |
---|
434 | { |
---|
435 | x1 = CLHEP::RandFlat::shoot(anXmin, xMax); |
---|
436 | y = std::pow(x1, alpha); |
---|
437 | y *= std::pow( std::pow(1-x1-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea); |
---|
438 | y *= std::pow(1-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); |
---|
439 | x2 = ymax*G4UniformRand(); |
---|
440 | } |
---|
441 | while(x2>y); |
---|
442 | result = x1; |
---|
443 | return result; |
---|
444 | } |
---|