[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 | #include "G4QGSMSplitableHadron.hh" |
---|
| 27 | #include "G4ParticleTable.hh" |
---|
| 28 | #include "G4PionPlus.hh" |
---|
| 29 | #include "G4PionMinus.hh" |
---|
| 30 | #include "G4Gamma.hh" |
---|
| 31 | #include "G4PionZero.hh" |
---|
| 32 | #include "G4KaonPlus.hh" |
---|
| 33 | #include "G4KaonMinus.hh" |
---|
| 34 | |
---|
| 35 | // based on prototype by Maxim Komogorov |
---|
| 36 | // Splitting into methods, and centralizing of model parameters HPW Feb 1999 |
---|
| 37 | // restructuring HPW Feb 1999 |
---|
| 38 | // fixing bug in the sampling of 'x', HPW Feb 1999 |
---|
| 39 | // fixing bug in sampling pz, HPW Feb 1999. |
---|
| 40 | // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999. |
---|
| 41 | // Using Parton more directly, HPW Feb 1999. |
---|
| 42 | // Shortening the algorithm for sampling x, HPW Feb 1999. |
---|
| 43 | // sampling of x replaced by formula, taking X_min into account in the correlated sampling. HPW, Feb 1999. |
---|
| 44 | // logic much clearer now. HPW Feb 1999 |
---|
| 45 | // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99. |
---|
| 46 | // Fixing p-t distributions for scattering of nuclei. |
---|
| 47 | // Separating out parameters. |
---|
| 48 | |
---|
| 49 | void G4QGSMSplitableHadron::InitParameters() |
---|
| 50 | { |
---|
| 51 | // changing rapidity distribution for all |
---|
| 52 | alpha = -0.5; // Note that this number is still assumed in the algorithm |
---|
| 53 | // needs to be generalized. |
---|
| 54 | // changing rapidity distribution for projectile like |
---|
| 55 | beta = 2.5;// Note that this number is still assumed in the algorithm |
---|
| 56 | // needs to be generalized. |
---|
| 57 | theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass(); |
---|
| 58 | // theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass(); |
---|
| 59 | // theMinPz = G4PionMinus::PionMinus()->GetPDGMass(); |
---|
| 60 | // as low as possible, otherwise, we have unphysical boundary conditions in the sampling. |
---|
| 61 | StrangeSuppress = 0.48; |
---|
| 62 | sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7, |
---|
| 63 | // but Maxim's algorithm breaks energy conservation |
---|
| 64 | // to be revised. |
---|
| 65 | widthOfPtSquare = 0.01*GeV*GeV; |
---|
| 66 | Direction = FALSE; |
---|
| 67 | minTransverseMass = 1*keV; |
---|
| 68 | } |
---|
| 69 | |
---|
| 70 | G4QGSMSplitableHadron::G4QGSMSplitableHadron() |
---|
| 71 | { |
---|
| 72 | InitParameters(); |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection) |
---|
| 76 | :G4VSplitableHadron(aPrimary) |
---|
| 77 | { |
---|
| 78 | InitParameters(); |
---|
| 79 | Direction = aDirection; |
---|
| 80 | } |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary) |
---|
| 84 | : G4VSplitableHadron(aPrimary) |
---|
| 85 | { |
---|
| 86 | InitParameters(); |
---|
| 87 | } |
---|
| 88 | |
---|
| 89 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon) |
---|
| 90 | : G4VSplitableHadron(aNucleon) |
---|
| 91 | { |
---|
| 92 | InitParameters(); |
---|
| 93 | } |
---|
| 94 | |
---|
| 95 | G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection) |
---|
| 96 | : G4VSplitableHadron(aNucleon) |
---|
| 97 | { |
---|
| 98 | InitParameters(); |
---|
| 99 | Direction = aDirection; |
---|
| 100 | } |
---|
| 101 | |
---|
| 102 | G4QGSMSplitableHadron::~G4QGSMSplitableHadron(){} |
---|
| 103 | |
---|
| 104 | const G4QGSMSplitableHadron & G4QGSMSplitableHadron::operator=(const G4QGSMSplitableHadron &) |
---|
| 105 | { |
---|
| 106 | throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron::operator= meant to not be accessable"); |
---|
| 107 | return *this; |
---|
| 108 | } |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | //************************************************************************************************************************** |
---|
| 112 | |
---|
| 113 | void G4QGSMSplitableHadron::SplitUp() |
---|
| 114 | { |
---|
| 115 | if (IsSplit()) return; |
---|
| 116 | Splitting(); |
---|
| 117 | if (Color.size()!=0) return; |
---|
| 118 | if (GetSoftCollisionCount() == 0) |
---|
| 119 | { |
---|
| 120 | DiffractiveSplitUp(); |
---|
| 121 | } |
---|
| 122 | else |
---|
| 123 | { |
---|
| 124 | SoftSplitUp(); |
---|
| 125 | } |
---|
| 126 | } |
---|
| 127 | |
---|
| 128 | void G4QGSMSplitableHadron::DiffractiveSplitUp() |
---|
| 129 | { |
---|
| 130 | // take the particle definitions and get the partons HPW |
---|
| 131 | G4Parton * Left = NULL; |
---|
| 132 | G4Parton * Right = NULL; |
---|
| 133 | GetValenceQuarkFlavors(GetDefinition(), Left, Right); |
---|
| 134 | Left->SetPosition(GetPosition()); |
---|
| 135 | Right->SetPosition(GetPosition()); |
---|
| 136 | |
---|
| 137 | G4LorentzVector HadronMom = Get4Momentum(); |
---|
| 138 | //std::cout << "DSU 1 - "<<HadronMom<<std::endl; |
---|
| 139 | |
---|
| 140 | // momenta of string ends |
---|
| 141 | G4double pt2 = HadronMom.perp2(); |
---|
| 142 | G4double transverseMass2 = HadronMom.plus()*HadronMom.minus(); |
---|
| 143 | G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2)); |
---|
| 144 | G4ThreeVector pt(minTransverseMass, minTransverseMass, 0); |
---|
| 145 | if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2); |
---|
| 146 | //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl; |
---|
| 147 | |
---|
| 148 | G4LorentzVector LeftMom(pt, 0.); |
---|
| 149 | G4LorentzVector RightMom; |
---|
| 150 | RightMom.setPx(HadronMom.px() - pt.x()); |
---|
| 151 | RightMom.setPy(HadronMom.py() - pt.y()); |
---|
| 152 | //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl; |
---|
| 153 | |
---|
| 154 | G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus(); |
---|
| 155 | G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus())); |
---|
| 156 | //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl; |
---|
| 157 | if (Direction) Local2 = -Local2; |
---|
| 158 | G4double RightMinus = 0.5*(Local1 + Local2); |
---|
| 159 | G4double LeftMinus = HadronMom.minus() - RightMinus; |
---|
| 160 | //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl; |
---|
| 161 | |
---|
| 162 | G4double LeftPlus = LeftMom.perp2()/LeftMinus; |
---|
| 163 | G4double RightPlus = HadronMom.plus() - LeftPlus; |
---|
| 164 | //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl; |
---|
| 165 | LeftMom.setPz(0.5*(LeftPlus - LeftMinus)); |
---|
| 166 | LeftMom.setE (0.5*(LeftPlus + LeftMinus)); |
---|
| 167 | RightMom.setPz(0.5*(RightPlus - RightMinus)); |
---|
| 168 | RightMom.setE (0.5*(RightPlus + RightMinus)); |
---|
| 169 | //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl; |
---|
| 170 | Left->Set4Momentum(LeftMom); |
---|
| 171 | Right->Set4Momentum(RightMom); |
---|
| 172 | Color.push_back(Left); |
---|
| 173 | AntiColor.push_back(Right); |
---|
| 174 | } |
---|
| 175 | |
---|
| 176 | |
---|
| 177 | void G4QGSMSplitableHadron::SoftSplitUp() |
---|
| 178 | { |
---|
| 179 | //... sample transversal momenta for sea and valence quarks |
---|
| 180 | G4double phi, pts; |
---|
| 181 | G4double SumPy = 0.; |
---|
| 182 | G4double SumPx = 0.; |
---|
| 183 | G4ThreeVector Pos = GetPosition(); |
---|
| 184 | G4int nSeaPair = GetSoftCollisionCount()-1; |
---|
| 185 | |
---|
| 186 | // here the condition,to ensure viability of splitting, also in cases |
---|
| 187 | // where difractive excitation occured together with soft scattering. |
---|
| 188 | // G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus(); |
---|
| 189 | // G4double Xmin = theMinPz/LightConeMomentum; |
---|
| 190 | G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() ); |
---|
| 191 | while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95; |
---|
| 192 | |
---|
| 193 | G4int aSeaPair; |
---|
| 194 | for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++) |
---|
| 195 | { |
---|
| 196 | // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2) |
---|
| 197 | |
---|
| 198 | G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); |
---|
| 199 | |
---|
| 200 | // BuildSeaQuark() determines quark spin, isospin and colour |
---|
| 201 | // via parton-constructor G4Parton(aPDGCode) |
---|
| 202 | |
---|
| 203 | G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair); |
---|
| 204 | |
---|
| 205 | // G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl; |
---|
| 206 | |
---|
| 207 | // G4cerr << "Parton 1: " |
---|
| 208 | // << " PDGcode: " << aPDGCode |
---|
| 209 | // << " - Name: " << aParton->GetDefinition()->GetParticleName() |
---|
| 210 | // << " - Type: " << aParton->GetDefinition()->GetParticleType() |
---|
| 211 | // << " - Spin-3: " << aParton->GetSpinZ() |
---|
| 212 | // << " - Colour: " << aParton->GetColour() << G4endl; |
---|
| 213 | |
---|
| 214 | // save colour a spin-3 for anti-quark |
---|
| 215 | |
---|
| 216 | G4int firstPartonColour = aParton->GetColour(); |
---|
| 217 | G4double firstPartonSpinZ = aParton->GetSpinZ(); |
---|
| 218 | |
---|
| 219 | SumPx += aParton->Get4Momentum().px(); |
---|
| 220 | SumPy += aParton->Get4Momentum().py(); |
---|
| 221 | Color.push_back(aParton); |
---|
| 222 | |
---|
| 223 | // create anti-quark |
---|
| 224 | |
---|
| 225 | aParton = BuildSeaQuark(true, aPDGCode, nSeaPair); |
---|
| 226 | aParton->SetSpinZ(-firstPartonSpinZ); |
---|
| 227 | aParton->SetColour(-firstPartonColour); |
---|
| 228 | |
---|
| 229 | // G4cerr << "Parton 2: " |
---|
| 230 | // << " PDGcode: " << -aPDGCode |
---|
| 231 | // << " - Name: " << aParton->GetDefinition()->GetParticleName() |
---|
| 232 | // << " - Type: " << aParton->GetDefinition()->GetParticleType() |
---|
| 233 | // << " - Spin-3: " << aParton->GetSpinZ() |
---|
| 234 | // << " - Colour: " << aParton->GetColour() << G4endl; |
---|
| 235 | // G4cerr << "------------" << G4endl; |
---|
| 236 | |
---|
| 237 | SumPx += aParton->Get4Momentum().px(); |
---|
| 238 | SumPy += aParton->Get4Momentum().py(); |
---|
| 239 | AntiColor.push_back(aParton); |
---|
| 240 | } |
---|
| 241 | // Valence quark |
---|
| 242 | G4Parton* pColorParton = NULL; |
---|
| 243 | G4Parton* pAntiColorParton = NULL; |
---|
| 244 | GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton); |
---|
| 245 | G4int ColorEncoding = pColorParton->GetPDGcode(); |
---|
| 246 | G4int AntiColorEncoding = pAntiColorParton->GetPDGcode(); |
---|
| 247 | |
---|
| 248 | pts = sigmaPt*std::sqrt(-std::log(G4UniformRand())); |
---|
| 249 | phi = 2.*pi*G4UniformRand(); |
---|
| 250 | G4double Px = pts*std::cos(phi); |
---|
| 251 | G4double Py = pts*std::sin(phi); |
---|
| 252 | SumPx += Px; |
---|
| 253 | SumPy += Py; |
---|
| 254 | |
---|
| 255 | if (ColorEncoding < 0) // use particle definition |
---|
| 256 | { |
---|
| 257 | G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0); |
---|
| 258 | pColorParton->Set4Momentum(ColorMom); |
---|
| 259 | G4LorentzVector AntiColorMom(Px, Py, 0, 0); |
---|
| 260 | pAntiColorParton->Set4Momentum(AntiColorMom); |
---|
| 261 | } |
---|
| 262 | else |
---|
| 263 | { |
---|
| 264 | G4LorentzVector ColorMom(Px, Py, 0, 0); |
---|
| 265 | pColorParton->Set4Momentum(ColorMom); |
---|
| 266 | G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0); |
---|
| 267 | pAntiColorParton->Set4Momentum(AntiColorMom); |
---|
| 268 | } |
---|
| 269 | Color.push_back(pColorParton); |
---|
| 270 | AntiColor.push_back(pAntiColorParton); |
---|
| 271 | |
---|
| 272 | // Sample X |
---|
| 273 | G4int nAttempt = 0; |
---|
| 274 | G4double SumX = 0; |
---|
| 275 | G4double aBeta = beta; |
---|
| 276 | G4double ColorX, AntiColorX; |
---|
| 277 | G4double HPWtest = 0; |
---|
| 278 | if (GetDefinition() == G4PionMinus::PionMinusDefinition()) aBeta = 1.; |
---|
| 279 | if (GetDefinition() == G4Gamma::GammaDefinition()) aBeta = 1.; |
---|
| 280 | if (GetDefinition() == G4PionPlus::PionPlusDefinition()) aBeta = 1.; |
---|
| 281 | if (GetDefinition() == G4PionZero::PionZeroDefinition()) aBeta = 1.; |
---|
| 282 | if (GetDefinition() == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.; |
---|
| 283 | if (GetDefinition() == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.; |
---|
| 284 | do |
---|
| 285 | { |
---|
| 286 | SumX = 0; |
---|
| 287 | nAttempt++; |
---|
| 288 | G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair; |
---|
| 289 | G4double beta1 = beta; |
---|
| 290 | if (std::abs(ColorEncoding) <= 1000 && std::abs(AntiColorEncoding) <= 1000) beta1 = 1.; //... in a meson |
---|
| 291 | ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
| 292 | HPWtest = ColorX; |
---|
| 293 | while (ColorX < Xmin || ColorX > 1.|| 1. - ColorX <= Xmin) {;} |
---|
| 294 | Color.back()->SetX(SumX = ColorX);// this is the valenz quark. |
---|
| 295 | for(G4int aPair = 0; aPair < nSeaPair; aPair++) |
---|
| 296 | { |
---|
| 297 | NumberOfUnsampledSeaQuarks--; |
---|
| 298 | ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
| 299 | Color[aPair]->SetX(ColorX); |
---|
| 300 | SumX += ColorX; |
---|
| 301 | NumberOfUnsampledSeaQuarks--; |
---|
| 302 | AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta); |
---|
| 303 | AntiColor[aPair]->SetX(AntiColorX); // the 'sea' partons |
---|
| 304 | SumX += AntiColorX; |
---|
| 305 | if (1. - SumX <= Xmin) break; |
---|
| 306 | } |
---|
| 307 | } |
---|
| 308 | while (1. - SumX <= Xmin); |
---|
| 309 | (*(AntiColor.end()-1))->SetX(1. - SumX); // the di-quark takes the rest, then go to momentum |
---|
| 310 | /// and here is the bug ;-) @@@@@@@@@@@@@ |
---|
| 311 | if(getenv("debug_QGSMSplitableHadron") )G4cout << "particle energy at split = "<<Get4Momentum().t()<<G4endl; |
---|
| 312 | G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus()); |
---|
| 313 | G4double lightCone2 = ((!Direction) ? Get4Momentum().plus() : Get4Momentum().minus()); |
---|
| 314 | // lightCone -= 0.5*Get4Momentum().m(); |
---|
| 315 | // hpw testing @@@@@ lightCone = 2.*Get4Momentum().t(); |
---|
| 316 | if(getenv("debug_QGSMSplitableHadron") )G4cout << "Light cone = "<<lightCone<<G4endl; |
---|
| 317 | for(aSeaPair = 0; aSeaPair < nSeaPair+1; aSeaPair++) |
---|
| 318 | { |
---|
| 319 | G4Parton* aParton = Color[aSeaPair]; |
---|
| 320 | aParton->DefineMomentumInZ(lightCone, lightCone2, Direction); |
---|
| 321 | |
---|
| 322 | aParton = AntiColor[aSeaPair]; |
---|
| 323 | aParton->DefineMomentumInZ(lightCone, lightCone2, Direction); |
---|
| 324 | } |
---|
| 325 | //--DEBUG-- cout <<G4endl<<"XSAMPLE "<<HPWtest<<G4endl; |
---|
| 326 | return; |
---|
| 327 | } |
---|
| 328 | |
---|
| 329 | void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart, G4Parton *& Parton1, G4Parton *& Parton2) |
---|
| 330 | { |
---|
| 331 | // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq. |
---|
| 332 | G4int aEnd; |
---|
| 333 | G4int bEnd; |
---|
| 334 | G4int HadronEncoding = aPart->GetPDGEncoding(); |
---|
| 335 | if (aPart->GetBaryonNumber() == 0) |
---|
| 336 | { |
---|
| 337 | theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd); |
---|
| 338 | } |
---|
| 339 | else |
---|
| 340 | { |
---|
| 341 | theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd); |
---|
| 342 | } |
---|
| 343 | |
---|
| 344 | Parton1 = new G4Parton(aEnd); |
---|
| 345 | Parton1->SetPosition(GetPosition()); |
---|
| 346 | |
---|
| 347 | // G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl; |
---|
| 348 | // G4cerr << "Parton 1: " |
---|
| 349 | // << " PDGcode: " << aEnd |
---|
| 350 | // << " - Name: " << Parton1->GetDefinition()->GetParticleName() |
---|
| 351 | // << " - Type: " << Parton1->GetDefinition()->GetParticleType() |
---|
| 352 | // << " - Spin-3: " << Parton1->GetSpinZ() |
---|
| 353 | // << " - Colour: " << Parton1->GetColour() << G4endl; |
---|
| 354 | |
---|
| 355 | Parton2 = new G4Parton(bEnd); |
---|
| 356 | Parton2->SetPosition(GetPosition()); |
---|
| 357 | |
---|
| 358 | // G4cerr << "Parton 2: " |
---|
| 359 | // << " PDGcode: " << bEnd |
---|
| 360 | // << " - Name: " << Parton2->GetDefinition()->GetParticleName() |
---|
| 361 | // << " - Type: " << Parton2->GetDefinition()->GetParticleType() |
---|
| 362 | // << " - Spin-3: " << Parton2->GetSpinZ() |
---|
| 363 | // << " - Colour: " << Parton2->GetColour() << G4endl; |
---|
| 364 | // G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl; |
---|
| 365 | |
---|
| 366 | // colour of parton 1 choosen at random by G4Parton(aEnd) |
---|
| 367 | // colour of parton 2 is the opposite: |
---|
| 368 | |
---|
| 369 | Parton2->SetColour(-(Parton1->GetColour())); |
---|
| 370 | |
---|
| 371 | // isospin-3 of both partons is handled by G4Parton(PDGCode) |
---|
| 372 | |
---|
| 373 | // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd) |
---|
| 374 | // spin-3 of parton 2 may be constrained by spin of original particle: |
---|
| 375 | |
---|
| 376 | if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin()) |
---|
| 377 | { |
---|
| 378 | Parton2->SetSpinZ(-(Parton2->GetSpinZ())); |
---|
| 379 | } |
---|
| 380 | |
---|
| 381 | // G4cerr << "Parton 2: " |
---|
| 382 | // << " PDGcode: " << bEnd |
---|
| 383 | // << " - Name: " << Parton2->GetDefinition()->GetParticleName() |
---|
| 384 | // << " - Type: " << Parton2->GetDefinition()->GetParticleType() |
---|
| 385 | // << " - Spin-3: " << Parton2->GetSpinZ() |
---|
| 386 | // << " - Colour: " << Parton2->GetColour() << G4endl; |
---|
| 387 | // G4cerr << "------------" << G4endl; |
---|
| 388 | |
---|
| 389 | } |
---|
| 390 | |
---|
| 391 | |
---|
| 392 | G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare) |
---|
| 393 | { |
---|
| 394 | G4double R; |
---|
| 395 | while((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare) {;} |
---|
| 396 | R = std::sqrt(R); |
---|
| 397 | G4double phi = twopi*G4UniformRand(); |
---|
| 398 | return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.); |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | G4Parton * G4QGSMSplitableHadron:: |
---|
| 402 | BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/) |
---|
| 403 | { |
---|
| 404 | if (isAntiQuark) aPDGCode*=-1; |
---|
| 405 | G4Parton* result = new G4Parton(aPDGCode); |
---|
| 406 | result->SetPosition(GetPosition()); |
---|
| 407 | G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX); |
---|
| 408 | G4LorentzVector a4Momentum(aPtVector, 0); |
---|
| 409 | result->Set4Momentum(a4Momentum); |
---|
| 410 | return result; |
---|
| 411 | } |
---|
| 412 | |
---|
| 413 | G4double G4QGSMSplitableHadron:: |
---|
| 414 | SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta) |
---|
| 415 | { |
---|
| 416 | G4double result; |
---|
| 417 | G4double x1, x2; |
---|
| 418 | G4double ymax = 0; |
---|
| 419 | for(G4int ii=1; ii<100; ii++) |
---|
| 420 | { |
---|
| 421 | G4double y = std::pow(1./G4double(ii), alpha); |
---|
| 422 | y *= std::pow( std::pow(1-anXmin-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea); |
---|
| 423 | y *= std::pow(1-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); |
---|
| 424 | if(y>ymax) ymax = y; |
---|
| 425 | } |
---|
| 426 | G4double y; |
---|
| 427 | G4double xMax=1-(totalSea+1)*anXmin; |
---|
| 428 | if(anXmin > xMax) |
---|
| 429 | { |
---|
| 430 | G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl; |
---|
| 431 | throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints."); |
---|
| 432 | } |
---|
| 433 | do |
---|
| 434 | { |
---|
| 435 | x1 = CLHEP::RandFlat::shoot(anXmin, xMax); |
---|
| 436 | y = std::pow(x1, alpha); |
---|
| 437 | y *= std::pow( std::pow(1-x1-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea); |
---|
| 438 | y *= std::pow(1-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); |
---|
| 439 | x2 = ymax*G4UniformRand(); |
---|
| 440 | } |
---|
| 441 | while(x2>y); |
---|
| 442 | result = x1; |
---|
| 443 | return result; |
---|
| 444 | } |
---|