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