Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QHadron.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 //
27 // $Id$
28 //
29 // ---------------- G4QHadron ----------------
30 // by Mikhail Kossov, Sept 1999.
31 // class for Quasmon initiated Hadrons generated by CHIPS Model
32 // -------------------------------------------------------------------
33 // Short description: In CHIPS all particles are G4QHadrons, while they
34 // can be leptons, gammas or nuclei. The G4QPDGCode makes the difference.
35 // In addition the 4-momentum is a basic value, so the mass can be
36 // different from the GS mass (e.g. for the virtual gamma).
37 // -------------------------------------------------------------------
38 //
39 //#define debug
40 //#define edebug
41 //#define pdebug
42 //#define sdebug
43 //#define ppdebug
44 
45 #include <cmath>
46 
47 #include "G4QHadron.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 
51 using namespace std;
52 
53 G4double G4QHadron::StrangeSuppress = 0.48; // ? M.K.
54 G4double G4QHadron::sigmaPt = 1.7*GeV; // Can be 0 ?
55 G4double G4QHadron::widthOfPtSquare = 0.01*GeV*GeV; // ? M.K.
56 
57 G4QHadron::G4QHadron(): theMomentum(0.,0.,0.,0.), theQPDG(0), valQ(0,0,0,0,0,0), nFragm(0),
58  thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
59  Color(), AntiColor(), bindE(0.), formTime(0.) {}
60 
61 G4QHadron::G4QHadron(G4LorentzVector p): theMomentum(p), theQPDG(0), valQ(0,0,0,0,0,0),
62  nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
63  Color(), AntiColor(), bindE(0.), formTime(0.) {}
64 
65 // For Chipolino or Quasmon doesn't make any sense
66 G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p): theMomentum(p), theQPDG(PDGCode),
67  nFragm(0),thePosition(0.,0.,0.),theCollisionCount(0),isSplit(false),Direction(true),
68  Color(), AntiColor(), bindE(0.), formTime(0.)
69 {
70 #ifdef debug
71  G4cout<<"G4QHadron must be created with PDG="<<PDGCode<<", 4M="<<p<<G4endl;
72 #endif
73  if(GetQCode()>-1)
74  {
75  if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass());
76  valQ=theQPDG.GetQuarkContent();
77  }
78  else if(PDGCode>80000000) DefineQC(PDGCode);
79  else G4cerr<<"***G4QHadron:(P) PDG="<<PDGCode<<", use other constructor"<<G4endl;
80 #ifdef debug
81  G4cout<<"G4QHadron is created with QCode="<<GetQCode()<<", QC="<<valQ<<G4endl;
82 #endif
83 }
84 
85 // For Chipolino or Quasmon doesn't make any sense
86 G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p): theMomentum(p), theQPDG(QPDG),
87  nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
88  Color(), AntiColor(), bindE(0.), formTime(0.)
89 {
90  if(theQPDG.GetQCode()>-1)
91  {
92  if(theMomentum.e()==0.) theMomentum.setE(theQPDG.GetMass());
93  valQ=theQPDG.GetQuarkContent();
94  }
95  else
96  {
97  G4int cPDG=theQPDG.GetPDGCode();
98  if(cPDG>80000000) DefineQC(cPDG);
99  else G4cerr<<"***G4QHadr:(QP) PDG="<<cPDG<<" use other constructor"<<G4endl;
100  }
101 }
102 
103 // Make sense Chipolino or Quasmon
104 G4QHadron::G4QHadron(G4QContent QC, G4LorentzVector p): theMomentum(p),theQPDG(0),valQ(QC),
105  nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0), isSplit(false), Direction(true),
106  Color(), AntiColor(), bindE(0.), formTime(0.)
107 {
108  G4int curPDG=valQ.GetSPDGCode();
109  if(curPDG==10&&valQ.GetBaryonNumber()>0) curPDG=valQ.GetZNSPDGCode();
110  if(curPDG&&curPDG!=10) theQPDG.SetPDGCode(curPDG);
111  else theQPDG.InitByQCont(QC);
112 }
113 
114 G4QHadron::G4QHadron(G4int PDGCode, G4double aMass, G4QContent QC) :
115  theMomentum(0.,0.,0.,aMass), theQPDG(PDGCode), valQ(QC), nFragm(0),thePosition(0.,0.,0.),
116  theCollisionCount(0), isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.),
117  formTime(0.)
118 {}
119 
120 G4QHadron::G4QHadron(G4QPDGCode QPDG, G4double aMass, G4QContent QC) :
121  theMomentum(0.,0.,0.,aMass), theQPDG(QPDG), valQ(QC), nFragm(0), thePosition(0.,0.,0.),
122  theCollisionCount(0), isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.),
123  formTime(0.)
124 {}
125 
126 G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p, G4QContent QC) : theMomentum(p),
127  theQPDG(PDGCode), valQ(QC), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
128  isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
129 {}
130 
131 G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p, G4QContent QC) : theMomentum(p),
132  theQPDG(QPDG), valQ(QC), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
133  isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
134 {}
135 
136 G4QHadron::G4QHadron(G4QParticle* pPart, G4double maxM) : theMomentum(0.,0.,0.,0.),
137  theQPDG(pPart->GetQPDG()), nFragm(0), thePosition(0.,0.,0.), theCollisionCount(0),
138  isSplit(false), Direction(true), Color(), AntiColor(), bindE(0.), formTime(0.)
139 {
140 #ifdef debug
141  G4cout<<"G4QHadron is created & randomized with maxM="<<maxM<<G4endl;
142 #endif
143  G4int PDGCode = theQPDG.GetPDGCode();
144  if(PDGCode<2)G4cerr<<"***G4QHadron:(M) PDGC="<<PDGCode<<" use other constructor"<<G4endl;
145  valQ=theQPDG.GetQuarkContent();
146  theMomentum.setE(RandomizeMass(pPart, maxM));
147 }
148 
149 G4QHadron::G4QHadron(const G4QHadron& right)
150 {
151  theMomentum = right.theMomentum;
152  theQPDG = right.theQPDG;
153  valQ = right.valQ;
154  nFragm = right.nFragm;
155  thePosition = right.thePosition;
156  theCollisionCount = 0;
157  isSplit = false;
158  Direction = right.Direction;
159  bindE = right.bindE;
160  formTime = right.formTime;
161 }
162 
163 G4QHadron::G4QHadron(const G4QHadron* right)
164 {
165  theMomentum = right->theMomentum;
166  theQPDG = right->theQPDG;
167  valQ = right->valQ;
168  nFragm = right->nFragm;
169  thePosition = right->thePosition;
170  theCollisionCount = 0;
171  isSplit = false;
172  Direction = right->Direction;
173  bindE = right->bindE;
174  formTime = right->formTime;
175 }
176 
177 G4QHadron::G4QHadron(const G4QHadron* right, G4int C, G4ThreeVector P, G4LorentzVector M)
178 {
179  theMomentum = M;
180  theQPDG = right->theQPDG;
181  valQ = right->valQ;
182  nFragm = right->nFragm;
183  thePosition = P;
184  theCollisionCount = C;
185  isSplit = false;
186  Direction = right->Direction;
187  bindE = right->bindE;
188  formTime = right->formTime;
189 }
190 
191 const G4QHadron& G4QHadron::operator=(const G4QHadron &right)
192 {
193  if(this != &right) // Beware of self assignment
194  {
195  theMomentum = right.theMomentum;
196  theQPDG = right.theQPDG;
197  valQ = right.valQ;
198  nFragm = right.nFragm;
199  thePosition = right.thePosition;
200  theCollisionCount = 0;
201  isSplit = false;
202  Direction = right.Direction;
203  bindE = right.bindE;
204  }
205  return *this;
206 }
207 
208 G4QHadron::~G4QHadron()
209 {
210  std::list<G4QParton*>::iterator ipos = Color.begin();
211  std::list<G4QParton*>::iterator epos = Color.end();
212  for( ; ipos != epos; ipos++) {delete [] *ipos;}
213  Color.clear();
214 
215  ipos = AntiColor.begin();
216  epos = AntiColor.end();
217  for( ; ipos != epos; ipos++) {delete [] *ipos;}
218  AntiColor.clear();
219 }
220 
221 // Define quark content of the particle with a particular PDG Code
222 void G4QHadron::DefineQC(G4int PDGCode)
223 {
224  //G4cout<<"G4QHadron::DefineQC is called with PDGCode="<<PDGCode<<G4endl;
225  G4int szn=PDGCode-90000000;
226  G4int ds=0;
227  G4int dz=0;
228  G4int dn=0;
229  if(szn<-100000)
230  {
231  G4int ns_value=(-szn)/1000000+1;
232  szn+=ns_value*1000000;
233  ds+=ns_value;
234  }
235  else if(szn<-100)
236  {
237  G4int nz=(-szn)/1000+1;
238  szn+=nz*1000;
239  dz+=nz;
240  }
241  else if(szn<0)
242  {
243  G4int nn=-szn;
244  szn=0;
245  dn+=nn;
246  }
247  G4int sz =szn/1000;
248  G4int n =szn%1000;
249  if(n>700)
250  {
251  n-=1000;
252  dz--;
253  }
254  G4int z =sz%1000-dz;
255  if(z>700)
256  {
257  z-=1000;
258  ds--;
259  }
260  G4int Sq =sz/1000-ds;
261  G4int zns=z+n+Sq;
262  G4int Dq=n+zns;
263  G4int Uq=z+zns;
264  if (Dq<0&&Uq<0&&Sq<0)valQ=G4QContent(0 ,0 ,0 ,-Dq,-Uq,-Sq);
265  else if (Uq<0&&Sq<0) valQ=G4QContent(Dq,0 ,0 ,0 ,-Uq,-Sq);
266  else if (Dq<0&&Sq<0) valQ=G4QContent(0 ,Uq,0 ,-Dq,0 ,-Sq);
267  else if (Dq<0&&Uq<0) valQ=G4QContent(0 ,0 ,Sq,-Dq,-Uq,0 );
268  else if (Uq<0) valQ=G4QContent(Dq,0 ,Sq,0 ,-Uq,0 );
269  else if (Sq<0) valQ=G4QContent(Dq,Uq,0 ,0 ,0 ,-Sq);
270  else if (Dq<0) valQ=G4QContent(0 ,Uq,Sq,-Dq,0 ,0 );
271  else valQ=G4QContent(Dq,Uq,Sq,0 ,0 ,0 );
272 }
273 
274 // Redefine a Hadron with a new PDGCode
275 void G4QHadron::SetQPDG(const G4QPDGCode& newQPDG)
276 {
277  theQPDG = newQPDG;
278  G4int PDG= newQPDG.GetPDGCode();
279  G4int Q = newQPDG.GetQCode();
280 #ifdef debug
281  G4cout<<"G4QHadron::SetQPDG is called with PDGCode="<<PDG<<", QCode="<<Q<<G4endl;
282 #endif
283  if (Q>-1) valQ=theQPDG.GetQuarkContent();
284  else if(PDG>80000000) DefineQC(PDG);
285  else
286  {
287  // G4cerr<<"***G4QHadron::SetQPDG: QPDG="<<newQPDG<<G4endl;
288  // throw G4QException("***G4QHadron::SetQPDG: Impossible QPDG Probably a Chipolino");
289  G4ExceptionDescription ed;
290  ed << "Impossible QPDG Probably a Chipolino: QPDG=" << newQPDG << G4endl;
291  G4Exception("G4QHadron::SetQPDG()", "HAD_CHPS_0000", FatalException, ed);
292  }
293 }
294 
295 // Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir
296 G4bool G4QHadron::RelDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
297  G4LorentzVector& dir, G4double maxCost, G4double minCost)
298 {
299  G4double fM2 = f4Mom.m2();
300  G4double fM = sqrt(fM2); // Mass of the 1st Hadron
301  G4double sM2 = s4Mom.m2();
302  G4double sM = sqrt(sM2); // Mass of the 2nd Hadron
303  G4double iM2 = theMomentum.m2();
304  G4double iM = sqrt(iM2); // Mass of the decaying hadron
305  G4double vP = theMomentum.rho(); // Momentum of the decaying hadron
306  G4double dE = theMomentum.e(); // Energy of the decaying hadron
307  if(dE<vP)
308  {
309  G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
310  G4double accuracy=.000001*vP;
311  G4double emodif=std::fabs(dE-vP);
312  //if(emodif<accuracy)
313  //{
314  G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
315  theMomentum.setE(vP+emodif+.01*accuracy);
316  //}
317  }
318  G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
319  G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans.
320  G4LorentzVector cdir = dir; // A copy to make a transformation to CMS
321 #ifdef ppdebug
322  if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p="
323  <<cdir.e()-cdir.rho()<<G4endl;
324 #endif
325  cdir.boost(ltf); // Direction transpormed to CMS of the Momentum
326  G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle
327 #ifdef ppdebug
328  G4cout<<"G4QHad::RelDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl;
329 #endif
330  G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle
331  G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction
332  G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction
333  if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS
334  {
335  vx = vdir.unit(); // Ort in the direction of the reference particle
336 #ifdef ppdebug
337  G4cout<<"G4QH::RelDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl;
338 #endif
339  G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
340  vy = vv.unit(); // First ort orthogonal to the direction
341  vz = vx.cross(vy); // Second ort orthoganal to the direction
342  }
343 #ifdef ppdebug
344  G4cout<<"G4QHad::RelDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl;
345 #endif
346  if(maxCost> 1.) maxCost= 1.;
347  if(minCost<-1.) minCost=-1.;
348  if(maxCost<-1.) maxCost=-1.;
349  if(minCost> 1.) minCost= 1.;
350  if(minCost> maxCost) minCost=maxCost;
351  if(fabs(iM-fM-sM)<.00000001)
352  {
353  G4double fR=fM/iM;
354  G4double sR=sM/iM;
355  f4Mom=fR*theMomentum;
356  s4Mom=sR*theMomentum;
357  return true;
358  }
359  else if (iM+.001<fM+sM || iM==0.)
360  {//@@ Later on make a quark content check for the decay
361  G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
362  return false;
363  }
364  G4double d2 = iM2-fM2-sM2;
365  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon
366  if(p2<0.)
367  {
368 #ifdef ppdebug
369  G4cout<<"**G4QH:RDIn2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
370 #endif
371  p2=0.;
372  }
373  G4double p = sqrt(p2);
374  G4double ct = maxCost;
375  if(maxCost>minCost)
376  {
377  G4double dcost=maxCost-minCost;
378  ct = minCost+dcost*G4UniformRand();
379  }
380  G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?)
381  G4double ps=0.;
382  if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct);
383  else
384  {
385 #ifdef ppdebug
386  G4cout<<"**G4QH::RDIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl;
387  //throw G4QException("***G4QHadron::RDIn2: bad cos(theta)");
388 #endif
389  if(ct>1.) ct=1.;
390  if(ct<-1.) ct=-1.;
391  }
392  G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx;
393 #ifdef ppdebug
394  G4cout<<"G4QH::RelDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl;
395 #endif
396 
397  f4Mom.setVect(pVect);
398  f4Mom.setE(sqrt(fM2+p2));
399  s4Mom.setVect((-1)*pVect);
400  s4Mom.setE(sqrt(sM2+p2));
401 
402 #ifdef ppdebug
403  G4cout<<"G4QHadr::RelDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = "
404  <<f4Mom+s4Mom<<", M="<<iM<<G4endl;
405 #endif
406  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
407  <<f4Mom.e()-f4Mom.rho()<<G4endl;
408  f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS
409  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
410  <<s4Mom.e()-s4Mom.rho()<<G4endl;
411  s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS
412 #ifdef ppdebug
413  G4cout<<"G4QHadron::RelDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = "
414  <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl;
415 #endif
416  return true;
417 } // End of "RelDecayIn2"
418 
419 // Decay of Hadron In2Particles f&s, f w/r/to dN/dO [cp>0: ~cost^cp, cp<0: ~(1-cost)^(-cp)]
420 G4bool G4QHadron::CopDecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
421  G4LorentzVector& dir, G4double cosp)
422 {
423  G4double fM2 = f4Mom.m2();
424  G4double fM = sqrt(fM2); // Mass of the 1st Hadron
425  G4double sM2 = s4Mom.m2();
426  G4double sM = sqrt(sM2); // Mass of the 2nd Hadron
427  G4double iM2 = theMomentum.m2();
428  G4double iM = sqrt(iM2); // Mass of the decaying hadron
429  G4double vP = theMomentum.rho(); // Momentum of the decaying hadron
430  G4double dE = theMomentum.e(); // Energy of the decaying hadron
431  G4bool neg=false; // Negative (backward) distribution of t
432  if(cosp<0)
433  {
434  cosp=-cosp;
435  neg=true;
436  }
437  if(dE<vP)
438  {
439  G4cerr<<"***G4QHad::CopDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
440  G4double accuracy=.000001*vP;
441  G4double emodif=std::fabs(dE-vP);
442  //if(emodif<accuracy)
443  //{
444  G4cerr<<"G4QHadron::CopDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
445  theMomentum.setE(vP+emodif+.01*accuracy);
446  //}
447  }
448  G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans.
449  G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans.
450  G4LorentzVector cdir = dir; // A copy to make a transformation to CMS
451 #ifdef ppdebug
452  if(cdir.e()+.001<cdir.rho()) G4cerr<<"*G4QH::RDIn2:*Boost* cd4M="<<cdir<<",e-p="
453  <<cdir.e()-cdir.rho()<<G4endl;
454 #endif
455  cdir.boost(ltf); // Direction transpormed to CMS of the Momentum
456  G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle
457 #ifdef ppdebug
458  G4cout<<"G4QHad::CopDI2:dir="<<dir<<",ltf="<<ltf<<",cdir="<<cdir<<",vdir="<<vdir<<G4endl;
459 #endif
460  G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle
461  G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction
462  G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction
463  if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS
464  {
465  vx = vdir.unit(); // Ort in the direction of the reference particle
466 #ifdef ppdebug
467  G4cout<<"G4QH::CopDecIn2:Vx="<<vx<<",M="<<theMomentum<<",d="<<dir<<",c="<<cdir<<G4endl;
468 #endif
469  G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!)
470  vy = vv.unit(); // First ort orthogonal to the direction
471  vz = vx.cross(vy); // Second ort orthoganal to the direction
472  }
473 #ifdef ppdebug
474  G4cout<<"G4QHad::CopDecIn2:iM="<<iM<<"=>fM="<<fM<<"+sM="<<sM<<",ob="<<vx<<vy<<vz<<G4endl;
475 #endif
476  if(fabs(iM-fM-sM)<.00000001)
477  {
478  G4double fR=fM/iM;
479  G4double sR=sM/iM;
480  f4Mom=fR*theMomentum;
481  s4Mom=sR*theMomentum;
482  return true;
483  }
484  else if (iM+.001<fM+sM || iM==0.)
485  {//@@ Later on make a quark content check for the decay
486  G4cerr<<"***G4QH::CopDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl;
487  return false;
488  }
489  G4double d2 = iM2-fM2-sM2;
490  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon
491  if(p2<0.)
492  {
493 #ifdef ppdebug
494  G4cout<<"*G4QH:CopDI2:p2="<<p2<<"<0,d4/4="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
495 #endif
496  p2=0.;
497  }
498  G4double p = sqrt(p2);
499  G4double ct = 0;
500  G4double rn = pow(G4UniformRand(),cosp+1.);
501  if(neg) ct = rn+rn-1.; // More backward than forward
502  else ct = 1.-rn-rn; // More forward than backward
503  //
504  G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?)
505  G4double ps=0.;
506  if(fabs(ct)<1.) ps = p * sqrt(1.-ct*ct);
507  else
508  {
509 #ifdef ppdebug
510  G4cout<<"**G4QH::CopDecayIn2:ct="<<ct<<",mac="<<maxCost<<",mic="<<minCost<<G4endl;
511  //throw G4QException("***G4QHadron::RDIn2: bad cos(theta)");
512 #endif
513  if(ct>1.) ct=1.;
514  if(ct<-1.) ct=-1.;
515  }
516  G4ThreeVector pVect=(ps*sin(phi))*vz+(ps*cos(phi))*vy+p*ct*vx;
517 #ifdef ppdebug
518  G4cout<<"G4QH::CopDIn2:ct="<<ct<<",p="<<p<<",ps="<<ps<<",ph="<<phi<<",v="<<pVect<<G4endl;
519 #endif
520 
521  f4Mom.setVect(pVect);
522  f4Mom.setE(sqrt(fM2+p2));
523  s4Mom.setVect((-1)*pVect);
524  s4Mom.setE(sqrt(sM2+p2));
525 
526 #ifdef ppdebug
527  G4cout<<"G4QHadr::CopDecIn2:p2="<<p2<<",v="<<ltb<<",f4M="<<f4Mom<<" + s4M="<<s4Mom<<" = "
528  <<f4Mom+s4Mom<<", M="<<iM<<G4endl;
529 #endif
530  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p="
531  <<f4Mom.e()-f4Mom.rho()<<G4endl;
532  f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS
533  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p="
534  <<s4Mom.e()-s4Mom.rho()<<G4endl;
535  s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS
536 #ifdef ppdebug
537  G4cout<<"G4QHadron::CopDecayIn2:Output, f4Mom="<<f4Mom<<" + s4Mom="<<s4Mom<<" = "
538  <<f4Mom+s4Mom<<", d4M="<<theMomentum-f4Mom-s4Mom<<G4endl;
539 #endif
540  return true;
541 } // End of "CopDecayIn2"
542 
543 // Decay of the Hadron in 2 particles (f + s)
544 G4bool G4QHadron::DecayIn2(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom)
545 {
546  G4double fM2 = f4Mom.m2();
547  if(fM2<0.) fM2=0.;
548  G4double fM = sqrt(fM2); // Mass of the 1st Hadron
549  G4double sM2 = s4Mom.m2();
550  if(sM2<0.) sM2=0.;
551  G4double sM = sqrt(sM2); // Mass of the 2nd Hadron
552  G4double iM2 = theMomentum.m2();
553  if(iM2<0.) iM2=0.;
554  G4double iM = sqrt(iM2); // Mass of the decaying hadron
555 #ifdef debug
556  G4cout<<"G4QHadron::DecIn2: iM="<<iM<<" => fM="<<fM<<" + sM="<<sM<<" = "<<fM+sM<<G4endl;
557 #endif
558  //@@ Later on make a quark content check for the decay
559  if (fabs(iM-fM-sM)<.0000001)
560  {
561  G4double fR=fM/iM;
562  G4double sR=sM/iM;
563  f4Mom=fR*theMomentum;
564  s4Mom=sR*theMomentum;
565  return true;
566  }
567  else if (iM+.001<fM+sM || iM==0.)
568  {
569 #ifdef debug
570  G4cerr<<"***G4QHadron::DecayIn2*** fM="<<fM<<" + sM="<<sM<<"="<<fM+sM<<" > iM="<<iM
571  <<", d="<<iM-fM-sM<<G4endl;
572 #endif
573  return false;
574  }
575 
576  G4double d2 = iM2-fM2-sM2;
577  G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon
578  if (p2<0.)
579  {
580 #ifdef debug
581  G4cerr<<"***G4QH::DI2:p2="<<p2<<"<0,d2^2="<<d2*d2/4.<<"<4*fM2*sM2="<<4*fM2*sM2<<G4endl;
582 #endif
583  p2=0.;
584  }
585  G4double p = sqrt(p2);
586  G4double ct = 1.-2*G4UniformRand();
587 #ifdef debug
588  G4cout<<"G4QHadron::DecayIn2: ct="<<ct<<", p="<<p<<G4endl;
589 #endif
590  G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?)
591  G4double ps = p * sqrt(1.-ct*ct);
592  G4ThreeVector pVect(ps*sin(phi),ps*cos(phi),p*ct);
593 
594  f4Mom.setVect(pVect);
595  f4Mom.setE(sqrt(fM2+p2));
596  s4Mom.setVect((-1)*pVect);
597  s4Mom.setE(sqrt(sM2+p2));
598 
599  if(theMomentum.e()<theMomentum.rho())
600  {
601  G4cerr<<"*G4QH::DecIn2:*Boost* 4M="<<theMomentum<<",e-p="
602  <<theMomentum.e()-theMomentum.rho()<<G4endl;
603  //throw G4QException("G4QHadron::DecayIn2: Decay of particle with zero mass")
604  //theMomentum.setE(1.0000001*theMomentum.rho()); // Lead to TeV error !
605  return false;
606  }
607  G4double vP = theMomentum.rho(); // Momentum of the decaying hadron
608  G4double dE = theMomentum.e(); // Energy of the decaying hadron
609  if(dE<vP)
610  {
611  G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
612  G4double accuracy=.000001*vP;
613  G4double emodif=std::fabs(dE-vP);
614  if(emodif<accuracy)
615  {
616  G4cerr<<"G4QHadron::DecayIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
617  theMomentum.setE(vP+emodif+.01*accuracy);
618  }
619  }
620  G4ThreeVector ltb = theMomentum.boostVector(); // Boost vector for backward Lor.Trans.
621 #ifdef debug
622  G4cout<<"G4QHadron::DecIn2:LorTrans v="<<ltb<<",f4Mom="<<f4Mom<<",s4Mom="<<s4Mom<<G4endl;
623 #endif
624  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* f4M="<<f4Mom<<G4endl;
625  f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS
626  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* s4M="<<s4Mom<<G4endl;
627  s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS
628 #ifdef debug
629  G4cout<<"G4QHadron::DecayIn2: ROOT OUTPUT f4Mom="<<f4Mom<<", s4Mom="<<s4Mom<<G4endl;
630 #endif
631  return true;
632 } // End of "DecayIn2"
633 
634 // Correction for the Hadron + fr decay in case of the new corM mass of the Hadron
635 G4bool G4QHadron::CorMDecayIn2(G4double corM, G4LorentzVector& fr4Mom)
636 {
637  G4double fM = fr4Mom.m(); // Mass of the Fragment
638  G4LorentzVector comp=theMomentum+fr4Mom; // 4Mom of the decaying compound system
639  G4double iM = comp.m(); // mass of the decaying compound system
640 #ifdef debug
641  G4cout<<"G4QH::CMDIn2: iM="<<iM<<comp<<"=>fM="<<fM<<"+corM="<<corM<<"="<<fM+corM<<G4endl;
642 #endif
643  G4double dE=iM-fM-corM;
644  //@@ Later on make a quark content check for the decay
645  if (fabs(dE)<.001)
646  {
647  G4double fR=fM/iM;
648  G4double cR=corM/iM;
649  fr4Mom=fR*comp;
650  theMomentum=cR*comp;
651  return true;
652  }
653  else if (dE<-.001 || iM==0.)
654  {
655  G4cerr<<"***G4QH::CorMDIn2***fM="<<fM<<" + cM="<<corM<<" > iM="<<iM<<",d="<<dE<<G4endl;
656  return false;
657  }
658  G4double corM2= corM*corM;
659  G4double fM2 = fM*fM;
660  G4double iM2 = iM*iM;
661  G4double d2 = iM2-fM2-corM2;
662  G4double p2 = (d2*d2/4.-fM2*corM2)/iM2; // Decay momentum(^2) in CMS of Quasmon
663  if (p2<0.)
664  {
665 #ifdef debug
666  G4cerr<<"**G4QH::CMDI2:p2="<<p2<<"<0,d="<<d2*d2/4.<<"<4*fM2*hM2="<<4*fM2*corM2<<G4endl;
667 #endif
668  p2=0.;
669  }
670  G4double p = sqrt(p2);
671  if(comp.e()<comp.rho())G4cerr<<"*G4QH::CorMDecayIn2:*Boost* comp4M="<<comp<<",e-p="
672  <<comp.e()-comp.rho()<<G4endl;
673  G4ThreeVector ltb = comp.boostVector(); // Boost vector for backward Lor.Trans.
674  G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans.
675  G4LorentzVector cm4Mom=fr4Mom; // Copy of fragment 4Mom to transform to CMS
676  if(cm4Mom.e()<cm4Mom.rho())
677  {
678  G4cerr<<"*G4QH::CorMDecIn2:*Boost* c4M="<<cm4Mom<<G4endl;
679  //cm4Mom.setE(1.0000001*cm4Mom.rho());
680  return false;
681  }
682  cm4Mom.boost(ltf); // Now it is in CMS (Forward Lor.Trans.)
683  G4double pfx= cm4Mom.px();
684  G4double pfy= cm4Mom.py();
685  G4double pfz= cm4Mom.pz();
686  G4double pt2= pfx*pfx+pfy*pfy;
687  G4double tx=0.;
688  G4double ty=0.;
689  if(pt2<=0.)
690  {
691  G4double phi= 360.*deg*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?)
692  tx=sin(phi);
693  ty=cos(phi);
694  }
695  else
696  {
697  G4double pt=sqrt(pt2);
698  tx=pfx/pt;
699  ty=pfy/pt;
700  }
701  G4double pc2=pt2+pfz*pfz;
702  G4double ct=0.;
703  if(pc2<=0.)
704  {
705  G4double rnd= G4UniformRand();
706  ct=1.-rnd-rnd;
707  }
708  else
709  {
710  G4double pc_value=sqrt(pc2);
711  ct=pfz/pc_value;
712  }
713 #ifdef debug
714  G4cout<<"G4QHadron::CorMDecayIn2: ct="<<ct<<", p="<<p<<G4endl;
715 #endif
716  G4double ps = p * sqrt(1.-ct*ct);
717  G4ThreeVector pVect(ps*tx,ps*ty,p*ct);
718  fr4Mom.setVect(pVect);
719  fr4Mom.setE(sqrt(fM2+p2));
720  theMomentum.setVect((-1)*pVect);
721  theMomentum.setE(sqrt(corM2+p2));
722 #ifdef debug
723  G4LorentzVector dif2=comp-fr4Mom-theMomentum;
724  G4cout<<"G4QH::CorMDIn2:c="<<comp<<"-f="<<fr4Mom<<"-4M="<<theMomentum<<"="<<dif2<<G4endl;
725 #endif
726  if(fr4Mom.e()+.001<fr4Mom.rho())G4cerr<<"*G4QH::CorMDecIn2:*Boost*fr4M="<<fr4Mom<<G4endl;
727  fr4Mom.boost(ltb); // Lor.Trans. of the Fragment back to LS
728  if(theMomentum.e()<theMomentum.rho())
729  {
730  G4cerr<<"*G4QH::CMDI2:4="<<theMomentum<<G4endl;
731  theMomentum.setE(1.0000001*theMomentum.rho());
732  }
733  theMomentum.boost(ltb); // Lor.Trans. of the Hadron back to LS
734 #ifdef debug
735  G4LorentzVector dif3=comp-fr4Mom-theMomentum;
736  G4cout<<"G4QH::CorMDecIn2:OUTPUT:f4M="<<fr4Mom<<",h4M="<<theMomentum<<"d="<<dif3<<G4endl;
737 #endif
738  return true;
739 } // End of "CorMDecayIn2"
740 
741 
742 // Fragment fr4Mom louse energy corE and transfer it to This Hadron
743 G4bool G4QHadron::CorEDecayIn2(G4double corE, G4LorentzVector& fr4Mom)
744 {
745  G4double fE = fr4Mom.m(); // Energy of the Fragment
746 #ifdef debug
747  G4cout<<"G4QH::CorEDecIn2:fE="<<fE<<fr4Mom<<">corE="<<corE<<",h4M="<<theMomentum<<G4endl;
748 #endif
749  if (fE+.001<=corE)
750  {
751 #ifdef debug
752  G4cerr<<"***G4QHadron::CorEDecIn2*** fE="<<fE<<"<corE="<<corE<<", d="<<corE-fE<<G4endl;
753 #endif
754  return false;
755  }
756  G4double fM2=fr4Mom.m2(); // Squared Mass of the Fragment
757  if(fM2<0.) fM2=0.;
758  G4double iPx=fr4Mom.px(); // Initial Px of the Fragment
759  G4double iPy=fr4Mom.py(); // Initial Py of the Fragment
760  G4double iPz=fr4Mom.pz(); // Initial Pz of the Fragment
761  G4double fP2=iPx*iPx+iPy*iPy+iPz*iPz; // Initial Squared 3-momentum of the Fragment
762  G4double finE = fE - corE; // Final energy of the fragment
763  G4double rP = sqrt((finE*finE-fM2)/fP2); // Reduction factor for the momentum
764  G4double fPx=iPx*rP;
765  G4double fPy=iPy*rP;
766  G4double fPz=iPz*rP;
767  fr4Mom= G4LorentzVector(fPx,fPy,fPz,finE);
768  G4double Px=theMomentum.px()+iPx-fPx;
769  G4double Py=theMomentum.py()+iPy-fPy;
770  G4double Pz=theMomentum.pz()+iPz-fPz;
771  G4double mE=theMomentum.e();
773  theMomentum= G4LorentzVector(Px,Py,Pz,mE+corE);
774 #ifdef debug
775  G4double difF=fr4Mom.m2()-fM2;
776  G4cout<<"G4QH::CorEDecIn2: dF="<<difF<<",out:"<<theMomentum<<fr4Mom<<G4endl;
777 #endif
778  return true;
779 } // End of "CorEDecayIn2"
780 
781 // Decay of the hadron in 3 particles i=>r+s+t
782 G4bool G4QHadron::DecayIn3
783  (G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, G4LorentzVector& t4Mom)
784 {
785 #ifdef debug
786  G4cout<<"G4QH::DIn3:"<<theMomentum<<"=>pf="<<f4Mom<<"+ps="<<s4Mom<<"+pt="<<t4Mom<<G4endl;
787 #endif
788  G4double iM = theMomentum.m(); // Mass of the decaying hadron
789  G4double fM = f4Mom.m(); // Mass of the 1st hadron
790  G4double sM = s4Mom.m(); // Mass of the 2nd hadron
791  G4double tM = t4Mom.m(); // Mass of the 3rd hadron
792  G4double eps = 0.001; // Accuracy of the split condition
793  if (fabs(iM-fM-sM-tM)<=eps)
794  {
795  G4double fR=fM/iM;
796  G4double sR=sM/iM;
797  G4double tR=tM/iM;
798  f4Mom=fR*theMomentum;
799  s4Mom=sR*theMomentum;
800  t4Mom=tR*theMomentum;
801  return true;
802  }
803  if (iM+eps<fM+sM+tM)
804  {
805  G4cout<<"***G4QHadron::DecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
806  <<",d="<<iM-fM-sM-tM<<G4endl;
807  return false;
808  }
809  G4double fM2 = fM*fM;
810  G4double sM2 = sM*sM;
811  G4double tM2 = tM*tM;
812  G4double iM2 = iM*iM;
813  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
814  G4double m12sMin =(fM+sM)*(fM+sM);
815  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
816  G4double rR = 0.;
817  G4double rnd= 1.;
818 #ifdef debug
819  G4int tr = 0; //@@ Comment if "cout" below is skiped @@
820 #endif
821  G4double m12s = 0.; // Fake definition before the Loop
822  while (rnd > rR)
823  {
824  m12s = m12sMin + m12sBase*G4UniformRand();
825  G4double e1=m12s+fM2-sM2;
826  G4double e2=iM2-m12s-tM2;
827  G4double four12=4*m12s;
828  G4double m13sRange=0.;
829  G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
830  if(dif<0.)
831  {
832 #ifdef debug
833  if(dif<-.01) G4cerr<<"*G4QHadron::DecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
834  <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
835 #endif
836  }
837  else m13sRange=sqrt(dif)/m12s;
838  rR = m13sRange/m13sBase;
839  rnd= G4UniformRand();
840 #ifdef debug
841  G4cout<<"G4QHadron::DecayIn3: try to decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
842 #endif
843  }
844  G4double m12 = sqrt(m12s); // Mass of the H1+H2 system
845  G4LorentzVector dh4Mom(0.,0.,0.,m12);
846 
847  if(!DecayIn2(t4Mom,dh4Mom))
848  {
849  G4cerr<<"***G4QHadron::DecayIn3: Exception1"<<G4endl;
850  //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
851  return false;
852  }
853 #ifdef debug
854  G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
855 #endif
856  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
857  {
858  G4cerr<<"***G4QHadron::DecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
859  //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
860  return false;
861  }
862  return true;
863 } // End of DecayIn3
864 
865 // Relative Decay of the hadron in 3 particles i=>f+s+t (t is with respect to minC<ct<maxC)
866 G4bool G4QHadron::RelDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
867  G4LorentzVector& t4Mom, G4LorentzVector& dir,
868  G4double maxCost, G4double minCost)
869 {
870 #ifdef debug
871  G4cout<<"G4QH::RelDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
872 #endif
873  G4double iM = theMomentum.m(); // Mass of the decaying hadron
874  G4double fM = f4Mom.m(); // Mass of the 1st hadron
875  G4double sM = s4Mom.m(); // Mass of the 2nd hadron
876  G4double tM = t4Mom.m(); // Mass of the 3rd hadron
877  G4double eps = 0.001; // Accuracy of the split condition
878  if (fabs(iM-fM-sM-tM)<=eps)
879  {
880  G4double fR=fM/iM;
881  G4double sR=sM/iM;
882  G4double tR=tM/iM;
883  f4Mom=fR*theMomentum;
884  s4Mom=sR*theMomentum;
885  t4Mom=tR*theMomentum;
886  return true;
887  }
888  if (iM+eps<fM+sM+tM)
889  {
890  G4cout<<"***G4QHadron::RelDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
891  <<",d="<<iM-fM-sM-tM<<G4endl;
892  return false;
893  }
894  G4double fM2 = fM*fM;
895  G4double sM2 = sM*sM;
896  G4double tM2 = tM*tM;
897  G4double iM2 = iM*iM;
898  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
899  G4double m12sMin =(fM+sM)*(fM+sM);
900  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
901  G4double rR = 0.;
902  G4double rnd= 1.;
903 #ifdef debug
904  G4int tr = 0; //@@ Comment if "cout" below is skiped @@
905 #endif
906  G4double m12s = 0.; // Fake definition before the Loop
907  while (rnd > rR)
908  {
909  m12s = m12sMin + m12sBase*G4UniformRand();
910  G4double e1=m12s+fM2-sM2;
911  G4double e2=iM2-m12s-tM2;
912  G4double four12=4*m12s;
913  G4double m13sRange=0.;
914  G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
915  if(dif<0.)
916  {
917 #ifdef debug
918  if(dif<-.01) G4cerr<<"G4QHadron::RelDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
919  <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
920 #endif
921  }
922  else m13sRange=sqrt(dif)/m12s;
923  rR = m13sRange/m13sBase;
924  rnd= G4UniformRand();
925 #ifdef debug
926  G4cout<<"G4QHadron::RelDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
927 #endif
928  }
929  G4double m12 = sqrt(m12s); // Mass of the H1+H2 system
930  G4LorentzVector dh4Mom(0.,0.,0.,m12);
931 
932  if(!RelDecayIn2(t4Mom,dh4Mom,dir,maxCost,minCost))
933  {
934  G4cerr<<"***G4QHadron::RelDecayIn3: Exception1"<<G4endl;
935  //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
936  return false;
937  }
938 #ifdef debug
939  G4cout<<"G4QHadron::RelDecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
940 #endif
941  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
942  {
943  G4cerr<<"***G4QHadron::RelDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
944  //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
945  return false;
946  }
947  return true;
948 } // End of RelDecayIn3
949 
950 // Relative Decay of hadron in 3: i=>f+s+t. dN/dO [cp>0:~cost^cp, cp<0:~(1-cost)^(-cp)]
951 G4bool G4QHadron::CopDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
952  G4LorentzVector& t4Mom, G4LorentzVector& dir, G4double cosp)
953 {
954 #ifdef debug
955  G4cout<<"G4QH::CopDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
956 #endif
957  G4double iM = theMomentum.m(); // Mass of the decaying hadron
958  G4double fM = f4Mom.m(); // Mass of the 1st hadron
959  G4double sM = s4Mom.m(); // Mass of the 2nd hadron
960  G4double tM = t4Mom.m(); // Mass of the 3rd hadron
961  G4double eps = 0.001; // Accuracy of the split condition
962  if (fabs(iM-fM-sM-tM)<=eps)
963  {
964  G4double fR=fM/iM;
965  G4double sR=sM/iM;
966  G4double tR=tM/iM;
967  f4Mom=fR*theMomentum;
968  s4Mom=sR*theMomentum;
969  t4Mom=tR*theMomentum;
970  return true;
971  }
972  if (iM+eps<fM+sM+tM)
973  {
974  G4cout<<"***G4QHadron::CopDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
975  <<",d="<<iM-fM-sM-tM<<G4endl;
976  return false;
977  }
978  G4double fM2 = fM*fM;
979  G4double sM2 = sM*sM;
980  G4double tM2 = tM*tM;
981  G4double iM2 = iM*iM;
982  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
983  G4double m12sMin =(fM+sM)*(fM+sM);
984  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
985  G4double rR = 0.;
986  G4double rnd= 1.;
987 #ifdef debug
988  G4int tr = 0; //@@ Comment if "cout" below is skiped @@
989 #endif
990  G4double m12s = 0.; // Fake definition before the Loop
991  while (rnd > rR)
992  {
993  m12s = m12sMin + m12sBase*G4UniformRand();
994  G4double e1=m12s+fM2-sM2;
995  G4double e2=iM2-m12s-tM2;
996  G4double four12=4*m12s;
997  G4double m13sRange=0.;
998  G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
999  if(dif<0.)
1000  {
1001 #ifdef debug
1002  if(dif<-.01) G4cerr<<"G4QHadron::CopDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
1003  <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
1004 #endif
1005  }
1006  else m13sRange=sqrt(dif)/m12s;
1007  rR = m13sRange/m13sBase;
1008  rnd= G4UniformRand();
1009 #ifdef debug
1010  G4cout<<"G4QHadron::CopDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
1011 #endif
1012  }
1013  G4double m12 = sqrt(m12s); // Mass of the H1+H2 system
1014  G4LorentzVector dh4Mom(0.,0.,0.,m12);
1015 
1016  if(!CopDecayIn2(t4Mom,dh4Mom,dir,cosp))
1017  {
1018  G4cerr<<"***G4QHadron::CopDecayIn3: Exception1"<<G4endl;
1019  //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
1020  return false;
1021  }
1022 #ifdef debug
1023  G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
1024 #endif
1025  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
1026  {
1027  G4cerr<<"***G4QHadron::CopDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
1028  //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
1029  return false;
1030  }
1031  return true;
1032 } // End of CopDecayIn3
1033 
1034 // Randomize particle mass taking into account the width
1035 G4double G4QHadron::RandomizeMass(G4QParticle* pPart, G4double maxM)
1036 {
1037  G4double meanM = theQPDG.GetMass();
1038  G4double width = theQPDG.GetWidth()/2.;
1039 #ifdef debug
1040  G4cout<<"G4QHadron::RandomizeMass: meanM="<<meanM<<", halfWidth="<<width<<G4endl;
1041 #endif
1042  if(maxM<meanM-3*width)
1043  {
1044 #ifdef debug
1045  G4cout<<"***G4QH::RandM:m=0 maxM="<<maxM<<"<meanM="<<meanM<<"-3*halfW="<<width<<G4endl;
1046 #endif
1047  return 0.;
1048  }
1050  if(width==0.)
1051  {
1052 #ifdef debug
1053  if(meanM>maxM) G4cerr<<"***G4QHadron::RandM:Stable m="<<meanM<<">maxM="<<maxM<<G4endl;
1054 #endif
1055  return meanM;
1056  //return 0.;
1057  }
1058  else if(width<0.)
1059  {
1060  // G4cerr<<"***G4QHadron::RandM: width="<<width<<"<0,PDGC="<<theQPDG.GetPDGCode()<<G4endl;
1061  // throw G4QException("G4QHadron::RandomizeMass: with the width of the Hadron < 0.");
1062  G4ExceptionDescription ed;
1063  ed << "width of the Hadron < 0. : width=" << width << "<0,PDGC="
1064  << theQPDG.GetPDGCode() << G4endl;
1065  G4Exception("G4QHadron::RandomizeMass()", "HAD_CHPS_0000", FatalException, ed);
1066  }
1067  G4double minM = pPart->MinMassOfFragm();
1068  if(minM>maxM)
1069  {
1070 #ifdef debug
1071  G4cout<<"***G4QHadron::RandomizeMass:for PDG="<<theQPDG.GetPDGCode()<<" minM="<<minM
1072  <<" > maxM="<<maxM<<G4endl;
1073 #endif
1074  return 0.;
1075  }
1076  //Now calculate the Breit-Wigner distribution with two cuts
1077  G4double v1=atan((minM-meanM)/width);
1078  G4double v2=atan((maxM-meanM)/width);
1079  G4double dv=v2-v1;
1080 #ifdef debug
1081  G4cout<<"G4QHadr::RandM:Mi="<<minM<<",i="<<v1<<",Ma="<<maxM<<",a="<<v2<<","<<dv<<G4endl;
1082 #endif
1083  return meanM+width*tan(v1+dv*G4UniformRand());
1084 }
1085 
1086 // Split the hadron in two partons ( quark = anti-diquark v.s. anti-quark = diquark)
1087 void G4QHadron::SplitUp()
1088 {
1089  if (isSplit) return;
1090 #ifdef pdebug
1091  G4cout<<"G4QHadron::SplitUp ***IsCalled***, before Splitting nC="<<Color.size()
1092  <<", SoftColCount="<<theCollisionCount<<G4endl;
1093 #endif
1094  isSplit=true; // PutUp isSplit flag to avoid remake
1095  if (!Color.empty()) return; // Do not split if it is already split
1096  if (!theCollisionCount) // Diffractive splitting from particDef
1097  {
1098  G4QParton* Left = 0;
1099  G4QParton* Right = 0;
1100  GetValenceQuarkFlavors(Left, Right);
1101  G4ThreeVector Pos=GetPosition();
1102  Left->SetPosition(Pos);
1103  Right->SetPosition(Pos);
1104 
1105  G4double theMomPlus = theMomentum.plus(); // E+pz
1106  G4double theMomMinus = theMomentum.minus(); // E-pz
1107 #ifdef pdebug
1108  G4cout<<"G4QHadron::SplitUp: *Dif* possition="<<Pos<<", 4M="<<theMomentum<<G4endl;
1109 #endif
1110 
1111  // momenta of string ends
1112  G4double pt2 = theMomentum.perp2();
1113  G4double transverseMass2 = theMomPlus*theMomMinus;
1114  if(transverseMass2<0.) transverseMass2=0.;
1115  G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
1116  G4ThreeVector pt(0., 0., 0.); // Prototype
1117  if(maxAvailMomentum2/widthOfPtSquare > 0.01)
1118  pt=GaussianPt(widthOfPtSquare, maxAvailMomentum2);
1119 #ifdef pdebug
1120  G4cout<<"G4QHadron::SplitUp: *Dif* maxMom2="<<maxAvailMomentum2<<", pt="<<pt<<G4endl;
1121 #endif
1122  G4LorentzVector LeftMom(pt, 0.);
1123  G4LorentzVector RightMom;
1124  RightMom.setPx(theMomentum.px() - pt.x());
1125  RightMom.setPy(theMomentum.py() - pt.y());
1126 #ifdef pdebug
1127  G4cout<<"G4QHadron::SplitUp: *Dif* right4m="<<RightMom<<", left4M="<<LeftMom<<G4endl;
1128 #endif
1129  G4double Local1 = theMomMinus + (RightMom.perp2() - LeftMom.perp2()) / theMomPlus;
1130  G4double Local2 = std::sqrt(std::max(0., Local1*Local1 -
1131  4*RightMom.perp2()*theMomMinus / theMomPlus));
1132 #ifdef pdebug
1133  G4cout<<"G4QHadron::SplitUp:Dif,L1="<<Local1<<",L2="<<Local2<<",D="<<Direction<<G4endl;
1134 #endif
1135  if (Direction) Local2 = -Local2;
1136  G4double RightMinus = 0.5*(Local1 + Local2);
1137  G4double LeftMinus = theMomentum.minus() - RightMinus;
1138 #ifdef pdebug
1139  G4cout<<"G4QHadron::SplitUp: *Dif* Rminus="<<RightMinus<<",Lminus="<<LeftMinus<<",hmm="
1140  <<theMomentum.minus()<<G4endl;
1141 #endif
1142  G4double LeftPlus = 0.;
1143  if(LeftMinus) LeftPlus = LeftMom.perp2()/LeftMinus;
1144  G4double RightPlus = theMomentum.plus() - LeftPlus;
1145 #ifdef pdebug
1146  G4cout<<"G4QHadron::SplitUp: *Dif* Rplus="<<RightPlus<<", Lplus="<<LeftPlus<<G4endl;
1147 #endif
1148  LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
1149  LeftMom.setE (0.5*(LeftPlus + LeftMinus));
1150  RightMom.setPz(0.5*(RightPlus - RightMinus));
1151  RightMom.setE (0.5*(RightPlus + RightMinus));
1152  //G4cout<<"DSU 6: Left4M="<<LeftMom<<", Right4M="<<RightMom<<G4endl;
1153 #ifdef pdebug
1154  G4cout<<"G4QHadron::SplitUp: *Dif* -final- R4m="<<RightMom<<", L4M="<<LeftMom<<", L+R="
1155  <<RightMom+LeftMom<<", D4M="<<theMomentum-RightMom+LeftMom<<G4endl;
1156 #endif
1157  Left->Set4Momentum(LeftMom);
1158  Right->Set4Momentum(RightMom);
1159  Color.push_back(Left);
1160  AntiColor.push_back(Right);
1161  }
1162  else
1163  {
1164  // Soft hadronization splitting: sample transversal momenta for sea and valence quarks
1165  //G4double phi, pts;
1166  G4ThreeVector SumP(0.,0.,0.); // Remember the hadron position
1167  G4ThreeVector Pos = GetPosition(); // Remember the hadron position
1168  G4int nSeaPair = theCollisionCount-1; // a#of sea-pairs
1169 #ifdef pdebug
1170  G4cout<<"G4QHadron::SplitUp:*Soft* Pos="<<Pos<<", nSeaPair="<<nSeaPair<<G4endl;
1171 #endif
1172  // here the condition,to ensure viability of splitting, also in cases
1173  // where difractive excitation occured together with soft scattering.
1174  for (G4int aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++) // If the sea pairs exist!
1175  {
1176  // choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
1177  G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress);
1178 
1179  // BuildSeaQuark() determines quark spin, isospin and colour
1180  // via parton-constructor G4QParton(aPDGCode)
1181  G4QParton* aParton = BuildSeaQuark(false, aPDGCode); // quark/anti-diquark creation
1182 
1183  // save colour a spin-3 for anti-quark
1184  G4int firstPartonColour = aParton->GetColour();
1185  G4double firstPartonSpinZ = aParton->GetSpinZ();
1186 #ifdef pdebug
1187  G4cout<<"G4QHadron::SplitUp:*Soft* Part1 PDG="<<aPDGCode<<", Col="<<firstPartonColour
1188  <<", SpinZ="<<firstPartonSpinZ<<", 4M="<<aParton->Get4Momentum()<<G4endl;
1189 #endif
1190  SumP+=aParton->Get4Momentum();
1191  Color.push_back(aParton); // Quark/anti-diquark is filled
1192 
1193  // create anti-quark
1194  aParton = BuildSeaQuark(true, aPDGCode); // Redefine "aParton" (working pointer)
1195  aParton->SetSpinZ(-firstPartonSpinZ);
1196  aParton->SetColour(-firstPartonColour);
1197 #ifdef pdebug
1198  G4cout<<"G4QHadron::SplUp:Sft,P2="<<aParton->Get4Momentum()<<",i="<<aSeaPair<<G4endl;
1199 #endif
1200 
1201  SumP+=aParton->Get4Momentum();
1202  AntiColor.push_back(aParton); // Anti-quark/diquark is filled
1203 #ifdef pdebug
1204  G4cout<<"G4QHadron::SplUp:*Sft* Antiquark is filled, i="<<aSeaPair<<G4endl;
1205 #endif
1206  }
1207  // ---- Create valence quarks/diquarks
1208  G4QParton* pColorParton = 0;
1209  G4QParton* pAntiColorParton = 0;
1210  GetValenceQuarkFlavors(pColorParton, pAntiColorParton);
1211  G4int ColorEncoding = pColorParton->GetPDGCode();
1212 #ifdef pdebug
1213  G4int AntiColorEncoding = pAntiColorParton->GetPDGCode();
1214  G4cout<<"G4QHadron::SplUp:*Sft*,C="<<ColorEncoding<<", AC="<<AntiColorEncoding<<G4endl;
1215 #endif
1216  G4ThreeVector ptr = GaussianPt(sigmaPt, DBL_MAX);
1217  SumP += ptr;
1218 #ifdef pdebug
1219  G4cout<<"G4QHadron::SplitUp: *Sft*, ptr="<<ptr<<G4endl;
1220 #endif
1221 
1222  if (ColorEncoding < 0) // use particle definition
1223  {
1224  G4LorentzVector ColorMom(-SumP, 0);
1225  pColorParton->Set4Momentum(ColorMom);
1226  G4LorentzVector AntiColorMom(ptr, 0.);
1227  pAntiColorParton->Set4Momentum(AntiColorMom);
1228  }
1229  else
1230  {
1231  G4LorentzVector ColorMom(ptr, 0);
1232  pColorParton->Set4Momentum(ColorMom);
1233  G4LorentzVector AntiColorMom(-SumP, 0);
1234  pAntiColorParton->Set4Momentum(AntiColorMom);
1235  }
1236  Color.push_back(pColorParton);
1237  AntiColor.push_back(pAntiColorParton);
1238 #ifdef pdebug
1239  G4cout<<"G4QHadron::SplitUp: *Soft* Col&Anticol are filled PDG="<<GetPDGCode()<<G4endl;
1240 #endif
1241  // Sample X
1242  G4int nColor=Color.size();
1243  G4int nAntiColor=AntiColor.size();
1244  if(nColor!=nAntiColor || nColor != nSeaPair+1)
1245  {
1246  G4cerr<<"***G4QHadron::SplitUp: nA="<<nAntiColor<<",nAC="<<nColor<<",nSea="<<nSeaPair
1247  <<G4endl;
1248  G4Exception("G4QHadron::SplitUp:","72",FatalException,"Colours&AntiColours notSinc");
1249  }
1250 #ifdef pdebug
1251  G4cout<<"G4QHad::SpUp:,nPartons="<<nColor+nColor<<G4endl;
1252 #endif
1253  G4int dnCol=nColor+nColor;
1254  // From here two algorithm of splitting can be used (All(default): New, OBO: Olg, Bad)
1255  G4double* xs=RandomX(dnCol); // All-Non-iterative CHIPS algorithm of splitting
1256  // Instead one can try one-by-one CHIPS algorithm (faster? but not exact). OBO comment.
1257  //G4double Xmax=1.; // OBO
1258 #ifdef pdebug
1259  G4cout<<"G4QHadron::SplitUp:*Sft* Loop ColorX="<<ColorX<<G4endl;
1260 #endif
1261  std::list<G4QParton*>::iterator icolor = Color.begin();
1262  std::list<G4QParton*>::iterator ecolor = Color.end();
1263  std::list<G4QParton*>::iterator ianticolor = AntiColor.begin();
1264  std::list<G4QParton*>::iterator eanticolor = AntiColor.end();
1265  G4int xi=-1; // XIndex for All-Non-interactive CHIPS algorithm
1266  //G4double X=0.; // OBO
1267  for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor)
1268  {
1269  (*icolor)->SetX(xs[++xi]); // All-Non-iterative CHIPS algorithm of splitting
1270  //X=SampleCHIPSX(Xmax, dnCol); // OBO
1271  //Xmax-=X; // OBO
1272  //--dnCol; // OBO
1273  //(*icolor)->SetX(X); // OBO
1274  // ----
1275  (*icolor)->DefineEPz(theMomentum);
1276  (*ianticolor)->SetX(xs[++xi]); // All-Non-iterative CHIPS algorithm of splitting
1277  //X=SampleCHIPSX(Xmax, dnCol); // OBO
1278  //Xmax-=X; // OBO
1279  //--dnCol; // OBO
1280  //(*ianticolor)->SetX(X); // OBO
1281  // ----
1282  (*ianticolor)->DefineEPz(theMomentum);
1283  }
1284  delete[] xs; // The calculated array must be deleted (All)
1285 #ifdef pdebug
1286  G4cout<<"G4QHadron::SplitUp: *Soft* ===> End, ColSize="<<Color.size()<<G4endl;
1287 #endif
1288  return;
1289  }
1290 } // End of SplitUp
1291 
1292 // Boost hadron 4-momentum, using Boost Lorentz vector
1293 void G4QHadron::Boost(const G4LorentzVector& boost4M)
1294 {
1295  // see CERNLIB short writeup U101 for the algorithm
1296  G4double bm=boost4M.mag();
1297  G4double factor=(theMomentum.vect()*boost4M.vect()/(boost4M.e()+bm)-theMomentum.e())/bm;
1298  theMomentum.setE(theMomentum.dot(boost4M)/bm);
1299  theMomentum.setVect(factor*boost4M.vect() + theMomentum.vect());
1300 } // End of Boost
1301 
1302 // Build one (?M.K.) sea-quark
1303 G4QParton* G4QHadron::BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode)
1304 {
1305  if (isAntiQuark) aPDGCode*=-1;
1306  G4QParton* result = new G4QParton(aPDGCode);
1307  result->SetPosition(GetPosition());
1308  G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
1309  G4LorentzVector a4Momentum(aPtVector, 0);
1310  result->Set4Momentum(a4Momentum);
1311  return result;
1312 } // End of BuildSeaQuark
1313 
1314 // Fast non-iterative CHIPS algorithm
1315 G4double* G4QHadron::RandomX(G4int nPart)
1316 {
1317  G4double* x = 0;
1318  if(nPart<2)
1319  {
1320  G4cout<<"-Warning-G4QHadron::RandomX: nPart="<<nPart<<" < 2"<<G4endl;
1321  return x;
1322  }
1323  x = new G4double[nPart];
1324  G4int nP1=nPart-1;
1325  x[0]=G4UniformRand();
1326  for(G4int i=1; i<nP1; ++i)
1327  {
1328  G4double r=G4UniformRand();
1329  G4int j=0;
1330  for( ; j<i; ++j) if(r < x[j])
1331  {
1332  for(G4int k=i; k>j; --k) x[k]=x[k-1];
1333  x[j]=r;
1334  break;
1335  }
1336  if(j==i) x[i]=r;
1337  }
1338  x[nP1]=1.;
1339  for(G4int i=nP1; i>0; --i) x[i]-=x[i-1];
1340  return x;
1341 }
1342 
1343 // Non-iterative recursive phase-space CHIPS algorthm
1344 G4double G4QHadron::SampleCHIPSX(G4double anXtot, G4int nSea)
1345 {
1346  G4double ns_value=nSea;
1347  if(nSea<1 || anXtot<=0.) G4cout<<"-Warning-G4QHad::SCX:N="<<nSea<<",tX="<<anXtot<<G4endl;
1348  if(nSea<2) return anXtot;
1349  return anXtot*(1.-std::pow(G4UniformRand(),1./ns_value));
1350 }
1351 
1352 // Get flavors for the valence quarks of this hadron
1353 void G4QHadron::GetValenceQuarkFlavors(G4QParton* &Parton1, G4QParton* &Parton2)
1354 {
1355  // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
1356  G4int aEnd=0;
1357  G4int bEnd=0;
1358  G4int HadronEncoding = GetPDGCode();
1359  if(!(GetBaryonNumber())) SplitMeson(HadronEncoding, &aEnd, &bEnd);
1360  else SplitBaryon(HadronEncoding, &aEnd, &bEnd);
1361 
1362  Parton1 = new G4QParton(aEnd);
1363  Parton1->SetPosition(GetPosition());
1364 
1365 // G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
1366 // G4cerr << "Parton 1: "
1367 // << " PDGcode: " << aEnd
1368 // << " - Name: " << Parton1->GetDefinition()->GetParticleName()
1369 // << " - Type: " << Parton1->GetDefinition()->GetParticleType()
1370 // << " - Spin-3: " << Parton1->GetSpinZ()
1371 // << " - Colour: " << Parton1->GetColour() << G4endl;
1372 
1373  Parton2 = new G4QParton(bEnd);
1374  Parton2->SetPosition(GetPosition());
1375 
1376 // G4cerr << "Parton 2: "
1377 // << " PDGcode: " << bEnd
1378 // << " - Name: " << Parton2->GetDefinition()->GetParticleName()
1379 // << " - Type: " << Parton2->GetDefinition()->GetParticleType()
1380 // << " - Spin-3: " << Parton2->GetSpinZ()
1381 // << " - Colour: " << Parton2->GetColour() << G4endl;
1382 // G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;
1383 
1384  // colour of parton 1 choosen at random by G4QParton(aEnd)
1385  // colour of parton 2 is the opposite:
1386 
1387  Parton2->SetColour(-(Parton1->GetColour()));
1388 
1389  // isospin-3 of both partons is handled by G4Parton(PDGCode)
1390 
1391  // spin-3 of parton 1 and 2 choosen at random by G4QParton(aEnd)
1392  // spin-3 of parton 2 may be constrained by spin of original particle:
1393 
1394  if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > GetSpin())
1395  Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
1396 
1397 // G4cerr << "Parton 2: "
1398 // << " PDGcode: " << bEnd
1399 // << " - Name: " << Parton2->GetDefinition()->GetParticleName()
1400 // << " - Type: " << Parton2->GetDefinition()->GetParticleType()
1401 // << " - Spin-3: " << Parton2->GetSpinZ()
1402 // << " - Colour: " << Parton2->GetColour() << G4endl;
1403 // G4cerr << "------------" << G4endl;
1404 
1405 } // End of GetValenceQuarkFlavors
1406 
1407 G4bool G4QHadron::SplitMeson(G4int PDGcode, G4int* aEnd, G4int* bEnd)
1408 {
1409  G4bool result = true;
1410  G4int absPDGcode = std::abs(PDGcode);
1411  if (absPDGcode >= 1000) return false;
1412  if(absPDGcode == 22 || absPDGcode == 111) // only u-ubar, d-dbar configurations
1413  {
1414  G4int it=1;
1415  if(G4UniformRand()<.5) it++;
1416  *aEnd = it;
1417  *bEnd = -it;
1418  }
1419  else if(absPDGcode == 130 || absPDGcode == 310) // K0-K0bar mixing
1420  {
1421  G4int it=1;
1422  if(G4UniformRand()<.5) it=-1;
1423  *aEnd = it;
1424  if(it>0) *bEnd = -3;
1425  else *bEnd = 3;
1426  }
1427  else
1428  {
1429  G4int heavy = absPDGcode/100;
1430  G4int light = (absPDGcode%100)/10;
1431  G4int anti = 1 - 2*(std::max(heavy, light)%2);
1432  if (PDGcode < 0 ) anti = -anti;
1433  heavy *= anti;
1434  light *= -anti;
1435  if ( anti < 0) G4SwapObj(&heavy, &light);
1436  *aEnd = heavy;
1437  *bEnd = light;
1438  }
1439  return result;
1440 }
1441 
1442 G4bool G4QHadron::SplitBaryon(G4int PDGcode, G4int* quark, G4int* diQuark)
1443 {
1444  static const G4double r2=.5;
1445  static const G4double r3=1./3.;
1446  static const G4double d3=2./3.;
1447  static const G4double r4=1./4.;
1448  static const G4double r6=1./6.;
1449  static const G4double r12=1./12.;
1450  //
1451  std::pair<G4int,G4int> qdq[5];
1452  G4double prb[5];
1453  G4int nc=0;
1454  G4int aPDGcode=std::abs(PDGcode);
1455  if(aPDGcode==2212) // ==> Proton
1456  {
1457  nc=3;
1458  qdq[0]=make_pair(2203, 1); prb[0]=r3; // uu_1, d
1459  qdq[1]=make_pair(2103, 2); prb[1]=r6; // ud_1, u
1460  qdq[2]=make_pair(2101, 2); prb[2]=r2; // ud_0, u
1461  }
1462  else if(aPDGcode==2112) // ==> Neutron
1463  {
1464  nc=3;
1465  qdq[0]=make_pair(2103, 1); prb[0]=r6; // ud_1, d
1466  qdq[1]=make_pair(2101, 1); prb[1]=r2; // ud_0, d
1467  qdq[2]=make_pair(1103, 2); prb[2]=r3; // dd_1, u
1468  }
1469  else if(aPDGcode%10<3) // ==> Spin 1/2 Hyperons
1470  {
1471  if(aPDGcode==3122) // Lambda
1472  {
1473  nc=5;
1474  qdq[0]=make_pair(2103, 3); prb[0]=r3; // ud_1, s
1475  qdq[1]=make_pair(3203, 1); prb[1]=r4; // su_1, d
1476  qdq[2]=make_pair(3201, 1); prb[2]=r12; // su_0, d
1477  qdq[3]=make_pair(3103, 2); prb[3]=r4; // sd_1, u
1478  qdq[4]=make_pair(3101, 2); prb[4]=r12; // sd_0, u
1479  }
1480  else if(aPDGcode==3222) // Sigma+
1481  {
1482  nc=3;
1483  qdq[0]=make_pair(2203, 3); prb[0]=r3; // uu_1, s
1484  qdq[1]=make_pair(3203, 2); prb[1]=r6; // su_1, d
1485  qdq[2]=make_pair(3201, 2); prb[2]=r2; // su_0, d
1486  }
1487  else if(aPDGcode==3212) // Sigma0
1488  {
1489  nc=5;
1490  qdq[0]=make_pair(2103, 3); prb[0]=r3; // ud_1, s
1491  qdq[1]=make_pair(3203, 1); prb[1]=r12; // su_1, d
1492  qdq[2]=make_pair(3201, 1); prb[2]=r4; // su_0, d
1493  qdq[3]=make_pair(3103, 2); prb[3]=r12; // sd_1, u
1494  qdq[4]=make_pair(3101, 2); prb[4]=r4; // sd_0, u
1495  }
1496  else if(aPDGcode==3112) // Sigma-
1497  {
1498  nc=3;
1499  qdq[0]=make_pair(1103, 3); prb[0]=r3; // dd_1, s
1500  qdq[1]=make_pair(3103, 1); prb[1]=r6; // sd_1, d
1501  qdq[2]=make_pair(3101, 1); prb[2]=r2; // sd_0, d
1502  }
1503  else if(aPDGcode==3312) // Xi-
1504  {
1505  nc=3;
1506  qdq[0]=make_pair(3103, 3); prb[0]=r6; // sd_1, s
1507  qdq[1]=make_pair(3101, 3); prb[1]=r2; // sd_0, s
1508  qdq[2]=make_pair(3303, 1); prb[2]=r3; // ss_1, d
1509  }
1510  else if(aPDGcode==3322) // Xi0
1511  {
1512  nc=3;
1513  qdq[0]=make_pair(3203, 3); prb[0]=r6; // su_1, s
1514  qdq[1]=make_pair(3201, 3); prb[1]=r2; // su_0, s
1515  qdq[2]=make_pair(3303, 2); prb[2]=r3; // ss_1, u
1516  }
1517  else return false;
1518  }
1519  else // ==> Spin 3/2 Baryons (Spin>3/2 is ERROR)
1520  {
1521  if(aPDGcode==3334)
1522  {
1523  nc=1;
1524  qdq[0]=make_pair(3303, 3); prb[0]=1.; // ss_1, s
1525  }
1526  else if(aPDGcode==2224)
1527  {
1528  nc=1;
1529  qdq[0]=make_pair(2203, 2); prb[0]=1.; // uu_1, s
1530  }
1531  else if(aPDGcode==2214)
1532  {
1533  nc=2;
1534  qdq[0]=make_pair(2203, 1); prb[0]=r3; // uu_1, d
1535  qdq[1]=make_pair(2103, 2); prb[1]=d3; // ud_1, u
1536  }
1537  else if(aPDGcode==2114)
1538  {
1539  nc=2;
1540  qdq[0]=make_pair(1103, 2); prb[0]=r3; // dd_1, u
1541  qdq[1]=make_pair(2103, 1); prb[1]=d3; // ud_1, d
1542  }
1543  else if(aPDGcode==1114)
1544  {
1545  nc=1;
1546  qdq[0]=make_pair(1103, 1); prb[0]=1.; // uu_1, s
1547  }
1548  else if(aPDGcode==3224)
1549  {
1550  nc=2;
1551  qdq[0]=make_pair(2203, 3); prb[0]=r3; // uu_1, s
1552  qdq[1]=make_pair(3203, 2); prb[1]=d3; // su_1, u
1553  }
1554  else if(aPDGcode==3214) // @@ SU(3) is broken because of the s-quark mass
1555  {
1556  nc=3;
1557  qdq[0]=make_pair(2103, 3); prb[0]=r3; // ud_1, s
1558  qdq[1]=make_pair(3203, 1); prb[1]=r3; // su_1, d
1559  qdq[2]=make_pair(3103, 2); prb[2]=r3; // sd_1, u
1560  }
1561  else if(aPDGcode==3114)
1562  {
1563  nc=2;
1564  qdq[0]=make_pair(1103, 3); prb[0]=r3; // dd_1, s
1565  qdq[1]=make_pair(3103, 1); prb[1]=d3; // sd_1, d
1566  }
1567  else if(aPDGcode==3324)
1568  {
1569  nc=2;
1570  qdq[0]=make_pair(3203, 3); prb[0]=r3; // su_1, s
1571  qdq[1]=make_pair(3303, 2); prb[1]=d3; // ss_1, u
1572  }
1573  else if(aPDGcode==3314)
1574  {
1575  nc=2;
1576  qdq[0]=make_pair(3103, 3); prb[0]=d3; // sd_1, s
1577  qdq[1]=make_pair(3303, 1); prb[1]=r3; // ss_1, d
1578  }
1579  else return false;
1580  }
1581  G4double random = G4UniformRand();
1582  G4double sum = 0.;
1583  for(G4int i=0; i<nc; i++)
1584  {
1585  sum += prb[i];
1586  if(sum>random)
1587  {
1588  if(PDGcode>0)
1589  {
1590  *diQuark= qdq[i].first;
1591  *quark = qdq[i].second;
1592  }
1593  else
1594  {
1595  *diQuark= -qdq[i].second;
1596  *quark = -qdq[i].first;
1597  }
1598  break;
1599  }
1600  }
1601  return true;
1602 }
1603 
1604 // This is not usual Gaussian, in fact it is dN/d(pt) ~ pt * exp(-pt^2/pt0^2)
1605 G4ThreeVector G4QHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
1606 {
1607  G4double R=0.;
1608  while ((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare){}
1609  R = std::sqrt(R);
1610  G4double phi = twopi*G4UniformRand();
1611  return G4ThreeVector(R*std::cos(phi), R*std::sin(phi), 0.);
1612 }
1613 
1614 G4QParton* G4QHadron::GetNextParton()
1615 {
1616  if(Color.size()==0) return 0;
1617  G4QParton* result = Color.back();
1618  Color.pop_back();
1619  return result;
1620 }
1621 
1622 G4QParton* G4QHadron::GetNextAntiParton()
1623 {
1624  if(AntiColor.size() == 0) return 0;
1625  G4QParton* result = AntiColor.front();
1626  AntiColor.pop_front();
1627  return result;
1628 }
1629 
1630 // Random Split of the Hadron in 2 Partons (caller is responsible for G4QPartonPair delete)
1631 G4QPartonPair* G4QHadron::SplitInTwoPartons() // If result=0: impossible to split (?)
1632 {
1633  if(std::abs(GetBaryonNumber())>1) // Not Baryons or Mesons or Anti-Baryons
1634  {
1635  G4cerr<<"***G4QHadron::SplitInTwoPartons: Can not split QC="<<valQ<< G4endl;
1636  G4Exception("G4QFragmentation::ChooseX:","72",FatalException,"NeitherMesonNorBaryon");
1637  }
1638  std::pair<G4int,G4int> PP = valQ.MakePartonPair();
1639  return new G4QPartonPair(new G4QParton(PP.first), new G4QParton(PP.second));
1640 }
Generated on Sat May 25 2013 14:33:46 for Geant4 by   doxygen 1.8.4