Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4HadronBuilder.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//
28// -----------------------------------------------------------------------------
29// GEANT 4 class implementation file
30//
31// History:
32// Gunter Folger, August/September 2001
33// Create class; algorithm previously in G4VLongitudinalStringDecay.
34// -----------------------------------------------------------------------------
35
36#include "G4HadronBuilder.hh"
37#include "G4SystemOfUnits.hh"
38#include "Randomize.hh"
39#include "G4HadronicException.hh"
40#include "G4ParticleTable.hh"
41
42//#define debug_Hbuilder
43//#define debug_heavyHadrons
44
45G4HadronBuilder::G4HadronBuilder(const std::vector<G4double> & mesonMix, const G4double barionMix,
46 const std::vector<G4double> & scalarMesonMix,
47 const std::vector<G4double> & vectorMesonMix,
48 const G4double Eta_cProb, const G4double Eta_bProb)
49{
50 mesonSpinMix = mesonMix;
51 barionSpinMix = barionMix;
52 scalarMesonMixings = scalarMesonMix;
53 vectorMesonMixings = vectorMesonMix;
54 ProbEta_c = Eta_cProb;
55 ProbEta_b = Eta_bProb;
56}
57
58//-------------------------------------------------------------------------
59
60G4ParticleDefinition * G4HadronBuilder::Build(G4ParticleDefinition * black, G4ParticleDefinition * white)
61{
62 if (black->GetParticleSubType()== "di_quark" || white->GetParticleSubType()== "di_quark" ) {
63 // Barion
64 Spin spin = (G4UniformRand() < barionSpinMix) ? SpinHalf : SpinThreeHalf;
65 return Barion(black,white,spin);
66 } else {
67 // Meson
68 G4int StrangeQ = 0;
69 if( std::abs(black->GetPDGEncoding()) >= 3 ) StrangeQ++;
70 if( std::abs(white->GetPDGEncoding()) >= 3 ) StrangeQ++;
71 Spin spin = (G4UniformRand() < mesonSpinMix[StrangeQ]) ? SpinZero : SpinOne;
72 return Meson(black,white,spin);
73 }
74}
75
76//-------------------------------------------------------------------------
77
78G4ParticleDefinition * G4HadronBuilder::BuildLowSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
79{
80 if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
81 return Meson(black,white, SpinZero);
82 } else {
83 // will return a SpinThreeHalf Barion if all quarks the same
84 return Barion(black,white, SpinHalf);
85 }
86}
87
88//-------------------------------------------------------------------------
89
90G4ParticleDefinition * G4HadronBuilder::BuildHighSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
91{
92 if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
93 return Meson(black,white, SpinOne);
94 } else {
95 return Barion(black,white,SpinThreeHalf);
96 }
97}
98
99//-------------------------------------------------------------------------
100
101G4ParticleDefinition * G4HadronBuilder::Meson(G4ParticleDefinition * black,
102 G4ParticleDefinition * white, Spin theSpin)
103{
104 #ifdef debug_Hbuilder
105 // Verify Input Charge
106 G4double charge = black->GetPDGCharge() + white->GetPDGCharge();
107 if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent ) // 1.001 to avoid int(.9999) -> 0
108 {
109 G4cerr << " G4HadronBuilder::Build()" << G4endl;
110 G4cerr << " Invalid total charge found for on input: "
111 << charge<< G4endl;
112 G4cerr << " PGDcode input quark1/quark2 : " <<
113 black->GetPDGEncoding() << " / "<<
114 white->GetPDGEncoding() << G4endl;
115 G4cerr << G4endl;
116 }
117 #endif
118
119 G4int id1 = black->GetPDGEncoding();
120 G4int id2 = white->GetPDGEncoding();
121
122 // G4int ifl1= std::max(std::abs(id1), std::abs(id2));
123 if ( std::abs(id1) < std::abs(id2) )
124 {
125 G4int xchg = id1;
126 id1 = id2;
127 id2 = xchg;
128 }
129
130 G4int abs_id1 = std::abs(id1);
131
132 if ( abs_id1 > 5 )
133 throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Meson : Illegal Quark content as input");
134
135 G4int PDGEncoding=0;
136
137 if (id1 + id2 == 0) {
138 if ( abs_id1 < 4) { // light quarks: u, d or s
139 G4double rmix = G4UniformRand();
140 G4int imix = 2*std::abs(id1) - 1;
141 if (theSpin == SpinZero) {
142 PDGEncoding = 110*(1 + (G4int)(rmix + scalarMesonMixings[imix - 1])
143 + (G4int)(rmix + scalarMesonMixings[imix])
144 ) + theSpin;
145 } else {
146 PDGEncoding = 110*(1 + (G4int)(rmix + vectorMesonMixings[imix - 1])
147 + (G4int)(rmix + vectorMesonMixings[imix])
148 ) + theSpin;
149 }
150 } else { // for c and b quarks
151
152 PDGEncoding = abs_id1*100 + abs_id1*10;
153
154 if (PDGEncoding == 440) {
155 if ( G4UniformRand() < ProbEta_c ) {
156 PDGEncoding +=1;
157 } else {
158 PDGEncoding +=3;
159 }
160 }
161 if (PDGEncoding == 550) {
162 if ( G4UniformRand() < ProbEta_b ) {
163 PDGEncoding +=1;
164 } else {
165 PDGEncoding +=3;
166 }
167 }
168 }
169 } else {
170 PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) + theSpin;
171 G4bool IsUp = (std::abs(id1)&1) == 0; // quark 1 up type quark (u or c)
172 G4bool IsAnti = id1 < 0; // quark 1 is antiquark?
173 if ( (IsUp && IsAnti ) || (!IsUp && !IsAnti ) ) PDGEncoding = - PDGEncoding;
174 }
175
176 // ---------------------------------------------------------------------
177 // Special treatment for charmed and bottom mesons : in Geant4 there are
178 // no excited charmed or bottom mesons, therefore we need to transform these
179 // into existing charmed and bottom mesons in Geant4. Whenever possible,
180 // we use the corresponding ground state mesons with the same quantum numbers;
181 // else, we prefer to conserve the electric charge rather than other flavor numbers.
182 #ifdef debug_heavyHadrons
183 G4int initialPDGEncoding = PDGEncoding;
184 #endif
185 if ( std::abs( PDGEncoding ) == 10411 ) // D*0(2400)+ -> D+
186 ( PDGEncoding > 0 ? PDGEncoding = 411 : PDGEncoding = -411 );
187 else if ( std::abs( PDGEncoding ) == 10421 ) // D*0(2400)0 -> D0
188 ( PDGEncoding > 0 ? PDGEncoding = 421 : PDGEncoding = -421 );
189 else if ( std::abs( PDGEncoding ) == 413 ) // D*(2010)+ -> D+
190 ( PDGEncoding > 0 ? PDGEncoding = 411 : PDGEncoding = -411 );
191 else if ( std::abs( PDGEncoding ) == 423 ) // D*(2007)0 -> D0
192 ( PDGEncoding > 0 ? PDGEncoding = 421 : PDGEncoding = -421 );
193 else if ( std::abs( PDGEncoding ) == 10413 ) // D1(2420)+ -> D+
194 ( PDGEncoding > 0 ? PDGEncoding = 411 : PDGEncoding = -411 );
195 else if ( std::abs( PDGEncoding ) == 10423 ) // D1(2420)0 -> D0
196 ( PDGEncoding > 0 ? PDGEncoding = 421 : PDGEncoding = -421 );
197 else if ( std::abs( PDGEncoding ) == 20413 ) // D1(H)+ -> D+
198 ( PDGEncoding > 0 ? PDGEncoding = 411 : PDGEncoding = -411 );
199 else if ( std::abs( PDGEncoding ) == 20423 ) // D1(2430)0 -> D0
200 ( PDGEncoding > 0 ? PDGEncoding = 421 : PDGEncoding = -421 );
201 else if ( std::abs( PDGEncoding ) == 415 ) // D2*(2460)+ -> D+
202 ( PDGEncoding > 0 ? PDGEncoding = 411 : PDGEncoding = -411 );
203 else if ( std::abs( PDGEncoding ) == 425 ) // D2*(2460)0 -> D0
204 ( PDGEncoding > 0 ? PDGEncoding = 421 : PDGEncoding = -421 );
205 else if ( std::abs( PDGEncoding ) == 10431 ) // Ds0*(2317)+ -> Ds+
206 ( PDGEncoding > 0 ? PDGEncoding = 431 : PDGEncoding = -431 );
207 else if ( std::abs( PDGEncoding ) == 433 ) // Ds*+ -> Ds+
208 ( PDGEncoding > 0 ? PDGEncoding = 431 : PDGEncoding = -431 );
209 else if ( std::abs( PDGEncoding ) == 10433 ) // Ds1(2536)+ -> Ds+
210 ( PDGEncoding > 0 ? PDGEncoding = 431 : PDGEncoding = -431 );
211 else if ( std::abs( PDGEncoding ) == 20433 ) // Ds1(2460)+ -> Ds+
212 ( PDGEncoding > 0 ? PDGEncoding = 431 : PDGEncoding = -431 );
213 else if ( std::abs( PDGEncoding ) == 435 ) // Ds2*(2573)+ -> Ds+
214 ( PDGEncoding > 0 ? PDGEncoding = 431 : PDGEncoding = -431 );
215 else if ( std::abs( PDGEncoding ) == 10441 ) PDGEncoding = 441; // chi_c0(1P) -> eta_c
216 else if ( std::abs( PDGEncoding ) == 100441 ) PDGEncoding = 441; // eta_c(2S) -> eta_c
217 else if ( std::abs( PDGEncoding ) == 10443 ) PDGEncoding = 443; // h_c(1P) -> J/psi
218 else if ( std::abs( PDGEncoding ) == 20443 ) PDGEncoding = 443; // chi_c1(1P) -> J/psi
219 else if ( std::abs( PDGEncoding ) == 100443 ) PDGEncoding = 443; // psi(2S) -> J/psi
220 else if ( std::abs( PDGEncoding ) == 30443 ) PDGEncoding = 443; // psi(3770) -> J/psi
221 else if ( std::abs( PDGEncoding ) == 9000443 ) PDGEncoding = 443; // psi(4040) -> J/psi
222 else if ( std::abs( PDGEncoding ) == 9010443 ) PDGEncoding = 443; // psi(4160) -> J/psi
223 else if ( std::abs( PDGEncoding ) == 9020443 ) PDGEncoding = 443; // psi(4415) -> J/psi
224 else if ( std::abs( PDGEncoding ) == 445 ) PDGEncoding = 443; // chi_c2(1P) -> J/psi
225 else if ( std::abs( PDGEncoding ) == 100445 ) PDGEncoding = 443; // chi_c2(2P) -> J/psi
226 // Bottom mesons
227 else if ( std::abs( PDGEncoding ) == 10511 ) // B0*0 -> B0
228 ( PDGEncoding > 0 ? PDGEncoding = 511 : PDGEncoding = -511 );
229 else if ( std::abs( PDGEncoding ) == 10521 ) // B0*+ -> B+
230 ( PDGEncoding > 0 ? PDGEncoding = 521 : PDGEncoding = -521 );
231 else if ( std::abs( PDGEncoding ) == 513 ) // B*0 -> B0
232 ( PDGEncoding > 0 ? PDGEncoding = 511 : PDGEncoding = -511 );
233 else if ( std::abs( PDGEncoding ) == 523 ) // B*+ -> B+
234 ( PDGEncoding > 0 ? PDGEncoding = 521 : PDGEncoding = -521 );
235 else if ( std::abs( PDGEncoding ) == 10513 ) // B1(L)0 -> B0
236 ( PDGEncoding > 0 ? PDGEncoding = 511 : PDGEncoding = -511 );
237 else if ( std::abs( PDGEncoding ) == 10523 ) // B1(L)+ -> B+
238 ( PDGEncoding > 0 ? PDGEncoding = 521 : PDGEncoding = -521 );
239 else if ( std::abs( PDGEncoding ) == 20513 ) // B1(H)0 -> B0
240 ( PDGEncoding > 0 ? PDGEncoding = 511 : PDGEncoding = -511 );
241 else if ( std::abs( PDGEncoding ) == 20523 ) // B1(H)+ -> B+
242 ( PDGEncoding > 0 ? PDGEncoding = 521 : PDGEncoding = -521 );
243 else if ( std::abs( PDGEncoding ) == 515 ) // B2*0 -> B0
244 ( PDGEncoding > 0 ? PDGEncoding = 511 : PDGEncoding = -511 );
245 else if ( std::abs( PDGEncoding ) == 525 ) // B2*+ -> B+
246 ( PDGEncoding > 0 ? PDGEncoding = 521 : PDGEncoding = -521 );
247 else if ( std::abs( PDGEncoding ) == 10531 ) // Bs0*0 -> Bs0
248 ( PDGEncoding > 0 ? PDGEncoding = 531 : PDGEncoding = -531 );
249 else if ( std::abs( PDGEncoding ) == 533 ) // Bs*0 -> Bs0
250 ( PDGEncoding > 0 ? PDGEncoding = 531 : PDGEncoding = -531 );
251 else if ( std::abs( PDGEncoding ) == 10533 ) // Bs1(L)0 -> Bs0
252 ( PDGEncoding > 0 ? PDGEncoding = 531 : PDGEncoding = -531 );
253 else if ( std::abs( PDGEncoding ) == 20533 ) // Bs1(H)0 -> Bs0
254 ( PDGEncoding > 0 ? PDGEncoding = 531 : PDGEncoding = -531 );
255 else if ( std::abs( PDGEncoding ) == 535 ) // Bs2*0 -> Bs0
256 ( PDGEncoding > 0 ? PDGEncoding = 531 : PDGEncoding = -531 );
257 else if ( std::abs( PDGEncoding ) == 10541 ) // Bc0*+ -> Bc+
258 ( PDGEncoding > 0 ? PDGEncoding = 541 : PDGEncoding = -541 );
259 else if ( std::abs( PDGEncoding ) == 543 ) // Bc*+ -> Bc+
260 ( PDGEncoding > 0 ? PDGEncoding = 541 : PDGEncoding = -541 );
261 else if ( std::abs( PDGEncoding ) == 10543 ) // Bc1(L)+ -> Bc+
262 ( PDGEncoding > 0 ? PDGEncoding = 541 : PDGEncoding = -541 );
263 else if ( std::abs( PDGEncoding ) == 20543 ) // Bc1(H)+ -> Bc+
264 ( PDGEncoding > 0 ? PDGEncoding = 541 : PDGEncoding = -541 );
265 else if ( std::abs( PDGEncoding ) == 545 ) // Bc2*+ -> Bc+
266 ( PDGEncoding > 0 ? PDGEncoding = 541 : PDGEncoding = -541 );
267 else if ( std::abs( PDGEncoding ) == 551 ) PDGEncoding = 553; // eta_b(1S) -> Upsilon
268 else if ( std::abs( PDGEncoding ) == 10551 ) PDGEncoding = 553; // chi_b0(1P) -> Upsilon
269 else if ( std::abs( PDGEncoding ) == 100551 ) PDGEncoding = 553; // eta_b(2S) -> Upsilon
270 else if ( std::abs( PDGEncoding ) == 110551 ) PDGEncoding = 553; // chi_b0(2P) -> Upsilon
271 else if ( std::abs( PDGEncoding ) == 200551 ) PDGEncoding = 553; // eta_b(3S) -> Upsilon
272 else if ( std::abs( PDGEncoding ) == 210551 ) PDGEncoding = 553; // chi_b0(3P) -> Upsilon
273 else if ( std::abs( PDGEncoding ) == 10553 ) PDGEncoding = 553; // h_b(1P) -> Upsilon
274 else if ( std::abs( PDGEncoding ) == 20553 ) PDGEncoding = 553; // chi_b1(1P) -> Upsilon
275 else if ( std::abs( PDGEncoding ) == 30553 ) PDGEncoding = 553; // Upsilon_1(1D) -> Upsilon
276 else if ( std::abs( PDGEncoding ) == 100553 ) PDGEncoding = 553; // Upsilon(2S) -> Upsilon
277 else if ( std::abs( PDGEncoding ) == 110553 ) PDGEncoding = 553; // h_b(2P) -> Upsilon
278 else if ( std::abs( PDGEncoding ) == 120553 ) PDGEncoding = 553; // chi_b1(2P) -> Upsilon
279 else if ( std::abs( PDGEncoding ) == 130553 ) PDGEncoding = 553; // Upsilon_1(2D) -> Upsilon
280 else if ( std::abs( PDGEncoding ) == 200553 ) PDGEncoding = 553; // Upsilon(3S) -> Upsilon
281 else if ( std::abs( PDGEncoding ) == 210553 ) PDGEncoding = 553; // h_b(3P) -> Upsilon
282 else if ( std::abs( PDGEncoding ) == 220553 ) PDGEncoding = 553; // chi_b1(3P) -> Upsilon
283 else if ( std::abs( PDGEncoding ) == 300553 ) PDGEncoding = 553; // Upsilon(4S) -> Upsilon
284 else if ( std::abs( PDGEncoding ) == 9000553 ) PDGEncoding = 553; // Upsilon(10860) -> Upsilon
285 else if ( std::abs( PDGEncoding ) == 9010553 ) PDGEncoding = 553; // Upsilon(11020) -> Upsilon
286 else if ( std::abs( PDGEncoding ) == 555 ) PDGEncoding = 553; // chi_b2(1P) -> Upsilon
287 else if ( std::abs( PDGEncoding ) == 10555 ) PDGEncoding = 553; // eta_b2(1D) -> Upsilon
288 else if ( std::abs( PDGEncoding ) == 20555 ) PDGEncoding = 553; // Upsilon_2(1D) -> Upsilon
289 else if ( std::abs( PDGEncoding ) == 100555 ) PDGEncoding = 553; // chi_b2(2P) -> Upsilon
290 else if ( std::abs( PDGEncoding ) == 110555 ) PDGEncoding = 553; // eta_b2(2D) -> Upsilon
291 else if ( std::abs( PDGEncoding ) == 120555 ) PDGEncoding = 553; // Upsilon_2(2D) -> Upsilon
292 else if ( std::abs( PDGEncoding ) == 200555 ) PDGEncoding = 553; // chi_b2(3P) -> Upsilon
293 else if ( std::abs( PDGEncoding ) == 557 ) PDGEncoding = 553; // Upsilon_3(1D) -> Upsilon
294 else if ( std::abs( PDGEncoding ) == 100557 ) PDGEncoding = 553; // Upsilon_3(2D) -> Upsilon
295 #ifdef debug_heavyHadrons
296 if ( initialPDGEncoding != PDGEncoding ) {
297 G4cout << "G4HadronBuilder::Meson : forcing (inexisting in G4) heavy meson with pdgCode="
298 << initialPDGEncoding << " into pdgCode=" << PDGEncoding << G4endl;
299 }
300 #endif
301 // ---------------------------------------------------------------------
302
303 G4ParticleDefinition * MesonDef=
304 G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
305
306 #ifdef debug_Hbuilder
307 if (MesonDef == 0 ) {
308 G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
309 << PDGEncoding << G4endl;
310 } else if ( ( black->GetPDGCharge() + white->GetPDGCharge()
311 - MesonDef->GetPDGCharge() ) > perCent ) {
312 G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
313 << " Quark1/2 = "
314 << black->GetParticleName() << " / "
315 << white->GetParticleName()
316 << " resulting Hadron " << MesonDef->GetParticleName()
317 << G4endl;
318 }
319 #endif
320
321 return MesonDef;
322}
323
324//-------------------------------------------------------------------------
325
326G4ParticleDefinition * G4HadronBuilder::Barion(G4ParticleDefinition * black,
327 G4ParticleDefinition * white,Spin theSpin)
328{
329 #ifdef debug_Hbuilder
330 // Verify Input Charge
331 G4double charge = black->GetPDGCharge() + white->GetPDGCharge();
332 if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent )
333 {
334 G4cerr << " G4HadronBuilder::Build()" << G4endl;
335 G4cerr << " Invalid total charge found for on input: "
336 << charge<< G4endl;
337 G4cerr << " PGDcode input quark1/quark2 : " <<
338 black->GetPDGEncoding() << " / "<<
339 white->GetPDGEncoding() << G4endl;
340 G4cerr << G4endl;
341 }
342 #endif
343
344 G4int id1 = black->GetPDGEncoding();
345 G4int id2 = white->GetPDGEncoding();
346
347 if ( std::abs(id1) < std::abs(id2) )
348 {
349 G4int xchg = id1;
350 id1 = id2;
351 id2 = xchg;
352 }
353
354 if (std::abs(id1) < 1000 || std::abs(id2) > 5 )
355 throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Barion: Illegal quark content as input");
356
357 G4int ifl1= std::abs(id1)/1000;
358 G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100;
359 G4int diquarkSpin = std::abs(id1)%10;
360 G4int ifl3 = id2;
361 if (id1 < 0)
362 {
363 ifl1 = - ifl1;
364 ifl2 = - ifl2;
365 }
366 //... Construct barion, distinguish Lambda and Sigma barions.
367 G4int kfla = std::abs(ifl1);
368 G4int kflb = std::abs(ifl2);
369 G4int kflc = std::abs(ifl3);
370
371 G4int kfld = std::max(kfla,kflb);
372 kfld = std::max(kfld,kflc);
373 G4int kflf = std::min(kfla,kflb);
374 kflf = std::min(kflf,kflc);
375
376 G4int kfle = kfla + kflb + kflc - kfld - kflf;
377
378 //... barion with content uuu or ddd or sss has always spin = 3/2
379 theSpin = (kfla == kflb && kflb == kflc)? SpinThreeHalf : theSpin;
380
381 G4int kfll = 0;
382 if (kfld < 6) {
383 if (theSpin == SpinHalf && kfld > kfle && kfle > kflf) {
384 // Spin J=1/2 and all three quarks different
385 // Two states exist: (uds -> lambda or sigma0)
386 // - lambda: s(ud)0 s : 3122; ie. reverse the two lighter quarks
387 // - sigma0: s(ud)1 s : 3212
388 if (diquarkSpin == 1 ) {
389 if ( kfla == kfld) { // heaviest quark in diquark
390 kfll = 1;
391 } else {
392 kfll = (G4int)(0.25 + G4UniformRand());
393 }
394 }
395 if (diquarkSpin == 3 && kfla != kfld)
396 kfll = (G4int)(0.75 + G4UniformRand());
397 }
398 }
399
400 G4int PDGEncoding;
401 if (kfll == 1)
402 PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin;
403 else
404 PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin;
405
406 if (id1 < 0)
407 PDGEncoding = -PDGEncoding;
408
409 // ---------------------------------------------------------------------
410 // Special treatment for charmed and bottom baryons : in Geant4 there are
411 // neither excited charmed or bottom baryons, nor baryons with two or three
412 // heavy (c, b) constitutent quarks:
413 // Sigma_c* , Xi_c' , Xi_c* , Omega_c* ,
414 // Xi_cc , Xi_cc* , Omega_cc , Omega_cc* , Omega_ccc ;
415 // Sigma_b* , Xi_b' , Xi_b* , Omega_b*,
416 // Xi_bc , Xi_bc' , Xi_bc* , Omega_bc , Omega_bc' , Omega_bc* ,
417 // Omega_bcc , Omega_bcc* , Xi_bb, Xi_bb* , Omega_bb, Omega_bb* ,
418 // Omega_bbc , Omega_bbc* , Omega_bbb
419 // therefore we need to transform these into existing charmed and bottom
420 // baryons in Geant4. Whenever possible, we use the corresponding ground state
421 // baryons with the same quantum numbers; else, we prefer to conserve the
422 // electric charge rather than other flavor numbers.
423 #ifdef debug_heavyHadrons
424 G4int charmViolation = 0, bottomViolation = 0; // Only positive
425 G4int initialPDGEncoding = PDGEncoding;
426 #endif
427 if ( std::abs( PDGEncoding ) == 4224 ) { // Sigma_c*++ -> Sigma_c++
428 ( PDGEncoding > 0 ? PDGEncoding = 4222 : PDGEncoding = -4222 );
429 } else if ( std::abs( PDGEncoding ) == 4214 ) { // Sigma_c*+ -> Sigma_c+
430 ( PDGEncoding > 0 ? PDGEncoding = 4212 : PDGEncoding = -4212 );
431 } else if ( std::abs( PDGEncoding ) == 4114 ) { // Sigma_c*0 -> Sigma_c0
432 ( PDGEncoding > 0 ? PDGEncoding = 4112 : PDGEncoding = -4112 );
433 } else if ( std::abs( PDGEncoding ) == 4322 ) { // Xi_c'+ -> Xi_c+
434 ( PDGEncoding > 0 ? PDGEncoding = 4232 : PDGEncoding = -4232 );
435 } else if ( std::abs( PDGEncoding ) == 4312 ) { // Xi_c'0 -> Xi_c0
436 ( PDGEncoding > 0 ? PDGEncoding = 4132 : PDGEncoding = -4132 );
437 } else if ( std::abs( PDGEncoding ) == 4324 ) { // Xi_c*+ -> Xi_c+
438 ( PDGEncoding > 0 ? PDGEncoding = 4232 : PDGEncoding = -4232 );
439 } else if ( std::abs( PDGEncoding ) == 4314 ) { // Xi_c*0 -> Xi_c0
440 ( PDGEncoding > 0 ? PDGEncoding = 4132 : PDGEncoding = -4132 );
441 } else if ( std::abs( PDGEncoding ) == 4334 ) { // Omega_c*0 -> Omega_c0
442 ( PDGEncoding > 0 ? PDGEncoding = 4332 : PDGEncoding = -4332 );
443 } else if ( std::abs( PDGEncoding ) == 4412 ) { // Xi_cc+ -> Xi_c+
444 ( PDGEncoding > 0 ? PDGEncoding = 4232 : PDGEncoding = -4232 );
445 #ifdef debug_heavyHadrons
446 charmViolation = 1;
447 #endif
448 } else if ( std::abs( PDGEncoding ) == 4422 ) { // Xi_cc++ -> Sigma_c++ (use Sigma to conserve charge)
449 ( PDGEncoding > 0 ? PDGEncoding = 4222 : PDGEncoding = -4222 );
450 #ifdef debug_heavyHadrons
451 charmViolation = 1;
452 #endif
453 } else if ( std::abs( PDGEncoding ) == 4414 ) { // Xi_cc*+ -> Xi_c+
454 ( PDGEncoding > 0 ? PDGEncoding = 4232 : PDGEncoding = -4232 );
455 #ifdef debug_heavyHadrons
456 charmViolation = 1;
457 #endif
458 } else if ( std::abs( PDGEncoding ) == 4424 ) { // Xi_cc*++ -> Sigma_c++ (use Sigma to conserve charge)
459 ( PDGEncoding > 0 ? PDGEncoding = 4222 : PDGEncoding = -4222 );
460 #ifdef debug_heavyHadrons
461 charmViolation = 1;
462 #endif
463 } else if ( std::abs( PDGEncoding ) == 4432 ) { // Omega_cc+ -> Xi_c+ (use Xi to conserve charge)
464 ( PDGEncoding > 0 ? PDGEncoding = 4232 : PDGEncoding = -4232 );
465 #ifdef debug_heavyHadrons
466 charmViolation = 1;
467 #endif
468 } else if ( std::abs( PDGEncoding ) == 4434 ) { // Omega_cc*+ -> Xi_c+ (use Xi to conserve charge)
469 ( PDGEncoding > 0 ? PDGEncoding = 4232 : PDGEncoding = -4232 );
470 #ifdef debug_heavyHadrons
471 charmViolation = 1;
472 #endif
473 } else if ( std::abs( PDGEncoding ) == 4444 ) { // Omega_ccc++ -> Sigma_c++ (use Sigma to conserve charge)
474 ( PDGEncoding > 0 ? PDGEncoding = 4222 : PDGEncoding = -4222 );
475 #ifdef debug_heavyHadrons
476 charmViolation = 2;
477 #endif
478 // Bottom baryons
479 } else if ( std::abs( PDGEncoding ) == 5114 ) { // Sigma_b*- -> Sigma_b-
480 ( PDGEncoding > 0 ? PDGEncoding = 5112 : PDGEncoding = -5112 );
481 } else if ( std::abs( PDGEncoding ) == 5214 ) { // Sigma_b*0 -> Sigma_b0
482 ( PDGEncoding > 0 ? PDGEncoding = 5212 : PDGEncoding = -5212 );
483 } else if ( std::abs( PDGEncoding ) == 5224 ) { // Sigma_b*+ -> Sigma_b+
484 ( PDGEncoding > 0 ? PDGEncoding = 5222 : PDGEncoding = -5222 );
485 } else if ( std::abs( PDGEncoding ) == 5312 ) { // Xi_b'- -> Xi_b-
486 ( PDGEncoding > 0 ? PDGEncoding = 5132 : PDGEncoding = -5132 );
487 } else if ( std::abs( PDGEncoding ) == 5322 ) { // Xi_b'0 -> Xi_b0
488 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
489 } else if ( std::abs( PDGEncoding ) == 5314 ) { // Xi_b*- -> Xi_b-
490 ( PDGEncoding > 0 ? PDGEncoding = 5132 : PDGEncoding = -5132 );
491 } else if ( std::abs( PDGEncoding ) == 5324 ) { // Xi_b*0 -> Xi_b0
492 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
493 } else if ( std::abs( PDGEncoding ) == 5334 ) { // Omega_b*- -> Omega_b-
494 ( PDGEncoding > 0 ? PDGEncoding = 5332 : PDGEncoding = -5332 );
495 } else if ( std::abs( PDGEncoding ) == 5142 ) { // Xi_bc0 -> Xi_b0
496 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
497 #ifdef debug_heavyHadrons
498 charmViolation = 1;
499 #endif
500 } else if ( std::abs( PDGEncoding ) == 5242 ) { // Xi_bc+ -> Sigma_b+ (use Sigma to conserve charge)
501 ( PDGEncoding > 0 ? PDGEncoding = 5222 : PDGEncoding = -5222 );
502 #ifdef debug_heavyHadrons
503 charmViolation = 1;
504 #endif
505 } else if ( std::abs( PDGEncoding ) == 5412 ) { // Xi_bc'0 -> Xi_b0
506 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
507 #ifdef debug_heavyHadrons
508 charmViolation = 1;
509 #endif
510 } else if ( std::abs( PDGEncoding ) == 5422 ) { // Xi_bc'+ -> Sigma_b+ (use Sigma to conserve charge)
511 ( PDGEncoding > 0 ? PDGEncoding = 5222 : PDGEncoding = -5222 );
512 #ifdef debug_heavyHadrons
513 charmViolation = 1;
514 #endif
515 } else if ( std::abs( PDGEncoding ) == 5414 ) { // Xi_bc*0 -> Xi_b0
516 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
517 #ifdef debug_heavyHadrons
518 charmViolation = 1;
519 #endif
520 } else if ( std::abs( PDGEncoding ) == 5424 ) { // Xi_bc*+ -> Sigma_b+ (use Sigma to conserve charge)
521 ( PDGEncoding > 0 ? PDGEncoding = 5222 : PDGEncoding = -5222 );
522 #ifdef debug_heavyHadrons
523 charmViolation = 1;
524 #endif
525 } else if ( std::abs( PDGEncoding ) == 5342 ) { // Omega_bc0 -> Xi_b0 (use Xi to conserve charge)
526 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
527 #ifdef debug_heavyHadrons
528 charmViolation = 1;
529 #endif
530 } else if ( std::abs( PDGEncoding ) == 5432 ) { // Omega_bc'0 -> Xi_b0 (use Xi to conserve charge)
531 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
532 #ifdef debug_heavyHadrons
533 charmViolation = 1;
534 #endif
535 } else if ( std::abs( PDGEncoding ) == 5434 ) { // Omega_bc*0 -> Xi_b0 (use Xi to conserve charge)
536 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
537 #ifdef debug_heavyHadrons
538 charmViolation = 1;
539 #endif
540 } else if ( std::abs( PDGEncoding ) == 5442 ) { // Omega_bcc+ -> Sigma_b+ (use Sigma to conserve charge)
541 ( PDGEncoding > 0 ? PDGEncoding = 5222 : PDGEncoding = -5222 );
542 #ifdef debug_heavyHadrons
543 charmViolation = 2;
544 #endif
545 } else if ( std::abs( PDGEncoding ) == 5444 ) { // Omega_bcc*+ -> Sigma_b+ (use Sigma to conserve charge)
546 ( PDGEncoding > 0 ? PDGEncoding = 5222 : PDGEncoding = -5222 );
547 #ifdef debug_heavyHadrons
548 charmViolation = 2;
549 #endif
550 } else if ( std::abs( PDGEncoding ) == 5512 ) { // Xi_bb- -> Xi_b-
551 ( PDGEncoding > 0 ? PDGEncoding = 5132 : PDGEncoding = -5132 );
552 #ifdef debug_heavyHadrons
553 bottomViolation = 1;
554 #endif
555 } else if ( std::abs( PDGEncoding ) == 5522 ) { // Xi_bb0 -> Xi_b0
556 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
557 #ifdef debug_heavyHadrons
558 bottomViolation = 1;
559 #endif
560 } else if ( std::abs( PDGEncoding ) == 5514 ) { // Xi_bb*- -> Xi_b-
561 ( PDGEncoding > 0 ? PDGEncoding = 5132 : PDGEncoding = -5132 );
562 #ifdef debug_heavyHadrons
563 bottomViolation = 1;
564 #endif
565 } else if ( std::abs( PDGEncoding ) == 5524 ) { // Xi_bb*0 -> Xi_b0
566 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
567 #ifdef debug_heavyHadrons
568 bottomViolation = 1;
569 #endif
570 } else if ( std::abs( PDGEncoding ) == 5532 ) { // Omega_bb- -> Omega_b-
571 ( PDGEncoding > 0 ? PDGEncoding = 5332 : PDGEncoding = -5332 );
572 #ifdef debug_heavyHadrons
573 bottomViolation = 1;
574 #endif
575 } else if ( std::abs( PDGEncoding ) == 5534 ) { // Omega_bb*- -> Omega_b-
576 ( PDGEncoding > 0 ? PDGEncoding = 5332 : PDGEncoding = -5332 );
577 #ifdef debug_heavyHadrons
578 bottomViolation = 1;
579 #endif
580 } else if ( std::abs( PDGEncoding ) == 5542 ) { // Omega_bbc0 -> Xi_b0 (use Xi to conserve charge)
581 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
582 #ifdef debug_heavyHadrons
583 charmViolation = 1; bottomViolation = 1;
584 #endif
585 } else if ( std::abs( PDGEncoding ) == 5544 ) { // Omega_bbc*0 -> Xi_b0 (use Xi to conserve charge)
586 ( PDGEncoding > 0 ? PDGEncoding = 5232 : PDGEncoding = -5232 );
587 #ifdef debug_heavyHadrons
588 charmViolation = 1; bottomViolation = 1;
589 #endif
590 } else if ( std::abs( PDGEncoding ) == 5554 ) { // Omega_bbb- -> Omega_b-
591 ( PDGEncoding > 0 ? PDGEncoding = 5332 : PDGEncoding = -5332 );
592 #ifdef debug_heavyHadrons
593 bottomViolation = 2;
594 #endif
595 }
596 #ifdef debug_heavyHadrons
597 if ( initialPDGEncoding != PDGEncoding ) {
598 G4cout << "G4HadronBuilder::Barion : forcing (inexisting in G4) heavy baryon with pdgCode="
599 << initialPDGEncoding << " into pdgCode=" << PDGEncoding << G4endl;
600 if ( charmViolation != 0 || bottomViolation != 0 ) {
601 G4cout << "\t --> VIOLATION of " << ( charmViolation != 0 ? " CHARM " : " " )
602 << ( charmViolation != 0 && bottomViolation != 0 ? " and " : " " )
603 << ( bottomViolation != 0 ? " BOTTOM " : " " ) << " quantum number ! " << G4endl;
604 }
605 }
606 #endif
607 // ---------------------------------------------------------------------
608
609 G4ParticleDefinition * BarionDef=
610 G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
611
612 #ifdef debug_Hbuilder
613 if (BarionDef == 0 ) {
614 G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
615 << PDGEncoding << G4endl;
616 } else if ( ( black->GetPDGCharge() + white->GetPDGCharge()
617 - BarionDef->GetPDGCharge() ) > perCent ) {
618 G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
619 << " DiQuark/Quark = "
620 << black->GetParticleName() << " / "
621 << white->GetParticleName()
622 << " resulting Hadron " << BarionDef->GetParticleName()
623 << G4endl;
624 }
625 #endif
626
627 return BarionDef;
628}
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4cerr
G4GLOB_DLL std::ostream G4cerr
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4HadronBuilder::Build
G4ParticleDefinition * Build(G4ParticleDefinition *black, G4ParticleDefinition *white)
Definition G4HadronBuilder.cc:60
G4HadronBuilder::BuildLowSpin
G4ParticleDefinition * BuildLowSpin(G4ParticleDefinition *black, G4ParticleDefinition *white)
Definition G4HadronBuilder.cc:78
G4HadronBuilder::G4HadronBuilder
G4HadronBuilder(const std::vector< G4double > &mesonMix, const G4double barionMix, const std::vector< G4double > &scalarMesonMix, const std::vector< G4double > &vectorMesonMix, const G4double Eta_cProb, const G4double Eta_bProb)
Definition G4HadronBuilder.cc:45
G4HadronBuilder::BuildHighSpin
G4ParticleDefinition * BuildHighSpin(G4ParticleDefinition *black, G4ParticleDefinition *white)
Definition G4HadronBuilder.cc:90
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4ParticleDefinition::GetParticleSubType
const G4String & GetParticleSubType() const
Definition G4ParticleDefinition.hh:119
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition G4ParticleTable.cc:498
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition G4ParticleTable.cc:82