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 ]

Diff markup

Differences between /processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc (Version 11.3.0) and /processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc (Version 6.2.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 //                                                 23 //
                                                   >>  24 // $Id: G4LundStringFragmentation.cc,v 1.2 2003/11/03 17:54:53 hpw Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-06-00-patch-01 $
 27 //                                                 26 //
 28 // -------------------------------------------     27 // -----------------------------------------------------------------------------
 29 //      GEANT 4 class implementation file          28 //      GEANT 4 class implementation file
 30 //                                                 29 //
 31 //      History: first implementation, Maxim K     30 //      History: first implementation, Maxim Komogorov, 10-Jul-1998
 32 // -------------------------------------------     31 // -----------------------------------------------------------------------------
 33 #include "G4LundStringFragmentation.hh"            32 #include "G4LundStringFragmentation.hh"
 34 #include "G4PhysicalConstants.hh"              << 
 35 #include "G4SystemOfUnits.hh"                  << 
 36 #include "Randomize.hh"                            33 #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                                                    34 
 46 // Class G4LundStringFragmentation                 35 // Class G4LundStringFragmentation 
 47 //******************************************** <<  36 //****************************************************************************************
 48                                                    37 
 49 G4LundStringFragmentation::G4LundStringFragmen     38 G4LundStringFragmentation::G4LundStringFragmentation()
 50   : G4VLongitudinalStringDecay("LundStringFrag <<  39    {
 51 {                                              << 
 52     SetMassCut(210.*MeV);   //  Mpi + Delta    << 
 53                             // For ProduceOneH << 
 54                             // that no one pi- << 
 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 << 
 65     // set the non-zero probabilities for c-cb << 
 66     // else set them to 0.0. If these probabil << 
 67     // hadrons can't/can be created during the << 
 68     // (i.e. not heavy) projectile hadron nucl << 
 69     if ( G4HadronicParameters::Instance()->Ena << 
 70       SetProbCCbar(0.0002);  // According to O << 
 71       SetProbBBbar(5.0e-5);  // According to O << 
 72     } else {                                   << 
 73       SetProbCCbar(0.0);                       << 
 74       SetProbBBbar(0.0);                       << 
 75     }                                          << 
 76                                                << 
 77     SetMinMasses();  // For treating of small  << 
 78 }                                              << 
 79                                                << 
 80 //-------------------------------------------- << 
 81                                                << 
 82 G4KineticTrackVector* G4LundStringFragmentatio << 
 83 {                                              << 
 84   // Can no longer modify Parameters for Fragm << 
 85                                                << 
 86   PastInitPhase=true;                          << 
 87                                                << 
 88   G4FragmentingString  aString(theString);     << 
 89   SetMinimalStringMass(&aString);              << 
 90                                                << 
 91         #ifdef debug_LUNDfragmentation         << 
 92         G4cout<<G4endl<<"LUND StringFragmentat << 
 93         G4cout<<G4endl<<"LUND StringFragm: Str << 
 94                       <<theString.Get4Momentum << 
 95                       <<"4Mom "<<theString.Get << 
 96                       <<"--------------------- << 
 97         G4cout<<"String ends Direct "<<theStri << 
 98                                      <<theStri << 
 99                                      <<theStri << 
100         G4cout<<"Left  mom "<<theString.GetLef << 
101         G4cout<<"Right mom "<<theString.GetRig << 
102         G4cout<<"Check for Fragmentation "<<G4 << 
103         #endif                                 << 
104                                                << 
105   G4KineticTrackVector * LeftVector(0);        << 
106                                                << 
107   if (!aString.IsAFourQuarkString() && !IsItFr << 
108   {                                            << 
109                 #ifdef debug_LUNDfragmentation << 
110                 G4cout<<"Non fragmentable - th << 
111                 #endif                         << 
112                 // SetMassCut(210.*MeV);  // F << 
113                                           // t << 
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)->SetForm << 
125             LeftVector->operator[](0)->SetPosi << 
126                   }                            << 
127       if (LeftVector->size() > 1)              << 
128                   {                            << 
129             // 2 hadrons created from qq-qqbar << 
130       LeftVector->operator[](1)->SetFormationT << 
131       LeftVector->operator[](1)->SetPosition(t << 
132       }                                        << 
133     }                                          << 
134     return LeftVector;                         << 
135   }                                            << 
136                                                << 
137         #ifdef debug_LUNDfragmentation         << 
138         G4cout<<"The string will be fragmented << 
139         #endif                                 << 
140                                                << 
141   // The string can fragment. At least two par << 
142              LeftVector =new G4KineticTrackVec << 
143   G4KineticTrackVector * RightVector=new G4Kin << 
144                                                << 
145         G4bool success = Loop_toFragmentString << 
146                                                << 
147   if ( ! success )                             << 
148   {                                            << 
149     std::for_each(LeftVector->begin(), LeftVec << 
150     LeftVector->clear();                       << 
151     std::for_each(RightVector->begin(), RightV << 
152     delete RightVector;                        << 
153     return LeftVector;                         << 
154   }                                            << 
155                                                << 
156   // Join Left- and RightVector into LeftVecto << 
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::IsItFragment << 
170 {                                              << 
171   SetMinimalStringMass(string);                << 
172         //G4cout<<"MinM StrM "<<MinimalStringM << 
173                                                << 
174   return std::abs(MinimalStringMass) < string- << 
175                                                << 
176         //MinimalStringMass is negative and la << 
177 }                                              << 
178                                                << 
179 //-------------------------------------------- << 
180                                                << 
181 G4bool G4LundStringFragmentation::Loop_toFragm << 
182                                                << 
183                                                << 
184 {                                              << 
185         #ifdef debug_LUNDfragmentation         << 
186         G4cout<<"Loop_toFrag "<<theString.GetL << 
187                               <<theString.GetL << 
188               <<"            "<<theString.GetR << 
189                               <<theString.GetR << 
190               <<"Direction   "<<theString.GetD << 
191         #endif                                 << 
192                                                << 
193         G4LorentzRotation toCmsI, toObserverFr << 
194                                                << 
195   G4bool final_success=false;                  << 
196   G4bool inner_success=true;                   << 
197                                                << 
198   G4int attempt=0;                             << 
199                                                << 
200   while ( ! final_success && attempt++ < Strin << 
201   {       // If the string fragmentation does  << 
202           // repeat the fragmentation.         << 
203                                                << 
204                 G4FragmentingString *currentSt << 
205                 toCmsI = currentString->Transf << 
206                 toObserverFrameI = toCmsI.inve << 
207                                                << 
208     G4LorentzRotation toCms, toObserverFrame;  << 
209                                                << 
210     //G4cout<<"Main loop start whilecounter "< << 
211                                                << 
212     // Cleaning up the previously produced had << 
213     std::for_each(LeftVector->begin(), LeftVec << 
214     LeftVector->clear();                       << 
215     std::for_each(RightVector->begin(), RightV << 
216     RightVector->clear();                      << 
217                                                << 
218     // Main fragmentation loop until the strin << 
219     inner_success=true;  // set false on failu << 
220                 const G4int maxNumberOfLoops = << 
221                 G4int loopCounter = -1;        << 
222                                                << 
223     while ( (! StopFragmenting(currentString)) << 
224     {       // Split current string into hadro << 
225                         #ifdef debug_LUNDfragm << 
226                         G4cout<<"The string wi << 
227                         //G4cout<<"1 "<<curren << 
228                         #endif                 << 
229       G4FragmentingString *newString=0;  // us << 
230                                                << 
231       toCms=currentString->TransformToAlignedC << 
232                         toObserverFrame= toCms << 
233                                                << 
234                         #ifdef debug_LUNDfragm << 
235                         //G4cout<<"CMS Left  m << 
236                         //G4cout<<"CMS Right m << 
237                         //G4cout<<"CMS String  << 
238                         #endif                 << 
239                                                << 
240       G4KineticTrack * Hadron=Splitup(currentS << 
241                                                << 
242       if ( Hadron != 0 )  // Store the hadron  << 
243       {                                        << 
244                                 #ifdef debug_L << 
245                                 G4cout<<"Hadro << 
246                                 //G4cout<<"2 " << 
247                                 #endif         << 
248                                                << 
249         Hadron->Set4Momentum(toObserverFrame*H << 
250                                                << 
251         G4double TimeOftheStringCreation=theSt << 
252         G4ThreeVector PositionOftheStringCreat << 
253                                                << 
254         G4LorentzVector Coordinate(Hadron->Get << 
255         G4LorentzVector Momentum = toObserverF << 
256         Hadron->SetFormationTime(TimeOftheStri << 
257         G4ThreeVector aPosition(Momentum.vect( << 
258         Hadron->SetPosition(PositionOftheStrin << 
259                                                << 
260                                 // Open to pro << 
261         if ( currentString->GetDecayDirection( << 
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 ) del << 
272                         }                      << 
273                                                << 
274                         currentString->Lorentz << 
275     };                                         << 
276                                                << 
277                 if ( loopCounter >= maxNumberO << 
278                   inner_success=false;         << 
279                 }                              << 
280                                                << 
281     // Split remaining string into 2 final had << 
282                 #ifdef debug_LUNDfragmentation << 
283                 if (inner_success) G4cout<<"Sp << 
284                 #endif                         << 
285                                                << 
286     if ( inner_success && SplitLast(currentStr << 
287     {                                          << 
288       final_success = true;                    << 
289     }                                          << 
290                                                << 
291     delete currentString;                      << 
292   }  // End of the loop where we try to fragme << 
293                                                << 
294         G4int sign = +1;                       << 
295         if ( theString.GetDirection() < 0 ) si << 
296         for ( unsigned int hadronI = 0; hadron << 
297            G4LorentzVector Tmp = LeftVector->o << 
298            Tmp.setZ(sign*Tmp.getZ());          << 
299            Tmp *= toObserverFrameI;            << 
300            LeftVector->operator[](hadronI)->Se << 
301         }                                      << 
302         for ( unsigned int hadronI = 0; hadron << 
303            G4LorentzVector Tmp = RightVector-> << 
304            Tmp.setZ(sign*Tmp.getZ());          << 
305            Tmp *= toObserverFrameI;            << 
306            RightVector->operator[](hadronI)->S << 
307         }                                      << 
308                                                << 
309   return final_success;                        << 
310 }                                              << 
311                                                << 
312 //-------------------------------------------- << 
313                                                << 
314 G4bool G4LundStringFragmentation::StopFragment << 
315 {                                              << 
316   SetMinimalStringMass(string);                << 
317                                                << 
318   if ( MinimalStringMass < 0.) return true;    << 
319                                                << 
320   if (string->IsAFourQuarkString())            << 
321   {                                            << 
322     return G4UniformRand() < G4Exp(-0.0005*(st << 
323   } else {                                     << 
324                                                << 
325                 if (MinimalStringMass < 0.0 )  << 
326                                                << 
327     G4bool Result = G4UniformRand() <          << 
328         G4Exp(-0.66e-6*(string->Mass()*string- << 
329                 // G4bool Result = string->Mas << 
330                                                << 
331                 #ifdef debug_LUNDfragmentation << 
332                 G4cout<<"StopFragmenting Minim << 
333                       <<" "<<string->Mass()<<G << 
334                 G4cout<<"StopFragmenting - Yes << 
335                 #endif                         << 
336     return Result;                             << 
337   }                                            << 
338 }                                              << 
339                                                << 
340 //-------------------------------------------- << 
341                                                << 
342 G4KineticTrack * G4LundStringFragmentation::Sp << 
343                                   G4Fragmentin << 
344 {                                              << 
345        #ifdef debug_LUNDfragmentation          << 
346        G4cout<<G4endl;                         << 
347        G4cout<<"Start SplitUP ================ << 
348        G4cout<<"String partons: " <<string->Ge << 
349                                   <<string->Ge << 
350              <<"Direction "       <<string->Ge << 
351        #endif                                  << 
352                                                << 
353        //... random choice of string end to us << 
354        G4int SideOfDecay = (G4UniformRand() <  << 
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 * GetStr << 
370                                                << 
371        #ifdef debug_LUNDfragmentation          << 
372        G4cout<<"StrMass DiquarkSuppression     << 
373        #endif                                  << 
374                                                << 
375        G4int NumberOfpossibleBaryons = 2;      << 
376                                                << 
377        if (string->GetLeftParton()->GetParticl << 
378        if (string->GetRightParton()->GetPartic << 
379                                                << 
380        G4double ActualProb  = ProbDqADq ;      << 
381        ActualProb *= (1.0-G4Pow::GetInstance() << 
382        if(ActualProb <0.0) ActualProb = 0.;    << 
383                                                << 
384        SetDiquarkSuppression(ActualProb);      << 
385                                                << 
386        G4double Mth = 1250.0;                  << 
387        if ( NumberOfpossibleBaryons == 3 ){Mth << 
388        else if ( NumberOfpossibleBaryons == 4  << 
389        else {}                                 << 
390                                                << 
391        ActualProb = ProbSaS;                   << 
392        ActualProb *= (1.0 - G4Pow::GetInstance << 
393        if ( ActualProb < 0.0 ) ActualProb = 0. << 
394        SetStrangenessSuppression((1.0-ActualPr << 
395                                                << 
396        #ifdef debug_LUNDfragmentation          << 
397        G4cout<<"StrMass DiquarkSuppression cor << 
398        #endif                                  << 
399                                                << 
400        if (string->DecayIsQuark())             << 
401        {                                       << 
402           HadronDefinition= QuarkSplitup(strin << 
403        } else {                                << 
404           HadronDefinition= DiQuarkSplitup(str << 
405        }                                       << 
406                                                << 
407        SetDiquarkSuppression(ProbDqADq);       << 
408        SetStrangenessSuppression((1.0-ProbSaS) << 
409                                                << 
410        if ( HadronDefinition == NULL ) { G4Kin << 
411                                                << 
412        #ifdef debug_LUNDfragmentation          << 
413        G4cout<<"The parton "<<string->GetDecay << 
414              <<" produces hadron "<<HadronDefi << 
415              <<" and is transformed to "<<newS << 
416        G4cout<<"The side of the string decay L << 
417        #endif                                  << 
418        // create new String from old, ie. keep << 
419                                                << 
420        if ( newString ) delete newString;      << 
421                                                << 
422        newString=new G4FragmentingString(*stri << 
423                                                << 
424        #ifdef debug_LUNDfragmentation          << 
425        G4cout<<"An attempt to determine its en << 
426        #endif                                  << 
427        G4LorentzVector* HadronMomentum=SplitEa << 
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 << 
436            #endif                              << 
437      G4ThreeVector   Pos;                      << 
438      Hadron = new G4KineticTrack(HadronDefinit << 
439                                                << 
440            if ( newString ) delete newString;  << 
441                                                << 
442      newString=new G4FragmentingString(*string << 
443             HadronMomentum);                   << 
444      delete HadronMomentum;                    << 
445        }                                       << 
446        else                                    << 
447        {                                       << 
448            #ifdef debug_LUNDfragmentation      << 
449            G4cout<<"The attempt was not succes << 
450            #endif                              << 
451        }                                       << 
452                                                << 
453        #ifdef debug_LUNDfragmentation          << 
454        G4cout<<"End SplitUP (G4VLongitudinalSt << 
455        #endif                                  << 
456                                                << 
457        return Hadron;                          << 
458 }                                              << 
459                                                << 
460 //-------------------------------------------- << 
461                                                << 
462 G4ParticleDefinition * G4LundStringFragmentati << 
463                                                << 
464 {                                              << 
465    G4double StrSup=GetStrangeSuppress();       << 
466    G4double ProbQQbar = (1.0 - 2.0*StrSup)*1.2 << 
467                                                << 
468    //... can Diquark break or not?             << 
469    if (G4UniformRand() < DiquarkBreakProb ){   << 
470                                                << 
471       //... Diquark break                      << 
472       G4int stableQuarkEncoding = decay->GetPD << 
473       G4int decayQuarkEncoding = (decay->GetPD << 
474       if (G4UniformRand() < 0.5)               << 
475       {                                        << 
476          G4int Swap = stableQuarkEncoding;     << 
477          stableQuarkEncoding = decayQuarkEncod << 
478          decayQuarkEncoding = Swap;            << 
479       }                                        << 
480                                                << 
481       G4int IsParticle=(decayQuarkEncoding>0)  << 
482                                                << 
483       SetStrangenessSuppression((1.0-ProbQQbar << 
484       pDefPair QuarkPair = CreatePartonPair(Is << 
485       SetStrangenessSuppression((1.0-StrSup)/2 << 
486                                                << 
487       //... Build new Diquark                  << 
488       G4int QuarkEncoding=QuarkPair.second->Ge << 
489       G4int i10  = std::max(std::abs(QuarkEnco << 
490       G4int i20  = std::min(std::abs(QuarkEnco << 
491       G4int spin = (i10 != i20 && G4UniformRan << 
492       G4int NewDecayEncoding = -1*IsParticle*( << 
493       created = FindParticle(NewDecayEncoding) << 
494       G4ParticleDefinition * decayQuark=FindPa << 
495       G4ParticleDefinition * had=hadronizer->B << 
496       StrangeSuppress=StrSup;                  << 
497                                                << 
498       return had;                              << 
499                                                << 
500    } else {                                    << 
501       //... Diquark does not break             << 
502                                                << 
503       G4int IsParticle=(decay->GetPDGEncoding( << 
504                                                << 
505       StrangeSuppress=(1.0 - ProbQQbar)/2.0;   << 
506       pDefPair QuarkPair = CreatePartonPair(Is << 
507                                                << 
508       created = QuarkPair.second;              << 
509                                                << 
510       G4ParticleDefinition * had=hadronizer->B << 
511       StrangeSuppress=StrSup;                  << 
512                                                << 
513       return had;                              << 
514    }                                               40    }
515 }                                              << 
516                                                << 
517 //-------------------------------------------- << 
518                                                << 
519 G4LorentzVector * G4LundStringFragmentation::S << 
520                                                << 
521                                                << 
522 {                                              << 
523   G4LorentzVector String4Momentum=string->Get4 << 
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 n << 
531                                                << 
532         #ifdef debug_LUNDfragmentation         << 
533         G4cout<<G4endl<<"Start LUND SplitEandP << 
534         G4cout<<"String 4 mom, String M and Mt << 
535               <<" "<<std::sqrt(StringMT2)<<G4e << 
536         G4cout<<"Hadron "<<pHadron->GetParticl << 
537         G4cout<<"HadM MinimalStringMassLeft St << 
538               <<String4Momentum.mag()<<" "<<Ha << 
539         #endif                                 << 
540                                                << 
541   if ((HadronMass + MinimalStringMass > string << 
542   {                                            << 
543           #ifdef debug_LUNDfragmentation       << 
544           G4cout<<"Mass of the string is not s << 
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 mo << 
554   G4ThreeVector HadronPt    , RemSysPt;        << 
555   G4double      HadronMassT2, ResidualMassT2;  << 
556   G4double HadronMt, Pt, Pt2, phi;             << 
557                                                << 
558         G4double TmtCur = Tmt;                 << 
559                                                << 
560         if ( (string->GetDecayParton()->GetPar << 
561              (pHadron->GetBaryonNumber() != 0) << 
562           TmtCur = Tmt*0.37;              // q << 
563         } else if ( (string->GetDecayParton()- << 
564                     (pHadron->GetBaryonNumber( << 
565           //TmtCur = Tmt;                      << 
566   } else if ( (string->GetDecayParton()->GetPa << 
567                     (pHadron->GetBaryonNumber( << 
568           //TmtCur = Tmt*0.89;                 << 
569         } else if ( (string->GetDecayParton()- << 
570                     (pHadron->GetBaryonNumber( << 
571           TmtCur = Tmt*1.35;                   << 
572         }                                      << 
573                                                << 
574         //...  sample Pt of the hadron         << 
575         G4int attempt=0;                       << 
576         do                                     << 
577         {                                      << 
578           attempt++; if (attempt > StringLoopI << 
579                                                << 
580           HadronMt = HadronMass - TmtCur*G4Log << 
581     Pt2 = sqr(HadronMt)-sqr(HadronMass); Pt=st << 
582     phi = 2.*pi*G4UniformRand();               << 
583           HadronPt = G4ThreeVector( Pt*std::co << 
584           RemSysPt = StringPt - HadronPt;      << 
585           HadronMassT2 = sqr(HadronMass) + Had << 
586           ResidualMassT2=sqr(MinimalStringMass << 
587                                                << 
588         } while (std::sqrt(HadronMassT2) + std << 
589                                                << 
590   //...  sample z to define hadron longitudina << 
591   //... but first check the available phase sp << 
592                                                << 
593   G4double Pz2 = (sqr(StringMT2 - HadronMassT2 << 
594       4*HadronMassT2 * ResidualMassT2)/4./Stri << 
595                                                << 
596   if (Pz2 < 0 ) {return 0;}          // have t << 
597                                                << 
598   //... then compute allowed z region  z_min < << 
599                                                << 
600   G4double Pz = std::sqrt(Pz2);                << 
601   G4double zMin = (std::sqrt(HadronMassT2+Pz2) << 
602         // G4double zMin = (std::sqrt(HadronMa << 
603   G4double zMax = (std::sqrt(HadronMassT2+Pz2) << 
604                                                << 
605   if (zMin >= zMax) return 0;   // have to sta << 
606                                                << 
607   G4double z = GetLightConeZ(zMin, zMax,       << 
608                  string->GetDecayParton()->Get << 
609                  HadronPt.x(), HadronPt.y());  << 
610                                                << 
611   //... now compute hadron longitudinal moment << 
612   // longitudinal hadron momentum component Ha << 
613                                                << 
614   HadronPt.setZ(0.5* string->GetDecayDirection << 
615            (z * string->LightConeDecay() -     << 
616           HadronMassT2/(z * string->LightConeD << 
617   G4double HadronE  = 0.5* (z * string->LightC << 
618           HadronMassT2/(z * string->LightConeD << 
619                                                << 
620   G4LorentzVector * a4Momentum= new G4LorentzV << 
621                                                << 
622         #ifdef debug_LUNDfragmentation         << 
623         G4cout<<G4endl<<" string->GetDecayDire << 
624         G4cout<<"string->LightConeDecay() "<<s << 
625         G4cout<<"HadronPt,HadronE "<<HadronPt< << 
626         G4cout<<"String4Momentum "<<String4Mom << 
627         G4cout<<"Out of LUND SplitEandP "<<G4e << 
628         #endif                                 << 
629                                                    41 
630   return a4Momentum;                           <<  42 // G4LundStringFragmentation::G4LundStringFragmentation(G4double sigmaPt)
631 }                                              <<  43 // : G4VLongitudinalStringDecay(sigmaPt)
                                                   >>  44 //    {
                                                   >>  45 //    }
632                                                    46 
633 //-------------------------------------------- <<  47 G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
634                                                <<  48    {
635 G4double G4LundStringFragmentation::GetLightCo <<  49    }
636                                       G4int PD << 
637                                       G4Partic << 
638                                       G4double << 
639 {                                              << 
640   G4double Mass = pHadron->GetPDGMass();       << 
641         G4int HadronEncoding=std::abs(pHadron- << 
642                                                << 
643   G4double Mt2 = Px*Px + Py*Py + Mass*Mass;    << 
644                                                << 
645   G4double  Alund, Blund;                      << 
646   G4double zOfMaxyf(0.), maxYf(1.), z(0.), yf( << 
647                                                << 
648   if (!((std::abs(PDGEncodingOfDecayParton) >  << 
649   {    // ---------------- Quark fragmentation << 
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(sq << 
658      }                                         << 
659                                                << 
660      if (zOfMaxyf < zmin) {zOfMaxyf=zmin;}     << 
661      if (zOfMaxyf > zmax) {zOfMaxyf=zmax;}     << 
662      maxYf=(1-zOfMaxyf)/zOfMaxyf * G4Exp(-Blun << 
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* << 
670     yf = G4Pow::GetInstance()->powA(1.0-z,Alun << 
671      }                                         << 
672      while ( (G4UniformRand()*maxYf > yf) && + << 
673            if ( loopCounter >= maxNumberOfLoop << 
674              z = 0.5*(zmin + zmax);  // Just a << 
675            }                                   << 
676      return z;                                 << 
677         }                                      << 
678                                                << 
679   if (std::abs(PDGEncodingOfDecayParton) > 100 << 
680   {                                            << 
681                 G4double an = 2.5;             << 
682                 an +=(sqr(Px)+sqr(Py))/sqr(GeV << 
683                 z=zmin + (zmax-zmin)*G4Pow::Ge << 
684                 if( PDGEncodingOfDecayParton > << 
685   }                                            << 
686                                                << 
687   return z;                                    << 
688 }                                              << 
689                                                << 
690 //-------------------------------------------- << 
691                                                << 
692 G4bool G4LundStringFragmentation::SplitLast(G4 << 
693                                             G4 << 
694                                             G4 << 
695 {                                              << 
696   //... perform last cluster decay             << 
697         SetMinimalStringMass( string);         << 
698         if ( MinimalStringMass < 0.) return fa << 
699         #ifdef debug_LUNDfragmentation         << 
700         G4cout<<G4endl<<"Split last----------- << 
701         G4cout<<"MinimalStringMass "<<MinimalS << 
702         G4cout<<"Left  "<<string->GetLeftParto << 
703         G4cout<<"Right "<<string->GetRightPart << 
704         G4cout<<"String4mom "<<string->GetPstr << 
705         #endif                                 << 
706                                                << 
707         G4LorentzVector Str4Mom=string->Get4Mo << 
708         G4LorentzRotation toCms(-1*Str4Mom.boo << 
709         G4LorentzVector Pleft = toCms * string << 
710         toCms.rotateZ(-1*Pleft.phi());         << 
711         toCms.rotateY(-1*Pleft.theta());       << 
712                                                << 
713         G4LorentzRotation toObserverFrame= toC << 
714                                                << 
715   G4double StringMass=string->Mass();          << 
716                                                << 
717   G4ParticleDefinition * LeftHadron(0), * Righ << 
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()->GetPa << 
727               <<string->GetRightParton()->GetP << 
728         #endif                                 << 
729                                                << 
730   string->SetLeftPartonStable(); // to query q << 
731                                                << 
732   if (string->IsAFourQuarkString() )           << 
733   {                                            << 
734           G4int IDleft =std::abs(string->GetLe << 
735           G4int IDright=std::abs(string->GetRi << 
736                                                << 
737           if ( (IDleft > 3000) || (IDright > 3 << 
738             if ( ! Diquark_AntiDiquark_belowTh << 
739             {                                  << 
740               return false;                    << 
741             }                                  << 
742           } else {                             << 
743     // The string is qq-qqbar type. Diquarks a << 
744           if (StringMass-MinimalStringMass < 0 << 
745     {                                          << 
746       if (! Diquark_AntiDiquark_belowThreshold << 
747                         {                      << 
748         return false;                          << 
749                         }                      << 
750     } else                                     << 
751     {                                          << 
752       Diquark_AntiDiquark_aboveThreshold_lastS << 
753                                                << 
754       if (NumberOf_FS == 0) return false;      << 
755                                                << 
756                         sampledState = SampleS << 
757       if (string->GetLeftParton()->GetPDGEncod << 
758       {                                        << 
759         LeftHadron =FS_LeftHadron[sampledState << 
760         RightHadron=FS_RightHadron[sampledStat << 
761       } else                                   << 
762       {                                        << 
763         LeftHadron =FS_RightHadron[sampledStat << 
764         RightHadron=FS_LeftHadron[sampledState << 
765       }                                        << 
766     }                                          << 
767           }  // ID > 3300                      << 
768         } else {                               << 
769     if (string->DecayIsQuark() && string->Stab << 
770     {       //... there are quarks on cluster  << 
771                         #ifdef debug_LUNDfragm << 
772                         G4cout<<"Q Q string La << 
773                         #endif                 << 
774                                                << 
775       Quark_AntiQuark_lastSplitting(string, Le << 
776                                                << 
777       if (NumberOf_FS == 0) return false;      << 
778                   sampledState = SampleState() << 
779       if (string->GetLeftParton()->GetPDGEncod << 
780       {                                        << 
781         LeftHadron =FS_RightHadron[sampledStat << 
782         RightHadron=FS_LeftHadron[sampledState << 
783       } else                                   << 
784       {                                        << 
785         LeftHadron =FS_LeftHadron[sampledState << 
786         RightHadron=FS_RightHadron[sampledStat << 
787       }                                        << 
788     } else                                     << 
789     {       //... there is a Diquark on one of << 
790                         #ifdef debug_LUNDfragm << 
791                         G4cout<<"DiQ Q string  << 
792                         #endif                 << 
793                                                << 
794       Quark_Diquark_lastSplitting(string, Left << 
795                                                << 
796       if (NumberOf_FS == 0) return false;      << 
797                   sampledState = SampleState() << 
798                                                << 
799       if (string->GetLeftParton()->GetParticle << 
800       {                                        << 
801         LeftHadron =FS_LeftHadron[sampledState << 
802         RightHadron=FS_RightHadron[sampledStat << 
803       } else                                   << 
804       {                                        << 
805         LeftHadron =FS_RightHadron[sampledStat << 
806         RightHadron=FS_LeftHadron[sampledState << 
807       }                                        << 
808     }                                          << 
809                                                << 
810   }                                            << 
811                                                << 
812         #ifdef debug_LUNDfragmentation         << 
813         G4cout<<"Sampled hadrons: "<<LeftHadro << 
814         #endif                                 << 
815                                                << 
816   G4LorentzVector P_left  =string->GetPleft(), << 
817                                                << 
818   G4LorentzVector  LeftMom, RightMom;          << 
819   G4ThreeVector    Pos;                        << 
820                                                << 
821   Sample4Momentum(&LeftMom,  LeftHadron->GetPD << 
822       &RightMom, RightHadron->GetPDGMass(),    << 
823       StringMass);                             << 
824                                                << 
825         // Sample4Momentum ascribes LeftMom.pz << 
826         // It must be negative in case when th << 
827                                                << 
828   if (!(string->DecayIsQuark() && string->Stab << 
829   { // Only for qq - q, q - qq, and qq - qqbar << 
830                                                << 
831           if ( G4UniformRand() <= 0.5 )        << 
832     {                                          << 
833       if (P_left.z() <= 0.) {G4LorentzVector t << 
834     }                                          << 
835     else                                       << 
836     {                                          << 
837       if (P_right.z() >= 0.) {G4LorentzVector  << 
838     }                                          << 
839   }                                            << 
840                                                << 
841   LeftMom *=toObserverFrame;                   << 
842   RightMom*=toObserverFrame;                   << 
843                                                << 
844   LeftVector->push_back(new G4KineticTrack(Lef << 
845   RightVector->push_back(new G4KineticTrack(Ri << 
846                                                << 
847   string->LorentzRotate(toObserverFrame);      << 
848   return true;                                 << 
849 }                                              << 
850                                                << 
851 //-------------------------------------------- << 
852                                                << 
853 G4bool G4LundStringFragmentation::             << 
854 Diquark_AntiDiquark_belowThreshold_lastSplitti << 
855                                                << 
856                                                << 
857 {                                              << 
858   G4double StringMass   = string->Mass();      << 
859                                                << 
860   G4int cClusterInterrupt = 0;                 << 
861         G4bool isOK = false;                   << 
862   do                                           << 
863   {                                            << 
864     G4int LeftQuark1= string->GetLeftParton()- << 
865     G4int LeftQuark2=(string->GetLeftParton()- << 
866                                                << 
867     G4int RightQuark1= string->GetRightParton( << 
868     G4int RightQuark2=(string->GetRightParton( << 
869                                                << 
870     if (G4UniformRand()<0.5)                   << 
871     {                                          << 
872       LeftHadron =hadronizer->Build(FindPartic << 
873                   FindParticle(RightQuark1));  << 
874       RightHadron= (LeftHadron == nullptr) ? n << 
875                                                << 
876                   FindParticle(RightQuark2));  << 
877     } else                                     << 
878     {                                          << 
879       LeftHadron =hadronizer->Build(FindPartic << 
880                   FindParticle(RightQuark2));  << 
881       RightHadron=(LeftHadron == nullptr) ? nu << 
882                                     hadronizer << 
883                   FindParticle(RightQuark1));  << 
884     }                                          << 
885                                                << 
886     isOK = (LeftHadron != nullptr) && (RightHa << 
887                                                << 
888     if(isOK) { isOK = (StringMass > LeftHadron << 
889     ++cClusterInterrupt;                       << 
890     //... repeat procedure, if mass of cluster << 
891     //... ClusterMassCut = 0.15*GeV model para << 
892   }                                            << 
893   while (isOK == false && cClusterInterrupt <  << 
894   /* Loop checking, 07.08.2015, A.Ribon */     << 
895     return isOK;                               << 
896 }                                              << 
897                                                << 
898 //-------------------------------------------- << 
899                                                << 
900 G4bool G4LundStringFragmentation::             << 
901 Diquark_AntiDiquark_aboveThreshold_lastSplitti << 
902                                                << 
903                                                << 
904 {                                              << 
905   // StringMass-MinimalStringMass > 0. Creatio << 
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( << 
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->Get << 
923   G4int AbsIDAnti_di_quark =std::abs(IDAnti_di << 
924   G4int IDdi_quark         =Di_Quark->GetPDGEn << 
925   G4int AbsIDdi_quark      =std::abs(IDdi_quar << 
926                                                << 
927   G4int ADi_q1=AbsIDAnti_di_quark/1000;        << 
928   G4int ADi_q2=(AbsIDAnti_di_quark-ADi_q1*1000 << 
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 = << 
938                 G4int loopCounter = 0;         << 
939     do  // while(Meson[AbsIDquark-1][ProdQ-1][ << 
940     {                                          << 
941       LeftHadron=G4ParticleTable::GetParticleT << 
942               -Baryon[ADi_q1-1][ADi_q2-1][Prod << 
943                                                << 
944       if (LeftHadron == NULL) continue;        << 
945       G4double LeftHadronMass=LeftHadron->GetP << 
946                                                << 
947       G4int StateDiQ=0;                        << 
948                         const G4int maxNumberO << 
949                         G4int internalLoopCoun << 
950       do // while(Baryon[Di_q1-1][Di_q2-1][Pro << 
951       {                                        << 
952         RightHadron=G4ParticleTable::GetPartic << 
953                 +Baryon[Di_q1-1][Di_q2-1][Prod << 
954                                                << 
955         if (RightHadron == NULL) continue;     << 
956         G4double RightHadronMass=RightHadron-> << 
957                                                << 
958         if (StringMass > LeftHadronMass + Righ << 
959         {                                      << 
960                                         if ( N << 
961                                           G4Ex << 
962                                           ed < << 
963                                           G4Ex << 
964                                                << 
965                                           Numb << 
966                                         }      << 
967                                                << 
968           G4double FS_Psqr=lambda(StringMassSq << 
969                 sqr(RightHadronMass));         << 
970           //FS_Psqr=1.;                        << 
971           FS_Weight[NumberOf_FS]=std::sqrt(FS_ << 
972                      BaryonWeight[ADi_q1-1][AD << 
973                      BaryonWeight[Di_q1-1][Di_ << 
974                      Prob_QQbar[ProdQ-1];      << 
975                                                << 
976           FS_LeftHadron[NumberOf_FS] = LeftHad << 
977           FS_RightHadron[NumberOf_FS]= RightHa << 
978           NumberOf_FS++;                       << 
979         } // End of if (StringMass > LeftHadro << 
980                                                << 
981         StateDiQ++;                            << 
982                                                << 
983       } while( (Baryon[Di_q1-1][Di_q2-1][ProdQ << 
984                                  ++internalLoo << 
985                         if ( internalLoopCount << 
986                           return false;        << 
987                         }                      << 
988                                                << 
989       StateADiQ++;                             << 
990     } while( (Baryon[ADi_q1-1][ADi_q2-1][ProdQ << 
991                          ++loopCounter < maxNu << 
992                 if ( loopCounter >= maxNumberO << 
993                   return false;                << 
994                 }                              << 
995   } // End of for (G4int ProdQ=1; ProdQ < 4; P << 
996                                                << 
997     return true;                               << 
998 }                                              << 
999                                                << 
1000 //------------------------------------------- << 
1001                                               << 
1002 G4bool G4LundStringFragmentation::Quark_Diqua << 
1003                                               << 
1004                                               << 
1005 {                                             << 
1006   G4double StringMass   = string->Mass();     << 
1007   G4double StringMassSqr= sqr(StringMass);    << 
1008                                               << 
1009   G4ParticleDefinition * Di_Quark;            << 
1010   G4ParticleDefinition * Quark;               << 
1011                                               << 
1012   if (string->GetLeftParton()->GetParticleSub << 
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->GetPDGEncodin << 
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++)  // << 
1033   {                                        // << 
1034     G4int SignQ;                              << 
1035     if (IDquark > 0)                          << 
1036     {                                         << 
1037             SignQ=-1;                         << 
1038             if (IDquark == 2)                 << 
1039       if ((IDquark == 1) && (ProdQ == 3)) Sig << 
1040       if ((IDquark == 3) && (ProdQ == 1)) Sig << 
1041             if (IDquark == 4)                 << 
1042       if (IDquark == 5)                   Sig << 
1043     } else                                    << 
1044     {                                         << 
1045       SignQ= 1;                               << 
1046       if (IDquark == -2)                  Sig << 
1047       if ((IDquark ==-1) && (ProdQ == 3)) Sig << 
1048       if ((IDquark ==-3) && (ProdQ == 1)) Sig << 
1049       if (IDquark == -4)                  Sig << 
1050       if (IDquark == -5)                  Sig << 
1051     }                                         << 
1052                                               << 
1053     if (AbsIDquark == ProdQ)            SignQ << 
1054                                               << 
1055     G4int StateQ=0;                           << 
1056                 const G4int maxNumberOfLoops  << 
1057                 G4int loopCounter = 0;        << 
1058     do  // while(Meson[AbsIDquark-1][ProdQ-1] << 
1059     {                                         << 
1060       LeftHadron=G4ParticleTable::GetParticle << 
1061               Meson[AbsIDquark-1][ProdQ-1][St << 
1062       if (LeftHadron == NULL) continue;       << 
1063       G4double LeftHadronMass=LeftHadron->Get << 
1064                                               << 
1065       G4int StateDiQ=0;                       << 
1066                         const G4int maxNumber << 
1067                         G4int internalLoopCou << 
1068       do // while(Baryon[Di_q1-1][Di_q2-1][Pr << 
1069       {                                       << 
1070         RightHadron=G4ParticleTable::GetParti << 
1071                 Baryon[Di_q1-1][Di_q2-1][Prod << 
1072         if (RightHadron == NULL) continue;    << 
1073         G4double RightHadronMass=RightHadron- << 
1074                                               << 
1075         if (StringMass > LeftHadronMass + Rig << 
1076         {                                     << 
1077                                         if (  << 
1078                                           G4E << 
1079                                           ed  << 
1080                                           G4E << 
1081                                               << 
1082                                           Num << 
1083                                         }     << 
1084                                               << 
1085           G4double FS_Psqr=lambda(StringMassS << 
1086                 sqr(RightHadronMass));        << 
1087           FS_Weight[NumberOf_FS]=std::sqrt(FS << 
1088                      MesonWeight[AbsIDquark-1 << 
1089                      BaryonWeight[Di_q1-1][Di << 
1090                      Prob_QQbar[ProdQ-1];     << 
1091                                               << 
1092           FS_LeftHadron[NumberOf_FS] = LeftHa << 
1093           FS_RightHadron[NumberOf_FS]= RightH << 
1094                                               << 
1095           NumberOf_FS++;                      << 
1096         } // End of if (StringMass > LeftHadr << 
1097                                               << 
1098         StateDiQ++;                           << 
1099                                               << 
1100       } while( (Baryon[Di_q1-1][Di_q2-1][Prod << 
1101                                  ++internalLo << 
1102                         if ( internalLoopCoun << 
1103                           return false;       << 
1104                         }                     << 
1105                                               << 
1106       StateQ++;                               << 
1107     } while( (Meson[AbsIDquark-1][ProdQ-1][St << 
1108                          ++loopCounter < maxN << 
1109                                               << 
1110                   if ( loopCounter >= maxNumb << 
1111                     return false;             << 
1112                   }                           << 
1113   }                                           << 
1114                                               << 
1115   return true;                                << 
1116 }                                             << 
1117                                               << 
1118 //------------------------------------------- << 
1119                                               << 
1120 G4bool G4LundStringFragmentation::Quark_AntiQ << 
1121                                               << 
1122                                               << 
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 << 
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->GetPDGEncodin << 
1141   G4int AbsIDquark      =std::abs(IDquark);   << 
1142         G4int QuarkCharge     =Qcharge[IDquar << 
1143                                               << 
1144   G4int IDanti_quark    =Anti_Quark->GetPDGEn << 
1145   G4int AbsIDanti_quark =std::abs(IDanti_quar << 
1146         G4int AntiQuarkCharge =-Qcharge[AbsID << 
1147                                               << 
1148         G4int LeftHadronCharge(0), RightHadro << 
1149                                               << 
1150         //G4cout<<"Q Qbar "<<IDquark<<" "<<ID << 
1151                                               << 
1152   NumberOf_FS=0;                              << 
1153   for (G4int ProdQ=1; ProdQ < 4; ProdQ++)  // << 
1154   {                                        // << 
1155                 //G4cout<<"NumberOf_FS ProdQ  << 
1156     LeftHadronCharge = QuarkCharge - Qcharge[ << 
1157     G4int SignQ = LeftHadronCharge/3; if (Sig << 
1158                                               << 
1159     if ((IDquark == 1) && (ProdQ == 3)) SignQ << 
1160     if ((IDquark == 3) && (ProdQ == 1)) SignQ << 
1161                 if ((IDquark == 4) && (ProdQ  << 
1162                 if ((IDquark == 5) && (ProdQ  << 
1163                 if ((IDquark == 5) && (ProdQ  << 
1164                                               << 
1165                 RightHadronCharge = AntiQuark << 
1166     G4int SignAQ = RightHadronCharge/3; if (S << 
1167                                               << 
1168     if ((IDanti_quark ==-1) && (ProdQ == 3))  << 
1169     if ((IDanti_quark ==-3) && (ProdQ == 1))  << 
1170     if ((IDanti_quark ==-4) && (ProdQ == 2))  << 
1171     if ((IDanti_quark ==-5) && (ProdQ == 1))  << 
1172     if ((IDanti_quark ==-5) && (ProdQ == 3))  << 
1173                                               << 
1174                 //G4cout<<"ProQ signs "<<Prod << 
1175                                               << 
1176     G4int StateQ=0;                           << 
1177                 const G4int maxNumberOfLoops  << 
1178                 G4int loopCounter = 0;        << 
1179     do                                        << 
1180     {                                         << 
1181                         //G4cout<<"[AbsIDquar << 
1182                         //<<ProdQ-1<<" "<<Sta << 
1183       LeftHadron=G4ParticleTable::GetParticle << 
1184                    Meson[AbsIDquark-1][ProdQ- << 
1185                         //G4cout<<"LeftHadron << 
1186       if (LeftHadron == NULL) { StateQ++; con << 
1187                         //G4cout<<"LeftHadron << 
1188       G4double LeftHadronMass=LeftHadron->Get << 
1189                                               << 
1190       G4int StateAQ=0;                        << 
1191                         const G4int maxNumber << 
1192                         G4int internalLoopCou << 
1193       do                                      << 
1194       {                                       << 
1195                                 //G4cout<<"   << 
1196                                 //      <<Pro << 
1197         RightHadron=G4ParticleTable::GetParti << 
1198                 Meson[AbsIDanti_quark-1][Prod << 
1199                                 //G4cout<<"Ri << 
1200         if(RightHadron == NULL) { StateAQ++;  << 
1201                                 //G4cout<<"Ri << 
1202         G4double RightHadronMass=RightHadron- << 
1203                                               << 
1204         if (StringMass > LeftHadronMass + Rig << 
1205         {                                     << 
1206                                         if (  << 
1207                                           G4E << 
1208                                           ed  << 
1209                                           G4E << 
1210                                               << 
1211                                           Num << 
1212                                         }     << 
1213                                               << 
1214                                         G4dou << 
1215                 sqr(RightHadronMass));        << 
1216           //FS_Psqr=1.;                       << 
1217           FS_Weight[NumberOf_FS]=std::sqrt(FS << 
1218                      MesonWeight[AbsIDquark-1 << 
1219                      MesonWeight[AbsIDanti_qu << 
1220                      Prob_QQbar[ProdQ-1];     << 
1221           if (string->GetLeftParton()->GetPDG << 
1222           {                                   << 
1223             FS_LeftHadron[NumberOf_FS] = Righ << 
1224             FS_RightHadron[NumberOf_FS]= Left << 
1225           } else                              << 
1226           {                                   << 
1227             FS_LeftHadron[NumberOf_FS] = Left << 
1228             FS_RightHadron[NumberOf_FS]= Righ << 
1229           }                                   << 
1230                                               << 
1231           NumberOf_FS++;                      << 
1232                                               << 
1233         }                                     << 
1234                                               << 
1235         StateAQ++;                            << 
1236                                 //G4cout<<"   << 
1237                                 //      <<Mes << 
1238       } while ( (Meson[AbsIDanti_quark-1][Pro << 
1239                                   ++internalL << 
1240                           if ( internalLoopCo << 
1241                             return false;     << 
1242                           }                   << 
1243                                               << 
1244       StateQ++;                               << 
1245                         //G4cout<<"StateQ Mes << 
1246                         //      <<Meson[AbsID << 
1247                                               << 
1248     } while ( (Meson[AbsIDquark-1][ProdQ-1][S << 
1249                          ++loopCounter < maxN << 
1250                   if ( loopCounter >= maxNumb << 
1251                     return false;             << 
1252                   }                           << 
1253   } // End of for (G4int ProdQ=1; ProdQ < 4;  << 
1254                                               << 
1255   return true;                                << 
1256 }                                             << 
1257                                               << 
1258 //------------------------------------------- << 
1259                                               << 
1260 G4int G4LundStringFragmentation::SampleState( << 
1261 {                                             << 
1262         if ( NumberOf_FS > 349 ) {            << 
1263           G4ExceptionDescription ed;          << 
1264           ed << " NumberOf_FS exceeds its lim << 
1265           G4Exception( "G4LundStringFragmenta << 
1266           NumberOf_FS = 349;                  << 
1267         }                                     << 
1268                                               << 
1269   G4double SumWeights=0.;                     << 
1270   for (G4int i=0; i<NumberOf_FS; i++) {SumWei << 
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::Sample4Moment << 
1288                                               << 
1289                                               << 
1290 {                                             << 
1291   // ------ Sampling of momenta of 2 last pro << 
1292   G4ThreeVector Pt;                           << 
1293   G4double MassMt, AntiMassMt;                << 
1294   G4double AvailablePz, AvailablePz2;         << 
1295         #ifdef debug_LUNDfragmentation        << 
1296         G4cout<<"Sampling of momenta of 2 las << 
1297         G4cout<<"Init Mass "<<InitialMass<<"  << 
1298         #endif                                << 
1299                                               << 
1300   G4double r_val = sqr(InitialMass*InitialMas << 
1301       sqr(2.*Mass*AntiMass);                  << 
1302   G4double Pabs = (r_val > 0.)? std::sqrt(r_v << 
1303                                               << 
1304         const G4int maxNumberOfLoops = 1000;  << 
1305   G4double SigmaQTw = SigmaQT;                << 
1306   if ( Mass > 930. || AntiMass > 930. ) {     << 
1307     SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+AntiMa << 
1308         }                                     << 
1309         if ( Mass < 930. && AntiMass < 930. ) << 
1310         if ( ( Mass < 930. && AntiMass > 930. << 
1311        ( Mass > 930. && AntiMass < 930. ) ) { << 
1312           //SigmaQT = -1.;  // isotropical de << 
1313         }                                     << 
1314         if ( Mass > 930. && AntiMass > 930. ) << 
1315           SigmaQT *= ( 1.0 - 0.55*sqr( (Mass+ << 
1316         }                                     << 
1317                                               << 
1318         G4int loopCounter = 0;                << 
1319   do                                          << 
1320   {                                           << 
1321     Pt=SampleQuarkPt(Pabs); Pt.setZ(0); G4dou << 
1322     MassMt    = std::sqrt(    Mass *     Mass << 
1323     AntiMassMt= std::sqrt(AntiMass * AntiMass << 
1324   }                                           << 
1325   while ( (InitialMass < MassMt + AntiMassMt) << 
1326                                                   50 
1327         SigmaQT = SigmaQTw;                   << 
1328                                                   51 
1329   AvailablePz2= sqr(InitialMass*InitialMass - <<  52 G4LundStringFragmentation::~G4LundStringFragmentation()
1330         4.*sqr(MassMt*AntiMassMt);            <<  53    { 
                                                   >>  54    }
1331                                                   55 
1332   AvailablePz2 /=(4.*InitialMass*InitialMass) <<  56 //****************************************************************************************
1333   AvailablePz = std::sqrt(AvailablePz2);      << 
1334                                                   57 
1335   G4double Px=Pt.getX();                      <<  58 const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
1336   G4double Py=Pt.getY();                      <<  59    {
                                                   >>  60      throw G4HadronicException(__FILE__, __LINE__, "G4LundStringFragmentation::operator= meant to not be accessable");
                                                   >>  61      return *this;
                                                   >>  62    }
1337                                                   63 
1338   Mom->setPx(Px); Mom->setPy(Py); Mom->setPz( <<  64 int G4LundStringFragmentation::operator==(const G4LundStringFragmentation &right) const
1339   Mom->setE(std::sqrt(sqr(MassMt)+AvailablePz <<  65    {
                                                   >>  66    return !memcmp(this, &right, sizeof(G4LundStringFragmentation));
                                                   >>  67    }
1340                                                   68 
1341   AntiMom->setPx(-Px); AntiMom->setPy(-Py); A <<  69 int G4LundStringFragmentation::operator!=(const G4LundStringFragmentation &right) const
1342   AntiMom->setE (std::sqrt(sqr(AntiMassMt)+Av <<  70    {
                                                   >>  71    return memcmp(this, &right, sizeof(G4LundStringFragmentation));
                                                   >>  72    }
1343                                                   73 
1344         #ifdef debug_LUNDfragmentation        <<  74 //****************************************************************************************
1345         G4cout<<"Fmass Mom "<<Mom->getX()<<"  << 
1346         G4cout<<"Smass Mom "<<AntiMom->getX() << 
1347               <<" "<<AntiMom->getT()<<G4endl; << 
1348         #endif                                << 
1349 }                                             << 
1350                                                   75 
1351 //------------------------------------------- <<  76 G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int ,  G4ParticleDefinition* pHadron, G4double Px, G4double Py)
                                                   >>  77     {
                                                   >>  78     const G4double  alund = 0.7/GeV/GeV; 
1352                                                   79 
1353 G4double G4LundStringFragmentation::lambda(G4 <<  80 //    If blund get restored, you MUST adapt the calculation of zOfMaxyf.
1354 {                                             <<  81 //    const G4double  blund = 1;
1355   G4double lam = sqr(S - m1_Sqr - m2_Sqr) - 4 << 
1356   return lam;                                 << 
1357 }                                             << 
1358                                                   82 
1359 // ------------------------------------------ <<  83     G4double z, yf;
1360 G4LundStringFragmentation::~G4LundStringFragm <<  84     G4double Mass = pHadron->GetPDGMass();
1361 {}                                            <<  85     
                                                   >>  86     G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
                                                   >>  87     G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
                                                   >>  88     G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * exp(-alund*Mt2/zOfMaxyf);
                                                   >>  89     do
                                                   >>  90        {
                                                   >>  91        z = zmin + G4UniformRand()*(zmax-zmin);
                                                   >>  92 //       yf = pow(1. - z, blund)/z*exp(-alund*Mt2/z);
                                                   >>  93    yf = (1-z)/z * exp(-alund*Mt2/z);
                                                   >>  94        } 
                                                   >>  95     while (G4UniformRand()*maxYf > yf); 
                                                   >>  96     return z;
                                                   >>  97     }
1362                                                   98 
                                                   >>  99 //****************************************************************************************
1363                                                  100