source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QHadron.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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