Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.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/hadronization/src/G4QGSMFragmentation.cc (Version 11.3.0) and /processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.cc (Version 11.1.3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 //                                                 27 //
 28 // -------------------------------------------     28 // -----------------------------------------------------------------------------
 29 //      GEANT 4 class implementation file          29 //      GEANT 4 class implementation file
 30 //                                                 30 //
 31 //      History: first implementation, Maxim K     31 //      History: first implementation, Maxim Komogorov, 10-Jul-1998
 32 // -------------------------------------------     32 // -----------------------------------------------------------------------------
 33 #include "G4QGSMFragmentation.hh"                  33 #include "G4QGSMFragmentation.hh"
 34 #include "G4PhysicalConstants.hh"                  34 #include "G4PhysicalConstants.hh"
 35 #include "G4SystemOfUnits.hh"                      35 #include "G4SystemOfUnits.hh"
 36 #include "Randomize.hh"                            36 #include "Randomize.hh"
 37 #include "G4ios.hh"                                37 #include "G4ios.hh"
 38 #include "G4FragmentingString.hh"                  38 #include "G4FragmentingString.hh"
 39 #include "G4DiQuarks.hh"                           39 #include "G4DiQuarks.hh"
 40 #include "G4Quarks.hh"                             40 #include "G4Quarks.hh"
 41 #include "G4HadronicParameters.hh"                 41 #include "G4HadronicParameters.hh"
 42 #include "G4Pow.hh"                                42 #include "G4Pow.hh"
 43                                                    43 
 44 //#define debug_QGSMfragmentation                  44 //#define debug_QGSMfragmentation 
 45                                                    45 
 46 // Class G4QGSMFragmentation                       46 // Class G4QGSMFragmentation 
 47 //********************************************     47 //****************************************************************************************
 48                                                    48  
 49 G4QGSMFragmentation::G4QGSMFragmentation()         49 G4QGSMFragmentation::G4QGSMFragmentation() 
 50 {                                                  50 {
 51     SigmaQT = 0.45 * GeV;                          51     SigmaQT = 0.45 * GeV;
 52                                                    52 
 53     MassCut = 0.35*GeV;                            53     MassCut = 0.35*GeV; 
 54                                                    54 
 55     SetStrangenessSuppression((1.0 - 0.16)/2.)     55     SetStrangenessSuppression((1.0 - 0.16)/2.);
 56                                                    56 
 57     // Check if charmed and bottom hadrons are     57     // Check if charmed and bottom hadrons are enabled: if this is the case, then
 58     // set the non-zero probabilities for c-cb     58     // set the non-zero probabilities for c-cbar and b-bbar creation from the vacuum,
 59     // else set them to 0.0. If these probabil     59     // else set them to 0.0. If these probabilities are/aren't zero then charmed or bottom
 60     // hadrons can't/can be created during the     60     // hadrons can't/can be created during the string fragmentation of ordinary
 61     // (i.e. not heavy) projectile hadron nucl     61     // (i.e. not heavy) projectile hadron nuclear reactions.
 62     if ( G4HadronicParameters::Instance()->Ena     62     if ( G4HadronicParameters::Instance()->EnableBCParticles() ) {
 63       SetProbCCbar(0.0002);  // According to O     63       SetProbCCbar(0.0002);  // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094; tuned by Uzhi Oct. 2022
 64       SetProbBBbar(5.0e-5);  // According to O     64       SetProbBBbar(5.0e-5);  // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
 65     } else {                                       65     } else {
 66       SetProbCCbar(0.0);                           66       SetProbCCbar(0.0);
 67       SetProbBBbar(0.0);                           67       SetProbBBbar(0.0);
 68     }                                              68     }
 69                                                    69     
 70     SetDiquarkSuppression(0.32);                   70     SetDiquarkSuppression(0.32);
 71     SetDiquarkBreakProbability(0.7);               71     SetDiquarkBreakProbability(0.7);
 72                                                    72 
 73     SetMinMasses();                                73     SetMinMasses();
 74                                                    74 
 75     arho = 0.5;    // alpha_rho0                   75     arho = 0.5;    // alpha_rho0
 76     aphi = 0.0;    // alpha_fi                     76     aphi = 0.0;    // alpha_fi
 77     aJPs =-2.2;    // alpha_J/Psi                  77     aJPs =-2.2;    // alpha_J/Psi
 78     aUps =-8.0;    // alpha_Y      ??? O. Pisk     78     aUps =-8.0;    // alpha_Y      ??? O. Piskunova Yad. Phys. 56 (1993) 1094.
 79                                                    79 
 80     aksi =-1.0;                                    80     aksi =-1.0; 
 81     alft = 0.5;    // 2 * alpha'_R *<Pt^2>         81     alft = 0.5;    // 2 * alpha'_R *<Pt^2>
 82                                                    82 
 83     an    = -0.5 ;                                 83     an    = -0.5 ; 
 84     ala   = -0.75; // an - arho/2 + aphi/2         84     ala   = -0.75; // an - arho/2 + aphi/2 
 85     alaC  =  an - arho/2.0 + aJPs/2.0;             85     alaC  =  an - arho/2.0 + aJPs/2.0;  
 86     alaB  =  an - arho/2.0 + aUps/2.0;             86     alaB  =  an - arho/2.0 + aUps/2.0;  
 87     aXi   =  0.0;  // ??                           87     aXi   =  0.0;  // ??
 88     aXiC  =  0.0;  // ??                           88     aXiC  =  0.0;  // ??
 89     aXiB  =  0.0;  // ??                           89     aXiB  =  0.0;  // ??
 90     aXiCC =  0.0;  // ??                           90     aXiCC =  0.0;  // ??
 91     aXiCB =  0.0;  // ??                           91     aXiCB =  0.0;  // ??
 92     aXiBB =  0.0;  // ??                           92     aXiBB =  0.0;  // ??
 93                                                    93 
 94     SetFFq2q();                                    94     SetFFq2q();
 95     SetFFq2qq();                                   95     SetFFq2qq();
 96     SetFFqq2q();                                   96     SetFFqq2q();
 97     SetFFqq2qq();                                  97     SetFFqq2qq();
 98                          // d  u   s   c   b       98                          // d  u   s   c   b
 99     G4int Index[5][5] = { { 0, 1,  2,  3,  4 }     99     G4int Index[5][5] = { { 0, 1,  2,  3,  4 },    // d
100                           { 1, 5,  6,  7,  8 }    100                           { 1, 5,  6,  7,  8 },    // u
101                           { 2, 6,  9, 10, 11 }    101                           { 2, 6,  9, 10, 11 },    // s
102                           { 3, 7, 10, 12, 13 }    102                           { 3, 7, 10, 12, 13 },    // c
103                           { 4, 8, 11, 13, 14 }    103                           { 4, 8, 11, 13, 14 } };  // b
104     for (G4int i = 0; i < 5; i++ ) {              104     for (G4int i = 0; i < 5; i++ ) {
105       for ( G4int j = 0; j < 5; j++ ) {           105       for ( G4int j = 0; j < 5; j++ ) { 
106         IndexDiQ[i][j] = Index[i][j];             106         IndexDiQ[i][j] = Index[i][j];
107       }                                           107       } 
108     };                                            108     };
109 }                                                 109 }
110                                                   110 
111 G4QGSMFragmentation::~G4QGSMFragmentation()       111 G4QGSMFragmentation::~G4QGSMFragmentation()
112 {}                                                112 {}
113                                                   113 
114 //--------------------------------------------    114 //----------------------------------------------------------------------------------------------------------
115                                                   115 
116 G4KineticTrackVector* G4QGSMFragmentation::Fra    116 G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
117 {                                                 117 {
118                                                   118 
119   G4FragmentingString  aString(theString);        119   G4FragmentingString  aString(theString);
120   SetMinimalStringMass(&aString);                 120   SetMinimalStringMass(&aString);
121                                                   121 
122   #ifdef debug_QGSMfragmentation                  122   #ifdef debug_QGSMfragmentation
123   G4cout<<G4endl<<"QGSM StringFragm: String Ma    123   G4cout<<G4endl<<"QGSM StringFragm: String Mass "
124                              <<theString.Get4M    124                              <<theString.Get4Momentum().mag()<<" Pz "
125                              <<theString.Get4M    125                              <<theString.Get4Momentum().pz()
126                              <<"--------------    126                              <<"------------------------------------"<<G4endl;
127   G4cout<<"String ends Direct "<<theString.Get    127   G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
128                                <<theString.Get    128                                <<theString.GetRightParton()->GetPDGcode()<<" "
129                                <<theString.Get    129                                <<theString.GetDirection()<< G4endl;
130   G4cout<<"Left  mom "<<theString.GetLeftParto    130   G4cout<<"Left  mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
131   G4cout<<"Right mom "<<theString.GetRightPart    131   G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl;
132   G4cout<<"Check for Fragmentation "<<G4endl;     132   G4cout<<"Check for Fragmentation "<<G4endl;
133   #endif                                          133   #endif
134                                                   134 
135   // Can no longer modify Parameters for Fragm    135   // Can no longer modify Parameters for Fragmentation.
136   PastInitPhase=true;                             136   PastInitPhase=true;
137                                                   137   
138   // Check if string has enough mass to fragme    138   // Check if string has enough mass to fragment...
139   G4KineticTrackVector * LeftVector=NULL;         139   G4KineticTrackVector * LeftVector=NULL;
140                                                   140 
141   if ( !IsItFragmentable(&aString) ) {            141   if ( !IsItFragmentable(&aString) ) {
142      LeftVector=ProduceOneHadron(&theString);     142      LeftVector=ProduceOneHadron(&theString);
143                                                   143 
144      #ifdef debug_QGSMfragmentation               144      #ifdef debug_QGSMfragmentation
145      if ( LeftVector != 0 ) G4cout<<"Non fragm    145      if ( LeftVector != 0 ) G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
146      #endif                                       146      #endif
147                                                   147 
148      if ( LeftVector == nullptr ) LeftVector =    148      if ( LeftVector == nullptr ) LeftVector = new G4KineticTrackVector;
149      return LeftVector;                           149      return LeftVector;
150   }                                               150   }
151                                                   151 
152   #ifdef debug_QGSMfragmentation                  152   #ifdef debug_QGSMfragmentation
153   G4cout<<"The string will be fragmented. "<<G    153   G4cout<<"The string will be fragmented. "<<G4endl;
154   #endif                                          154   #endif
155                                                   155   
156   LeftVector = new G4KineticTrackVector;          156   LeftVector = new G4KineticTrackVector;
157   G4KineticTrackVector * RightVector=new G4Kin    157   G4KineticTrackVector * RightVector=new G4KineticTrackVector;
158                                                   158 
159   G4ExcitedString *theStringInCMS=CopyExcited(    159   G4ExcitedString *theStringInCMS=CopyExcited(theString);
160   G4LorentzRotation toCms=theStringInCMS->Tran    160   G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
161                                                   161 
162   G4bool success=false, inner_sucess=true;        162   G4bool success=false, inner_sucess=true;
163   G4int attempt=0;                                163   G4int attempt=0;
164   while ( !success && attempt++ < StringLoopIn    164   while ( !success && attempt++ < StringLoopInterrupt )  /* Loop checking, 07.08.2015, A.Ribon */
165   {                                               165   {
166                 #ifdef debug_QGSMfragmentation    166                 #ifdef debug_QGSMfragmentation
167                 G4cout<<"Loop_toFrag "<<theStr    167                 G4cout<<"Loop_toFrag "<<theStringInCMS->GetLeftParton()->GetPDGcode()<<" "
168                                       <<theStr    168                                       <<theStringInCMS->GetRightParton()->GetPDGcode()<<" "
169                                       <<theStr    169                                       <<theStringInCMS->GetDirection()<< G4endl;
170                 #endif                            170                 #endif
171                                                   171 
172     G4FragmentingString *currentString=new G4F    172     G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
173                                                   173 
174     std::for_each(LeftVector->begin(), LeftVec    174     std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
175     LeftVector->clear();                          175     LeftVector->clear();
176     std::for_each(RightVector->begin(), RightV    176     std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
177     RightVector->clear();                         177     RightVector->clear();
178                                                   178     
179     inner_sucess=true;  // set false on failur    179     inner_sucess=true;  // set false on failure..
180                 const G4int maxNumberOfLoops =    180                 const G4int maxNumberOfLoops = 1000;
181                 G4int loopCounter = -1;           181                 G4int loopCounter = -1;
182     while (! StopFragmenting(currentString) &&    182     while (! StopFragmenting(currentString) && ++loopCounter < maxNumberOfLoops )   /* Loop checking, 07.08.2015, A.Ribon */
183     {  // Split current string into hadron + n    183     {  // Split current string into hadron + new string
184                                                   184 
185                         #ifdef debug_QGSMfragm    185                         #ifdef debug_QGSMfragmentation
186                         G4cout<<"The string ca    186                         G4cout<<"The string can fragment. "<<G4endl;;
187                         #endif                    187                         #endif
188       G4FragmentingString *newString=0;  // us    188       G4FragmentingString *newString=0;  // used as output from SplitUp...
189       G4KineticTrack * Hadron=Splitup(currentS    189       G4KineticTrack * Hadron=Splitup(currentString,newString);
190                                                   190 
191       if ( Hadron != 0 )                          191       if ( Hadron != 0 ) 
192       {                                           192       {
193                            #ifdef debug_QGSMfr    193                            #ifdef debug_QGSMfragmentation
194                            G4cout<<"Hadron pro    194                            G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
195                            #endif                 195                            #endif
196                            // To close the pro    196                            // To close the production of hadrons at fragmentation stage
197          if ( currentString->GetDecayDirection    197          if ( currentString->GetDecayDirection() > 0 )
198            LeftVector->push_back(Hadron);         198            LeftVector->push_back(Hadron);
199                else                               199                else
200              RightVector->push_back(Hadron);      200              RightVector->push_back(Hadron);
201                                                   201 
202          delete currentString;                    202          delete currentString;
203          currentString=newString;                 203          currentString=newString;
204                                                   204 
205       } else {                                    205       } else {
206                                                   206 
207                            #ifdef debug_QGSMfr    207                            #ifdef debug_QGSMfragmentation
208                            G4cout<<"abandon ..    208                            G4cout<<"abandon ... start from the beginning ---------------"<<G4endl;
209                            #endif                 209                            #endif
210                                                   210 
211          // Abandon ... start from the beginni    211          // Abandon ... start from the beginning
212          if (newString) delete newString;         212          if (newString) delete newString;
213          inner_sucess=false;                      213          inner_sucess=false;
214          break;                                   214          break;
215       }                                           215       }
216     }                                             216     }
217                 if ( loopCounter >= maxNumberO    217                 if ( loopCounter >= maxNumberOfLoops ) {
218                   inner_sucess=false;             218                   inner_sucess=false;
219                 }                                 219                 }
220                                                   220 
221     // Split current string into 2 final Hadro    221     // Split current string into 2 final Hadrons
222                 #ifdef debug_QGSMfragmentation    222                 #ifdef debug_QGSMfragmentation
223                 if( inner_sucess ) {              223                 if( inner_sucess ) {
224                   G4cout<<"Split remaining str    224                   G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
225                 } else {                          225                 } else {
226       G4cout<<" New attempt to fragment string    226       G4cout<<" New attempt to fragment string"<<G4endl;
227     }                                             227     }
228                 #endif                            228                 #endif
229                 // To the close production of     229                 // To the close production of hadrons at last string decay
230     if ( inner_sucess &&                          230     if ( inner_sucess && 
231          SplitLast(currentString,LeftVector, R    231          SplitLast(currentString,LeftVector, RightVector) ) 
232     {                                             232     {
233       success=true;                               233       success=true;
234     }                                             234     }
235     delete currentString;                         235     delete currentString;
236   }                                               236   }
237                                                   237   
238   delete theStringInCMS;                          238   delete theStringInCMS;
239                                                   239   
240   if ( ! success )                                240   if ( ! success )
241   {                                               241   {
242     std::for_each(LeftVector->begin(), LeftVec    242     std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
243     LeftVector->clear();                          243     LeftVector->clear();
244     std::for_each(RightVector->begin(), RightV    244     std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
245     delete RightVector;                           245     delete RightVector;
246     return LeftVector;                            246     return LeftVector;
247   }                                               247   }
248                                                   248     
249   // Join Left- and RightVector into LeftVecto    249   // Join Left- and RightVector into LeftVector in correct order.
250   while(!RightVector->empty())  /* Loop checki    250   while(!RightVector->empty())  /* Loop checking, 07.08.2015, A.Ribon */
251   {                                               251   {
252       LeftVector->push_back(RightVector->back(    252       LeftVector->push_back(RightVector->back());
253       RightVector->erase(RightVector->end()-1)    253       RightVector->erase(RightVector->end()-1);
254   }                                               254   }
255   delete RightVector;                             255   delete RightVector;
256                                                   256 
257   CalculateHadronTimePosition(theString.Get4Mo    257   CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
258                                                   258 
259   G4LorentzRotation toObserverFrame(toCms.inve    259   G4LorentzRotation toObserverFrame(toCms.inverse());
260                                                   260 
261   for (size_t C1 = 0; C1 < LeftVector->size();    261   for (size_t C1 = 0; C1 < LeftVector->size(); C1++)
262   {                                               262   {
263      G4KineticTrack* Hadron = LeftVector->oper    263      G4KineticTrack* Hadron = LeftVector->operator[](C1);
264      G4LorentzVector Momentum = Hadron->Get4Mo    264      G4LorentzVector Momentum = Hadron->Get4Momentum();
265      Momentum = toObserverFrame*Momentum;         265      Momentum = toObserverFrame*Momentum;
266      Hadron->Set4Momentum(Momentum);              266      Hadron->Set4Momentum(Momentum);
267      G4LorentzVector Coordinate(Hadron->GetPos    267      G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
268      Momentum = toObserverFrame*Coordinate;       268      Momentum = toObserverFrame*Coordinate;
269      Hadron->SetFormationTime(Momentum.e());      269      Hadron->SetFormationTime(Momentum.e());
270      G4ThreeVector aPosition(Momentum.vect());    270      G4ThreeVector aPosition(Momentum.vect());
271      Hadron->SetPosition(theString.GetPosition    271      Hadron->SetPosition(theString.GetPosition()+aPosition);
272   }                                               272   }
273   return LeftVector;                              273   return LeftVector;
274 }                                                 274 }
275                                                   275 
276 //--------------------------------------------    276 //----------------------------------------------------------------------------------------------------------
277                                                   277 
278 G4bool G4QGSMFragmentation::IsItFragmentable(c    278 G4bool G4QGSMFragmentation::IsItFragmentable(const G4FragmentingString * string)
279 {                                                 279 {
280   return sqr( MinimalStringMass + MassCut ) <     280   return sqr( MinimalStringMass + MassCut ) < string->Mass2();
281 }                                                 281 }
282                                                   282 
283 //--------------------------------------------    283 //----------------------------------------------------------------------------------------------------------
284                                                   284 
285 G4bool G4QGSMFragmentation::StopFragmenting(co    285 G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * string)
286 {                                                 286 {
287   SetMinimalStringMass(string);                   287   SetMinimalStringMass(string);
288         if ( MinimalStringMass < 0.0 ) return     288         if ( MinimalStringMass < 0.0 ) return true;
289                                                   289 
290         G4double smass = string->Mass();          290         G4double smass = string->Mass();
291   G4double x = (string->IsAFourQuarkString())     291   G4double x = (string->IsAFourQuarkString()) ? 0.005*(smass - MinimalStringMass)
292     : 0.66e-6*(smass - MinimalStringMass)*(sma    292     : 0.66e-6*(smass - MinimalStringMass)*(smass + MinimalStringMass);
293                                                   293 
294         G4bool res = true;                        294         G4bool res = true;
295         if(x > 0.0) {                             295         if(x > 0.0) {
296           res = (x < 200.) ? (G4UniformRand()     296           res = (x < 200.) ? (G4UniformRand() < G4Exp(-x)) : false;
297   }                                               297   }
298   return res;                                     298   return res;
299 }                                                 299 }
300                                                   300 
301 //--------------------------------------------    301 //-----------------------------------------------------------------------------
302                                                   302 
303 G4KineticTrack * G4QGSMFragmentation::Splitup(    303 G4KineticTrack * G4QGSMFragmentation::Splitup( G4FragmentingString *string, 
304                              G4FragmentingStri    304                              G4FragmentingString *&newString )
305 {                                                 305 {
306        #ifdef debug_QGSMfragmentation             306        #ifdef debug_QGSMfragmentation
307        G4cout<<G4endl;                            307        G4cout<<G4endl;
308        G4cout<<"Start SplitUP (G4VLongitudinal    308        G4cout<<"Start SplitUP (G4VLongitudinalStringDecay) ========================="<<G4endl;
309        G4cout<<"String partons: " <<string->Ge    309        G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
310                                   <<string->Ge    310                                   <<string->GetRightParton()->GetPDGEncoding()<<" "
311              <<"Direction "       <<string->Ge    311              <<"Direction "       <<string->GetDecayDirection()<<G4endl;
312        #endif                                     312        #endif
313                                                   313 
314        //... random choice of string end to us    314        //... random choice of string end to use for creating the hadron (decay)   
315        G4int SideOfDecay = (G4UniformRand() <     315        G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
316        if (SideOfDecay < 0)                       316        if (SideOfDecay < 0)
317        {                                          317        {
318     string->SetLeftPartonStable();                318     string->SetLeftPartonStable();
319        } else                                     319        } else
320        {                                          320        {
321           string->SetRightPartonStable();         321           string->SetRightPartonStable();
322        }                                          322        }
323                                                   323 
324        G4ParticleDefinition *newStringEnd;        324        G4ParticleDefinition *newStringEnd;
325        G4ParticleDefinition * HadronDefinition    325        G4ParticleDefinition * HadronDefinition;
326        if (string->DecayIsQuark())                326        if (string->DecayIsQuark())
327        {                                          327        {
328     G4double ProbDqADq = GetDiquarkSuppress();    328     G4double ProbDqADq = GetDiquarkSuppress();
329                                                   329 
330     G4int NumberOfpossibleBaryons = 2;            330     G4int NumberOfpossibleBaryons = 2;
331                                                   331 
332     if (string->GetLeftParton()->GetParticleSu    332     if (string->GetLeftParton()->GetParticleSubType()  != "quark") NumberOfpossibleBaryons++; 
333     if (string->GetRightParton()->GetParticleS    333     if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
334                                                   334 
335     G4double ActualProb  = ProbDqADq ;            335     G4double ActualProb  = ProbDqADq ;
336           ActualProb *= (1.0-G4Exp(2.0*(1.0 -     336           ActualProb *= (1.0-G4Exp(2.0*(1.0 - string->Mass()/(NumberOfpossibleBaryons*1400.0))));
337                                                   337 
338     SetDiquarkSuppression(ActualProb);            338     SetDiquarkSuppression(ActualProb); 
339                                                   339 
340           HadronDefinition= QuarkSplitup(strin    340           HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
341                                                   341 
342     SetDiquarkSuppression(ProbDqADq);             342     SetDiquarkSuppression(ProbDqADq);
343        } else {                                   343        } else {
344           HadronDefinition= DiQuarkSplitup(str    344           HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
345        }                                          345        }      
346                                                   346 
347        if ( HadronDefinition == NULL ) return     347        if ( HadronDefinition == NULL ) return NULL;
348                                                   348 
349        #ifdef debug_QGSMfragmentation             349        #ifdef debug_QGSMfragmentation
350        G4cout<<"The parton "<<string->GetDecay    350        G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
351              <<" produces hadron "<<HadronDefi    351              <<" produces hadron "<<HadronDefinition->GetParticleName()
352              <<" and is transformed to "<<newS    352              <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
353        G4cout<<"The side of the string decay L    353        G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
354        #endif                                     354        #endif
355        // create new String from old, ie. keep    355        // create new String from old, ie. keep Left and Right order, but replace decay
356                                                   356 
357        newString=new G4FragmentingString(*stri    357        newString=new G4FragmentingString(*string,newStringEnd); // To store possible
358                                                   358                                                                 // quark containt of new string
359                                                   359 
360        #ifdef debug_QGSMfragmentation             360        #ifdef debug_QGSMfragmentation
361        G4cout<<"An attempt to determine its en    361        G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
362        #endif                                     362        #endif
363        G4LorentzVector* HadronMomentum=SplitEa    363        G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
364                                                   364 
365        delete newString; newString=0;             365        delete newString; newString=0;
366                                                   366   
367        G4KineticTrack * Hadron =0;                367        G4KineticTrack * Hadron =0;
368        if ( HadronMomentum != 0 ) {               368        if ( HadronMomentum != 0 ) {
369                                                   369 
370            #ifdef debug_QGSMfragmentation         370            #ifdef debug_QGSMfragmentation                     
371            G4cout<<"The attempt was successful    371            G4cout<<"The attempt was successful"<<G4endl;
372            #endif                                 372            #endif
373      G4ThreeVector   Pos;                         373      G4ThreeVector   Pos;
374      Hadron = new G4KineticTrack(HadronDefinit    374      Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
375                                                   375 
376          newString=new G4FragmentingString(*st    376          newString=new G4FragmentingString(*string,newStringEnd,HadronMomentum);
377                                                   377 
378      delete HadronMomentum;                       378      delete HadronMomentum;
379        }                                          379        }
380        else                                       380        else
381        {                                          381        {
382          #ifdef debug_QGSMfragmentation           382          #ifdef debug_QGSMfragmentation
383          G4cout<<"The attempt was not successf    383          G4cout<<"The attempt was not successful !!!"<<G4endl;
384          #endif                                   384          #endif
385        }                                          385        }
386                                                   386 
387        #ifdef debug_VStringDecay                  387        #ifdef debug_VStringDecay
388        G4cout<<"End SplitUP (G4VLongitudinalSt    388        G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
389        #endif                                     389        #endif
390                                                   390 
391        return Hadron;                             391        return Hadron;
392 }                                                 392 }
393                                                   393 
394 //--------------------------------------------    394 //-----------------------------------------------------------------------------
395                                                   395 
396 G4ParticleDefinition *G4QGSMFragmentation::DiQ    396 G4ParticleDefinition *G4QGSMFragmentation::DiQuarkSplitup( G4ParticleDefinition* decay,
397                                                   397                                                            G4ParticleDefinition *&created )
398 {                                                 398 {
399    //... can Diquark break or not?                399    //... can Diquark break or not?
400    if (G4UniformRand() < DiquarkBreakProb )  /    400    if (G4UniformRand() < DiquarkBreakProb )  //... Diquark break
401    {                                              401    {
402       G4int stableQuarkEncoding = decay->GetPD    402       G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
403       G4int decayQuarkEncoding = (decay->GetPD    403       G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
404                                                   404 
405       if (G4UniformRand() < 0.5)                  405       if (G4UniformRand() < 0.5)
406       {                                           406       {
407          G4int Swap = stableQuarkEncoding;        407          G4int Swap = stableQuarkEncoding;
408          stableQuarkEncoding = decayQuarkEncod    408          stableQuarkEncoding = decayQuarkEncoding;
409          decayQuarkEncoding = Swap;               409          decayQuarkEncoding = Swap;
410       }                                           410       }
411                                                   411 
412       G4int IsParticle=(decayQuarkEncoding>0)     412       G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;  // if we have a quark, we need antiquark
413                                                   413 
414       G4double StrSup=GetStrangeSuppress();       414       G4double StrSup=GetStrangeSuppress();
415       SetStrangenessSuppression((1.0 - 0.07)/2    415       SetStrangenessSuppression((1.0 - 0.07)/2.);  // Prob qq->K qq' 0.07
416       pDefPair QuarkPair = CreatePartonPair(Is    416       pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
417       SetStrangenessSuppression(StrSup);          417       SetStrangenessSuppression(StrSup);
418                                                   418 
419       //... Build new Diquark                     419       //... Build new Diquark
420       G4int QuarkEncoding=QuarkPair.second->Ge    420       G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
421       G4int i10  = std::max(std::abs(QuarkEnco    421       G4int i10  = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
422       G4int i20  = std::min(std::abs(QuarkEnco    422       G4int i20  = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
423       G4int spin = (i10 != i20 && G4UniformRan    423       G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
424       G4int NewDecayEncoding = -1*IsParticle*(    424       G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
425                                                   425       
426       created = FindParticle(NewDecayEncoding)    426       created = FindParticle(NewDecayEncoding);
427       G4ParticleDefinition * decayQuark=FindPa    427       G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
428       G4ParticleDefinition * had=hadronizer->B    428       G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
429                                                   429 
430       DecayQuark = decay->GetPDGEncoding();       430       DecayQuark = decay->GetPDGEncoding();
431       NewQuark   = NewDecayEncoding;              431       NewQuark   = NewDecayEncoding;
432                                                   432 
433       return had;                                 433       return had;
434                                                   434 
435    } else {  //... Diquark does not break         435    } else {  //... Diquark does not break
436                                                   436 
437       G4int IsParticle=(decay->GetPDGEncoding(    437       G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;  // if we have a diquark, we need quark)
438       G4double StrSup=GetStrangeSuppress();  /    438       G4double StrSup=GetStrangeSuppress();  // for changing s-sbar production
439       SetStrangenessSuppression((1.0 - 0.07)/2    439       SetStrangenessSuppression((1.0 - 0.07)/2.);
440       pDefPair QuarkPair = CreatePartonPair(Is    440       pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
441       SetStrangenessSuppression(StrSup);          441       SetStrangenessSuppression(StrSup);
442                                                   442 
443       created = QuarkPair.second;                 443       created = QuarkPair.second;
444                                                   444 
445       DecayQuark = decay->GetPDGEncoding();       445       DecayQuark = decay->GetPDGEncoding();
446       NewQuark   = created->GetPDGEncoding();     446       NewQuark   = created->GetPDGEncoding();
447                                                   447 
448       G4ParticleDefinition * had=hadronizer->B    448       G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
449       return had;                                 449       return had;
450    }                                              450    }
451 }                                                 451 }
452                                                   452 
453 //--------------------------------------------    453 //-----------------------------------------------------------------------------------------
454                                                   454 
455 G4LorentzVector * G4QGSMFragmentation::SplitEa    455 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
456                                                   456                                                   G4FragmentingString * string,  
457                                                   457                                                   G4FragmentingString * NewString)
458 {                                                 458 {
459        G4double HadronMass = pHadron->GetPDGMa    459        G4double HadronMass = pHadron->GetPDGMass();
460                                                   460 
461        SetMinimalStringMass(NewString);           461        SetMinimalStringMass(NewString);
462                                                   462 
463        if ( MinimalStringMass < 0.0 ) return n    463        if ( MinimalStringMass < 0.0 ) return nullptr;
464                                                   464 
465        #ifdef debug_QGSMfragmentation             465        #ifdef debug_QGSMfragmentation
466        G4cout<<"G4QGSMFragmentation::SplitEand    466        G4cout<<"G4QGSMFragmentation::SplitEandP "<<pHadron->GetParticleName()<<G4endl;
467        G4cout<<"String 4 mom, String M "<<stri    467        G4cout<<"String 4 mom, String M "<<string->Get4Momentum()<<" "<<string->Mass()<<G4endl;
468        G4cout<<"HadM MinimalStringMassLeft Str    468        G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
469              <<string->Mass()<<" "<<HadronMass    469              <<string->Mass()<<" "<<HadronMass+MinimalStringMass<<G4endl;
470        #endif                                     470        #endif
471                                                   471 
472        if (HadronMass + MinimalStringMass > st    472        if (HadronMass + MinimalStringMass > string->Mass())
473        {                                          473        {
474          #ifdef debug_QGSMfragmentation           474          #ifdef debug_QGSMfragmentation
475          G4cout<<"Mass of the string is not su    475          G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
476          #endif                                   476          #endif
477    return 0;                                      477    return 0;
478        }  // have to start all over!              478        }  // have to start all over!
479                                                   479 
480        // calculate and assign hadron transver    480        // calculate and assign hadron transverse momentum component HadronPx andHadronPy
481        G4double StringMT2 = string->MassT2();     481        G4double StringMT2 = string->MassT2();
482        G4double StringMT  = std::sqrt(StringMT    482        G4double StringMT  = std::sqrt(StringMT2);
483                                                   483 
484        G4LorentzVector String4Momentum = strin    484        G4LorentzVector String4Momentum = string->Get4Momentum();
485        String4Momentum.setPz(0.);                 485        String4Momentum.setPz(0.);
486        G4ThreeVector StringPt = String4Momentu    486        G4ThreeVector StringPt = String4Momentum.vect();
487                                                   487 
488        G4ThreeVector HadronPt    , RemSysPt;      488        G4ThreeVector HadronPt    , RemSysPt;
489        G4double      HadronMassT2, ResidualMas    489        G4double      HadronMassT2, ResidualMassT2;
490                                                   490 
491        // Mt distribution is implemented          491        // Mt distribution is implemented
492        G4double HadronMt, Pt, Pt2, phi;           492        G4double HadronMt, Pt, Pt2, phi;
493                                                   493 
494        //...  sample Pt of the hadron             494        //...  sample Pt of the hadron
495        G4int attempt=0;                           495        G4int attempt=0;
496        do                                         496        do
497        {                                          497        {
498          attempt++; if (attempt > StringLoopIn    498          attempt++; if (attempt > StringLoopInterrupt) return 0;
499                                                   499 
500          HadronMt = HadronMass - 200.0*G4Log(G    500          HadronMt = HadronMass - 200.0*G4Log(G4UniformRand());    // 200.0 must be tuned
501          Pt2 = sqr(HadronMt)-sqr(HadronMass);     501          Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2);
502          phi = 2.*pi*G4UniformRand();             502          phi = 2.*pi*G4UniformRand();
503          G4ThreeVector SampleQuarkPtw= G4Three    503          G4ThreeVector SampleQuarkPtw= G4ThreeVector(Pt*std::cos(phi), Pt*std::sin(phi), 0);
504          HadronPt =SampleQuarkPtw  + string->D    504          HadronPt =SampleQuarkPtw  + string->DecayPt();
505          HadronPt.setZ(0);                        505          HadronPt.setZ(0);
506          RemSysPt = StringPt - HadronPt;          506          RemSysPt = StringPt - HadronPt;
507                                                   507 
508          HadronMassT2 = sqr(HadronMass) + Hadr    508          HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
509          ResidualMassT2=sqr(MinimalStringMass)    509          ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
510                                                   510 
511        } while (std::sqrt(HadronMassT2) + std:    511        } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);  /* Loop checking, 07.08.2015, A.Ribon */
512                                                   512 
513        //...  sample z to define hadron longit    513        //...  sample z to define hadron longitudinal momentum and energy
514        //... but first check the available pha    514        //... but first check the available phase space
515                                                   515 
516        G4double Pz2 = (sqr(StringMT2 - HadronM    516        G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
517           4*HadronMassT2 * ResidualMassT2)/4./    517           4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
518                                                   518 
519        if ( Pz2 < 0 ) {return 0;}          //     519        if ( Pz2 < 0 ) {return 0;}          // have to start all over!
520                                                   520 
521        //... then compute allowed z region  z_    521        //... then compute allowed z region  z_min <= z <= z_max
522                                                   522 
523        G4double Pz = std::sqrt(Pz2);              523        G4double Pz = std::sqrt(Pz2);
524        G4double zMin = (std::sqrt(HadronMassT2    524        G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
525        G4double zMax = (std::sqrt(HadronMassT2    525        G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
526                                                   526 
527        if (zMin >= zMax) return 0;    // have     527        if (zMin >= zMax) return 0;    // have to start all over!
528                                                   528   
529        G4double z = GetLightConeZ(zMin, zMax,     529        G4double z = GetLightConeZ(zMin, zMax,
530                       string->GetDecayParton()    530                       string->GetDecayParton()->GetPDGEncoding(), pHadron,
531                       HadronPt.x(), HadronPt.y    531                       HadronPt.x(), HadronPt.y());      
532                                                   532 
533        //... now compute hadron longitudinal m    533        //... now compute hadron longitudinal momentum and energy
534        // longitudinal hadron momentum compone    534        // longitudinal hadron momentum component HadronPz
535                                                   535 
536        HadronPt.setZ( 0.5* string->GetDecayDir    536        HadronPt.setZ( 0.5* string->GetDecayDirection() *
537           (z * string->LightConeDecay() -         537           (z * string->LightConeDecay() - 
538            HadronMassT2/(z * string->LightCone    538            HadronMassT2/(z * string->LightConeDecay())) );
539        G4double HadronE  = 0.5* (z * string->L    539        G4double HadronE  = 0.5* (z * string->LightConeDecay() + 
540          HadronMassT2/(z * string->LightConeDe    540          HadronMassT2/(z * string->LightConeDecay()) );
541                                                   541 
542        G4LorentzVector * a4Momentum= new G4Lor    542        G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
543                                                   543 
544        #ifdef debug_QGSMfragmentation             544        #ifdef debug_QGSMfragmentation
545        G4cout<<"string->GetDecayDirection() st    545        G4cout<<"string->GetDecayDirection() string->LightConeDecay() "
546              <<string->GetDecayDirection()<<"     546              <<string->GetDecayDirection()<<" "<<string->LightConeDecay()<<G4endl;
547        G4cout<<"HadronPt,HadronE "<<HadronPt<<    547        G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
548        G4cout<<"Out of QGSM SplitEandP "<<G4en    548        G4cout<<"Out of QGSM SplitEandP "<<G4endl;
549        #endif                                     549        #endif
550                                                   550 
551        return a4Momentum;                         551        return a4Momentum;
552 }                                                 552 }
553                                                   553 
554 //--------------------------------------------    554 //----------------------------------------------------------------------------------------------------------
555                                                   555 
556 G4double G4QGSMFragmentation::GetLightConeZ(G4    556 G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int /* PartonEncoding */ ,  
557                                             G4    557                                             G4ParticleDefinition* /* pHadron */, G4double ptx , G4double pty)
558 {                                                 558 {    
559   G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq    559   G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sqr(GeV);
560                                                   560 
561   #ifdef debug_QGSMfragmentation                  561   #ifdef debug_QGSMfragmentation
562   G4cout<<"GetLightConeZ zmin zmax Parton pHad    562   G4cout<<"GetLightConeZ zmin zmax Parton pHadron "<<zmin<<" "<<zmax<<" "/*<< PartonEncoding */
563         <<" "/*<< pHadron->GetParticleName() *    563         <<" "/*<< pHadron->GetParticleName() */ <<G4endl;
564   #endif                                          564   #endif
565                                                   565 
566   G4double z(0.);                                 566   G4double z(0.);    
567   G4int DiQold(0), DiQnew(0);                     567   G4int DiQold(0), DiQnew(0);
568   G4double d1(-1.0), d2(0.);                      568   G4double d1(-1.0), d2(0.);
569   G4double invD1(0.),invD2(0.), r1(0.),r2(0.),    569   G4double invD1(0.),invD2(0.), r1(0.),r2(0.),r12(0.);
570                                                   570 
571   G4int absDecayQuarkCode = std::abs( DecayQua    571   G4int absDecayQuarkCode = std::abs( DecayQuark );
572   G4int absNewQuarkCode   = std::abs( NewQuark    572   G4int absNewQuarkCode   = std::abs( NewQuark   );
573                                                   573 
574   G4int q1(0), q2(0);                             574   G4int q1(0), q2(0);
575   //  q1 = absDecayQuarkCode/1000; q2 = (absDe    575   //  q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100;
576                                                   576 
577   G4int qA(0), qB(0);                             577   G4int qA(0), qB(0);
578   //  qA = absNewQuarkCode/1000;   qB = (absNe    578   //  qA = absNewQuarkCode/1000;   qB = (absNewQuarkCode % 1000)/100;
579                                                   579 
580   if ( (absDecayQuarkCode < 6) && (absNewQuark    580   if ( (absDecayQuarkCode < 6) && (absNewQuarkCode < 6) ) {
581     d1 = FFq2q[absDecayQuarkCode-1][absNewQuar    581     d1 = FFq2q[absDecayQuarkCode-1][absNewQuarkCode-1][0]; d2 = FFq2q[absDecayQuarkCode-1][absNewQuarkCode-1][1];
582   }                                               582   }
583                                                   583 
584   if ( (absDecayQuarkCode < 6) && (absNewQuark    584   if ( (absDecayQuarkCode < 6) && (absNewQuarkCode > 6) ) {
585    qA = absNewQuarkCode/1000;   qB = (absNewQu    585    qA = absNewQuarkCode/1000;   qB = (absNewQuarkCode % 1000)/100;   DiQnew = IndexDiQ[qA-1][qB-1];
586    d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0]    586    d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0]; d2 = FFq2q[absDecayQuarkCode-1][DiQnew][1];
587   }                                               587   }
588                                                   588 
589   if ( (absDecayQuarkCode > 6) && (absNewQuark    589   if ( (absDecayQuarkCode > 6) && (absNewQuarkCode < 6) ) {
590     q1 = absDecayQuarkCode/1000; q2 = (absDeca    590     q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; DiQold = IndexDiQ[q1-1][q2-1];
591     d1 = FFqq2q[DiQold][absNewQuarkCode-1][0];    591     d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; d2 = FFqq2q[DiQold][absNewQuarkCode-1][1];
592   }                                               592   }
593                                                   593 
594   if ( d1 < 0. ) {                                594   if ( d1 < 0. ) {
595     q1 = absDecayQuarkCode/1000; q2 = (absDeca    595     q1 = absDecayQuarkCode/1000; q2 = (absDecayQuarkCode % 1000)/100; DiQold = IndexDiQ[q1-1][q2-1];
596     qA = absNewQuarkCode/1000;   qB = (absNewQ    596     qA = absNewQuarkCode/1000;   qB = (absNewQuarkCode % 1000)/100;   DiQnew = IndexDiQ[qA-1][qB-1];
597     d1 = FFqq2qq[DiQold][DiQnew][0]; d2 = FFqq    597     d1 = FFqq2qq[DiQold][DiQnew][0]; d2 = FFqq2qq[DiQold][DiQnew][1];
598   }                                               598   }
599                                                   599 
600   d2 +=lambda;                                    600   d2 +=lambda;
601   d1+=1.0; d2+=1.0;                               601   d1+=1.0; d2+=1.0;
602                                                   602 
603   invD1=1./d1; invD2=1./d2;                       603   invD1=1./d1; invD2=1./d2;
604                                                   604 
605   const G4int maxNumberOfLoops = 10000;           605   const G4int maxNumberOfLoops = 10000;
606   G4int loopCounter = 0;                          606   G4int loopCounter = 0;
607   do  // Jong's algorithm                         607   do  // Jong's algorithm
608   {                                               608   {
609     r1=G4Pow::GetInstance()->powA(G4UniformRan    609     r1=G4Pow::GetInstance()->powA(G4UniformRand(),invD1);
610     r2=G4Pow::GetInstance()->powA(G4UniformRan    610     r2=G4Pow::GetInstance()->powA(G4UniformRand(),invD2);
611     r12=r1+r2;                                    611     r12=r1+r2;
612     z=r1/r12;                                     612     z=r1/r12;
613   } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z    613   } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z <= zmax))) && 
614             ++loopCounter < maxNumberOfLoops )    614             ++loopCounter < maxNumberOfLoops );  /* Loop checking, 07.08.2015, A.Ribon */
615                                                   615 
616   if ( loopCounter >= maxNumberOfLoops ) {        616   if ( loopCounter >= maxNumberOfLoops ) {
617     z = 0.5*(zmin + zmax);  // Just a value be    617     z = 0.5*(zmin + zmax);  // Just a value between zmin and zmax, no physics considerations at all! 
618   }                                               618   }
619                                                   619 
620   return z;                                       620   return z;
621 }                                                 621 }
622                                                   622 
623 //--------------------------------------------    623 //-----------------------------------------------------------------------------------------
624                                                   624 
625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme    625 G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
626                     G4KineticTrackVector * Lef    626                     G4KineticTrackVector * LeftVector,
627                   G4KineticTrackVector * Right    627                   G4KineticTrackVector * RightVector)
628 {                                                 628 {
629     //... perform last cluster decay              629     //... perform last cluster decay
630                                                   630 
631     G4ThreeVector ClusterVel =string->Get4Mome    631     G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
632     G4double ResidualMass    =string->Mass();     632     G4double ResidualMass    =string->Mass();
633                                                   633 
634     #ifdef debug_QGSMfragmentation                634     #ifdef debug_QGSMfragmentation
635     G4cout<<"Split last-----------------------    635     G4cout<<"Split last-----------------------------------------"<<G4endl;
636     G4cout<<"StrMass "<<ResidualMass<<" q's "     636     G4cout<<"StrMass "<<ResidualMass<<" q's "
637           <<string->GetLeftParton()->GetPartic    637           <<string->GetLeftParton()->GetParticleName()<<" "
638           <<string->GetRightParton()->GetParti    638           <<string->GetRightParton()->GetParticleName()<<G4endl;
639     #endif                                        639     #endif
640                                                   640 
641     G4int cClusterInterrupt = 0;                  641     G4int cClusterInterrupt = 0;
642     G4ParticleDefinition *LeftHadron = nullptr    642     G4ParticleDefinition *LeftHadron = nullptr;
643     G4ParticleDefinition *RightHadron = nullpt    643     G4ParticleDefinition *RightHadron = nullptr;
644     const G4int maxNumberOfLoops = 1000;          644     const G4int maxNumberOfLoops = 1000;
645     G4int loopCounter = 0;                        645     G4int loopCounter = 0;
646                                                   646 
647     G4double LeftHadronMass(0.); G4double Righ    647     G4double LeftHadronMass(0.); G4double RightHadronMass(0.);
648     do                                            648     do
649     {                                             649     {
650         if (cClusterInterrupt++ >= ClusterLoop    650         if (cClusterInterrupt++ >= ClusterLoopInterrupt) return false;
651         LeftHadronMass = -MaxMass; RightHadron    651         LeftHadronMass = -MaxMass; RightHadronMass = -MaxMass;
652                                                   652 
653   G4ParticleDefinition * quark = nullptr;         653   G4ParticleDefinition * quark = nullptr;
654   string->SetLeftPartonStable(); // to query q    654   string->SetLeftPartonStable(); // to query quark contents..
655                                                   655 
656   if (string->DecayIsQuark() && string->Stable    656   if (string->DecayIsQuark() && string->StableIsQuark() ) 
657   {                                               657   {
658     //... there are quarks on cluster ends        658     //... there are quarks on cluster ends
659                                                   659 
660     G4int IsParticle=(string->GetLeftParton()-    660     G4int IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 
661                 // if we have a quark, we need    661                 // if we have a quark, we need antiquark or diquark
662                                                   662 
663     pDefPair QuarkPair = CreatePartonPair(IsPa    663     pDefPair QuarkPair = CreatePartonPair(IsParticle);
664     quark = QuarkPair.second;                     664     quark = QuarkPair.second;
665                                                   665 
666     LeftHadron= hadronizer->Build(QuarkPair.fi    666     LeftHadron= hadronizer->Build(QuarkPair.first, string->GetLeftParton());
667           if ( LeftHadron == NULL ) continue;     667           if ( LeftHadron == NULL ) continue;
668           RightHadron = hadronizer->Build(stri    668           RightHadron = hadronizer->Build(string->GetRightParton(), quark);
669           if ( RightHadron == NULL ) continue;    669           if ( RightHadron == NULL ) continue;    
670   } else if( (!string->DecayIsQuark() &&  stri    670   } else if( (!string->DecayIsQuark() &&  string->StableIsQuark() ) ||   
671              ( string->DecayIsQuark() && !stri    671              ( string->DecayIsQuark() && !string->StableIsQuark() )   ) {
672     //... there is a Diquark on one of cluster    672     //... there is a Diquark on one of cluster ends
673     G4int IsParticle;                             673     G4int IsParticle;
674     if ( string->StableIsQuark() ) {              674     if ( string->StableIsQuark() ) {
675       IsParticle=(string->GetLeftParton()->Get    675       IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 
676     } else {                                      676     } else {
677       IsParticle=(string->GetLeftParton()->Get    677       IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
678     }                                             678     }
679                                                   679 
680           //G4double ProbSaS   = 1.0 - 2.0 * G    680           //G4double ProbSaS   = 1.0 - 2.0 * GetStrangeSuppress();
681           //G4double ActualProb = ProbSaS * 1.    681           //G4double ActualProb = ProbSaS * 1.4;
682           //SetStrangenessSuppression((1.0-Act    682           //SetStrangenessSuppression((1.0-ActualProb)/2.0);
683                                                   683 
684           pDefPair QuarkPair = CreatePartonPai    684           pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
685           //SetStrangenessSuppression((1.0-Pro    685           //SetStrangenessSuppression((1.0-ProbSaS)/2.0);
686           quark = QuarkPair.second;               686           quark = QuarkPair.second;
687           LeftHadron=hadronizer->Build(QuarkPa    687           LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
688           if ( LeftHadron == NULL ) continue;     688           if ( LeftHadron == NULL ) continue;
689           RightHadron = hadronizer->Build(stri    689           RightHadron = hadronizer->Build(string->GetRightParton(), quark);
690           if ( RightHadron == NULL ) continue;    690           if ( RightHadron == NULL ) continue;
691         } else {  // Diquark and anti-diquark     691         } else {  // Diquark and anti-diquark are on the string ends 
692           if (cClusterInterrupt++ >= ClusterLo    692           if (cClusterInterrupt++ >= ClusterLoopInterrupt) return false;
693           G4int LeftQuark1= string->GetLeftPar    693           G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
694           G4int LeftQuark2=(string->GetLeftPar    694           G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
695           G4int RightQuark1= string->GetRightP    695           G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
696           G4int RightQuark2=(string->GetRightP    696           G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
697           if (G4UniformRand()<0.5) {              697           if (G4UniformRand()<0.5) {
698             LeftHadron  =hadronizer->Build(Fin    698             LeftHadron  =hadronizer->Build(FindParticle( LeftQuark1), FindParticle(RightQuark1));
699             RightHadron =hadronizer->Build(Fin    699             RightHadron =hadronizer->Build(FindParticle( LeftQuark2), FindParticle(RightQuark2));
700           } else {                                700           } else {
701             LeftHadron  =hadronizer->Build(Fin    701             LeftHadron  =hadronizer->Build(FindParticle( LeftQuark1), FindParticle(RightQuark2));
702             RightHadron =hadronizer->Build(Fin    702             RightHadron =hadronizer->Build(FindParticle( LeftQuark2), FindParticle(RightQuark1));
703           }                                       703           }
704     if ( (LeftHadron == NULL) || (RightHadron     704     if ( (LeftHadron == NULL) || (RightHadron == NULL) ) continue;
705         }                                         705         }
706         LeftHadronMass  = LeftHadron->GetPDGMa    706         LeftHadronMass  = LeftHadron->GetPDGMass();
707         RightHadronMass = RightHadron->GetPDGM    707         RightHadronMass = RightHadron->GetPDGMass();
708         //... repeat procedure, if mass of clu    708         //... repeat procedure, if mass of cluster is too low to produce hadrons
709     } while ( ( ResidualMass <= LeftHadronMass    709     } while ( ( ResidualMass <= LeftHadronMass + RightHadronMass )
710               && ++loopCounter < maxNumberOfLo    710               && ++loopCounter < maxNumberOfLoops );  /* Loop checking, 07.08.2015, A.Ribon */
711                                                   711 
712     if ( loopCounter >= maxNumberOfLoops ) {      712     if ( loopCounter >= maxNumberOfLoops ) {
713       return false;                               713       return false;
714     }                                             714     }
715                                                   715 
716     //... compute hadron momenta and energies     716     //... compute hadron momenta and energies   
717     G4LorentzVector  LeftMom, RightMom;           717     G4LorentzVector  LeftMom, RightMom;
718     G4ThreeVector    Pos;                         718     G4ThreeVector    Pos;
719     Sample4Momentum(&LeftMom , LeftHadron->Get    719     Sample4Momentum(&LeftMom , LeftHadron->GetPDGMass() , 
720                     &RightMom, RightHadron->Ge    720                     &RightMom, RightHadron->GetPDGMass(), ResidualMass);
721     LeftMom.boost(ClusterVel);                    721     LeftMom.boost(ClusterVel);
722     RightMom.boost(ClusterVel);                   722     RightMom.boost(ClusterVel);
723                                                   723 
724     #ifdef debug_QGSMfragmentation                724     #ifdef debug_QGSMfragmentation
725     G4cout<<LeftHadron->GetParticleName()<<" "    725     G4cout<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
726     G4cout<<"Left  Hadrom P M "<<LeftMom<<" "<    726     G4cout<<"Left  Hadrom P M "<<LeftMom<<" "<<LeftMom.mag()<<G4endl;
727     G4cout<<"Right Hadrom P M "<<RightMom<<" "    727     G4cout<<"Right Hadrom P M "<<RightMom<<" "<<RightMom.mag()<<G4endl;
728     #endif                                        728     #endif
729                                                   729 
730     LeftVector->push_back(new G4KineticTrack(L    730     LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
731     RightVector->push_back(new G4KineticTrack(    731     RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
732                                                   732 
733     return true;                                  733     return true;
734 }                                                 734 }
735                                                   735 
736 //--------------------------------------------    736 //----------------------------------------------------------------------------------------------------------
737                                                   737 
738 void G4QGSMFragmentation::Sample4Momentum(G4Lo    738 void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom    , G4double Mass    , 
739                                           G4Lo    739                                           G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) 
740 {                                                 740 {
741     #ifdef debug_QGSMfragmentation                741     #ifdef debug_QGSMfragmentation
742     G4cout<<"Sample4Momentum Last-------------    742     G4cout<<"Sample4Momentum Last-----------------------------------------"<<G4endl;
743     G4cout<<"  StrMass "<<InitialMass<<" Mass1    743     G4cout<<"  StrMass "<<InitialMass<<" Mass1 "<<Mass<<" Mass2 "<<AntiMass<<G4endl;
744     G4cout<<"  SumMass "<<Mass+AntiMass<<G4end    744     G4cout<<"  SumMass "<<Mass+AntiMass<<G4endl;
745     #endif                                        745     #endif
746                                                   746 
747     G4double r_val = sqr(InitialMass*InitialMa    747     G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
748     G4double Pabs = (r_val > 0.)? std::sqrt(r_    748     G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
749                                                   749 
750     //... sample unit vector                      750     //... sample unit vector       
751     G4double pz = 1. - 2.*G4UniformRand();        751     G4double pz = 1. - 2.*G4UniformRand();  
752     G4double st     = std::sqrt(1. - pz * pz)*    752     G4double st     = std::sqrt(1. - pz * pz)*Pabs;
753     G4double phi    = 2.*pi*G4UniformRand();      753     G4double phi    = 2.*pi*G4UniformRand();
754     G4double px = st*std::cos(phi);               754     G4double px = st*std::cos(phi);
755     G4double py = st*std::sin(phi);               755     G4double py = st*std::sin(phi);
756     pz *= Pabs;                                   756     pz *= Pabs;
757                                                   757     
758     Mom->setPx(px); Mom->setPy(py); Mom->setPz    758     Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
759     Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass)    759     Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
760                                                   760 
761     AntiMom->setPx(-px); AntiMom->setPy(-py);     761     AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
762     AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM    762     AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
763 }                                                 763 }
764                                                   764 
765 //--------------------------------------------    765 //----------------------------------------------------------------------------------------------------------
766                                                   766 
767 void G4QGSMFragmentation::SetFFq2q()  // q-> q    767 void G4QGSMFragmentation::SetFFq2q()  // q-> q' + Meson (q anti q')
768 {                                                 768 {
769   for (G4int i=0; i < 5; i++) {                   769   for (G4int i=0; i < 5; i++) {
770     FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -a    770     FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -arho + alft;  // q->d + (q dbar) Pi0, Eta, Eta', Rho0, omega
771     FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -a    771     FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -arho + alft;  // q->u + (q ubar) Pi-, Rho-
772     FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -a    772     FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -aphi + alft;  // q->s + (q sbar) K0, K*0
773     FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -a    773     FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -aJPs + alft;  // q->c + (q+cbar) D-, D*-
774     FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -a    774     FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -aUps + alft;  // q->b + (q bbar) EtaB, Upsilon
775   }                                               775   }
776 }                                                 776 }
777                                                   777 
778 //--------------------------------------------    778 //----------------------------------------------------------------------------------------------------------
779                                                   779 
780 void G4QGSMFragmentation::SetFFq2qq()  // q->     780 void G4QGSMFragmentation::SetFFq2qq()  // q-> anti (q1'q2') + Baryon (q + q1 + q2)
781 {                                                 781 {
782   for (G4int i=0; i < 5; i++) {                   782   for (G4int i=0; i < 5; i++) {
783     FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1]     783     FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] = arho - 2.0*an    + alft;  // q->dd bar + (q dd)
784     FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1]     784     FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] = arho - 2.0*an    + alft;  // q->ud bar + (q ud)
785     FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1]     785     FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] = arho - 2.0*ala   + alft;  // q->sd bar + (q sd) 
786     FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1]     786     FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] = arho - 2.0*alaC  + alft;  // q->cd bar + (q cd)
787     FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1]     787     FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] = arho - 2.0*alaB  + alft;  // q->bd bar + (q bd)
788     FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1]     788     FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] = arho - 2.0*an    + alft;  // q->uu bar + (q uu)
789     FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1]     789     FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] = arho - 2.0*ala   + alft;  // q->su bar + (q su)
790     FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1]     790     FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] = arho - 2.0*alaC  + alft;  // q->cu bar + (q cu)
791     FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1]     791     FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] = arho - 2.0*alaB  + alft;  // q->bu bar + (q bu)
792     FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1]     792     FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] = arho - 2.0*aXi   + alft;  // q->ss bar + (q ss)    
793     FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1]     793     FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] = arho - 2.0*aXiC  + alft;  // q->cs bar + (q cs)
794     FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1]     794     FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] = arho - 2.0*aXiB  + alft;  // q->bs bar + (q bc)
795     FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1]     795     FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] = arho - 2.0*aXiCC + alft;  // q->cc bar + (q cc) 
796     FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1]     796     FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] = arho - 2.0*aXiCB + alft;  // q->bc bar + (q bc)
797     FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1]     797     FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] = arho - 2.0*aXiBB + alft;  // q->bb bar + (q bb)
798   }                                               798   }
799 }                                                 799 }
800                                                   800 
801 //--------------------------------------------    801 //----------------------------------------------------------------------------------------------------------
802                                                   802 
803 void G4QGSMFragmentation::SetFFqq2q()  // q1q2    803 void G4QGSMFragmentation::SetFFqq2q()  // q1q2-> anti(q') + Baryon (q1 + q2 + q')
804 {                                                 804 {
805   for (G4int i=0; i < 15; i++) {                  805   for (G4int i=0; i < 15; i++) {   
806     FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[    806     FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[i][0][1] = -arho + alft;
807     FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[    807     FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[i][1][1] = -arho + alft;
808     FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[    808     FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[i][2][1] = -aphi + alft;
809     FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[    809     FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[i][3][1] = -aJPs + alft;
810     FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[    810     FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[i][4][1] = -aUps + alft;
811   }                                               811   }
812 }                                                 812 }
813                                                   813 
814 //--------------------------------------------    814 //----------------------------------------------------------------------------------------------------------
815                                                   815 
816 void G4QGSMFragmentation::SetFFqq2qq()  // q1(    816 void G4QGSMFragmentation::SetFFqq2qq()  // q1(q2)-> q'(q2) + Meson(q1 anti q')
817 {                                                 817 {
818   for (G4int i=0; i < 15; i++) {                  818   for (G4int i=0; i < 15; i++) {
819     FFqq2qq[i][0][0] = 0.  ;  FFqq2qq[i][0][1]    819     FFqq2qq[i][0][0] = 0.  ;  FFqq2qq[i][0][1] = 2.0*arho - 2.0*an -arho + alft;  // dd -> dd + Pi0 (d d bar)
820     FFqq2qq[i][1][0] = 0.  ;  FFqq2qq[i][1][1]    820     FFqq2qq[i][1][0] = 0.  ;  FFqq2qq[i][1][1] = 2.0*arho - 2.0*an -arho + alft;  // dd -> ud + Pi- (d u bar)
821     FFqq2qq[i][2][0] = 0.  ;  FFqq2qq[i][2][1]    821     FFqq2qq[i][2][0] = 0.  ;  FFqq2qq[i][2][1] = 2.0*arho - 2.0*an -aphi + alft;  // dd -> sd + K0  (d s bar)
822     FFqq2qq[i][3][0] = 0.  ;  FFqq2qq[i][3][1]    822     FFqq2qq[i][3][0] = 0.  ;  FFqq2qq[i][3][1] = 2.0*arho - 2.0*an -aJPs + alft;  // dd -> cd + D-  (d c bar)
823     FFqq2qq[i][4][0] = 0.  ;  FFqq2qq[i][4][1]    823     FFqq2qq[i][4][0] = 0.  ;  FFqq2qq[i][4][1] = 2.0*arho - 2.0*an -aUps + alft;  // dd -> bd + B0  (d b bar)
824   }                                               824   }
825 }                                                 825 }
826                                                   826 
827                                                   827