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