| 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" <<
28 #include "G4SystemOfUnits.hh" <<
29 #include "G4ParticleTable.hh" 27 #include "G4ParticleTable.hh"
30 #include "G4PionPlus.hh" 28 #include "G4PionPlus.hh"
31 #include "G4PionMinus.hh" 29 #include "G4PionMinus.hh"
32 #include "G4Gamma.hh" 30 #include "G4Gamma.hh"
33 #include "G4PionZero.hh" 31 #include "G4PionZero.hh"
34 #include "G4KaonPlus.hh" 32 #include "G4KaonPlus.hh"
35 #include "G4KaonMinus.hh" 33 #include "G4KaonMinus.hh"
36 34
37 #include "G4Log.hh" <<
38 #include "G4Pow.hh" <<
39 <<
40 // based on prototype by Maxim Komogorov 35 // based on prototype by Maxim Komogorov
41 // Splitting into methods, and centralizing of 36 // Splitting into methods, and centralizing of model parameters HPW Feb 1999
42 // restructuring HPW Feb 1999 37 // restructuring HPW Feb 1999
43 // fixing bug in the sampling of 'x', HPW Feb 38 // fixing bug in the sampling of 'x', HPW Feb 1999
44 // fixing bug in sampling pz, HPW Feb 1999. 39 // fixing bug in sampling pz, HPW Feb 1999.
45 // Code now also good for p-nucleus scattering 40 // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999.
46 // Using Parton more directly, HPW Feb 1999. 41 // Using Parton more directly, HPW Feb 1999.
47 // Shortening the algorithm for sampling x, HP 42 // Shortening the algorithm for sampling x, HPW Feb 1999.
48 // sampling of x replaced by formula, taking X 43 // 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 44 // logic much clearer now. HPW Feb 1999
50 // Removed the ordering problem. No Direction 45 // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99.
51 // Fixing p-t distributions for scattering of 46 // Fixing p-t distributions for scattering of nuclei.
52 // Separating out parameters. 47 // Separating out parameters.
53 48
54 void G4QGSMSplitableHadron::InitParameters() 49 void G4QGSMSplitableHadron::InitParameters()
55 { 50 {
56 // changing rapidity distribution for all 51 // changing rapidity distribution for all
57 alpha = -0.5; // Note that this number is st 52 alpha = -0.5; // Note that this number is still assumed in the algorithm
58 // needs to be generalized. << 53 // needs to be generalized.
59 // changing rapidity distribution for projec 54 // changing rapidity distribution for projectile like
60 beta = 2.5;// Note that this number is still 55 beta = 2.5;// Note that this number is still assumed in the algorithm
61 // needs to be generalized. << 56 // needs to be generalized.
62 theMinPz = 0.5*G4PionMinus::PionMinus()->Get << 57 theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass();
63 // theMinPz = 0.1*G4PionMinus::PionMinus()- << 58 // theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass();
64 // theMinPz = G4PionMinus::PionMinus()->Get << 59 // theMinPz = G4PionMinus::PionMinus()->GetPDGMass();
65 // as low as possible, otherwise, we have un 60 // as low as possible, otherwise, we have unphysical boundary conditions in the sampling.
66 StrangeSuppress = 0.48; 61 StrangeSuppress = 0.48;
67 sigmaPt = 0.*GeV; // widens eta slightly, if << 62 sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7,
68 // but Maxim's algorithm breaks energy conse << 63 // but Maxim's algorithm breaks energy conservation
69 widthOfPtSquare = 0.5*sqr(GeV); << 64 // to be revised.
>> 65 widthOfPtSquare = 0.01*GeV*GeV;
70 Direction = FALSE; 66 Direction = FALSE;
71 minTransverseMass = 1*keV; 67 minTransverseMass = 1*keV;
72 iP =0;// Color.begin(); <<
73 iAP =0;// AntiColor.begin(); <<
74 } 68 }
75 69
76 G4QGSMSplitableHadron::G4QGSMSplitableHadron() 70 G4QGSMSplitableHadron::G4QGSMSplitableHadron()
77 { 71 {
78 InitParameters(); 72 InitParameters();
79 } 73 }
80 74
81 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 75 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection)
82 :G4VSplitableHadron(aPrimary) << 76 :G4VSplitableHadron(aPrimary)
83 { 77 {
84 InitParameters(); 78 InitParameters();
85 Direction = aDirection; 79 Direction = aDirection;
86 } 80 }
87 81
88 << 82
89 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 83 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary)
90 : G4VSplitableHadron(aPrimary) << 84 : G4VSplitableHadron(aPrimary)
91 { 85 {
92 InitParameters(); 86 InitParameters();
93 } 87 }
94 88
95 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 89 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon)
96 : G4VSplitableHadron(aNucleon) << 90 : G4VSplitableHadron(aNucleon)
97 { 91 {
98 InitParameters(); 92 InitParameters();
99 } 93 }
100 94
101 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c 95 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection)
102 : G4VSplitableHadron(aNucleon) << 96 : G4VSplitableHadron(aNucleon)
103 { 97 {
104 InitParameters(); 98 InitParameters();
105 Direction = aDirection; 99 Direction = aDirection;
106 } 100 }
107 101
108 G4QGSMSplitableHadron::~G4QGSMSplitableHadron( << 102 G4QGSMSplitableHadron::~G4QGSMSplitableHadron(){}
109 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 }
110 109
111 //******************************************** <<
112 110
>> 111 //**************************************************************************************************************************
>> 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();
117 if (Color.size()!=0) return; 117 if (Color.size()!=0) return;
118 if (GetSoftCollisionCount() == 0) // GetSof << 118 if (GetSoftCollisionCount() == 0)
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 //std::cout << "DSU 1 - "<<HadronMom<<std::endl;
138 139
139 G4double maxAvailMomentum2 = sqr(HadronMom.m << 140 // momenta of string ends
140 << 141 G4double pt2 = HadronMom.perp2();
>> 142 G4double transverseMass2 = HadronMom.plus()*HadronMom.minus();
>> 143 G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
141 G4ThreeVector pt(minTransverseMass, minTrans 144 G4ThreeVector pt(minTransverseMass, minTransverseMass, 0);
142 if (maxAvailMomentum2/widthOfPtSquare>0.01) << 145 if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2);
>> 146 //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl;
143 147
144 G4LorentzVector LeftMom(pt, 0.); 148 G4LorentzVector LeftMom(pt, 0.);
145 G4LorentzVector RightMom; 149 G4LorentzVector RightMom;
146 RightMom.setPx(HadronMom.px() - pt.x()); 150 RightMom.setPx(HadronMom.px() - pt.x());
147 RightMom.setPy(HadronMom.py() - pt.y()); 151 RightMom.setPy(HadronMom.py() - pt.y());
>> 152 //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl;
148 153
149 G4double Local1 = HadronMom.minus() + (Right 154 G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus();
150 G4double Local2 = std::sqrt(std::max(0., sqr 155 G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus()));
151 << 156 //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl;
152 if (Direction) Local2 = -Local2; << 157 if (Direction) Local2 = -Local2;
153 G4double RightMinus = 0.5*(Local1 + Local2 158 G4double RightMinus = 0.5*(Local1 + Local2);
154 G4double LeftMinus = HadronMom.minus() - Rig 159 G4double LeftMinus = HadronMom.minus() - RightMinus;
155 << 160 //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl;
156 if (LeftMinus <= 0.) { <<
157 RightMinus = 0.5*(Local1 - Local2); <<
158 LeftMinus = HadronMom.minus() - RightMinus <<
159 } <<
160 161
161 G4double LeftPlus = LeftMom.perp2()/LeftMin 162 G4double LeftPlus = LeftMom.perp2()/LeftMinus;
162 G4double RightPlus = HadronMom.plus() - Left 163 G4double RightPlus = HadronMom.plus() - LeftPlus;
163 << 164 //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl;
164 LeftMom.setPz(0.5*(LeftPlus - LeftMinus)); 165 LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
165 LeftMom.setE (0.5*(LeftPlus + LeftMinus)); 166 LeftMom.setE (0.5*(LeftPlus + LeftMinus));
166 RightMom.setPz(0.5*(RightPlus - RightMinus)) 167 RightMom.setPz(0.5*(RightPlus - RightMinus));
167 RightMom.setE (0.5*(RightPlus + RightMinus)) 168 RightMom.setE (0.5*(RightPlus + RightMinus));
168 << 169 //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl;
169 Left->Set4Momentum(LeftMom); 170 Left->Set4Momentum(LeftMom);
170 Right->Set4Momentum(RightMom); 171 Right->Set4Momentum(RightMom);
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; <<
175 } 174 }
176 175
177 176
178 void G4QGSMSplitableHadron::SoftSplitUp() 177 void G4QGSMSplitableHadron::SoftSplitUp()
179 { 178 {
180 G4int nSeaPair = GetSoftCollisionCount()-1; << 179 //... sample transversal momenta for sea and valence quarks
181 << 180 G4double phi, pts;
182 G4LorentzVector tmp(0., 0., 0., 0.); << 181 G4double SumPy = 0.;
183 << 182 G4double SumPx = 0.;
184 G4int aSeaPair; << 183 G4ThreeVector Pos = GetPosition();
185 for (aSeaPair = 0; aSeaPair < nSeaPair; aSea << 184 G4int nSeaPair = GetSoftCollisionCount()-1;
186 { << 185
187 // choose quark flavour, d:u:s = 1:1:(1/S << 186 // here the condition,to ensure viability of splitting, also in cases
188 G4int aPDGCode = 1 + (G4int)(G4UniformRand << 187 // where difractive excitation occured together with soft scattering.
189 << 188 // G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus();
190 // BuildSeaQuark() determines quark spin, << 189 // G4double Xmin = theMinPz/LightConeMomentum;
191 // via parton-constructor G4Parton(aPDGCo << 190 G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() );
192 G4Parton * aParton = BuildSeaQuark(false, << 191 while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95;
193 << 192
194 G4int firstPartonColour = aParton->GetColo << 193 G4int aSeaPair;
195 G4double firstPartonSpinZ = aParton->GetSp << 194 for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
196 << 195 {
197 aParton->Set4Momentum(tmp); << 196 // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
198 Color.push_back(aParton); << 197
199 << 198 G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
200 // create anti-quark << 199
201 aParton = BuildSeaQuark(true, aPDGCode, nS << 200 // BuildSeaQuark() determines quark spin, isospin and colour
202 aParton->SetSpinZ(-firstPartonSpinZ); << 201 // via parton-constructor G4Parton(aPDGCode)
203 aParton->SetColour(-firstPartonColour); << 202
204 AntiColor.push_back(aParton); << 203 G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
205 } << 204
206 << 205 // G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl;
207 // Valence quark << 206
208 G4Parton* pColorParton = NULL; << 207 // G4cerr << "Parton 1: "
209 G4Parton* pAntiColorParton = NULL; << 208 // << " PDGcode: " << aPDGCode
210 GetValenceQuarkFlavors(GetDefinition(), pCol << 209 // << " - Name: " << aParton->GetDefinition()->GetParticleName()
211 << 210 // << " - Type: " << aParton->GetDefinition()->GetParticleType()
212 pColorParton->Set4Momentum(tmp); << 211 // << " - Spin-3: " << aParton->GetSpinZ()
213 pAntiColorParton->Set4Momentum(tmp); << 212 // << " - Colour: " << aParton->GetColour() << G4endl;
214 << 213
215 Color.push_back(pColorParton); << 214 // save colour a spin-3 for anti-quark
216 AntiColor.push_back(pAntiColorParton); << 215
217 << 216 G4int firstPartonColour = aParton->GetColour();
218 iP=0; iAP=0; << 217 G4double firstPartonSpinZ = aParton->GetSpinZ();
219 << 218
220 return; << 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;
221 } 327 }
222 328
223 << 329 void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart, G4Parton *& Parton1, G4Parton *& Parton2)
224 void G4QGSMSplitableHadron::GetValenceQuarkFla <<
225 <<
226 { 330 {
227 // Note! convention aEnd = q or (qq)bar and << 331 // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
228 G4int aEnd=0; << 332 G4int aEnd;
229 G4int bEnd=0; << 333 G4int bEnd;
230 G4int HadronEncoding = aPart->GetPDGEncoding 334 G4int HadronEncoding = aPart->GetPDGEncoding();
231 if (aPart->GetBaryonNumber() == 0) << 335 if (aPart->GetBaryonNumber() == 0)
232 { 336 {
233 theMesonSplitter.SplitMeson(HadronEncoding << 337 theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd);
234 } 338 }
235 else 339 else
236 { 340 {
237 theBaryonSplitter.SplitBarion(HadronEncodi << 341 theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd);
238 } << 342 }
239 343
240 Parton1 = new G4Parton(aEnd); 344 Parton1 = new G4Parton(aEnd);
241 Parton1->SetPosition(GetPosition()); 345 Parton1->SetPosition(GetPosition());
242 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
243 Parton2 = new G4Parton(bEnd); 355 Parton2 = new G4Parton(bEnd);
244 Parton2->SetPosition(GetPosition()); 356 Parton2->SetPosition(GetPosition());
245 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
246 // colour of parton 1 choosen at random by G 366 // colour of parton 1 choosen at random by G4Parton(aEnd)
247 // colour of parton 2 is the opposite: 367 // colour of parton 2 is the opposite:
248 368
249 Parton2->SetColour(-(Parton1->GetColour())); 369 Parton2->SetColour(-(Parton1->GetColour()));
250 370
251 // isospin-3 of both partons is handled by G << 371 // isospin-3 of both partons is handled by G4Parton(PDGCode)
252 372
253 // spin-3 of parton 1 and 2 choosen at rando << 373 // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd)
254 // spin-3 of parton 2 may be constrained by << 374 // spin-3 of parton 2 may be constrained by spin of original particle:
255 375
256 if ( std::abs(Parton1->GetSpinZ() + Parton2- << 376 if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin())
257 { 377 {
258 Parton2->SetSpinZ(-(Parton2->GetSpinZ())); << 378 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
259 } << 379 }
260 } <<
261 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;
262 388
>> 389 }
>> 390
>> 391
263 G4ThreeVector G4QGSMSplitableHadron::GaussianP 392 G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
264 { 393 {
265 G4double R; 394 G4double R;
266 const G4int maxNumberOfLoops = 1000; << 395 while((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare) {;}
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); 396 R = std::sqrt(R);
274 G4double phi = twopi*G4UniformRand(); 397 G4double phi = twopi*G4UniformRand();
275 return G4ThreeVector (R*std::cos(phi), R*std << 398 return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.);
276 } 399 }
277 400
278 <<
279 G4Parton * G4QGSMSplitableHadron:: 401 G4Parton * G4QGSMSplitableHadron::
280 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCo 402 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/)
281 { 403 {
282 if (isAntiQuark) aPDGCode*=-1; 404 if (isAntiQuark) aPDGCode*=-1;
283 G4Parton* result = new G4Parton(aPDGCode); << 405 G4Parton* result = new G4Parton(aPDGCode);
284 result->SetPosition(GetPosition()); 406 result->SetPosition(GetPosition());
285 G4ThreeVector aPtVector = GaussianPt(sigmaPt 407 G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
286 G4LorentzVector a4Momentum(aPtVector, 0); 408 G4LorentzVector a4Momentum(aPtVector, 0);
287 result->Set4Momentum(a4Momentum); 409 result->Set4Momentum(a4Momentum);
288 return result; 410 return result;
289 } 411 }
290 412
291 <<
292 G4double G4QGSMSplitableHadron:: 413 G4double G4QGSMSplitableHadron::
293 SampleX(G4double anXmin, G4int nSea, G4int tot 414 SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
294 { 415 {
295 G4double result; 416 G4double result;
296 G4double x1, x2; 417 G4double x1, x2;
297 G4double ymax = 0; 418 G4double ymax = 0;
298 for(G4int ii=1; ii<100; ii++) 419 for(G4int ii=1; ii<100; ii++)
299 { 420 {
300 G4double y = G4Pow::GetInstance()->powA(1. << 421 G4double y = std::pow(1./G4double(ii), alpha);
301 y *= G4Pow::GetInstance()->powN( G4Pow::Ge << 422 y *= std::pow( std::pow(1-anXmin-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
302 G4Pow::Ge << 423 y *= std::pow(1-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
303 y *= G4Pow::GetInstance()->powA(1-anXmin-t << 424 if(y>ymax) ymax = y;
304 G4Pow::GetInstance()->powA(anXmin, aB <<
305 if (y>ymax) ymax = y; <<
306 } 425 }
307 G4double y; 426 G4double y;
308 G4double xMax=1-(totalSea+1)*anXmin; 427 G4double xMax=1-(totalSea+1)*anXmin;
309 if (anXmin > xMax) << 428 if(anXmin > xMax)
310 { 429 {
311 throw G4HadronicException(__FILE__, __LINE << 430 G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl;
312 "G4QGSMSplitableHadron - Fatal: Ca << 431 throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints.");
313 } 432 }
314 const G4int maxNumberOfLoops = 1000; <<
315 G4int loopCounter = 0; <<
316 do 433 do
317 { 434 {
318 x1 = G4RandFlat::shoot(anXmin, xMax); << 435 x1 = CLHEP::RandFlat::shoot(anXmin, xMax);
319 y = G4Pow::GetInstance()->powA(x1, alpha); << 436 y = std::pow(x1, alpha);
320 y *= G4Pow::GetInstance()->powN( G4Pow::Ge << 437 y *= std::pow( std::pow(1-x1-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
321 G4Pow::Ge << 438 y *= std::pow(1-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
322 y *= G4Pow::GetInstance()->powA(1-x1-total <<
323 G4Pow::GetInstance()->powA(anXmin, aB <<
324 x2 = ymax*G4UniformRand(); 439 x2 = ymax*G4UniformRand();
325 } while( (x2>y) && ++loopCounter < maxNumber <<
326 if ( loopCounter >= maxNumberOfLoops ) { <<
327 x1 = 0.5*( anXmin + xMax ); // Just an ac <<
328 } 440 }
>> 441 while(x2>y);
329 result = x1; 442 result = x1;
330 return result; << 443 return result;
331 } 444 }
332 445