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