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

Last change on this file since 836 was 819, checked in by garnier, 16 years ago

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