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