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