Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 // -----------------------------------------------------------------------------
 29 //      GEANT 4 class implementation file
 30 //
 31 //      History: first implementation, Maxim Komogorov, 10-Jul-1998
 32 // -----------------------------------------------------------------------------
 33 #include "G4LundStringFragmentation.hh"
 34 #include "G4PhysicalConstants.hh"
 35 #include "G4SystemOfUnits.hh"
 36 #include "Randomize.hh"
 37 #include "G4FragmentingString.hh"
 38 #include "G4DiQuarks.hh"
 39 #include "G4Quarks.hh"
 40 #include "G4HadronicParameters.hh"
 41 #include "G4Exp.hh"
 42 #include "G4Pow.hh"
 43 
 44 //#define debug_LUNDfragmentation
 45 
 46 // Class G4LundStringFragmentation 
 47 //*************************************************************************************
 48 
 49 G4LundStringFragmentation::G4LundStringFragmentation()
 50   : G4VLongitudinalStringDecay("LundStringFragmentation")
 51 {
 52     SetMassCut(210.*MeV);   //  Mpi + Delta
 53                             // For ProduceOneHadron it is required
 54                             // that no one pi-meson can be produced.
 55     SigmaQT = 0.435 * GeV;
 56     Tmt = 190.0 * MeV;       
 57 
 58     SetStringTensionParameter(1.*GeV/fermi);
 59     SetDiquarkBreakProbability(0.3);
 60 
 61     SetStrangenessSuppression((1.0 - 0.12)/2.0);
 62     SetDiquarkSuppression(0.07);
 63 
 64     // Check if charmed and bottom hadrons are enabled: if this is the case, then
 65     // set the non-zero probabilities for c-cbar and b-bbar creation from the vacuum,
 66     // else set them to 0.0. If these probabilities are/aren't zero then charmed or bottom
 67     // hadrons can't/can be created during the string fragmentation of ordinary
 68     // (i.e. not heavy) projectile hadron nuclear reactions.
 69     if ( G4HadronicParameters::Instance()->EnableBCParticles() ) {
 70       SetProbCCbar(0.0002);  // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094; tuned by Uzhi Oct. 2022
 71       SetProbBBbar(5.0e-5);  // According to O.I. Piskunova Yad. Fiz. 56 (1993) 1094
 72     } else {
 73       SetProbCCbar(0.0);
 74       SetProbBBbar(0.0);
 75     }
 76     
 77     SetMinMasses();  // For treating of small string decays
 78 }
 79 
 80 //--------------------------------------------------------------------------------------
 81 
 82 G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
 83 {
 84   // Can no longer modify Parameters for Fragmentation.
 85 
 86   PastInitPhase=true;
 87 
 88   G4FragmentingString  aString(theString);
 89   SetMinimalStringMass(&aString);
 90 
 91         #ifdef debug_LUNDfragmentation
 92         G4cout<<G4endl<<"LUND StringFragmentation ------------------------------------"<<G4endl;
 93         G4cout<<G4endl<<"LUND StringFragm: String Mass "
 94                       <<theString.Get4Momentum().mag()<<G4endl
 95                       <<"4Mom "<<theString.Get4Momentum()<<G4endl
 96                       <<"------------------------------------"<<G4endl;
 97         G4cout<<"String ends Direct "<<theString.GetLeftParton()->GetPDGcode()<<" "
 98                                      <<theString.GetRightParton()->GetPDGcode()<<" "
 99                                      <<theString.GetDirection()<< G4endl;
100         G4cout<<"Left  mom "<<theString.GetLeftParton()->Get4Momentum()<<G4endl;
101         G4cout<<"Right mom "<<theString.GetRightParton()->Get4Momentum()<<G4endl<<G4endl;
102         G4cout<<"Check for Fragmentation "<<G4endl;
103         #endif
104 
105   G4KineticTrackVector * LeftVector(0);
106 
107   if (!aString.IsAFourQuarkString() && !IsItFragmentable(&aString))
108   {
109                 #ifdef debug_LUNDfragmentation
110                 G4cout<<"Non fragmentable - the string is converted to one hadron "<<G4endl;
111                 #endif
112                 // SetMassCut(210.*MeV);  // For ProduceOneHadron it is required
113                                           // that no one pi-meson can be produced.
114 
115     G4double Mcut = GetMassCut();
116     SetMassCut(10000.*MeV);
117     LeftVector=ProduceOneHadron(&theString);
118     SetMassCut(Mcut);
119 
120     if ( LeftVector )
121     {
122       if ( LeftVector->size() > 0)
123                   {
124             LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
125             LeftVector->operator[](0)->SetPosition(theString.GetPosition());
126                   }
127       if (LeftVector->size() > 1)
128                   {
129             // 2 hadrons created from qq-qqbar are stored
130       LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
131       LeftVector->operator[](1)->SetPosition(theString.GetPosition());
132       }
133     }   
134     return LeftVector;
135   }
136 
137         #ifdef debug_LUNDfragmentation
138         G4cout<<"The string will be fragmented. "<<G4endl;
139         #endif
140 
141   // The string can fragment. At least two particles can be produced.
142              LeftVector =new G4KineticTrackVector;
143   G4KineticTrackVector * RightVector=new G4KineticTrackVector;
144 
145         G4bool success = Loop_toFragmentString(theString, LeftVector, RightVector);
146 
147   if ( ! success )
148   {
149     std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
150     LeftVector->clear();
151     std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
152     delete RightVector;
153     return LeftVector;
154   }
155 
156   // Join Left- and RightVector into LeftVector in correct order.
157   while (!RightVector->empty())
158   {
159     LeftVector->push_back(RightVector->back());
160     RightVector->erase(RightVector->end()-1);
161   }
162   delete RightVector;
163 
164   return LeftVector;
165 }
166 
167 //----------------------------------------------------------------------------------
168 
169 G4bool G4LundStringFragmentation::IsItFragmentable(const G4FragmentingString * const string)
170 {
171   SetMinimalStringMass(string);
172         //G4cout<<"MinM StrM "<<MinimalStringMass<<" "<< string->Get4Momentum().mag()<<G4endl;
173 
174   return std::abs(MinimalStringMass) < string->Get4Momentum().mag();
175 
176         //MinimalStringMass is negative and large for a string with unknown particles in a final 2-particle decay.
177 }
178 
179 //----------------------------------------------------------------------------------------
180 
181 G4bool G4LundStringFragmentation::Loop_toFragmentString( const G4ExcitedString  &theString,
182                                                          G4KineticTrackVector * & LeftVector,
183                                                          G4KineticTrackVector * & RightVector )
184 {
185         #ifdef debug_LUNDfragmentation
186         G4cout<<"Loop_toFrag "<<theString.GetLeftParton()->GetPDGcode()<<" "
187                               <<theString.GetLeftParton()->Get4Momentum()<<G4endl
188               <<"            "<<theString.GetRightParton()->GetPDGcode()<<" "
189                               <<theString.GetRightParton()->Get4Momentum()<<G4endl
190               <<"Direction   "<<theString.GetDirection()<< G4endl;
191         #endif
192 
193         G4LorentzRotation toCmsI, toObserverFrameI;
194 
195   G4bool final_success=false;
196   G4bool inner_success=true;
197 
198   G4int attempt=0;
199 
200   while ( ! final_success && attempt++ < StringLoopInterrupt )
201   {       // If the string fragmentation does not be happend, 
202           // repeat the fragmentation.
203 
204                 G4FragmentingString *currentString = new G4FragmentingString(theString);
205                 toCmsI = currentString->TransformToAlignedCms();
206                 toObserverFrameI = toCmsI.inverse();
207 
208     G4LorentzRotation toCms, toObserverFrame;
209 
210     //G4cout<<"Main loop start whilecounter "<<attempt<<G4endl;
211 
212     // Cleaning up the previously produced hadrons
213     std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
214     LeftVector->clear();
215     std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
216     RightVector->clear();
217 
218     // Main fragmentation loop until the string will not be able to fragment
219     inner_success=true;  // set false on failure.
220                 const G4int maxNumberOfLoops = 1000;
221                 G4int loopCounter = -1;
222     
223     while ( (! StopFragmenting(currentString)) && ++loopCounter < maxNumberOfLoops )
224     {       // Split current string into hadron + new string
225                         #ifdef debug_LUNDfragmentation
226                         G4cout<<"The string will fragment. "<<G4endl;;
227                         //G4cout<<"1 "<<currentString->GetDecayDirection()<<G4endl;
228                         #endif
229       G4FragmentingString *newString=0;  // used as output from SplitUp.
230 
231       toCms=currentString->TransformToAlignedCms();
232                         toObserverFrame= toCms.inverse();
233       
234                         #ifdef debug_LUNDfragmentation
235                         //G4cout<<"CMS Left  mom "<<currentString->GetPleft()<<G4endl;
236                         //G4cout<<"CMS Right mom "<<currentString->GetPright()<<G4endl;
237                         //G4cout<<"CMS String M  "<<currentString->GetPstring()<<G4endl;
238                         #endif
239 
240       G4KineticTrack * Hadron=Splitup(currentString,newString);
241 
242       if ( Hadron != 0 )  // Store the hadron                               
243       {
244                                 #ifdef debug_LUNDfragmentation
245                                 G4cout<<"Hadron prod at fragm. "<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
246                                 //G4cout<<"2 "<<currentString->GetDecayDirection()<<G4endl;
247                                 #endif
248 
249         Hadron->Set4Momentum(toObserverFrame*Hadron->Get4Momentum());
250 
251         G4double TimeOftheStringCreation=theString.GetTimeOfCreation();
252         G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
253 
254         G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
255         G4LorentzVector Momentum = toObserverFrame*Coordinate;
256         Hadron->SetFormationTime(TimeOftheStringCreation + Momentum.e() - fermi/c_light);
257         G4ThreeVector aPosition(Momentum.vect());
258         Hadron->SetPosition(PositionOftheStringCreation+aPosition);
259  
260                                 // Open to protect hadron production at fragmentation 
261         if ( currentString->GetDecayDirection() > 0 )
262                                 {
263           LeftVector->push_back(Hadron); 
264                                 } else
265                                 {
266           RightVector->push_back(Hadron);
267                                 }
268         delete currentString;
269         currentString=newString;
270       } else {
271                           if ( newString ) delete newString;
272                         }
273 
274                         currentString->LorentzRotate(toObserverFrame);
275     };
276 
277                 if ( loopCounter >= maxNumberOfLoops ) {
278                   inner_success=false;
279                 }
280 
281     // Split remaining string into 2 final hadrons.
282                 #ifdef debug_LUNDfragmentation
283                 if (inner_success) G4cout<<"Split remaining string into 2 final hadrons."<<G4endl;
284                 #endif
285 
286     if ( inner_success && SplitLast(currentString, LeftVector, RightVector) )  // Close to protect Last Str. Decay
287     {
288       final_success = true;
289     }
290 
291     delete currentString;
292   }  // End of the loop where we try to fragment the string.
293 
294         G4int sign = +1;
295         if ( theString.GetDirection() < 0 ) sign = -1;
296         for ( unsigned int hadronI = 0; hadronI < LeftVector->size(); ++hadronI ) {
297            G4LorentzVector Tmp = LeftVector->operator[](hadronI)->Get4Momentum();
298            Tmp.setZ(sign*Tmp.getZ());
299            Tmp *= toObserverFrameI;
300            LeftVector->operator[](hadronI)->Set4Momentum(Tmp);
301         }
302         for ( unsigned int hadronI = 0; hadronI < RightVector->size(); ++hadronI ) {
303            G4LorentzVector Tmp = RightVector->operator[](hadronI)->Get4Momentum();
304            Tmp.setZ(sign*Tmp.getZ());
305            Tmp *= toObserverFrameI;
306            RightVector->operator[](hadronI)->Set4Momentum(Tmp);
307         }
308 
309   return final_success;
310 }
311 
312 //----------------------------------------------------------------------------------------
313 
314 G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
315 {
316   SetMinimalStringMass(string); 
317 
318   if ( MinimalStringMass < 0.) return true;
319 
320   if (string->IsAFourQuarkString())
321   {
322     return G4UniformRand() < G4Exp(-0.0005*(string->Mass() - MinimalStringMass));
323   } else {
324            
325                 if (MinimalStringMass < 0.0 ) return false;  // For a string with di-quark having c or b quarks and s, c, b quarks
326 
327     G4bool Result = G4UniformRand() < 
328         G4Exp(-0.66e-6*(string->Mass()*string->Mass() - MinimalStringMass*MinimalStringMass));
329                 // G4bool Result = string->Mass() < MinimalStringMass + 150.*MeV*G4UniformRand();     // a'la LUND
330 
331                 #ifdef debug_LUNDfragmentation 
332                 G4cout<<"StopFragmenting MinimalStringMass string->Mass() "<<MinimalStringMass
333                       <<" "<<string->Mass()<<G4endl;
334                 G4cout<<"StopFragmenting - Yes/No "<<Result<<G4endl;
335                 #endif  
336     return Result;
337   }
338 }
339 
340 //-----------------------------------------------------------------------------
341 
342 G4KineticTrack * G4LundStringFragmentation::Splitup(G4FragmentingString *string, 
343                                   G4FragmentingString *&newString)
344 {
345        #ifdef debug_LUNDfragmentation
346        G4cout<<G4endl;
347        G4cout<<"Start SplitUP ========================="<<G4endl;
348        G4cout<<"String partons: " <<string->GetLeftParton()->GetPDGEncoding()<<" "
349                                   <<string->GetRightParton()->GetPDGEncoding()<<" "
350              <<"Direction "       <<string->GetDecayDirection()<<G4endl;
351        #endif
352 
353        //... random choice of string end to use for creating the hadron (decay)   
354        G4int SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
355        if (SideOfDecay < 0)
356        {
357     string->SetLeftPartonStable();
358        } else
359        {
360           string->SetRightPartonStable();
361        }
362 
363        G4ParticleDefinition *newStringEnd;
364        G4ParticleDefinition * HadronDefinition;
365 
366        G4double StringMass=string->Mass();
367 
368        G4double ProbDqADq = GetDiquarkSuppress();
369        G4double ProbSaS   = 1.0 - 2.0 * GetStrangeSuppress();
370 
371        #ifdef debug_LUNDfragmentation
372        G4cout<<"StrMass DiquarkSuppression           "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl; 
373        #endif   
374 
375        G4int NumberOfpossibleBaryons = 2;
376 
377        if (string->GetLeftParton()->GetParticleSubType()  != "quark") NumberOfpossibleBaryons++; 
378        if (string->GetRightParton()->GetParticleSubType() != "quark") NumberOfpossibleBaryons++;
379 
380        G4double ActualProb  = ProbDqADq ;
381        ActualProb *= (1.0-G4Pow::GetInstance()->powA(NumberOfpossibleBaryons*1400.0/StringMass, 8.0));
382        if(ActualProb <0.0) ActualProb = 0.;
383 
384        SetDiquarkSuppression(ActualProb); 
385 
386        G4double Mth = 1250.0;                                  // 2 Mk + Mpi
387        if ( NumberOfpossibleBaryons == 3 ){Mth = 2520.0;}       // Mlambda/Msigma + Mk + Mpi
388        else if ( NumberOfpossibleBaryons == 4 ){Mth = 2380.0;}  // 2 Mlambda/Msigma + Mk + Mpi
389        else {}
390 
391        ActualProb = ProbSaS;
392        ActualProb *= (1.0 - G4Pow::GetInstance()->powA( Mth/StringMass, 2.5 ));
393        if ( ActualProb < 0.0 ) ActualProb = 0.0;
394        SetStrangenessSuppression((1.0-ActualProb)/2.0);
395 
396        #ifdef debug_LUNDfragmentation
397        G4cout<<"StrMass DiquarkSuppression corrected "<<StringMass<<" "<<GetDiquarkSuppress()<<G4endl; 
398        #endif
399 
400        if (string->DecayIsQuark())
401        {
402           HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
403        } else {
404           HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
405        }
406 
407        SetDiquarkSuppression(ProbDqADq);
408        SetStrangenessSuppression((1.0-ProbSaS)/2.0);
409 
410        if ( HadronDefinition == NULL ) { G4KineticTrack * Hadron =0; return Hadron; }
411 
412        #ifdef debug_LUNDfragmentation
413        G4cout<<"The parton "<<string->GetDecayParton()->GetPDGEncoding()<<" "
414              <<" produces hadron "<<HadronDefinition->GetParticleName()
415              <<" and is transformed to "<<newStringEnd->GetPDGEncoding()<<G4endl;
416        G4cout<<"The side of the string decay Left/Right (1/-1) "<<SideOfDecay<<G4endl;
417        #endif
418        // create new String from old, ie. keep Left and Right order, but replace decay
419 
420        if ( newString ) delete newString;
421 
422        newString=new G4FragmentingString(*string,newStringEnd);  // To store possible quark containt of new string
423 
424        #ifdef debug_LUNDfragmentation
425        G4cout<<"An attempt to determine its energy (SplitEandP)"<<G4endl;
426        #endif
427        G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
428 
429        delete newString; newString=0;
430   
431        G4KineticTrack * Hadron =0;
432        if ( HadronMomentum != 0 ) {
433 
434            #ifdef debug_LUNDfragmentation
435            G4cout<<"The attempt was successful"<<G4endl;
436            #endif
437      G4ThreeVector   Pos;
438      Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
439      
440            if ( newString ) delete newString;
441 
442      newString=new G4FragmentingString(*string,newStringEnd,
443             HadronMomentum);
444      delete HadronMomentum;
445        }
446        else
447        {
448            #ifdef debug_LUNDfragmentation
449            G4cout<<"The attempt was not successful !!!"<<G4endl;
450            #endif
451        }
452 
453        #ifdef debug_LUNDfragmentation
454        G4cout<<"End SplitUP (G4VLongitudinalStringDecay) ====================="<<G4endl;
455        #endif
456 
457        return Hadron;
458 }
459 
460 //-----------------------------------------------------------------------------
461 
462 G4ParticleDefinition * G4LundStringFragmentation::DiQuarkSplitup(G4ParticleDefinition* decay,
463                                                                  G4ParticleDefinition *&created)
464 {
465    G4double StrSup=GetStrangeSuppress();
466    G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.25;
467 
468    //... can Diquark break or not?
469    if (G4UniformRand() < DiquarkBreakProb ){
470 
471       //... Diquark break
472       G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
473       G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
474       if (G4UniformRand() < 0.5)
475       {
476          G4int Swap = stableQuarkEncoding;
477          stableQuarkEncoding = decayQuarkEncoding;
478          decayQuarkEncoding = Swap;
479       }
480 
481       G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;  // if we have a quark, we need antiquark
482 
483       SetStrangenessSuppression((1.0-ProbQQbar)/2.0);
484       pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
485       SetStrangenessSuppression((1.0-StrSup)/2.0);
486 
487       //... Build new Diquark
488       G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
489       G4int i10  = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
490       G4int i20  = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
491       G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
492       G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
493       created = FindParticle(NewDecayEncoding);
494       G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
495       G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
496       StrangeSuppress=StrSup;
497 
498       return had;
499 
500    } else {
501       //... Diquark does not break
502 
503       G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;  // if we have a diquark, we need quark
504 
505       StrangeSuppress=(1.0 - ProbQQbar)/2.0;
506       pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
507 
508       created = QuarkPair.second;
509 
510       G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
511       StrangeSuppress=StrSup;
512 
513       return had;
514    }
515 }
516 
517 //-----------------------------------------------------------------------------
518 
519 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
520                                                         G4FragmentingString *   string, 
521                                                         G4FragmentingString * newString)
522 { 
523   G4LorentzVector String4Momentum=string->Get4Momentum();
524   G4double StringMT2=string->MassT2();
525   G4double StringMT =std::sqrt(StringMT2);
526 
527   G4double HadronMass = pHadron->GetPDGMass();
528   SetMinimalStringMass(newString);
529 
530        if ( MinimalStringMass < 0.0 ) return nullptr;
531 
532         #ifdef debug_LUNDfragmentation
533         G4cout<<G4endl<<"Start LUND SplitEandP "<<G4endl;
534         G4cout<<"String 4 mom, String M and Mt "<<String4Momentum<<" "<<String4Momentum.mag()
535               <<" "<<std::sqrt(StringMT2)<<G4endl;
536         G4cout<<"Hadron "<<pHadron->GetParticleName()<<G4endl;
537         G4cout<<"HadM MinimalStringMassLeft StringM hM+sM "<<HadronMass<<" "<<MinimalStringMass<<" "
538               <<String4Momentum.mag()<<" "<<HadronMass+MinimalStringMass<<G4endl;
539         #endif
540 
541   if ((HadronMass + MinimalStringMass > string->Mass()) || MinimalStringMass < 0.) 
542   {
543           #ifdef debug_LUNDfragmentation
544           G4cout<<"Mass of the string is not sufficient to produce the hadron!"<<G4endl;
545           #endif
546     return 0;
547   }  // have to start all over!
548 
549   String4Momentum.setPz(0.);
550   G4ThreeVector StringPt=String4Momentum.vect();
551         StringPt.setZ(0.);
552   
553   // calculate and assign hadron transverse momentum component HadronPx and HadronPy
554   G4ThreeVector HadronPt    , RemSysPt; 
555   G4double      HadronMassT2, ResidualMassT2;
556   G4double HadronMt, Pt, Pt2, phi;
557 
558         G4double TmtCur = Tmt;
559 
560         if ( (string->GetDecayParton()->GetParticleSubType()== "quark") &&
561              (pHadron->GetBaryonNumber() != 0) ) {
562           TmtCur = Tmt*0.37;              // q->B     
563         } else if ( (string->GetDecayParton()->GetParticleSubType()== "quark") &&
564                     (pHadron->GetBaryonNumber() == 0) ) {
565           //TmtCur = Tmt;                               // q->M
566   } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") && 
567                     (pHadron->GetBaryonNumber() == 0) ) {
568           //TmtCur = Tmt*0.89;                          // qq -> M
569         } else if ( (string->GetDecayParton()->GetParticleSubType()== "di_quark") &&
570                     (pHadron->GetBaryonNumber() != 0) ) {
571           TmtCur = Tmt*1.35;                            // qq -> B
572         }
573 
574         //...  sample Pt of the hadron
575         G4int attempt=0;
576         do
577         {
578           attempt++; if (attempt > StringLoopInterrupt) {return 0;}
579 
580           HadronMt = HadronMass - TmtCur*G4Log(G4UniformRand());
581     Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=std::sqrt(Pt2);
582     phi = 2.*pi*G4UniformRand();
583           HadronPt = G4ThreeVector( Pt*std::cos(phi), Pt*std::sin(phi), 0. );
584           RemSysPt = StringPt - HadronPt;
585           HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
586           ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
587 
588         } while (std::sqrt(HadronMassT2) + std::sqrt(ResidualMassT2) > StringMT);
589 
590   //...  sample z to define hadron longitudinal momentum and energy
591   //... but first check the available phase space
592 
593   G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
594       4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
595 
596   if (Pz2 < 0 ) {return 0;}          // have to start all over!
597 
598   //... then compute allowed z region  z_min <= z <= z_max
599 
600   G4double Pz = std::sqrt(Pz2);
601   G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
602         // G4double zMin = (std::sqrt(HadronMassT2+Pz2) - 0.)/std::sqrt(StringMT2); // For testing purposes
603   G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
604 
605   if (zMin >= zMax) return 0;   // have to start all over!
606 
607   G4double z = GetLightConeZ(zMin, zMax, 
608                  string->GetDecayParton()->GetPDGEncoding(), pHadron,
609                  HadronPt.x(), HadronPt.y());
610 
611   //... now compute hadron longitudinal momentum and energy
612   // longitudinal hadron momentum component HadronPz
613 
614   HadronPt.setZ(0.5* string->GetDecayDirection() *
615            (z * string->LightConeDecay() - 
616           HadronMassT2/(z * string->LightConeDecay())));
617   G4double HadronE  = 0.5* (z * string->LightConeDecay() +
618           HadronMassT2/(z * string->LightConeDecay()));
619 
620   G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
621 
622         #ifdef debug_LUNDfragmentation
623         G4cout<<G4endl<<" string->GetDecayDirection() "<<string->GetDecayDirection()<<G4endl<<G4endl;
624         G4cout<<"string->LightConeDecay() "<<string->LightConeDecay()<<G4endl;
625         G4cout<<"HadronPt,HadronE "<<HadronPt<<" "<<HadronE<<G4endl;
626         G4cout<<"String4Momentum "<<String4Momentum<<G4endl;
627         G4cout<<"Out of LUND SplitEandP "<<G4endl<<G4endl;
628         #endif
629 
630   return a4Momentum;
631 }
632 
633 //-----------------------------------------------------------------------------------------
634 
635 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, 
636                                       G4int PDGEncodingOfDecayParton,
637                                       G4ParticleDefinition* pHadron,
638                                       G4double Px, G4double Py)
639 {
640   G4double Mass = pHadron->GetPDGMass();
641         G4int HadronEncoding=std::abs(pHadron->GetPDGEncoding());
642 
643   G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
644 
645   G4double  Alund, Blund;
646   G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf(1.);
647 
648   if (!((std::abs(PDGEncodingOfDecayParton) > 1000) && (HadronEncoding > 1000)) )
649   {    // ---------------- Quark fragmentation  and qq-> meson ----------------------
650            Alund=1.;
651            Blund=0.7/GeV/GeV;
652 
653      G4double BMt2 = Blund*Mt2;
654      if (Alund == 1.0) {
655        zOfMaxyf=BMt2/(Blund*Mt2 + 1.);}
656      else {  
657          zOfMaxyf = ((1.0+BMt2) - std::sqrt(sqr(1.0-BMt2) + 4.0*BMt2*Alund))/2.0/(1.-Alund);
658      }
659 
660      if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}
661      if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}
662      maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blund*Mt2/zOfMaxyf);
663 
664            const G4int maxNumberOfLoops = 1000;
665            G4int loopCounter = 0;
666      do
667      {
668     z = zmin + G4UniformRand()*(zmax-zmin);
669                 //yf = (1-z)/z * G4Exp(-Blund*Mt2/z);
670     yf = G4Pow::GetInstance()->powA(1.0-z,Alund)/z*G4Exp(-BMt2/z);
671      }
672      while ( (G4UniformRand()*maxYf > yf) && ++loopCounter < maxNumberOfLoops );
673            if ( loopCounter >= maxNumberOfLoops ) {
674              z = 0.5*(zmin + zmax);  // Just a value between zmin and zmax, no physics considerations at all! 
675            }
676      return z;
677         }
678 
679   if (std::abs(PDGEncodingOfDecayParton) > 1000)     
680   {
681                 G4double an = 2.5;
682                 an +=(sqr(Px)+sqr(Py))/sqr(GeV)-0.5;
683                 z=zmin + (zmax-zmin)*G4Pow::GetInstance()->powA(G4UniformRand(),1./an);
684                 if( PDGEncodingOfDecayParton > 3000 ) z=zmin+zmax-z;
685   }
686 
687   return z;
688 }
689 
690 //----------------------------------------------------------------------------------------------------------
691 
692 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
693                                             G4KineticTrackVector * LeftVector,
694                                             G4KineticTrackVector * RightVector)
695 {
696   //... perform last cluster decay
697         SetMinimalStringMass( string);
698         if ( MinimalStringMass < 0.) return false;
699         #ifdef debug_LUNDfragmentation
700         G4cout<<G4endl<<"Split last-----------------------------------------"<<G4endl;
701         G4cout<<"MinimalStringMass "<<MinimalStringMass<<G4endl;
702         G4cout<<"Left  "<<string->GetLeftParton()->GetPDGEncoding()<<" "<<string->GetPleft()<<G4endl;
703         G4cout<<"Right "<<string->GetRightParton()->GetPDGEncoding()<<" "<<string->GetPright()<<G4endl;
704         G4cout<<"String4mom "<<string->GetPstring()<<" "<<string->GetPstring().mag()<<G4endl;
705         #endif
706 
707         G4LorentzVector Str4Mom=string->Get4Momentum();
708         G4LorentzRotation toCms(-1*Str4Mom.boostVector());
709         G4LorentzVector Pleft = toCms * string->GetPleft();
710         toCms.rotateZ(-1*Pleft.phi());
711         toCms.rotateY(-1*Pleft.theta());
712   
713         G4LorentzRotation toObserverFrame= toCms.inverse();
714 
715   G4double StringMass=string->Mass();
716 
717   G4ParticleDefinition * LeftHadron(0), * RightHadron(0);
718 
719         NumberOf_FS=0;
720   for (G4int i=0; i<350; i++) {FS_Weight[i]=0.;}
721 
722   G4int sampledState = 0;
723 
724         #ifdef debug_LUNDfragmentation
725         G4cout<<"StrMass "<<StringMass<<" q's "
726               <<string->GetLeftParton()->GetParticleName()<<" "
727               <<string->GetRightParton()->GetParticleName()<<G4endl;
728         #endif
729 
730   string->SetLeftPartonStable(); // to query quark contents..
731 
732   if (string->IsAFourQuarkString() )
733   {
734           G4int IDleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
735           G4int IDright=std::abs(string->GetRightParton()->GetPDGEncoding());
736 
737           if ( (IDleft > 3000) || (IDright > 3000) ) {
738             if ( ! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
739             {
740               return false;
741             } 
742           } else {
743     // The string is qq-qqbar type. Diquarks are on the string ends
744           if (StringMass-MinimalStringMass < 0.)
745     {
746       if (! Diquark_AntiDiquark_belowThreshold_lastSplitting(string, LeftHadron, RightHadron) )
747                         {
748         return false;
749                         }
750     } else
751     {
752       Diquark_AntiDiquark_aboveThreshold_lastSplitting(string, LeftHadron, RightHadron);
753 
754       if (NumberOf_FS == 0) return false;
755 
756                         sampledState = SampleState();
757       if (string->GetLeftParton()->GetPDGEncoding() < 0)
758       {
759         LeftHadron =FS_LeftHadron[sampledState];
760         RightHadron=FS_RightHadron[sampledState];
761       } else
762       {
763         LeftHadron =FS_RightHadron[sampledState];
764         RightHadron=FS_LeftHadron[sampledState];
765       }
766     }
767           }  // ID > 3300
768         } else {
769     if (string->DecayIsQuark() && string->StableIsQuark() )
770     {       //... there are quarks on cluster ends
771                         #ifdef debug_LUNDfragmentation
772                         G4cout<<"Q Q string LastSplit"<<G4endl;
773                         #endif
774 
775       Quark_AntiQuark_lastSplitting(string, LeftHadron, RightHadron);
776     
777       if (NumberOf_FS == 0) return false;
778                   sampledState = SampleState();
779       if (string->GetLeftParton()->GetPDGEncoding() < 0)
780       {
781         LeftHadron =FS_RightHadron[sampledState];
782         RightHadron=FS_LeftHadron[sampledState];
783       } else
784       {
785         LeftHadron =FS_LeftHadron[sampledState];
786         RightHadron=FS_RightHadron[sampledState];
787       }
788     } else 
789     {       //... there is a Diquark on one of the cluster ends
790                         #ifdef debug_LUNDfragmentation
791                         G4cout<<"DiQ Q string Last Split"<<G4endl;
792                         #endif
793 
794       Quark_Diquark_lastSplitting(string, LeftHadron, RightHadron);
795 
796       if (NumberOf_FS == 0) return false;                           
797                   sampledState = SampleState();
798 
799       if (string->GetLeftParton()->GetParticleSubType() == "quark")
800       {
801         LeftHadron =FS_LeftHadron[sampledState];
802         RightHadron=FS_RightHadron[sampledState];
803       } else 
804       {
805         LeftHadron =FS_RightHadron[sampledState];
806         RightHadron=FS_LeftHadron[sampledState];
807       }     
808     }
809 
810   }
811 
812         #ifdef debug_LUNDfragmentation
813         G4cout<<"Sampled hadrons: "<<LeftHadron->GetParticleName()<<" "<<RightHadron->GetParticleName()<<G4endl;
814         #endif
815 
816   G4LorentzVector P_left  =string->GetPleft(), P_right = string->GetPright();
817 
818   G4LorentzVector  LeftMom, RightMom;
819   G4ThreeVector    Pos;
820 
821   Sample4Momentum(&LeftMom,  LeftHadron->GetPDGMass(),
822       &RightMom, RightHadron->GetPDGMass(),
823       StringMass);
824 
825         // Sample4Momentum ascribes LeftMom.pz() along positive Z axis for baryons in many cases.
826         // It must be negative in case when the system is moving against Z axis.
827 
828   if (!(string->DecayIsQuark() && string->StableIsQuark() ))
829   { // Only for qq - q, q - qq, and qq - qqbar -------------------
830 
831           if ( G4UniformRand() <= 0.5 )
832     {
833       if (P_left.z() <= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
834     } 
835     else  
836     {
837       if (P_right.z() >= 0.) {G4LorentzVector tmp = LeftMom; LeftMom=RightMom; RightMom=tmp;}
838     }
839   }
840 
841   LeftMom *=toObserverFrame;
842   RightMom*=toObserverFrame;
843 
844   LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
845   RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
846 
847   string->LorentzRotate(toObserverFrame);
848   return true;
849 }
850 
851 //----------------------------------------------------------------------------------------
852 
853 G4bool G4LundStringFragmentation::
854 Diquark_AntiDiquark_belowThreshold_lastSplitting(G4FragmentingString * & string,
855                                                  G4ParticleDefinition * & LeftHadron,
856                                                  G4ParticleDefinition * & RightHadron)
857 {
858   G4double StringMass   = string->Mass();
859 
860   G4int cClusterInterrupt = 0;
861         G4bool isOK = false;
862   do
863   {
864     G4int LeftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
865     G4int LeftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
866 
867     G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
868     G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
869 
870     if (G4UniformRand()<0.5)
871     {
872       LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
873                   FindParticle(RightQuark1));
874       RightHadron= (LeftHadron == nullptr) ? nullptr :
875                                                       hadronizer->Build(FindParticle( LeftQuark2),
876                   FindParticle(RightQuark2));
877     } else
878     {
879       LeftHadron =hadronizer->Build(FindParticle( LeftQuark1),
880                   FindParticle(RightQuark2));
881       RightHadron=(LeftHadron == nullptr) ? nullptr :
882                                     hadronizer->Build(FindParticle( LeftQuark2),
883                   FindParticle(RightQuark1));
884     }
885 
886     isOK = (LeftHadron != nullptr) && (RightHadron != nullptr);
887 
888     if(isOK) { isOK = (StringMass > LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()); }
889     ++cClusterInterrupt;
890     //... repeat procedure, if mass of cluster is too low to produce hadrons
891     //... ClusterMassCut = 0.15*GeV model parameter
892   }
893   while (isOK == false && cClusterInterrupt < ClusterLoopInterrupt);
894   /* Loop checking, 07.08.2015, A.Ribon */
895     return isOK;
896 }
897 
898 //----------------------------------------------------------------------------------------
899 
900 G4bool G4LundStringFragmentation::
901 Diquark_AntiDiquark_aboveThreshold_lastSplitting(G4FragmentingString * & string,
902                                                  G4ParticleDefinition * & LeftHadron,
903                                                  G4ParticleDefinition * & RightHadron)
904 {
905   // StringMass-MinimalStringMass > 0. Creation of 2 baryons is possible ----
906 
907   G4double StringMass   = string->Mass();
908   G4double StringMassSqr= sqr(StringMass); 
909   G4ParticleDefinition * Di_Quark;
910   G4ParticleDefinition * Anti_Di_Quark;
911 
912   if (string->GetLeftParton()->GetPDGEncoding() < 0)
913   {
914     Anti_Di_Quark   =string->GetLeftParton();
915     Di_Quark=string->GetRightParton();
916   } else
917   {
918     Anti_Di_Quark   =string->GetRightParton();
919     Di_Quark=string->GetLeftParton();
920   }
921 
922   G4int IDAnti_di_quark    =Anti_Di_Quark->GetPDGEncoding();
923   G4int AbsIDAnti_di_quark =std::abs(IDAnti_di_quark);
924   G4int IDdi_quark         =Di_Quark->GetPDGEncoding();
925   G4int AbsIDdi_quark      =std::abs(IDdi_quark);
926 
927   G4int ADi_q1=AbsIDAnti_di_quark/1000;
928   G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000)/100;
929 
930   G4int Di_q1=AbsIDdi_quark/1000;
931   G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
932 
933   NumberOf_FS=0;
934   for (G4int ProdQ=1; ProdQ < 6; ProdQ++)
935   {
936     G4int StateADiQ=0;
937                 const G4int maxNumberOfLoops = 1000;
938                 G4int loopCounter = 0;
939     do  // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
940     {
941       LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(
942               -Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]);
943 
944       if (LeftHadron == NULL) continue;
945       G4double LeftHadronMass=LeftHadron->GetPDGMass();
946 
947       G4int StateDiQ=0;
948                         const G4int maxNumberOfInternalLoops = 1000;
949                         G4int internalLoopCounter = 0;
950       do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
951       {
952         RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(
953                 +Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
954 
955         if (RightHadron == NULL) continue;
956         G4double RightHadronMass=RightHadron->GetPDGMass();
957 
958         if (StringMass > LeftHadronMass + RightHadronMass)
959         {
960                                         if ( NumberOf_FS > 349 ) {
961                                           G4ExceptionDescription ed;
962                                           ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
963                                           G4Exception( "G4LundStringFragmentation::Diquark_AntiDiquark_aboveThreshold_lastSplitting ",
964                                                        "HAD_LUND_001", JustWarning, ed );
965                                           NumberOf_FS = 349;
966                                         }
967 
968           G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
969                 sqr(RightHadronMass));
970           //FS_Psqr=1.;
971           FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*FS_Psqr*
972                      BaryonWeight[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]*
973                      BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
974                      Prob_QQbar[ProdQ-1];
975 
976           FS_LeftHadron[NumberOf_FS] = LeftHadron;
977           FS_RightHadron[NumberOf_FS]= RightHadron;
978           NumberOf_FS++;
979         } // End of if (StringMass > LeftHadronMass + RightHadronMass)
980 
981         StateDiQ++;
982 
983       } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) && 
984                                  ++internalLoopCounter < maxNumberOfInternalLoops );
985                         if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
986                           return false;
987                         }
988  
989       StateADiQ++;
990     } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ-1][StateADiQ]!=0) &&
991                          ++loopCounter < maxNumberOfLoops ); 
992                 if ( loopCounter >= maxNumberOfLoops ) {
993                   return false;
994                 }
995   } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
996 
997     return true;
998 }
999 
1000 //----------------------------------------------------------------------------------------
1001 
1002 G4bool G4LundStringFragmentation::Quark_Diquark_lastSplitting(G4FragmentingString * & string,
1003                                                               G4ParticleDefinition * & LeftHadron,
1004                                                               G4ParticleDefinition * & RightHadron)
1005 {
1006   G4double StringMass   = string->Mass();
1007   G4double StringMassSqr= sqr(StringMass);
1008 
1009   G4ParticleDefinition * Di_Quark;
1010   G4ParticleDefinition * Quark;
1011 
1012   if (string->GetLeftParton()->GetParticleSubType()== "quark")
1013   {
1014     Quark   =string->GetLeftParton();
1015     Di_Quark=string->GetRightParton();
1016   } else
1017   {
1018     Quark   =string->GetRightParton();
1019     Di_Quark=string->GetLeftParton();
1020   }
1021 
1022   G4int IDquark        =Quark->GetPDGEncoding();
1023   G4int AbsIDquark     =std::abs(IDquark);
1024   G4int IDdi_quark   =Di_Quark->GetPDGEncoding();
1025   G4int AbsIDdi_quark=std::abs(IDdi_quark);
1026   G4int Di_q1=AbsIDdi_quark/1000;
1027   G4int Di_q2=(AbsIDdi_quark-Di_q1*1000)/100;
1028   G4int SignDiQ= 1;
1029   if (IDdi_quark < 0) SignDiQ=-1;
1030 
1031   NumberOf_FS=0;
1032   for (G4int ProdQ=1; ProdQ < 4; ProdQ++)  // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar 
1033   {                                        // (as last splitting, do not consider c-cbar and b-bbar cases)
1034     G4int SignQ;
1035     if (IDquark > 0)
1036     {
1037             SignQ=-1;
1038             if (IDquark == 2)                   SignQ= 1;
1039       if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0
1040       if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar
1041             if (IDquark == 4)                   SignQ= 1; // D+, D0, Ds+
1042       if (IDquark == 5)                   SignQ=-1; // B-, anti_B0, anti_Bs0
1043     } else
1044     {
1045       SignQ= 1;
1046       if (IDquark == -2)                  SignQ=-1;
1047       if ((IDquark ==-1) && (ProdQ == 3)) SignQ=-1; // K0bar
1048       if ((IDquark ==-3) && (ProdQ == 1)) SignQ= 1; // K0
1049       if (IDquark == -4)                  SignQ=-1; // D-, anti_D0, anti_Ds+
1050       if (IDquark == -5)                  SignQ= 1; // B+, B0, Bs0
1051     }
1052 
1053     if (AbsIDquark == ProdQ)            SignQ= 1;
1054 
1055     G4int StateQ=0;
1056                 const G4int maxNumberOfLoops = 1000;
1057                 G4int loopCounter = 0;
1058     do  // while(Meson[AbsIDquark-1][ProdQ-1][StateQ]<>0);
1059     {
1060       LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1061               Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1062       if (LeftHadron == NULL) continue;
1063       G4double LeftHadronMass=LeftHadron->GetPDGMass();
1064 
1065       G4int StateDiQ=0;
1066                         const G4int maxNumberOfInternalLoops = 1000;
1067                         G4int internalLoopCounter = 0;
1068       do // while(Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]<>0);
1069       {
1070         RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignDiQ*
1071                 Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]);
1072         if (RightHadron == NULL) continue;
1073         G4double RightHadronMass=RightHadron->GetPDGMass();
1074 
1075         if (StringMass > LeftHadronMass + RightHadronMass)
1076         {
1077                                         if ( NumberOf_FS > 349 ) {
1078                                           G4ExceptionDescription ed;
1079                                           ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1080                                           G4Exception( "G4LundStringFragmentation::Quark_Diquark_lastSplitting ",
1081                                                        "HAD_LUND_002", JustWarning, ed );
1082                                           NumberOf_FS = 349;
1083                                         }
1084 
1085           G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1086                 sqr(RightHadronMass));
1087           FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1088                      MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1089                      BaryonWeight[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]*
1090                      Prob_QQbar[ProdQ-1];
1091 
1092           FS_LeftHadron[NumberOf_FS] = LeftHadron;
1093           FS_RightHadron[NumberOf_FS]= RightHadron;
1094 
1095           NumberOf_FS++;
1096         } // End of if (StringMass > LeftHadronMass + RightHadronMass)
1097 
1098         StateDiQ++;
1099 
1100       } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ-1][StateDiQ]!=0) &&
1101                                  ++internalLoopCounter < maxNumberOfInternalLoops );
1102                         if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1103                           return false;
1104                         }
1105 
1106       StateQ++;
1107     } while( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1108                          ++loopCounter < maxNumberOfLoops );  /* Loop checking, 07.08.2015, A.Ribon */
1109 
1110                   if ( loopCounter >= maxNumberOfLoops ) {
1111                     return false;
1112                   }
1113   }
1114 
1115   return true;
1116 }
1117 
1118 //----------------------------------------------------------------------------------------
1119 
1120 G4bool G4LundStringFragmentation::Quark_AntiQuark_lastSplitting(G4FragmentingString * & string,
1121                                                                 G4ParticleDefinition * & LeftHadron,
1122                                                                 G4ParticleDefinition * & RightHadron)
1123 {
1124   G4double StringMass   = string->Mass();
1125   G4double StringMassSqr= sqr(StringMass);
1126 
1127   G4ParticleDefinition * Quark;
1128   G4ParticleDefinition * Anti_Quark;
1129 
1130   if (string->GetLeftParton()->GetPDGEncoding()>0)
1131   {
1132     Quark     =string->GetLeftParton();
1133     Anti_Quark=string->GetRightParton();
1134   } else
1135   {
1136     Quark     =string->GetRightParton();
1137     Anti_Quark=string->GetLeftParton();
1138   }
1139 
1140   G4int IDquark         =Quark->GetPDGEncoding();
1141   G4int AbsIDquark      =std::abs(IDquark);
1142         G4int QuarkCharge     =Qcharge[IDquark-1];
1143 
1144   G4int IDanti_quark    =Anti_Quark->GetPDGEncoding();
1145   G4int AbsIDanti_quark =std::abs(IDanti_quark);
1146         G4int AntiQuarkCharge =-Qcharge[AbsIDanti_quark-1];
1147 
1148         G4int LeftHadronCharge(0), RightHadronCharge(0);
1149 
1150         //G4cout<<"Q Qbar "<<IDquark<<" "<<IDanti_quark<<G4endl;
1151 
1152   NumberOf_FS=0;
1153   for (G4int ProdQ=1; ProdQ < 4; ProdQ++)  // Loop over quark-antiquark cases: u-ubar, d-dbar, s-sbar 
1154   {                                        // (as last splitting, do not consider c-cbar and b-bbar cases)
1155                 //G4cout<<"NumberOf_FS ProdQ "<<NumberOf_FS<<" "<<ProdQ<<G4endl;
1156     LeftHadronCharge = QuarkCharge - Qcharge[ProdQ-1];
1157     G4int SignQ = LeftHadronCharge/3; if (SignQ == 0) SignQ = 1;
1158 
1159     if ((IDquark == 1) && (ProdQ == 3)) SignQ= 1; // K0       (d,sbar)
1160     if ((IDquark == 3) && (ProdQ == 1)) SignQ=-1; // K0bar    (s,dbar)
1161                 if ((IDquark == 4) && (ProdQ == 2)) SignQ= 1; // D0       (c,ubar)
1162                 if ((IDquark == 5) && (ProdQ == 1)) SignQ=-1; // anti_B0  (b,dbar)
1163                 if ((IDquark == 5) && (ProdQ == 3)) SignQ=-1; // anti_Bs0 (b,sbar)
1164     
1165                 RightHadronCharge = AntiQuarkCharge + Qcharge[ProdQ-1];
1166     G4int SignAQ = RightHadronCharge/3; if (SignAQ == 0) SignAQ = 1;
1167 
1168     if ((IDanti_quark ==-1) && (ProdQ == 3)) SignAQ=-1; // K0bar   (dbar,s)
1169     if ((IDanti_quark ==-3) && (ProdQ == 1)) SignAQ= 1; // K0      (sbar,d)
1170     if ((IDanti_quark ==-4) && (ProdQ == 2)) SignAQ=-1; // anti_D0 (cbar,u)
1171     if ((IDanti_quark ==-5) && (ProdQ == 1)) SignAQ= 1; // B0      (bbar,d)
1172     if ((IDanti_quark ==-5) && (ProdQ == 3)) SignAQ= 1; // Bs0     (bbar,s)
1173 
1174                 //G4cout<<"ProQ signs "<<ProdQ<<" "<<SignQ<<" "<<SignAQ<<G4endl;
1175 
1176     G4int StateQ=0;
1177                 const G4int maxNumberOfLoops = 1000;
1178                 G4int loopCounter = 0;
1179     do
1180     {
1181                         //G4cout<<"[AbsIDquark-1][ProdQ-1][StateQ "<<AbsIDquark-1<<" "
1182                         //<<ProdQ-1<<" "<<StateQ<<" "<<SignQ*Meson[AbsIDquark-1][ProdQ-1][StateQ]<<G4endl;
1183       LeftHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignQ*
1184                    Meson[AbsIDquark-1][ProdQ-1][StateQ]);
1185                         //G4cout<<"LeftHadron "<<LeftHadron<<G4endl; 
1186       if (LeftHadron == NULL) { StateQ++; continue; }
1187                         //G4cout<<"LeftHadron "<<LeftHadron->GetParticleName()<<G4endl;
1188       G4double LeftHadronMass=LeftHadron->GetPDGMass();
1189 
1190       G4int StateAQ=0;
1191                         const G4int maxNumberOfInternalLoops = 1000;
1192                         G4int internalLoopCounter = 0;
1193       do
1194       {
1195                                 //G4cout<<"           [AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<AbsIDanti_quark-1<<" "
1196                                 //      <<ProdQ-1<<" "<<StateAQ<<" "<<SignAQ*Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<G4endl;
1197         RightHadron=G4ParticleTable::GetParticleTable()->FindParticle(SignAQ*
1198                 Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]);
1199                                 //G4cout<<"RightHadron "<<RightHadron<<G4endl;
1200         if(RightHadron == NULL) { StateAQ++; continue; }
1201                                 //G4cout<<"RightHadron "<<RightHadron->GetParticleName()<<G4endl;
1202         G4double RightHadronMass=RightHadron->GetPDGMass();
1203 
1204         if (StringMass > LeftHadronMass + RightHadronMass)
1205         {
1206                                         if ( NumberOf_FS > 349 ) {
1207                                           G4ExceptionDescription ed;
1208                                           ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1209                                           G4Exception( "G4LundStringFragmentation::Quark_AntiQuark_lastSplitting ",
1210                                                        "HAD_LUND_003", JustWarning, ed );
1211                                           NumberOf_FS = 349;
1212                                         }
1213 
1214                                         G4double FS_Psqr=lambda(StringMassSqr,sqr(LeftHadronMass),
1215                 sqr(RightHadronMass));
1216           //FS_Psqr=1.;
1217           FS_Weight[NumberOf_FS]=std::sqrt(FS_Psqr)*
1218                      MesonWeight[AbsIDquark-1][ProdQ-1][StateQ]*
1219                      MesonWeight[AbsIDanti_quark-1][ProdQ-1][StateAQ]*
1220                      Prob_QQbar[ProdQ-1];
1221           if (string->GetLeftParton()->GetPDGEncoding()>0)
1222           {
1223             FS_LeftHadron[NumberOf_FS] = RightHadron;
1224             FS_RightHadron[NumberOf_FS]= LeftHadron;
1225           } else
1226           {
1227             FS_LeftHadron[NumberOf_FS] = LeftHadron;
1228             FS_RightHadron[NumberOf_FS]= RightHadron;
1229           }
1230 
1231           NumberOf_FS++;
1232 
1233         }
1234 
1235         StateAQ++;
1236                                 //G4cout<<"               StateAQ Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ] "<<StateAQ<<" "
1237                                 //      <<Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]<<" "<<internalLoopCounter<<G4endl;
1238       } while ( (Meson[AbsIDanti_quark-1][ProdQ-1][StateAQ]!=0) &&
1239                                   ++internalLoopCounter < maxNumberOfInternalLoops );
1240                           if ( internalLoopCounter >= maxNumberOfInternalLoops ) {
1241                             return false;
1242                           }
1243 
1244       StateQ++;
1245                         //G4cout<<"StateQ Meson[AbsIDquark-1][ProdQ-1][StateQ] "<<StateQ<<" "
1246                         //      <<Meson[AbsIDquark-1][ProdQ-1][StateQ]<<" "<<loopCounter<<G4endl;
1247 
1248     } while ( (Meson[AbsIDquark-1][ProdQ-1][StateQ]!=0) &&
1249                          ++loopCounter < maxNumberOfLoops );
1250                   if ( loopCounter >= maxNumberOfLoops ) {
1251                     return false;
1252                   }
1253   } // End of for (G4int ProdQ=1; ProdQ < 4; ProdQ++)
1254 
1255   return true;
1256 }
1257 
1258 //----------------------------------------------------------------------------------------------------------
1259 
1260 G4int G4LundStringFragmentation::SampleState(void) 
1261 {
1262         if ( NumberOf_FS > 349 ) {
1263           G4ExceptionDescription ed;
1264           ed << " NumberOf_FS exceeds its limit: NumberOf_FS=" << NumberOf_FS << G4endl;
1265           G4Exception( "G4LundStringFragmentation::SampleState ", "HAD_LUND_004", JustWarning, ed );
1266           NumberOf_FS = 349;
1267         }
1268 
1269   G4double SumWeights=0.;
1270   for (G4int i=0; i<NumberOf_FS; i++) {SumWeights+=FS_Weight[i];}
1271 
1272   G4double ksi=G4UniformRand();
1273   G4double Sum=0.;
1274   G4int indexPosition = 0;
1275 
1276   for (G4int i=0; i<NumberOf_FS; i++)
1277   {
1278     Sum+=(FS_Weight[i]/SumWeights);
1279     indexPosition=i;
1280     if (Sum >= ksi) break;
1281   }
1282   return indexPosition;
1283 }
1284 
1285 //----------------------------------------------------------------------------------------------------------
1286 
1287 void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector*     Mom, G4double     Mass, 
1288                                                 G4LorentzVector* AntiMom, G4double AntiMass, 
1289                                                 G4double InitialMass) 
1290 {
1291   // ------ Sampling of momenta of 2 last produced hadrons --------------------
1292   G4ThreeVector Pt;
1293   G4double MassMt, AntiMassMt;
1294   G4double AvailablePz, AvailablePz2;
1295         #ifdef debug_LUNDfragmentation
1296         G4cout<<"Sampling of momenta of 2 last produced hadrons ----------------"<<G4endl;
1297         G4cout<<"Init Mass "<<InitialMass<<" FirstM "<<Mass<<" SecondM "<<AntiMass<<" ProbIsotropy "<<G4endl;  
1298         #endif
1299 
1300   G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -   
1301       sqr(2.*Mass*AntiMass);
1302   G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
1303 
1304         const G4int maxNumberOfLoops = 1000;
1305   G4double SigmaQTw = SigmaQT;
1306   if ( Mass > 930. || AntiMass > 930. ) {
1307     SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) );
1308         }
1309         if ( Mass < 930. && AntiMass < 930. ) {}     // q-qbar string
1310         if ( ( Mass < 930. && AntiMass > 930. ) ||
1311        ( Mass > 930. && AntiMass < 930. ) ) {  // q-di_q string
1312           //SigmaQT = -1.;  // isotropical decay
1313         }
1314         if ( Mass > 930. && AntiMass > 930. ) {      // qq-qqbar string
1315           SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMass)/InitialMass ) );
1316         }
1317 
1318         G4int loopCounter = 0;
1319   do
1320   {
1321     Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4double Pt2=Pt.mag2();
1322     MassMt    = std::sqrt(    Mass *     Mass + Pt2);
1323     AntiMassMt= std::sqrt(AntiMass * AntiMass + Pt2);
1324   }
1325   while ( (InitialMass < MassMt + AntiMassMt) && ++loopCounter < maxNumberOfLoops ); 
1326 
1327         SigmaQT = SigmaQTw;
1328 
1329   AvailablePz2= sqr(InitialMass*InitialMass - sqr(MassMt) - sqr(AntiMassMt)) -
1330         4.*sqr(MassMt*AntiMassMt);
1331 
1332   AvailablePz2 /=(4.*InitialMass*InitialMass);
1333   AvailablePz = std::sqrt(AvailablePz2);
1334 
1335   G4double Px=Pt.getX();
1336   G4double Py=Pt.getY();
1337 
1338   Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
1339   Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz2));
1340 
1341   AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
1342   AntiMom->setE (std::sqrt(sqr(AntiMassMt)+AvailablePz2));
1343 
1344         #ifdef debug_LUNDfragmentation
1345         G4cout<<"Fmass Mom "<<Mom->getX()<<" "<<Mom->getY()<<" "<<Mom->getZ()<<" "<<Mom->getT()<<G4endl;
1346         G4cout<<"Smass Mom "<<AntiMom->getX()<<" "<<AntiMom->getY()<<" "<<AntiMom->getZ()
1347               <<" "<<AntiMom->getT()<<G4endl;
1348         #endif
1349 }
1350 
1351 //------------------------------------------------------------------------
1352 
1353 G4double G4LundStringFragmentation::lambda(G4double S, G4double m1_Sqr, G4double m2_Sqr)
1354 { 
1355   G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4.*m1_Sqr*m2_Sqr;
1356   return lam;
1357 }
1358 
1359 // --------------------------------------------------------------
1360 G4LundStringFragmentation::~G4LundStringFragmentation()
1361 {}
1362 
1363