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 | // $Id: G4QHadronBuilder.cc,v 1.4 2009/04/29 07:53:18 mkossov Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
29 | // |
---|
30 | // ----------------------------------------------------------------------------- |
---|
31 | // GEANT4 class implementation file |
---|
32 | // |
---|
33 | // Created by Mikhail Kosov, October 2006 |
---|
34 | // Simple algorithm making a hadron out of two partons |
---|
35 | // For comparison mirror member functions are taken from G4 class: |
---|
36 | // G4HadronBuilder |
---|
37 | // ----------------------------------------------------------------------------- |
---|
38 | // Short description: A CHIPS class for the builder of the G4QHadron, which is a |
---|
39 | // resulting object for the string fragmentation. the G4QHadron has specific |
---|
40 | // parameters, which are not included in the G4QParticle from the CHIPS World, |
---|
41 | // but necessary for the string fragmentation. When the G4QHadron is created |
---|
42 | // (builded), it is converted to the CHIPS particle. |
---|
43 | // ----------------------------------------------------------------------------- |
---|
44 | |
---|
45 | |
---|
46 | //#define debug |
---|
47 | |
---|
48 | #include "G4QHadronBuilder.hh" |
---|
49 | #include "Randomize.hh" |
---|
50 | #include "G4ParticleTable.hh" |
---|
51 | |
---|
52 | G4QHadronBuilder::G4QHadronBuilder() |
---|
53 | { |
---|
54 | mesonSpinMix = 0.5; // probability to create vector meson |
---|
55 | baryonSpinMix= 0.5; // probability to create 3/2 baryon |
---|
56 | scalarMesonMixings.resize(6); |
---|
57 | scalarMesonMixings[0] = 0.5; |
---|
58 | scalarMesonMixings[1] = 0.25; |
---|
59 | scalarMesonMixings[2] = 0.5; |
---|
60 | scalarMesonMixings[3] = 0.25; |
---|
61 | scalarMesonMixings[4] = 1.0; |
---|
62 | scalarMesonMixings[5] = 0.5; |
---|
63 | vectorMesonMixings.resize(6); |
---|
64 | vectorMesonMixings[0] = 0.5; |
---|
65 | vectorMesonMixings[1] = 0.0; |
---|
66 | vectorMesonMixings[2] = 0.5; |
---|
67 | vectorMesonMixings[3] = 0.0; |
---|
68 | vectorMesonMixings[4] = 1.0; |
---|
69 | vectorMesonMixings[5] = 1.0; |
---|
70 | } |
---|
71 | |
---|
72 | G4QHadron* G4QHadronBuilder::Build(G4QParton* black, G4QParton* white) |
---|
73 | { |
---|
74 | if(black->GetParticleSubType()=="di_quark" || white->GetParticleSubType()=="di_quark") |
---|
75 | { |
---|
76 | // Baryon consists of quark and at least one di-quark |
---|
77 | Spin spin = (G4UniformRand() < baryonSpinMix) ? SpinHalf : SpinThreeHalf; |
---|
78 | return Baryon(black,white,spin); |
---|
79 | } |
---|
80 | else |
---|
81 | { |
---|
82 | // Meson consists of quark and abnti-quark |
---|
83 | Spin spin = (G4UniformRand() < mesonSpinMix) ? SpinZero : SpinOne; |
---|
84 | return Meson(black,white,spin); |
---|
85 | } |
---|
86 | } |
---|
87 | |
---|
88 | G4QHadron* G4QHadronBuilder::BuildLowSpin(G4QParton* black, G4QParton* white) |
---|
89 | { |
---|
90 | if(black->GetParticleSubType() == "quark" && white->GetParticleSubType() == "quark") |
---|
91 | return Meson(black,white, SpinZero); |
---|
92 | // returns a SpinThreeHalf Baryon if all quarks are the same |
---|
93 | else return Baryon(black,white, SpinHalf); |
---|
94 | } |
---|
95 | |
---|
96 | G4QHadron* G4QHadronBuilder::BuildHighSpin(G4QParton* black, G4QParton* white) |
---|
97 | { |
---|
98 | if(black->GetParticleSubType() == "quark" && white->GetParticleSubType() == "quark") |
---|
99 | return Meson(black,white, SpinOne); |
---|
100 | else return Baryon(black,white,SpinThreeHalf); |
---|
101 | } |
---|
102 | |
---|
103 | G4QHadron* G4QHadronBuilder::Meson(G4QParton* black, G4QParton* white, Spin theSpin) |
---|
104 | { |
---|
105 | #ifdef debug |
---|
106 | // Verify Input Charge |
---|
107 | G4double charge = black->GetPDGCharge() + white->GetPDGCharge(); |
---|
108 | if (std::abs(charge)>2 || std::abs(3*(charge-G4int(charge*1.001))) > perCent) |
---|
109 | G4cerr<<"-Warning-G4QHadronBuilder::Meson: invalid TotalCharge="<<charge<<", PDG_q1=" |
---|
110 | <<black->GetPDGEncoding()<< ", PDG_q2="<<white->GetPDGEncoding()<<G4endl; |
---|
111 | #endif |
---|
112 | |
---|
113 | G4int id1= black->GetPDGCode(); |
---|
114 | G4int id2= white->GetPDGCode(); |
---|
115 | if (std::abs(id1) < std::abs(id2)) // exchange black and white |
---|
116 | { |
---|
117 | G4int xchg = id1; |
---|
118 | id1 = id2; |
---|
119 | id2 = xchg; |
---|
120 | } |
---|
121 | if(std::abs(id1)>3) |
---|
122 | { |
---|
123 | G4cerr<<"***G4QHadronBuilder::Meson: q1="<<id1<<", q2="<<id2 |
---|
124 | <<" while CHIPS is only SU(3)"<<G4endl; |
---|
125 | G4Exception("G4QHadronBuilder::Meson:","72",FatalException,"HeavyQuarkFound"); |
---|
126 | } |
---|
127 | G4int PDGEncoding=0; |
---|
128 | if(!(id1+id2)) // annihilation case (neutral) |
---|
129 | { |
---|
130 | G4double rmix = G4UniformRand(); |
---|
131 | G4int imix = 2*std::abs(id1) - 1; |
---|
132 | if(theSpin == SpinZero) |
---|
133 | PDGEncoding = 110*(1 + G4int(rmix + scalarMesonMixings[imix - 1]) |
---|
134 | + G4int(rmix + scalarMesonMixings[imix] ) ) + theSpin; |
---|
135 | else |
---|
136 | PDGEncoding = 110*(1 + G4int(rmix + vectorMesonMixings[imix - 1]) |
---|
137 | + G4int(rmix + vectorMesonMixings[imix] ) ) + theSpin; |
---|
138 | } |
---|
139 | else |
---|
140 | { |
---|
141 | PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) + theSpin; |
---|
142 | G4bool IsUp = (std::abs(id1)&1) == 0; // quark 1 is up type quark (u or c?) |
---|
143 | G4bool IsAnti = id1 < 0; // quark 1 is an antiquark? |
---|
144 | if( (IsUp && IsAnti) || (!IsUp && !IsAnti) ) PDGEncoding = - PDGEncoding; |
---|
145 | // Correction for the true neutral mesons |
---|
146 | if( PDGEncoding == -111 || PDGEncoding == -113 || PDGEncoding == -223 || |
---|
147 | PDGEncoding == -221 || PDGEncoding == -331 || PDGEncoding == -333 ) |
---|
148 | PDGEncoding = - PDGEncoding; |
---|
149 | } |
---|
150 | G4QHadron* Meson= new G4QHadron(PDGEncoding); |
---|
151 | #ifdef debug |
---|
152 | if(std::abs(black->GetPDGCharge() + white->GetPDGCharge() - Meson->GetCharge()) > .001) |
---|
153 | G4cout<<"-Warning-G4QHadronBuilder::Meson:wrongCharge, q1="<<black->GetPDGCharge()<<"(" |
---|
154 | <<black->->GetParticleName()<<"), q2="<<white->GetPDGCharge()<<"(" |
---|
155 | <<white->->GetParticleName()<<"), qM="<<Meson->GetCharge()<<"/"<<PDGEncoding |
---|
156 | <<G4endl; |
---|
157 | #endif |
---|
158 | return Meson; |
---|
159 | } |
---|
160 | |
---|
161 | G4QHadron* G4QHadronBuilder::Baryon(G4QParton* black, G4QParton* white, Spin theSpin) |
---|
162 | { |
---|
163 | #ifdef debug |
---|
164 | // Verify Input Charge |
---|
165 | G4double charge = black->GetPDGCharge() + white->GetPDGCharge(); |
---|
166 | if(std::abs(charge) > 2 || std::abs(3*(charge - G4int(charge*1.001))) > perCent ) |
---|
167 | G4cerr<<"-Warning-G4QHadronBuilder::Baryon: invalid TotalCharge="<<charge<<", PDG_q1=" |
---|
168 | <<black->GetPDGEncoding()<< ", PDG_q2="<<white->GetPDGEncoding()<<G4endl; |
---|
169 | #endif |
---|
170 | G4int id1= black->GetPDGCode(); |
---|
171 | G4int id2= white->GetPDGCode(); |
---|
172 | if(std::abs(id1) < std::abs(id2)) |
---|
173 | { |
---|
174 | G4int xchg = id1; |
---|
175 | id1 = id2; |
---|
176 | id2 = xchg; |
---|
177 | } |
---|
178 | if(std::abs(id1)<1000 || std::abs(id2)> 3) |
---|
179 | { |
---|
180 | G4cerr<<"***G4QHadronBuilder::Baryon: q1="<<id1<<", q2="<<id2 |
---|
181 | <<" can't create a Baryon"<<G4endl; |
---|
182 | G4Exception("G4QHadronBuilder::Baryon:","72",FatalException,"WrongQdQSequence"); |
---|
183 | } |
---|
184 | G4int ifl1= std::abs(id1)/1000; |
---|
185 | G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100; |
---|
186 | G4int diquarkSpin = std::abs(id1)%10; |
---|
187 | G4int ifl3 = id2; |
---|
188 | if (id1 < 0) {ifl1 = - ifl1; ifl2 = - ifl2;} |
---|
189 | //... Construct baryon, distinguish Lambda and Sigma baryons. |
---|
190 | G4int kfla = std::abs(ifl1); |
---|
191 | G4int kflb = std::abs(ifl2); |
---|
192 | G4int kflc = std::abs(ifl3); |
---|
193 | G4int kfld = std::max(kfla,kflb); |
---|
194 | kfld = std::max(kfld,kflc); |
---|
195 | G4int kflf = std::min(kfla,kflb); |
---|
196 | kflf = std::min(kflf,kflc); |
---|
197 | G4int kfle = kfla + kflb + kflc - kfld - kflf; |
---|
198 | //... baryon with content uuu or ddd or sss has always spin = 3/2 |
---|
199 | if(kfla==kflb && kflb==kflc) theSpin=SpinThreeHalf; |
---|
200 | |
---|
201 | G4int kfll = 0; |
---|
202 | if(theSpin == SpinHalf && kfld > kfle && kfle > kflf) |
---|
203 | { |
---|
204 | // Spin J=1/2 and all three quarks different |
---|
205 | // Two states exist: (uds -> lambda or sigma0) |
---|
206 | // - lambda: s(ud)0 s : 3122; ie. reverse the two lighter quarks |
---|
207 | // - sigma0: s(ud)1 s : 3212 |
---|
208 | if(diquarkSpin == 1 ) |
---|
209 | { |
---|
210 | if ( kfla == kfld) kfll = 1; // heaviest quark in diquark |
---|
211 | else kfll = G4int(0.25 + G4UniformRand()); |
---|
212 | } |
---|
213 | if(diquarkSpin==3 && kfla!=kfld) kfll = G4int(0.75+G4UniformRand()); |
---|
214 | } |
---|
215 | G4int PDGEncoding=0; |
---|
216 | if (kfll == 1) PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin; |
---|
217 | else PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin; |
---|
218 | if (id1 < 0) PDGEncoding = -PDGEncoding; |
---|
219 | G4QHadron* Baryon= new G4QHadron(PDGEncoding); |
---|
220 | #ifdef debug |
---|
221 | if(std::abs(black->GetPDGCharge() + white->GetPDGCharge()- Baryon->GetCharge()) > .001) |
---|
222 | G4cout<<"-Warning-G4QHadronBuilder::Baryon:wrongCharge,dq="<<black->GetPDGCharge()<<"(" |
---|
223 | <<black->->GetParticleName()<<"), q="<<white->GetPDGCharge()<<"(" |
---|
224 | <<white->->GetParticleName()<<"), qB="<<Baryon->GetCharge()<<"/"<<PDGEncoding |
---|
225 | <<G4endl; |
---|
226 | #endif |
---|
227 | return Baryon; |
---|
228 | } |
---|