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

Last change on this file since 1071 was 1055, checked in by garnier, 15 years ago

maj sur la beta de geant 4.9.3
File size: 56.7 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.54 2009/02/23 09:49:24 mkossov Exp $
28// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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 pdebug
42//#define sdebug
43//#define ppdebug
44
45#include "G4QHadron.hh"
46#include <cmath>
47using namespace std;
48
49G4double G4QHadron::alpha = -0.5;   // changing rapidity distribution for all
50G4double G4QHadron::beta  = 2.5;    // changing rapidity distribution for projectile region
51G4double G4QHadron::theMinPz = 70.*MeV;             // Can be from 14 to 140 MeV
52G4double G4QHadron::StrangeSuppress = 0.48;         // ? M.K.
53G4double G4QHadron::sigmaPt = 1.7*GeV;              // Can be 0
54G4double G4QHadron::widthOfPtSquare = 0.01*GeV*GeV; // ? M.K.
55G4double G4QHadron::minTransverseMass = 1.*keV;     // ? M.K.
56
57G4QHadron::G4QHadron() : theQPDG(0),theMomentum(0.,0.,0.,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
61G4QHadron::G4QHadron(G4LorentzVector p) : theQPDG(0),theMomentum(p),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
66G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p) : theQPDG(PDGCode),theMomentum(p),
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
86G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p) : theQPDG(QPDG),theMomentum(p),
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
104G4QHadron::G4QHadron(G4QContent QC, G4LorentzVector p): theQPDG(0),theMomentum(p),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
114G4QHadron::G4QHadron(G4int PDGCode, G4double aMass, G4QContent QC) :
115  theQPDG(PDGCode),theMomentum(0.,0.,0., aMass),valQ(QC),nFragm(0),thePosition(0.,0.,0.),
116  theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.),
117  formTime(0.)
118{}
119
120G4QHadron::G4QHadron(G4QPDGCode QPDG, G4double aMass, G4QContent QC) :
121  theQPDG(QPDG),theMomentum(0.,0.,0., aMass),valQ(QC),nFragm(0),thePosition(0.,0.,0.),
122  theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.),
123  formTime(0.)
124{}
125
126G4QHadron::G4QHadron(G4int PDGCode, G4LorentzVector p, G4QContent QC) :
127  theQPDG(PDGCode),theMomentum(p),valQ(QC),nFragm(0),thePosition(0.,0.,0.),
128  theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.),
129  formTime(0.)
130{}
131
132G4QHadron::G4QHadron(G4QPDGCode QPDG, G4LorentzVector p, G4QContent QC) :
133  theQPDG(QPDG),theMomentum(p),valQ(QC),nFragm(0),thePosition(0.,0.,0.),
134  theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.),
135  formTime(0.)
136{}
137
138G4QHadron::G4QHadron(G4QParticle* pPart, G4double maxM) :
139  theQPDG(pPart->GetQPDG()),theMomentum(0.,0.,0.,0.),nFragm(0),thePosition(0.,0.,0.),
140  theCollisionCount(0),isSplit(false),Direction(true),Color(),AntiColor(),bindE(0.),
141  formTime(0.)
142{
143#ifdef debug
144  G4cout<<"G4QHadron is created & randomized with maxM="<<maxM<<G4endl;
145#endif
146  G4int PDGCode = theQPDG.GetPDGCode();
147  if(PDGCode<2)G4cerr<<"***G4QHadron:(M) PDGC="<<PDGCode<<" use other constructor"<<G4endl;
148  valQ=theQPDG.GetQuarkContent();
149  theMomentum.setE(RandomizeMass(pPart, maxM));
150}
151
152G4QHadron::G4QHadron(const G4QHadron& right)
153{
154  theMomentum         = right.theMomentum;
155  theQPDG             = right.theQPDG;
156  valQ                = right.valQ;
157  nFragm              = right.nFragm;
158  thePosition         = right.thePosition;     
159  theCollisionCount   = 0;
160  isSplit             = false;
161  Direction           = right.Direction;
162  bindE               = right.bindE;
163  formTime            = right.formTime;
164}
165
166G4QHadron::G4QHadron(const G4QHadron* right)
167{
168  theMomentum         = right->theMomentum;
169  theQPDG             = right->theQPDG;
170  valQ                = right->valQ;
171  nFragm              = right->nFragm;
172  thePosition         = right->thePosition;     
173  theCollisionCount   = 0;
174  isSplit             = false;
175  Direction           = right->Direction;
176  bindE               = right->bindE;
177  formTime            = right->formTime;
178}
179
180G4QHadron::G4QHadron(const G4QHadron* right, G4int C, G4ThreeVector P, G4LorentzVector M)
181{
182  theMomentum         = M;
183  theQPDG             = right->theQPDG;
184  valQ                = right->valQ;
185  nFragm              = right->nFragm;
186  thePosition         = P;     
187  theCollisionCount   = C;
188  isSplit             = false;
189  Direction           = right->Direction;
190  bindE               = right->bindE;
191  formTime            = right->formTime;
192}
193
194const G4QHadron& G4QHadron::operator=(const G4QHadron &right)
195{
196  if(this != &right)                          // Beware of self assignment
197  {
198    theMomentum         = right.theMomentum;
199    theQPDG             = right.theQPDG;
200    valQ                = right.valQ;
201    nFragm              = right.nFragm;
202    thePosition         = right.thePosition;     
203    theCollisionCount   = 0;
204    isSplit             = false;
205    Direction           = right.Direction;
206    bindE               = right.bindE;
207  }
208  return *this;
209}
210
211G4QHadron::~G4QHadron()
212{
213  std::list<G4QParton*>::iterator ipos = Color.begin();
214  std::list<G4QParton*>::iterator epos = Color.end();
215  for( ; ipos != epos; ipos++) {delete [] *ipos;}
216  Color.clear();
217
218  ipos = AntiColor.begin();
219  epos = AntiColor.end();
220  for( ; ipos != epos; ipos++) {delete [] *ipos;}
221  AntiColor.clear();
222}
223
224// Define quark content of the particle with a particular PDG Code
225void G4QHadron::DefineQC(G4int PDGCode)
226{
227  //G4cout<<"G4QHadron::DefineQC is called with PDGCode="<<PDGCode<<G4endl;
228  G4int szn=PDGCode-90000000;
229  G4int ds=0;
230  G4int dz=0;
231  G4int dn=0;
232  if(szn<-100000)
233  {
234    G4int ns=(-szn)/1000000+1;
235    szn+=ns*1000000;
236    ds+=ns;
237  }
238  else if(szn<-100)
239  {
240    G4int nz=(-szn)/1000+1;
241    szn+=nz*1000;
242    dz+=nz;
243  }
244  else if(szn<0)
245  {
246    G4int nn=-szn;
247    szn=0;
248    dn+=nn;
249  }
250  G4int sz =szn/1000;
251  G4int n  =szn%1000;
252  if(n>700)
253  {
254    n-=1000;
255    dz--;
256  }
257  G4int z  =sz%1000-dz;
258  if(z>700)
259  {
260    z-=1000;
261    ds--;
262  }
263  G4int Sq =sz/1000-ds;
264  G4int zns=z+n+Sq;
265  G4int Dq=n+zns;
266  G4int Uq=z+zns;
267  if      (Dq<0&&Uq<0&&Sq<0)valQ=G4QContent(0 ,0 ,0 ,-Dq,-Uq,-Sq);
268  else if (Uq<0&&Sq<0)      valQ=G4QContent(Dq,0 ,0 ,0  ,-Uq,-Sq);
269  else if (Dq<0&&Sq<0)      valQ=G4QContent(0 ,Uq,0 ,-Dq,0  ,-Sq);
270  else if (Dq<0&&Uq<0)      valQ=G4QContent(0 ,0 ,Sq,-Dq,-Uq,0  );
271  else if (Uq<0)            valQ=G4QContent(Dq,0 ,Sq,0  ,-Uq,0  );
272  else if (Sq<0)            valQ=G4QContent(Dq,Uq,0 ,0  ,0  ,-Sq);
273  else if (Dq<0)            valQ=G4QContent(0 ,Uq,Sq,-Dq,0  ,0  );
274  else                      valQ=G4QContent(Dq,Uq,Sq,0  ,0  ,0  );
275}
276
277// Redefine a Hadron with a new PDGCode
278void G4QHadron::SetQPDG(const G4QPDGCode& newQPDG)
279{
280  theQPDG  = newQPDG;
281  G4int PDG= newQPDG.GetPDGCode();
282  G4int Q  = newQPDG.GetQCode();
283#ifdef debug
284  G4cout<<"G4QHadron::SetQPDG is called with PDGCode="<<PDG<<", QCode="<<Q<<G4endl;
285#endif
286  if     (Q>-1) valQ=theQPDG.GetQuarkContent();
287  else if(PDG>80000000) DefineQC(PDG);
288  else
289  {
290    G4cerr<<"***G4QHadron::SetQPDG: QPDG="<<newQPDG<<G4endl;
291    throw G4QException("***G4QHadron::SetQPDG: Impossible QPDG Probably a Chipolino");
292  }
293}
294
295// Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir
296G4bool 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)]
420G4bool 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)
544G4bool 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 pdebug
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());
605  }
606  G4double vP  = theMomentum.rho();      // Momentum of the decaying hadron
607  G4double dE  = theMomentum.e();        // Energy of the decaying hadron
608  if(dE<vP)
609  {
610    G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl;
611    G4double accuracy=.000001*vP;
612    G4double emodif=std::fabs(dE-vP);
613    if(emodif<accuracy)
614    {
615      G4cerr<<"G4QHadron::DecayIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl;
616      theMomentum.setE(vP+emodif+.01*accuracy);
617    }
618  }
619  G4ThreeVector ltb = theMomentum.boostVector(); // Boost vector for backward Lor.Trans.
620#ifdef pdebug
621  G4cout<<"G4QHadron::DecIn2:LorTrans v="<<ltb<<",f4Mom="<<f4Mom<<",s4Mom="<<s4Mom<<G4endl;
622#endif
623  if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* f4M="<<f4Mom<<G4endl;
624  f4Mom.boost(ltb);                        // Lor.Trans. of 1st hadron back to LS
625  if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::DecIn2:*Boost* s4M="<<s4Mom<<G4endl; 
626  s4Mom.boost(ltb);                        // Lor.Trans. of 2nd hadron back to LS
627#ifdef pdebug
628  G4cout<<"G4QHadron::DecayIn2: ROOT OUTPUT f4Mom="<<f4Mom<<", s4Mom="<<s4Mom<<G4endl;
629#endif
630  return true;
631} // End of "DecayIn2"
632
633// Correction for the Hadron + fr decay in case of the new corM mass of the Hadron
634G4bool G4QHadron::CorMDecayIn2(G4double corM, G4LorentzVector& fr4Mom)
635{//    ===============================================================
636  G4double fM  = fr4Mom.m();                // Mass of the Fragment
637  G4LorentzVector comp=theMomentum+fr4Mom;  // 4Mom of the decaying compound system
638  G4double iM  = comp.m();                  // mass of the decaying compound system
639#ifdef pdebug
640  G4cout<<"G4QH::CMDIn2: iM="<<iM<<comp<<"=>fM="<<fM<<"+corM="<<corM<<"="<<fM+corM<<G4endl;
641#endif
642  G4double dE=iM-fM-corM;
643  //@@ Later on make a quark content check for the decay
644  if (fabs(dE)<.001)
645  {
646    G4double fR=fM/iM;
647    G4double cR=corM/iM;
648    fr4Mom=fR*comp;
649    theMomentum=cR*comp;
650    return true;
651  }
652  else if (dE<-.001 || iM==0.)
653  {
654    G4cerr<<"***G4QH::CorMDIn2***fM="<<fM<<" + cM="<<corM<<" > iM="<<iM<<",d="<<dE<<G4endl;
655    return false;
656  }
657  G4double corM2= corM*corM;
658  G4double fM2 = fM*fM;
659  G4double iM2 = iM*iM;
660  G4double d2 = iM2-fM2-corM2;
661  G4double p2 = (d2*d2/4.-fM2*corM2)/iM2;    // Decay momentum(^2) in CMS of Quasmon
662  if (p2<0.)
663  {
664#ifdef pdebug
665    G4cerr<<"**G4QH::CMDI2:p2="<<p2<<"<0,d="<<d2*d2/4.<<"<4*fM2*hM2="<<4*fM2*corM2<<G4endl;
666#endif
667    p2=0.;
668  }
669  G4double p  = sqrt(p2);
670  if(comp.e()<comp.rho())G4cerr<<"*G4QH::CorMDecayIn2:*Boost* comp4M="<<comp<<",e-p="
671                               <<comp.e()-comp.rho()<<G4endl;
672  G4ThreeVector ltb = comp.boostVector();      // Boost vector for backward Lor.Trans.
673  G4ThreeVector ltf = -ltb;                    // Boost vector for forward Lorentz Trans.
674  G4LorentzVector cm4Mom=fr4Mom;               // Copy of fragment 4Mom to transform to CMS
675  if(cm4Mom.e()<cm4Mom.rho())
676  {
677    G4cerr<<"*G4QH::CorMDecIn2:*Boost* c4M="<<cm4Mom<<G4endl;
678    cm4Mom.setE(1.0000001*cm4Mom.rho());
679  }
680  cm4Mom.boost(ltf);                           // Now it is in CMS (Forward Lor.Trans.)
681  G4double pfx= cm4Mom.px();
682  G4double pfy= cm4Mom.py();
683  G4double pfz= cm4Mom.pz();
684  G4double pt2= pfx*pfx+pfy*pfy;
685  G4double tx=0.;
686  G4double ty=0.;
687  if(pt2<=0.)
688  {
689    G4double phi= 360.*deg*G4UniformRand();  // @@ Change 360.*deg to M_TWOPI (?)
690    tx=sin(phi);
691    ty=cos(phi);
692  }
693  else
694  {
695    G4double pt=sqrt(pt2);
696    tx=pfx/pt;
697    ty=pfy/pt;
698  }
699  G4double pc2=pt2+pfz*pfz;
700  G4double ct=0.;
701  if(pc2<=0.)
702  {
703    G4double rnd= G4UniformRand();
704    ct=1.-rnd-rnd;
705  }
706  else
707  {
708    G4double pc=sqrt(pc2);
709    ct=pfz/pc;
710  }
711#ifdef debug
712  G4cout<<"G4QHadron::CorMDecayIn2: ct="<<ct<<", p="<<p<<G4endl;
713#endif
714  G4double ps = p * sqrt(1.-ct*ct);
715  G4ThreeVector pVect(ps*tx,ps*ty,p*ct);
716  fr4Mom.setVect(pVect);
717  fr4Mom.setE(sqrt(fM2+p2));
718  theMomentum.setVect((-1)*pVect);
719  theMomentum.setE(sqrt(corM2+p2));
720#ifdef pdebug
721  G4LorentzVector dif2=comp-fr4Mom-theMomentum;
722  G4cout<<"G4QH::CorMDIn2:c="<<comp<<"-f="<<fr4Mom<<"-4M="<<theMomentum<<"="<<dif2<<G4endl;
723#endif
724  if(fr4Mom.e()+.001<fr4Mom.rho())G4cerr<<"*G4QH::CorMDecIn2:*Boost*fr4M="<<fr4Mom<<G4endl;
725  fr4Mom.boost(ltb);                        // Lor.Trans. of the Fragment back to LS
726  if(theMomentum.e()<theMomentum.rho())
727  {
728    G4cerr<<"*G4QH::CMDI2:4="<<theMomentum<<G4endl;
729    theMomentum.setE(1.0000001*theMomentum.rho());
730  }
731  theMomentum.boost(ltb);                  // Lor.Trans. of the Hadron back to LS
732#ifdef pdebug
733  G4LorentzVector dif3=comp-fr4Mom-theMomentum;
734  G4cout<<"G4QH::CorMDecIn2:OUTPUT:f4M="<<fr4Mom<<",h4M="<<theMomentum<<"d="<<dif3<<G4endl;
735#endif
736  return true;
737} // End of "CorMDecayIn2"
738
739
740// Fragment fr4Mom louse energy corE and transfer it to This Hadron
741G4bool G4QHadron::CorEDecayIn2(G4double corE, G4LorentzVector& fr4Mom)
742{//    ===============================================================
743  G4double fE  = fr4Mom.m();                // Energy of the Fragment
744#ifdef pdebug
745  G4cout<<"G4QH::CorEDecIn2:fE="<<fE<<fr4Mom<<">corE="<<corE<<",h4M="<<theMomentum<<G4endl;
746#endif
747  if (fE+.001<=corE)
748  {
749#ifdef pdebug
750    G4cerr<<"***G4QHadron::CorEDecIn2*** fE="<<fE<<"<corE="<<corE<<", d="<<corE-fE<<G4endl;
751#endif
752    return false;
753  }
754  G4double fM2=fr4Mom.m2();                 // Squared Mass of the Fragment
755  if(fM2<0.) fM2=0.;
756  G4double iPx=fr4Mom.px();                 // Initial Px of the Fragment
757  G4double iPy=fr4Mom.py();                 // Initial Py of the Fragment
758  G4double iPz=fr4Mom.pz();                 // Initial Pz of the Fragment
759  G4double fP2=iPx*iPx+iPy*iPy+iPz*iPz;     // Initial Squared 3-momentum of the Fragment
760  G4double finE = fE - corE;                // Final energy of the fragment
761  G4double rP = sqrt((finE*finE-fM2)/fP2);  // Reduction factor for the momentum
762  G4double fPx=iPx*rP;
763  G4double fPy=iPy*rP;
764  G4double fPz=iPz*rP;
765  fr4Mom= G4LorentzVector(fPx,fPy,fPz,finE);
766  G4double Px=theMomentum.px()+iPx-fPx;
767  G4double Py=theMomentum.py()+iPy-fPy;
768  G4double Pz=theMomentum.pz()+iPz-fPz;
769  G4double mE=theMomentum.e();
770  ///////////G4double mM2=theMomentum.m2();
771  theMomentum= G4LorentzVector(Px,Py,Pz,mE+corE);
772#ifdef pdebug
773  G4double difF=fr4Mom.m2()-fM2;
774  G4cout<<"G4QH::CorEDecIn2: dF="<<difF<<",out:"<<theMomentum<<fr4Mom<<G4endl;
775#endif
776  return true;
777} // End of "CorEDecayIn2"
778
779// Decay of the hadron in 3 particles i=>r+s+t
780G4bool G4QHadron::DecayIn3
781                   (G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, G4LorentzVector& t4Mom)
782{//    ====================================================================================
783#ifdef debug
784  G4cout<<"G4QH::DIn3:"<<theMomentum<<"=>pf="<<f4Mom<<"+ps="<<s4Mom<<"+pt="<<t4Mom<<G4endl;
785#endif
786  G4double iM  = theMomentum.m();  // Mass of the decaying hadron
787  G4double fM  = f4Mom.m();        // Mass of the 1st hadron
788  G4double sM  = s4Mom.m();        // Mass of the 2nd hadron
789  G4double tM  = t4Mom.m();        // Mass of the 3rd hadron
790  G4double eps = 0.001;            // Accuracy of the split condition
791  if (fabs(iM-fM-sM-tM)<=eps)
792  {
793    G4double fR=fM/iM;
794    G4double sR=sM/iM;
795    G4double tR=tM/iM;
796    f4Mom=fR*theMomentum;
797    s4Mom=sR*theMomentum;
798    t4Mom=tR*theMomentum;
799    return true;
800  }
801  if (iM+eps<fM+sM+tM)
802  {
803    G4cout<<"***G4QHadron::DecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
804          <<",d="<<iM-fM-sM-tM<<G4endl;
805    return false;
806  }
807  G4double fM2 = fM*fM;
808  G4double sM2 = sM*sM;
809  G4double tM2 = tM*tM;
810  G4double iM2 = iM*iM;
811  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
812  G4double m12sMin =(fM+sM)*(fM+sM);
813  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
814  G4double rR = 0.;
815  G4double rnd= 1.;
816#ifdef debug
817  G4int    tr = 0;                 //@@ Comment if "cout" below is skiped @@
818#endif
819  G4double m12s = 0.;              // Fake definition before the Loop
820  while (rnd > rR)
821  {
822    m12s = m12sMin + m12sBase*G4UniformRand();
823    G4double e1=m12s+fM2-sM2;
824    G4double e2=iM2-m12s-tM2;
825    G4double four12=4*m12s;
826    G4double m13sRange=0.;
827    G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
828    if(dif<0.)
829    {
830#ifdef debug
831      if(dif<-.01) G4cerr<<"*G4QHadron::DecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
832                         <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
833#endif
834    }
835    else m13sRange=sqrt(dif)/m12s;
836    rR = m13sRange/m13sBase;
837    rnd= G4UniformRand();
838#ifdef debug
839    G4cout<<"G4QHadron::DecayIn3: try to decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
840#endif
841  }
842  G4double m12 = sqrt(m12s);       // Mass of the H1+H2 system
843  G4LorentzVector dh4Mom(0.,0.,0.,m12);
844 
845  if(!DecayIn2(t4Mom,dh4Mom))
846  {
847    G4cerr<<"***G4QHadron::DecayIn3: Exception1"<<G4endl;
848    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
849    return false;
850  }
851#ifdef debug
852  G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
853#endif
854  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
855  {
856    G4cerr<<"***G4QHadron::DecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
857    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
858    return false;
859  }
860  return true;
861} // End of DecayIn3
862
863// Relative Decay of the hadron in 3 particles i=>f+s+t (t is with respect to minC<ct<maxC)
864G4bool G4QHadron::RelDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
865                              G4LorentzVector& t4Mom, G4LorentzVector& dir,
866                              G4double maxCost, G4double minCost)
867{//    ====================================================================================
868#ifdef debug
869  G4cout<<"G4QH::RelDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
870#endif
871  G4double iM  = theMomentum.m();  // Mass of the decaying hadron
872  G4double fM  = f4Mom.m();        // Mass of the 1st hadron
873  G4double sM  = s4Mom.m();        // Mass of the 2nd hadron
874  G4double tM  = t4Mom.m();        // Mass of the 3rd hadron
875  G4double eps = 0.001;            // Accuracy of the split condition
876  if (fabs(iM-fM-sM-tM)<=eps)
877  {
878    G4double fR=fM/iM;
879    G4double sR=sM/iM;
880    G4double tR=tM/iM;
881    f4Mom=fR*theMomentum;
882    s4Mom=sR*theMomentum;
883    t4Mom=tR*theMomentum;
884    return true;
885  }
886  if (iM+eps<fM+sM+tM)
887  {
888    G4cout<<"***G4QHadron::RelDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
889          <<",d="<<iM-fM-sM-tM<<G4endl;
890    return false;
891  }
892  G4double fM2 = fM*fM;
893  G4double sM2 = sM*sM;
894  G4double tM2 = tM*tM;
895  G4double iM2 = iM*iM;
896  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
897  G4double m12sMin =(fM+sM)*(fM+sM);
898  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
899  G4double rR = 0.;
900  G4double rnd= 1.;
901#ifdef debug
902  G4int    tr = 0;                 //@@ Comment if "cout" below is skiped @@
903#endif
904  G4double m12s = 0.;              // Fake definition before the Loop
905  while (rnd > rR)
906  {
907    m12s = m12sMin + m12sBase*G4UniformRand();
908    G4double e1=m12s+fM2-sM2;
909    G4double e2=iM2-m12s-tM2;
910    G4double four12=4*m12s;
911    G4double m13sRange=0.;
912    G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
913    if(dif<0.)
914    {
915#ifdef debug
916      if(dif<-.01) G4cerr<<"G4QHadron::RelDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
917                         <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
918#endif
919    }
920    else m13sRange=sqrt(dif)/m12s;
921    rR = m13sRange/m13sBase;
922    rnd= G4UniformRand();
923#ifdef debug
924    G4cout<<"G4QHadron::RelDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
925#endif
926  }
927  G4double m12 = sqrt(m12s);       // Mass of the H1+H2 system
928  G4LorentzVector dh4Mom(0.,0.,0.,m12);
929 
930  if(!RelDecayIn2(t4Mom,dh4Mom,dir,maxCost,minCost))
931  {
932    G4cerr<<"***G4QHadron::RelDecayIn3: Exception1"<<G4endl;
933    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
934    return false;
935  }
936#ifdef debug
937  G4cout<<"G4QHadron::RelDecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
938#endif
939  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
940  {
941    G4cerr<<"***G4QHadron::RelDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
942    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
943    return false;
944  }
945  return true;
946} // End of RelDecayIn3
947
948// Relative Decay of hadron in 3: i=>f+s+t.  dN/dO [cp>0:~cost^cp, cp<0:~(1-cost)^(-cp)]
949G4bool G4QHadron::CopDecayIn3(G4LorentzVector& f4Mom, G4LorentzVector& s4Mom,
950                              G4LorentzVector& t4Mom, G4LorentzVector& dir, G4double cosp)
951{//    ====================================================================================
952#ifdef debug
953  G4cout<<"G4QH::CopDIn3:"<<theMomentum<<"=>f="<<f4Mom<<"+s="<<s4Mom<<"+t="<<t4Mom<<G4endl;
954#endif
955  G4double iM  = theMomentum.m();  // Mass of the decaying hadron
956  G4double fM  = f4Mom.m();        // Mass of the 1st hadron
957  G4double sM  = s4Mom.m();        // Mass of the 2nd hadron
958  G4double tM  = t4Mom.m();        // Mass of the 3rd hadron
959  G4double eps = 0.001;            // Accuracy of the split condition
960  if (fabs(iM-fM-sM-tM)<=eps)
961  {
962    G4double fR=fM/iM;
963    G4double sR=sM/iM;
964    G4double tR=tM/iM;
965    f4Mom=fR*theMomentum;
966    s4Mom=sR*theMomentum;
967    t4Mom=tR*theMomentum;
968    return true;
969  }
970  if (iM+eps<fM+sM+tM)
971  {
972    G4cout<<"***G4QHadron::CopDecayIn3:fM="<<fM<<" + sM="<<sM<<" + tM="<<tM<<" > iM="<<iM
973          <<",d="<<iM-fM-sM-tM<<G4endl;
974    return false;
975  }
976  G4double fM2 = fM*fM;
977  G4double sM2 = sM*sM;
978  G4double tM2 = tM*tM;
979  G4double iM2 = iM*iM;
980  G4double m13sBase=(iM-sM)*(iM-sM)-(fM+tM)*(fM+tM);
981  G4double m12sMin =(fM+sM)*(fM+sM);
982  G4double m12sBase=(iM-tM)*(iM-tM)-m12sMin;
983  G4double rR = 0.;
984  G4double rnd= 1.;
985#ifdef debug
986  G4int    tr = 0;                 //@@ Comment if "cout" below is skiped @@
987#endif
988  G4double m12s = 0.;              // Fake definition before the Loop
989  while (rnd > rR)
990  {
991    m12s = m12sMin + m12sBase*G4UniformRand();
992    G4double e1=m12s+fM2-sM2;
993    G4double e2=iM2-m12s-tM2;
994    G4double four12=4*m12s;
995    G4double m13sRange=0.;
996    G4double dif=(e1*e1-four12*fM2)*(e2*e2-four12*tM2);
997    if(dif<0.)
998    {
999#ifdef debug
1000      if(dif<-.01) G4cerr<<"G4QHadron::CopDecayIn3:iM="<<iM<<",tM="<<tM<<",sM="<<sM<<",fM="
1001                         <<fM<<",m12(s+f)="<<sqrt(m12s)<<", d="<<iM-fM-sM-tM<<G4endl;
1002#endif
1003    }
1004    else m13sRange=sqrt(dif)/m12s;
1005    rR = m13sRange/m13sBase;
1006    rnd= G4UniformRand();
1007#ifdef debug
1008    G4cout<<"G4QHadron::CopDecayIn3: try decay #"<<++tr<<", rR="<<rR<<",rnd="<<rnd<<G4endl;
1009#endif
1010  }
1011  G4double m12 = sqrt(m12s);       // Mass of the H1+H2 system
1012  G4LorentzVector dh4Mom(0.,0.,0.,m12);
1013 
1014  if(!CopDecayIn2(t4Mom,dh4Mom,dir,cosp))
1015  {
1016    G4cerr<<"***G4QHadron::CopDecayIn3: Exception1"<<G4endl;
1017    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
1018    return false;
1019  }
1020#ifdef debug
1021  G4cout<<"G4QHadron::DecayIn3: Now the last decay of m12="<<dh4Mom.m()<<G4endl;
1022#endif
1023  if(!G4QHadron(dh4Mom).DecayIn2(f4Mom,s4Mom))
1024  {
1025    G4cerr<<"***G4QHadron::CopDecayIn3: Error in DecayIn2 -> Exception2"<<G4endl;
1026    //throw G4QException("G4QHadron::DecayIn3(): DecayIn2 did not succeed");
1027    return false;
1028  }
1029  return true;
1030} // End of CopDecayIn3
1031
1032// Randomize particle mass taking into account the width
1033G4double G4QHadron::RandomizeMass(G4QParticle* pPart, G4double maxM)
1034//       ===========================================================
1035{
1036  G4double meanM = theQPDG.GetMass();
1037  G4double width = theQPDG.GetWidth()/2.;
1038#ifdef debug
1039  G4cout<<"G4QHadron::RandomizeMass: meanM="<<meanM<<", halfWidth="<<width<<G4endl;
1040#endif
1041  if(maxM<meanM-3*width) 
1042  {
1043#ifdef debug
1044    G4cout<<"***G4QH::RandM:m=0 maxM="<<maxM<<"<meanM="<<meanM<<"-3*halfW="<<width<<G4endl;
1045#endif
1046    return 0.;
1047  }
1048  ///////////////G4double theMass  = 0.;
1049  if(width==0.)
1050  {
1051#ifdef debug
1052    if(meanM>maxM) G4cerr<<"***G4QHadron::RandM:Stable m="<<meanM<<">maxM="<<maxM<<G4endl;
1053#endif
1054    return meanM;
1055    //return 0.;
1056  }
1057  else if(width<0.)
1058  {
1059    G4cerr<<"***G4QHadron::RandM: width="<<width<<"<0,PDGC="<<theQPDG.GetPDGCode()<<G4endl;
1060    throw G4QException("G4QHadron::RandomizeMass: with the width of the Hadron < 0.");
1061  }
1062  G4double minM = pPart->MinMassOfFragm();
1063  if(minM>maxM)
1064  {
1065#ifdef debug
1066    G4cout<<"***G4QHadron::RandomizeMass:for PDG="<<theQPDG.GetPDGCode()<<" minM="<<minM
1067          <<" > maxM="<<maxM<<G4endl;
1068#endif
1069    return 0.;
1070  }
1071  //Now calculate the Breit-Wigner distribution with two cuts
1072  G4double v1=atan((minM-meanM)/width);
1073  G4double v2=atan((maxM-meanM)/width);
1074  G4double dv=v2-v1;
1075#ifdef debug
1076  G4cout<<"G4QHadr::RandM:Mi="<<minM<<",i="<<v1<<",Ma="<<maxM<<",a="<<v2<<","<<dv<<G4endl;
1077#endif
1078  return meanM+width*tan(v1+dv*G4UniformRand());
1079}
1080
1081// Split hadron in partons
1082void G4QHadron::SplitUp()
1083{ 
1084  if (IsSplit()) return;
1085  Splitting();
1086  if (Color.empty()) return;
1087  if (GetSoftCollisionCount() == 0)
1088  {
1089    // Diffractive splitting: take the particle definition and get the partons
1090    G4QParton* Left = 0;
1091    G4QParton* Right = 0;
1092    GetValenceQuarkFlavors(Left, Right);
1093    Left->SetPosition(GetPosition());
1094    Right->SetPosition(GetPosition());
1095 
1096    G4LorentzVector HadronMom = Get4Momentum();
1097    //G4cout<<"DSU 1 - "<<HadronMom<<G4endl;
1098
1099    // momenta of string ends
1100    G4double pt2 = HadronMom.perp2();
1101    G4double transverseMass2 = HadronMom.plus()*HadronMom.minus();
1102    G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
1103    G4ThreeVector pt(minTransverseMass, minTransverseMass, 0);
1104    if(maxAvailMomentum2/widthOfPtSquare>0.01)
1105           pt=GaussianPt(widthOfPtSquare, maxAvailMomentum2);
1106    //G4cout<<"DSU 1.1 - "<<maxAvailMomentum2<<", pt="<<pt<<G4endl;
1107
1108    G4LorentzVector LeftMom(pt, 0.);
1109    G4LorentzVector RightMom;
1110    RightMom.setPx(HadronMom.px() - pt.x());
1111    RightMom.setPy(HadronMom.py() - pt.y());
1112    //G4cout<<"DSU 2: Right4M="<<RightMom<<", Left4M= "<<LeftMom<<G4endl;
1113
1114    G4double Local1 = HadronMom.minus() +
1115                      (RightMom.perp2() - LeftMom.perp2()) / HadronMom.plus();
1116    G4double Local2 = std::sqrt(std::max(0., Local1*Local1 -
1117                                4*RightMom.perp2()*HadronMom.minus() / HadronMom.plus()));
1118    //G4cout<<"DSU 3: L1="<< Local1 <<", L2="<<Local2<<G4endl;
1119
1120    if (Direction) Local2 = -Local2;
1121    G4double RightMinus   = 0.5*(Local1 + Local2);
1122    G4double LeftMinus = HadronMom.minus() - RightMinus;
1123    //G4cout<<"DSU 4: Rm="<<RightMinus<<", Lm="<<LeftMinus<<" "<<HadronMom.minus()<<G4endl;
1124
1125    G4double LeftPlus  = LeftMom.perp2()/LeftMinus;
1126    G4double RightPlus = HadronMom.plus() - LeftPlus;
1127    //G4cout<<"DSU 5: Rp="<<RightPlus<<", Lp="<<LeftPlus<<G4endl;
1128    LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
1129    LeftMom.setE (0.5*(LeftPlus + LeftMinus));
1130    RightMom.setPz(0.5*(RightPlus - RightMinus));
1131    RightMom.setE (0.5*(RightPlus + RightMinus));
1132    //G4cout<<"DSU 6: Left4M="<<LeftMom<<", Right4M="<<RightMom<<G4endl;
1133    Left->Set4Momentum(LeftMom);
1134    Right->Set4Momentum(RightMom);
1135    Color.push_back(Left);
1136    AntiColor.push_back(Right);
1137  }
1138  else
1139  {
1140    // Soft hadronization splitting: sample transversal momenta for sea and valence quarks
1141    G4double phi, pts;
1142    G4double SumPy = 0.;
1143    G4double SumPx = 0.;
1144    G4ThreeVector Pos = GetPosition();
1145    G4int nSeaPair = GetSoftCollisionCount()-1; 
1146   
1147    // here the condition,to ensure viability of splitting, also in cases
1148    // where difractive excitation occured together with soft scattering.
1149    //G4double LightConeMomentum=(Direction)? Get4Momentum().plus():Get4Momentum().minus();
1150    //G4double Xmin = theMinPz/LightConeMomentum;
1151    G4double Xmin = theMinPz/( Get4Momentum().e() - GetMass() );
1152    while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95;
1153
1154    G4int aSeaPair;
1155    for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
1156    {
1157      //  choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
1158      G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); 
1159
1160      //  BuildSeaQuark() determines quark spin, isospin and colour
1161      //  via parton-constructor G4QParton(aPDGCode)
1162      G4QParton* aParton = BuildSeaQuark(false, aPDGCode);
1163
1164      //  G4cout << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl;
1165
1166      //  G4cout << "Parton 1: "
1167      //         << " PDGcode: "  << aPDGCode
1168      //         << " - Name: "   << aParton->GetDefinition()->GetParticleName()
1169      //         << " - Type: "   << aParton->GetDefinition()->GetParticleType()
1170      //         << " - Spin-3: " << aParton->GetSpinZ()
1171      //         << " - Colour: " << aParton->GetColour() << G4endl;
1172
1173      // save colour a spin-3 for anti-quark
1174      G4int firstPartonColour = aParton->GetColour();
1175      G4double firstPartonSpinZ = aParton->GetSpinZ();
1176
1177      SumPx += aParton->Get4Momentum().px();
1178      SumPy += aParton->Get4Momentum().py();
1179      Color.push_back(aParton);
1180
1181      // create anti-quark
1182      aParton = BuildSeaQuark(true, aPDGCode);
1183      aParton->SetSpinZ(-firstPartonSpinZ);
1184      aParton->SetColour(-firstPartonColour);
1185
1186      //  G4cout << "Parton 2: "
1187      //         << " PDGcode: "  << -aPDGCode
1188      //         << " - Name: "   << aParton->GetDefinition()->GetParticleName()
1189      //         << " - Type: "   << aParton->GetDefinition()->GetParticleType()
1190      //         << " - Spin-3: " << aParton->GetSpinZ()
1191      //         << " - Colour: " << aParton->GetColour() << G4endl;
1192      //  G4cerr << "------------" << G4endl;
1193
1194      SumPx += aParton->Get4Momentum().px();
1195      SumPy += aParton->Get4Momentum().py();
1196      AntiColor.push_back(aParton);
1197    }
1198    // Valence quark   
1199    G4QParton* pColorParton = 0;   
1200    G4QParton* pAntiColorParton = 0;   
1201    GetValenceQuarkFlavors(pColorParton, pAntiColorParton);
1202    G4int ColorEncoding = pColorParton->GetPDGCode();
1203    G4int AntiColorEncoding = pAntiColorParton->GetPDGCode();
1204   
1205    pts   =  sigmaPt*std::sqrt(-std::log(G4UniformRand()));
1206    phi   = twopi*G4UniformRand();
1207    G4double Px = pts*std::cos(phi);
1208    G4double Py = pts*std::sin(phi);
1209    SumPx += Px;
1210    SumPy += Py;
1211
1212    if (ColorEncoding < 0) // use particle definition
1213    {
1214      G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0);
1215      pColorParton->Set4Momentum(ColorMom);
1216      G4LorentzVector AntiColorMom(Px, Py, 0, 0);
1217      pAntiColorParton->Set4Momentum(AntiColorMom);
1218    }
1219    else
1220    {
1221      G4LorentzVector ColorMom(Px, Py, 0, 0);
1222      pColorParton->Set4Momentum(ColorMom);
1223      G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0);
1224      pAntiColorParton->Set4Momentum(AntiColorMom);
1225    }
1226    Color.push_back(pColorParton);
1227    AntiColor.push_back(pAntiColorParton);
1228
1229    // Sample X
1230    G4int nAttempt = 0;
1231    G4double SumX = 0;
1232    G4double aBeta = beta;
1233    G4double ColorX, AntiColorX;
1234    G4double HPWtest = 0;
1235    G4int aPDG=std::abs(GetPDGCode());
1236    if (aPDG ==211 || aPDG == 22 || aPDG == 111) aBeta = 1.;       
1237    else if (aPDG == 321) aBeta = 0.;       
1238    else G4cout<<"-Warning-G4QHadron::SplitUp: wrong PDG="<<GetPDGCode()<<G4endl;
1239    do
1240    {
1241      SumX = 0;
1242      nAttempt++;
1243      G4int    NumberOfUnsampledSeaQuarks = 2*nSeaPair;
1244      G4double beta1 = beta;
1245      if (std::abs(ColorEncoding) <= 1000 && std::abs(AntiColorEncoding) <= 1000) beta1 = 1.; //...  in a meson       
1246      ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
1247      HPWtest = ColorX;
1248      //while (ColorX < Xmin || ColorX > 1. || 1. -  ColorX <= Xmin)
1249      //{
1250      //}
1251      Color.back()->SetX(SumX = ColorX); // this is the valenz quark.
1252
1253      std::list<G4QParton*>::iterator icolor = Color.begin();
1254      std::list<G4QParton*>::iterator ecolor = Color.end();
1255      std::list<G4QParton*>::iterator ianticolor = AntiColor.begin();
1256      std::list<G4QParton*>::iterator eanticolor = AntiColor.end();
1257      for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor)
1258      {
1259        NumberOfUnsampledSeaQuarks--;
1260        ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
1261        (*icolor)->SetX(ColorX);
1262        SumX += ColorX; 
1263        NumberOfUnsampledSeaQuarks--;
1264        AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
1265        (*ianticolor)->SetX(AntiColorX); // the 'sea' partons
1266        SumX += AntiColorX;
1267        if (1. - SumX <= Xmin)  break;
1268      }
1269    } while (1. - SumX <= Xmin); 
1270    AntiColor.back()->SetX(1.0 - SumX); // the di-quark takes the rest, then go to momentum
1271    // and here is the bug ;-) @@@@@@@@@@@@@
1272    if(getenv("debug_QGSMSplitableHadron"))
1273                         G4cout<<"particle energy at split = "<<Get4Momentum().t()<<G4endl;
1274    G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus());
1275    // lightCone -= 0.5*Get4Momentum().m();
1276    // hpw testing @@@@@ lightCone = 2.*Get4Momentum().t();
1277    if(getenv("debug_QGSMSplitableHadron") )G4cout << "Light cone = "<<lightCone<<G4endl;
1278    std::list<G4QParton*>::iterator icolor = Color.begin();
1279    std::list<G4QParton*>::iterator ecolor = Color.end();
1280    std::list<G4QParton*>::iterator ianticolor = AntiColor.begin();
1281    std::list<G4QParton*>::iterator eanticolor = AntiColor.end();
1282    for ( ; icolor != ecolor && ianticolor != eanticolor; ++icolor, ++ianticolor)
1283    {
1284      (*icolor)->DefineMomentumInZ(lightCone, Direction);
1285      (*ianticolor)->DefineMomentumInZ(lightCone, Direction);
1286    }
1287    //G4cout <<G4endl<<"XSAMPLE "<<HPWtest<<G4endl;
1288    return;
1289  }
1290} // End of SplitUp
1291
1292// Boost hadron 4-momentum, using Boost Lorentz vector
1293void 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
1303G4QParton* 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
1314G4double G4QHadron::SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
1315{
1316  G4double result;
1317  G4double x1, x2;
1318  G4double ymax = 0;
1319  for(G4int ii=0; ii<100; ii++)                    // @@ 100 is hardwired ? M.K.
1320  {
1321    G4double y = std::pow(1./G4double(ii), alpha);
1322    y*=std::pow(std::pow(1.-anXmin-totalSea*anXmin,alpha+1)-std::pow(anXmin,alpha+1),nSea);
1323    y*=std::pow(1.-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
1324    if(y>ymax) ymax = y;
1325  }
1326  G4double y;
1327  G4double xMax=1.-(totalSea+1.)*anXmin;
1328  if(anXmin > xMax) 
1329  {
1330    G4cerr<<"***G4QHadron::SampleX: anXmin="<<anXmin<<" > xMax="<<xMax<<", nSea="<<nSea
1331          <<", totSea="<<totalSea<<G4endl;
1332    G4Exception("G4QHadron::SampleX:","72",FatalException,"TooBigXValue");
1333  }
1334  do
1335  {
1336    x1 = CLHEP::RandFlat::shoot(anXmin, xMax);
1337    y = std::pow(x1, alpha);
1338    y*=std::pow(std::pow(1.-x1-totalSea*anXmin,alpha+1) - std::pow(anXmin, alpha+1), nSea);
1339    y*=std::pow(1.-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1); 
1340    x2 = ymax*G4UniformRand();
1341  } while(x2>y);
1342  result = x1;
1343  return result; 
1344} // End of SampleX
1345
1346
1347void G4QHadron::GetValenceQuarkFlavors(G4QParton* &Parton1, G4QParton* &Parton2)
1348{
1349  // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
1350  G4int aEnd=0;
1351  G4int bEnd=0;
1352  G4int HadronEncoding = GetPDGCode();
1353  if(!(GetBaryonNumber())) SplitMeson(HadronEncoding,&aEnd,&bEnd);
1354  else SplitBaryon(HadronEncoding, &aEnd, &bEnd);
1355
1356  Parton1 = new G4QParton(aEnd);
1357  Parton1->SetPosition(GetPosition());
1358
1359// G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
1360// G4cerr << "Parton 1: "
1361//        << " PDGcode: "  << aEnd
1362//        << " - Name: "   << Parton1->GetDefinition()->GetParticleName()
1363//        << " - Type: "   << Parton1->GetDefinition()->GetParticleType()
1364//        << " - Spin-3: " << Parton1->GetSpinZ()
1365//        << " - Colour: " << Parton1->GetColour() << G4endl;
1366
1367  Parton2 = new G4QParton(bEnd);
1368  Parton2->SetPosition(GetPosition());
1369
1370// G4cerr << "Parton 2: "
1371//        << " PDGcode: "  << bEnd
1372//        << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
1373//        << " - Type: "   << Parton2->GetDefinition()->GetParticleType()
1374//        << " - Spin-3: " << Parton2->GetSpinZ()
1375//        << " - Colour: " << Parton2->GetColour() << G4endl;
1376// G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;
1377
1378  // colour of parton 1 choosen at random by G4QParton(aEnd)
1379  // colour of parton 2 is the opposite:
1380
1381  Parton2->SetColour(-(Parton1->GetColour()));
1382
1383 // isospin-3 of both partons is handled by G4Parton(PDGCode)
1384
1385 // spin-3 of parton 1 and 2 choosen at random by G4QParton(aEnd)
1386 // spin-3 of parton 2 may be constrained by spin of original particle:
1387
1388  if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > GetSpin()) 
1389                                                 Parton2->SetSpinZ(-(Parton2->GetSpinZ()));
1390
1391// G4cerr << "Parton 2: "
1392//        << " PDGcode: "  << bEnd
1393//        << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
1394//        << " - Type: "   << Parton2->GetDefinition()->GetParticleType()
1395//        << " - Spin-3: " << Parton2->GetSpinZ()
1396//        << " - Colour: " << Parton2->GetColour() << G4endl;
1397// G4cerr << "------------" << G4endl;
1398
1399} // End of GetValenceQuarkFlavors
1400
1401G4bool G4QHadron::SplitMeson(G4int PDGcode, G4int* aEnd, G4int* bEnd)
1402{
1403  G4bool result = true;
1404  G4int absPDGcode = std::abs(PDGcode);
1405  if (absPDGcode >= 1000) return false;
1406  if(absPDGcode == 22)
1407  {
1408    G4int it=1;
1409    if(G4UniformRand()<.5) it++;
1410    *aEnd = it;
1411    *bEnd = -it;
1412  }
1413  else
1414  {
1415    G4int heavy =  absPDGcode/100;
1416    G4int light = (absPDGcode%100)/10;
1417    G4int anti  = 1 - 2*(std::max(heavy, light)%2);
1418    if (PDGcode < 0 ) anti = -anti;
1419    heavy *=  anti;
1420    light *= -anti;
1421    if ( anti < 0) G4SwapObj(&heavy, &light);
1422    *aEnd = heavy;
1423    *bEnd = light;
1424  }
1425  return result;
1426}
1427
1428G4bool G4QHadron::SplitBaryon(G4int PDGcode, G4int* quark, G4int* diQuark)
1429{
1430  static const G4double r2=.5;
1431  static const G4double r3=1./3.;
1432  static const G4double d3=2./3.;
1433  static const G4double r4=1./4.;
1434  static const G4double r6=1./6.;
1435  static const G4double r12=1./12.;
1436  //
1437  std::pair<G4int,G4int> qdq[5];
1438  G4double               prb[5];
1439  G4int                  nc=0;
1440  G4int            aPDGcode=std::abs(PDGcode);
1441  if(aPDGcode==2212)                         // ==> Proton
1442  {
1443    nc=3;
1444    qdq[0]=make_pair(2203, 1); prb[0]=r3;    // uu_1, d
1445    qdq[1]=make_pair(2103, 2); prb[1]=r6;    // ud_1, u
1446    qdq[2]=make_pair(2101, 2); prb[2]=r2;    // ud_0, u
1447  }
1448  else if(aPDGcode==2112)                    // ==> Neutron
1449  {
1450    nc=3;
1451    qdq[0]=make_pair(2103, 1); prb[0]=r6;    // ud_1, d
1452    qdq[1]=make_pair(2101, 1); prb[1]=r2;    // ud_0, d
1453    qdq[2]=make_pair(1103, 2); prb[2]=r3;    // dd_1, u
1454  }
1455  else if(aPDGcode%10<3)                     // ==> Spin 1/2 Hyperons
1456  {
1457    if(aPDGcode==3122)
1458    {
1459      nc=5;
1460      qdq[0]=make_pair(2103, 3); prb[0]=r3;  // ud_1, s
1461      qdq[1]=make_pair(3203, 1); prb[1]=r4;  // su_1, d
1462      qdq[2]=make_pair(3201, 1); prb[2]=r12; // su_0, d
1463      qdq[3]=make_pair(3103, 2); prb[3]=r4;  // sd_1, u
1464      qdq[4]=make_pair(3101, 2); prb[4]=r12; // sd_0, u
1465    }
1466    else if(aPDGcode==3222)
1467    {
1468      nc=3;
1469      qdq[0]=make_pair(2203, 3); prb[0]=r3;
1470      qdq[1]=make_pair(3203, 2); prb[1]=r6;
1471      qdq[2]=make_pair(3201, 2); prb[2]=r2;
1472    }
1473    else if(aPDGcode==3212)
1474    {
1475      nc=5;
1476      qdq[0]=make_pair(2103, 3); prb[0]=r3;
1477      qdq[1]=make_pair(3203, 1); prb[1]=r12;
1478      qdq[2]=make_pair(3201, 1); prb[2]=r4;
1479      qdq[3]=make_pair(3103, 2); prb[3]=r12;
1480      qdq[4]=make_pair(3101, 2); prb[4]=r4;
1481    }
1482    else if(aPDGcode==3112)
1483    {
1484      nc=3;
1485      qdq[0]=make_pair(1103, 3); prb[0]=r3;
1486      qdq[1]=make_pair(3103, 1); prb[1]=r6;
1487      qdq[2]=make_pair(3101, 1); prb[2]=r2;
1488    }
1489    else if(aPDGcode==3312)
1490    {
1491      nc=3;
1492      qdq[0]=make_pair(3103, 3); prb[0]=r6;
1493      qdq[1]=make_pair(3101, 3); prb[1]=r2;
1494      qdq[2]=make_pair(3303, 1); prb[2]=r3;
1495    }
1496    else if(aPDGcode==3322)
1497    {
1498      nc=3;
1499      qdq[0]=make_pair(3203, 3); prb[0]=r6;
1500      qdq[1]=make_pair(3201, 3); prb[1]=r2;
1501      qdq[2]=make_pair(3303, 2); prb[2]=r3;
1502    }
1503    else return false;
1504  }
1505  else                                       // ==> Spin 3/2 Baryons (Spin>3/2 is ERROR)
1506  {
1507    if(aPDGcode==3334)
1508    {
1509      nc=1;
1510      qdq[0]=make_pair(3303, 3); prb[0]=1.;
1511    }
1512    else if(aPDGcode==2224)
1513    {
1514      nc=1;
1515      qdq[0]=make_pair(2203, 2); prb[0]=1.;
1516    }
1517    else if(aPDGcode==2214)
1518    {
1519      nc=2;
1520      qdq[0]=make_pair(2203, 1); prb[0]=r3;
1521      qdq[1]=make_pair(2103, 2); prb[1]=d3;
1522    }
1523    else if(aPDGcode==2114)
1524    {
1525      nc=2;
1526      qdq[0]=make_pair(2103, 1); prb[0]=d3;
1527      qdq[1]=make_pair(2103, 2); prb[1]=r3;
1528    }
1529    else if(aPDGcode==1114)
1530    {
1531      nc=1;
1532      qdq[0]=make_pair(1103, 1); prb[0]=1.;
1533    }
1534    else if(aPDGcode==3224)
1535    {
1536      nc=2;
1537      qdq[0]=make_pair(2203, 3); prb[0]=r3;
1538      qdq[1]=make_pair(3203, 2); prb[1]=d3;
1539    }
1540    else if(aPDGcode==3214)
1541    {
1542      nc=3;
1543      qdq[0]=make_pair(2103, 3); prb[0]=r3;
1544      qdq[1]=make_pair(3203, 1); prb[1]=r3;
1545      qdq[2]=make_pair(3103, 2); prb[2]=r3;
1546    }
1547    else if(aPDGcode==3114)
1548    {
1549      nc=2;
1550      qdq[0]=make_pair(1103, 3); prb[0]=r3;
1551      qdq[1]=make_pair(3103, 1); prb[1]=d3;
1552    }
1553    else if(aPDGcode==3324)
1554    {
1555      nc=2;
1556      qdq[0]=make_pair(3203, 3); prb[0]=r3;
1557      qdq[1]=make_pair(3303, 2); prb[1]=d3;
1558    }
1559    else if(aPDGcode==3314)
1560    {
1561      nc=2;
1562      qdq[0]=make_pair(3103, 3); prb[0]=d3;
1563      qdq[1]=make_pair(3303, 1); prb[1]=r3;
1564    }
1565    else return false;
1566  }
1567  G4double random = G4UniformRand();
1568  G4double sum = 0.;
1569  for(G4int i=0; i<nc; i++)
1570  {
1571    sum += prb[i];
1572    if(sum>random)
1573    {
1574      if(PDGcode>0)
1575      {
1576        *diQuark= qdq[i].first;
1577        *quark  = qdq[i].second;
1578      }
1579      else
1580      {
1581        *diQuark= -qdq[i].first;
1582        *quark  = -qdq[i].second;
1583      }
1584      break;
1585    }
1586  }
1587  return true;
1588}
1589
1590G4ThreeVector G4QHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
1591{
1592  G4double R=0.;
1593  while ((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare){}
1594  R = std::sqrt(R);
1595  G4double phi = twopi*G4UniformRand();
1596  return G4ThreeVector(R*std::cos(phi), R*std::sin(phi), 0.);   
1597}
Note: See TracBrowser for help on using the repository browser.