Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/qgsm/src/G4QGSMSplitableHadron.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/parton_string/qgsm/src/G4QGSMSplitableHadron.cc (Version 11.3.0) and /processes/hadronic/models/parton_string/qgsm/src/G4QGSMSplitableHadron.cc (Version 8.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 #include "G4QGSMSplitableHadron.hh"                23 #include "G4QGSMSplitableHadron.hh"
 27 #include "G4PhysicalConstants.hh"              << 
 28 #include "G4SystemOfUnits.hh"                  << 
 29 #include "G4ParticleTable.hh"                      24 #include "G4ParticleTable.hh" 
 30 #include "G4PionPlus.hh"                           25 #include "G4PionPlus.hh"
 31 #include "G4PionMinus.hh"                          26 #include "G4PionMinus.hh"
 32 #include "G4Gamma.hh"                              27 #include "G4Gamma.hh"
 33 #include "G4PionZero.hh"                           28 #include "G4PionZero.hh"
 34 #include "G4KaonPlus.hh"                           29 #include "G4KaonPlus.hh"
 35 #include "G4KaonMinus.hh"                          30 #include "G4KaonMinus.hh"
 36                                                    31 
 37 #include "G4Log.hh"                            << 
 38 #include "G4Pow.hh"                            << 
 39                                                << 
 40 // based on prototype by Maxim Komogorov           32 // based on prototype by Maxim Komogorov
 41 // Splitting into methods, and centralizing of     33 // Splitting into methods, and centralizing of model parameters HPW Feb 1999
 42 // restructuring HPW Feb 1999                      34 // restructuring HPW Feb 1999
 43 // fixing bug in the sampling of 'x', HPW Feb      35 // fixing bug in the sampling of 'x', HPW Feb 1999
 44 // fixing bug in sampling pz, HPW Feb 1999.        36 // fixing bug in sampling pz, HPW Feb 1999. 
 45 // Code now also good for p-nucleus scattering     37 // Code now also good for p-nucleus scattering (before only p-p), HPW Feb 1999.
 46 // Using Parton more directly, HPW Feb 1999.       38 // Using Parton more directly, HPW Feb 1999.
 47 // Shortening the algorithm for sampling x, HP     39 // Shortening the algorithm for sampling x, HPW Feb 1999.
 48 // sampling of x replaced by formula, taking X     40 // sampling of x replaced by formula, taking X_min into account in the correlated sampling. HPW, Feb 1999.
 49 // logic much clearer now. HPW Feb 1999            41 // logic much clearer now. HPW Feb 1999
 50 // Removed the ordering problem. No Direction      42 // Removed the ordering problem. No Direction needed in selection of valence quark types. HPW Mar'99.
 51 // Fixing p-t distributions for scattering of      43 // Fixing p-t distributions for scattering of nuclei.
 52 // Separating out parameters.                      44 // Separating out parameters.
 53                                                    45 
 54 void G4QGSMSplitableHadron::InitParameters()       46 void G4QGSMSplitableHadron::InitParameters()
 55 {                                                  47 {
 56   // changing rapidity distribution for all        48   // changing rapidity distribution for all
 57   alpha = -0.5; // Note that this number is st     49   alpha = -0.5; // Note that this number is still assumed in the algorithm
 58   // needs to be generalized.                  <<  50                 // needs to be generalized.
 59   // changing rapidity distribution for projec     51   // changing rapidity distribution for projectile like
 60   beta = 2.5;// Note that this number is still     52   beta = 2.5;// Note that this number is still assumed in the algorithm
 61   // needs to be generalized.                  <<  53                 // needs to be generalized.
 62   theMinPz = 0.5*G4PionMinus::PionMinus()->Get <<  54   theMinPz = 0.5*G4PionMinus::PionMinus()->GetPDGMass(); 
 63   //  theMinPz = 0.1*G4PionMinus::PionMinus()- <<  55 //  theMinPz = 0.1*G4PionMinus::PionMinus()->GetPDGMass(); 
 64   //  theMinPz = G4PionMinus::PionMinus()->Get <<  56 //  theMinPz = G4PionMinus::PionMinus()->GetPDGMass(); 
 65   // as low as possible, otherwise, we have un     57   // as low as possible, otherwise, we have unphysical boundary conditions in the sampling.
 66   StrangeSuppress = 0.48;                          58   StrangeSuppress = 0.48;
 67   sigmaPt = 0.*GeV; // widens eta slightly, if <<  59   sigmaPt = 0.*GeV; // widens eta slightly, if increased to 1.7, 
 68   // but Maxim's algorithm breaks energy conse <<  60                     // but Maxim's algorithm breaks energy conservation
 69   widthOfPtSquare = 0.5*sqr(GeV);              <<  61         // to be revised.
                                                   >>  62   widthOfPtSquare = 0.01*GeV*GeV;
 70   Direction = FALSE;                               63   Direction = FALSE;
 71   minTransverseMass = 1*keV;                       64   minTransverseMass = 1*keV;
 72   iP  =0;//     Color.begin();                 << 
 73   iAP =0;// AntiColor.begin();                 << 
 74 }                                                  65 } 
 75                                                    66 
 76 G4QGSMSplitableHadron::G4QGSMSplitableHadron()     67 G4QGSMSplitableHadron::G4QGSMSplitableHadron() 
 77 {                                                  68 {
 78   InitParameters();                                69   InitParameters();
 79 }                                                  70 }
 80                                                    71 
 81 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c     72 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary, G4bool aDirection)
 82 :G4VSplitableHadron(aPrimary)                  <<  73  :G4VSplitableHadron(aPrimary)
 83 {                                                  74 {
 84   InitParameters();                                75   InitParameters();
 85   Direction = aDirection;                          76   Direction = aDirection;
 86 }                                                  77 }
 87                                                    78 
 88                                                <<  79  
 89 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c     80 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4ReactionProduct & aPrimary)
 90 :  G4VSplitableHadron(aPrimary)                <<  81       :  G4VSplitableHadron(aPrimary)
 91 {                                                  82 {
 92   InitParameters();                                83   InitParameters();
 93 }                                                  84 }
 94                                                    85 
 95 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c     86 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon)
 96 :  G4VSplitableHadron(aNucleon)                <<  87       :  G4VSplitableHadron(aNucleon)
 97 {                                                  88 {
 98   InitParameters();                                89   InitParameters();
 99 }                                                  90 }
