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