Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4QGSMSplitableHadron.cc
이 파일의 문서화 페이지로 가기
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 "G4PhysicalConstants.hh"
28 #include "G4SystemOfUnits.hh"
29 #include "G4ParticleTable.hh"
30 #include "G4PionPlus.hh"
31 #include "G4PionMinus.hh"
32 #include "G4Gamma.hh"
33 #include "G4PionZero.hh"
34 #include "G4KaonPlus.hh"
35 #include "G4KaonMinus.hh"
36 
37 #include "G4Log.hh"
38 #include "G4Pow.hh"
39 
40 // based on prototype by Maxim Komogorov
41 // Splitting into methods, and centralizing of model parameters HPW Feb 1999
42 // restructuring HPW Feb 1999
43 // fixing bug in the sampling of 'x', HPW Feb 1999
44 // fixing bug in sampling pz, HPW Feb 1999.
45 // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999.
46 // Using Parton more directly, HPW Feb 1999.
47 // Shortening the algorithm for sampling x, HPW Feb 1999.
48 // 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
50 // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99.
51 // Fixing p-t distributions for scattering of nuclei.
52 // Separating out parameters.
53 
55 {
56  // changing rapidity distribution for all
57  alpha = -0.5; // Note that this number is still assumed in the algorithm
58  // needs to be generalized.
59  // changing rapidity distribution for projectile like
60  beta = 2.5;// Note that this number is still assumed in the algorithm
61  // needs to be generalized.
63  // theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass();
64  // theMinPz = G4PionMinus::PionMinus()->GetPDGMass();
65  // as low as possible, otherwise, we have unphysical boundary conditions in the sampling.
66  StrangeSuppress = 0.48;
67  sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7,
68  // but Maxim's algorithm breaks energy conservation
69  // to be revised.
70  widthOfPtSquare = 0.5*sqr(GeV); // 0.01*GeV*GeV; // Uzhi Apr. 2016
71  Direction = FALSE;
73  iP =0;// Color.begin();
74  iAP =0;// AntiColor.begin();
75 }
76 
78 {
80 }
81 
83 :G4VSplitableHadron(aPrimary)
84 {
86  Direction = aDirection;
87 }
88 
89 
91 : G4VSplitableHadron(aPrimary)
92 {
94 }
95 
97 : G4VSplitableHadron(aNucleon)
98 {
100 }
101 
103 : G4VSplitableHadron(aNucleon)
104 {
105  InitParameters();
106  Direction = aDirection;
107 }
108 
110 {
111 /*
112 G4cout<<"Destructor "<<Color.size()<<" "<<AntiColor.size()<<G4endl;
113  for(unsigned int i=0; i<Color.size();i++) {
114 G4cout<<"i "<<i<<G4endl;
115  delete Color.operator[](i);
116  delete AntiColor.operator[](i);
117  }
118 G4cout<<"empty"<<G4endl;
119  while(!Color.empty()) {Color.pop_back();}
120  while(!AntiColor.empty()) {AntiColor.pop_back();}
121 G4cout<<"clear"<<G4endl;
122  Color.clear(); AntiColor.clear();
123 */
124 }
125 
126 
127 
128 //**************************************************************************************************************************
129 
131 {
132 //G4cout<<G4endl<<"SplitUp() this "<<this<<" IsSplit() "<<IsSplit()<<G4endl;
133  if (IsSplit()) return;
134  Splitting(); // Uzhi To mark that a hadron is split
135 //G4cout<<"Color.size() "<<Color.size()<<G4endl;
136  if (Color.size()!=0) return;
137 //G4cout<<"GetSoftCollisionCount() "<<GetSoftCollisionCount()<<G4endl;
138  if (GetSoftCollisionCount() == 0) // GetSoftCollisionCount() from G4VSplitableHadron
139  {
141  }
142  else
143  {
144  SoftSplitUp();
145  }
146 //G4cout<<"Color.size() "<<Color.size()<<G4endl;
147 }
148 
150 {
151 /*
152 G4cout<<G4endl<<"G4QGSMSplitableHadron::DiffractiveSplitUp() "<<GetDefinition()->GetParticleName()<<G4endl;
153 G4cout<<" GetSoftCollisionCount() "<<GetSoftCollisionCount()<<G4endl;
154 G4cout<<"Mom M "<<Get4Momentum()<<" "<<Get4Momentum().mag()<<G4endl;
155 G4cout<<"Direction "<<Direction<<"---"<<G4endl;
156 */
157  // take the particle definitions and get the partons HPW
158  G4Parton * Left = NULL;
159  G4Parton * Right = NULL;
160  GetValenceQuarkFlavors(GetDefinition(), Left, Right);
161  Left->SetPosition(GetPosition());
162  Right->SetPosition(GetPosition());
163 
164 //G4cout<<"Partons Left Right "<<Left->GetDefinition()->GetParticleName()<<" "<<Right->GetDefinition()->GetParticleName()<<G4endl;
165 /*
166 G4LorentzVector tmp(0., 0., 0., 0.);
167 Left->Set4Momentum(tmp);
168 Right->Set4Momentum(tmp);
169 Color.push_back(Left);
170 AntiColor.push_back(Right);
171 */ // Uzhi
172  G4LorentzVector HadronMom = Get4Momentum();
173 //std::cout << "DSU 1 - "<<HadronMom<<std::endl;
174 
175  // momenta of string ends
176 // G4double pt2 = HadronMom.perp2();
177 // G4double transverseMass2 = HadronMom.plus()*HadronMom.minus();
178 // G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2)); // It is wrong! Uzhi
179 G4double maxAvailMomentum2 = sqr(HadronMom.mag()/2.); // Uzhi
180 //G4cout<<"Hadron M M estimated Pt "<<HadronMom.mag()<<" "<<std::sqrt(transverseMass2) - std::sqrt(pt2)<<" "<<std::sqrt(pt2)<<G4endl;
182 //G4cout<<"maxAvailMomentum2 widthOfPtSquare "<<maxAvailMomentum2<<" "<<widthOfPtSquare<<G4endl;
183  if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2);
184  //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl;
185 
186  G4LorentzVector LeftMom(pt, 0.);
187  G4LorentzVector RightMom;
188  RightMom.setPx(HadronMom.px() - pt.x());
189  RightMom.setPy(HadronMom.py() - pt.y());
190 //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl;
191 
192  G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus();
193  G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus()));
194 
195 //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl;
196  if (Direction) Local2 = -Local2;
197  G4double RightMinus = 0.5*(Local1 + Local2);
198  G4double LeftMinus = HadronMom.minus() - RightMinus;
199 //
200  if(LeftMinus <= 0.) { // Uzhi
201  RightMinus = 0.5*(Local1 - Local2); // Uzhi
202  LeftMinus = HadronMom.minus() - RightMinus; // Uzhi
203  } // Uzhi
204 //
205 //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl;
206 
207  G4double LeftPlus = LeftMom.perp2()/LeftMinus;
208  G4double RightPlus = HadronMom.plus() - LeftPlus;
209  //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl;
210  LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
211  LeftMom.setE (0.5*(LeftPlus + LeftMinus));
212  RightMom.setPz(0.5*(RightPlus - RightMinus));
213  RightMom.setE (0.5*(RightPlus + RightMinus));
214  //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl;
215  Left->Set4Momentum(LeftMom);
216  Right->Set4Momentum(RightMom);
217 //G4cout<<"Momenta H q AntiQ"<<G4endl;
218 //G4cout<<Get4Momentum()<<G4endl<<Left->Get4Momentum()<<G4endl<<Right->Get4Momentum()<<G4endl;
219 //G4cout<<"Color AntiColor "<<Left<<" "<<Right<<G4endl;
220  Color.push_back(Left);
221  AntiColor.push_back(Right);
222  iP=0; iAP=0;
223 
224 }
225 
226 
228 {
229 /*
230 G4cout<<"G4QGSMSplitableHadron::SoftSplitUp()"<<G4endl;
231 G4cout<<" GetSoftCollisionCount() "<<GetSoftCollisionCount()<<G4endl;
232 */
233  //... sample transversal momenta for sea and valence quarks
234 /* Uzhi
235  G4double phi, pts;
236  G4double SumPy = 0.;
237  G4double SumPx = 0.;
238  G4ThreeVector Pos = GetPosition();
239 */ // Uzhi
240  G4int nSeaPair = GetSoftCollisionCount()-1;
241 
242  G4LorentzVector tmp(0., 0., 0., 0.);
243 
244  G4int aSeaPair;
245  for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
246  {
247  // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
248  G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
249 
250  // BuildSeaQuark() determines quark spin, isospin and colour
251  // via parton-constructor G4Parton(aPDGCode)
252  G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
253 
254  G4int firstPartonColour = aParton->GetColour();
255  G4double firstPartonSpinZ = aParton->GetSpinZ();
256 
257  aParton->Set4Momentum(tmp);
258  Color.push_back(aParton);
259 
260  // create anti-quark
261 
262  aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
263  aParton->SetSpinZ(-firstPartonSpinZ);
264  aParton->SetColour(-firstPartonColour);
265  AntiColor.push_back(aParton);
266  }
267 
268  // Valence quark
269  G4Parton* pColorParton = NULL;
270  G4Parton* pAntiColorParton = NULL;
271  GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton);
272 // G4int ColorEncoding = pColorParton->GetPDGcode();
273 
274  pColorParton->Set4Momentum(tmp);
275  pAntiColorParton->Set4Momentum(tmp);
276 
277 //G4cout<<"Color AntiColor "<<pColorParton<<" "<<pAntiColorParton<<G4endl;
278  Color.push_back(pColorParton);
279  AntiColor.push_back(pAntiColorParton);
280 
281  iP=0; iAP=0;
282 
283 /* Uzhi
284  // here the condition,to ensure viability of splitting, also in cases
285  // where difractive excitation occured together with soft scattering.
286  // G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus();
287  // G4double Xmin = theMinPz/LightConeMomentum;
288  G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() );
289  while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95;
290 
291  G4int aSeaPair;
292  for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
293  {
294  // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
295 
296  G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
297 
298  // BuildSeaQuark() determines quark spin, isospin and colour
299  // via parton-constructor G4Parton(aPDGCode)
300 
301  G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
302 
303  // G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl;
304 
305  // G4cerr << "Parton 1: "
306  // << " PDGcode: " << aPDGCode
307  // << " - Name: " << aParton->GetDefinition()->GetParticleName()
308  // << " - Type: " << aParton->GetDefinition()->GetParticleType()
309  // << " - Spin-3: " << aParton->GetSpinZ()
310  // << " - Colour: " << aParton->GetColour() << G4endl;
311 
312  // save colour a spin-3 for anti-quark
313 
314  G4int firstPartonColour = aParton->GetColour();
315  G4double firstPartonSpinZ = aParton->GetSpinZ();
316 
317  SumPx += aParton->Get4Momentum().px();
318  SumPy += aParton->Get4Momentum().py();
319  Color.push_back(aParton);
320 
321  // create anti-quark
322 
323  aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
324  aParton->SetSpinZ(-firstPartonSpinZ);
325  aParton->SetColour(-firstPartonColour);
326 
327  // G4cerr << "Parton 2: "
328  // << " PDGcode: " << -aPDGCode
329  // << " - Name: " << aParton->GetDefinition()->GetParticleName()
330  // << " - Type: " << aParton->GetDefinition()->GetParticleType()
331  // << " - Spin-3: " << aParton->GetSpinZ()
332  // << " - Colour: " << aParton->GetColour() << G4endl;
333  // G4cerr << "------------" << G4endl;
334 
335  SumPx += aParton->Get4Momentum().px();
336  SumPy += aParton->Get4Momentum().py();
337  AntiColor.push_back(aParton);
338  }
339 */ // Uzhi
340 /* Uzhi
341  // Valence quark
342  G4Parton* pColorParton = NULL;
343  G4Parton* pAntiColorParton = NULL;
344  GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton);
345  G4int ColorEncoding = pColorParton->GetPDGcode();
346 
347  pts = sigmaPt*std::sqrt(-G4Log(G4UniformRand()));
348  phi = 2.*pi*G4UniformRand();
349  G4double Px = pts*std::cos(phi);
350  G4double Py = pts*std::sin(phi);
351  SumPx += Px;
352  SumPy += Py;
353 
354  if (ColorEncoding < 0) // use particle definition
355  {
356  G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0);
357  pColorParton->Set4Momentum(ColorMom);
358  G4LorentzVector AntiColorMom(Px, Py, 0, 0);
359  pAntiColorParton->Set4Momentum(AntiColorMom);
360  }
361  else
362  {
363  G4LorentzVector ColorMom(Px, Py, 0, 0);
364  pColorParton->Set4Momentum(ColorMom);
365  G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0);
366  pAntiColorParton->Set4Momentum(AntiColorMom);
367  }
368  Color.push_back(pColorParton);
369  AntiColor.push_back(pAntiColorParton);
370 
371  // Sample X
372  G4int nAttempt = 0;
373  G4double SumX = 0;
374  G4double aBeta = beta;
375  G4double ColorX, AntiColorX;
376  if (GetDefinition() == G4PionMinus::PionMinusDefinition()) aBeta = 1.;
377  if (GetDefinition() == G4Gamma::GammaDefinition()) aBeta = 1.;
378  if (GetDefinition() == G4PionPlus::PionPlusDefinition()) aBeta = 1.;
379  if (GetDefinition() == G4PionZero::PionZeroDefinition()) aBeta = 1.;
380  if (GetDefinition() == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.;
381  if (GetDefinition() == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.;
382  do
383  {
384  SumX = 0;
385  nAttempt++;
386  G4int NumberOfUnsampledSeaQuarks = 2*nSeaPair;
387  ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
388  Color.back()->SetX(SumX = ColorX);// this is the valenz quark.
389  for(G4int aPair = 0; aPair < nSeaPair; aPair++)
390  {
391  NumberOfUnsampledSeaQuarks--;
392  ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
393  Color[aPair]->SetX(ColorX);
394  SumX += ColorX;
395  NumberOfUnsampledSeaQuarks--;
396  AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
397  AntiColor[aPair]->SetX(AntiColorX); // the 'sea' partons
398  SumX += AntiColorX;
399  if (1. - SumX <= Xmin) break;
400  }
401  }
402  while (1. - SumX <= Xmin);
403 
404  (*(AntiColor.end()-1))->SetX(1. - SumX); // the di-quark takes the rest, then go to momentum
405  G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus());
406  G4double lightCone2 = ((!Direction) ? Get4Momentum().plus() : Get4Momentum().minus());
407  for(aSeaPair = 0; aSeaPair < nSeaPair+1; aSeaPair++)
408  {
409  G4Parton* aParton = Color[aSeaPair];
410  aParton->DefineMomentumInZ(lightCone, lightCone2, Direction);
411 
412  aParton = AntiColor[aSeaPair];
413  aParton->DefineMomentumInZ(lightCone, lightCone2, Direction);
414  }
415 */ // Uzhi
416  return;
417 }
418 
420 {
421  // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
422  G4int aEnd;
423  G4int bEnd;
424  G4int HadronEncoding = aPart->GetPDGEncoding();
425  if (aPart->GetBaryonNumber() == 0)
426  {
427  theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd);
428  }
429  else
430  {
431  theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd);
432  }
433 
434  Parton1 = new G4Parton(aEnd);
435  Parton1->SetPosition(GetPosition());
436 
437  // G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
438  // G4cerr << "Parton 1: "
439  // << " PDGcode: " << aEnd
440  // << " - Name: " << Parton1->GetDefinition()->GetParticleName()
441  // << " - Type: " << Parton1->GetDefinition()->GetParticleType()
442  // << " - Spin-3: " << Parton1->GetSpinZ()
443  // << " - Colour: " << Parton1->GetColour() << G4endl;
444 
445  Parton2 = new G4Parton(bEnd);
446  Parton2->SetPosition(GetPosition());
447 
448  // G4cerr << "Parton 2: "
449  // << " PDGcode: " << bEnd
450  // << " - Name: " << Parton2->GetDefinition()->GetParticleName()
451  // << " - Type: " << Parton2->GetDefinition()->GetParticleType()
452  // << " - Spin-3: " << Parton2->GetSpinZ()
453  // << " - Colour: " << Parton2->GetColour() << G4endl;
454  // G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;
455 
456  // colour of parton 1 choosen at random by G4Parton(aEnd)
457  // colour of parton 2 is the opposite:
458 
459  Parton2->SetColour(-(Parton1->GetColour()));
460 
461  // isospin-3 of both partons is handled by G4Parton(PDGCode)
462 
463  // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd)
464  // spin-3 of parton 2 may be constrained by spin of original particle:
465 
466  if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin())
467  {
468  Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
469  }
470 
471  // G4cerr << "Parton 2: "
472  // << " PDGcode: " << bEnd
473  // << " - Name: " << Parton2->GetDefinition()->GetParticleName()
474  // << " - Type: " << Parton2->GetDefinition()->GetParticleType()
475  // << " - Spin-3: " << Parton2->GetSpinZ()
476  // << " - Colour: " << Parton2->GetColour() << G4endl;
477  // G4cerr << "------------" << G4endl;
478 
479 }
480 
481 
483 {
484  G4double R;
485  const G4int maxNumberOfLoops = 1000;
486  G4int loopCounter = -1;
487  while( ((R = -widthSquare*G4Log(G4UniformRand())) > maxPtSquare) &&
488  ++loopCounter < maxNumberOfLoops ) {;} /* Loop checking, 07.08.2015, A.Ribon */
489  if ( loopCounter >= maxNumberOfLoops ) {
490  R = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration.
491  }
492  R = std::sqrt(R);
493  G4double phi = twopi*G4UniformRand();
494  return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.);
495 }
496 
498 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/)
499 {
500  if (isAntiQuark) aPDGCode*=-1;
501  G4Parton* result = new G4Parton(aPDGCode);
502  result->SetPosition(GetPosition());
503  G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
504  G4LorentzVector a4Momentum(aPtVector, 0);
505  result->Set4Momentum(a4Momentum);
506  return result;
507 }
508 
510 SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
511 {
513  G4double x1, x2;
514  G4double ymax = 0;
515  for(G4int ii=1; ii<100; ii++)
516  {
518  y *= G4Pow::GetInstance()->powN( G4Pow::GetInstance()->powA(1-anXmin-totalSea*anXmin, alpha+1) -
519  G4Pow::GetInstance()->powA(anXmin, alpha+1), nSea);
520  y *= G4Pow::GetInstance()->powA(1-anXmin-totalSea*anXmin, aBeta+1) -
521  G4Pow::GetInstance()->powA(anXmin, aBeta+1);
522  if(y>ymax) ymax = y;
523  }
524  G4double y;
525  G4double xMax=1-(totalSea+1)*anXmin;
526  if(anXmin > xMax)
527  {
528 // G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl;
529  throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints.");
530  }
531  const G4int maxNumberOfLoops = 1000;
532  G4int loopCounter = 0;
533  do
534  {
535  x1 = G4RandFlat::shoot(anXmin, xMax);
536  y = G4Pow::GetInstance()->powA(x1, alpha);
537  y *= G4Pow::GetInstance()->powN( G4Pow::GetInstance()->powA(1-x1-totalSea*anXmin, alpha+1) -
538  G4Pow::GetInstance()->powA(anXmin, alpha+1), nSea);
539  y *= G4Pow::GetInstance()->powA(1-x1-totalSea*anXmin, aBeta+1) -
540  G4Pow::GetInstance()->powA(anXmin, aBeta+1);
541  x2 = ymax*G4UniformRand();
542  }
543  while( (x2>y) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
544  if ( loopCounter >= maxNumberOfLoops ) {
545  x1 = 0.5*( anXmin + xMax ); // Just an acceptable value, without any physics consideration.
546  }
547  result = x1;
548  return result;
549 }
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4ParticleDefinition * GetDefinition() const
void SetSpinZ(G4double aSpinZ)
Definition: G4Parton.hh:95
CLHEP::Hep3Vector G4ThreeVector
double plus() const
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare)
double py() const
static constexpr double keV
Definition: G4SIunits.hh:216
Float_t x1[n_points_granero]
Definition: compare.C:5
void GetValenceQuarkFlavors(const G4ParticleDefinition *aPart, G4Parton *&Parton1, G4Parton *&Parton2)
Float_t y
Definition: compare.C:6
double px() const
G4int GetColour()
Definition: G4Parton.hh:90
G4bool SplitMeson(G4int PDGcode, G4int *aEnd, G4int *bEnd)
Float_t tmp
G4double GetPDGMass() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:177
double perp2() const
void SetColour(G4int aColour)
Definition: G4Parton.hh:89
G4double GetPDGSpin() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
const G4ThreeVector & GetPosition() const
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
double minus() const
#define FALSE
Definition: globals.hh:52
#define G4UniformRand()
Definition: Randomize.hh:53
G4double SampleX(G4double anXmin, G4int nSea, G4int theTotalSea, G4double aBeta)
static constexpr double twopi
Definition: G4SIunits.hh:76
TMarker * pt
Definition: egs.C:25
Double_t R
G4bool SplitBarion(G4int PDGCode, G4int *q_or_qqbar, G4int *qbar_or_qq)
G4Parton * BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int nSeaPair)
G4BaryonSplitter theBaryonSplitter
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double G4ParticleHPJENDLHEData::G4double result
void SetPosition(const G4ThreeVector &aPosition)
Definition: G4Parton.hh:137
std::deque< G4Parton * > AntiColor
Color
Definition: test07.cc:36
double mag() const
int G4int
Definition: G4Types.hh:78
G4MesonSplitter theMesonSplitter
double x() const
const G4LorentzVector & Get4Momentum() const
T sqr(const T &x)
Definition: templates.hh:145
double y() const
Float_t x2[n_points_geant4]
Definition: compare.C:26
#define DBL_MAX
Definition: templates.hh:83
G4double GetSpinZ()
Definition: G4Parton.hh:96
static constexpr double GeV
Definition: G4SIunits.hh:217
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148