100                                                    91 
101 G4QGSMSplitableHadron::G4QGSMSplitableHadron(c     92 G4QGSMSplitableHadron::G4QGSMSplitableHadron(const G4Nucleon & aNucleon, G4bool aDirection)
102 :  G4VSplitableHadron(aNucleon)                <<  93       :  G4VSplitableHadron(aNucleon)
103 {                                                  94 {
104   InitParameters();                                95   InitParameters();
105   Direction = aDirection;                          96   Direction = aDirection;
106 }                                                  97 }
107                                                    98 
108 G4QGSMSplitableHadron::~G4QGSMSplitableHadron( <<  99 G4QGSMSplitableHadron::~G4QGSMSplitableHadron(){}
109                                                   100 
                                                   >> 101 const G4QGSMSplitableHadron & G4QGSMSplitableHadron::operator=(const G4QGSMSplitableHadron &)
                                                   >> 102 {
                                                   >> 103   throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron::operator= meant to not be accessable");
                                                   >> 104   return *this;
                                                   >> 105 }
110                                                   106 
111 //******************************************** << 
112                                                   107 
                                                   >> 108 //**************************************************************************************************************************
                                                   >> 109     
113 void G4QGSMSplitableHadron::SplitUp()             110 void G4QGSMSplitableHadron::SplitUp()
114 {                                                 111 {  
115   if (IsSplit()) return;                          112   if (IsSplit()) return;
116   Splitting();  // To mark that a hadron is sp << 113   Splitting();
117   if (Color.size()!=0) return;                    114   if (Color.size()!=0) return;
118   if (GetSoftCollisionCount() == 0)  // GetSof << 115   if (GetSoftCollisionCount() == 0)
119   {                                               116   {
120     DiffractiveSplitUp();                         117     DiffractiveSplitUp();
121   }                                               118   }
122   else                                            119   else
123   {                                               120   {
124     SoftSplitUp();                                121     SoftSplitUp();
125   }                                               122   }
126 }                                                 123 }
127                                                << 124    
128 void G4QGSMSplitableHadron::DiffractiveSplitUp    125 void G4QGSMSplitableHadron::DiffractiveSplitUp()
129 {                                                 126 {   
130   // take the particle definitions and get the    127   // take the particle definitions and get the partons HPW
131   G4Parton * Left = NULL;                         128   G4Parton * Left = NULL;
132   G4Parton * Right = NULL;                        129   G4Parton * Right = NULL;
133   GetValenceQuarkFlavors(GetDefinition(), Left    130   GetValenceQuarkFlavors(GetDefinition(), Left, Right);
134   Left->SetPosition(GetPosition());               131   Left->SetPosition(GetPosition());
135   Right->SetPosition(GetPosition());              132   Right->SetPosition(GetPosition());
136                                                << 133   
137   G4LorentzVector HadronMom = Get4Momentum();     134   G4LorentzVector HadronMom = Get4Momentum();
                                                   >> 135   //std::cout << "DSU 1 - "<<HadronMom<<std::endl;
138                                                   136 
139   G4double maxAvailMomentum2 = sqr(HadronMom.m << 137   // momenta of string ends 
140                                                << 138   G4double pt2 = HadronMom.perp2();
                                                   >> 139   G4double transverseMass2 = HadronMom.plus()*HadronMom.minus();
                                                   >> 140   G4double maxAvailMomentum2 = sqr(std::sqrt(transverseMass2) - std::sqrt(pt2));
141   G4ThreeVector pt(minTransverseMass, minTrans    141   G4ThreeVector pt(minTransverseMass, minTransverseMass, 0);
142   if (maxAvailMomentum2/widthOfPtSquare>0.01)  << 142   if(maxAvailMomentum2/widthOfPtSquare>0.01) pt = GaussianPt(widthOfPtSquare, maxAvailMomentum2);
                                                   >> 143   //std::cout << "DSU 1.1 - "<< maxAvailMomentum2<< pt <<std::endl;
143                                                   144 
144   G4LorentzVector LeftMom(pt, 0.);                145   G4LorentzVector LeftMom(pt, 0.);
145   G4LorentzVector RightMom;                       146   G4LorentzVector RightMom;
146   RightMom.setPx(HadronMom.px() - pt.x());        147   RightMom.setPx(HadronMom.px() - pt.x());
147   RightMom.setPy(HadronMom.py() - pt.y());        148   RightMom.setPy(HadronMom.py() - pt.y());
                                                   >> 149   //std::cout << "DSU 2 - "<<RightMom<<" "<< LeftMom <<std::endl;
148                                                   150 
149   G4double Local1 = HadronMom.minus() + (Right    151   G4double Local1 = HadronMom.minus() + (RightMom.perp2() - LeftMom.perp2())/HadronMom.plus();
150   G4double Local2 = std::sqrt(std::max(0., sqr    152   G4double Local2 = std::sqrt(std::max(0., sqr(Local1) - 4.*RightMom.perp2()*HadronMom.minus()/HadronMom.plus()));
151                                                << 153   //std::cout << "DSU 3 - "<< Local1 <<" "<< Local2 <<std::endl;
152   if (Direction) Local2 = -Local2;             << 154   if (Direction) Local2 = -Local2;
153   G4double RightMinus   = 0.5*(Local1 + Local2    155   G4double RightMinus   = 0.5*(Local1 + Local2);
154   G4double LeftMinus = HadronMom.minus() - Rig    156   G4double LeftMinus = HadronMom.minus() - RightMinus;
155                                                << 157   //std::cout << "DSU 4 - "<< RightMinus <<" "<< LeftMinus << " "<<HadronMom.minus() <<std::endl;
156   if (LeftMinus <= 0.) {                       << 
157     RightMinus   = 0.5*(Local1 - Local2);      << 
158     LeftMinus = HadronMom.minus() - RightMinus << 
159   }                                            << 
160                                                   158 
161   G4double LeftPlus  = LeftMom.perp2()/LeftMin    159   G4double LeftPlus  = LeftMom.perp2()/LeftMinus;
162   G4double RightPlus = HadronMom.plus() - Left    160   G4double RightPlus = HadronMom.plus() - LeftPlus;
163                                                << 161   //std::cout << "DSU 5 - "<< RightPlus <<" "<< LeftPlus <<std::endl;
164   LeftMom.setPz(0.5*(LeftPlus - LeftMinus));      162   LeftMom.setPz(0.5*(LeftPlus - LeftMinus));
165   LeftMom.setE (0.5*(LeftPlus + LeftMinus));      163   LeftMom.setE (0.5*(LeftPlus + LeftMinus));
166   RightMom.setPz(0.5*(RightPlus - RightMinus))    164   RightMom.setPz(0.5*(RightPlus - RightMinus));
167   RightMom.setE (0.5*(RightPlus + RightMinus))    165   RightMom.setE (0.5*(RightPlus + RightMinus));
168                                                << 166   //std::cout << "DSU 6 - "<< LeftMom <<" "<< RightMom <<std::endl;
169   Left->Set4Momentum(LeftMom);                    167   Left->Set4Momentum(LeftMom);
170   Right->Set4Momentum(RightMom);                  168   Right->Set4Momentum(RightMom);
171                                                << 
172   Color.push_back(Left);                          169   Color.push_back(Left);
173   AntiColor.push_back(Right);                     170   AntiColor.push_back(Right);
174   iP=0; iAP=0;                                 << 
175 }                                                 171 }
176                                                   172 
177                                                   173 
178 void G4QGSMSplitableHadron::SoftSplitUp()         174 void G4QGSMSplitableHadron::SoftSplitUp()
179 {                                                 175 {
180   G4int nSeaPair = GetSoftCollisionCount()-1;  << 176    //... sample transversal momenta for sea and valence quarks
181                                                << 177    G4double phi, pts;
182   G4LorentzVector tmp(0., 0., 0., 0.);         << 178    G4double SumPy = 0.;
183                                                << 179    G4double SumPx = 0.;
184   G4int aSeaPair;                              << 180    G4ThreeVector Pos    = GetPosition();
185   for (aSeaPair = 0; aSeaPair < nSeaPair; aSea << 181    G4int nSeaPair = GetSoftCollisionCount()-1; 
186   {                                            << 182    
187     //  choose quark flavour, d:u:s = 1:1:(1/S << 183    // here the condition,to ensure viability of splitting, also in cases
188     G4int aPDGCode = 1 + (G4int)(G4UniformRand << 184    // where difractive excitation occured together with soft scattering.
189                                                << 185    //   G4double LightConeMomentum = (Direction)? Get4Momentum().plus() : Get4Momentum().minus();
190     //  BuildSeaQuark() determines quark spin, << 186    //   G4double Xmin = theMinPz/LightConeMomentum;
191     //  via parton-constructor G4Parton(aPDGCo << 187    G4double Xmin = theMinPz/( Get4Momentum().e() - GetDefinition()->GetPDGMass() );
192     G4Parton * aParton = BuildSeaQuark(false,  << 188    while(Xmin>=1-(2*nSeaPair+1)*Xmin) Xmin*=0.95;
193                                                << 189 
194     G4int firstPartonColour = aParton->GetColo << 190    G4int aSeaPair;
195     G4double firstPartonSpinZ = aParton->GetSp << 191    for (aSeaPair = 0; aSeaPair < nSeaPair; aSeaPair++)
196                                                << 192    {
197     aParton->Set4Momentum(tmp);                << 193      //  choose quark flavour, d:u:s = 1:1:(1/StrangeSuppress-2)
198     Color.push_back(aParton);                  << 194 
199                                                << 195      G4int aPDGCode = 1 + (G4int)(G4UniformRand()/StrangeSuppress); 
200     // create anti-quark                       << 196 
201     aParton = BuildSeaQuark(true, aPDGCode, nS << 197      //  BuildSeaQuark() determines quark spin, isospin and colour 
202     aParton->SetSpinZ(-firstPartonSpinZ);      << 198      //  via parton-constructor G4Parton(aPDGCode) 
203     aParton->SetColour(-firstPartonColour);    << 199 
204     AntiColor.push_back(aParton);              << 200      G4Parton * aParton = BuildSeaQuark(false, aPDGCode, nSeaPair);
205   }                                            << 201 
206                                                << 202 //    G4cerr << "G4QGSMSplitableHadron::SoftSplitUp()" << G4endl;
207   // Valence quark                             << 203 
208   G4Parton* pColorParton = NULL;               << 204 //    G4cerr << "Parton 1: " 
209   G4Parton* pAntiColorParton = NULL;           << 205 //           << " PDGcode: "  << aPDGCode
210   GetValenceQuarkFlavors(GetDefinition(), pCol << 206 //           << " - Name: "   << aParton->GetDefinition()->GetParticleName()
211                                                << 207 //           << " - Type: "   << aParton->GetDefinition()->GetParticleType() 
212   pColorParton->Set4Momentum(tmp);             << 208 //           << " - Spin-3: " << aParton->GetSpinZ() 
213   pAntiColorParton->Set4Momentum(tmp);         << 209 //           << " - Colour: " << aParton->GetColour() << G4endl;
214                                                << 210 
215   Color.push_back(pColorParton);               << 211      // save colour a spin-3 for anti-quark
216   AntiColor.push_back(pAntiColorParton);       << 212 
217                                                << 213      G4int firstPartonColour = aParton->GetColour();
218   iP=0; iAP=0;                                 << 214      G4double firstPartonSpinZ = aParton->GetSpinZ();
219                                                << 215 
220   return;                                      << 216      SumPx += aParton->Get4Momentum().px();
                                                   >> 217      SumPy += aParton->Get4Momentum().py();
                                                   >> 218      Color.push_back(aParton);
                                                   >> 219 
                                                   >> 220      // create anti-quark
                                                   >> 221 
                                                   >> 222      aParton = BuildSeaQuark(true, aPDGCode, nSeaPair);
                                                   >> 223      aParton->SetSpinZ(-firstPartonSpinZ);
                                                   >> 224      aParton->SetColour(-firstPartonColour);
                                                   >> 225 
                                                   >> 226 //    G4cerr << "Parton 2: " 
                                                   >> 227 //           << " PDGcode: "  << -aPDGCode
                                                   >> 228 //           << " - Name: "   << aParton->GetDefinition()->GetParticleName()
                                                   >> 229 //           << " - Type: "   << aParton->GetDefinition()->GetParticleType() 
                                                   >> 230 //           << " - Spin-3: " << aParton->GetSpinZ() 
                                                   >> 231 //           << " - Colour: " << aParton->GetColour() << G4endl;
                                                   >> 232 //    G4cerr << "------------" << G4endl;
                                                   >> 233 
                                                   >> 234      SumPx += aParton->Get4Momentum().px();
                                                   >> 235      SumPy += aParton->Get4Momentum().py();
                                                   >> 236      AntiColor.push_back(aParton);
                                                   >> 237    }
                                                   >> 238    // Valence quark    
                                                   >> 239    G4Parton* pColorParton = NULL;   
                                                   >> 240    G4Parton* pAntiColorParton = NULL;   
                                                   >> 241    GetValenceQuarkFlavors(GetDefinition(), pColorParton, pAntiColorParton);
                                                   >> 242    G4int ColorEncoding = pColorParton->GetPDGcode();
                                                   >> 243    G4int AntiColorEncoding = pAntiColorParton->GetPDGcode();
                                                   >> 244    
                                                   >> 245    pts   =  sigmaPt*std::sqrt(-std::log(G4UniformRand()));
                                                   >> 246    phi   = 2.*pi*G4UniformRand();
                                                   >> 247    G4double Px = pts*std::cos(phi);
                                                   >> 248    G4double Py = pts*std::sin(phi);
                                                   >> 249    SumPx += Px;
                                                   >> 250    SumPy += Py;
                                                   >> 251 
                                                   >> 252    if (ColorEncoding < 0) // use particle definition
                                                   >> 253    {
                                                   >> 254      G4LorentzVector ColorMom(-SumPx, -SumPy, 0, 0);
                                                   >> 255      pColorParton->Set4Momentum(ColorMom);
                                                   >> 256      G4LorentzVector AntiColorMom(Px, Py, 0, 0);
                                                   >> 257      pAntiColorParton->Set4Momentum(AntiColorMom);
                                                   >> 258    }
                                                   >> 259    else
                                                   >> 260    {
                                                   >> 261      G4LorentzVector ColorMom(Px, Py, 0, 0);
                                                   >> 262      pColorParton->Set4Momentum(ColorMom);
                                                   >> 263      G4LorentzVector AntiColorMom(-SumPx, -SumPy, 0, 0);
                                                   >> 264      pAntiColorParton->Set4Momentum(AntiColorMom);
                                                   >> 265    }
                                                   >> 266    Color.push_back(pColorParton);
                                                   >> 267    AntiColor.push_back(pAntiColorParton);
                                                   >> 268 
                                                   >> 269    // Sample X
                                                   >> 270    G4int nAttempt = 0;
                                                   >> 271    G4double SumX = 0;
                                                   >> 272    G4double aBeta = beta;
                                                   >> 273    G4double ColorX, AntiColorX;
                                                   >> 274    G4double HPWtest = 0;
                                                   >> 275    if (GetDefinition() == G4PionMinus::PionMinusDefinition()) aBeta = 1.;        
                                                   >> 276    if (GetDefinition() == G4Gamma::GammaDefinition()) aBeta = 1.;        
                                                   >> 277    if (GetDefinition() == G4PionPlus::PionPlusDefinition()) aBeta = 1.;     
                                                   >> 278    if (GetDefinition() == G4PionZero::PionZeroDefinition()) aBeta = 1.;       
                                                   >> 279    if (GetDefinition() == G4KaonPlus::KaonPlusDefinition()) aBeta = 0.;       
                                                   >> 280    if (GetDefinition() == G4KaonMinus::KaonMinusDefinition()) aBeta = 0.;       
                                                   >> 281    do
                                                   >> 282    {
                                                   >> 283      SumX = 0;
                                                   >> 284      nAttempt++;
                                                   >> 285      G4int    NumberOfUnsampledSeaQuarks = 2*nSeaPair;
                                                   >> 286      G4double beta1 = beta;
                                                   >> 287      if (std::abs(ColorEncoding) <= 1000 && std::abs(AntiColorEncoding) <= 1000) beta1 = 1.; //...  in a meson        
                                                   >> 288      ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
                                                   >> 289      HPWtest = ColorX;
                                                   >> 290      while (ColorX < Xmin || ColorX > 1.|| 1. -  ColorX <= Xmin); 
                                                   >> 291      Color.back()->SetX(SumX = ColorX);// this is the valenz quark.
                                                   >> 292      for(G4int aPair = 0; aPair < nSeaPair; aPair++) 
                                                   >> 293      {
                                                   >> 294        NumberOfUnsampledSeaQuarks--;
                                                   >> 295        ColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
                                                   >> 296        Color[aPair]->SetX(ColorX);
                                                   >> 297        SumX += ColorX; 
                                                   >> 298        NumberOfUnsampledSeaQuarks--;
                                                   >> 299        AntiColorX = SampleX(Xmin, NumberOfUnsampledSeaQuarks, 2*nSeaPair, aBeta);
                                                   >> 300        AntiColor[aPair]->SetX(AntiColorX); // the 'sea' partons
                                                   >> 301        SumX += AntiColorX;
                                                   >> 302        if (1. - SumX <= Xmin)  break;
                                                   >> 303      }
                                                   >> 304    } 
                                                   >> 305    while (1. - SumX <= Xmin); 
                                                   >> 306    (*(AntiColor.end()-1))->SetX(1. - SumX); // the di-quark takes the rest, then go to momentum
                                                   >> 307    /// and here is the bug ;-) @@@@@@@@@@@@@
                                                   >> 308    if(getenv("debug_QGSMSplitableHadron") )G4cout << "particle energy at split = "<<Get4Momentum().t()<<G4endl;
                                                   >> 309    G4double lightCone = ((!Direction) ? Get4Momentum().minus() : Get4Momentum().plus());
                                                   >> 310    // lightCone -= 0.5*Get4Momentum().m();
                                                   >> 311    // hpw testing @@@@@ lightCone = 2.*Get4Momentum().t();
                                                   >> 312    if(getenv("debug_QGSMSplitableHadron") )G4cout << "Light cone = "<<lightCone<<G4endl;
                                                   >> 313    for(aSeaPair = 0; aSeaPair < nSeaPair+1; aSeaPair++) 
                                                   >> 314    {
                                                   >> 315      G4Parton* aParton = Color[aSeaPair];
                                                   >> 316      aParton->DefineMomentumInZ(lightCone, Direction);
                                                   >> 317 
                                                   >> 318      aParton = AntiColor[aSeaPair]; 
                                                   >> 319      aParton->DefineMomentumInZ(lightCone, Direction);
                                                   >> 320    }  
                                                   >> 321 //--DEBUG--   cout <<G4endl<<"XSAMPLE "<<HPWtest<<G4endl;
                                                   >> 322    return;
221 }                                                 323 } 
222                                                   324 
223                                                << 325 void G4QGSMSplitableHadron::GetValenceQuarkFlavors(const G4ParticleDefinition * aPart, G4Parton *& Parton1, G4Parton *& Parton2)
224 void G4QGSMSplitableHadron::GetValenceQuarkFla << 
225                                                << 
226 {                                                 326 {
227   // Note! convention aEnd = q or (qq)bar and  << 327    // Note! convention aEnd = q or (qq)bar and bEnd = qbar or qq.
228   G4int aEnd=0;                                << 328   G4int aEnd;
229   G4int bEnd=0;                                << 329   G4int bEnd;
230   G4int HadronEncoding = aPart->GetPDGEncoding    330   G4int HadronEncoding = aPart->GetPDGEncoding();
231   if (aPart->GetBaryonNumber() == 0)           << 331   if (aPart->GetBaryonNumber() == 0)  
232   {                                               332   {
233     theMesonSplitter.SplitMeson(HadronEncoding << 333     theMesonSplitter.SplitMeson(HadronEncoding, &aEnd, &bEnd);
234   }                                               334   }
235   else                                            335   else
236   {                                               336   {
237     theBaryonSplitter.SplitBarion(HadronEncodi << 337     theBaryonSplitter.SplitBarion(HadronEncoding, &aEnd, &bEnd);
238   }                                            << 338   } 
239                                                   339 
240   Parton1 = new G4Parton(aEnd);                   340   Parton1 = new G4Parton(aEnd);
241   Parton1->SetPosition(GetPosition());            341   Parton1->SetPosition(GetPosition());
242                                                   342 
                                                   >> 343 //  G4cerr << "G4QGSMSplitableHadron::GetValenceQuarkFlavors()" << G4endl;
                                                   >> 344 //  G4cerr << "Parton 1: " 
                                                   >> 345 //         << " PDGcode: "  << aEnd
                                                   >> 346 //         << " - Name: "   << Parton1->GetDefinition()->GetParticleName()
                                                   >> 347 //         << " - Type: "   << Parton1->GetDefinition()->GetParticleType() 
                                                   >> 348 //         << " - Spin-3: " << Parton1->GetSpinZ() 
                                                   >> 349 //         << " - Colour: " << Parton1->GetColour() << G4endl;
                                                   >> 350 
243   Parton2 = new G4Parton(bEnd);                   351   Parton2 = new G4Parton(bEnd);
244   Parton2->SetPosition(GetPosition());            352   Parton2->SetPosition(GetPosition());
245                                                   353 
                                                   >> 354 //  G4cerr << "Parton 2: " 
                                                   >> 355 //         << " PDGcode: "  << bEnd
                                                   >> 356 //         << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
                                                   >> 357 //         << " - Type: "   << Parton2->GetDefinition()->GetParticleType() 
                                                   >> 358 //         << " - Spin-3: " << Parton2->GetSpinZ() 
                                                   >> 359 //         << " - Colour: " << Parton2->GetColour() << G4endl;
                                                   >> 360 //  G4cerr << "... now checking for color and spin conservation - yielding: " << G4endl;
                                                   >> 361 
246   // colour of parton 1 choosen at random by G    362   // colour of parton 1 choosen at random by G4Parton(aEnd)
247   // colour of parton 2 is the opposite:          363   // colour of parton 2 is the opposite:
248                                                   364 
249   Parton2->SetColour(-(Parton1->GetColour()));    365   Parton2->SetColour(-(Parton1->GetColour()));
250                                                   366 
251   // isospin-3 of both partons is handled by G << 367   // isospin-3 of both partons is handled by G4Parton(PDGCode)
252                                                   368 
253   // spin-3 of parton 1 and 2 choosen at rando << 369   // spin-3 of parton 1 and 2 choosen at random by G4Parton(aEnd)
254   // spin-3 of parton 2 may be constrained by  << 370   // spin-3 of parton 2 may be constrained by spin of original particle:
255                                                   371 
256   if ( std::abs(Parton1->GetSpinZ() + Parton2- << 372   if ( std::abs(Parton1->GetSpinZ() + Parton2->GetSpinZ()) > aPart->GetPDGSpin()) 
257   {                                               373   {
258     Parton2->SetSpinZ(-(Parton2->GetSpinZ())); << 374     Parton2->SetSpinZ(-(Parton2->GetSpinZ()));    
259   }                                            << 375   } 
260 }                                              << 
261                                                   376 
                                                   >> 377 //  G4cerr << "Parton 2: " 
                                                   >> 378 //         << " PDGcode: "  << bEnd
                                                   >> 379 //         << " - Name: "   << Parton2->GetDefinition()->GetParticleName()
                                                   >> 380 //         << " - Type: "   << Parton2->GetDefinition()->GetParticleType() 
                                                   >> 381 //         << " - Spin-3: " << Parton2->GetSpinZ() 
                                                   >> 382 //         << " - Colour: " << Parton2->GetColour() << G4endl;
                                                   >> 383 //  G4cerr << "------------" << G4endl;
262                                                   384 
                                                   >> 385 }
                                                   >> 386  
                                                   >> 387     
263 G4ThreeVector G4QGSMSplitableHadron::GaussianP    388 G4ThreeVector G4QGSMSplitableHadron::GaussianPt(G4double widthSquare, G4double maxPtSquare)
264 {                                                 389 {
265   G4double R;                                     390   G4double R;
266   const G4int maxNumberOfLoops = 1000;         << 391   while((R = -widthSquare*std::log(G4UniformRand())) > maxPtSquare);
267   G4int loopCounter = -1;                      << 
268   while( ((R = -widthSquare*G4Log(G4UniformRan << 
269          ++loopCounter < maxNumberOfLoops ) {; << 
270   if ( loopCounter >= maxNumberOfLoops ) {     << 
271     R = 0.99*maxPtSquare;  // Just an acceptab << 
272   }                                            << 
273   R = std::sqrt(R);                               392   R = std::sqrt(R);
274   G4double phi = twopi*G4UniformRand();           393   G4double phi = twopi*G4UniformRand();
275   return G4ThreeVector (R*std::cos(phi), R*std << 394   return G4ThreeVector (R*std::cos(phi), R*std::sin(phi), 0.);    
276 }                                                 395 }
277                                                   396 
278                                                << 
279 G4Parton * G4QGSMSplitableHadron::                397 G4Parton * G4QGSMSplitableHadron::
280 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCo    398 BuildSeaQuark(G4bool isAntiQuark, G4int aPDGCode, G4int /* nSeaPair*/)
281 {                                                 399 {
282   if (isAntiQuark) aPDGCode*=-1;                  400   if (isAntiQuark) aPDGCode*=-1;
283   G4Parton* result = new G4Parton(aPDGCode);   << 401   G4Parton* result = new G4Parton(aPDGCode);   
284   result->SetPosition(GetPosition());             402   result->SetPosition(GetPosition());
285   G4ThreeVector aPtVector = GaussianPt(sigmaPt    403   G4ThreeVector aPtVector = GaussianPt(sigmaPt, DBL_MAX);
286   G4LorentzVector a4Momentum(aPtVector, 0);       404   G4LorentzVector a4Momentum(aPtVector, 0);
287   result->Set4Momentum(a4Momentum);               405   result->Set4Momentum(a4Momentum);
288   return result;                                  406   return result;
289 }                                                 407 }
290                                                   408 
291                                                << 
292 G4double G4QGSMSplitableHadron::                  409 G4double G4QGSMSplitableHadron::
293 SampleX(G4double anXmin, G4int nSea, G4int tot    410 SampleX(G4double anXmin, G4int nSea, G4int totalSea, G4double aBeta)
294 {                                                 411 {
295   G4double result;                                412   G4double result;
296   G4double x1, x2;                                413   G4double x1, x2;
297   G4double ymax = 0;                              414   G4double ymax = 0;
298   for(G4int ii=1; ii<100; ii++)                << 415   for(G4int ii=0; ii<100; ii++)
299   {                                               416   {
300     G4double y = G4Pow::GetInstance()->powA(1. << 417     G4double y = std::pow(1./G4double(ii), alpha);
301     y *= G4Pow::GetInstance()->powN( G4Pow::Ge << 418     y *= std::pow( std::pow(1-anXmin-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
302                                      G4Pow::Ge << 419     y *= std::pow(1-anXmin-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);
303     y *= G4Pow::GetInstance()->powA(1-anXmin-t << 420     if(y>ymax) ymax = y;
304          G4Pow::GetInstance()->powA(anXmin, aB << 
305    if (y>ymax) ymax = y;                       << 
306   }                                               421   }
307   G4double y;                                     422   G4double y;
308   G4double xMax=1-(totalSea+1)*anXmin;            423   G4double xMax=1-(totalSea+1)*anXmin;
309   if (anXmin > xMax)                           << 424   if(anXmin > xMax) 
310   {                                               425   {
311     throw G4HadronicException(__FILE__, __LINE << 426     G4cout << "anXmin = "<<anXmin<<" nSea = "<<nSea<<" totalSea = "<< totalSea<<G4endl;
312             "G4QGSMSplitableHadron - Fatal: Ca << 427     throw G4HadronicException(__FILE__, __LINE__, "G4QGSMSplitableHadron - Fatal: Cannot sample parton densities under these constraints.");
313   }                                               428   }
314   const G4int maxNumberOfLoops = 1000;         << 
315   G4int loopCounter = 0;                       << 
316   do                                              429   do
317   {                                               430   {
318     x1 = G4RandFlat::shoot(anXmin, xMax);      << 431     x1 = CLHEP::RandFlat::shoot(anXmin, xMax);
319     y = G4Pow::GetInstance()->powA(x1, alpha); << 432     y = std::pow(x1, alpha);
320     y *= G4Pow::GetInstance()->powN( G4Pow::Ge << 433     y *= std::pow( std::pow(1-x1-totalSea*anXmin, alpha+1) - std::pow(anXmin, alpha+1), nSea);
321                                      G4Pow::Ge << 434     y *= std::pow(1-x1-totalSea*anXmin, aBeta+1) - std::pow(anXmin, aBeta+1);  
322     y *= G4Pow::GetInstance()->powA(1-x1-total << 
323          G4Pow::GetInstance()->powA(anXmin, aB << 
324     x2 = ymax*G4UniformRand();                    435     x2 = ymax*G4UniformRand();
325   } while( (x2>y) && ++loopCounter < maxNumber << 
326   if ( loopCounter >= maxNumberOfLoops ) {     << 
327     x1 = 0.5*( anXmin + xMax );  // Just an ac << 
328   }                                               436   }
                                                   >> 437   while(x2>y);
329   result = x1;                                    438   result = x1;
330   return result;                               << 439   return result;  
331 }                                                 440 }
332                                                   441