| 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: G4HadronBuilder.cc,v 1.12 2010/11/03 17:11:42 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-04 $
27 // 29 //
28 // ------------------------------------------- 30 // -----------------------------------------------------------------------------
29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
30 // 32 //
31 // History: 33 // History:
32 // Gunter Folger, August/September 34 // Gunter Folger, August/September 2001
33 // Create class; algorithm previ 35 // Create class; algorithm previously in G4VLongitudinalStringDecay.
34 // ------------------------------------------- 36 // -----------------------------------------------------------------------------
35 37
36 #include "G4HadronBuilder.hh" 38 #include "G4HadronBuilder.hh"
37 #include "G4SystemOfUnits.hh" <<
38 #include "Randomize.hh" <<
39 #include "G4HadronicException.hh" 39 #include "G4HadronicException.hh"
>> 40 #include "Randomize.hh"
40 #include "G4ParticleTable.hh" 41 #include "G4ParticleTable.hh"
41 42
42 //#define debug_Hbuilder << 43 G4HadronBuilder::G4HadronBuilder(G4double mesonMix, G4double barionMix,
43 //#define debug_heavyHadrons << 44 std::vector<double> scalarMesonMix,
44 << 45 std::vector<double> vectorMesonMix)
45 G4HadronBuilder::G4HadronBuilder(const std::ve << 46 {
46 const std::vector<G4doubl << 47 mesonSpinMix=mesonMix;
47 const std::vector<G4doubl << 48 barionSpinMix=barionMix;
48 const G4doubl << 49 scalarMesonMixings=scalarMesonMix;
49 { << 50 vectorMesonMixings=vectorMesonMix;
50 mesonSpinMix = mesonMix; <<
51 barionSpinMix = barionMix; <<
52 scalarMesonMixings = scalarMesonMix; <<
53 vectorMesonMixings = vectorMesonMix; <<
54 ProbEta_c = Eta_cProb; <<
55 ProbEta_b = Eta_bProb; <<
56 } 51 }
57 52
58 //-------------------------------------------- <<
59 <<
60 G4ParticleDefinition * G4HadronBuilder::Build( 53 G4ParticleDefinition * G4HadronBuilder::Build(G4ParticleDefinition * black, G4ParticleDefinition * white)
61 { 54 {
>> 55
62 if (black->GetParticleSubType()== "di_quark" 56 if (black->GetParticleSubType()== "di_quark" || white->GetParticleSubType()== "di_quark" ) {
63 // Barion << 57
>> 58 // Barion
64 Spin spin = (G4UniformRand() < barionSpin 59 Spin spin = (G4UniformRand() < barionSpinMix) ? SpinHalf : SpinThreeHalf;
65 return Barion(black,white,spin); 60 return Barion(black,white,spin);
66 } else { << 61
67 // Meson << 62 } else {
68 G4int StrangeQ = 0; << 63
69 if( std::abs(black->GetPDGEncoding()) >= << 64 // Meson
70 if( std::abs(white->GetPDGEncoding()) >= << 65 Spin spin = (G4UniformRand() < mesonSpinMix) ? SpinZero : SpinOne;
71 Spin spin = (G4UniformRand() < meso <<
72 return Meson(black,white,spin); 66 return Meson(black,white,spin);
>> 67
73 } 68 }
74 } 69 }
75 70
76 //-------------------------------------------- 71 //-------------------------------------------------------------------------
77 72
78 G4ParticleDefinition * G4HadronBuilder::BuildL 73 G4ParticleDefinition * G4HadronBuilder::BuildLowSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
79 { 74 {
80 if ( black->GetParticleSubType()== "quark" & 75 if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
81 return Meson(black,white, SpinZero); 76 return Meson(black,white, SpinZero);
82 } else { 77 } else {
83 // will return a SpinThreeHalf << 78 // will return a SpinThreeHalf Barion if all quarks the same
84 return Barion(black,white, SpinHalf); 79 return Barion(black,white, SpinHalf);
85 } 80 }
86 } 81 }
87 82
88 //-------------------------------------------- 83 //-------------------------------------------------------------------------
89 84
90 G4ParticleDefinition * G4HadronBuilder::BuildH 85 G4ParticleDefinition * G4HadronBuilder::BuildHighSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
91 { 86 {
92 if ( black->GetParticleSubType()== "quark" & 87 if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
93 return Meson(black,white, SpinOne); 88 return Meson(black,white, SpinOne);
94 } else { 89 } else {
95 return Barion(black,white,SpinThreeHalf); 90 return Barion(black,white,SpinThreeHalf);
96 } 91 }
97 } 92 }
98 93
99 //-------------------------------------------- 94 //-------------------------------------------------------------------------
100 95
101 G4ParticleDefinition * G4HadronBuilder::Meson( 96 G4ParticleDefinition * G4HadronBuilder::Meson(G4ParticleDefinition * black,
102 G4ParticleDefinition * white, 97 G4ParticleDefinition * white, Spin theSpin)
103 { 98 {
104 #ifdef debug_Hbuilder << 99 #ifdef G4VERBOSE
105 // Verify Input Charge << 100 // Verify Input Charge
106 G4double charge = black->GetPDGCharge( << 101
107 if (std::abs(charge) > 2 || std::abs(3. << 102 G4double charge = black->GetPDGCharge()
108 { << 103 + white->GetPDGCharge();
>> 104 if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent ) // 1.001 to avoid int(.9999) -> 0
>> 105 {
109 G4cerr << " G4HadronBuilder::Build()" << 106 G4cerr << " G4HadronBuilder::Build()" << G4endl;
110 G4cerr << " Invalid total charge foun 107 G4cerr << " Invalid total charge found for on input: "
111 << charge<< G4endl; 108 << charge<< G4endl;
112 G4cerr << " PGDcode input quark1/quar 109 G4cerr << " PGDcode input quark1/quark2 : " <<
113 black->GetPDGEncoding() << " / "<< 110 black->GetPDGEncoding() << " / "<<
114 white->GetPDGEncoding() << G4endl; 111 white->GetPDGEncoding() << G4endl;
115 G4cerr << G4endl; 112 G4cerr << G4endl;
116 } 113 }
117 #endif << 114 #endif
118 << 115
119 G4int id1 = black->GetPDGEncoding(); << 116 G4int id1= black->GetPDGEncoding();
120 G4int id2 = white->GetPDGEncoding(); << 117 G4int id2= white->GetPDGEncoding();
121 << 118 // G4int ifl1= std::max(std::abs(id1), std::abs(id2));
122 // G4int ifl1= std::max(std::abs(id1), <<
123 if ( std::abs(id1) < std::abs(id2) ) 119 if ( std::abs(id1) < std::abs(id2) )
124 { << 120 {
125 G4int xchg = id1; 121 G4int xchg = id1;
126 id1 = id2; 122 id1 = id2;
127 id2 = xchg; 123 id2 = xchg;
128 } << 124 }
129 <<
130 G4int abs_id1 = std::abs(id1); <<
131 125
132 if ( abs_id1 > 5 ) << 126 if (std::abs(id1) > 3 )
133 throw G4HadronicException(__FILE__, __LIN 127 throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Meson : Illegal Quark content as input");
134 128
135 G4int PDGEncoding=0; 129 G4int PDGEncoding=0;
136 130
137 if (id1 + id2 == 0) { 131 if (id1 + id2 == 0) {
138 if ( abs_id1 < 4) { // light quark << 132 G4double rmix = G4UniformRand();
139 G4double rmix = G4UniformRand(); << 133 G4int imix = 2*std::abs(id1) - 1;
140 G4int imix = 2*std::abs(id1) - 1; << 134 if(theSpin == SpinZero) {
141 if (theSpin == SpinZero) { << 135 PDGEncoding = 110*(1 + (G4int)(rmix + scalarMesonMixings[imix - 1])
142 PDGEncoding = 110*(1 + (G4int)(rmix << 136 + (G4int)(rmix + scalarMesonMixings[imix])
143 + (G4int)(rmix << 137 ) + theSpin;
144 ) + theSpin; << 138 } else {
145 } else { << 139 PDGEncoding = 110*(1 + (G4int)(rmix + vectorMesonMixings[imix - 1])
146 PDGEncoding = 110*(1 + (G4int)(rmix << 140 + (G4int)(rmix + vectorMesonMixings[imix])
147 + (G4int)(rmix + vectorMes << 141 ) + theSpin;
148 ) + theSpin; << 142 }
149 } <<
150 } else { // for c and b quarks <<
151 <<
152 PDGEncoding = abs_id1*100 + abs_ <<
153 <<
154 if (PDGEncoding == 440) { <<
155 if ( G4UniformRand() < ProbEta <<
156 PDGEncoding +=1; <<
157 } else { <<
158 PDGEncoding +=3; <<
159 } <<
160 } <<
161 if (PDGEncoding == 550) { <<
162 if ( G4UniformRand() < ProbEta <<
163 PDGEncoding +=1; <<
164 } else { <<
165 PDGEncoding +=3; <<
166 } <<
167 } <<
168 } <<
169 } else { 143 } else {
170 PDGEncoding = 100 * std::abs(id1) + 10 * 144 PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) + theSpin;
171 G4bool IsUp = (std::abs(id1)&1) == 0; // 145 G4bool IsUp = (std::abs(id1)&1) == 0; // quark 1 up type quark (u or c)
172 G4bool IsAnti = id1 < 0; // q << 146 G4bool IsAnti = id1 < 0; // quark 1 is antiquark?
173 if ( (IsUp && IsAnti ) || (!IsUp && !IsAn << 147 if( (IsUp && IsAnti ) || (!IsUp && !IsAnti ) )
>> 148 PDGEncoding = - PDGEncoding;
174 } 149 }
175 << 150
176 // ----------------------------------- << 151
177 // Special treatment for charmed and b <<
178 // no excited charmed or bottom mesons, ther <<
179 // into existing charmed and bottom mesons i <<
180 // we use the corresponding ground state mes <<
181 // else, we prefer to conserve the electric <<
182 #ifdef debug_heavyHadrons <<
183 G4int initialPDGEncoding = PDGEncoding; <<
184 #endif <<
185 if ( std::abs( PDGEncoding ) == 1 <<
186 ( PDGEncoding > 0 ? PDGEncoding = 411 : PD <<
187 else if ( std::abs( PDGEncoding ) == 1 <<
188 ( PDGEncoding > 0 ? PDGEncoding = 421 : PD <<
189 else if ( std::abs( PDGEncoding ) == 4 <<
190 ( PDGEncoding > 0 ? PDGEncoding = 411 : PD <<
191 else if ( std::abs( PDGEncoding ) == 4 <<
192 ( PDGEncoding > 0 ? PDGEncoding = 421 : PD <<
193 else if ( std::abs( PDGEncoding ) == 1 <<
194 ( PDGEncoding > 0 ? PDGEncoding = 411 : PD <<
195 else if ( std::abs( PDGEncoding ) == 1 <<
196 ( PDGEncoding > 0 ? PDGEncoding = 421 : PD <<
197 else if ( std::abs( PDGEncoding ) == 2 <<
198 ( PDGEncoding > 0 ? PDGEncoding = 411 : PD <<
199 else if ( std::abs( PDGEncoding ) == 2 <<
200 ( PDGEncoding > 0 ? PDGEncoding = 421 : PD <<
201 else if ( std::abs( PDGEncoding ) == 4 <<
202 ( PDGEncoding > 0 ? PDGEncoding = 411 : PD <<
203 else if ( std::abs( PDGEncoding ) == 4 <<
204 ( PDGEncoding > 0 ? PDGEncoding = 421 : PD <<
205 else if ( std::abs( PDGEncoding ) == 1 <<
206 ( PDGEncoding > 0 ? PDGEncoding = 431 : PD <<
207 else if ( std::abs( PDGEncoding ) == 4 <<
208 ( PDGEncoding > 0 ? PDGEncoding = 431 : PD <<
209 else if ( std::abs( PDGEncoding ) == 1 <<
210 ( PDGEncoding > 0 ? PDGEncoding = 431 : PD <<
211 else if ( std::abs( PDGEncoding ) == 2 <<
212 ( PDGEncoding > 0 ? PDGEncoding = 431 : PD <<
213 else if ( std::abs( PDGEncoding ) == 4 <<
214 ( PDGEncoding > 0 ? PDGEncoding = 431 : PD <<
215 else if ( std::abs( PDGEncoding ) == <<
216 else if ( std::abs( PDGEncoding ) == <<
217 else if ( std::abs( PDGEncoding ) == <<
218 else if ( std::abs( PDGEncoding ) == <<
219 else if ( std::abs( PDGEncoding ) == <<
220 else if ( std::abs( PDGEncoding ) == <<
221 else if ( std::abs( PDGEncoding ) == 9 <<
222 else if ( std::abs( PDGEncoding ) == 9 <<
223 else if ( std::abs( PDGEncoding ) == 9 <<
224 else if ( std::abs( PDGEncoding ) == <<
225 else if ( std::abs( PDGEncoding ) == <<
226 // Bottom mesons <<
227 else if ( std::abs( PDGEncoding ) == 1 <<
228 ( PDGEncoding > 0 ? PDGEncoding = 511 : PD <<
229 else if ( std::abs( PDGEncoding ) == 1 <<
230 ( PDGEncoding > 0 ? PDGEncoding = 521 : PD <<
231 else if ( std::abs( PDGEncoding ) == 5 <<
232 ( PDGEncoding > 0 ? PDGEncoding = 511 : PD <<
233 else if ( std::abs( PDGEncoding ) == 5 <<
234 ( PDGEncoding > 0 ? PDGEncoding = 521 : PD <<
235 else if ( std::abs( PDGEncoding ) == 1 <<
236 ( PDGEncoding > 0 ? PDGEncoding = 511 : PD <<
237 else if ( std::abs( PDGEncoding ) == 1 <<
238 ( PDGEncoding > 0 ? PDGEncoding = 521 : PD <<
239 else if ( std::abs( PDGEncoding ) == 2 <<
240 ( PDGEncoding > 0 ? PDGEncoding = 511 : PD <<
241 else if ( std::abs( PDGEncoding ) == 2 <<
242 ( PDGEncoding > 0 ? PDGEncoding = 521 : PD <<
243 else if ( std::abs( PDGEncoding ) == 5 <<
244 ( PDGEncoding > 0 ? PDGEncoding = 511 : PD <<
245 else if ( std::abs( PDGEncoding ) == 5 <<
246 ( PDGEncoding > 0 ? PDGEncoding = 521 : PD <<
247 else if ( std::abs( PDGEncoding ) == 1 <<
248 ( PDGEncoding > 0 ? PDGEncoding = 531 : PD <<
249 else if ( std::abs( PDGEncoding ) == 5 <<
250 ( PDGEncoding > 0 ? PDGEncoding = 531 : PD <<
251 else if ( std::abs( PDGEncoding ) == 1 <<
252 ( PDGEncoding > 0 ? PDGEncoding = 531 : PD <<
253 else if ( std::abs( PDGEncoding ) == 2 <<
254 ( PDGEncoding > 0 ? PDGEncoding = 531 : PD <<
255 else if ( std::abs( PDGEncoding ) == 5 <<
256 ( PDGEncoding > 0 ? PDGEncoding = 531 : PD <<
257 else if ( std::abs( PDGEncoding ) == 1 <<
258 ( PDGEncoding > 0 ? PDGEncoding = 541 : PD <<
259 else if ( std::abs( PDGEncoding ) == 5 <<
260 ( PDGEncoding > 0 ? PDGEncoding = 541 : PD <<
261 else if ( std::abs( PDGEncoding ) == 1 <<
262 ( PDGEncoding > 0 ? PDGEncoding = 541 : PD <<
263 else if ( std::abs( PDGEncoding ) == 20543 ) <<
264 ( PDGEncoding > 0 ? PDGEncoding = 541 : PD <<
265 else if ( std::abs( PDGEncoding ) == 5 <<
266 ( PDGEncoding > 0 ? PDGEncoding = 541 : PD <<
267 else if ( std::abs( PDGEncoding ) == <<
268 else if ( std::abs( PDGEncoding ) == <<
269 else if ( std::abs( PDGEncoding ) == <<
270 else if ( std::abs( PDGEncoding ) == <<
271 else if ( std::abs( PDGEncoding ) == <<
272 else if ( std::abs( PDGEncoding ) == <<
273 else if ( std::abs( PDGEncoding ) == <<
274 else if ( std::abs( PDGEncoding ) == <<
275 else if ( std::abs( PDGEncoding ) == <<
276 else if ( std::abs( PDGEncoding ) == <<
277 else if ( std::abs( PDGEncoding ) == <<
278 else if ( std::abs( PDGEncoding ) == <<
279 else if ( std::abs( PDGEncoding ) == <<
280 else if ( std::abs( PDGEncoding ) == <<
281 else if ( std::abs( PDGEncoding ) == <<
282 else if ( std::abs( PDGEncoding ) == <<
283 else if ( std::abs( PDGEncoding ) == <<
284 else if ( std::abs( PDGEncoding ) == 9 <<
285 else if ( std::abs( PDGEncoding ) == 9 <<
286 else if ( std::abs( PDGEncoding ) == <<
287 else if ( std::abs( PDGEncoding ) == <<
288 else if ( std::abs( PDGEncoding ) == <<
289 else if ( std::abs( PDGEncoding ) == <<
290 else if ( std::abs( PDGEncoding ) == <<
291 else if ( std::abs( PDGEncoding ) == <<
292 else if ( std::abs( PDGEncoding ) == <<
293 else if ( std::abs( PDGEncoding ) == <<
294 else if ( std::abs( PDGEncoding ) == <<
295 #ifdef debug_heavyHadrons <<
296 if ( initialPDGEncoding != PDGEncoding ) { <<
297 G4cout << "G4HadronBuilder::Meson : forcin <<
298 << initialPDGEncoding << " into pdgCode=" <<
299 } <<
300 #endif <<
301 // ----------------------------------- <<
302 <<
303 G4ParticleDefinition * MesonDef= 152 G4ParticleDefinition * MesonDef=
304 G4ParticleTable::GetParticleTable()->FindP 153 G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
305 << 154 #ifdef G4VERBOSE
306 #ifdef debug_Hbuilder <<
307 if (MesonDef == 0 ) { 155 if (MesonDef == 0 ) {
308 G4cerr << " G4HadronBuilder - Warning: No 156 G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
309 << PDGEncoding << G4endl; 157 << PDGEncoding << G4endl;
310 } else if ( ( black->GetPDGCharge() + whit 158 } else if ( ( black->GetPDGCharge() + white->GetPDGCharge()
311 - MesonDef->GetPDGCharge() ) > per 159 - MesonDef->GetPDGCharge() ) > perCent ) {
312 G4cerr << " G4HadronBuilder - Warnin 160 G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
313 << " Quark1/2 = " 161 << " Quark1/2 = "
314 << black->GetParticleName() << " / " 162 << black->GetParticleName() << " / "
315 << white->GetParticleName() 163 << white->GetParticleName()
316 << " resulting Hadron " << MesonDef->Get 164 << " resulting Hadron " << MesonDef->GetParticleName()
317 << G4endl; 165 << G4endl;
318 } 166 }
319 #endif << 167 #endif
320 168
321 return MesonDef; 169 return MesonDef;
322 } 170 }
323 171
324 //-------------------------------------------- <<
325 172
326 G4ParticleDefinition * G4HadronBuilder::Barion 173 G4ParticleDefinition * G4HadronBuilder::Barion(G4ParticleDefinition * black,
327 G4ParticleDefinition * white, << 174 G4ParticleDefinition * white,Spin theSpin)
328 { 175 {
329 #ifdef debug_Hbuilder << 176
330 // Verify Input Charge << 177 #ifdef G4VERBOSE
331 G4double charge = black->GetPDGCharge << 178 // Verify Input Charge
332 if (std::abs(charge) > 2 || std::abs(3 << 179 G4double charge = black->GetPDGCharge()
>> 180 + white->GetPDGCharge();
>> 181 if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent )
333 { 182 {
334 G4cerr << " G4HadronBuilder::Build()" << 183 G4cerr << " G4HadronBuilder::Build()" << G4endl;
335 G4cerr << " Invalid total charge foun 184 G4cerr << " Invalid total charge found for on input: "
336 << charge<< G4endl; 185 << charge<< G4endl;
337 G4cerr << " PGDcode input quark1/quar 186 G4cerr << " PGDcode input quark1/quark2 : " <<
338 black->GetPDGEncoding() << " / "<< 187 black->GetPDGEncoding() << " / "<<
339 white->GetPDGEncoding() << G4endl; 188 white->GetPDGEncoding() << G4endl;
340 G4cerr << G4endl; 189 G4cerr << G4endl;
341 } 190 }
342 #endif << 191 #endif
343 << 192 G4int id1= black->GetPDGEncoding();
344 G4int id1 = black->GetPDGEncoding(); << 193 G4int id2= white->GetPDGEncoding();
345 G4int id2 = white->GetPDGEncoding(); <<
346 <<
347 if ( std::abs(id1) < std::abs(id2) ) 194 if ( std::abs(id1) < std::abs(id2) )
348 { << 195 {
349 G4int xchg = id1; 196 G4int xchg = id1;
350 id1 = id2; 197 id1 = id2;
351 id2 = xchg; 198 id2 = xchg;
352 } << 199 }
353 200
354 if (std::abs(id1) < 1000 || std::abs(id2) > << 201 if (std::abs(id1) < 1000 || std::abs(id2) > 3 )
355 throw G4HadronicException(__FILE__, __LIN 202 throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Barion: Illegal quark content as input");
356 203
357 G4int ifl1= std::abs(id1)/1000; 204 G4int ifl1= std::abs(id1)/1000;
358 G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/1 205 G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100;
359 G4int diquarkSpin = std::abs(id1)%10; 206 G4int diquarkSpin = std::abs(id1)%10;
360 G4int ifl3 = id2; 207 G4int ifl3 = id2;
361 if (id1 < 0) 208 if (id1 < 0)
362 { << 209 {
363 ifl1 = - ifl1; 210 ifl1 = - ifl1;
364 ifl2 = - ifl2; 211 ifl2 = - ifl2;
365 } << 212 }
366 //... Construct barion, distinguish Lambda a 213 //... Construct barion, distinguish Lambda and Sigma barions.
367 G4int kfla = std::abs(ifl1); 214 G4int kfla = std::abs(ifl1);
368 G4int kflb = std::abs(ifl2); 215 G4int kflb = std::abs(ifl2);
369 G4int kflc = std::abs(ifl3); 216 G4int kflc = std::abs(ifl3);
370 217
371 G4int kfld = std::max(kfla,kflb); 218 G4int kfld = std::max(kfla,kflb);
372 kfld = std::max(kfld,kflc); 219 kfld = std::max(kfld,kflc);
373 G4int kflf = std::min(kfla,kflb); 220 G4int kflf = std::min(kfla,kflb);
374 kflf = std::min(kflf,kflc); 221 kflf = std::min(kflf,kflc);
375 222
376 G4int kfle = kfla + kflb + kflc - kfld - kfl 223 G4int kfle = kfla + kflb + kflc - kfld - kflf;
377 224
378 //... barion with content uuu or ddd or sss 225 //... barion with content uuu or ddd or sss has always spin = 3/2
379 theSpin = (kfla == kflb && kflb == kflc)? Sp 226 theSpin = (kfla == kflb && kflb == kflc)? SpinThreeHalf : theSpin;
380 227
381 G4int kfll = 0; 228 G4int kfll = 0;
382 if (kfld < 6) { << 229 if(theSpin == SpinHalf && kfld > kfle && kfle > kflf) {
383 if (theSpin == SpinHalf && kfld > kfle && << 230 // Spin J=1/2 and all three quarks different
384 // Spin J=1/2 and all three quar << 231 // Two states exist: (uds -> lambda or sigma0)
385 // Two states exist: (uds -> lam << 232 // - lambda: s(ud)0 s : 3122; ie. reverse the two lighter quarks
386 // - lambda: s(ud)0 s : 3122; << 233 // - sigma0: s(ud)1 s : 3212
387 // - sigma0: s(ud)1 s : 3212 << 234 if(diquarkSpin == 1 ) {
388 if (diquarkSpin == 1 ) { << 235 if ( kfla == kfld) { // heaviest quark in diquark
389 if ( kfla == kfld) { // heaviest << 236 kfll = 1;
390 kfll = 1; << 237 } else {
391 } else { << 238 kfll = (G4int)(0.25 + G4UniformRand());
392 kfll = (G4int)(0.25 + G4UniformR << 239 }
393 } << 240 }
394 } << 241 if(diquarkSpin == 3 && kfla != kfld)
395 if (diquarkSpin == 3 && kfla != kfld) << 242 kfll = (G4int)(0.75 + G4UniformRand());
396 kfll = (G4int)(0.75 + G4UniformRan << 243 }
397 } <<
398 } <<
399 244
400 G4int PDGEncoding; 245 G4int PDGEncoding;
401 if (kfll == 1) 246 if (kfll == 1)
402 PDGEncoding = 1000 * kfld + 100 * kflf + 247 PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin;
403 else 248 else
404 PDGEncoding = 1000 * kfld + 100 * kfle + 249 PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin;
405 250
406 if (id1 < 0) 251 if (id1 < 0)
407 PDGEncoding = -PDGEncoding; 252 PDGEncoding = -PDGEncoding;
408 253
409 // ----------------------------------- << 254
410 // Special treatment for charmed and b <<
411 // neither excited charmed or bottom baryons <<
412 // heavy (c, b) constitutent quarks: <<
413 // Sigma_c* , Xi_c' , Xi_c* , Omega_c* , <<
414 // Xi_cc , Xi_cc* , Omega_cc , Omega_cc* , <<
415 // Sigma_b* , Xi_b' , Xi_b* , Omega_b*, <<
416 // Xi_bc , Xi_bc' , Xi_bc* , Omega_bc , Om <<
417 // Omega_bcc , Omega_bcc* , Xi_bb, Xi_bb* <<
418 // Omega_bbc , Omega_bbc* , Omega_bbb <<
419 // therefore we need to transform these into <<
420 // baryons in Geant4. Whenever possible, we <<
421 // baryons with the same quantum numbers; el <<
422 // electric charge rather than other flavor <<
423 #ifdef debug_heavyHadrons <<
424 G4int charmViolation = 0, bottomViolation = <<
425 G4int initialPDGEncoding = PDGEncoding; <<
426 #endif <<
427 if ( std::abs( PDGEncoding ) == 4 <<
428 ( PDGEncoding > 0 ? PDGEncoding = 4222 : P <<
429 } else if ( std::abs( PDGEncoding ) == 4214 <<
430 ( PDGEncoding > 0 ? PDGEncoding = 4212 : P <<
431 } else if ( std::abs( PDGEncoding ) == 4114 <<
432 ( PDGEncoding > 0 ? PDGEncoding = 4112 : P <<
433 } else if ( std::abs( PDGEncoding ) == 4322 <<
434 ( PDGEncoding > 0 ? PDGEncoding = 4232 : P <<
435 } else if ( std::abs( PDGEncoding ) == 4312 <<
436 ( PDGEncoding > 0 ? PDGEncoding = 4132 : P <<
437 } else if ( std::abs( PDGEncoding ) == 4324 <<
438 ( PDGEncoding > 0 ? PDGEncoding = 4232 : P <<
439 } else if ( std::abs( PDGEncoding ) == 4314 <<
440 ( PDGEncoding > 0 ? PDGEncoding = 4132 : P <<
441 } else if ( std::abs( PDGEncoding ) == 4334 <<
442 ( PDGEncoding > 0 ? PDGEncoding = 4332 : P <<
443 } else if ( std::abs( PDGEncoding ) == 4412 <<
444 ( PDGEncoding > 0 ? PDGEncoding = 4232 : P <<
445 #ifdef debug_heavyHadrons <<
446 charmViolation = 1; <<
447 #endif <<
448 } else if ( std::abs( PDGEncoding ) == <<
449 ( PDGEncoding > 0 ? PDGEncoding = 4222 : P <<
450 #ifdef debug_heavyHadrons <<
451 charmViolation = 1; <<
452 #endif <<
453 } else if ( std::abs( PDGEncoding ) == 4414 <<
454 ( PDGEncoding > 0 ? PDGEncoding = 4232 : P <<
455 #ifdef debug_heavyHadrons <<
456 charmViolation = 1; <<
457 #endif <<
458 } else if ( std::abs( PDGEncoding ) == <<
459 ( PDGEncoding > 0 ? PDGEncoding = 4222 : P <<
460 #ifdef debug_heavyHadrons <<
461 charmViolation = 1; <<
462 #endif <<
463 } else if ( std::abs( PDGEncoding ) == 4432 <<
464 ( PDGEncoding > 0 ? PDGEncoding = 4232 : P <<
465 #ifdef debug_heavyHadrons <<
466 charmViolation = 1; <<
467 #endif <<
468 } else if ( std::abs( PDGEncoding ) == 4434 <<
469 ( PDGEncoding > 0 ? PDGEncoding = 4232 : P <<
470 #ifdef debug_heavyHadrons <<
471 charmViolation = 1; <<
472 #endif <<
473 } else if ( std::abs( PDGEncoding ) == <<
474 ( PDGEncoding > 0 ? PDGEncoding = 4222 : P <<
475 #ifdef debug_heavyHadrons <<
476 charmViolation = 2; <<
477 #endif <<
478 // Bottom baryons <<
479 } else if ( std::abs( PDGEncoding ) == <<
480 ( PDGEncoding > 0 ? PDGEncoding = 5112 : P <<
481 } else if ( std::abs( PDGEncoding ) == <<
482 ( PDGEncoding > 0 ? PDGEncoding = 5212 : P <<
483 } else if ( std::abs( PDGEncoding ) == <<
484 ( PDGEncoding > 0 ? PDGEncoding = 5222 : P <<
485 } else if ( std::abs( PDGEncoding ) == <<
486 ( PDGEncoding > 0 ? PDGEncoding = 5132 : P <<
487 } else if ( std::abs( PDGEncoding ) == <<
488 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
489 } else if ( std::abs( PDGEncoding ) == <<
490 ( PDGEncoding > 0 ? PDGEncoding = 5132 : P <<
491 } else if ( std::abs( PDGEncoding ) == <<
492 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
493 } else if ( std::abs( PDGEncoding ) == <<
494 ( PDGEncoding > 0 ? PDGEncoding = 5332 : P <<
495 } else if ( std::abs( PDGEncoding ) == <<
496 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
497 #ifdef debug_heavyHadrons <<
498 charmViolation = 1; <<
499 #endif <<
500 } else if ( std::abs( PDGEncoding ) == <<
501 ( PDGEncoding > 0 ? PDGEncoding = 5222 : P <<
502 #ifdef debug_heavyHadrons <<
503 charmViolation = 1; <<
504 #endif <<
505 } else if ( std::abs( PDGEncoding ) == <<
506 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
507 #ifdef debug_heavyHadrons <<
508 charmViolation = 1; <<
509 #endif <<
510 } else if ( std::abs( PDGEncoding ) == <<
511 ( PDGEncoding > 0 ? PDGEncoding = 5222 : P <<
512 #ifdef debug_heavyHadrons <<
513 charmViolation = 1; <<
514 #endif <<
515 } else if ( std::abs( PDGEncoding ) == <<
516 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
517 #ifdef debug_heavyHadrons <<
518 charmViolation = 1; <<
519 #endif <<
520 } else if ( std::abs( PDGEncoding ) == <<
521 ( PDGEncoding > 0 ? PDGEncoding = 5222 : P <<
522 #ifdef debug_heavyHadrons <<
523 charmViolation = 1; <<
524 #endif <<
525 } else if ( std::abs( PDGEncoding ) == <<
526 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
527 #ifdef debug_heavyHadrons <<
528 charmViolation = 1; <<
529 #endif <<
530 } else if ( std::abs( PDGEncoding ) == <<
531 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
532 #ifdef debug_heavyHadrons <<
533 charmViolation = 1; <<
534 #endif <<
535 } else if ( std::abs( PDGEncoding ) == <<
536 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
537 #ifdef debug_heavyHadrons <<
538 charmViolation = 1; <<
539 #endif <<
540 } else if ( std::abs( PDGEncoding ) == <<
541 ( PDGEncoding > 0 ? PDGEncoding = 5222 : P <<
542 #ifdef debug_heavyHadrons <<
543 charmViolation = 2; <<
544 #endif <<
545 } else if ( std::abs( PDGEncoding ) == <<
546 ( PDGEncoding > 0 ? PDGEncoding = 5222 : P <<
547 #ifdef debug_heavyHadrons <<
548 charmViolation = 2; <<
549 #endif <<
550 } else if ( std::abs( PDGEncoding ) == <<
551 ( PDGEncoding > 0 ? PDGEncoding = 5132 : P <<
552 #ifdef debug_heavyHadrons <<
553 bottomViolation = 1; <<
554 #endif <<
555 } else if ( std::abs( PDGEncoding ) == <<
556 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
557 #ifdef debug_heavyHadrons <<
558 bottomViolation = 1; <<
559 #endif <<
560 } else if ( std::abs( PDGEncoding ) == <<
561 ( PDGEncoding > 0 ? PDGEncoding = 5132 : P <<
562 #ifdef debug_heavyHadrons <<
563 bottomViolation = 1; <<
564 #endif <<
565 } else if ( std::abs( PDGEncoding ) == <<
566 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
567 #ifdef debug_heavyHadrons <<
568 bottomViolation = 1; <<
569 #endif <<
570 } else if ( std::abs( PDGEncoding ) == <<
571 ( PDGEncoding > 0 ? PDGEncoding = 5332 : P <<
572 #ifdef debug_heavyHadrons <<
573 bottomViolation = 1; <<
574 #endif <<
575 } else if ( std::abs( PDGEncoding ) == <<
576 ( PDGEncoding > 0 ? PDGEncoding = 5332 : P <<
577 #ifdef debug_heavyHadrons <<
578 bottomViolation = 1; <<
579 #endif <<
580 } else if ( std::abs( PDGEncoding ) == <<
581 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
582 #ifdef debug_heavyHadrons <<
583 charmViolation = 1; bottomViolation = 1; <<
584 #endif <<
585 } else if ( std::abs( PDGEncoding ) == <<
586 ( PDGEncoding > 0 ? PDGEncoding = 5232 : P <<
587 #ifdef debug_heavyHadrons <<
588 charmViolation = 1; bottomViolation = 1; <<
589 #endif <<
590 } else if ( std::abs( PDGEncoding ) == <<
591 ( PDGEncoding > 0 ? PDGEncoding = 5332 : P <<
592 #ifdef debug_heavyHadrons <<
593 bottomViolation = 2; <<
594 #endif <<
595 } <<
596 #ifdef debug_heavyHadrons <<
597 if ( initialPDGEncoding != PDGEncoding ) { <<
598 G4cout << "G4HadronBuilder::Barion : forci <<
599 << initialPDGEncoding << " into pdgCode=" <<
600 if ( charmViolation != 0 || bottomViolat <<
601 G4cout << "\t --> VIOLATION of " << ( ch <<
602 << ( charmViolation != 0 && bot <<
603 << ( bottomViolation != 0 ? " BOT <<
604 } <<
605 } <<
606 #endif <<
607 // ----------------------------------- <<
608 <<
609 G4ParticleDefinition * BarionDef= 255 G4ParticleDefinition * BarionDef=
610 G4ParticleTable::GetParticleTable()->FindP 256 G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
611 << 257 #ifdef G4VERBOSE
612 #ifdef debug_Hbuilder <<
613 if (BarionDef == 0 ) { 258 if (BarionDef == 0 ) {
614 G4cerr << " G4HadronBuilder - Warning: No 259 G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
615 << PDGEncoding << G4endl; 260 << PDGEncoding << G4endl;
616 } else if ( ( black->GetPDGCharge() + whit 261 } else if ( ( black->GetPDGCharge() + white->GetPDGCharge()
617 - BarionDef->GetPDGCharge() ) > pe 262 - BarionDef->GetPDGCharge() ) > perCent ) {
618 G4cerr << " G4HadronBuilder - Warnin 263 G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
619 << " DiQuark/Quark = " 264 << " DiQuark/Quark = "
620 << black->GetParticleName() << " / " 265 << black->GetParticleName() << " / "
621 << white->GetParticleName() 266 << white->GetParticleName()
622 << " resulting Hadron " << BarionDef->Ge 267 << " resulting Hadron " << BarionDef->GetParticleName()
623 << G4endl; 268 << G4endl;
624 } 269 }
625 #endif << 270 #endif
626 271
627 return BarionDef; 272 return BarionDef;
628 } 273 }
629 274