00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "G4HadronBuilder.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "Randomize.hh"
00040 #include "G4HadronicException.hh"
00041 #include "G4ParticleTable.hh"
00042
00043 G4HadronBuilder::G4HadronBuilder(G4double mesonMix, G4double barionMix,
00044 std::vector<double> scalarMesonMix,
00045 std::vector<double> vectorMesonMix)
00046 {
00047 mesonSpinMix=mesonMix;
00048 barionSpinMix=barionMix;
00049 scalarMesonMixings=scalarMesonMix;
00050 vectorMesonMixings=vectorMesonMix;
00051 }
00052
00053 G4ParticleDefinition * G4HadronBuilder::Build(G4ParticleDefinition * black, G4ParticleDefinition * white)
00054 {
00055
00056 if (black->GetParticleSubType()== "di_quark" || white->GetParticleSubType()== "di_quark" ) {
00057
00058
00059 Spin spin = (G4UniformRand() < barionSpinMix) ? SpinHalf : SpinThreeHalf;
00060 return Barion(black,white,spin);
00061
00062 } else {
00063
00064
00065 Spin spin = (G4UniformRand() < mesonSpinMix) ? SpinZero : SpinOne;
00066 return Meson(black,white,spin);
00067
00068 }
00069 }
00070
00071
00072
00073 G4ParticleDefinition * G4HadronBuilder::BuildLowSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
00074 {
00075 if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
00076 return Meson(black,white, SpinZero);
00077 } else {
00078
00079 return Barion(black,white, SpinHalf);
00080 }
00081 }
00082
00083
00084
00085 G4ParticleDefinition * G4HadronBuilder::BuildHighSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
00086 {
00087 if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
00088 return Meson(black,white, SpinOne);
00089 } else {
00090 return Barion(black,white,SpinThreeHalf);
00091 }
00092 }
00093
00094
00095
00096 G4ParticleDefinition * G4HadronBuilder::Meson(G4ParticleDefinition * black,
00097 G4ParticleDefinition * white, Spin theSpin)
00098 {
00099 #ifdef G4VERBOSE
00100
00101
00102 G4double charge = black->GetPDGCharge()
00103 + white->GetPDGCharge();
00104 if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent )
00105 {
00106 G4cerr << " G4HadronBuilder::Build()" << G4endl;
00107 G4cerr << " Invalid total charge found for on input: "
00108 << charge<< G4endl;
00109 G4cerr << " PGDcode input quark1/quark2 : " <<
00110 black->GetPDGEncoding() << " / "<<
00111 white->GetPDGEncoding() << G4endl;
00112 G4cerr << G4endl;
00113 }
00114 #endif
00115
00116 G4int id1= black->GetPDGEncoding();
00117 G4int id2= white->GetPDGEncoding();
00118
00119 if ( std::abs(id1) < std::abs(id2) )
00120 {
00121 G4int xchg = id1;
00122 id1 = id2;
00123 id2 = xchg;
00124 }
00125
00126 if (std::abs(id1) > 3 )
00127 throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Meson : Illegal Quark content as input");
00128
00129 G4int PDGEncoding=0;
00130
00131 if (id1 + id2 == 0) {
00132 G4double rmix = G4UniformRand();
00133 G4int imix = 2*std::abs(id1) - 1;
00134 if(theSpin == SpinZero) {
00135 PDGEncoding = 110*(1 + (G4int)(rmix + scalarMesonMixings[imix - 1])
00136 + (G4int)(rmix + scalarMesonMixings[imix])
00137 ) + theSpin;
00138 } else {
00139 PDGEncoding = 110*(1 + (G4int)(rmix + vectorMesonMixings[imix - 1])
00140 + (G4int)(rmix + vectorMesonMixings[imix])
00141 ) + theSpin;
00142 }
00143 } else {
00144 PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) + theSpin;
00145 G4bool IsUp = (std::abs(id1)&1) == 0;
00146 G4bool IsAnti = id1 < 0;
00147 if( (IsUp && IsAnti ) || (!IsUp && !IsAnti ) )
00148 PDGEncoding = - PDGEncoding;
00149 }
00150
00151
00152 G4ParticleDefinition * MesonDef=
00153 G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
00154 #ifdef G4VERBOSE
00155 if (MesonDef == 0 ) {
00156 G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
00157 << PDGEncoding << G4endl;
00158 } else if ( ( black->GetPDGCharge() + white->GetPDGCharge()
00159 - MesonDef->GetPDGCharge() ) > perCent ) {
00160 G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
00161 << " Quark1/2 = "
00162 << black->GetParticleName() << " / "
00163 << white->GetParticleName()
00164 << " resulting Hadron " << MesonDef->GetParticleName()
00165 << G4endl;
00166 }
00167 #endif
00168
00169 return MesonDef;
00170 }
00171
00172
00173 G4ParticleDefinition * G4HadronBuilder::Barion(G4ParticleDefinition * black,
00174 G4ParticleDefinition * white,Spin theSpin)
00175 {
00176
00177 #ifdef G4VERBOSE
00178
00179 G4double charge = black->GetPDGCharge()
00180 + white->GetPDGCharge();
00181 if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent )
00182 {
00183 G4cerr << " G4HadronBuilder::Build()" << G4endl;
00184 G4cerr << " Invalid total charge found for on input: "
00185 << charge<< G4endl;
00186 G4cerr << " PGDcode input quark1/quark2 : " <<
00187 black->GetPDGEncoding() << " / "<<
00188 white->GetPDGEncoding() << G4endl;
00189 G4cerr << G4endl;
00190 }
00191 #endif
00192 G4int id1= black->GetPDGEncoding();
00193 G4int id2= white->GetPDGEncoding();
00194 if ( std::abs(id1) < std::abs(id2) )
00195 {
00196 G4int xchg = id1;
00197 id1 = id2;
00198 id2 = xchg;
00199 }
00200
00201 if (std::abs(id1) < 1000 || std::abs(id2) > 3 )
00202 throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Barion: Illegal quark content as input");
00203
00204 G4int ifl1= std::abs(id1)/1000;
00205 G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100;
00206 G4int diquarkSpin = std::abs(id1)%10;
00207 G4int ifl3 = id2;
00208 if (id1 < 0)
00209 {
00210 ifl1 = - ifl1;
00211 ifl2 = - ifl2;
00212 }
00213
00214 G4int kfla = std::abs(ifl1);
00215 G4int kflb = std::abs(ifl2);
00216 G4int kflc = std::abs(ifl3);
00217
00218 G4int kfld = std::max(kfla,kflb);
00219 kfld = std::max(kfld,kflc);
00220 G4int kflf = std::min(kfla,kflb);
00221 kflf = std::min(kflf,kflc);
00222
00223 G4int kfle = kfla + kflb + kflc - kfld - kflf;
00224
00225
00226 theSpin = (kfla == kflb && kflb == kflc)? SpinThreeHalf : theSpin;
00227
00228 G4int kfll = 0;
00229 if(theSpin == SpinHalf && kfld > kfle && kfle > kflf) {
00230
00231
00232
00233
00234 if(diquarkSpin == 1 ) {
00235 if ( kfla == kfld) {
00236 kfll = 1;
00237 } else {
00238 kfll = (G4int)(0.25 + G4UniformRand());
00239 }
00240 }
00241 if(diquarkSpin == 3 && kfla != kfld)
00242 kfll = (G4int)(0.75 + G4UniformRand());
00243 }
00244
00245 G4int PDGEncoding;
00246 if (kfll == 1)
00247 PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin;
00248 else
00249 PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin;
00250
00251 if (id1 < 0)
00252 PDGEncoding = -PDGEncoding;
00253
00254
00255 G4ParticleDefinition * BarionDef=
00256 G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
00257 #ifdef G4VERBOSE
00258 if (BarionDef == 0 ) {
00259 G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
00260 << PDGEncoding << G4endl;
00261 } else if ( ( black->GetPDGCharge() + white->GetPDGCharge()
00262 - BarionDef->GetPDGCharge() ) > perCent ) {
00263 G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
00264 << " DiQuark/Quark = "
00265 << black->GetParticleName() << " / "
00266 << white->GetParticleName()
00267 << " resulting Hadron " << BarionDef->GetParticleName()
00268 << G4endl;
00269 }
00270 #endif
00271
00272 return BarionDef;
00273 }