| 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 // $Id: G4DiffractiveSplitableHadron.cc 74627 2013-10-17 07:04:38Z gcosmo $
>> 28 // GEANT4 tag $Name: $
27 // 29 //
28 30
29 // ------------------------------------------- 31 // ------------------------------------------------------------
30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file
31 // 33 //
32 // ---------------- G4DiffractiveSplitabl 34 // ---------------- G4DiffractiveSplitableHadron----------------
33 // by Gunter Folger, August 1998. 35 // by Gunter Folger, August 1998.
34 // class splitting an interacting partic 36 // class splitting an interacting particle. Used by FTF String Model.
35 // ------------------------------------------- 37 // ------------------------------------------------------------
36 38
37 #include "G4DiffractiveSplitableHadron.hh" 39 #include "G4DiffractiveSplitableHadron.hh"
38 40
39 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh"
40 #include "Randomize.hh" 42 #include "Randomize.hh"
41 43
42 44
43 //============================================ 45 //============================================================================
44 46
45 G4DiffractiveSplitableHadron::G4DiffractiveSpl << 47 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron() {
46 { <<
47 PartonIndex = -1; 48 PartonIndex = -1;
48 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. <<
49 Parton[0] = new G4Parton( 1 ); 49 Parton[0] = new G4Parton( 1 );
50 Parton[1] = new G4Parton(-1 ); 50 Parton[1] = new G4Parton(-1 );
51 <<
52 Parton[0]->Set4Momentum(tmp); Parton[1]->Set <<
53 } 51 }
54 52
55 53
56 //============================================ 54 //============================================================================
57 55
58 G4DiffractiveSplitableHadron::G4DiffractiveSpl 56 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron( const G4ReactionProduct& aPrimary ) :
59 G4VSplitableHadron( aPrimary ) 57 G4VSplitableHadron( aPrimary )
60 { 58 {
61 PartonIndex = -1; << 59 PartonIndex = -2;
62 Parton[0] = nullptr; << 60 Parton[0] = NULL;
63 Parton[1] = nullptr; <<
64 } 61 }
65 62
66 63
67 //============================================ 64 //============================================================================
68 65
69 G4DiffractiveSplitableHadron::G4DiffractiveSpl 66 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron( const G4Nucleon& aNucleon ) :
70 G4VSplitableHadron( aNucleon ) 67 G4VSplitableHadron( aNucleon )
71 { 68 {
72 PartonIndex = -1; << 69 PartonIndex = -2;
73 Parton[0] = nullptr; << 70 Parton[0] = NULL;
74 Parton[1] = nullptr; <<
75 } 71 }
76 72
77 73
78 //============================================ 74 //============================================================================
79 75
80 G4DiffractiveSplitableHadron::G4DiffractiveSpl 76 G4DiffractiveSplitableHadron::G4DiffractiveSplitableHadron( const G4VKineticNucleon* aNucleon ) :
81 G4VSplitableHadron( aNucleon ) 77 G4VSplitableHadron( aNucleon )
82 { 78 {
83 PartonIndex = -1; << 79 PartonIndex = -2;
84 Parton[0] = nullptr; << 80 Parton[0] = NULL;
85 Parton[1] = nullptr; <<
86 } 81 }
87 82
88 83
89 //============================================ 84 //============================================================================
90 85
91 G4DiffractiveSplitableHadron::~G4DiffractiveSp << 86 G4DiffractiveSplitableHadron::~G4DiffractiveSplitableHadron() {
>> 87 //G4cout << "Destruct G4DiffractiveSplitableHadron" << Parton[0] << " " << Parton[1] << G4endl;
>> 88 //if ( Parton[0] != NULL ) { delete Parton[0]; delete Parton[1]; }
>> 89 }
92 90
93 91
94 //============================================ 92 //============================================================================
95 93
96 void G4DiffractiveSplitableHadron::SplitUp() { 94 void G4DiffractiveSplitableHadron::SplitUp() {
97 << 95 //G4cout << "SplitUp() IsSplit() Parton[0] " << IsSplit() << " " << Parton[0] << G4endl;
98 if ( IsSplit() ) return; 96 if ( IsSplit() ) return;
99 Splitting(); 97 Splitting();
100 // Split once only... 98 // Split once only...
101 if ( Parton[0] != nullptr ) return; << 99 if ( Parton[0] != NULL ) return;
102 100
103 // flavours of quark ends 101 // flavours of quark ends
104 G4int PDGcode = GetDefinition()->GetPDGEncod 102 G4int PDGcode = GetDefinition()->GetPDGEncoding();
105 G4int stringStart, stringEnd; 103 G4int stringStart, stringEnd;
106 ChooseStringEnds( PDGcode, &stringStart, &st 104 ChooseStringEnds( PDGcode, &stringStart, &stringEnd );
107 105
108 Parton[0] = new G4Parton( stringStart ); 106 Parton[0] = new G4Parton( stringStart );
109 Parton[1] = new G4Parton( stringEnd ); 107 Parton[1] = new G4Parton( stringEnd );
110 <<
111 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. <<
112 Parton[0]->Set4Momentum(tmp); Parton[1]->Set <<
113 <<
114 /* // <<
115 if ( G4UniformRand() < 1.75 ) { //0.75 <<
116 Parton[0] = new G4Parton( stringStart ); <<
117 Parton[1] = new G4Parton( stringEnd ); <<
118 } else { <<
119 Parton[0] = new G4Parton( stringEnd ); <<
120 Parton[1] = new G4Parton( stringStart ); <<
121 } <<
122 */ <<
123 <<
124 PartonIndex = -1; 108 PartonIndex = -1;
125 } 109 }
126 110
127 111
128 //============================================ 112 //============================================================================
129 113
130 G4Parton* G4DiffractiveSplitableHadron::GetNex 114 G4Parton* G4DiffractiveSplitableHadron::GetNextParton() {
131 ++PartonIndex; 115 ++PartonIndex;
132 if ( PartonIndex > 1 || PartonIndex < 0 ) << 116 if ( PartonIndex > 1 || PartonIndex < 0 ) return NULL;
133 G4int PartonInd( PartonIndex ); 117 G4int PartonInd( PartonIndex );
134 if ( PartonIndex == 1 ) PartonIndex = -1; 118 if ( PartonIndex == 1 ) PartonIndex = -1;
135 return Parton[ PartonInd ]; 119 return Parton[ PartonInd ];
136 } 120 }
137 121
138 122
139 //============================================ 123 //============================================================================
140 124
141 G4Parton* G4DiffractiveSplitableHadron::GetNex 125 G4Parton* G4DiffractiveSplitableHadron::GetNextAntiParton() {
142 ++PartonIndex; 126 ++PartonIndex;
143 if ( PartonIndex > 1 || PartonIndex < 0 ) << 127 if ( PartonIndex > 1 || PartonIndex < 0 ) return NULL;
144 G4int PartonInd( PartonIndex ); 128 G4int PartonInd( PartonIndex );
145 if ( PartonIndex == 1 ) PartonIndex = -1; 129 if ( PartonIndex == 1 ) PartonIndex = -1;
146 return Parton[ PartonInd ]; 130 return Parton[ PartonInd ];
147 } 131 }
148 132
149 133
150 //============================================ 134 //============================================================================
151 135
152 void G4DiffractiveSplitableHadron::SetFirstPar 136 void G4DiffractiveSplitableHadron::SetFirstParton( G4int PDGcode ) {
153 delete Parton[0]; 137 delete Parton[0];
154 Parton[0] = new G4Parton( PDGcode ); 138 Parton[0] = new G4Parton( PDGcode );
155 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. <<
156 Parton[0]->Set4Momentum(tmp); <<
157 } 139 }
158 140
159 141
160 //============================================ 142 //============================================================================
161 143
162 void G4DiffractiveSplitableHadron::SetSecondPa 144 void G4DiffractiveSplitableHadron::SetSecondParton( G4int PDGcode ) {
163 delete Parton[1]; 145 delete Parton[1];
164 Parton[1] = new G4Parton( PDGcode ); 146 Parton[1] = new G4Parton( PDGcode );
165 G4LorentzVector tmp=G4LorentzVector(0.,0.,0. <<
166 Parton[1]->Set4Momentum(tmp); <<
167 } 147 }
168 148
169 149
170 //============================================ 150 //============================================================================
171 151
172 void G4DiffractiveSplitableHadron::ChooseStrin 152 void G4DiffractiveSplitableHadron::ChooseStringEnds( G4int PDGcode, G4int* aEnd,
173 153 G4int* bEnd ) const {
>> 154 const G4double udspin1 = 1.0/6.0;
>> 155 const G4double uuspin1 = 1.0/3.0;
>> 156 //const G4double udspin0 = 1.0/2.0;
>> 157
174 G4int absPDGcode = std::abs( PDGcode ); 158 G4int absPDGcode = std::abs( PDGcode );
175 159
176 if ( absPDGcode < 1000 ) { //-------------- 160 if ( absPDGcode < 1000 ) { //-------------------- Meson -------------
177 G4int heavy(0), light(0); << 161 G4int heavy = absPDGcode/100;
178 if (!((absPDGcode == 111)||(absPDGcode == << 162 G4int light = (absPDGcode % 100)/10;
179 { // Ordinary mes << 163 //G4int anti = std::pow( -1 , std::max( heavy, light ) );
180 heavy = absPDGcode/100; << 164 G4int anti = 1 - 2*( std::max( heavy, light ) % 2 );
181 light = (absPDGcode % 100)/10; << 165 if (PDGcode < 0 ) anti *= -1;
182 //G4int anti = std::pow( -1 , std::max( h << 166 heavy *= anti;
183 G4int anti = 1 - 2*( std::max( heavy, lig << 167 light *= -1 * anti;
184 if (PDGcode < 0 ) anti *= -1; <<
185 heavy *= anti; <<
186 light *= -1 * anti; <<
187 } <<
188 else <<
189 { // Pi0, Eta, Eta <<
190 if ( G4UniformRand() < 0.5 ) {heavy = 1; <<
191 else {heavy = 2; <<
192 } <<
193 if ( G4UniformRand() < 0.5 ) { 168 if ( G4UniformRand() < 0.5 ) {
194 *aEnd = heavy; 169 *aEnd = heavy;
195 *bEnd = light; 170 *bEnd = light;
196 } else { 171 } else {
197 *aEnd = light; 172 *aEnd = light;
198 *bEnd = heavy; 173 *bEnd = heavy;
199 } 174 }
200 } else { //-------------- 175 } else { //-------------------- Baryon --------------
201 G4int j1000 = PDGcode/1000; 176 G4int j1000 = PDGcode/1000;
202 G4int j100 = (PDGcode % 1000)/100; 177 G4int j100 = (PDGcode % 1000)/100;
203 G4int j10 = (PDGcode % 100)/10; 178 G4int j10 = (PDGcode % 100)/10;
204 << 179 G4double random = G4UniformRand();
205 if ( absPDGcode > 4000 ) { << 180 if ( std::abs( j100 ) >= std::abs( j10 ) ) {
206 *aEnd = j10; << 181 if ( random < udspin1 ) {
207 if ( G4UniformRand() > 0.25 ) { << 182 *aEnd = j1000;
208 *bEnd = Diquark( j1000, j100, 0 ); << 183 *bEnd = Diquark( j100, j10, 1 );
209 } else { << 184 } else if ( random < udspin1 + uuspin1 ) {
>> 185 *aEnd = j10;
210 *bEnd = Diquark( j1000, j100, 1 ); 186 *bEnd = Diquark( j1000, j100, 1 );
>> 187 } else {
>> 188 *aEnd = j100;
>> 189 // Careful, there is no diquark q1=q2, (q1 q2)0, possible for Omega-
>> 190 *bEnd = Diquark( j1000, j10, j1000 != j100 ? 0 : 1 );
211 } 191 }
212 return; << 192 } else {
213 } << 193 // Lambda-like hadrons have two lightest quarks in spin 0
214 << 194 if ( random < udspin1 ) {
215 G4double SuppresUUDDSS=1.0/2.0; <<
216 if ((j1000 == j100) && (j1000 == j10)) Sup <<
217 <<
218 const G4int maxNumberOfLoops = 1000; <<
219 G4int loopCounter = 0; <<
220 do <<
221 { <<
222 G4double random = G4UniformRand(); <<
223 <<
224 if (random < 0.33333) <<
225 { <<
226 if (( j100 == j10 ) && ( G4UniformRand <<
227 *aEnd = j1000; 195 *aEnd = j1000;
228 if ( j100 == j10 ) {*bEnd << 196 // as above, but with charmed baryons
229 else << 197 *bEnd = Diquark( j100, j10, j100 != j10 ? 0 : 10 );
230 if ( G4UniformRand() > 0.25) {*bEnd << 198 } else if ( random < udspin1 + uuspin1 ) {
231 else {*bEnd = <<
232 break; <<
233 } <<
234 else if (random < 0.66667) <<
235 { <<
236 if (( j1000 == j10 ) && ( G4UniformRan <<
237 *aEnd = j100; <<
238 if ( j1000 == j10 ) {*bEnd <<
239 else <<
240 if ( G4UniformRand() > 0.25) {*bEnd <<
241 else {*bEnd = <<
242 break; <<
243 } <<
244 else <<
245 { <<
246 if (( j1000 == j100 ) && ( G4UniformRa <<
247 *aEnd = j10; 199 *aEnd = j10;
248 if ( j1000 == j100 ) {*bEnd << 200 *bEnd = Diquark( j1000, j100, 1 );
249 else << 201 } else {
250 if ( G4UniformRand() > 0.25) {*bEnd << 202 *aEnd = j100;
251 else {*bEnd = << 203 *bEnd = Diquark( j1000, j10, 1 );
252 break; << 204 }
253 } <<
254 } while ( (true) && <<
255 ++loopCounter < maxNumberOfLoops <<
256 if ( loopCounter >= maxNumberOfLoops ) { <<
257 *aEnd = j10; *bEnd = Diquark( j1000, j10 <<
258 } 205 }
259 <<
260 } 206 }
>> 207
261 } 208 }
262 209
263 210
264 //============================================ 211 //============================================================================
265 212
266 G4int G4DiffractiveSplitableHadron::Diquark( G 213 G4int G4DiffractiveSplitableHadron::Diquark( G4int aquark, G4int bquark, G4int Spin) const {
267 G4int diquarkPDG = std::max( std::abs( aquar 214 G4int diquarkPDG = std::max( std::abs( aquark ), std::abs( bquark ) ) * 1000 +
268 std::min( std::abs( aquar 215 std::min( std::abs( aquark ), std::abs( bquark ) ) * 100 +
269 2*Spin + 1; 216 2*Spin + 1;
270 return ( aquark > 0 && bquark > 0 ) ? diqu 217 return ( aquark > 0 && bquark > 0 ) ? diquarkPDG : -1*diquarkPDG;
271 } 218 }
272 <<
273 219